« back to all changes in this revision
Viewing changes to R/family.rrr.q
- browse files at revision 9
- compare with another revision
- download diff
- download tarball
- view history from revision 9
- Committer: Bazaar Package Importer
- Author(s): Chris Lawrence
- Date: 2009-07-10 21:27:58 UTC
- mfrom: (1.1.7 upstream)
- Revision ID: james.westby@ubuntu.com-20090710212758-q06yias6vzajwe54
Tags: 0.7-9-1
* New upstream release.
* Update debian/watch.
* Update debian/watch.
- files added:
- R/formula.vlm.q
- files removed:
- R/build.terms.s
- files modified:
- DESCRIPTION
added
removed
R/family.rrr.q
5
5
6
6
7
7
8
replace.constraints <- function(Blist, cm, index)
9
{
8
replace.constraints = function(Blist, cm, index) {
10
9
11
for(i in index)
12
Blist[[i]] = cm
10
for(iii in index)
11
Blist[[iii]] = cm
13
12
Blist
14
13
}
15
14
18
17
Alphavec=c(2, 4, 6, 9, 12, 16, 20, 25, 30, 40, 50,
19
18
60, 80, 100, 125, 2^(8:12)),
20
19
Criterion = c("rss", "coefficients"),
21
Linesearch= FALSE, Maxit=7,
20
Linesearch = FALSE, Maxit=7,
22
21
Suppress.warning=TRUE,
23
22
Tolerance=1e-7, ...)
24
23
{
68
67
Maxit=20,
69
68
Structural.zero=NULL,
70
69
SD.Cinit=0.02,
71
Suppress.warning= FALSE,
70
Suppress.warning=FALSE,
72
71
Tolerance=1e-6,
73
trace= FALSE,
72
trace=FALSE,
74
73
xij=NULL)
75
74
{
76
75
110
109
if(length(Structural.zero))
111
110
diag(M)[,-Structural.zero,drop=FALSE] else diag(M), 1:Rank)
112
111
if(p1)
113
for(k in 1:p1)
114
cmat2[[Rank+k]] <- Blist[[colx1.index[k]]]
112
for(kk in 1:p1)
113
cmat2[[Rank+kk]] <- Blist[[colx1.index[kk]]]
115
114
116
115
if(is.null(Cinit))
117
116
Cinit <- matrix(rnorm(p2*Rank, sd=SD.Cinit), p2, Rank)
127
126
128
127
lv.mat <- x[,colx2.index,drop=FALSE] %*% C
129
128
new.lv.model.matrix = cbind(lv.mat, if(p1) x[,colx1.index] else NULL)
130
fit = vlm.wfit(new.lv.model.matrix, z, Blist=cmat2, U=U,
131
matrix.out=TRUE, XBIG=FALSE, rss=FALSE, qr=FALSE, xij=xij)
129
fit = vlm.wfit(xmat=new.lv.model.matrix, z, Blist=cmat2, U=U,
130
matrix.out=TRUE, is.vlmX=FALSE, rss=FALSE, qr=FALSE, xij=xij)
132
131
A <- t(fit$mat.coef[1:Rank,,drop=FALSE])
133
132
134
133
cmat1 = replace.constraints(Blist, A, colx2.index)
135
fit = vlm.wfit(x, z, Blist=cmat1, U=U,
136
matrix.out=TRUE, XBIG= FALSE, rss=TRUE, qr= FALSE, xij=xij)
134
fit = vlm.wfit(xmat=x, z, Blist=cmat1, U=U,
135
matrix.out=TRUE, is.vlmX=FALSE, rss=TRUE, qr=FALSE, xij=xij)
137
136
C = fit$mat.coef[colx2.index,,drop=FALSE] %*% A %*% solve(t(A) %*% A)
138
137
139
138
numat = x[,colx2.index,drop=FALSE] %*% C
156
155
" ratio =", format(ratio), "\n")
157
156
if(!is.null(fit$rss))
158
157
cat(" rss =", fit$rss, "\n")
159
if(exists("flush.console"))
160
flush.console()
158
flush.console()
161
159
}
162
160
163
161
if(ratio < Tolerance) {
180
178
try.new.lv.model.matrix = cbind(try.lv.mat,
181
179
if(p1) x[,colx1.index] else NULL)
182
180
183
try = vlm.wfit(try.new.lv.model.matrix, z, Blist=cmat2, U=U,
184
matrix.out=TRUE, XBIG= FALSE, rss=TRUE, qr= FALSE,
185
xij=xij)
181
try = vlm.wfit(xmat=try.new.lv.model.matrix, z, Blist=cmat2,
182
U=U, matrix.out=TRUE, is.vlmX=FALSE,
183
rss=TRUE, qr=FALSE, xij=xij)
186
184
if(try$rss < ftemp) {
187
185
use.alpha <- Alphavec[itter]
188
186
fit <- try
195
193
if(trace && use.alpha>0) {
196
194
cat(" Finished line search using Alpha =",
197
195
use.alpha, "\n")
198
if(exists("flush.console"))
199
flush.console()
196
flush.console()
200
197
}
201
198
fini.linesearch = TRUE
202
199
}
247
244
if(length(Dzero)) diag(M)[,-Dzero,drop=FALSE] else diag(M)},
248
245
Aoffset + (1:Qoffset))
249
246
if(p1)
250
for(k in 1:p1)
251
cmat2[[Aoffset+Qoffset+k]] <- Blist[[colx1.index[k]]]
247
for(kk in 1:p1)
248
cmat2[[Aoffset+Qoffset+kk]] <- Blist[[colx1.index[kk]]]
252
249
if(!no.thrills) {
253
250
i63 = iam(NA, NA, M=Rank, both=TRUE)
254
251
names(cmat2) = c(
274
271
asx = attr(x, "assign")
275
272
asx = vector("list", ncol(new.lv.model.matrix))
276
273
names(asx) = names(cmat2)
277
for(i in 1:length(names(asx))) {
278
asx[[i]] = i
274
for(ii in 1:length(names(asx))) {
275
asx[[ii]] = ii
279
276
}
280
277
attr(new.lv.model.matrix, "assign") = asx
281
278
}
292
289
293
290
294
291
cmat1 = replace.constraints(Blist, A, control$colx2.index)
295
fit <- vlm.wfit(x, z, Blist=cmat1, U=U, matrix.out=TRUE,
296
XBIG= FALSE, rss=TRUE, qr= FALSE, xij=control$xij)
292
fit <- vlm.wfit(xmat=x, z, Blist=cmat1, U=U, matrix.out=TRUE,
293
is.vlmX=FALSE, rss=TRUE, qr=FALSE, xij=control$xij)
297
294
C = fit$mat.coef[control$colx2.index,,drop=FALSE] %*% A %*% solve(t(A) %*% A)
298
295
299
296
list(A=A, C=C, fitted=fit$fitted, new.coeffs = fit$coef,
331
328
for(ii in 1:NOS) {
332
329
i5 = i5 + 1:MSratio
333
330
334
tmp100 = vlm.wfit(new.lv.model.matrix, zedd[,i5,drop=FALSE],
331
tmp100 = vlm.wfit(xmat=new.lv.model.matrix, zedd[,i5,drop=FALSE],
335
332
Blist=cmat2, U=U[i5,,drop=FALSE],
336
matrix.out=TRUE, XBIG=FALSE, rss=TRUE, qr=FALSE,
337
Eta.range = control$Eta.range,
333
matrix.out=TRUE, is.vlmX=FALSE, rss=TRUE,
334
qr=FALSE, Eta.range = control$Eta.range,
338
335
xij=control$xij, lp.names=lp.names[i5])
339
336
fit$rss = fit$rss + tmp100$rss
340
337
fit$mat.coef = cbind(fit$mat.coef, tmp100$mat.coef)
341
338
fit$fitted.values = cbind(fit$fitted.values, tmp100$fitted.values)
342
339
}
343
340
} else {
344
fit = vlm.wfit(new.lv.model.matrix, zedd, Blist=cmat2, U=U,
345
matrix.out=TRUE, XBIG= FALSE, rss=TRUE, qr= FALSE,
341
fit = vlm.wfit(xmat=new.lv.model.matrix, zedd, Blist=cmat2, U=U,
342
matrix.out=TRUE, is.vlmX=FALSE, rss=TRUE, qr=FALSE,
346
343
Eta.range = control$Eta.range,
347
344
xij=control$xij, lp.names=lp.names)
348
345
}
376
373
377
374
378
375
rrr.init.expression <- expression({
379
if(backchat || control$Quadratic)
380
copyxbig <- TRUE
376
if(control$Quadratic)
377
copy_X_vlm <- TRUE
381
378
382
379
modelno = switch(family@vfamily[1], "poissonff"=2,
383
380
"quasipoissonff"=2, "quasipoisson"=2,
384
381
"binomialff"=1, "quasibinomialff"=1,
385
382
"quasibinomial"=1, "negbinomial"=3,
386
383
"gamma2"=5, "gaussianff"=8,
387
0) # stop("can't fit this model using fast algorithm")
384
0) # stop("cannot fit this model using fast algorithm")
388
385
if(modelno == 1) modelno = get("modelno", envir = VGAMenv)
389
386
rrcontrol$modelno = control$modelno = modelno
390
387
if(modelno==3 || modelno==5) {
507
504
508
505
rrr.end.expression = expression({
509
506
510
if(is.R()) {
511
if(exists(".VGAM.etamat", envir = VGAMenv))
512
rm(".VGAM.etamat", envir = VGAMenv)
513
} else {
514
while(exists(".VGAM.etamat", inherits=TRUE))
515
rm(".VGAM.etamat", inherits=TRUE)
516
}
507
if(exists(".VGAM.etamat", envir = VGAMenv))
508
rm(".VGAM.etamat", envir = VGAMenv)
517
509
518
510
519
511
if(control$Quadratic) {
525
517
Blist = replace.constraints(Blist.save, Amat, colx2.index)
526
518
}
527
519
528
xbig.save = if(control$Quadratic) {
520
X_vlm_save = if(control$Quadratic) {
529
521
tmp300 = lm2qrrvlm.model.matrix(x=x, Blist=Blist.save,
530
522
C=Cmat, control=control)
531
523
lv.mat = tmp300$lv.mat # Needed at the top of new.s.call
567
559
568
560
569
561
570
which.optimizer = if(is.R()) {
571
if(control$Quadratic && control$FastAlgorithm) {
572
"BFGS"
573
} else {
574
if(iter <= rrcontrol$Switch.optimizer) "Nelder-Mead" else "BFGS"
575
}
576
} else "Quasi-Newton"
562
which.optimizer = if(control$Quadratic && control$FastAlgorithm) {
563
"BFGS"
564
} else {
565
if(iter <= rrcontrol$Switch.optimizer) "Nelder-Mead" else "BFGS"
566
}
577
567
if(trace && control$OptimizeWrtC) {
578
568
cat("\n\n")
579
569
cat("Using", which.optimizer, "\n")
580
if(exists("flush.console"))
581
flush.console()
570
flush.console()
582
571
}
583
572
584
573
constraints=replace.constraints(constraints,diag(M),rrcontrol$colx2.index)
586
575
all(trivial.constraints(constraints) == 1)
587
576
588
577
theta0 <- c(Cmat)
589
if(is.R()) assign(".VGAM.dot.counter", 0, envir = VGAMenv) else
590
.VGAM.dot.counter <<- 0
591
if(control$OptimizeWrtC) {
592
if(is.R()) {
578
assign(".VGAM.dot.counter", 0, envir = VGAMenv)
579
if(control$OptimizeWrtC) {
593
580
if(control$Quadratic && control$FastAlgorithm) {
594
581
if(iter == 2) {
595
if(is.R()) {
596
if(exists(".VGAM.etamat", envir = VGAMenv))
597
rm(".VGAM.etamat", envir = VGAMenv)
598
} else {
599
if(exists(".VGAM.etamat", inherits=TRUE))
600
rm(".VGAM.etamat", inherits=TRUE)
601
}
582
if(exists(".VGAM.etamat", envir = VGAMenv))
583
rm(".VGAM.etamat", envir = VGAMenv)
602
584
}
603
585
if(iter > 2 && !quasi.newton$convergence) {
604
if(is.R()) {
605
if(zthere <- exists(".VGAM.z", envir = VGAMenv)) {
606
..VGAM.z = get(".VGAM.z", envir = VGAMenv)
607
..VGAM.U = get(".VGAM.U", envir = VGAMenv)
608
..VGAM.beta = get(".VGAM.beta", envir = VGAMenv)
609
}
610
} else {
611
if(zthere <- exists(".VGAM.z")) {
612
..VGAM.z = .VGAM.z
613
..VGAM.U = .VGAM.U
614
..VGAM.beta = .VGAM.beta
615
}
586
if(zthere <- exists(".VGAM.z", envir = VGAMenv)) {
587
..VGAM.z = get(".VGAM.z", envir = VGAMenv)
588
..VGAM.U = get(".VGAM.U", envir = VGAMenv)
589
..VGAM.beta = get(".VGAM.beta", envir = VGAMenv)
616
590
}
617
591
if(zthere) {
618
592
z = matrix(..VGAM.z, n, M) # minus any offset
624
598
if(iter == 2 || quasi.newton$convergence) {
625
599
NOS = ifelse(modelno==3 || modelno==5, M/2, M)
626
600
627
canfitok = if(is.R())
628
(exists("CQO.FastAlgorithm", envir=VGAMenv) &&
629
get("CQO.FastAlgorithm", envir = VGAMenv)) else
630
(exists("CQO.FastAlgorithm",inherits=TRUE) && CQO.FastAlgorithm)
601
canfitok = (exists("CQO.FastAlgorithm", envir=VGAMenv) &&
602
get("CQO.FastAlgorithm", envir = VGAMenv))
631
603
if(!canfitok)
632
stop("can't fit this model using fast algorithm")
604
stop("cannot fit this model using fast algorithm")
633
605
p2star = if(nice31)
634
ifelse(control$IToleran, Rank, Rank+0.5*Rank*(Rank+1)) else
635
(NOS*Rank + Rank*(Rank+1)/2 * ifelse(control$EqualTol,1,NOS))
606
ifelse(control$IToleran, Rank, Rank+0.5*Rank*(Rank+1)) else
607
(NOS*Rank + Rank*(Rank+1)/2 * ifelse(control$EqualTol,1,NOS))
636
608
p1star = if(nice31) p1 * ifelse(modelno==3 || modelno==5,2,1) else
637
(ncol(xbig.save) - p2star)
638
xbig.save1 = if(p1star > 0) xbig.save[,-(1:p2star)] else NULL
609
(ncol(X_vlm_save) - p2star)
610
X_vlm_1save = if(p1star > 0) X_vlm_save[,-(1:p2star)] else NULL
639
611
quasi.newton = optim(par=Cmat, fn=callcqof,
640
612
gr=if(control$GradientFunction) calldcqof else NULL,
641
613
method=which.optimizer,
643
615
parscale=rep(control$Parscale, len=length(Cmat)),
644
616
maxit=250),
645
617
etamat=eta, xmat=x, ymat=y, wvec=w,
646
xbig.save1 = if(nice31) NULL else xbig.save1,
618
X_vlm_1save = if(nice31) NULL else X_vlm_1save,
647
619
modelno=modelno, Control=control,
648
620
n=n, M=M, p1star=p1star, p2star=p2star, nice31=nice31)
649
621
650
622
651
if(is.R()) {
652
if(zthere <- exists(".VGAM.z", envir = VGAMenv)) {
653
..VGAM.z = get(".VGAM.z", envir = VGAMenv)
654
..VGAM.U = get(".VGAM.U", envir = VGAMenv)
655
..VGAM.beta = get(".VGAM.beta", envir = VGAMenv)
656
}
657
} else {
658
if(zthere <- exists(".VGAM.z")) {
659
..VGAM.z = .VGAM.z
660
..VGAM.U = .VGAM.U
661
..VGAM.beta = .VGAM.beta
662
}
623
if(zthere <- exists(".VGAM.z", envir = VGAMenv)) {
624
..VGAM.z = get(".VGAM.z", envir = VGAMenv)
625
..VGAM.U = get(".VGAM.U", envir = VGAMenv)
626
..VGAM.beta = get(".VGAM.beta", envir = VGAMenv)
663
627
}
664
628
if(zthere) {
665
629
z = matrix(..VGAM.z, n, M) # minus any offset
666
630
U = matrix(..VGAM.U, M, n)
667
631
}
668
632
} else {
669
if(is.R()) {
670
if(exists(".VGAM.offset", envir = VGAMenv))
671
rm(".VGAM.offset", envir = VGAMenv)
672
} else {
673
while(exists(".VGAM.offset", inherits=TRUE))
674
rm(".VGAM.offset", inherits=TRUE)
675
}
633
if(exists(".VGAM.offset", envir = VGAMenv))
634
rm(".VGAM.offset", envir = VGAMenv)
676
635
}
677
636
} else {
678
637
use.reltol = if(length(rrcontrol$Reltol) >= iter)
689
648
M=M, xmat=x, # varbix2=varbix2,
690
649
Blist=Blist, rrcontrol=rrcontrol)
691
650
}
692
} else {
693
quasi.newton <-
694
nlminb(start=theta0,
695
objective=rrr.derivC.rss,
696
control = nlminb.control(x.tol = rrcontrol$X.tol,
697
eval.max = rrcontrol$Eval.max,
698
iter.max = rrcontrol$Iter.max,
699
abs.tol = rrcontrol$Abs.tol,
700
rel.tol = rrcontrol$Rel.tol,
701
step.min = rrcontrol$Step.min,
702
rel.err = rrcontrol$Rel.err),
703
U=U, z= if(control$ITolerances) z+offset else z,
704
M=M, xmat=x, # varbix2=varbix2,
705
Blist=Blist, rrcontrol=rrcontrol)
706
}
651
707
652
708
653
709
654
714
659
evnu = eigen(var(numat))
715
660
Cmat = Cmat %*% evnu$vector
716
661
numat = x[,rrcontrol$colx2.index,drop=FALSE] %*% Cmat
717
offset = if(Rank > 1) -0.5*apply(numat^2, 1, sum) else -0.5*numat^2
662
offset = if(Rank > 1) -0.5*rowSums(numat^2) else -0.5*numat^2
718
663
}
719
664
}
720
665
735
680
736
681
if(trace && control$OptimizeWrtC) {
737
682
cat("\n")
738
cat(which.optimizer, "using",
739
if(is.R()) "optim():" else "nlminb():", "\n")
740
cat("Objective =", if(is.R())
741
quasi.newton$value else format(quasi.newton$objective), "\n")
683
cat(which.optimizer, "using optim():\n")
684
cat("Objective =", quasi.newton$value, "\n")
742
685
cat("Parameters (= c(C)) = ", if(length(quasi.newton$par) < 5)
743
686
"" else "\n")
744
cat(if(is.R()) alt$Cmat else format(alt$Cmat), fill=TRUE)
687
cat(alt$Cmat, fill=TRUE)
745
688
cat("\n")
746
if(!is.R())
747
cat("Gradient norm =", format(quasi.newton$grad.norm), "\n")
748
cat("Number of function evaluations =", if(is.R())
749
quasi.newton$count[1] else quasi.newton$f.evals, "\n")
750
if(!is.R())
751
cat("Number of gradient evaluations =", quasi.newton$g.evals, "\n")
689
cat("Number of function evaluations =", quasi.newton$count[1], "\n")
752
690
if(length(quasi.newton$message))
753
691
cat("Message =", quasi.newton$message, "\n")
754
692
cat("\n")
755
if(exists("flush.console"))
756
flush.console()
693
flush.console()
757
694
}
758
695
759
696
771
708
772
709
if(rrcontrol$trace) {
773
710
cat(".")
774
if(exists("flush.console"))
775
flush.console()
711
flush.console()
776
712
}
777
alreadyThere = if(is.R())
778
exists(".VGAM.dot.counter", envir = VGAMenv) else
779
exists(".VGAM.dot.counter")
713
alreadyThere = exists(".VGAM.dot.counter", envir = VGAMenv)
780
714
if(alreadyThere) {
781
if(is.R()) {
782
VGAM.dot.counter = get(".VGAM.dot.counter", envir = VGAMenv)
783
VGAM.dot.counter = VGAM.dot.counter + 1
784
assign(".VGAM.dot.counter", VGAM.dot.counter, envir = VGAMenv)
785
} else {
786
.VGAM.dot.counter <<- .VGAM.dot.counter + 1
787
}
715
VGAM.dot.counter = get(".VGAM.dot.counter", envir = VGAMenv)
716
VGAM.dot.counter = VGAM.dot.counter + 1
717
assign(".VGAM.dot.counter", VGAM.dot.counter, envir = VGAMenv)
788
718
if(VGAM.dot.counter > max(50, options()$width - 5)) {
789
719
if(rrcontrol$trace) {
790
720
cat("\n")
791
if(exists("flush.console"))
792
flush.console()
721
flush.console()
793
722
}
794
if(is.R()) assign(".VGAM.dot.counter", 0, envir = VGAMenv) else
795
.VGAM.dot.counter <<- 0
723
assign(".VGAM.dot.counter", 0, envir = VGAMenv)
796
724
}
797
725
}
798
726
801
729
802
730
tmp700 = lm2qrrvlm.model.matrix(x=xmat, Blist=Blist,
803
731
no.thrills = !rrcontrol$Corner,
804
C=Cmat, control=rrcontrol, assign= FALSE)
732
C=Cmat, control=rrcontrol, assign=FALSE)
805
733
Blist = tmp700$constraints # Doesn't contain \bI_{Rank} \bnu
806
734
807
735
if(rrcontrol$Corner) {
812
740
if(length(tmp700$offset)) z = z - tmp700$offset
813
741
814
742
815
vlm.wfit(x=tmp700$new.lv.model.matrix, z=z,
743
vlm.wfit(xmat=tmp700$new.lv.model.matrix, zmat=z,
816
744
Blist=Blist, ncolx=ncol(xmat), U=U, only.rss=TRUE,
817
matrix.out= FALSE, XBIG= FALSE, rss= TRUE, qr= FALSE,
745
matrix.out=FALSE, is.vlmX=FALSE, rss= TRUE, qr=FALSE,
818
746
Eta.range = rrcontrol$Eta.range,
819
747
xij=rrcontrol$xij)$rss
820
748
}
841
769
842
770
843
771
844
if(is.R())
845
772
nlminbcontrol = function(Abs.tol = 10^(-6),
846
773
Eval.max=91,
847
774
Iter.max=91,
869
796
870
797
871
798
if(length(varlvI) != 1 || !is.logical(varlvI))
872
stop("\"varlvI\" must be TRUE or FALSE")
873
if(length(reference) > 1) stop("\"reference\" must be of length 0 or 1")
799
stop("'varlvI' must be TRUE or FALSE")
800
if(length(reference) > 1) stop("'reference' must be of length 0 or 1")
874
801
if(length(reference) && is.Numeric(reference))
875
802
if(!is.Numeric(reference, allow=1, integ=TRUE))
876
stop("bad input for argument \"reference\"")
803
stop("bad input for argument 'reference'")
877
804
if(!is.logical(ConstrainedQO <- object@control$ConstrainedQO))
878
stop("can't determine whether the model is constrained or not")
805
stop("cannot determine whether the model is constrained or not")
879
806
ocontrol = object@control
880
807
coef.object = object@coefficients
881
808
Rank = ocontrol$Rank
913
840
914
841
td.expression = expression({
915
842
Tolerance = Darray = m2adefault(Dmat, M=Rank)
916
for(i in 1:M)
917
if(length(Dzero) && any(Dzero == i)) {
918
Tolerance[,,i] = NA # Darray[,,i] == O
919
bellshaped[i] = FALSE
843
for(ii in 1:M)
844
if(length(Dzero) && any(Dzero == ii)) {
845
Tolerance[,,ii] = NA # Darray[,,ii] == O
846
bellshaped[ii] = FALSE
920
847
} else {
921
Tolerance[,,i] = -0.5 * solve(Darray[,,i])
922
bellshaped[i] = all(eigen(Tolerance[,,i])$values > 0)
848
Tolerance[,,ii] = -0.5 * solve(Darray[,,ii])
849
bellshaped[ii] = all(eigen(Tolerance[,,ii])$values > 0)
923
850
}
924
851
optimum = matrix(as.numeric(NA),Rank,M) # dimnames=list(lv.names,ynames)
925
for(i in 1:M)
926
if(bellshaped[i])
927
optimum[,i] = Tolerance[,,i] %*% cbind(Amat[i,])
852
for(ii in 1:M)
853
if(bellshaped[ii])
854
optimum[,ii] = Tolerance[,,ii] %*% cbind(Amat[ii,])
928
855
})
929
856
Amat = object@extra$Amat # M x Rank
930
857
Cmat = object@extra$Cmat # p2 x Rank
935
862
if(is.character(reference)) {
936
863
reference = (1:NOS)[reference == ynames]
937
864
if(length(reference) != 1)
938
stop("could not match argument \"reference\" with any response")
865
stop("could not match argument 'reference' with any response")
939
866
}
940
867
ptr1 = 1
941
868
candidates = if(length(reference)) reference else {
971
898
eval(td.expression)
972
899
} else {
973
900
if(length(reference) == 1)
974
stop(paste("tolerance matrix specified by \"reference\"",
975
"is not positive-definite")) else
976
warning(paste("could not find any positive-definite",
977
"tolerance matrix"))
901
stop("tolerance matrix specified by 'reference' ",
902
"is not positive-definite") else
903
warning("could not find any positive-definite ",
904
"tolerance matrix")
978
905
}
979
906
980
907
1013
940
cx1i = ocontrol$colx1.index
1014
941
maximum = if(length(cx1i)==1 && names(cx1i)=="(Intercept)") {
1015
942
eta.temp = B1
1016
for(i in 1:M)
1017
eta.temp[i] = eta.temp[i] +
1018
Amat[i,,drop=FALSE] %*% optimum[,i,drop=FALSE] +
1019
t(optimum[,i,drop=FALSE]) %*%
1020
Darray[,,i,drop= TRUE] %*% optimum[,i,drop=FALSE]
943
for(ii in 1:M)
944
eta.temp[ii] = eta.temp[ii] +
945
Amat[ii,,drop=FALSE] %*% optimum[,ii,drop=FALSE] +
946
t(optimum[,ii,drop=FALSE]) %*%
947
Darray[,,ii,drop= TRUE] %*% optimum[,ii,drop=FALSE]
1021
948
mymax = object@family@inverse(rbind(eta.temp), extra=object@extra)
1022
949
c(mymax) # Convert from matrix to vector
1023
950
} else {
1051
978
Tolerance=Tolerance)
1052
979
if(ConstrainedQO) {ans@C = Cmat} else {Cmat = NULL}
1053
980
1054
for(r in 1:Rank)
1055
ans@OptimumOrder[r,] = order(ans@Optimum[r,])
1056
for(r in 1:Rank)
1057
ans@lvOrder[,r] = order(ans@lv[,r])
981
for(rrr in 1:Rank)
982
ans@OptimumOrder[rrr,] = order(ans@Optimum[rrr,])
983
for(rrr in 1:Rank)
984
ans@lvOrder[,rrr] = order(ans@lv[,rrr])
1058
985
1059
986
if(length(object@misc$estimated.dispersion) &&
1060
987
object@misc$estimated.dispersion) {
1112
1039
"Dzero" = "logical",
1113
1040
"Tolerance" = "array"))
1114
1041
1115
setClass(Class="Coef.qrrvglm", representation("Coef.uqo",
1116
"C" = "matrix"))
1042
setClass(Class="Coef.qrrvglm", representation(
1043
"C" = "matrix"),
1044
contains = "Coef.uqo")
1117
1045
1118
1046
printCoef.qrrvglm = function(x, ...) {
1119
1047
1121
1049
Rank = object@Rank
1122
1050
M = nrow(object@A)
1123
1051
NOS = object@NOS
1124
iii = matrix(as.numeric(NA), NOS, Rank)
1052
mymat = matrix(as.numeric(NA), NOS, Rank)
1125
1053
if(Rank == 1) { # || object@Diagonal
1126
for(i in 1:NOS) {
1127
fred = if(Rank>1) diag(object@Tolerance[,,i,drop=F]) else
1128
object@Tolerance[,,i]
1054
for(ii in 1:NOS) {
1055
fred = if(Rank>1) diag(object@Tolerance[,,ii,drop=F]) else
1056
object@Tolerance[,,ii]
1129
1057
if(all(fred > 0))
1130
iii[i,] = sqrt(fred)
1058
mymat[ii,] = sqrt(fred)
1131
1059
}
1132
dimnames(iii) = list(dimnames(object@Tolerance)[[3]],
1133
if(Rank==1) "lv" else
1134
paste("Tolerance", dimnames(iii)[[2]], sep=""))
1060
dimnames(mymat) = list(dimnames(object@Tolerance)[[3]],
1061
if(Rank==1) "lv" else
1062
paste("Tolerance", dimnames(mymat)[[2]], sep=""))
1135
1063
} else {
1136
for(i in 1:NOS) {
1137
fred = eigen(object@Tolerance[,,i])
1064
for(ii in 1:NOS) {
1065
fred = eigen(object@Tolerance[,,ii])
1138
1066
if(all(fred$value > 0))
1139
iii[i,] = sqrt(fred$value)
1067
mymat[ii,] = sqrt(fred$value)
1140
1068
}
1141
dimnames(iii) = list(dimnames(object@Tolerance)[[3]],
1142
paste("tol", 1:Rank, sep=""))
1069
dimnames(mymat) = list(dimnames(object@Tolerance)[[3]],
1070
paste("tol", 1:Rank, sep=""))
1143
1071
}
1144
1072
1145
1073
dimnames(object@A) = list(dimnames(object@A)[[1]],
1146
1074
if(Rank > 1) paste("A", dimnames(object@A)[[2]], sep=".") else "A")
1147
1075
1148
1076
Maximum = if(length(object@Maximum)) cbind(Maximum=object@Maximum) else NULL
1149
if(length(Maximum) && length(iii) && Rank==1)
1150
Maximum[is.na(iii),] = NA
1077
if(length(Maximum) && length(mymat) && Rank==1)
1078
Maximum[is.na(mymat),] = NA
1151
1079
1152
1080
optmat = cbind(t(object@Optimum))
1153
1081
dimnames(optmat) = list(dimnames(optmat)[[1]],
1154
1082
if(Rank > 1) paste("Optimum", dimnames(optmat)[[2]], sep=".")
1155
1083
else "Optimum")
1156
if(length(optmat) && length(iii) && Rank==1)
1157
optmat[is.na(iii),] = NA
1084
if(length(optmat) && length(mymat) && Rank==1)
1085
optmat[is.na(mymat),] = NA
1158
1086
1159
1087
if( object@Constrained ) {
1160
1088
cat("\nC matrix (constrained/canonical coefficients)\n")
1169
1097
if(Rank > 1) { # !object@Diagonal && Rank > 1
1170
1098
cat("\nTolerances\n") } else
1171
1099
cat("\nTolerance\n")
1172
print(iii, ...)
1100
print(mymat, ...)
1173
1101
1174
1102
cat("\nStandard deviation of the latent variables (site scores)\n")
1175
1103
print(sd(object@lv))
1177
1105
}
1178
1106
1179
1107
1180
setMethod("show", "Coef.qrrvglm", function(object)
1181
printCoef.qrrvglm(object))
1182
setMethod("print", "Coef.qrrvglm", function(x, ...)
1183
printCoef.qrrvglm(x, ...))
1184
setMethod("summary", "qrrvglm", function(object, ...)
1185
summary.qrrvglm(object, ...))
1108
setMethod("show", "Coef.qrrvglm", function(object)
1109
printCoef.qrrvglm(object))
1110
setMethod("print", "Coef.qrrvglm", function(x, ...)
1111
printCoef.qrrvglm(x, ...))
1112
setMethod("summary", "qrrvglm", function(object, ...)
1113
summary.qrrvglm(object, ...))
1186
1114
1187
1115
1188
1116
predict.qrrvglm <- function(object,
1195
1123
varlvI = FALSE, reference = NULL, ...)
1196
1124
{
1197
1125
if(se.fit)
1198
stop("can't handle se.fit==TRUE yet")
1126
stop("cannot handle se.fit==TRUE yet")
1199
1127
if(deriv != 0)
1200
1128
stop("derivative is not equal to 0")
1201
1129
1203
1131
type <- as.character(substitute(type))
1204
1132
type <- match.arg(type, c("link", "response", "lv", "terms"))[1]
1205
1133
if(type=="lv")
1206
stop("can't handle type='lv' yet")
1134
stop("cannot handle type='lv' yet")
1207
1135
if(type=="terms")
1208
stop("can't handle type='terms' yet")
1136
stop("cannot handle type='terms' yet")
1209
1137
1210
1138
M = object@misc$M
1211
1139
Rank = object@control$Rank
1221
1149
}
1222
1150
}
1223
1151
1224
attrassignlm <- function(object, ...)
1225
attrassigndefault(model.matrix(object), object@terms)
1226
1227
attrassigndefault <- function(mmat, tt) {
1228
if (!inherits(tt, "terms"))
1229
stop("need terms object")
1230
aa <- attr(mmat, "assign")
1231
if (is.null(aa))
1232
stop("argument is not really a model matrix")
1233
ll <- attr(tt, "term.labels")
1234
if (attr(tt, "intercept") > 0)
1235
ll <- c("(Intercept)", ll)
1236
aaa <- factor(aa, labels = ll)
1237
split(order(aa), aaa)
1238
}
1239
1240
1152
if(!length(newdata)) {
1241
1153
X <- model.matrixvlm(object, type="lm", ...)
1242
1154
offset <- object@offset
1243
1155
tt <- object@terms$terms # terms(object)
1244
if(is.R() && !length(object@x))
1245
attr(X, "assign") <- attrassignlm(X, tt)
1156
if(!length(object@x))
1157
attr(X, "assign") = attrassignlm(X, tt)
1246
1158
} else {
1247
1159
if(is.smart(object) && length(object@smart.prediction)) {
1248
1160
setup.smart("read", smart.prediction=object@smart.prediction)
1253
1165
if(length(object@contrasts)) object@contrasts else NULL,
1254
1166
xlev = object@xlevels)
1255
1167
1256
if(is.R() && nrow(X)!=nrow(newdata)) {
1168
if(nrow(X) != nrow(newdata)) {
1257
1169
as.save = attr(X, "assign")
1258
1170
X = X[rep(1, nrow(newdata)),,drop=FALSE]
1259
1171
dimnames(X) = list(dimnames(newdata)[[1]], "(Intercept)")
1271
1183
wrapup.smart()
1272
1184
}
1273
1185
1274
if(is.R())
1275
attr(X, "assign") <- attrassigndefault(X, tt)
1186
attr(X, "assign") = attrassigndefault(X, tt)
1276
1187
}
1277
1188
1278
1189
ocontrol = object@control
1335
1246
predict.qrrvglm(object, ...))
1336
1247
1337
1248
coefqrrvglm = function(object, matrix.out = FALSE,
1338
label = TRUE, compress = TRUE) {
1249
label = TRUE) {
1339
1250
if(matrix.out)
1340
stop("currently can't handle matrix.out=TRUE")
1341
coefvlm(object, matrix.out = matrix.out, label = label, compress = compress)
1251
stop("currently cannot handle matrix.out=TRUE")
1252
coefvlm(object, matrix.out = matrix.out, label = label)
1342
1253
}
1343
1254
1344
1255
1359
1270
printrrvglm <- function(x, ...)
1360
1271
{
1361
1272
if(!is.null(cl <- x@call)) {
1362
cat("Call:\n")
1363
dput(cl)
1273
cat("Call:\n")
1274
dput(cl)
1364
1275
}
1365
coef <- x@coefficients
1366
if(any(nas <- is.na(coef))) {
1367
if(is.null(names(coef)))
1368
names(coef) <- paste("b", 1:length(coef), sep = "")
1369
cat("\nCoefficients: (", sum(nas),
1370
" not defined because of singularities)\n", sep = "")
1276
vecOfBetas <- x@coefficients
1277
if(any(nas <- is.na(vecOfBetas))) {
1278
if(is.null(names(vecOfBetas)))
1279
names(vecOfBetas) = paste("b", 1:length(vecOfBetas), sep="")
1280
cat("\nCoefficients: (", sum(nas),
1281
" not defined because of singularities)\n", sep = "")
1371
1282
} else
1372
1283
cat("\nCoefficients:\n")
1284
print.default(vecOfBetas, ...) # used to be print()
1373
1285
1374
1286
if(FALSE) {
1375
Rank <- x@Rank
1376
if(!length(Rank))
1377
Rank <- sum(!nas)
1287
Rank <- x@Rank
1288
if(!length(Rank))
1289
Rank <- sum(!nas)
1378
1290
}
1379
1291
1380
1292
if(FALSE) {
1392
1304
1393
1305
if(length(x@criterion)) {
1394
1306
ncrit <- names(x@criterion)
1395
for(i in ncrit)
1396
if(i!="loglikelihood" && i!="deviance")
1397
cat(paste(i, ":", sep=""), format(x@criterion[[i]]), "\n")
1307
for(iii in ncrit)
1308
if(iii != "loglikelihood" && iii != "deviance")
1309
cat(paste(iii, ":", sep=""), format(x@criterion[[iii]]), "\n")
1398
1310
}
1399
1311
1400
1312
invisible(x)
1405
1317
1406
1318
setMethod("print", "rrvglm", function(x, ...) printrrvglm(x, ...))
1407
1319
1408
setMethod("show", "rrvglm", function(object) printrrvglm(object))
1409
1410
1411
1412
1413
rrvglm.control.Gaussian <- function(backchat= FALSE, half.stepsizing= FALSE,
1320
setMethod("show", "rrvglm", function(object) printrrvglm(object))
1321
1322
1323
1324
1325
rrvglm.control.Gaussian <- function(half.stepsizing= FALSE,
1414
1326
save.weight= TRUE, ...)
1415
1327
{
1416
1328
1417
list(backchat= FALSE, half.stepsizing= FALSE,
1329
list(half.stepsizing= FALSE,
1418
1330
save.weight=as.logical(save.weight)[1])
1419
1331
}
1420
1332
1421
1333
1422
1334
1423
summary.rrvglm <- function(object, correlation= FALSE,
1335
summary.rrvglm <- function(object, correlation=FALSE,
1424
1336
dispersion=NULL, digits=NULL,
1425
1337
numerical= TRUE,
1426
1338
h.step = 0.0001,
1427
kill.all= FALSE, omit13= FALSE, fixA= FALSE, ...)
1339
kill.all=FALSE, omit13=FALSE, fixA=FALSE, ...)
1428
1340
{
1429
1341
1430
1342
1432
1344
1433
1345
1434
1346
if(!is.Numeric(h.step, allow=1) || abs(h.step)>1)
1435
stop("bad input for \"h.step\"")
1347
stop("bad input for 'h.step'")
1436
1348
1437
1349
if(!object@control$Corner)
1438
1350
stop("this function works with corner constraints only")
1440
1352
if(is.null(dispersion))
1441
1353
dispersion <- object@misc$dispersion
1442
1354
1443
newobject <- object
1444
class(newobject) <- "vglm" # 6/2/02; For Splus6
1355
newobject <- as(object, "vglm")
1356
1357
1445
1358
stuff <- summaryvglm(newobject, correlation=correlation,
1446
1359
dispersion=dispersion)
1447
1360
1515
1428
1516
1429
1517
1430
1518
get.rrvglm.se1 <- function(fit, omit13= FALSE, kill.all= FALSE,
1519
numerical= TRUE,
1520
fixA= FALSE, h.step=0.0001,
1521
trace.arg= FALSE, ...) {
1431
get.rrvglm.se1 = function(fit, omit13=FALSE, kill.all=FALSE,
1432
numerical=TRUE,
1433
fixA=FALSE, h.step=0.0001,
1434
trace.arg=FALSE, ...) {
1522
1435
1523
1436
1524
1437
1525
1438
1526
1439
if(length(fit@control$Nested) && fit@control$Nested)
1527
stop("sorry, can't handle nested models yet")
1440
stop("sorry, cannot handle nested models yet")
1528
1441
1529
1442
Structural.zero = fit@control$Structural.zero
1530
1443
1552
1465
1553
1466
wz <- weights(fit, type="w") # old: wweights(fit) #fit@weights
1554
1467
if(!length(wz))
1555
stop("can't get fit@weights")
1468
stop("cannot get fit@weights")
1556
1469
1557
1470
M <- fit@misc$M
1558
1471
n <- fit@misc$n
1580
1493
}
1581
1494
1582
1495
1583
newobject <- fit
1584
class(newobject) <- "vglm" # 6/2/02; For Splus6
1496
1497
1498
newobject <- as(fit, "vglm")
1499
1500
1501
1502
1585
1503
sfit2233 <- summaryvglm(newobject)
1586
1504
d8 <- dimnames(sfit2233@cov.unscaled)[[1]]
1587
1505
cov2233 <- solve(sfit2233@cov.unscaled) # Includes any intercepts
1591
1509
nassign = names(fit@constraints)
1592
1510
choose.from = varassign(fit@constraints, nassign)
1593
1511
for(ii in nassign)
1594
if(any(ii== names(colx2.index))) {
1512
if(any(ii == names(colx2.index))) {
1595
1513
log.vec33 = c(log.vec33, choose.from[[ii]])
1596
1514
}
1597
1515
cov33 = cov2233[ log.vec33, log.vec33, drop=FALSE] # r*p2 by r*p2
1603
1521
offs = matrix(0, n, M) # The "0" handles Structural.zero's
1604
1522
offs[,Index.corner] = lv.mat
1605
1523
if(M == (Rank+length(Structural.zero)))
1606
stop("can't handle full-rank models yet")
1524
stop("cannot handle full-rank models yet")
1607
1525
cm = matrix(0, M, M-Rank-length(Structural.zero))
1608
1526
cm[-c(Index.corner,Structural.zero),] = diag(M-Rank-length(Structural.zero))
1609
1527
1610
1528
Blist = vector("list", length(colx1.index)+1)
1611
1529
names(Blist) = c(names(colx1.index), "I(lv.mat)")
1612
1530
for(ii in names(colx1.index))
1613
Blist[[ii]] = fit@constraints[[ii]]
1531
Blist[[ii]] = fit@constraints[[ii]]
1614
1532
Blist[["I(lv.mat)"]] = cm
1615
1533
1616
1534
1637
1555
}
1638
1556
1639
1557
1640
if(( is.R() && fit@misc$dataname == "list") ||
1641
(!is.R() && fit@misc$dataname == "sys.parent")) {
1558
if(fit@misc$dataname == "list") {
1642
1559
dspec = FALSE
1643
1560
} else {
1644
if(is.R()) {
1645
mytext1 = "exists(x=fit@misc$dataname, envir = VGAMenv)"
1646
myexp1 = parse(text=mytext1)
1647
is.there = eval(myexp1)
1648
bbdata= if(is.there) get(fit@misc$dataname, envir=VGAMenv) else
1649
get(fit@misc$dataname)
1650
} else {
1651
bbdata = get(fit@misc$dataname)
1652
}
1561
mytext1 = "exists(x=fit@misc$dataname, envir = VGAMenv)"
1562
myexp1 = parse(text=mytext1)
1563
is.there = eval(myexp1)
1564
bbdata= if(is.there) get(fit@misc$dataname, envir=VGAMenv) else
1565
get(fit@misc$dataname)
1653
1566
dspec = TRUE
1654
1567
}
1655
1568
1656
if(!is.R()) {
1657
1658
stop("26-9-2007: uncomment out the following lines to run it in Splus")
1659
}
1660
1569
1661
1570
fit1122 <- if(dspec) vlm(bb,
1662
1571
constraint=Blist, crit="d", weight=wz, data=bbdata,
1663
save.weight= TRUE, smart= FALSE, trace=trace.arg, x= TRUE) else
1572
save.weight=TRUE, smart=FALSE, trace=trace.arg, x=TRUE) else
1664
1573
vlm(bb,
1665
1574
constraint=Blist, crit="d", weight=wz,
1666
save.weight= TRUE, smart= FALSE, trace=trace.arg, x= TRUE)
1575
save.weight=TRUE, smart=FALSE, trace=trace.arg, x=TRUE)
1667
1576
1668
1577
1669
1578
1737
1646
dct.da <- matrix(as.numeric(NA), (M-r-length(Structural.zero))*r, r*p2)
1738
1647
1739
1648
if((length(Index.corner) + length(Structural.zero)) == M)
1740
stop("can't handle full rank models yet")
1649
stop("cannot handle full rank models yet")
1741
1650
cbindex = (1:M)[-c(Index.corner, Structural.zero)]
1742
1651
1743
1652
ptr = 1
1744
for(s in 1:r)
1653
for(sss in 1:r)
1745
1654
for(tt in cbindex) {
1746
1655
small.Blist = vector("list", p2)
1747
1656
pAmat = Aimat
1748
pAmat[tt,s] = pAmat[tt,s] + h.step # Perturb it
1657
pAmat[tt,sss] = pAmat[tt,sss] + h.step # Perturb it
1749
1658
for(ii in 1:p2)
1750
1659
small.Blist[[ii]] = pAmat
1751
1660
1760
1669
1761
1670
fred = weights(fit, type="w", deriv= TRUE, ignore.slot= TRUE)
1762
1671
if(!length(fred))
1763
stop("can't get @weights and $deriv from object")
1672
stop("cannot get @weights and $deriv from object")
1764
1673
wz = fred$weights
1765
1674
deriv.mu <- fred$deriv
1766
1675
1768
1677
tvfor <- vforsub(U, as.matrix(deriv.mu), M=M, n=nn)
1769
1678
newzmat <- neweta + vbacksub(U, tvfor, M=M, n=nn) - offset
1770
1679
1771
newfit = vlm.wfit(x=x2mat, z=newzmat - x1mat %*% Bmat,
1772
Blist=small.Blist, U = U,
1773
matrix.out = FALSE, XBIG = FALSE,
1774
rss = TRUE, qr = FALSE, x.ret = FALSE, offset = NULL,
1775
xij=xij)
1680
newfit = vlm.wfit(xmat=x2mat, zmat=newzmat - x1mat %*% Bmat,
1681
Blist=small.Blist, U = U,
1682
matrix.out = FALSE, is.vlmX = FALSE,
1683
rss = TRUE, qr = FALSE, x.ret = FALSE,
1684
offset = NULL, xij=xij)
1776
1685
dct.da[ptr,] <- (newfit$coef - t(Cimat)) / h.step
1777
1686
ptr = ptr + 1
1778
1687
}
1783
1692
1784
1693
1785
1694
1786
dctda.fast.only <- function(theta, wz, U, zmat, M, r, x1mat, x2mat,
1787
p2, Index.corner, Aimat, Bmat, Cimat,
1788
xij=NULL,
1789
Structural.zero=NULL)
1695
dctda.fast.only = function(theta, wz, U, zmat, M, r, x1mat, x2mat,
1696
p2, Index.corner, Aimat, Bmat, Cimat,
1697
xij=NULL,
1698
Structural.zero=NULL)
1790
1699
{
1791
1700
1792
1701
1793
1702
if(length(Structural.zero))
1794
stop("can't handle Structural.zero in dctda.fast.only()")
1703
stop("cannot handle Structural.zero in dctda.fast.only()")
1795
1704
1796
nn <- nrow(x2mat)
1797
if(nrow(Cimat)!=p2 || ncol(Cimat)!=r)
1705
nn = nrow(x2mat)
1706
if(nrow(Cimat) != p2 || ncol(Cimat) != r)
1798
1707
stop("Cimat wrong shape")
1799
1708
1800
fred <- kronecker(matrix(1,1,r), x2mat)
1801
fred <- kronecker(fred, matrix(1,M,1))
1802
barney <- kronecker(Aimat, matrix(1,1,p2))
1803
barney <- kronecker(matrix(1,nn,1), barney)
1804
1805
temp <- array(t(barney*fred), c(p2*r, M, nn))
1806
temp <- aperm(temp, c(2,1,3)) # M by p2*r by nn
1807
temp <- mux5(wz, temp, M=M, matrix.arg= TRUE)
1808
temp <- m2adefault(temp, M=p2*r) # Note M != M here!
1809
G <- solve(apply(temp,1:2,sum)) # p2*r by p2*r
1810
1811
dc.da <- array(NA, c(p2, r, M, r)) # different from other functions
1709
fred = kronecker(matrix(1,1,r), x2mat)
1710
fred = kronecker(fred, matrix(1,M,1))
1711
barney = kronecker(Aimat, matrix(1,1,p2))
1712
barney = kronecker(matrix(1,nn,1), barney)
1713
1714
temp = array(t(barney*fred), c(p2*r, M, nn))
1715
temp = aperm(temp, c(2,1,3)) # M by p2*r by nn
1716
temp = mux5(wz, temp, M=M, matrix.arg= TRUE)
1717
temp = m2adefault(temp, M=p2*r) # Note M != M here!
1718
G = solve(rowSums(temp, dims=2)) # p2*r by p2*r
1719
1720
dc.da = array(NA, c(p2, r, M, r)) # different from other functions
1812
1721
if(length(Index.corner) == M)
1813
stop("can't handle full rank models yet")
1814
cbindex <- (1:M)[-Index.corner] # complement of Index.corner
1815
resid2 <- if(length(x1mat))
1816
mux22(t(wz), zmat - x1mat %*% Bmat, M=M, upper= FALSE, as.mat= TRUE) else
1817
mux22(t(wz), zmat , M=M, upper= FALSE, as.mat= TRUE)
1722
stop("cannot handle full rank models yet")
1723
cbindex = (1:M)[-Index.corner] # complement of Index.corner
1724
resid2 = if(length(x1mat))
1725
mux22(t(wz), zmat - x1mat %*% Bmat, M=M, upper=FALSE, as.mat=TRUE) else
1726
mux22(t(wz), zmat , M=M, upper=FALSE, as.mat=TRUE)
1818
1727
1819
for(s in 1:r)
1820
for(tt in cbindex) {
1821
fred <- t(x2mat) * matrix(resid2[,tt], p2, nn, byrow= TRUE) # p2 * nn
1822
temp2 <- kronecker(ei(s,r), apply(fred,1,sum))
1823
for(k in 1:r) {
1824
Wiak <- mux22(t(wz), matrix(Aimat[,k], nn, M, byrow= TRUE),
1825
M=M, upper= FALSE, as.mat= TRUE) # nn * M
1826
wxx <- Wiak[,tt] * x2mat
1827
blocki <- t(x2mat) %*% wxx
1828
temp4a <- blocki %*% Cimat[,k]
1829
if(k==1) {
1830
temp4b <- blocki %*% Cimat[,s]
1728
for(sss in 1:r)
1729
for(ttt in cbindex) {
1730
fred = t(x2mat) * matrix(resid2[,ttt], p2, nn, byrow= TRUE) # p2 * nn
1731
temp2 = kronecker(ei(sss,r), rowSums(fred))
1732
for(kkk in 1:r) {
1733
Wiak = mux22(t(wz), matrix(Aimat[,kkk], nn, M, byrow= TRUE),
1734
M=M, upper= FALSE, as.mat= TRUE) # nn * M
1735
wxx = Wiak[,ttt] * x2mat
1736
blocki = t(x2mat) %*% wxx
1737
temp4a = blocki %*% Cimat[,kkk]
1738
if(kkk==1) {
1739
temp4b = blocki %*% Cimat[,sss]
1831
1740
}
1832
temp2 = temp2 - kronecker(ei(s,r), temp4a) -
1833
kronecker(ei(k,r), temp4b)
1741
temp2 = temp2 - kronecker(ei(sss,r), temp4a) -
1742
kronecker(ei(kkk,r), temp4b)
1834
1743
}
1835
dc.da[,,tt,s] <- G %*% temp2
1744
dc.da[,,ttt,sss] = G %*% temp2
1836
1745
}
1837
ans1 <- dc.da[,,cbindex,,drop=FALSE] # p2 x r x (M-r) x r
1838
ans1 <- aperm(ans1, c(2,1,3,4)) # r x p2 x (M-r) x r
1746
ans1 = dc.da[,,cbindex,,drop=FALSE] # p2 x r x (M-r) x r
1747
ans1 = aperm(ans1, c(2,1,3,4)) # r x p2 x (M-r) x r
1839
1748
1840
ans1 <- matrix(c(ans1), r*p2, (M-r)*r)
1841
ans1 <- t(ans1)
1749
ans1 = matrix(c(ans1), r*p2, (M-r)*r)
1750
ans1 = t(ans1)
1842
1751
ans1
1843
1752
}
1844
1753
1845
1754
1846
1755
1847
dcda.fast <- function(theta, wz, U, z, M, r, xmat, pp, Index.corner,
1848
intercept= TRUE, xij=NULL)
1756
dcda.fast = function(theta, wz, U, z, M, r, xmat, pp, Index.corner,
1757
intercept= TRUE, xij=NULL)
1849
1758
{
1850
1759
1851
1760
1859
1768
if(intercept) {
1860
1769
Blist <- vector("list", pp+1)
1861
1770
Blist[[1]] <- diag(M)
1862
for(i in 2:(pp+1))
1863
Blist[[i]] <- Aimat
1771
for(ii in 2:(pp+1))
1772
Blist[[ii]] = Aimat
1864
1773
} else {
1865
1774
Blist <- vector("list", pp)
1866
for(i in 1:(pp))
1867
Blist[[i]] <- Aimat
1775
for(ii in 1:pp)
1776
Blist[[ii]] = Aimat
1868
1777
}
1869
1778
1870
coeffs <- vlm.wfit(xmat, z, Blist, U=U, matrix.out= TRUE,
1871
xij=xij)$mat.coef
1779
coeffs = vlm.wfit(xmat=xmat, z, Blist, U=U, matrix.out=TRUE,
1780
xij=xij)$mat.coef
1872
1781
c3 <- coeffs <- t(coeffs) # transpose to make M x (pp+1)
1873
1782
1874
1783
1887
1796
temp <- aperm(temp, c(2,1,3))
1888
1797
temp <- mux5(wz, temp, M=M, matrix.arg= TRUE)
1889
1798
temp <- m2adefault(temp, M=r*pp) # Note M != M here!
1890
G <- solve(apply(temp,1:2,sum))
1799
G = solve(rowSums(temp, dims=2))
1891
1800
1892
1801
dc.da <- array(NA, c(pp,r,M,r)) # different from other functions
1893
1802
cbindex <- (1:M)[-Index.corner]
1898
1807
for(tt in cbindex) {
1899
1808
fred <- (if(intercept) t(xmat[,-1,drop=FALSE]) else
1900
1809
t(xmat)) * matrix(resid2[,tt],pp,nn,byrow= TRUE)
1901
temp2 <- kronecker(ei(s,r), apply(fred,1,sum))
1810
temp2 <- kronecker(ei(s,r), rowSums(fred))
1902
1811
1903
1812
temp4 <- rep(0,pp)
1904
1813
for(k in 1:r) {
1928
1837
etastar <- (if(intercept) xmat[,-1,drop=FALSE] else xmat) %*% Cimat
1929
1838
eta <- matrix(int.vec, nn, M, byrow= TRUE) + etastar %*% t(Aimat)
1930
1839
1931
sumWinv <- solve((m2adefault(t(apply(wz, 2, sum)), M=M))[,,1])
1840
sumWinv <- solve((m2adefault(t(colSums(wz)), M=M))[,,1])
1932
1841
1933
1842
deta0.da <- array(0,c(M,M,r))
1934
1843
AtWi <- kronecker(matrix(1,nn,1), Aimat)
1935
1844
AtWi <- mux111(t(wz), AtWi, M=M, upper= FALSE) # matrix.arg= TRUE,
1936
1845
AtWi <- array(t(AtWi), c(r,M,nn))
1937
1846
for(ss in 1:r) {
1938
temp90 <- (m2adefault(t(apply(etastar[,ss]*wz,2,sum)), M=M))[,,1] #MxM
1847
temp90 <- (m2adefault(t(colSums(etastar[,ss]*wz)), M=M))[,,1] #MxM
1939
1848
temp92 <- array(detastar.da[,,ss,], c(M,r,nn))
1940
1849
temp93 <- mux7(temp92, AtWi)
1941
temp91 <- apply(temp93, 1:2, sum) # M x M
1850
temp91 = rowSums(temp93, dims=2) # M x M
1942
1851
deta0.da[,,ss] <- -(temp90 + temp91) %*% sumWinv
1943
1852
}
1944
1853
ans2 <- deta0.da[-(1:r),,,drop=FALSE] # (M-r) x M x r
1950
1859
1951
1860
1952
1861
1953
rrr.deriv.rss <- function(theta, wz, U, z, M, r, xmat,
1954
pp, Index.corner, intercept= TRUE,
1955
xij=NULL)
1862
rrr.deriv.rss = function(theta, wz, U, z, M, r, xmat,
1863
pp, Index.corner, intercept= TRUE,
1864
xij=NULL)
1956
1865
{
1957
1866
1958
Amat <- matrix(as.numeric(NA), M, r)
1959
Amat[Index.corner,] <- diag(r)
1960
Amat[-Index.corner,] <- theta # [-(1:M)]
1867
Amat = matrix(as.numeric(NA), M, r)
1868
Amat[Index.corner,] = diag(r)
1869
Amat[-Index.corner,] = theta # [-(1:M)]
1961
1870
1962
1871
if(intercept) {
1963
Blist <- vector("list", pp+1)
1964
Blist[[1]] <- diag(M)
1965
for(i in 2:(pp+1))
1966
Blist[[i]] <- Amat
1872
Blist = vector("list", pp+1)
1873
Blist[[1]] = diag(M)
1874
for(ii in 2:(pp+1))
1875
Blist[[ii]] = Amat
1967
1876
} else {
1968
Blist <- vector("list", pp)
1969
for(i in 1:(pp))
1970
Blist[[i]] <- Amat
1877
Blist = vector("list", pp)
1878
for(ii in 1:pp)
1879
Blist[[ii]] = Amat
1971
1880
}
1972
1881
1973
vlm.wfit(xmat, z, Blist, U=U, matrix.out= FALSE, rss= TRUE, xij=xij)$rss
1882
vlm.wfit(xmat=xmat, z, Blist, U=U, matrix.out=FALSE,
1883
rss=TRUE, xij=xij)$rss
1974
1884
}
1975
1885
1976
1886
1977
1887
1978
1888
1979
rrr.deriv.gradient.fast <- function(theta, wz, U, z, M, r, xmat,
1980
pp, Index.corner, intercept= TRUE)
1889
rrr.deriv.gradient.fast = function(theta, wz, U, z, M, r, xmat,
1890
pp, Index.corner, intercept= TRUE)
1981
1891
{
1982
1892
1983
1893
1984
1894
1985
1895
1986
nn <- nrow(xmat)
1896
nn = nrow(xmat)
1987
1897
1988
Aimat <- matrix(as.numeric(NA), M, r)
1989
Aimat[Index.corner,] <- diag(r)
1990
Aimat[-Index.corner,] <- theta # [-(1:M)]
1898
Aimat = matrix(as.numeric(NA), M, r)
1899
Aimat[Index.corner,] = diag(r)
1900
Aimat[-Index.corner,] = theta # [-(1:M)]
1991
1901
1992
1902
if(intercept) {
1993
Blist <- vector("list", pp+1)
1994
Blist[[1]] <- diag(M)
1903
Blist = vector("list", pp+1)
1904
Blist[[1]] = diag(M)
1995
1905
for(i in 2:(pp+1))
1996
Blist[[i]] <- Aimat
1906
Blist[[i]] = Aimat
1997
1907
} else {
1998
Blist <- vector("list", pp)
1908
Blist = vector("list", pp)
1999
1909
for(i in 1:(pp))
2000
Blist[[i]] <- Aimat
1910
Blist[[i]] = Aimat
2001
1911
}
2002
1912
2003
coeffs <- vlm.wfit(xmat, z, Blist, U=U, matrix.out= TRUE,
1913
coeffs = vlm.wfit(xmat, z, Blist, U=U, matrix.out= TRUE,
2004
1914
xij=NULL)$mat.coef
2005
c3 <- coeffs <- t(coeffs) # transpose to make M x (pp+1)
2006
2007
2008
int.vec <- if(intercept) c3[,1] else 0 # \boldeta_0
2009
Cimat <- if(intercept) t(c3[Index.corner,-1,drop=FALSE]) else
1915
c3 = coeffs = t(coeffs) # transpose to make M x (pp+1)
1916
1917
1918
int.vec = if(intercept) c3[,1] else 0 # \boldeta_0
1919
Cimat = if(intercept) t(c3[Index.corner,-1,drop=FALSE]) else
2010
1920
t(c3[Index.corner,,drop=FALSE])
2011
1921
if(nrow(Cimat)!=pp || ncol(Cimat)!=r)
2012
1922
stop("Cimat wrong shape")
2013
1923
2014
1924
fred = kronecker(matrix(1,1,r), if(intercept) xmat[,-1,drop=FALSE] else xmat)
2015
fred <- kronecker(fred, matrix(1,M,1))
2016
barney <- kronecker(Aimat, matrix(1,1,pp))
2017
barney <- kronecker(matrix(1,nn,1), barney)
2018
2019
temp <- array(t(barney*fred), c(r*pp,M,nn))
2020
temp <- aperm(temp, c(2,1,3))
2021
temp <- mux5(wz, temp, M=M, matrix.arg= TRUE)
2022
temp <- m2adefault(temp, M=r*pp) # Note M != M here!
2023
G <- solve(apply(temp,1:2,sum))
2024
2025
dc.da <- array(NA,c(pp,r,r,M))
2026
cbindex <- (1:M)[-Index.corner]
2027
resid2 <- mux22(t(wz), z - matrix(int.vec,nn,M,byrow= TRUE), M=M,
1925
fred = kronecker(fred, matrix(1,M,1))
1926
barney = kronecker(Aimat, matrix(1,1,pp))
1927
barney = kronecker(matrix(1,nn,1), barney)
1928
1929
temp = array(t(barney*fred), c(r*pp,M,nn))
1930
temp = aperm(temp, c(2,1,3))
1931
temp = mux5(wz, temp, M=M, matrix.arg= TRUE)
1932
temp = m2adefault(temp, M=r*pp) # Note M != M here!
1933
G = solve(rowSums(temp, dims=2))
1934
1935
dc.da = array(NA,c(pp,r,r,M))
1936
cbindex = (1:M)[-Index.corner]
1937
resid2 = mux22(t(wz), z - matrix(int.vec,nn,M,byrow= TRUE), M=M,
2028
1938
upper= FALSE, as.mat= TRUE) # mat= TRUE,
2029
1939
2030
1940
for(s in 1:r)
2031
1941
for(tt in cbindex) {
2032
fred <- (if(intercept) t(xmat[,-1,drop=FALSE]) else
1942
fred = (if(intercept) t(xmat[,-1,drop=FALSE]) else
2033
1943
t(xmat)) * matrix(resid2[,tt],pp,nn,byrow= TRUE)
2034
temp2 <- kronecker(ei(s,r), apply(fred,1,sum))
1944
temp2 = kronecker(ei(s,r), rowSums(fred))
2035
1945
2036
temp4 <- rep(0,pp)
1946
temp4 = rep(0,pp)
2037
1947
for(k in 1:r) {
2038
Wiak <- mux22(t(wz), matrix(Aimat[,k],nn,M,byrow= TRUE),
1948
Wiak = mux22(t(wz), matrix(Aimat[,k],nn,M,byrow= TRUE),
2039
1949
M=M, upper= FALSE, as.mat= TRUE) # mat= TRUE,
2040
wxx <- Wiak[,tt] * (if(intercept) xmat[,-1,drop=FALSE] else xmat)
2041
blocki <- (if(intercept) t(xmat[,-1,drop=FALSE]) else t(xmat)) %*% wxx
2042
temp4 <- temp4 + blocki %*% Cimat[,k]
1950
wxx = Wiak[,tt] * (if(intercept) xmat[,-1,drop=FALSE] else xmat)
1951
blocki = (if(intercept) t(xmat[,-1,drop=FALSE]) else t(xmat)) %*% wxx
1952
temp4 = temp4 + blocki %*% Cimat[,k]
2043
1953
}
2044
dc.da[,,s,tt] <- G %*% (temp2 - 2 * kronecker(ei(s,r),temp4))
1954
dc.da[,,s,tt] = G %*% (temp2 - 2 * kronecker(ei(s,r),temp4))
2045
1955
}
2046
1956
2047
detastar.da <- array(0,c(M,r,r,nn))
1957
detastar.da = array(0,c(M,r,r,nn))
2048
1958
for(s in 1:r)
2049
1959
for(j in 1:r) {
2050
t1 <- t(dc.da[,j,s,])
2051
t1 <- matrix(t1, M, pp)
2052
detastar.da[,j,s,] <- t1 %*% (if(intercept)
1960
t1 = t(dc.da[,j,s,])
1961
t1 = matrix(t1, M, pp)
1962
detastar.da[,j,s,] = t1 %*% (if(intercept)
2053
1963
t(xmat[,-1,drop=FALSE]) else t(xmat))
2054
1964
}
2055
1965
2056
etastar <- (if(intercept) xmat[,-1,drop=FALSE] else xmat) %*% Cimat
2057
eta <- matrix(int.vec, nn, M, byrow= TRUE) + etastar %*% t(Aimat)
2058
2059
sumWinv <- solve((m2adefault(t(apply(wz, 2, sum)), M=M))[,,1])
2060
2061
deta0.da <- array(0,c(M,M,r))
2062
2063
AtWi <- kronecker(matrix(1,nn,1), Aimat)
2064
AtWi <- mux111(t(wz), AtWi, M=M, upper= FALSE) # matrix.arg= TRUE,
2065
AtWi <- array(t(AtWi), c(r,M,nn))
1966
etastar = (if(intercept) xmat[,-1,drop=FALSE] else xmat) %*% Cimat
1967
eta = matrix(int.vec, nn, M, byrow= TRUE) + etastar %*% t(Aimat)
1968
1969
sumWinv = solve((m2adefault(t(colSums(wz)), M=M))[,,1])
1970
1971
deta0.da = array(0,c(M,M,r))
1972
1973
AtWi = kronecker(matrix(1,nn,1), Aimat)
1974
AtWi = mux111(t(wz), AtWi, M=M, upper= FALSE) # matrix.arg= TRUE,
1975
AtWi = array(t(AtWi), c(r,M,nn))
2066
1976
2067
1977
for(ss in 1:r) {
2068
temp90 <- (m2adefault(t(apply(etastar[,ss]*wz,2,sum)), M=M))[,,1] # M x M
2069
temp92 <- array(detastar.da[,,ss,],c(M,r,nn))
2070
temp93 <- mux7(temp92,AtWi)
2071
temp91 <- apply(temp93,1:2,sum) # M x M
2072
deta0.da[,,ss] <- -(temp90 + temp91) %*% sumWinv
1978
temp90 = (m2adefault(t(colSums(etastar[,ss]*wz)), M=M))[,,1] # M x M
1979
temp92 = array(detastar.da[,,ss,],c(M,r,nn))
1980
temp93 = mux7(temp92,AtWi)
1981
temp91 = rowSums(temp93, dims=2) # M x M
1982
deta0.da[,,ss] = -(temp90 + temp91) %*% sumWinv
2073
1983
}
2074
1984
2075
ans <- matrix(0,M,r)
2076
fred <- mux22(t(wz), z-eta, M=M, upper= FALSE, as.mat= TRUE) # mat= TRUE,
2077
fred.array <- array(t(fred %*% Aimat),c(r,1,nn))
1985
ans = matrix(0,M,r)
1986
fred = mux22(t(wz), z-eta, M=M, upper= FALSE, as.mat= TRUE) # mat= TRUE,
1987
fred.array = array(t(fred %*% Aimat),c(r,1,nn))
2078
1988
for(s in 1:r) {
2079
a1 <- apply(fred %*% t(deta0.da[,,s]),2,sum)
2080
a2 <- apply(fred * etastar[,s],2,sum)
2081
temp92 <- array(detastar.da[,,s,],c(M,r,nn))
2082
temp93 <- mux7(temp92, fred.array)
2083
a3 <- apply(temp93,1:2,sum)
2084
ans[,s] <- a1 + a2 + a3
1989
a1 = colSums(fred %*% t(deta0.da[,,s]))
1990
a2 = colSums(fred * etastar[,s])
1991
temp92 = array(detastar.da[,,s,],c(M,r,nn))
1992
temp93 = mux7(temp92, fred.array)
1993
a3 = rowSums(temp93, dims=2)
1994
ans[,s] = a1 + a2 + a3
2085
1995
}
2086
1996
2087
ans <- -2 * c(ans[cbindex,])
1997
ans = -2 * c(ans[cbindex,])
2088
1998
2089
1999
ans
2090
2000
}
2231
2141
if(length(ellipse)) {
2232
2142
ellipse.temp = if(ellipse > 0) ellipse else 0.95
2233
2143
if(ellipse < 0 && (!object@control$EqualTolerances || varlvI))
2234
stop(paste("an equal-tolerances assumption and varlvI=FALSE",
2235
"is needed for \"ellipse\" < 0"))
2144
stop("an equal-tolerances assumption and 'varlvI=FALSE' ",
2145
"is needed for 'ellipse' < 0")
2236
2146
if( check.ok ) {
2237
2147
colx1.index = object@control$colx1.index
2238
2148
if(!(length(colx1.index)==1 &&
2364
2274
Aadj = rep(Aadj, len=length(index.nosz))
2365
2275
Acex = rep(Acex, len=length(index.nosz))
2366
2276
Acol = rep(Acol, len=length(index.nosz))
2367
if(length(Alabels) != M) stop(paste("Alabels must be of length", M))
2277
if(length(Alabels) != M) stop("'Alabels' must be of length ", M)
2368
2278
if(length(Apch)) {
2369
2279
Apch = rep(Apch, len=length(index.nosz))
2370
2280
for(i in index.nosz)
2385
2295
Clwd = rep(Clwd, len=p2)
2386
2296
Clty = rep(Clty, len=p2)
2387
2297
if(length(Clabels) != p2)
2388
stop(paste("length(Clabels) must be equal to", p2))
2389
for(i in 1:p2) {
2390
if(is.R()) arrows(0, 0, Cmat[i,1], Cmat[i,2],
2391
lwd=Clwd[i], lty=Clty[i], col=Ccol[i]) else
2392
arrows(0,0,Cmat[i,1],Cmat[i,2],open=TRUE,
2393
lwd=Clwd[i], lty=Clty[i], col=Ccol[i])
2394
const = 1 + gapC[i] / sqrt(Cmat[i,1]^2 + Cmat[i,2]^2)
2395
text(const*Cmat[i,1], const*Cmat[i,2], Clabels[i], cex=Ccex[i],
2396
adj=Cadj[i], col=Ccol[i])
2298
stop("'length(Clabels)' must be equal to ", p2)
2299
for(ii in 1:p2) {
2300
arrows(0, 0, Cmat[ii,1], Cmat[ii,2],
2301
lwd=Clwd[ii], lty=Clty[ii], col=Ccol[ii])
2302
const = 1 + gapC[ii] / sqrt(Cmat[ii,1]^2 + Cmat[ii,2]^2)
2303
text(const*Cmat[ii,1], const*Cmat[ii,2], Clabels[ii], cex=Ccex[ii],
2304
adj=Cadj[ii], col=Ccol[ii])
2397
2305
}
2398
2306
}
2399
2307
2407
2315
spch = rep(spch, len=n)
2408
2316
scol = rep(scol, len=n)
2409
2317
scex = rep(scex, len=n)
2410
for(i in ugrp) {
2411
gp = groups==i
2318
for(ii in ugrp) {
2319
gp = groups == ii
2412
2320
if(nlev > 1 && (length(unique(spch[gp])) != 1 ||
2413
2321
length(unique(scol[gp])) != 1 ||
2414
2322
length(unique(scex[gp])) != 1))
2415
warning(paste("spch/scol/scex is different for individuals",
2416
"from the same group"))
2323
warning("spch/scol/scex is different for individuals ",
2324
"from the same group")
2417
2325
2418
2326
temp = nuhat[gp,,drop=FALSE]
2419
2327
if(length(spch)) {
2426
2334
if(chull.arg) {
2427
2335
hull = chull(temp[,1],temp[,2])
2428
2336
hull = c(hull, hull[1])
2429
lines(temp[hull,1], temp[hull,2], type="b", lty=clty[i],
2430
col=ccol[i], lwd=clwd[i], pch=" ")
2337
lines(temp[hull,1], temp[hull,2], type="b", lty=clty[ii],
2338
col=ccol[ii], lwd=clwd[ii], pch=" ")
2431
2339
}
2432
2340
}
2433
2341
}
2497
2405
}
2498
2406
2499
2407
2500
if(is.R()) {
2501
if(!isGeneric("biplot"))
2408
if(!isGeneric("biplot"))
2502
2409
setGeneric("biplot", function(x, ...) standardGeneric("biplot"))
2503
}
2504
2410
2505
2411
2506
2412
setMethod("Coef", "qrrvglm", function(object, ...) Coef.qrrvglm(object, ...))
2541
2447
answer@dispersion =
2542
2448
answer@misc$dispersion = (answer@post$Coef)@dispersion
2543
2449
2544
class(answer) = "summary.qrrvglm"
2545
answer
2450
as(answer, "summary.qrrvglm")
2546
2451
}
2547
2452
2548
2453
printsummary.qrrvglm = function(x, ...) {
2568
2473
2569
2474
}
2570
2475
2571
setClass(Class="summary.qrrvglm", representation("qrrvglm"))
2476
2477
setClass("summary.qrrvglm", contains = "qrrvglm")
2478
2479
2480
2572
2481
2573
2482
setMethod("summary", "qrrvglm",
2574
2483
function(object, ...)
2575
2484
summary.qrrvglm(object, ...))
2576
2485
2577
setMethod("print", "summary.qrrvglm",
2578
function(x, ...)
2579
invisible(printsummary.qrrvglm(x, ...)))
2486
setMethod("print", "summary.qrrvglm",
2487
function(x, ...)
2488
invisible(printsummary.qrrvglm(x, ...)))
2580
2489
2581
setMethod("show", "summary.qrrvglm",
2582
function(object)
2583
invisible(printsummary.qrrvglm(object)))
2490
setMethod("show", "summary.qrrvglm",
2491
function(object)
2492
invisible(printsummary.qrrvglm(object)))
2584
2493
2585
2494
setMethod("print", "Coef.rrvglm", function(x, ...)
2586
2495
invisible(printCoef.rrvglm(x, ...)))
2604
2513
y = object.save@y
2605
2514
} else {
2606
2515
y = as.matrix(y)
2607
class(y) = "matrix" # Needed in R
2516
y = as(y, "matrix")
2608
2517
}
2609
2518
if(length(dim(y)) != 2 || nrow(y) < 3 || ncol(y) < 3)
2610
2519
stop("y must be a matrix with >= 3 rows & columns, or a rrvglm() object")
2621
2530
options(warn=warn.save)
2622
2531
2623
2532
Row = factor(1:nrow(y))
2624
modmat.row = if(is.R()) model.matrix(~ Row) else {
2625
tmp3 = contrasts(Row)
2626
dimnames(tmp3) = list(dimnames(tmp3)[[1]], paste("Row", 2:nrow(y), sep=""))
2627
cbind("(Intercept)"=1, tmp3)
2628
}
2533
modmat.row = model.matrix( ~ Row)
2629
2534
Col = factor(1:ncol(y))
2630
modmat.col = if(is.R()) model.matrix(~ Col) else {
2631
tmp3 = contrasts(Col)
2632
dimnames(tmp3) = list(dimnames(tmp3)[[1]], paste("Col", 2:ncol(y), sep=""))
2633
cbind("(Intercept)"=1, tmp3)
2634
}
2535
modmat.col = model.matrix( ~ Col)
2635
2536
2636
2537
cms = list("(Intercept)" = matrix(1, ncol(y), 1))
2637
for(i in 2:nrow(y)) {
2638
cms[[paste("Row", i, sep="")]] = matrix(1, ncol(y), 1)
2639
.grc.df[[paste("Row", i, sep="")]] = modmat.row[,i]
2640
}
2641
for(i in 2:ncol(y)) {
2642
cms[[paste("Col", i, sep="")]] = modmat.col[,i,drop=FALSE]
2643
.grc.df[[paste("Col", i, sep="")]] = rep(1, nrow(y))
2644
}
2645
for(i in 2:nrow(y)) {
2646
cms[[yn1[i]]] = diag(ncol(y))
2647
.grc.df[[yn1[i]]] = ei(i, nrow(y))
2538
for(ii in 2:nrow(y)) {
2539
cms[[paste("Row", ii, sep="")]] = matrix(1, ncol(y), 1)
2540
.grc.df[[paste("Row", ii, sep="")]] = modmat.row[,ii]
2541
}
2542
for(ii in 2:ncol(y)) {
2543
cms[[paste("Col", ii, sep="")]] = modmat.col[,ii,drop=FALSE]
2544
.grc.df[[paste("Col", ii, sep="")]] = rep(1, nrow(y))
2545
}
2546
for(ii in 2:nrow(y)) {
2547
cms[[yn1[ii]]] = diag(ncol(y))
2548
.grc.df[[yn1[ii]]] = ei(ii, nrow(y))
2648
2549
}
2649
2550
2650
2551
dimnames(.grc.df) = list(if(length(dimnames(y)[[1]])) dimnames(y)[[1]] else
2653
2554
2654
2555
str1 = "~ Row2"
2655
2556
if(nrow(y)>2)
2656
for(i in 3:nrow(y))
2657
str1 = paste(str1, paste("Row", i, sep=""), sep=" + ")
2658
for(i in 2:ncol(y))
2659
str1 = paste(str1, paste("Col", i, sep=""), sep=" + ")
2557
for(ii in 3:nrow(y))
2558
str1 = paste(str1, paste("Row", ii, sep=""), sep=" + ")
2559
for(ii in 2:ncol(y))
2560
str1 = paste(str1, paste("Col", ii, sep=""), sep=" + ")
2660
2561
str2 = paste("y ", str1)
2661
for(i in 2:nrow(y))
2662
str2 = paste(str2, yn1[i], sep=" + ")
2562
for(ii in 2:nrow(y))
2563
str2 = paste(str2, yn1[ii], sep=" + ")
2663
2564
myrrcontrol$Norrr = as.formula(str1) # Overwrite this
2664
2565
2665
if(is.R()) assign(".grc.df", .grc.df, envir = VGAMenv) else
2666
.grc.df <<- .grc.df
2566
assign(".grc.df", .grc.df, envir = VGAMenv)
2667
2567
2668
2568
warn.save = options()$warn
2669
2569
options(warn=-3) # Suppress the warnings (hopefully, temporarily)
2673
2573
options(warn=warn.save)
2674
2574
2675
2575
if(summary.arg) {
2676
class(answer) = "rrvglm"
2576
answer = as(answer, "rrvglm")
2577
2677
2578
answer = summary.rrvglm(answer, h.step=h.step)
2678
2579
} else {
2679
class(answer) = "grc"
2580
answer = as(answer, "grc")
2680
2581
}
2681
2582
2682
if(is.R()) {
2683
if(exists(".grc.df", envir = VGAMenv))
2684
rm(".grc.df", envir = VGAMenv)
2685
} else {
2686
remove(".grc.df")
2687
}
2583
if(exists(".grc.df", envir = VGAMenv))
2584
rm(".grc.df", envir = VGAMenv)
2688
2585
2689
2586
answer
2690
2587
}
2818
2715
MaxScale <- as.character(substitute(MaxScale))
2819
2716
MaxScale <- match.arg(MaxScale, c("predictors", "response"))[1]
2820
2717
if(MaxScale != "predictors")
2821
stop("can currently only handle MaxScale=\"predictors\"")
2718
stop("can currently only handle MaxScale='predictors'")
2822
2719
2823
2720
sobj = summary(object)
2824
2721
cobj = Coef(object, ITolerances = ITolerances, ...)
3109
3006
3110
3007
3111
3008
cgo <- function(...) {
3112
stop("The function \"cgo\" has been renamed \"cqo\". Ouch! Sorry!")
3009
stop("The function 'cgo' has been renamed 'cqo'. Ouch! Sorry!")
3113
3010
}
3114
3011
3115
3012
clo <- function(...) {
3116
stop("Constrained linear ordination is fitted with the function \"rrvglm\"")
3013
stop("Constrained linear ordination is fitted with the function 'rrvglm'")
3117
3014
}
3118
3015
3119
3016
Loggerhead is a web-based interface for Breezy
Version: 2.0.1