1# These functions are 2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland. 3# All rights reserved. 4 5 6 7 8 9 10replace.constraints <- function(Hlist, cm, index) { 11 12 for (iii in index) 13 Hlist[[iii]] <- cm 14 Hlist 15} # replace.constraints 16 17 18 19 20 valt.control <- function( 21 Alphavec = c(2, 4, 6, 9, 12, 16, 20, 25, 30, 40, 50, 22 60, 80, 100, 125, 2^(8:12)), 23 Criterion = c("ResSS", "coefficients"), 24 Linesearch = FALSE, 25 Maxit = 7, 26 Suppress.warning = TRUE, 27 Tolerance = 1e-7, ...) { 28 29 if (mode(Criterion) != "character" && mode(Criterion) != "name") 30 Criterion <- as.character(substitute(Criterion)) 31 Criterion <- match.arg(Criterion, c("ResSS", "coefficients"))[1] 32 33 list(Alphavec = Alphavec, 34 Criterion = Criterion, 35 Linesearch = Linesearch, 36 Maxit = Maxit, 37 Suppress.warning = Suppress.warning, 38 Tolerance = Tolerance) 39} # valt.control 40 41 42 43 44 45qrrvglm.xprod <- function(numat, Aoffset, Quadratic, I.tolerances) { 46 Rank <- ncol(numat) 47 moff <- NULL 48 ans <- if (Quadratic) { 49 index <- iam(NA, NA, M = Rank, diag = TRUE, both = TRUE) 50 temp1 <- cbind(numat[, index$row] * numat[, index$col]) 51 if (I.tolerances) { 52 moff <- 0 53 for (ii in 1:Rank) 54 moff <- moff - 0.5 * temp1[, ii] 55 } 56 cbind(numat, if (I.tolerances) NULL else temp1) 57 } else { 58 as.matrix(numat) 59 } 60 list(matrix = if (Aoffset > 0) ans else ans[, -(1:Rank), drop = FALSE], 61 offset = moff) 62} # qrrvglm.xprod 63 64 65 66 67 68 69 valt <- function(x, z, U, Rank = 1, 70 Hlist = NULL, 71 Cinit = NULL, 72 Alphavec = c(2, 4, 6, 9, 12, 16, 20, 25, 30, 40, 50, 73 60, 80, 100, 125, 2^(8:12)), 74 Criterion = c("ResSS", "coefficients"), 75 Crow1positive = rep_len(TRUE, Rank), 76 colx1.index, 77 Linesearch = FALSE, 78 Maxit = 20, 79 str0 = NULL, 80 sd.Cinit = 0.02, 81 Suppress.warning = FALSE, 82 Tolerance = 1e-6, 83 trace = FALSE, 84 xij = NULL) { 85 86 87 88 89 90 91 92 if (mode(Criterion) != "character" && mode(Criterion) != "name") 93 Criterion <- as.character(substitute(Criterion)) 94 Criterion <- match.arg(Criterion, c("ResSS", "coefficients"))[1] 95 96 if (any(diff(Alphavec) <= 0)) 97 stop("'Alphavec' must be an increasing sequence") 98 99 if (!is.matrix(z)) 100 z <- as.matrix(z) 101 n <- nrow(z) 102 M <- ncol(z) 103 if (!is.matrix(x)) 104 x <- as.matrix(x) 105 106 colx2.index <- if (is.null(colx1.index)) 1:ncol(x) else 107 (1:ncol(x))[-colx1.index] 108 109 p1 <- length(colx1.index) 110 p2 <- length(colx2.index) 111 p <- p1 + p2 112 if (!p2) 113 stop("'p2', the number of variables for the ", 114 "reduced-rank regression, must be > 0") 115 116 if (!length(Hlist)) { 117 Hlist <- replace.constraints(vector("list", p), diag(M), 1:p) 118 } 119 120 dU <- dim(U) 121 if (dU[2] != n) 122 stop("input unconformable") 123 124 clist2 <- replace.constraints(vector("list", Rank+p1), 125 if (length(str0)) 126 diag(M)[, -str0, drop = FALSE] else 127 diag(M), 1:Rank) 128 if (p1) { 129 for (kk in 1:p1) 130 clist2[[Rank+kk]] <- Hlist[[colx1.index[kk]]] 131 } 132 133 if (is.null(Cinit)) 134 Cinit <- matrix(rnorm(p2*Rank, sd = sd.Cinit), p2, Rank) 135 136 fit <- list(ResSS = 0) # Only for initial old.crit below 137 138 C <- Cinit # This is input for the main iter loop 139 old.crit <- switch(Criterion, coefficients = C, ResSS = fit$ResSS) 140 141 recover <- 0 # Allow a few iterations between different line searches 142 for (iter in 1:Maxit) { 143 iter.save <- iter 144 145 latvar.mat <- x[, colx2.index, drop = FALSE] %*% C 146 new.latvar.model.matrix <- cbind(latvar.mat, 147 if (p1) x[, colx1.index] else NULL) 148 fit <- vlm.wfit(xmat = new.latvar.model.matrix, z, Hlist = clist2, 149 U = U, matrix.out = TRUE, is.vlmX = FALSE, 150 ResSS = FALSE, qr = FALSE, xij = xij) 151 A <- t(fit$mat.coef[1:Rank, , drop = FALSE]) 152 153 clist1 <- replace.constraints(Hlist, A, colx2.index) 154 fit <- vlm.wfit(xmat = x, z, Hlist = clist1, U = U, 155 matrix.out = TRUE, is.vlmX = FALSE, 156 ResSS = TRUE, qr = FALSE, xij = xij) 157 C <- fit$mat.coef[colx2.index, , drop = FALSE] %*% A %*% 158 solve(t(A) %*% A) 159 160 numat <- x[, colx2.index, drop = FALSE] %*% C 161 evnu <- eigen(var(numat), symmetric = TRUE) 162 temp7 <- if (Rank > 1) evnu$vector %*% diag(evnu$value^(-0.5)) else 163 evnu$vector %*% evnu$value^(-0.5) 164 C <- C %*% temp7 165 A <- A %*% t(solve(temp7)) 166 temp8 <- crow1C(cmat = C, Crow1positive, amat = A) 167 C <- temp8$cmat 168 A <- temp8$amat 169 170 171 ratio <- 172 switch(Criterion, 173 coefficients = max(abs(C - old.crit) / ( 174 Tolerance + abs(C))), 175 ResSS = max(abs(fit$ResSS - old.crit) / ( 176 Tolerance + fit$ResSS))) 177 178 if (trace) { 179 cat(" Alternating iteration", iter, 180 ", Convergence criterion = ", format(ratio), "\n") 181 if (!is.null(fit$ResSS)) 182 cat(" ResSS = ", fit$ResSS, "\n") 183 flush.console() 184 } 185 186 if (ratio < Tolerance) { 187 if (!Linesearch || (Linesearch && iter >= 3)) 188 break 189 } else if (iter == Maxit && !Suppress.warning) { 190 warning("did not converge") 191 } 192 193 fini.linesearch <- FALSE 194 if (Linesearch && iter - recover >= 2) { 195 xnew <- C 196 197 direction1 <- (xnew - xold) # / sqrt(1 + sum((xnew-xold)^2)) 198 ftemp <- fit$ResSS # Most recent objective function 199 use.alpha <- 0 # The current step relative to (xold, yold) 200 for (itter in seq_along(Alphavec)) { 201 CC <- xold + Alphavec[itter] * direction1 202 203 try.latvar.mat <- x[, colx2.index, drop = FALSE] %*% CC 204 try.new.latvar.model.matrix <- 205 cbind(try.latvar.mat, 206 if (p1) x[, colx1.index] else NULL) 207 208 try <- vlm.wfit(xmat = try.new.latvar.model.matrix, z, 209 Hlist = clist2, U = U, matrix.out = TRUE, 210 is.vlmX = FALSE, ResSS = TRUE, qr = FALSE, 211 xij = xij) 212 if (try$ResSS < ftemp) { 213 use.alpha <- Alphavec[itter] 214 fit <- try 215 ftemp <- try$ResSS 216 C <- CC 217 A <- t(fit$mat.coef[1:Rank, , drop = FALSE]) 218 latvar.mat <- x[, colx2.index, drop = FALSE] %*% C 219 recover <- iter # Give it some altg iters to recover 220 } else { 221 if (trace && use.alpha > 0) { 222 cat(" Finished line search using Alpha = ", 223 use.alpha, "\n") 224 flush.console() 225 } 226 fini.linesearch <- TRUE 227 } 228 if (fini.linesearch) break 229 } # End of itter loop 230 } 231 232 xold <- C # Do not take care of drift 233 old.crit <- switch(Criterion, 234 coefficients = C, 235 ResSS = fit$ResSS) 236 } # End of iter loop 237 238 list(A = A, 239 C = C, 240 fitted = fit$fitted, 241 new.coeffs = fit$coef, 242 ResSS = fit$ResSS) 243} # valt 244 245 246 247 248 249 250 lm2qrrvlm.model.matrix <- 251 function(x, Hlist, C, control, assign = TRUE, 252 no.thrills = FALSE) { 253 254 Rank <- control$Rank 255 colx1.index <- control$colx1.index 256 Quadratic <- control$Quadratic 257 Dzero <- control$Dzero 258 Corner <- control$Corner 259 I.tolerances <- control$I.tolerances 260 261 M <- nrow(Hlist[[1]]) 262 p1 <- length(colx1.index) 263 combine2 <- c(control$str0, 264 if (Corner) control$Index.corner else NULL) 265 266 Qoffset <- if (Quadratic) ifelse(I.tolerances, 0, sum(1:Rank)) else 0 267 NoA <- length(combine2) == M # No unknown parameters in A 268 clist2 <- if (NoA) { 269 Aoffset <- 0 270 vector("list", Aoffset+Qoffset+p1) 271 } else { 272 Aoffset <- Rank 273 replace.constraints(vector("list", Aoffset+Qoffset+p1), 274 if (length(combine2)) diag(M)[, -combine2, drop = FALSE] else 275 diag(M), 276 1:Rank) # If Corner then does not contain \bI_{Rank} 277 } 278 if (Quadratic && !I.tolerances) 279 clist2 <- replace.constraints(clist2, 280 if (control$eq.tolerances) 281 matrix(1, M, 1) - eijfun(Dzero, M) else { 282 if (length(Dzero)) diag(M)[,-Dzero, drop = FALSE] else diag(M)}, 283 Aoffset + (1:Qoffset)) 284 if (p1) 285 for (kk in 1:p1) 286 clist2[[Aoffset+Qoffset+kk]] <- Hlist[[colx1.index[kk]]] 287 if (!no.thrills) { 288 i63 <- iam(NA, NA, M = Rank, both = TRUE) 289 names(clist2) <- c( 290 if (NoA) NULL else paste("(latvar", 1:Rank, ")", sep = ""), 291 if (Quadratic && Rank == 1 && !I.tolerances) 292 "(latvar^2)" else 293 if (Quadratic && Rank>1 && !I.tolerances) 294 paste("(latvar", i63$row, ifelse(i63$row == i63$col, "^2", 295 paste("*latvar", i63$col, sep = "")), ")", sep = "") else 296 NULL, 297 if (p1) names(colx1.index) else NULL) 298 } 299 300 latvar.mat <- x[, control$colx2.index, drop = FALSE] %*% C 301 302 303 tmp900 <- qrrvglm.xprod(latvar.mat, Aoffset, Quadratic, I.tolerances) 304 new.latvar.model.matrix <- cbind(tmp900$matrix, 305 if (p1) x[,colx1.index] else NULL) 306 if (!no.thrills) 307 dimnames(new.latvar.model.matrix) <- list(dimnames(x)[[1]], 308 names(clist2)) 309 310 if (assign) { 311 asx <- attr(x, "assign") 312 asx <- vector("list", ncol(new.latvar.model.matrix)) 313 names(asx) <- names(clist2) 314 for (ii in seq_along(names(asx))) { 315 asx[[ii]] <- ii 316 } 317 attr(new.latvar.model.matrix, "assign") <- asx 318 } 319 320 if (no.thrills) 321 list(new.latvar.model.matrix = new.latvar.model.matrix, 322 constraints = clist2, 323 offset = tmp900$offset) else 324 list(new.latvar.model.matrix = new.latvar.model.matrix, 325 constraints = clist2, 326 NoA = NoA, 327 Aoffset = Aoffset, 328 latvar.mat = latvar.mat, 329 offset = tmp900$offset) 330} # lm2qrrvlm.model.matrix 331 332 333 334valt.2iter <- function(x, z, U, Hlist, A, control) { 335 336 337 clist1 <- replace.constraints(Hlist, A, control$colx2.index) 338 fit <- vlm.wfit(xmat = x, z, Hlist = clist1, U = U, matrix.out = TRUE, 339 is.vlmX = FALSE, ResSS = TRUE, 340 qr = FALSE, xij = control$xij) 341 C <- fit$mat.coef[control$colx2.index, , drop = FALSE] %*% 342 A %*% solve(t(A) %*% A) 343 344 list(A = A, C = C, 345 fitted = fit$fitted, new.coeffs = fit$coef, 346 Hlist = clist1, ResSS = fit$ResSS) 347} # valt.2iter 348 349 350 351valt.1iter <- function(x, z, U, Hlist, C, control, 352 lp.names = NULL, nice31 = FALSE, 353 MSratio = 1) { 354 355 Rank <- control$Rank 356 Quadratic <- control$Quadratic 357 Index.corner <- control$Index.corner 358 p1 <- length(control$colx1.index) 359 M <- ncol(zedd <- as.matrix(z)) 360 NOS <- M / MSratio 361 Corner <- control$Corner 362 I.tolerances <- control$I.tolerances 363 364 Qoffset <- if (Quadratic) ifelse(I.tolerances, 0, sum(1:Rank)) else 0 365 tmp833 <- lm2qrrvlm.model.matrix(x = x, Hlist = Hlist, C = C, 366 control = control) 367 new.latvar.model.matrix <- tmp833$new.latvar.model.matrix 368 clist2 <- tmp833$constraints # Does not contain \bI_{Rank} 369 latvar.mat <- tmp833$latvar.mat 370 if (Corner) 371 zedd[,Index.corner] <- zedd[,Index.corner] - latvar.mat 372 373 if (nice31 && MSratio == 1) { 374 fit <- list(mat.coef = NULL, fitted.values = NULL, ResSS = 0) 375 376 clist2 <- NULL # for vlm.wfit 377 378 i5 <- rep_len(0, MSratio) 379 for (ii in 1:NOS) { 380 i5 <- i5 + 1:MSratio 381 382 tmp100 <- vlm.wfit(xmat = new.latvar.model.matrix, 383 zedd[, i5, drop = FALSE], 384 Hlist = clist2, 385 U = U[i5,, drop = FALSE], 386 matrix.out = TRUE, 387 is.vlmX = FALSE, ResSS = TRUE, 388 qr = FALSE, 389 Eta.range = control$Eta.range, 390 xij = control$xij, 391 lp.names = lp.names[i5]) 392 fit$ResSS <- fit$ResSS + tmp100$ResSS 393 fit$mat.coef <- cbind(fit$mat.coef, tmp100$mat.coef) 394 fit$fitted.values <- cbind(fit$fitted.values, 395 tmp100$fitted.values) 396 } 397 } else { 398 fit <- vlm.wfit(xmat = new.latvar.model.matrix, 399 zedd, Hlist = clist2, U = U, 400 matrix.out = TRUE, 401 is.vlmX = FALSE, ResSS = TRUE, qr = FALSE, 402 Eta.range = control$Eta.range, 403 xij = control$xij, lp.names = lp.names) 404 } 405 A <- if (tmp833$NoA) matrix(0, M, Rank) else 406 t(fit$mat.coef[1:Rank,, drop = FALSE]) 407 if (Corner) 408 A[Index.corner,] <- diag(Rank) 409 410 B1 <- if (p1) 411 fit$mat.coef[-(1:(tmp833$Aoffset+Qoffset)),, drop = FALSE] else 412 NULL 413 fv <- as.matrix(fit$fitted.values) 414 if (Corner) 415 fv[,Index.corner] <- fv[,Index.corner] + latvar.mat 416 Dmat <- if (Quadratic) { 417 if (I.tolerances) { 418 tmp800 <- matrix(0, M, Rank*(Rank+1)/2) 419 tmp800[if (MSratio == 2) c(TRUE, FALSE) else 420 TRUE, 1:Rank] <- -0.5 421 tmp800 422 } else 423 t(fit$mat.coef[(tmp833$Aoffset+1): 424 (tmp833$Aoffset+Qoffset),, drop = FALSE]) 425 } else 426 NULL 427 428 list(Amat = A, B1 = B1, Cmat = C, Dmat = Dmat, 429 fitted = if (M == 1) c(fv) else fv, 430 new.coeffs = fit$coef, constraints = clist2, ResSS = fit$ResSS, 431 offset = if (length(tmp833$offset)) tmp833$offset else NULL) 432} # valt.1iter 433 434 435 436 437 438rrr.init.expression <- expression({ 439 if (length(control$Quadratic) && control$Quadratic) 440 copy.X.vlm <- TRUE 441 442 443 444 445 if (function.name %in% c("cqo", "cao")) { 446 447 modelno <- switch(family@vfamily[1], "poissonff" = 2, 448 "quasipoissonff" = 2, "quasipoisson" = 2, 449 "binomialff" = 1, "quasibinomialff" = 1, 450 "quasibinomial" = 1, "negbinomial" = 3, 451 "gamma2" = 5, "gaussianff" = 8, 452 0) # stop("cannot fit this model using fast algorithm") 453 if (modelno == 1) modelno = get("modelno", envir = VGAMenv) 454 rrcontrol$modelno = control$modelno = modelno 455 if (modelno == 3 || modelno == 5) { 456 457 458 M <- 2 * ifelse(is.matrix(y), ncol(y), 1) 459 control$str0 <- 460 rrcontrol$str0 <- seq(from = 2, to = M, by = 2) # Handles A 461 control$Dzero <- 462 rrcontrol$Dzero <- seq(from = 2, to = M, by = 2) # Handles D 463 464 465 } 466 } else { 467 modelno <- 0 # Any value will do as the variable is unused. 468 } 469 470 471}) # rrr.init.expression 472 473 474 475 476 477rrr.alternating.expression <- expression({ 478 479 alt <- valt(x, z, U, Rank = Rank, 480 Hlist = Hlist, 481 Cinit = rrcontrol$Cinit, 482 Criterion = rrcontrol$Criterion, 483 colx1.index = rrcontrol$colx1.index, 484 Linesearch = rrcontrol$Linesearch, 485 Maxit = rrcontrol$Maxit, 486 str0 = rrcontrol$str0, 487 sd.Cinit = rrcontrol$sd.Cinit, 488 Suppress.warning = rrcontrol$Suppress.warning, 489 Tolerance = rrcontrol$Tolerance, 490 trace = trace, 491 xij = control$xij) # This is subject to drift in A and C 492 493 ans2 <- rrr.normalize(rrcontrol = rrcontrol, A=alt$A, C=alt$C, x = x) 494 495 Amat <- ans2$A # Fed into Hlist below (in rrr.end.expression) 496 tmp.fitted <- alt$fitted # Also fed; was alt2$fitted 497 498 rrcontrol$Cinit <- ans2$C # For next valt() call 499 500 eval(rrr.end.expression) # Put Amat into Hlist, and create new z 501}) # rrr.alternating.expression 502 503 504 505 506 507 adjust.Dmat.expression <- function(Mmat, Rank, Dmat, M) { 508 509 if (length(Dmat)) { 510 ind0 <- iam(NA, NA, both = TRUE, M = Rank) 511 for (kay in 1:M) { 512 elts <- Dmat[kay, , drop = FALSE] # Manual recycling 513 if (length(elts) < Rank) 514 elts <- matrix(elts, 1, Rank) 515 Dk <- m2a(elts, M = Rank)[, , 1] 516 Dk <- matrix(Dk, Rank, Rank) 517 Dk <- t(Mmat) %*% Dk %*% Mmat # 20030822 Not diagonal in general 518 Dmat[kay, ] <- Dk[cbind(ind0$row.index[1:ncol(Dmat)], 519 ind0$col.index[1:ncol(Dmat)])] 520 } 521 } 522 Dmat 523 } # adjust.Dmat.expression 524 525 526 527 528 529rrr.normalize <- function(rrcontrol, A, C, x, Dmat = NULL) { 530 531 532 533 colx2.index <- rrcontrol$colx2.index 534 Rank <- rrcontrol$Rank 535 Index.corner <- rrcontrol$Index.corner 536 M <- nrow(A) 537 C.old <- C 538 539 if (rrcontrol$Corner) { 540 tmp87 <- A[Index.corner,, drop = FALSE] 541 Mmat <- solve(tmp87) # The normalizing matrix 542 C <- C %*% t(tmp87) 543 A <- A %*% Mmat 544 A[Index.corner,] <- diag(Rank) # Make sure 545 546 Dmat <- adjust.Dmat.expression(Mmat = Mmat, Rank = Rank, 547 Dmat = Dmat, M = M) 548 } 549 550 if (rrcontrol$Svd.arg) { 551 temp <- svd(C %*% t(A)) 552 if (!is.matrix(temp$v)) 553 temp$v <- as.matrix(temp$v) 554 C <- temp$u[, 1:Rank, drop = FALSE] %*% 555 diag(temp$d[1:Rank]^(1-rrcontrol$Alpha), nrow = Rank) 556 A <- diag(temp$d[1:Rank]^( rrcontrol$Alpha), nrow = Rank) %*% 557 t(temp$v[, 1:Rank, drop = FALSE]) 558 A <- t(A) 559 Mmat <- t(C.old) %*% C.old %*% solve(t(C) %*% C.old) 560 561 562 Dmat <- adjust.Dmat.expression(Mmat = Mmat, Rank = Rank, 563 Dmat = Dmat, M = M) 564 } 565 566 if (rrcontrol$Uncorrelated.latvar) { 567 latvar.mat <- x[, colx2.index, drop = FALSE] %*% C 568 var.latvar.mat <- var(latvar.mat) 569 UU <- chol(var.latvar.mat) 570 Ut <- solve(UU) 571 Mmat <- t(UU) 572 C <- C %*% Ut 573 A <- A %*% t(UU) 574 575 576 577 Dmat <- adjust.Dmat.expression(Mmat = Mmat, Rank = Rank, 578 Dmat = Dmat, M = M) 579 } 580 581 582 if (rrcontrol$Quadratic) { 583 Mmat <- diag(Rank) 584 for (LV in 1:Rank) 585 if (( rrcontrol$Crow1positive[LV] && C[1,LV] < 0) || 586 (!rrcontrol$Crow1positive[LV] && C[1,LV] > 0)) { 587 C[,LV] <- -C[,LV] 588 A[,LV] <- -A[,LV] 589 Mmat[LV,LV] <- -1 590 } 591 592 593 594 Dmat <- adjust.Dmat.expression(Mmat = Mmat, Rank = Rank, 595 Dmat = Dmat, M = M) 596 } 597 598 599 list(Amat = A, Cmat = C, Dmat = Dmat) 600} # rrr.normalize 601 602 603 604 605rrr.end.expression <- expression({ 606 607 if (exists(".VGAM.etamat", envir = VGAMenv)) 608 rm(".VGAM.etamat", envir = VGAMenv) 609 610 611 if (control$Quadratic) { 612 if (!length(extra)) 613 extra <- list() 614 extra$Cmat <- Cmat # Saves the latest iteration 615 extra$Dmat <- Dmat # Not the latest iteration 616 extra$B1 <- B1.save # Not the latest iteration (not good) 617 } else { 618 Hlist <- replace.constraints(Hlist.save, Amat, colx2.index) 619 } 620 621 X.vlm.save <- if (control$Quadratic) { 622 tmp300 <- lm2qrrvlm.model.matrix(x = x, Hlist = Hlist.save, 623 C = Cmat, control = control) 624 latvar.mat <- tmp300$latvar.mat # Needed at the top of new.s.call 625 626 lm2vlm.model.matrix(tmp300$new.latvar.model.matrix, 627 H.list, 628 xij = control$xij) 629 } else { 630 lm2vlm.model.matrix(x, Hlist, xij = control$xij) 631 } 632 633 634 fv <- tmp.fitted # Contains \bI \bnu 635 eta <- fv + offset 636 if (FALSE && control$Rank == 1) { 637 ooo <- order(latvar.mat[, 1]) 638 } 639 mu <- family@linkinv(eta, extra) 640 641 if (anyNA(mu)) 642 warning("there are NAs in mu") 643 644 deriv.mu <- eval(family@deriv) 645 wz <- eval(family@weight) 646 if (control$checkwz) 647 wz <- checkwz(wz, M = M, trace = trace, 648 wzepsilon = control$wzepsilon) 649 U <- vchol(wz, M = M, n = n, silent=!trace) 650 tvfor <- vforsub(U, as.matrix(deriv.mu), M = M, n = n) 651 z <- eta + vbacksub(U, tvfor, M = M, n = n) - 652 offset # Contains \bI \bnu 653 654 655 656}) # rrr.end.expression 657 658 659 660 661 662rrr.derivative.expression <- expression({ 663 664 665 666 667 668 669 which.optimizer <- if (control$Quadratic && control$FastAlgorithm) { 670 "BFGS" 671 } else { 672 if (iter <= rrcontrol$Switch.optimizer) "Nelder-Mead" else "BFGS" 673 } 674 if (trace && control$OptimizeWrtC) { 675 cat("\n\n") 676 cat("Using", which.optimizer, "\n") 677 flush.console() 678 } 679 680 constraints <- replace.constraints(constraints, diag(M), 681 rrcontrol$colx2.index) 682 nice31 <- (!control$eq.tol || control$I.tolerances) && 683 all(trivial.constraints(constraints) == 1) 684 685 theta0 <- c(Cmat) 686 assign(".VGAM.dot.counter", 0, envir = VGAMenv) 687 if (control$OptimizeWrtC) { 688 if (control$Quadratic && control$FastAlgorithm) { 689 if (iter == 2) { 690 if (exists(".VGAM.etamat", envir = VGAMenv)) 691 rm(".VGAM.etamat", envir = VGAMenv) 692 } 693 if (iter > 2 && !quasi.newton$convergence) { 694 if (zthere <- exists(".VGAM.z", envir = VGAMenv)) { 695 ..VGAM.z <- get(".VGAM.z", envir = VGAMenv) 696 ..VGAM.U <- get(".VGAM.U", envir = VGAMenv) 697 ..VGAM.beta <- get(".VGAM.beta", envir = VGAMenv) 698 } 699 if (zthere) { 700 z <- matrix(..VGAM.z, n, M) # minus any offset 701 U <- matrix(..VGAM.U, M, n) 702 } 703 704 } 705 706 if (iter == 2 || quasi.newton$convergence) { 707 NOS <- ifelse(modelno == 3 || modelno == 5, M/2, M) 708 709 canfitok <- 710 (exists("CQO.FastAlgorithm", envir=VGAMenv) && 711 get("CQO.FastAlgorithm", envir = VGAMenv)) 712 if (!canfitok) 713 stop("cannot fit this model using fast algorithm") 714 p2star <- if (nice31) 715 ifelse(control$I.toleran, 716 Rank, 717 Rank+0.5*Rank*(Rank+1)) else 718 (NOS*Rank + 719 Rank*(Rank+1)/2 * ifelse(control$eq.tol, 1,NOS)) 720 p1star <- if (nice31) p1 * 721 ifelse(modelno == 3 || modelno == 5, 2, 1) else 722 (ncol(X.vlm.save) - p2star) 723 X.vlm.1save <- if (p1star > 0) 724 X.vlm.save[,-(1:p2star)] else NULL 725 quasi.newton <- 726 optim(par = Cmat, fn = callcqof, 727 gr <- if (control$GradientFunction) calldcqo else NULL, 728 method = which.optimizer, 729 control = list(fnscale = 1, 730 trace = as.integer(control$trace), 731 parscale = rep_len(control$Parscale, 732 length(Cmat)), 733 maxit = 250), 734 etamat = eta, xmat = x, ymat = y, wvec = w, 735 X.vlm.1save = if (nice31) NULL else X.vlm.1save, 736 modelno = modelno, Control = control, 737 n = n, M = M, p1star = p1star, 738 p2star = p2star, nice31 = nice31) 739 740 741 if (zthere <- exists(".VGAM.z", envir = VGAMenv)) { 742 ..VGAM.z <- get(".VGAM.z", envir = VGAMenv) 743 ..VGAM.U <- get(".VGAM.U", envir = VGAMenv) 744 ..VGAM.beta <- get(".VGAM.beta", envir = VGAMenv) 745 } 746 if (zthere) { 747 z <- matrix(..VGAM.z, n, M) # minus any offset 748 U <- matrix(..VGAM.U, M, n) 749 } 750 } else { 751 if (exists(".VGAM.offset", envir = VGAMenv)) 752 rm(".VGAM.offset", envir = VGAMenv) 753 } 754 } else { 755 use.reltol <- if (length(rrcontrol$Reltol) >= iter) 756 rrcontrol$Reltol[iter] else rev(rrcontrol$Reltol)[1] 757 quasi.newton <- 758 optim(par = theta0, 759 fn = rrr.derivC.ResSS, 760 method = which.optimizer, 761 control = list(fnscale = rrcontrol$Fnscale, 762 maxit = rrcontrol$Maxit, 763 abstol = rrcontrol$Abstol, 764 reltol = use.reltol), 765 U = U, z = if (control$I.tolerances) z + offset else z, 766 M = M, xmat = x, # varbix2 = varbix2, 767 Hlist = Hlist, rrcontrol = rrcontrol) 768 } 769 770 771 772 773 Cmat <- matrix(quasi.newton$par, p2, Rank, byrow = FALSE) 774 775 if (Rank > 1 && rrcontrol$I.tolerances) { 776 numat <- x[, rrcontrol$colx2.index, drop = FALSE] %*% Cmat 777 evnu <- eigen(var(numat), symmetric = TRUE) 778 Cmat <- Cmat %*% evnu$vector 779 numat <- x[, rrcontrol$colx2.index, drop = FALSE] %*% Cmat 780 offset <- if (Rank > 1) -0.5*rowSums(numat^2) else -0.5*numat^2 781 } 782 } 783 784 785 alt <- valt.1iter(x = x, z = z, U = U, Hlist = Hlist, 786 C = Cmat, nice31 = nice31, 787 control = rrcontrol, 788 lp.names = predictors.names) 789 790 791 if (length(alt$offset)) 792 offset <- alt$offset 793 794 B1.save <- alt$B1 # Put later into extra 795 tmp.fitted <- alt$fitted # contains \bI_{Rank} \bnu if Corner 796 797 if (modelno != 33 && control$OptimizeWrtC) 798 alt <- rrr.normalize(rrc = rrcontrol, A = alt$Amat, C = alt$Cmat, 799 x = x, Dmat = alt$Dmat) 800 801 if (trace && control$OptimizeWrtC) { 802 cat("\n") 803 cat(which.optimizer, "using optim():\n") 804 cat("Objective = ", quasi.newton$value, "\n") 805 cat("Parameters (= c(C)) = ", if (length(quasi.newton$par) < 5) 806 "" else "\n") 807 cat(alt$Cmat, fill = TRUE) 808 cat("\n") 809 cat("Number of function evaluations = ", 810 quasi.newton$count[1], "\n") 811 if (length(quasi.newton$message)) 812 cat("Message = ", quasi.newton$message, "\n") 813 cat("\n") 814 flush.console() 815 } 816 817 818 819 Amat <- alt$Amat # Needed in rrr.end.expression 820 Cmat <- alt$Cmat # Needed in rrr.end.expression if Quadratic 821 Dmat <- alt$Dmat # Put later into extra 822 823 eval(rrr.end.expression) # Put Amat into Hlist, and create new z 824}) # rrr.derivative.expression 825 826 827 828 829 830rrr.derivC.ResSS <- function(theta, U, z, M, xmat, Hlist, rrcontrol, 831 omit.these = NULL) { 832 833 if (rrcontrol$trace) { 834 cat(".") 835 flush.console() 836 } 837 alreadyThere <- exists(".VGAM.dot.counter", envir = VGAMenv) 838 if (alreadyThere) { 839 VGAM.dot.counter <- get(".VGAM.dot.counter", envir = VGAMenv) 840 VGAM.dot.counter <- VGAM.dot.counter + 1 841 assign(".VGAM.dot.counter", VGAM.dot.counter, envir = VGAMenv) 842 if (VGAM.dot.counter > max(50, options()$width - 5)) { 843 if (rrcontrol$trace) { 844 cat("\n") 845 flush.console() 846 } 847 assign(".VGAM.dot.counter", 0, envir = VGAMenv) 848 } 849 } 850 851 Cmat <- matrix(theta, length(rrcontrol$colx2.index), rrcontrol$Rank) 852 853 854 tmp700 <- 855 lm2qrrvlm.model.matrix(x = xmat, Hlist = Hlist, 856 no.thrills = !rrcontrol$Corner, 857 C = Cmat, control = rrcontrol, assign = FALSE) 858 Hlist <- tmp700$constraints # Does not contain \bI_{Rank} \bnu 859 860 if (rrcontrol$Corner) { 861 z <- as.matrix(z) # should actually call this zedd 862 z[, rrcontrol$Index.corner] <- 863 z[, rrcontrol$Index.corner] - tmp700$latvar.mat 864 } 865 866 if (length(tmp700$offset)) z <- z - tmp700$offset 867 868 869 vlm.wfit(xmat = tmp700$new.latvar.model.matrix, zmat = z, 870 Hlist = Hlist, ncolx = ncol(xmat), U = U, only.ResSS = TRUE, 871 matrix.out = FALSE, is.vlmX = FALSE, ResSS = TRUE, 872 qr = FALSE, Eta.range = rrcontrol$Eta.range, 873 xij = rrcontrol$xij)$ResSS 874} # rrr.derivC.ResSS 875 876 877 878 879rrvglm.optim.control <- function(Fnscale = 1, 880 Maxit = 100, 881 Switch.optimizer = 3, 882 Abstol = -Inf, 883 Reltol = sqrt(.Machine$double.eps), 884 ...) { 885 886 887 888 889 list(Fnscale = Fnscale, 890 Maxit = Maxit, 891 Switch.optimizer = Switch.optimizer, 892 Abstol = Abstol, 893 Reltol = Reltol) 894} # rrvglm.optim.control 895 896 897 898nlminbcontrol <- function(Abs.tol = 10^(-6), 899 Eval.max = 91, 900 Iter.max = 91, 901 Rel.err = 10^(-6), 902 Rel.tol = 10^(-6), 903 Step.min = 10^(-6), 904 X.tol = 10^(-6), 905 ...) { 906 907 908 list(Abs.tol = Abs.tol, 909 Eval.max = Eval.max, 910 Iter.max = Iter.max, 911 Rel.err = Rel.err, 912 Rel.tol = Rel.tol, 913 Step.min = Step.min, 914 X.tol = X.tol) 915} # nlminbcontrol 916 917 918 919 920 921 922Coef.qrrvglm <- 923 function(object, 924 varI.latvar = FALSE, 925 refResponse = NULL, ...) { 926 927 928 929 930 if (length(varI.latvar) != 1 || !is.logical(varI.latvar)) 931 stop("argument 'varI.latvar' must be TRUE or FALSE") 932 if (length(refResponse) > 1) 933 stop("argument 'refResponse' must be of length 0 or 1") 934 if (length(refResponse) && 935 is.Numeric(refResponse)) 936 if (!is.Numeric(refResponse, length.arg = 1, 937 integer.valued = TRUE)) 938 stop("bad input for argument 'refResponse'") 939 if (!is.logical(ConstrainedQO <- object@control$ConstrainedQO)) 940 stop("cannot determine whether the model is constrained or not") 941 942 ocontrol <- object@control 943 coef.object <- object@coefficients # Unscaled 944 Rank <- ocontrol$Rank 945 M <- object@misc$M 946 NOS <- if (length(object@y)) NCOL(object@y) else M 947 MSratio <- M / NOS # First value is g(mean) = quadratic form in latvar 948 Quadratic <- if (ConstrainedQO) ocontrol$Quadratic else TRUE 949 if (!Quadratic) stop("object is not a quadratic ordination object") 950 p1 <- length(ocontrol$colx1.index) 951 p2 <- length(ocontrol$colx2.index) 952 Index.corner <- ocontrol$Index.corner 953 str0 <- ocontrol$str0 954 eq.tolerances <- ocontrol$eq.tolerances 955 Dzero <- ocontrol$Dzero 956 Corner <- if (ConstrainedQO) ocontrol$Corner else FALSE 957 958 estI.tol <- if (ConstrainedQO) object@control$I.tolerances else FALSE 959 modelno <- object@control$modelno # 1, 2, 3, 4, 5, 6, 7 or 0 960 combine2 <- c(str0, if (Corner) Index.corner else NULL) 961 NoA <- length(combine2) == M # A is fully known. 962 963 Qoffset <- if (Quadratic) ifelse(estI.tol, 0, sum(1:Rank)) else 0 964 965 ynames <- object@misc$ynames 966 if (!length(ynames)) ynames <- object@misc$predictors.names 967 if (!length(ynames)) ynames <- object@misc$ynames 968 if (!length(ynames)) ynames <- param.names("Y", NOS) 969 lp.names <- object@misc$predictors.names 970 if (!length(lp.names)) lp.names <- NULL 971 972 dzero.vector <- rep_len(FALSE, M) 973 if (length(Dzero)) 974 dzero.vector[Dzero] <- TRUE 975 names(dzero.vector) <- ynames 976 latvar.names <- param.names("latvar", Rank, skip1 = TRUE) 977 978 979 980 td.expression <- function(Dmat, Amat, M, Dzero, Rank, bellshaped) { 981 982 Tolerance <- Darray <- m2a(Dmat, M = Rank) 983 for (ii in 1:M) 984 if (length(Dzero) && any(Dzero == ii)) { 985 Tolerance[, , ii] <- NA_real_ # Darray[,,ii] == O 986 bellshaped[ii] <- FALSE 987 } else { 988 Tolerance[, , ii] <- -0.5 * solve(Darray[, , ii]) 989 bellshaped[ii] <- 990 all(eigen(Tolerance[, , ii], symmetric = TRUE)$values > 0) 991 } 992 optimum <- matrix(NA_real_, Rank, M) 993 for (ii in 1:M) 994 if (bellshaped[ii]) 995 optimum[, ii] <- Tolerance[, , ii] %*% cbind(Amat[ii, ]) 996 997 list(optimum = optimum, 998 Tolerance = Tolerance, 999 Darray = Darray, 1000 bellshaped = bellshaped) 1001 } # td.expression 1002 1003 1004 1005 Amat <- object@extra$Amat # M x Rank 1006 Cmat <- object@extra$Cmat # p2 x Rank 1007 Dmat <- object@extra$Dmat # 1008 B1 <- object@extra$B1 # 1009 bellshaped <- rep_len(FALSE, M) 1010 1011 if (is.character(refResponse)) { 1012 refResponse <- (1:NOS)[refResponse == ynames] 1013 if (length(refResponse) != 1) 1014 stop("could not match argument 'refResponse' with any response") 1015 } 1016 ptr1 <- 1 1017 candidates <- if (length(refResponse)) refResponse else { 1018 if (length(ocontrol$Dzero)) (1:M)[-ocontrol$Dzero] else (1:M)} 1019 repeat { 1020 if (ptr1 > 0) { 1021 this.spp <- candidates[ptr1] 1022 } 1023 elts <- Dmat[this.spp,, drop = FALSE] 1024 if (length(elts) < Rank) 1025 elts <- matrix(elts, 1, Rank) 1026 Dk <- m2a(elts, M = Rank)[, , 1] # Hopefully negative-def 1027 temp400 <- eigen(Dk, symmetric = TRUE) 1028 ptr1 <- ptr1 + 1 1029 if (all(temp400$value < 0)) 1030 break 1031 if (ptr1 > length(candidates)) 1032 break 1033 } # repeat 1034 if (all(temp400$value < 0)) { 1035 temp1tol <- -0.5 * solve(Dk) 1036 dim(temp1tol) <- c(Rank,Rank) 1037 Mmat <- t(chol(temp1tol)) 1038 if (ConstrainedQO) { 1039 temp900 <- solve(t(Mmat)) 1040 Cmat <- Cmat %*% temp900 1041 Amat <- Amat %*% Mmat 1042 } 1043 if (length(Cmat)) { 1044 temp800 <- crow1C(Cmat, ocontrol$Crow1positive, amat = Amat) 1045 Cmat <- temp800$cmat 1046 Amat <- temp800$amat 1047 } 1048 1049 1050 1051 Dmat <- adjust.Dmat.expression(Mmat = Mmat, Rank = Rank, 1052 Dmat = Dmat, M = M) 1053 1054 1055 1056 retlist <- td.expression(Dmat = Dmat, Amat = Amat, M = M, 1057 Dzero = Dzero, Rank = Rank, 1058 bellshaped = bellshaped) 1059 optimum <- retlist$optimum 1060 Tolerance <- retlist$Tolerance 1061 Darray <- retlist$Darray 1062 bellshaped <- retlist$bellshaped 1063 1064 1065 1066 1067 } else { 1068 if (length(refResponse) == 1) 1069 stop("tolerance matrix specified by 'refResponse' ", 1070 "is not positive-definite") else 1071 warning("could not find any positive-definite ", 1072 "tolerance matrix") 1073 } 1074 1075 1076 if (ConstrainedQO) 1077 if (Rank > 1) { 1078 if (!length(xmat <- object@x)) 1079 stop("cannot obtain the model matrix") 1080 numat <- xmat[,ocontrol$colx2.index, drop = FALSE] %*% Cmat 1081 evnu <- eigen(var(numat), symmetric = TRUE) 1082 Mmat <- solve(t(evnu$vector)) 1083 Cmat <- Cmat %*% evnu$vector # == Cmat %*% solve(t(Mmat)) 1084 Amat <- Amat %*% Mmat 1085 temp800 <- crow1C(Cmat, ocontrol$Crow1positive, amat = Amat) 1086 Cmat <- temp800$cmat 1087 Amat <- temp800$amat 1088 1089 1090 Dmat <- adjust.Dmat.expression(Mmat = Mmat, Rank = Rank, 1091 Dmat = Dmat, M = M) 1092 1093 1094 1095 retlist <- td.expression(Dmat = Dmat, Amat = Amat, M = M, 1096 Dzero = Dzero, Rank = Rank, 1097 bellshaped = bellshaped) 1098 optimum <- retlist$optimum 1099 Tolerance <- retlist$Tolerance 1100 Darray <- retlist$Darray 1101 bellshaped <- retlist$bellshaped 1102 } # Rank > 1 1103 1104 1105 if (ConstrainedQO) 1106 if (varI.latvar) { 1107 if (!length(xmat <- object@x)) 1108 stop("cannot obtain the model matrix") 1109 numat <- xmat[,ocontrol$colx2.index, drop = FALSE] %*% Cmat 1110 sdnumat <- apply(cbind(numat), 2, sd) 1111 Mmat <- if (Rank > 1) diag(sdnumat) else matrix(sdnumat, 1, 1) 1112 Cmat <- Cmat %*% solve(t(Mmat)) 1113 Amat <- Amat %*% Mmat 1114 temp800 <- crow1C(Cmat, ocontrol$Crow1positive, amat = Amat) 1115 Cmat <- temp800$cmat 1116 Amat <- temp800$amat 1117 Cmat # Not needed 1118 1119 Dmat <- adjust.Dmat.expression(Mmat = Mmat, Rank = Rank, 1120 Dmat = Dmat, M = M) 1121 1122 1123 retlist <- td.expression(Dmat = Dmat, Amat = Amat, M = M, 1124 Dzero = Dzero, Rank = Rank, 1125 bellshaped = bellshaped) 1126 optimum <- retlist$optimum 1127 Tolerance <- retlist$Tolerance 1128 Darray <- retlist$Darray 1129 bellshaped <- retlist$bellshaped 1130 } # (varI.latvar) 1131 1132 1133 cx1i <- ocontrol$colx1.index 1134 maximum <- if (length(cx1i) == 1 && names(cx1i) == "(Intercept)") { 1135 eta.temp <- B1 1136 for (ii in 1:M) 1137 eta.temp[ii] <- eta.temp[ii] + 1138 Amat[ii, , drop = FALSE] %*% optimum[, ii, drop = FALSE] + 1139 t(optimum[, ii, drop = FALSE]) %*% 1140 Darray[,, ii, drop = TRUE] %*% optimum[, ii, drop = FALSE] 1141 mymax <- object@family@linkinv(rbind(eta.temp), 1142 extra = object@extra) 1143 c(mymax) # Convert from matrix to vector 1144 } else { 1145 5 * rep_len(NA_real_, M) # Make "numeric" 1146 } 1147 names(maximum) <- ynames 1148 1149 latvar.mat <- if (ConstrainedQO) { 1150 object@x[, ocontrol$colx2.index, drop = FALSE] %*% Cmat 1151 } else { 1152 object@latvar 1153 } 1154 1155 dimnames(Amat) <- list(lp.names, latvar.names) 1156 if (ConstrainedQO) 1157 dimnames(Cmat) <- list(names(ocontrol$colx2.index), latvar.names) 1158 if (!length(xmat <- object@x)) stop("cannot obtain the model matrix") 1159 dimnames(latvar.mat) <- list(dimnames(xmat)[[1]], latvar.names) 1160 1161 ans <- 1162 new(Class <- if (ConstrainedQO) "Coef.qrrvglm" else "Coef.uqo", 1163 A = Amat, B1 = B1, Constrained = ConstrainedQO, D = Darray, 1164 NOS = NOS, Rank = Rank, 1165 latvar = latvar.mat, 1166 latvar.order = latvar.mat, 1167 Optimum = optimum, 1168 Optimum.order = optimum, 1169 bellshaped = bellshaped, 1170 Dzero = dzero.vector, 1171 Maximum = maximum, 1172 Tolerance = Tolerance) 1173 if (ConstrainedQO) {ans@C <- Cmat} else {Cmat <- NULL} 1174 1175 for (rrr in 1:Rank) 1176 ans@Optimum.order[rrr, ] <- order(ans@Optimum[rrr, ]) 1177 for (rrr in 1:Rank) 1178 ans@latvar.order[, rrr] <- order(ans@latvar[, rrr]) 1179 1180 if (length(object@misc$estimated.dispersion) && 1181 object@misc$estimated.dispersion) { 1182 p <- length(object@coefficients) 1183 n <- object@misc$n 1184 M <- object@misc$M 1185 NOS <- if (length(object@y)) ncol(object@y) else M 1186 pstar <- if (ConstrainedQO) (p + length(Cmat)) else 1187 p + n*Rank # Adjustment; not sure about UQO 1188 adjusted.dispersion <- object@misc$dispersion * (n*M - p) / 1189 (n*M - pstar) 1190 ans@dispersion <- adjusted.dispersion 1191 } 1192 1193 if (MSratio > 1) { 1194 keepIndex <- seq(from = 1, to = M, by = MSratio) 1195 ans@Dzero <- ans@Dzero[keepIndex] 1196 ans@Optimum <- ans@Optimum[,keepIndex, drop = FALSE] 1197 ans@Tolerance <- ans@Tolerance[,,keepIndex, drop = FALSE] 1198 ans@bellshaped <- ans@bellshaped[keepIndex] 1199 names(ans@Dzero) <- ynames 1200 } else { 1201 dimnames(ans@D) <- list(latvar.names, latvar.names, ynames) 1202 } 1203 names(ans@bellshaped) <- ynames 1204 dimnames(ans@Optimum) <- list(latvar.names, ynames) 1205 dimnames(ans@Tolerance) <- list(latvar.names, latvar.names, ynames) 1206 ans 1207} # Coef.qrrvglm 1208 1209 1210 1211setClass(Class = "Coef.rrvglm", representation( 1212 "A" = "matrix", 1213 "B1" = "matrix", # This may be unassigned if p1 = 0. 1214 "C" = "matrix", 1215 "Rank" = "numeric", 1216 "colx1.index" = "numeric", 1217 "colx2.index" = "numeric", 1218 "Atilde" = "matrix")) 1219 1220 1221setClass(Class = "Coef.uqo", representation( 1222 "A" = "matrix", 1223 "B1" = "matrix", 1224 "Constrained" = "logical", 1225 "D" = "array", 1226 "NOS" = "numeric", 1227 "Rank" = "numeric", 1228 "latvar" = "matrix", 1229 "latvar.order" = "matrix", 1230 "Maximum" = "numeric", 1231 "Optimum" = "matrix", 1232 "Optimum.order" = "matrix", 1233 "bellshaped" = "logical", 1234 "dispersion" = "numeric", 1235 "Dzero" = "logical", 1236 "Tolerance" = "array")) 1237 1238 1239setClass(Class = "Coef.qrrvglm", representation( 1240 "C" = "matrix"), 1241 contains = "Coef.uqo") 1242 1243 1244 1245 1246show.Coef.qrrvglm <- function(x, ...) { 1247 1248 object <- x 1249 Rank <- object@Rank 1250 M <- nrow(object@A) 1251 NOS <- object@NOS 1252 mymat <- matrix(NA_real_, NOS, Rank) 1253 if (Rank == 1) { # || object@Diagonal 1254 for (ii in 1:NOS) { 1255 fred <- if (Rank > 1) 1256 diag(object@Tolerance[, , ii, drop = FALSE]) else 1257 object@Tolerance[, , ii] 1258 if (all(fred > 0)) 1259 mymat[ii,] <- sqrt(fred) 1260 } 1261 dimnames(mymat) <- list(dimnames(object@Tolerance)[[3]], 1262 if (Rank == 1) "latvar" else 1263 paste("Tolerance", dimnames(mymat)[[2]], 1264 sep = "")) 1265 } else { 1266 for (ii in 1:NOS) { 1267 fred <- eigen(object@Tolerance[, , ii], symmetric = TRUE) 1268 if (all(fred$value > 0)) 1269 mymat[ii, ] <- sqrt(fred$value) 1270 } 1271 dimnames(mymat) <- list(dimnames(object@Tolerance)[[3]], 1272 param.names("tol", Rank)) 1273 } 1274 1275 dimnames(object@A) <- list(dimnames(object@A)[[1]], 1276 if (Rank > 1) paste("A", dimnames(object@A)[[2]], sep = ".") else 1277 "A") 1278 1279 Maximum <- if (length(object@Maximum)) 1280 cbind(Maximum = object@Maximum) else NULL 1281 if (length(Maximum) && length(mymat) && Rank == 1) 1282 Maximum[is.na(mymat),] <- NA 1283 1284 optmat <- cbind(t(object@Optimum)) 1285 dimnames(optmat) <- list(dimnames(optmat)[[1]], 1286 if (Rank > 1) 1287 paste("Optimum", dimnames(optmat)[[2]], sep = ".") else 1288 "Optimum") 1289 if (length(optmat) && length(mymat) && Rank == 1) 1290 optmat[is.na(mymat), ] <- NA 1291 1292 if ( object@Constrained ) { 1293 cat("\nC matrix (constrained/canonical coefficients)\n") 1294 print(object@C, ...) 1295 } 1296 cat("\nB1 and A matrices\n") 1297 print(cbind(t(object@B1), 1298 A = object@A), ...) 1299 cat("\nOptimums and maximums\n") 1300 print(cbind(Optimum = optmat, 1301 Maximum), ...) 1302 if (Rank > 1) { # !object@Diagonal && Rank > 1 1303 cat("\nTolerances\n") 1304 } else { 1305 cat("\nTolerance\n") 1306 } 1307 print(mymat, ...) 1308 1309 cat("\nStandard deviation of the latent variables (site scores)\n") 1310 print(apply(cbind(object@latvar), 2, sd)) 1311 invisible(object) 1312} # show.Coef.qrrvglm 1313 1314 1315 1316 1317 1318 1319 1320setMethod("show", "Coef.qrrvglm", function(object) 1321 show.Coef.qrrvglm(object)) 1322 1323 1324 1325 1326 1327 1328 1329 1330setMethod("summary", "qrrvglm", function(object, ...) 1331 summary.qrrvglm(object, ...)) 1332 1333 1334 1335predictqrrvglm <- 1336 function(object, 1337 newdata = NULL, 1338 type = c("link", "response", "latvar", "terms"), 1339 se.fit = FALSE, 1340 deriv = 0, 1341 dispersion = NULL, 1342 extra = object@extra, 1343 varI.latvar = FALSE, refResponse = NULL, ...) { 1344 if (se.fit) 1345 stop("cannot handle se.fit == TRUE yet") 1346 if (deriv != 0) 1347 stop("derivative is not equal to 0") 1348 1349 if (mode(type) != "character" && mode(type) != "name") 1350 type <- as.character(substitute(type)) 1351 type <- match.arg(type, c("link", "response", "latvar", "terms"))[1] 1352 if (type == "latvar") 1353 stop("cannot handle type='latvar' yet") 1354 if (type == "terms") 1355 stop("cannot handle type='terms' yet") 1356 1357 M <- object@misc$M 1358 Rank <- object@control$Rank 1359 1360 na.act <- object@na.action 1361 object@na.action <- list() 1362 1363 if (!length(newdata) && 1364 type == "response" && 1365 length(object@fitted.values)) { 1366 if (length(na.act)) { 1367 return(napredict(na.act[[1]], object@fitted.values)) 1368 } else { 1369 return(object@fitted.values) 1370 } 1371 } 1372 1373 if (!length(newdata)) { 1374 X <- model.matrixvlm(object, type = "lm", ...) 1375 offset <- object@offset 1376 tt <- object@terms$terms # terms(object) 1377 if (!length(object@x)) 1378 attr(X, "assign") <- attrassignlm(X, tt) 1379 } else { 1380 if (is.smart(object) && length(object@smart.prediction)) { 1381 setup.smart("read", smart.prediction = object@smart.prediction) 1382 } 1383 1384 tt <- object@terms$terms 1385 X <- model.matrix(delete.response(tt), newdata, 1386 contrasts = if (length(object@contrasts)) 1387 object@contrasts else NULL, 1388 xlev = object@xlevels) 1389 1390 if (nrow(X) != nrow(newdata)) { 1391 as.save <- attr(X, "assign") 1392 X <- X[rep_len(1, nrow(newdata)),, drop = FALSE] 1393 dimnames(X) <- list(dimnames(newdata)[[1]], "(Intercept)") 1394 attr(X, "assign") <- as.save # Restored 1395 } 1396 1397 offset <- if (!is.null(off.num<-attr(tt,"offset"))) { 1398 eval(attr(tt,"variables")[[off.num+1]], newdata) 1399 } else if (!is.null(object@offset)) 1400 eval(object@call$offset, newdata) 1401 1402 if (any(c(offset) != 0)) 1403 stop("currently cannot handle nonzero offsets") 1404 1405 if (is.smart(object) && length(object@smart.prediction)) { 1406 wrapup.smart() 1407 } 1408 1409 attr(X, "assign") <- attrassigndefault(X, tt) 1410 } 1411 1412 ocontrol <- object@control 1413 1414 Rank <- ocontrol$Rank 1415 NOS <- ncol(object@y) 1416 sppnames <- dimnames(object@y)[[2]] 1417 modelno <- ocontrol$modelno # 1, 2, 3, 5 or 0 1418 M <- if (any(slotNames(object) == "predictors") && 1419 is.matrix(object@predictors)) 1420 ncol(object@predictors) else 1421 object@misc$M 1422 MSratio <- M / NOS # 1st value is g(mean)=quadratic form in latvar 1423 if (MSratio != 1) stop("can only handle MSratio == 1 for now") 1424 1425 1426 if (length(newdata)) { 1427 Coefs <- Coef(object, varI.latvar = varI.latvar, 1428 refResponse = refResponse) 1429 X1mat <- X[, ocontrol$colx1.index, drop = FALSE] 1430 X2mat <- X[, ocontrol$colx2.index, drop = FALSE] 1431 latvarmat <- as.matrix(X2mat %*% Coefs@C) # n x Rank 1432 1433 etamat <- as.matrix(X1mat %*% Coefs@B1 + latvarmat %*% t(Coefs@A)) 1434 which.species <- 1:NOS # Do it all for all species 1435 for (sppno in seq_along(which.species)) { 1436 thisSpecies <- which.species[sppno] 1437 Dmat <- matrix(Coefs@D[,,thisSpecies], Rank, Rank) 1438 etamat[, thisSpecies] <- etamat[, thisSpecies] + 1439 mux34(latvarmat, Dmat, symmetric = TRUE) 1440 } 1441 } else { 1442 etamat <- object@predictors 1443 } 1444 1445 pred <- 1446 switch(type, 1447 response = { 1448 fv <- if (length(newdata)) 1449 object@family@linkinv(etamat, extra) else 1450 fitted(object) 1451 if (M > 1 && is.matrix(fv)) { 1452 dimnames(fv) <- list(dimnames(fv)[[1]], 1453 dimnames(object@fitted.values)[[2]]) 1454 } 1455 fv 1456 }, 1457 link = etamat, 1458 latvar = stop("failure here"), 1459 terms = stop("failure here")) 1460 1461 if (!length(newdata) && length(na.act)) { 1462 if (se.fit) { 1463 pred$fitted.values <- napredict(na.act[[1]], pred$fitted.values) 1464 pred$se.fit <- napredict(na.act[[1]], pred$se.fit) 1465 } else { 1466 pred <- napredict(na.act[[1]], pred) 1467 } 1468 } 1469 pred 1470} # predictqrrvglm 1471 1472 1473setMethod("predict", "qrrvglm", function(object, ...) 1474 predictqrrvglm(object, ...)) 1475 1476 1477 1478 1479coefqrrvglm <- 1480 function(object, matrix.out = FALSE, 1481 label = TRUE) { 1482 if (matrix.out) 1483 stop("currently cannot handle matrix.out = TRUE") 1484 1485 cobj <- coefvlm(object, matrix.out = matrix.out, label = label) 1486 1487 1488 M <- npred(object) 1489 M1 <- npred(object, type = "one.response") 1490 NOS <- M / M1 1491 Rank <- Rank(object) 1492 colx1.index <- object@control$colx1.index 1493 colx2.index <- object@control$colx2.index 1494 etol <- object@control$eq.tolerances 1495 Itol <- object@control$I.tolerances 1496 names.A <- param.names("A", M * Rank, skip1 = FALSE) 1497 if (Itol) 1498 names.D <- NULL # Because it was estimated by offsets 1499 if (etol && !Itol) 1500 names.D <- param.names("D", Rank * (Rank + 1) / 2, skip1 = FALSE) 1501 if (!etol) 1502 names.D <- param.names("D", 1503 NOS * Rank * (Rank + 1) / 2, skip1 = FALSE) 1504 names.B1 <- param.names("x1.", NOS * length(colx1.index), skip1 = FALSE) 1505 if (length(temp <- c(names.A, names.D, names.B1)) == length(cobj)) 1506 names(cobj) <- temp 1507 1508 cobj 1509} # coefqrrvglm 1510 1511 1512 1513 1514 1515 1516 1517residualsqrrvglm <- 1518 function(object, 1519 type = c("deviance", "pearson", "working", "response", "ldot"), 1520 matrix.arg = TRUE) { 1521 stop("this function has not been written yet") 1522} 1523 1524 1525setMethod("residuals", "qrrvglm", 1526 function(object, ...) 1527 residualsqrrvglm(object, ...)) 1528 1529 1530 1531 1532 1533show.rrvglm <- function(x, ...) { 1534 if (!is.null(cl <- x@call)) { 1535 cat("Call:\n") 1536 dput(cl) 1537 } 1538 vecOfBetas <- x@coefficients 1539 if (any(nas <- is.na(vecOfBetas))) { 1540 if (is.null(names(vecOfBetas))) 1541 names(vecOfBetas) <- param.names("b", length(vecOfBetas)) 1542 cat("\nCoefficients: (", sum(nas), 1543 " not defined because of singularities)\n", sep = "") 1544 } else 1545 cat("\nCoefficients:\n") 1546 print.default(vecOfBetas, ...) # used to be print() 1547 1548 if (FALSE) { 1549 Rank <- x@Rank 1550 if (!length(Rank)) 1551 Rank <- sum(!nas) 1552 } 1553 1554 if (FALSE) { 1555 nobs <- if (length(x@df.total)) x@df.total else length(x@residuals) 1556 rdf <- x@df.residual 1557 if (!length(rdf)) 1558 rdf <- nobs - Rank 1559 } 1560 cat("\n") 1561 1562 if (length(deviance(x))) 1563 cat("Residual deviance:", format(deviance(x)), "\n") 1564 if (length(vll <- logLik.vlm(x))) 1565 cat("Log-likelihood:", format(vll), "\n") 1566 1567 if (length(x@criterion)) { 1568 ncrit <- names(x@criterion) 1569 for (iii in ncrit) 1570 if (iii != "loglikelihood" && 1571 iii != "deviance") 1572 cat(paste(iii, ":", sep = ""), 1573 format(x@criterion[[iii]]), "\n") 1574 } 1575 1576 invisible(x) 1577} # show.rrvglm 1578 1579 1580 1581 1582 1583 1584setMethod("show", "rrvglm", function(object) show.rrvglm(object)) 1585 1586 1587 1588 1589 1590 1591 1592 1593summary.rrvglm <- 1594 function(object, correlation = FALSE, 1595 dispersion = NULL, digits = NULL, 1596 numerical = TRUE, 1597 h.step = 0.0001, 1598 kill.all = FALSE, omit13 = FALSE, 1599 fixA = FALSE, 1600 presid = TRUE, 1601 signif.stars = getOption("show.signif.stars"), 1602 nopredictors = FALSE, ...) { 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 if (!is.Numeric(h.step, length.arg = 1) || 1613 abs(h.step) > 1) 1614 stop("bad input for 'h.step'") 1615 1616 if (!object@control$Corner) 1617 stop("this function works with corner constraints only") 1618 1619 if (is.null(dispersion)) 1620 dispersion <- object@misc$dispersion 1621 1622 newobject <- as(object, "vglm") 1623 1624 1625 stuff <- summaryvglm(newobject, 1626 correlation = correlation, 1627 dispersion = dispersion, 1628 presid = presid) 1629 1630 answer <- 1631 new(Class = "summary.rrvglm", 1632 object, 1633 call = stuff@call, 1634 coef3 = stuff@coef3, 1635 cov.unscaled = stuff@cov.unscaled, 1636 correlation = stuff@correlation, 1637 df = stuff@df, 1638 sigma = stuff@sigma) 1639 1640 1641 if (is.numeric(stuff@dispersion)) 1642 slot(answer, "dispersion") <- stuff@dispersion 1643 1644 if (presid && length(stuff@pearson.resid)) 1645 slot(answer, "pearson.resid") <- stuff@pearson.resid 1646 1647 1648 1649 tmp5 <- get.rrvglm.se1(object, omit13 = omit13, 1650 numerical = numerical, h.step = h.step, 1651 kill.all = kill.all, fixA = fixA, ...) 1652 if (any(diag(tmp5$cov.unscaled) <= 0) || 1653 any(eigen(tmp5$cov.unscaled, symmetric = TRUE)$value <= 0)) { 1654 warning("cov.unscaled is not positive definite") 1655 } 1656 1657 answer@cov.unscaled <- tmp5$cov.unscaled 1658 1659 od <- if (is.numeric(object@misc$disper)) 1660 object@misc$disper else 1661 object@misc$default.disper 1662 if (is.numeric(dispersion)) { 1663 if (is.numeric(od) && dispersion != od) 1664 warning("dispersion != object@misc$dispersion; ", 1665 "using the former") 1666 } else { 1667 dispersion <- if (is.numeric(od)) od else 1 1668 } 1669 1670 tmp8 <- object@misc$M - object@control$Rank - 1671 length(object@control$str0) 1672 answer@df[1] <- answer@df[1] + tmp8 * object@control$Rank 1673 answer@df[2] <- answer@df[2] - tmp8 * object@control$Rank 1674 if (dispersion == 0) { 1675 dispersion <- tmp5$ResSS / answer@df[2] # Estimate 1676 } 1677 1678 answer@coef3 <- get.rrvglm.se2(answer@cov.unscaled, 1679 dispersion = dispersion, 1680 coefficients = tmp5$coefficients) 1681 1682 answer@dispersion <- dispersion 1683 answer@sigma <- dispersion^0.5 1684 1685 1686 answer@misc$signif.stars <- signif.stars # 20160629 1687 answer@misc$nopredictors <- nopredictors # 20150925 1688 1689 answer 1690} # summary.rrvglm 1691 1692 1693 1694 1695 1696 1697get.rrvglm.se1 <- 1698 function(fit, omit13 = FALSE, kill.all = FALSE, 1699 numerical = TRUE, 1700 fixA = FALSE, h.step = 0.0001, 1701 trace.arg = FALSE, ...) { 1702 1703 1704 1705 1706 if (length(fit@control$Nested) && fit@control$Nested) 1707 stop("sorry, cannot handle nested models yet") 1708 1709 str0 <- fit@control$str0 1710 1711 1712 if (!length(fit@x)) 1713 stop("fix@x is empty. Run rrvglm(... , x = TRUE)") 1714 1715 colx1.index <- fit@control$colx1.index # May be NULL 1716 colx2.index <- fit@control$colx2.index 1717 Hlist <- fit@constraints 1718 ncolHlist <- unlist(lapply(Hlist, ncol)) 1719 1720 p1 <- length(colx1.index) # May be 0 1721 p2 <- length(colx2.index) 1722 1723 Rank <- fit@control$Rank # fit@misc$Nested.Rank 1724 1725 Amat <- fit@constraints[[colx2.index[1]]] 1726 B1mat <- if (p1) 1727 coefvlm(fit, matrix.out = TRUE)[colx1.index, , drop = FALSE] else 1728 NULL 1729 C.try <- coefvlm(fit, matrix.out= TRUE)[colx2.index, , drop = FALSE] 1730 Cmat <- C.try %*% Amat %*% solve(t(Amat) %*% Amat) 1731 1732 x1mat <- if (p1) fit@x[, colx1.index, drop = FALSE] else NULL 1733 x2mat <- fit@x[, colx2.index, drop = FALSE] 1734 1735 wz <- weights(fit, type = "work") # old: wweights(fit) #fit@weights 1736 if (!length(wz)) 1737 stop("cannot get fit@weights") 1738 1739 M <- fit@misc$M 1740 n <- fit@misc$n 1741 Index.corner <- fit@control$Index.corner # used to be (1:Rank); 1742 zmat <- fit@predictors + fit@residuals 1743 theta <- c(Amat[-c(Index.corner,str0), ]) 1744 if (fit@control$checkwz) 1745 wz <- checkwz(wz, M = M, trace = trace, 1746 wzepsilon = fit@control$wzepsilon) 1747 U <- vchol(wz, M = M, n = n, silent= TRUE) 1748 1749 delct.da <- if (numerical) { 1750 num.deriv.rrr(fit, M = M, r = Rank, 1751 x1mat = x1mat, x2mat = x2mat, p2 = p2, 1752 Index.corner, Aimat = Amat, 1753 B1mat = B1mat, Cimat = Cmat, 1754 h.step = h.step, 1755 colx2.index = colx2.index, 1756 xij = fit@control$xij, 1757 str0 = str0) 1758 } else { 1759 dctda.fast.only(theta = theta, wz = wz, 1760 U = U, zmat, 1761 M = M, r = Rank, x1mat = x1mat, 1762 x2mat = x2mat, p2 = p2, 1763 Index.corner, Aimat = Amat, 1764 B1mat = B1mat, Cimat = Cmat, 1765 xij = fit@control$xij, 1766 str0 = str0) 1767 } 1768 1769 1770 1771 1772 newobject <- as(fit, "vglm") 1773 1774 1775 1776 1777 sfit2233 <- summaryvglm(newobject) 1778 d8 <- dimnames(sfit2233@cov.unscaled)[[1]] 1779 cov2233 <- solve(sfit2233@cov.unscaled) # Includes any intercepts 1780 dimnames(cov2233) <- list(d8, d8) 1781 1782 log.vec33 <- NULL 1783 nassign <- names(fit@constraints) 1784 choose.from <- varassign(fit@constraints, nassign) 1785 for (ii in nassign) 1786 if (any(ii == names(colx2.index))) { 1787 log.vec33 <- c(log.vec33, choose.from[[ii]]) 1788 } 1789 cov33 <- cov2233[ log.vec33, log.vec33, drop = FALSE] # r*p2 by r*p2 1790 cov23 <- cov2233[-log.vec33, log.vec33, drop = FALSE] 1791 cov22 <- cov2233[-log.vec33,-log.vec33, drop = FALSE] 1792 1793 1794 latvar.mat <- x2mat %*% Cmat 1795 offs <- matrix(0, n, M) # The "0" handles str0's 1796 offs[, Index.corner] <- latvar.mat 1797 if (M == (Rank + length(str0))) 1798 stop("cannot handle full-rank models yet") 1799 cm <- matrix(0, M, M - Rank - length(str0)) 1800 cm[-c(Index.corner, str0), ] <- diag(M - Rank - length(str0)) 1801 1802 Hlist <- vector("list", length(colx1.index)+1) 1803 names(Hlist) <- c(names(colx1.index), "I(latvar.mat)") 1804 for (ii in names(colx1.index)) 1805 Hlist[[ii]] <- fit@constraints[[ii]] 1806 Hlist[["I(latvar.mat)"]] <- cm 1807 1808 1809 if (p1) { 1810 ooo <- fit@assign 1811 bb <- NULL 1812 for (ii in seq_along(ooo)) { 1813 if (any(ooo[[ii]][1] == colx1.index)) 1814 bb <- c(bb, names(ooo)[ii]) 1815 } 1816 1817 has.intercept <- any(bb == "(Intercept)") 1818 bb[bb == "(Intercept)"] <- "1" 1819 if (p1 > 1) 1820 bb <- paste(bb, collapse = "+") 1821 if (has.intercept) { 1822 bb <- paste("zmat - offs ~ ", bb, " + I(latvar.mat)", 1823 collapse = " ") 1824 } else { 1825 bb <- paste("zmat - offs ~ -1 + ", bb, " + I(latvar.mat)", 1826 collapse = " ") 1827 } 1828 bb <- as.formula(bb) 1829 } else { 1830 bb <- as.formula("zmat - offs ~ -1 + I(latvar.mat)") 1831 } 1832 1833 1834 if (fit@misc$dataname == "list") { 1835 dspec <- FALSE 1836 } else { 1837 mytext1 <- "exists(x = fit@misc$dataname, envir = VGAMenv)" 1838 myexp1 <- parse(text = mytext1) 1839 is.there <- eval(myexp1) 1840 bbdata <- if (is.there) 1841 get(fit@misc$dataname, envir = VGAMenv) else 1842 get(fit@misc$dataname) 1843 dspec <- TRUE 1844 } 1845 1846 1847 fit1122 <- if (dspec) 1848 vlm(bb, 1849 constraints = Hlist, criterion = "d", weights = wz, 1850 data = bbdata, 1851 save.weights = TRUE, smart = FALSE, trace = trace.arg, 1852 x.arg = TRUE) else 1853 vlm(bb, 1854 constraints = Hlist, criterion = "d", weights = wz, 1855 save.weights = TRUE, smart = FALSE, trace = trace.arg, 1856 x.arg = TRUE) 1857 1858 1859 1860 sfit1122 <- summaryvlm(fit1122) 1861 d8 <- dimnames(sfit1122@cov.unscaled)[[1]] 1862 cov1122 <- solve(sfit1122@cov.unscaled) 1863 dimnames(cov1122) <- list(d8, d8) 1864 1865 lcs <- length(coefvlm(sfit1122)) 1866 log.vec11 <- (lcs-(M-Rank-length(str0))*Rank+1):lcs 1867 cov11 <- cov1122[log.vec11, log.vec11, drop = FALSE] 1868 cov12 <- cov1122[ log.vec11, -log.vec11, drop = FALSE] 1869 cov22 <- cov1122[-log.vec11, -log.vec11, drop = FALSE] 1870 cov13 <- delct.da %*% cov33 1871 1872 1873 if (omit13) 1874 cov13 <- cov13 * 0 # zero it 1875 1876 if (kill.all) { 1877 cov13 <- cov13 * 0 # zero it 1878 if (fixA) { 1879 cov12 <- cov12 * 0 # zero it 1880 } else { 1881 cov23 <- cov23 * 0 # zero it 1882 } 1883 } 1884 1885 cov13 <- -cov13 # Richards (1961) 1886 1887 if (fixA) { 1888 cov.unscaled <- rbind(cbind(cov1122, rbind(cov13, cov23)), 1889 cbind(t(cov13), t(cov23), cov33)) 1890 } else { 1891 cov.unscaled <- rbind(cbind(cov11, cov12, cov13), 1892 cbind(rbind(t(cov12), t(cov13)), cov2233)) 1893 } 1894 1895 ans <- solve(cov.unscaled) 1896 1897 acoefs <- c(fit1122@coefficients[log.vec11], fit@coefficients) 1898 dimnames(ans) <- list(names(acoefs), names(acoefs)) 1899 list(cov.unscaled = ans, 1900 coefficients = acoefs, 1901 ResSS = sfit1122@ResSS) 1902} # get.rrvglm.se1 1903 1904 1905 1906 1907 1908 1909get.rrvglm.se2 <- function(cov.unscaled, dispersion = 1, coefficients) { 1910 1911 d8 <- dimnames(cov.unscaled)[[1]] 1912 ans <- matrix(coefficients, length(coefficients), 4) 1913 ans[, 2] <- sqrt(dispersion) * sqrt(diag(cov.unscaled)) 1914 ans[, 3] <- ans[, 1] / ans[, 2] 1915 ans[, 4] <- pnorm(-abs(ans[, 3])) 1916 dimnames(ans) <- 1917 list(d8, c("Estimate", "Std. Error", "z value", "Pr(>|z|)")) 1918 ans 1919} # get.rrvglm.se2 1920 1921 1922 1923 1924 1925 1926num.deriv.rrr <- function(fit, M, r, x1mat, x2mat, 1927 p2, Index.corner, Aimat, B1mat, Cimat, 1928 h.step = 0.0001, colx2.index, 1929 xij = NULL, str0 = NULL) { 1930 1931 1932 nn <- nrow(x2mat) 1933 if (nrow(Cimat) != p2 || ncol(Cimat) != r) 1934 stop("'Cimat' wrong shape") 1935 1936 dct.da <- matrix(NA_real_, (M-r-length(str0))*r, r*p2) 1937 1938 if ((length(Index.corner) + length(str0)) == M) 1939 stop("cannot handle full rank models yet") 1940 cbindex <- (1:M)[-c(Index.corner, str0)] 1941 1942 ptr <- 1 1943 for (sss in 1:r) 1944 for (tt in cbindex) { 1945 small.Hlist <- vector("list", p2) 1946 pAmat <- Aimat 1947 pAmat[tt,sss] <- pAmat[tt,sss] + h.step # Perturb it 1948 for (ii in 1:p2) 1949 small.Hlist[[ii]] <- pAmat 1950 1951 offset <- if (length(fit@offset)) fit@offset else 0 1952 if (all(offset == 0)) 1953 offset <- 0 1954 neweta <- x2mat %*% Cimat %*% t(pAmat) 1955 if (is.numeric(x1mat)) 1956 neweta <- neweta + x1mat %*% B1mat 1957 fit@predictors <- neweta 1958 1959 1960 newmu <- fit@family@linkinv(neweta, fit@extra) 1961 fit@fitted.values <- as.matrix(newmu) # 20100909 1962 1963 fred <- weights(fit, type = "w", deriv = TRUE, ignore.slot = TRUE) 1964 if (!length(fred)) 1965 stop("cannot get @weights and @deriv from object") 1966 wz <- fred$weights 1967 deriv.mu <- fred$deriv 1968 1969 U <- vchol(wz, M = M, n = nn, silent = TRUE) 1970 tvfor <- vforsub(U, as.matrix(deriv.mu), M = M, n = nn) 1971 newzmat <- neweta + vbacksub(U, tvfor, M = M, n = nn) - offset 1972 if (is.numeric(x1mat)) 1973 newzmat <- newzmat - x1mat %*% B1mat 1974 1975 newfit <- vlm.wfit(xmat = x2mat, zmat = newzmat, 1976 Hlist = small.Hlist, U = U, 1977 matrix.out = FALSE, is.vlmX = FALSE, 1978 ResSS = TRUE, qr = FALSE, x.ret = FALSE, 1979 offset = NULL, xij = xij) 1980 dct.da[ptr, ] <- (newfit$coef - t(Cimat)) / h.step 1981 ptr <- ptr + 1 1982 } 1983 1984 dct.da 1985} # num.deriv.rrr 1986 1987 1988 1989 1990 1991dctda.fast.only <- function(theta, wz, U, zmat, M, r, x1mat, x2mat, 1992 p2, Index.corner, Aimat, B1mat, Cimat, 1993 xij = NULL, 1994 str0 = NULL) { 1995 1996 1997 if (length(str0)) 1998 stop("cannot handle 'str0' in dctda.fast.only()") 1999 2000 nn <- nrow(x2mat) 2001 if (nrow(Cimat) != p2 || ncol(Cimat) != r) 2002 stop("Cimat wrong shape") 2003 2004 fred <- kronecker(matrix(1, 1,r), x2mat) 2005 fred <- kronecker(fred, matrix(1,M, 1)) 2006 barney <- kronecker(Aimat, matrix(1, 1,p2)) 2007 barney <- kronecker(matrix(1, nn, 1), barney) 2008 2009 temp <- array(t(barney*fred), c(p2*r, M, nn)) 2010 temp <- aperm(temp, c(2, 1, 3)) # M by p2*r by nn 2011 temp <- mux5(wz, temp, M = M, matrix.arg= TRUE) 2012 temp <- m2a(temp, M = p2 * r) # Note M != M here! 2013 G <- solve(rowSums(temp, dims = 2)) # p2*r by p2*r 2014 2015 dc.da <- array(NA_real_, c(p2, r, M, r)) 2016 if (length(Index.corner) == M) 2017 stop("cannot handle full rank models yet") 2018 cbindex <- (1:M)[-Index.corner] # complement of Index.corner 2019 resid2 <- if (length(x1mat)) 2020 mux22(t(wz), zmat - x1mat %*% B1mat, M = M, 2021 upper = FALSE, as.matrix = TRUE) else 2022 mux22(t(wz), zmat , M = M, 2023 upper = FALSE, as.matrix = TRUE) 2024 2025 for (sss in 1:r) 2026 for (ttt in cbindex) { 2027 fred <- t(x2mat) * 2028 matrix(resid2[, ttt], p2, nn, byrow = TRUE) # p2 * nn 2029 temp2 <- kronecker(I.col(sss, r), rowSums(fred)) 2030 for (kkk in 1:r) { 2031 Wiak <- mux22(t(wz), matrix(Aimat[,kkk], nn, M, byrow = TRUE), 2032 M = M, upper = FALSE, 2033 as.matrix = TRUE) # nn * M 2034 wxx <- Wiak[,ttt] * x2mat 2035 blocki <- t(x2mat) %*% wxx 2036 temp4a <- blocki %*% Cimat[,kkk] 2037 if (kkk == 1) { 2038 temp4b <- blocki %*% Cimat[,sss] 2039 } 2040 temp2 <- temp2 - kronecker(I.col(sss, r), temp4a) - 2041 kronecker(I.col(kkk, r), temp4b) 2042 } 2043 dc.da[,,ttt,sss] <- G %*% temp2 2044 } 2045 ans1 <- dc.da[,,cbindex,, drop = FALSE] # p2 x r x (M-r) x r 2046 ans1 <- aperm(ans1, c(2, 1, 3, 4)) # r x p2 x (M-r) x r 2047 2048 ans1 <- matrix(c(ans1), r*p2, (M-r)*r) 2049 ans1 <- t(ans1) 2050 ans1 2051} # dcda.fast.only 2052 2053 2054 2055 2056 2057dcda.fast <- function(theta, wz, U, z, M, r, xmat, pp, Index.corner, 2058 intercept = TRUE, xij = NULL) { 2059 2060 2061 2062 nn <- nrow(xmat) 2063 2064 Aimat <- matrix(NA_real_, M, r) 2065 Aimat[Index.corner,] <- diag(r) 2066 Aimat[-Index.corner,] <- theta # [-(1:M)] 2067 2068 if (intercept) { 2069 Hlist <- vector("list", pp+1) 2070 Hlist[[1]] <- diag(M) 2071 for (ii in 2:(pp+1)) 2072 Hlist[[ii]] <- Aimat 2073 } else { 2074 Hlist <- vector("list", pp) 2075 for (ii in 1:pp) 2076 Hlist[[ii]] <- Aimat 2077 } 2078 2079 coeffs <- vlm.wfit(xmat = xmat, z, Hlist, U = U, matrix.out = TRUE, 2080 xij = xij)$mat.coef 2081 c3 <- coeffs <- t(coeffs) # transpose to make M x (pp+1) 2082 2083 2084 int.vec <- if (intercept) c3[, 1] else 0 # \boldeta_0 2085 Cimat <- if (intercept) t(c3[Index.corner,-1, drop = FALSE]) else 2086 t(c3[Index.corner,, drop = FALSE]) 2087 if (nrow(Cimat)!=pp || ncol(Cimat)!=r) 2088 stop("Cimat wrong shape") 2089 2090 fred <- kronecker(matrix(1, 1,r), 2091 if (intercept) xmat[,-1, drop = FALSE] else xmat) 2092 fred <- kronecker(fred, matrix(1,M, 1)) 2093 barney <- kronecker(Aimat, matrix(1, 1,pp)) 2094 barney <- kronecker(matrix(1, nn, 1), barney) 2095 2096 temp <- array(t(barney*fred), c(r*pp,M,nn)) 2097 temp <- aperm(temp, c(2, 1, 3)) 2098 temp <- mux5(wz, temp, M = M, matrix.arg = TRUE) 2099 temp <- m2a(temp, M = r * pp) # Note M != M here! 2100 G <- solve(rowSums(temp, dims = 2)) 2101 2102 dc.da <- array(NA_real_, c(pp, r, M, r)) 2103 cbindex <- (1:M)[-Index.corner] 2104 resid2 <- mux22(t(wz), 2105 z - matrix(int.vec, nn, M, byrow = TRUE), M = M, 2106 upper = FALSE, as.matrix = TRUE) # mat = TRUE, 2107 2108 for (s in 1:r) 2109 for (tt in cbindex) { 2110 fred <- (if (intercept) t(xmat[, -1, drop = FALSE]) else 2111 t(xmat)) * matrix(resid2[, tt], pp, nn, byrow = TRUE) 2112 temp2 <- kronecker(I.col(s, r), rowSums(fred)) 2113 2114 temp4 <- rep_len(0, pp) 2115 for (k in 1:r) { 2116 Wiak <- mux22(t(wz), 2117 matrix(Aimat[, k], nn, M, byrow = TRUE), 2118 M = M, upper = FALSE, as.matrix = TRUE) 2119 wxx <- Wiak[,tt] * (if (intercept) 2120 xmat[, -1, drop = FALSE] else 2121 xmat) 2122 blocki <- (if (intercept) 2123 t(xmat[, -1, drop = FALSE]) else 2124 t(xmat)) %*% wxx 2125 temp4 <- temp4 + blocki %*% Cimat[, k] 2126 } 2127 dc.da[,,tt,s] <- G %*% (temp2 - 2 * kronecker(I.col(s, r), temp4)) 2128 } 2129 ans1 <- dc.da[,,cbindex,, drop = FALSE] # pp x r x (M-r) x r 2130 ans1 <- aperm(ans1, c(2, 1, 3, 4)) # r x pp x (M-r) x r 2131 2132 ans1 <- matrix(c(ans1), (M-r)*r, r*pp, byrow = TRUE) 2133 2134 2135 detastar.da <- array(0,c(M,r,r,nn)) 2136 for (s in 1:r) 2137 for (j in 1:r) { 2138 t1 <- t(dc.da[,j,,s]) 2139 t1 <- matrix(t1, M, pp) 2140 detastar.da[,j,s,] <- t1 %*% (if (intercept) 2141 t(xmat[,-1, drop = FALSE]) else t(xmat)) 2142 } 2143 2144 etastar <- (if (intercept) xmat[,-1, drop = FALSE] else xmat) %*% Cimat 2145 eta <- matrix(int.vec, nn, M, byrow = TRUE) + etastar %*% t(Aimat) 2146 2147 sumWinv <- solve((m2a(t(colSums(wz)), M = M))[, , 1]) 2148 2149 deta0.da <- array(0,c(M,M,r)) 2150 AtWi <- kronecker(matrix(1, nn, 1), Aimat) 2151 AtWi <- mux111(t(wz), AtWi, M = M, upper= FALSE) # matrix.arg= TRUE, 2152 AtWi <- array(t(AtWi), c(r, M, nn)) 2153 for (ss in 1:r) { 2154 temp90 <- (m2a(t(colSums(etastar[, ss]*wz)), M = M))[, , 1] # MxM 2155 temp92 <- array(detastar.da[,,ss,], c(M, r, nn)) 2156 temp93 <- mux7(temp92, AtWi) 2157 temp91 <- rowSums(temp93, dims = 2) # M x M 2158 deta0.da[,,ss] <- -(temp90 + temp91) %*% sumWinv 2159 } 2160 ans2 <- deta0.da[-(1:r), , , drop = FALSE] # (M-r) x M x r 2161 ans2 <- aperm(ans2, c(1, 3, 2)) # (M-r) x r x M 2162 ans2 <- matrix(c(ans2), (M-r)*r, M) 2163 2164 list(dc.da = ans1, dint.da = ans2) 2165} # dcda.fast 2166 2167 2168 2169 2170 2171 2172rrr.deriv.ResSS <- function(theta, wz, U, z, M, r, xmat, 2173 pp, Index.corner, intercept = TRUE, 2174 xij = NULL) { 2175 2176 Amat <- matrix(NA_real_, M, r) 2177 Amat[Index.corner,] <- diag(r) 2178 Amat[-Index.corner,] <- theta # [-(1:M)] 2179 2180 if (intercept) { 2181 Hlist <- vector("list", pp+1) 2182 Hlist[[1]] <- diag(M) 2183 for (ii in 2:(pp+1)) 2184 Hlist[[ii]] <- Amat 2185 } else { 2186 Hlist <- vector("list", pp) 2187 for (ii in 1:pp) 2188 Hlist[[ii]] <- Amat 2189 } 2190 2191 vlm.wfit(xmat = xmat, z, Hlist, U = U, matrix.out = FALSE, 2192 ResSS = TRUE, xij = xij)$ResSS 2193} # rrr.deriv.ResSS 2194 2195 2196 2197 2198rrr.deriv.gradient.fast <- function(theta, wz, U, z, M, r, xmat, 2199 pp, Index.corner, 2200 intercept = TRUE) { 2201 2202 2203 2204 2205 nn <- nrow(xmat) 2206 2207 Aimat <- matrix(NA_real_, M, r) 2208 Aimat[Index.corner,] <- diag(r) 2209 Aimat[-Index.corner,] <- theta # [-(1:M)] 2210 2211 if (intercept) { 2212 Hlist <- vector("list", pp+1) 2213 Hlist[[1]] <- diag(M) 2214 for (i in 2:(pp+1)) 2215 Hlist[[i]] <- Aimat 2216 } else { 2217 Hlist <- vector("list", pp) 2218 for (i in 1:(pp)) 2219 Hlist[[i]] <- Aimat 2220 } 2221 2222 coeffs <- vlm.wfit(xmat, z, Hlist, U = U, matrix.out= TRUE, 2223 xij = NULL)$mat.coef 2224 c3 <- coeffs <- t(coeffs) # transpose to make M x (pp+1) 2225 2226 2227 int.vec <- if (intercept) c3[, 1] else 0 # \boldeta_0 2228 Cimat <- if (intercept) t(c3[Index.corner, -1, drop = FALSE]) else 2229 t(c3[Index.corner,, drop = FALSE]) 2230 if (nrow(Cimat) != pp || ncol(Cimat) != r) 2231 stop("Cimat wrong shape") 2232 2233 fred <- kronecker(matrix(1, 1,r), 2234 if (intercept) xmat[, -1, drop = FALSE] else xmat) 2235 fred <- kronecker(fred, matrix(1, M, 1)) 2236 barney <- kronecker(Aimat, matrix(1, 1, pp)) 2237 barney <- kronecker(matrix(1, nn, 1), barney) 2238 2239 temp <- array(t(barney*fred), c(r*pp, M, nn)) 2240 temp <- aperm(temp, c(2, 1, 3)) 2241 temp <- mux5(wz, temp, M = M, matrix.arg = TRUE) 2242 temp <- m2a(temp, M = r * pp) # Note M != M here! 2243 G <- solve(rowSums(temp, dims = 2)) 2244 2245 dc.da <- array(NA_real_, c(pp, r, r, M)) 2246 cbindex <- (1:M)[-Index.corner] 2247 resid2 <- mux22(t(wz), z - matrix(int.vec, nn, M, byrow = TRUE), 2248 M = M, 2249 upper = FALSE, as.matrix = TRUE) 2250 2251 for (s in 1:r) 2252 for (tt in cbindex) { 2253 fred <- (if (intercept) t(xmat[, -1, drop = FALSE]) else 2254 t(xmat)) * matrix(resid2[, tt], pp, nn, byrow = TRUE) 2255 temp2 <- kronecker(I.col(s, r), rowSums(fred)) 2256 2257 temp4 <- rep_len(0, pp) 2258 for (k in 1:r) { 2259 Wiak <- mux22(t(wz), 2260 matrix(Aimat[, k], nn, M, byrow = TRUE), 2261 M = M, upper = FALSE, as.matrix = TRUE) 2262 wxx <- Wiak[,tt] * (if (intercept) 2263 xmat[, -1, drop = FALSE] else xmat) 2264 blocki <- (if (intercept) t(xmat[, -1, drop = FALSE]) else 2265 t(xmat)) %*% wxx 2266 temp4 <- temp4 + blocki %*% Cimat[, k] 2267 } 2268 dc.da[,,s,tt] <- G %*% (temp2 - 2 * kronecker(I.col(s, r), temp4)) 2269 } 2270 2271 detastar.da <- array(0,c(M,r,r,nn)) 2272 for (s in 1:r) 2273 for (j in 1:r) { 2274 t1 <- t(dc.da[,j,s,]) 2275 t1 <- matrix(t1, M, pp) 2276 detastar.da[,j,s,] <- t1 %*% (if (intercept) 2277 t(xmat[, -1, drop = FALSE]) else t(xmat)) 2278 } 2279 2280 etastar <- (if (intercept) xmat[, -1, drop = FALSE] else 2281 xmat) %*% Cimat 2282 eta <- matrix(int.vec, nn, M, byrow = TRUE) + etastar %*% t(Aimat) 2283 2284 sumWinv <- solve((m2a(t(colSums(wz)), M = M))[, , 1]) 2285 2286 deta0.da <- array(0, c(M, M, r)) 2287 2288 AtWi <- kronecker(matrix(1, nn, 1), Aimat) 2289 AtWi <- mux111(t(wz), AtWi, M = M, upper = FALSE) # matrix.arg= TRUE, 2290 AtWi <- array(t(AtWi), c(r, M, nn)) 2291 2292 for (ss in 1:r) { 2293 temp90 <- (m2a(t(colSums(etastar[, ss] * wz)), M = M))[, , 1] 2294 temp92 <- array(detastar.da[, , ss, ], c(M, r, nn)) 2295 temp93 <- mux7(temp92,AtWi) 2296 temp91 <- apply(temp93, 1:2,sum) # M x M 2297 temp91 <- rowSums(temp93, dims = 2) # M x M 2298 deta0.da[,,ss] <- -(temp90 + temp91) %*% sumWinv 2299 } 2300 2301 ans <- matrix(0,M,r) 2302 fred <- mux22(t(wz), z - eta, M = M, 2303 upper = FALSE, as.matrix = TRUE) 2304 fred.array <- array(t(fred %*% Aimat),c(r, 1, nn)) 2305 for (s in 1:r) { 2306 a1 <- colSums(fred %*% t(deta0.da[,, s])) 2307 a2 <- colSums(fred * etastar[, s]) 2308 temp92 <- array(detastar.da[, , s, ],c(M, r, nn)) 2309 temp93 <- mux7(temp92, fred.array) 2310 a3 <- rowSums(temp93, dims = 2) 2311 ans[,s] <- a1 + a2 + a3 2312 } 2313 2314 ans <- -2 * c(ans[cbindex, ]) 2315 2316 ans 2317} # rrr.deriv.gradient.fast 2318 2319 2320 2321 2322 2323 2324 2325 2326vellipse <- function(R, ratio = 1, orientation = 0, 2327 center = c(0, 0), N = 300) { 2328 if (length(center) != 2) 2329 stop("argument 'center' must be of length 2") 2330 theta <- 2*pi*(0:N)/N 2331 x1 <- R*cos(theta) 2332 y1 <- ratio*R*sin(theta) 2333 x <- center[1] + cos(orientation)*x1 - sin(orientation)*y1 2334 y <- center[2] + sin(orientation)*x1 + cos(orientation)*y1 2335 cbind(x, y) 2336} # vellipse 2337 2338 2339biplot.qrrvglm <- function(x, ...) { 2340 stop("biplot.qrrvglm has been replaced by the function lvplot.qrrvglm") 2341} 2342 2343 2344 2345 2346 2347 2348 lvplot.qrrvglm <- 2349 function(object, varI.latvar = FALSE, refResponse = NULL, 2350 add = FALSE, show.plot = TRUE, rug = TRUE, y = FALSE, 2351 type = c("fitted.values", "predictors"), 2352 xlab = paste("Latent Variable", 2353 if (Rank == 1) "" else " 1", sep = ""), 2354 ylab = if (Rank == 1) switch(type, predictors = "Predictors", 2355 fitted.values = "Fitted values") else "Latent Variable 2", 2356 pcex = par()$cex, pcol = par()$col, pch = par()$pch, 2357 llty = par()$lty, lcol = par()$col, llwd = par()$lwd, 2358 label.arg = FALSE, adj.arg = -0.1, 2359 ellipse = 0.95, Absolute = FALSE, 2360 elty = par()$lty, ecol = par()$col, elwd = par()$lwd, 2361 egrid = 200, 2362 chull.arg = FALSE, clty = 2, ccol = par()$col, clwd = par()$lwd, 2363 cpch = " ", 2364 C = FALSE, 2365 OriginC = c("origin", "mean"), 2366 Clty = par()$lty, Ccol = par()$col, Clwd = par()$lwd, 2367 Ccex = par()$cex, Cadj.arg = -0.1, stretchC = 1, 2368 sites = FALSE, spch = NULL, scol = par()$col, scex = par()$cex, 2369 sfont = par()$font, 2370 check.ok = TRUE, 2371 jitter.y = FALSE, 2372 ...) { 2373 if (mode(type) != "character" && mode(type) != "name") 2374 type <- as.character(substitute(type)) 2375 type <- match.arg(type, c("fitted.values", "predictors"))[1] 2376 2377 if (is.numeric(OriginC)) 2378 OriginC <- rep_len(OriginC, 2) else { 2379 if (mode(OriginC) != "character" && mode(OriginC) != "name") 2380 OriginC <- as.character(substitute(OriginC)) 2381 OriginC <- match.arg(OriginC, c("origin","mean"))[1] 2382 } 2383 2384 if (length(ellipse) > 1) 2385 stop("ellipse must be of length 1 or 0") 2386 if (is.logical(ellipse)) {ellipse <- if (ellipse) 0.95 else NULL} 2387 2388 Rank <- object@control$Rank 2389 if (Rank > 2) 2390 stop("can only handle rank 1 or 2 models") 2391 M <- object@misc$M 2392 NOS <- ncol(object@y) 2393 MSratio <- M / NOS # 1st value is g(mean) = quadratic form in latvar 2394 n <- object@misc$n 2395 colx2.index <- object@control$colx2.index 2396 cx1i <- object@control$colx1.index # May be NULL 2397 if (check.ok) 2398 if (!(length(cx1i) == 1 && names(cx1i) == "(Intercept)")) 2399 stop("latent variable plots allowable only for ", 2400 "noRRR = ~ 1 models") 2401 2402 Coef.list <- Coef(object, varI.latvar = varI.latvar, 2403 refResponse = refResponse) 2404 if ( C) Cmat <- Coef.list@C 2405 nustar <- Coef.list@latvar # n x Rank 2406 2407 if (!show.plot) return(nustar) 2408 2409 r.curves <- slot(object, type) # n times M (\boldeta or \boldmu) 2410 if (!add) { 2411 if (Rank == 1) { 2412 matplot(nustar, 2413 if ( y && type == "fitted.values") 2414 (if (jitter.y) jitter(object@y) else object@y) else r.curves, 2415 type = "n", xlab = xlab, ylab = ylab, ...) 2416 } else { # Rank == 2 2417 matplot(c(Coef.list@Optimum[1, ], nustar[, 1]), 2418 c(Coef.list@Optimum[2, ], nustar[, 2]), 2419 type = "n", xlab = xlab, ylab = ylab, ...) 2420 } 2421 } 2422 2423 2424 2425 2426 pch <- rep_len(pch, ncol(r.curves)) 2427 pcol <- rep_len(pcol, ncol(r.curves)) 2428 pcex <- rep_len(pcex, ncol(r.curves)) 2429 llty <- rep_len(llty, ncol(r.curves)) 2430 lcol <- rep_len(lcol, ncol(r.curves)) 2431 llwd <- rep_len(llwd, ncol(r.curves)) 2432 elty <- rep_len(elty, ncol(r.curves)) 2433 ecol <- rep_len(ecol, ncol(r.curves)) 2434 elwd <- rep_len(elwd, ncol(r.curves)) 2435 adj.arg <- rep_len(adj.arg, ncol(r.curves)) 2436 if ( C ) { 2437 Clwd <- rep_len(Clwd, nrow(Cmat)) 2438 Clty <- rep_len(Clty, nrow(Cmat)) 2439 Ccol <- rep_len(Ccol, nrow(Cmat)) 2440 Ccex <- rep_len(Ccex, nrow(Cmat)) 2441 Cadj.arg <- rep_len(Cadj.arg, nrow(Cmat)) 2442 } 2443 2444 if (Rank == 1) { 2445 for (i in 1:ncol(r.curves)) { 2446 xx <- nustar 2447 yy <- r.curves[,i] 2448 o <- sort.list(xx) 2449 xx <- xx[o] 2450 yy <- yy[o] 2451 lines(xx, yy, col = lcol[i], lwd = llwd[i], lty = llty[i]) 2452 if ( y && type == "fitted.values") { 2453 ypts <- if (jitter.y) jitter(object@y) else object@y 2454 if (NCOL(ypts) == ncol(r.curves)) 2455 points(xx, ypts[o, i], col = pcol[i], 2456 cex = pcex[i], pch = pch[i]) 2457 } 2458 } 2459 if (rug) 2460 rug(xx) 2461 } else { 2462 for (i in 1:ncol(r.curves)) 2463 points(Coef.list@Optimum[1, i], Coef.list@Optimum[2, i], 2464 col = pcol[i], cex = pcex[i], pch = pch[i]) 2465 if (label.arg) { 2466 for (i in 1:ncol(r.curves)) 2467 text(Coef.list@Optimum[1, i], Coef.list@Optimum[2, i], 2468 labels = (dimnames(Coef.list@Optimum)[[2]])[i], 2469 adj = adj.arg[i], col = pcol[i], cex = pcex[i]) 2470 } 2471 if (chull.arg) { 2472 hull <- chull(nustar[, 1], nustar[, 2]) 2473 hull <- c(hull, hull[1]) 2474 lines(nustar[hull, 1], nustar[hull, 2], type = "b", pch = cpch, 2475 lty = clty, col = ccol, lwd = clwd) 2476 } 2477 if (length(ellipse)) { 2478 ellipse.temp <- if (ellipse > 0) ellipse else 0.95 2479 if (ellipse < 0 && (!object@control$eq.tolerances || varI.latvar)) 2480 stop("an equal-tolerances assumption and 'varI.latvar = FALSE' ", 2481 "is needed for 'ellipse' < 0") 2482 if ( check.ok ) { 2483 colx1.index <- object@control$colx1.index 2484 if (!(length(colx1.index) == 1 && 2485 names(colx1.index) == "(Intercept)")) 2486 stop("can only plot ellipses for intercept models only") 2487 } 2488 for (i in 1:ncol(r.curves)) { 2489 cutpoint <- object@family@linkfun( if (Absolute) ellipse.temp 2490 else Coef.list@Maximum[i] * ellipse.temp, 2491 extra = object@extra) 2492 if (MSratio > 1) 2493 cutpoint <- cutpoint[1, 1] 2494 2495 cutpoint <- object@family@linkfun(Coef.list@Maximum[i], 2496 extra = object@extra) - cutpoint 2497 if (is.finite(cutpoint) && cutpoint > 0) { 2498 Mmat <- diag(rep_len(ifelse(object@control$Crow1positive, 2499 1, -1), 2500 Rank)) 2501 etoli <- 2502 eigen(t(Mmat) %*% Coef.list@Tolerance[,,i] %*% Mmat, symmetric = TRUE) 2503 A <- ifelse(etoli$val[1]>0, sqrt(2*cutpoint*etoli$val[1]), Inf) 2504 B <- ifelse(etoli$val[2]>0, sqrt(2*cutpoint*etoli$val[2]), Inf) 2505 if (ellipse < 0) 2506 A <- B <- -ellipse / 2 2507 2508 theta.angle <- asin(etoli$vector[2, 1]) * 2509 ifelse(object@control$Crow1positive[2], 1, -1) 2510 if (object@control$Crow1positive[1]) 2511 theta.angle <- pi - theta.angle 2512 if (all(is.finite(c(A,B)))) 2513 lines(vellipse(R = 2*A, ratio = B/A, 2514 orientation = theta.angle, 2515 center = Coef.list@Optimum[, i], 2516 N = egrid), 2517 lwd = elwd[i], col =ecol[i], lty = elty[i]) 2518 } 2519 } 2520 } 2521 2522 if ( C ) { 2523 if (is.character(OriginC) && OriginC == "mean") 2524 OriginC <- c(mean(nustar[, 1]), mean(nustar[, 2])) 2525 if (is.character(OriginC) && OriginC == "origin") 2526 OriginC <- c(0,0) 2527 for (i in 1:nrow(Cmat)) 2528 arrows(x0 = OriginC[1], y0 = OriginC[2], 2529 x1 = OriginC[1] + stretchC * Cmat[i, 1], 2530 y1 = OriginC[2] + stretchC * Cmat[i, 2], 2531 lty = Clty[i], col = Ccol[i], lwd = Clwd[i]) 2532 if (label.arg) { 2533 temp200 <- dimnames(Cmat)[[1]] 2534 for (i in 1:nrow(Cmat)) 2535 text(OriginC[1] + stretchC * Cmat[i, 1], 2536 OriginC[2] + stretchC * Cmat[i, 2], col = Ccol[i], 2537 labels = temp200[i], adj = Cadj.arg[i], 2538 cex = Ccex[i]) 2539 } 2540 } 2541 if (sites) { 2542 text(nustar[, 1], nustar[, 2], adj = 0.5, 2543 labels = if (is.null(spch)) dimnames(nustar)[[1]] else 2544 rep_len(spch, nrow(nustar)), col = scol, 2545 cex = scex, font = sfont) 2546 } 2547 } 2548 invisible(nustar) 2549} # lvplot.qrrvglm 2550 2551 2552 2553 2554 2555 2556lvplot.rrvglm <- 2557 function(object, 2558 A = TRUE, 2559 C = TRUE, 2560 scores = FALSE, show.plot = TRUE, 2561 groups = rep(1, n), 2562 gapC = sqrt(sum(par()$cxy^2)), scaleA = 1, 2563 xlab = "Latent Variable 1", 2564 ylab = "Latent Variable 2", 2565 Alabels= if (length(object@misc$predictors.names)) 2566 object@misc$predictors.names else param.names("LP", M), 2567 Aadj = par()$adj, 2568 Acex = par()$cex, 2569 Acol = par()$col, 2570 Apch = NULL, 2571 Clabels = rownames(Cmat), 2572 Cadj = par()$adj, 2573 Ccex = par()$cex, 2574 Ccol = par()$col, 2575 Clty = par()$lty, 2576 Clwd = par()$lwd, 2577 chull.arg = FALSE, 2578 ccex = par()$cex, 2579 ccol = par()$col, 2580 clty = par()$lty, 2581 clwd = par()$lwd, 2582 spch = NULL, 2583 scex = par()$cex, 2584 scol = par()$col, 2585 slabels = rownames(x2mat), ...) { 2586 2587 2588 if (object@control$Rank != 2 && show.plot) 2589 stop("can only handle rank-2 models") 2590 M <- object@misc$M 2591 n <- object@misc$n 2592 colx2.index <- object@control$colx2.index 2593 Coef.list <- Coef(object) 2594 Amat <- Coef.list@A 2595 Cmat <- Coef.list@C 2596 2597 Amat <- Amat * scaleA 2598 dimnames(Amat) <- list(object@misc$predictors.names, NULL) 2599 Cmat <- Cmat / scaleA 2600 2601 if (!length(object@x)) { 2602 object@x <- model.matrixvlm(object, type = "lm") 2603 } 2604 x2mat <- object@x[, colx2.index, drop = FALSE] 2605 nuhat <- x2mat %*% Cmat 2606 if (!show.plot) return(as.matrix(nuhat)) 2607 2608 index.nosz <- 1:M 2609 allmat <- rbind(if (A) Amat else NULL, 2610 if (C) Cmat else NULL, 2611 if (scores) nuhat else NULL) 2612 2613 plot(allmat[, 1], allmat[, 2], type = "n", 2614 xlab=xlab, ylab=ylab, ...) # xlim etc. supplied through ... 2615 2616 if (A) { 2617 Aadj <- rep_len(Aadj, length(index.nosz)) 2618 Acex <- rep_len(Acex, length(index.nosz)) 2619 Acol <- rep_len(Acol, length(index.nosz)) 2620 if (length(Alabels) != M) 2621 stop("'Alabels' must be of length ", M) 2622 if (length(Apch)) { 2623 Apch <- rep_len(Apch, length(index.nosz)) 2624 for (i in index.nosz) 2625 points(Amat[i, 1], 2626 Amat[i, 2], 2627 pch=Apch[i],cex = Acex[i],col=Acol[i]) 2628 } else { 2629 for (i in index.nosz) 2630 text(Amat[i, 1], Amat[i, 2], 2631 Alabels[i], cex = Acex[i], 2632 col =Acol[i], adj=Aadj[i]) 2633 } 2634 } 2635 2636 if (C) { 2637 p2 <- nrow(Cmat) 2638 gapC <- rep_len(gapC, p2) 2639 Cadj <- rep_len(Cadj, p2) 2640 Ccex <- rep_len(Ccex, p2) 2641 Ccol <- rep_len(Ccol, p2) 2642 Clwd <- rep_len(Clwd, p2) 2643 Clty <- rep_len(Clty, p2) 2644 if (length(Clabels) != p2) 2645 stop("'length(Clabels)' must be equal to ", p2) 2646 for (ii in 1:p2) { 2647 arrows(0, 0, Cmat[ii, 1], Cmat[ii, 2], 2648 lwd = Clwd[ii], lty = Clty[ii], col = Ccol[ii]) 2649 const <- 1 + gapC[ii] / sqrt(Cmat[ii, 1]^2 + Cmat[ii, 2]^2) 2650 text(const*Cmat[ii, 1], const*Cmat[ii, 2], 2651 Clabels[ii], cex = Ccex[ii], 2652 adj = Cadj[ii], col = Ccol[ii]) 2653 } 2654 } 2655 2656 if (scores) { 2657 ugrp <- unique(groups) 2658 nlev <- length(ugrp) # number of groups 2659 clty <- rep_len(clty, nlev) 2660 clwd <- rep_len(clwd, nlev) 2661 ccol <- rep_len(ccol, nlev) 2662 if (length(spch)) spch <- rep_len(spch, n) 2663 scol <- rep_len(scol, n) 2664 scex <- rep_len(scex, n) 2665 for (ii in ugrp) { 2666 gp <- groups == ii 2667 if (nlev > 1 && (length(unique(spch[gp])) != 1 || 2668 length(unique(scol[gp])) != 1 || 2669 length(unique(scex[gp])) != 1)) 2670 warning("spch/scol/scex is different for individuals ", 2671 "from the same group") 2672 2673 temp <- nuhat[gp,, drop = FALSE] 2674 if (length(spch)) { 2675 points(temp[, 1], temp[, 2], cex = scex[gp], pch = spch[gp], 2676 col = scol[gp]) 2677 } else { 2678 text(temp[, 1], temp[, 2], label = slabels, cex = scex[gp], 2679 col = scol[gp]) 2680 } 2681 if (chull.arg) { 2682 hull <- chull(temp[, 1], temp[, 2]) 2683 hull <- c(hull, hull[1]) 2684 lines(temp[hull, 1], temp[hull, 2], 2685 type = "b", lty = clty[ii], 2686 col = ccol[ii], lwd = clwd[ii], pch = " ") 2687 } 2688 } 2689 } 2690 2691 invisible(nuhat) 2692} # lvplot.rrvglm 2693 2694 2695 2696 2697 2698 2699 Coef.rrvglm <- function(object, ...) { 2700 M <- object@misc$M 2701 n <- object@misc$n 2702 colx1.index <- object@control$colx1.index 2703 colx2.index <- object@control$colx2.index 2704 p1 <- length(colx1.index) # May be 0 2705 Amat <- object@constraints[[colx2.index[1]]] 2706 2707 B1mat <- if (p1) 2708 coefvlm(object, matrix.out = TRUE)[colx1.index,, drop = FALSE] else 2709 NULL 2710 2711 2712 C.try <- coefvlm(object, 2713 matrix.out = TRUE)[colx2.index, , drop = FALSE] 2714 2715 2716 Cmat <- C.try %*% Amat %*% solve(t(Amat) %*% Amat) 2717 2718 2719 Rank <- object@control$Rank 2720 latvar.names <- param.names("latvar", Rank, skip1 = TRUE) 2721 dimnames(Amat) <- list(object@misc$predictors.names, latvar.names) 2722 dimnames(Cmat) <- list(dimnames(Cmat)[[1]], latvar.names) 2723 2724 ans <- new(Class = "Coef.rrvglm", 2725 A = Amat, 2726 C = Cmat, 2727 Rank = Rank, 2728 colx2.index = colx2.index) 2729 2730 if (!is.null(colx1.index)) { 2731 ans@colx1.index <- colx1.index 2732 ans@B1 <- B1mat 2733 } 2734 2735 if (object@control$Corner) 2736 ans@Atilde <- Amat[-c(object@control$Index.corner, 2737 object@control$str0),, drop = FALSE] 2738 ans 2739} # Coef.rrvglm 2740 2741 2742 2743 2744setMethod("Coef", "rrvglm", 2745 function(object, ...) Coef.rrvglm(object, ...)) 2746 2747 2748 2749 2750 2751show.Coef.rrvglm <- function(x, ...) { 2752 2753 object <- x 2754 2755 cat("A matrix:\n") 2756 print(object@A, ...) 2757 2758 cat("\nC matrix:\n") 2759 print(object@C, ...) 2760 2761 p1 <- length(object@colx1.index) 2762 if (p1) { 2763 cat("\nB1 matrix:\n") 2764 print(object@B1, ...) 2765 } 2766 2767 invisible(object) 2768} # show.Coef.rrvglm 2769 2770 2771 if (!isGeneric("biplot")) 2772 setGeneric("biplot", function(x, ...) standardGeneric("biplot")) 2773 2774 2775setMethod("Coef", "qrrvglm", function(object, ...) 2776 Coef.qrrvglm(object, ...)) 2777 2778 2779 2780setMethod("biplot", "qrrvglm", 2781 function(x, ...) { 2782 biplot.qrrvglm(x, ...)}) 2783 2784setMethod("lvplot", "qrrvglm", 2785 function(object, ...) { 2786 invisible(lvplot.qrrvglm(object, ...))}) 2787 2788setMethod("lvplot", "rrvglm", 2789 function(object, ...) { 2790 invisible(lvplot.rrvglm(object, ...))}) 2791 2792 2793biplot.rrvglm <- function(x, ...) 2794 lvplot(object = x, ...) 2795 2796setMethod("biplot", "rrvglm", function(x, ...) 2797 invisible(biplot.rrvglm(x, ...))) 2798 2799 2800 2801 2802summary.qrrvglm <- 2803 function(object, 2804 varI.latvar = FALSE, refResponse = NULL, ...) { 2805 answer <- object 2806 answer@post$Coef <- Coef(object, 2807 varI.latvar = varI.latvar, 2808 refResponse = refResponse, 2809 ...) # Store it here; non-elegant 2810 2811 if (length((answer@post$Coef)@dispersion) && 2812 length(object@misc$estimated.dispersion) && 2813 object@misc$estimated.dispersion) 2814 answer@dispersion <- 2815 answer@misc$dispersion <- (answer@post$Coef)@dispersion 2816 2817 as(answer, "summary.qrrvglm") 2818} # summary.qrrvglm 2819 2820 2821 2822show.summary.qrrvglm <- function(x, ...) { 2823 2824 2825 2826 cat("\nCall:\n") 2827 dput(x@call) 2828 2829 print(x@post$Coef, ...) # non-elegant programming 2830 2831 if (length(x@dispersion) > 1) { 2832 cat("\nDispersion parameters:\n") 2833 if (length(x@misc$ynames)) { 2834 names(x@dispersion) <- x@misc$ynames 2835 print(x@dispersion, ...) 2836 } else { 2837 cat(x@dispersion, fill = TRUE) 2838 } 2839 cat("\n") 2840 } else if (length(x@dispersion) == 1) { 2841 cat("\nDispersion parameter: ", x@dispersion, "\n") 2842 } 2843 2844} # show.summary.qrrvglm 2845 2846 2847 setClass("summary.qrrvglm", contains = "qrrvglm") 2848 2849 2850 2851 2852 2853 2854 2855setMethod("summary", "qrrvglm", 2856 function(object, ...) 2857 summary.qrrvglm(object, ...)) 2858 2859 2860 2861 2862 2863 setMethod("show", "summary.qrrvglm", 2864 function(object) 2865 show.summary.qrrvglm(object)) 2866 2867 2868 2869 2870 2871setMethod("show", "Coef.rrvglm", function(object) 2872 show.Coef.rrvglm(object)) 2873 2874 2875 2876 2877 2878 grc <- function(y, Rank = 1, Index.corner = 2:(1+Rank), 2879 str0 = 1, 2880 summary.arg = FALSE, h.step = 0.0001, ...) { 2881 2882 2883 2884 myrrcontrol <- rrvglm.control(Rank = Rank, 2885 Index.corner = Index.corner, 2886 str0 = str0, ...) 2887 object.save <- y 2888 if (is(y, "rrvglm")) { 2889 y <- object.save@y 2890 } else { 2891 y <- as.matrix(y) 2892 y <- as(y, "matrix") 2893 } 2894 if (length(dim(y)) != 2 || nrow(y) < 3 || ncol(y) < 3) 2895 stop("y must be a matrix with >= 3 rows & columns, ", 2896 "or a rrvglm() object") 2897 2898 ei <- function(i, n) diag(n)[, i, drop = FALSE] 2899 .grc.df <- data.frame(Row.2 = I.col(2, nrow(y))) 2900 2901 yn1 <- if (length(dimnames(y)[[1]])) dimnames(y)[[1]] else 2902 param.names("X2.", nrow(y)) 2903 warn.save <- options()$warn 2904 options(warn = -3) # Suppress the warnings (hopefully, temporarily) 2905 if (any(!is.na(as.numeric(substring(yn1, 1, 1))))) 2906 yn1 <- param.names("X2.", nrow(y)) 2907 options(warn = warn.save) 2908 2909 Row. <- factor(1:nrow(y)) 2910 modmat.row <- model.matrix( ~ Row.) 2911 Col. <- factor(1:ncol(y)) 2912 modmat.col <- model.matrix( ~ Col.) 2913 2914 cms <- list("(Intercept)" = matrix(1, ncol(y), 1)) 2915 for (ii in 2:nrow(y)) { 2916 cms[[paste("Row.", ii, sep = "")]] <- matrix(1, ncol(y), 1) 2917 .grc.df[[paste("Row.", ii, sep = "")]] <- modmat.row[,ii] 2918 } 2919 for (ii in 2:ncol(y)) { 2920 cms[[paste("Col.", ii, sep = "")]] <- 2921 modmat.col[,ii, drop = FALSE] 2922 .grc.df[[paste("Col.", ii, sep = "")]] <- rep_len(1, nrow(y)) 2923 } 2924 for (ii in 2:nrow(y)) { 2925 cms[[yn1[ii]]] <- diag(ncol(y)) 2926 .grc.df[[yn1[ii]]] <- I.col(ii, nrow(y)) 2927 } 2928 2929 dimnames(.grc.df) <- list(if (length(dimnames(y)[[1]])) 2930 dimnames(y)[[1]] else 2931 as.character(1:nrow(y)), 2932 dimnames(.grc.df)[[2]]) 2933 2934 str1 <- "~ Row.2" 2935 if (nrow(y) > 2) 2936 for (ii in 3:nrow(y)) 2937 str1 <- paste(str1, paste("Row.", ii, sep = ""), sep = " + ") 2938 for (ii in 2:ncol(y)) 2939 str1 <- paste(str1, paste("Col.", ii, sep = ""), sep = " + ") 2940 str2 <- paste("y ", str1) 2941 for (ii in 2:nrow(y)) 2942 str2 <- paste(str2, yn1[ii], sep = " + ") 2943 myrrcontrol$noRRR <- as.formula(str1) # Overwrite this 2944 2945 assign(".grc.df", .grc.df, envir = VGAMenv) 2946 2947 warn.save <- options()$warn 2948 options(warn = -3) # Suppress the warnings (hopefully, temporarily) 2949 answer <- if (is(object.save, "rrvglm")) object.save else 2950 rrvglm(as.formula(str2), family = poissonff, 2951 constraints = cms, control = myrrcontrol, 2952 data = .grc.df) 2953 options(warn = warn.save) 2954 2955 if (summary.arg) { 2956 answer <- as(answer, "rrvglm") 2957 2958 answer <- summary.rrvglm(answer, h.step = h.step) 2959 } else { 2960 answer <- as(answer, "grc") 2961 } 2962 2963 if (exists(".grc.df", envir = VGAMenv)) 2964 rm(".grc.df", envir = VGAMenv) 2965 2966 answer 2967} # grc 2968 2969 2970summary.grc <- function(object, ...) { 2971 grc(object, summary.arg= TRUE, ...) 2972} 2973 2974 2975 2976 2977 2978 2979trplot.qrrvglm <- 2980 function(object, 2981 which.species = NULL, 2982 add = FALSE, show.plot = TRUE, 2983 label.sites = FALSE, 2984 sitenames = rownames(object@y), 2985 axes.equal = TRUE, 2986 cex = par()$cex, 2987 col = 1:(nos*(nos-1)/2), 2988 log = "", 2989 lty = rep_len(par()$lty, nos*(nos-1)/2), 2990 lwd = rep_len(par()$lwd, nos*(nos-1)/2), 2991 tcol = rep_len(par()$col, nos*(nos-1)/2), 2992 xlab = NULL, ylab = NULL, 2993 main = "", # "Trajectory plot", 2994 type = "b", 2995 check.ok = TRUE, ...) { 2996 coef.obj <- Coef(object) # use defaults for those two arguments 2997 if (coef.obj@Rank != 1) 2998 stop("object must be a rank-1 model") 2999 fv <- fitted(object) 3000 modelno <- object@control$modelno # 1, 2, 3, or 0 3001 NOS <- ncol(fv) # Number of species 3002 M <- object@misc$M # 3003 nn <- nrow(fv) # Number of sites 3004 if (length(sitenames)) 3005 sitenames <- rep_len(sitenames, nn) 3006 sppNames <- dimnames(object@y)[[2]] 3007 if (!length(which.species)) { 3008 which.species <- sppNames[1:NOS] 3009 which.species.numer <- 1:NOS 3010 } else 3011 if (is.numeric(which.species)) { 3012 which.species.numer <- which.species 3013 which.species <- sppNames[which.species.numer] # Convert to char 3014 } else { 3015 which.species.numer <- match(which.species, sppNames) 3016 } 3017 nos <- length(which.species) # nos = number of species to be plotted 3018 3019 if (length(which.species.numer) <= 1) 3020 stop("must have at least 2 species to be plotted") 3021 cx1i <- object@control$colx1.index 3022 if (check.ok) 3023 if (!(length(cx1i) == 1 && names(cx1i) == "(Intercept)")) 3024 stop("trajectory plots allowable only for noRRR = ~ 1 models") 3025 3026 first.spp <- iam(1, 1,M = M,both = TRUE,diag = FALSE)$row.index 3027 second.spp <- iam(1, 1,M = M,both = TRUE,diag = FALSE)$col.index 3028 myxlab <- if (length(which.species.numer) == 2) { 3029 paste("Fitted value for", 3030 if (is.character(which.species.numer)) 3031 which.species.numer[1] else 3032 sppNames[which.species.numer[1]]) 3033 } else "Fitted value for 'first' species" 3034 myxlab <- if (length(xlab)) xlab else myxlab 3035 myylab <- if (length(which.species.numer) == 2) { 3036 paste("Fitted value for", 3037 if (is.character(which.species.numer)) 3038 which.species.numer[2] else 3039 sppNames[which.species.numer[2]]) 3040 } else "Fitted value for 'second' species" 3041 myylab <- if (length(ylab)) ylab else myylab 3042 if (!add) { 3043 xxx <- if (axes.equal) fv[,which.species.numer] else 3044 fv[,which.species.numer[first.spp]] 3045 yyy <- if (axes.equal) fv[,which.species.numer] else 3046 fv[,which.species.numer[second.spp]] 3047 matplot(xxx, yyy, type = "n", log = log, xlab = myxlab, 3048 ylab = myylab, main = main, ...) 3049 } 3050 3051 lwd <- rep_len(lwd, nos*(nos-1)/2) 3052 col <- rep_len(col, nos*(nos-1)/2) 3053 lty <- rep_len(lty, nos*(nos-1)/2) 3054 tcol <- rep_len(tcol, nos*(nos-1)/2) 3055 3056 oo <- order(coef.obj@latvar) # Sort by the latent variable 3057 ii <- 0 3058 col <- rep_len(col, nos*(nos-1)/2) 3059 species.names <- NULL 3060 if (show.plot) 3061 for (i1 in seq(which.species.numer)) { 3062 for (i2 in seq(which.species.numer)) 3063 if (i1 < i2) { 3064 ii <- ii + 1 3065 species.names <- rbind(species.names, 3066 cbind(sppNames[i1], sppNames[i2])) 3067 matplot(fv[oo, which.species.numer[i1]], 3068 fv[oo, which.species.numer[i2]], 3069 type = type, add = TRUE, 3070 lty = lty[ii], lwd = lwd[ii], col = col[ii], 3071 pch = if (label.sites) " " else "*" ) 3072 if (label.sites && length(sitenames)) 3073 text(fv[oo, which.species.numer[i1]], 3074 fv[oo, which.species.numer[i2]], 3075 labels = sitenames[oo], cex = cex, col = tcol[ii]) 3076 } 3077 } 3078 invisible(list(species.names = species.names, 3079 sitenames = sitenames[oo])) 3080} # trplot.qrrvglm 3081 3082 3083 3084 if (!isGeneric("trplot")) 3085 setGeneric("trplot", 3086 function(object, ...) standardGeneric("trplot")) 3087setMethod("trplot", "qrrvglm", 3088 function(object, ...) trplot.qrrvglm(object, ...)) 3089 3090setMethod("trplot", "rrvgam", 3091 function(object, ...) trplot.qrrvglm(object, ...)) 3092 3093 3094 3095 3096 3097vcovrrvglm <- function(object, ...) { 3098 summary.rrvglm(object, ...)@cov.unscaled 3099} # vcovrrvglm 3100 3101 3102 3103vcovqrrvglm <- 3104 function(object, 3105 ...) { 3106 3107 3108 3109 object@control$trace <- FALSE # Suppress this for vglm.fit(). 3110 fit1 <- fnumat2R(object, refit.model = TRUE) 3111 RR <- fit1$R 3112 covun <- chol2inv(RR) 3113 dimnames(covun) <- list(names(coef(fit1)), names(coef(fit1))) 3114 return(covun) 3115 3116 3117 3118 stop("this function is not yet completed") 3119 3120 3121 3122I.tolerances = object@control$eq.tolerances 3123 3124 3125 3126 3127 3128 if (mode(MaxScale) != "character" && mode(MaxScale) != "name") 3129 MaxScale <- as.character(substitute(MaxScale)) 3130 MaxScale <- match.arg(MaxScale, c("predictors", "response"))[1] 3131 if (MaxScale != "predictors") 3132 stop("can currently only handle MaxScale='predictors'") 3133 3134 sobj <- summary(object) 3135 cobj <- Coef(object, I.tolerances = I.tolerances, ...) 3136 M <- nrow(cobj@A) 3137 dispersion <- rep_len(dispersion, M) 3138 if (cobj@Rank != 1) 3139 stop("object must be a rank 1 model") 3140 3141 dvecMax <- cbind(1, -0.5 * cobj@A / c(cobj@D), 3142 (cobj@A / c(2*cobj@D))^2) 3143 dvecTol <- cbind(0, 0, 1 / c(-2 * cobj@D)^1.5) 3144 dvecOpt <- cbind(0, -0.5 / c(cobj@D), 0.5 * cobj@A / c(cobj@D^2)) 3145 3146 if ((length(object@control$colx1.index) != 1) || 3147 (names(object@control$colx1.index) != "(Intercept)")) 3148 stop("Can only handle noRRR=~1 models") 3149 3150 okvals <- c(3*M, 2*M+1) 3151 if (all(length(coef(object)) != okvals)) 3152 stop("Can only handle intercepts-only model with ", 3153 "eq.tolerances = FALSE") 3154 3155 answer <- NULL 3156 Cov.unscaled <- array(NA_real_, c(3, 3, M), dimnames = list( 3157 c("(Intercept)", "latvar", "latvar^2"), 3158 c("(Intercept)", "latvar", "latvar^2"), dimnames(cobj@D)[[3]])) 3159 for (spp in 1:M) { 3160 index <- c(M + ifelse(object@control$eq.tolerances, 1, M) + spp, 3161 spp, 3162 M + ifelse(object@control$eq.tolerances, 1, spp)) 3163 vcov <- Cov.unscaled[,,spp] <- 3164 sobj@cov.unscaled[index, index] # Order is A, D, B1 3165 se2Max <- dvecMax[spp,, drop = FALSE] %*% vcov %*% cbind(dvecMax[spp,]) 3166 se2Tol <- dvecTol[spp,, drop = FALSE] %*% vcov %*% cbind(dvecTol[spp,]) 3167 se2Opt <- dvecOpt[spp,, drop = FALSE] %*% vcov %*% cbind(dvecOpt[spp,]) 3168 answer <- rbind(answer, dispersion[spp]^0.5 * 3169 c(se2Opt = se2Opt, se2Tol = se2Tol, se2Max = se2Max)) 3170 } 3171 3172 link.function <- if (MaxScale == "predictors") 3173 remove.arg(object@misc$predictors.names[1]) else "" 3174 dimnames(answer) <- list(dimnames(cobj@D)[[3]], 3175 c("Optimum", "Tolerance", 3176 if (nchar(link.function)) 3177 paste(link.function, "(Maximum)", sep = "") else 3178 "Maximum")) 3179 NAthere <- is.na(answer %*% rep_len(1, 3)) 3180 answer[NAthere,] <- NA # NA in tolerance means NA everywhere else 3181 new(Class = "vcov.qrrvglm", 3182 Cov.unscaled = Cov.unscaled, 3183 dispersion = dispersion, 3184 se = sqrt(answer)) 3185} # vcovqrrvglm 3186 3187 3188 3189setMethod("vcov", "rrvglm", function(object, ...) 3190 vcovrrvglm(object, ...)) 3191 3192 3193setMethod("vcov", "qrrvglm", function(object, ...) 3194 vcovqrrvglm(object, ...)) 3195 3196 3197setClass(Class = "vcov.qrrvglm", representation( 3198 Cov.unscaled = "array", # permuted cov.unscaled 3199 dispersion = "numeric", 3200 se = "matrix")) 3201 3202 3203 3204 3205model.matrixqrrvglm <- 3206 function(object, 3207 type = c("latvar", "lm", "vlm"), ...) { 3208 3209 3210 if (mode(type) != "character" && mode(type) != "name") 3211 type <- as.character(substitute(type)) 3212 type <- match.arg(type, c("latvar", "lm", "vlm"))[1] 3213 3214 switch(type, 3215 latvar = Coef(object, ...)@latvar, 3216 lm = object@x, # 20180516 was "vlm" 3217 vlm = fnumat2R(object) # 20180516 change 3218 ) 3219} # model.matrixqrrvglm 3220 3221 3222setMethod("model.matrix", "qrrvglm", function(object, ...) 3223 model.matrixqrrvglm(object, ...)) 3224 3225 3226 3227 3228 3229 3230 3231 3232perspqrrvglm <- 3233 function(x, varI.latvar = FALSE, refResponse = NULL, 3234 show.plot = TRUE, 3235 xlim = NULL, ylim = NULL, 3236 zlim = NULL, # zlim ignored if Rank == 1 3237 gridlength = if (Rank == 1) 301 else c(51, 51), 3238 which.species = NULL, 3239 xlab = if (Rank == 1) 3240 "Latent Variable" else "Latent Variable 1", 3241 ylab = if (Rank == 1) 3242 "Expected Value" else "Latent Variable 2", 3243 zlab = "Expected value", 3244 labelSpecies = FALSE, # For Rank == 1 only 3245 stretch = 1.05, # quick and dirty, Rank == 1 only 3246 main = "", 3247 ticktype = "detailed", 3248 col = if (Rank == 1) par()$col else "white", 3249 llty = par()$lty, llwd = par()$lwd, 3250 add1 = FALSE, 3251 ...) { 3252 oylim <- ylim 3253 object <- x # Do not like x as the primary argument 3254 coef.obj <- Coef(object, varI.latvar = varI.latvar, 3255 refResponse = refResponse) 3256 if ((Rank <- coef.obj@Rank) > 2) 3257 stop("object must be a rank-1 or rank-2 model") 3258 fv <- fitted(object) 3259 NOS <- ncol(fv) # Number of species 3260 M <- object@misc$M 3261 3262 xlim <- rep_len(if (length(xlim)) xlim else range(coef.obj@latvar[, 1]), 3263 2) 3264 if (!length(oylim)) { 3265 ylim <- if (Rank == 1) c(0, max(fv) * stretch) else 3266 rep_len(range(coef.obj@latvar[, 2]), 2) 3267 } 3268 gridlength <- rep_len(gridlength, Rank) 3269 latvar1 <- seq(xlim[1], xlim[2], length = gridlength[1]) 3270 if (Rank == 1) { 3271 m <- cbind(latvar1) 3272 } else { 3273 latvar2 <- seq(ylim[1], ylim[2], length = gridlength[2]) 3274 m <- expand.grid(latvar1,latvar2) 3275 } 3276 3277 if (dim(coef.obj@B1)[1] != 1 || 3278 dimnames(coef.obj@B1)[[1]] != "(Intercept)") 3279 stop("noRRR = ~ 1 is needed") 3280 LP <- coef.obj@A %*% t(cbind(m)) # M by n 3281 LP <- LP + c(coef.obj@B1) # Assumes \bix_1 = 1 (intercept only) 3282 3283 mm <- as.matrix(m) 3284 N <- ncol(LP) 3285 for (jay in 1:M) { 3286 for (ii in 1:N) { 3287 LP[jay, ii] <- LP[jay, ii] + 3288 mm[ii, , drop = FALSE] %*% 3289 coef.obj@D[,,jay] %*% 3290 t(mm[ii, , drop = FALSE]) 3291 } 3292 } 3293 LP <- t(LP) # n by M 3294 3295 3296 fitvals <- object@family@linkinv(LP, extra = object@extra) # n by NOS 3297 dimnames(fitvals) <- list(NULL, dimnames(fv)[[2]]) 3298 sppNames <- dimnames(object@y)[[2]] 3299 if (!length(which.species)) { 3300 which.species <- sppNames[1:NOS] 3301 which.species.numer <- 1:NOS 3302 } else 3303 if (is.numeric(which.species)) { 3304 which.species.numer <- which.species 3305 which.species <- sppNames[which.species.numer] # Convert to char 3306 } else { 3307 which.species.numer <- match(which.species, sppNames) 3308 } 3309 if (Rank == 1) { 3310 if (show.plot) { 3311 if (!length(oylim)) 3312 ylim <- c(0, max(fitvals[,which.species.numer]) * 3313 stretch) # A revision 3314 col <- rep_len(col, length(which.species.numer)) 3315 llty <- rep_len(llty, length(which.species.numer)) 3316 llwd <- rep_len(llwd, length(which.species.numer)) 3317 if (!add1) 3318 matplot(latvar1, fitvals, xlab = xlab, ylab = ylab, 3319 type = "n", 3320 main = main, xlim = xlim, ylim = ylim, ...) 3321 for (jloc in seq_along(which.species.numer)) { 3322 ptr2 <- which.species.numer[jloc] # points to species column 3323 lines(latvar1, fitvals[, ptr2], 3324 col = col[jloc], 3325 lty = llty[jloc], 3326 lwd = llwd[jloc], ...) 3327 if (labelSpecies) { 3328 ptr1 <- (1:nrow(fitvals))[max(fitvals[, ptr2]) == 3329 fitvals[, ptr2]] 3330 ptr1 <- ptr1[1] 3331 text(latvar1[ptr1], 3332 fitvals[ptr1, ptr2] + (stretch-1) * diff(range(ylim)), 3333 label = sppNames[jloc], col = col[jloc], ...) 3334 } 3335 } 3336 } 3337 } else { 3338 max.fitted <- matrix(fitvals[, which.species[1]], 3339 length(latvar1), length(latvar2)) 3340 if (length(which.species) > 1) 3341 for (jlocal in which.species[-1]) { 3342 max.fitted <- pmax(max.fitted, 3343 matrix(fitvals[, jlocal], 3344 length(latvar1), length(latvar2))) 3345 } 3346 if (!length(zlim)) 3347 zlim <- range(max.fitted, na.rm = TRUE) 3348 3349 3350 perspdefault <- getS3method("persp", "default") 3351 if (show.plot) 3352 perspdefault(latvar1, latvar2, max.fitted, 3353 zlim = zlim, 3354 xlab = xlab, ylab = ylab, zlab = zlab, 3355 ticktype = ticktype, col = col, main = main, ...) 3356 } 3357 3358 invisible(list(fitted = fitvals, 3359 latvar1grid = latvar1, 3360 latvar2grid = if (Rank == 2) latvar2 else NULL, 3361 max.fitted = if (Rank == 2) max.fitted else NULL)) 3362} # perspqrrvglm 3363 3364 3365 3366 if (!isGeneric("persp")) 3367 setGeneric("persp", function(x, ...) standardGeneric("persp"), 3368 package = "VGAM") 3369 3370setMethod("persp", "qrrvglm", 3371 function(x, ...) perspqrrvglm(x = x, ...)) 3372 3373 3374 3375 3376 3377 3378 3379 3380 3381Rank.rrvglm <- function(object, ...) { 3382 object@control$Rank 3383} 3384 3385 3386Rank.qrrvglm <- function(object, ...) { 3387 object@control$Rank 3388} 3389 3390 3391Rank.rrvgam <- function(object, ...) { 3392 object@control$Rank 3393} 3394 3395 3396 3397 3398concoef.qrrvglm <- function(object, varI.latvar = FALSE, 3399 refResponse = NULL, ...) { 3400 Coef(object, varI.latvar = varI.latvar, 3401 refResponse = refResponse, ...)@C 3402} 3403 3404 3405concoef.Coef.qrrvglm <- function(object, ...) { 3406 if (length(list(...))) 3407 warning("Too late! Ignoring the extra arguments") 3408 object@C 3409} 3410 3411 3412latvar.rrvglm <- function(object, ...) { 3413 ans <- lvplot(object, show.plot = FALSE) 3414 if (ncol(ans) == 1) 3415 dimnames(ans) <- list(dimnames(ans)[[1]], "lv") 3416 ans 3417} 3418 3419 3420 3421latvar.qrrvglm <- function(object, 3422 varI.latvar = FALSE, 3423 refResponse = NULL, ...) { 3424 Coef(object, 3425 varI.latvar = varI.latvar, 3426 refResponse = refResponse, ...)@latvar 3427} 3428 3429 3430latvar.Coef.qrrvglm <- function(object, ...) { 3431 if (length(list(...))) 3432 warning("Too late! Ignoring the extra arguments") 3433 object@latvar 3434} 3435 3436 3437Max.qrrvglm <- 3438 function(object, varI.latvar = FALSE, 3439 refResponse = NULL, ...) { 3440 Coef(object, varI.latvar = varI.latvar, 3441 refResponse = refResponse, ...)@Maximum 3442} 3443 3444 3445Max.Coef.qrrvglm <- function(object, ...) { 3446 if (length(list(...))) 3447 warning("Too late! Ignoring the extra arguments") 3448 if (any(slotNames(object) == "Maximum")) 3449 object@Maximum else 3450 Max(object, ...) 3451} 3452 3453 3454Opt.qrrvglm <- 3455 function(object, varI.latvar = FALSE, refResponse = NULL, ...) { 3456 Coef(object, varI.latvar = varI.latvar, 3457 refResponse = refResponse, ...)@Optimum 3458} 3459 3460 3461Opt.Coef.qrrvglm <- function(object, ...) { 3462 if (length(list(...))) 3463 warning("Too late! Ignoring the extra arguments") 3464 object@Optimum 3465} 3466 3467 3468Tol.qrrvglm <- 3469 function(object, varI.latvar = FALSE, refResponse = NULL, ...) { 3470 Coef(object, varI.latvar = varI.latvar, 3471 refResponse = refResponse, ...)@Tolerance 3472} 3473 3474 3475Tol.Coef.qrrvglm <- function(object, ...) { 3476 if (length(list(...))) 3477 warning("Too late! Ignoring the extra arguments") 3478 if (any(slotNames(object) == "Tolerance")) 3479 object@Tolerance else Tol(object, ...) 3480} 3481 3482 3483 3484 if (FALSE) { 3485 if (!isGeneric("ccoef")) 3486 setGeneric("ccoef", function(object, ...) { 3487 .Deprecated("concoef") 3488 3489 standardGeneric("ccoef") 3490 }) 3491 3492setMethod("ccoef", "rrvglm", 3493 function(object, ...) concoef.qrrvglm(object, ...)) 3494setMethod("ccoef", "qrrvglm", 3495 function(object, ...) concoef.qrrvglm(object, ...)) 3496setMethod("ccoef", "Coef.rrvglm", 3497 function(object, ...) concoef.Coef.qrrvglm(object, ...)) 3498setMethod("ccoef", "Coef.qrrvglm", 3499 function(object, ...) concoef.Coef.qrrvglm(object, ...)) 3500} 3501 3502 3503 if (!isGeneric("concoef")) 3504 setGeneric("concoef", function(object, ...) 3505 standardGeneric("concoef")) 3506 3507setMethod("concoef", "rrvglm", 3508 function(object, ...) concoef.qrrvglm(object, ...)) 3509setMethod("concoef", "qrrvglm", 3510 function(object, ...) concoef.qrrvglm(object, ...)) 3511setMethod("concoef", "Coef.rrvglm", 3512 function(object, ...) concoef.Coef.qrrvglm(object, ...)) 3513setMethod("concoef", "Coef.qrrvglm", 3514 function(object, ...) concoef.Coef.qrrvglm(object, ...)) 3515 3516 3517 3518 3519 3520 3521 3522 3523 3524setMethod("coef", "qrrvglm", 3525 function(object, ...) Coef.qrrvglm(object, ...)) 3526setMethod("coefficients", "qrrvglm", 3527 function(object, ...) Coef.qrrvglm(object, ...)) 3528 3529 3530 3531 if (!isGeneric("lv")) 3532 setGeneric("lv", 3533 function(object, ...) { 3534 .Deprecated("latvar") 3535 3536 standardGeneric("lv") 3537 }) 3538 3539 3540setMethod("lv", "rrvglm", 3541 function(object, ...) { 3542 3543 latvar.rrvglm(object, ...) 3544 }) 3545setMethod("lv", "qrrvglm", 3546 function(object, ...) { 3547 3548 latvar.qrrvglm(object, ...) 3549 }) 3550setMethod("lv", "Coef.rrvglm", 3551 function(object, ...) { 3552 3553 latvar.Coef.qrrvglm(object, ...) 3554 }) 3555setMethod("lv", "Coef.qrrvglm", 3556 function(object, ...) { 3557 3558 latvar.Coef.qrrvglm(object, ...) 3559 }) 3560 3561 3562 if (!isGeneric("latvar")) 3563 setGeneric("latvar", 3564 function(object, ...) standardGeneric("latvar")) 3565setMethod("latvar", "rrvglm", 3566 function(object, ...) latvar.rrvglm(object, ...)) 3567setMethod("latvar", "qrrvglm", 3568 function(object, ...) latvar.qrrvglm(object, ...)) 3569setMethod("latvar", "Coef.rrvglm", 3570 function(object, ...) latvar.Coef.qrrvglm(object, ...)) 3571setMethod("latvar", "Coef.qrrvglm", 3572 function(object, ...) latvar.Coef.qrrvglm(object, ...)) 3573 3574 3575 if (!isGeneric("Max")) 3576 setGeneric("Max", 3577 function(object, ...) standardGeneric("Max")) 3578setMethod("Max", "qrrvglm", 3579 function(object, ...) Max.qrrvglm(object, ...)) 3580setMethod("Max", "Coef.qrrvglm", 3581 function(object, ...) Max.Coef.qrrvglm(object, ...)) 3582 3583 3584 3585setMethod("Max", "rrvgam", 3586 function(object, ...) Coef(object, ...)@Maximum) 3587 3588 3589 3590 3591 if (!isGeneric("Opt")) 3592 setGeneric("Opt", 3593 function(object, ...) standardGeneric("Opt")) 3594setMethod("Opt", "qrrvglm", 3595 function(object, ...) Opt.qrrvglm(object, ...)) 3596setMethod("Opt", "Coef.qrrvglm", 3597 function(object, ...) Opt.Coef.qrrvglm(object, ...)) 3598 3599 3600setMethod("Opt", "rrvgam", 3601 function(object, ...) Coef(object, ...)@Optimum) 3602 3603 3604 3605 3606 if (!isGeneric("Tol")) 3607 setGeneric("Tol", 3608 function(object, ...) standardGeneric("Tol")) 3609setMethod("Tol", "qrrvglm", 3610 function(object, ...) Tol.qrrvglm(object, ...)) 3611setMethod("Tol", "Coef.qrrvglm", 3612 function(object, ...) Tol.Coef.qrrvglm(object, ...)) 3613 3614 3615 3616cgo <- function(...) { 3617 stop("The function 'cgo' has been renamed 'cqo'. ", 3618 "Ouch! Sorry!") 3619} 3620 3621clo <- function(...) { 3622 stop("Constrained linear ordination is fitted with ", 3623 "the function 'rrvglm'") 3624} 3625 3626 3627 3628 3629is.bell.vlm <- 3630is.bell.rrvglm <- function(object, ...) { 3631 M <- object@misc$M 3632 ynames <- object@misc$ynames 3633 ans <- rep_len(FALSE, M) 3634 if (length(ynames)) names(ans) <- ynames 3635 ans 3636} 3637 3638 3639is.bell.uqo <- 3640is.bell.qrrvglm <- function(object, ...) { 3641 is.finite(Max(object, ...)) 3642} 3643 3644 3645is.bell.rrvgam <- function(object, ...) { 3646 NA * Max(object, ...) 3647} 3648 3649 3650 if (!isGeneric("is.bell")) 3651 setGeneric("is.bell", 3652 function(object, ...) standardGeneric("is.bell")) 3653setMethod("is.bell","qrrvglm", 3654 function(object,...) is.bell.qrrvglm(object,...)) 3655setMethod("is.bell","rrvglm", 3656 function(object, ...) is.bell.rrvglm(object, ...)) 3657setMethod("is.bell","vlm", 3658 function(object, ...) is.bell.vlm(object, ...)) 3659setMethod("is.bell","rrvgam", 3660 function(object, ...) is.bell.rrvgam(object, ...)) 3661setMethod("is.bell","Coef.qrrvglm", 3662 function(object,...) is.bell.qrrvglm(object,...)) 3663 3664 3665 3666 3667 3668 3669 if (!isGeneric("Rank")) 3670 setGeneric("Rank", 3671 function(object, ...) standardGeneric("Rank"), package = "VGAM") 3672 3673setMethod("Rank", "rrvglm", 3674 function(object, ...) Rank.rrvglm(object, ...)) 3675setMethod("Rank", "qrrvglm", 3676 function(object, ...) Rank.qrrvglm(object, ...)) 3677setMethod("Rank", "rrvgam", 3678 function(object, ...) Rank.rrvgam(object, ...)) 3679 3680 3681 3682 3683 3684 3685 3686