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 10 11 12 13 14 15 16 17 18 19 rcim <- 20 function(y, 21 family = poissonff, 22 Rank = 0, 23 M1 = NULL, 24 weights = NULL, 25 which.linpred = 1, 26 Index.corner = ifelse(is.null(str0), 0, max(str0)) + 1:Rank, 27 rprefix = "Row.", 28 cprefix = "Col.", 29 iprefix = "X2.", 30 offset = 0, 31 32 str0 = if (Rank) 1 else NULL, # Ignored if Rank == 0 33 summary.arg = FALSE, h.step = 0.0001, 34 rbaseline = 1, cbaseline = 1, 35 36 has.intercept = TRUE, 37 38 M = NULL, 39 40 rindex = 2:nrow(y), # Row index 41 cindex = 2:ncol(y), # Col index 42 iindex = 2:nrow(y), # Interaction index 43 44 45 46 ...) { 47 48 49 50 51 52 rindex <- unique(sort(rindex)) 53 cindex <- unique(sort(cindex)) 54 iindex <- unique(sort(iindex)) 55 56 57 if (Rank == 0 && !has.intercept) 58 warning("probably 'has.intercept == TRUE' is better for ", 59 "a rank-0 model") 60 61 62 63 ncoly <- ncol(y) 64 65 66 noroweffects <- FALSE 67 nocoleffects <- FALSE 68 69 if (!is.Numeric(which.linpred, length.arg = 1, 70 integer.valued = TRUE, positive = TRUE)) 71 stop("bad input for argument 'which.linpred'") 72 73 if (!is.character(rprefix)) 74 stop("argument 'rprefix' must be character") 75 if (!is.character(cprefix)) 76 stop("argument 'cprefix' must be character") 77 78 if (is.character(family)) 79 family <- get(family) 80 if (is.function(family)) 81 family <- ((family)()) 82 if (!inherits(family, "vglmff")) { 83 stop("'family = ", family, "' is not a VGAM family function") 84 } 85 efamily <- family 86 87 88 if (!is.Numeric(M1)) { 89 iefamily <- efamily@infos 90 91 if (is.function(iefamily)) 92 M1 <- (iefamily())$M1 93 if (is.Numeric(M1)) 94 M1 <- abs(M1) 95 } 96 if (!is.Numeric(M1)) { 97 if (!is.Numeric(M)) 98 warning("cannot determine the value of 'M1'.", 99 "Assuming the value one.") 100 M1 <- 1 101 } 102 103 104 105 M <- if (is.null(M)) M1 * ncol(y) else M 106 107 special <- (M > 1) && (M1 == 1) 108 109 110 111 object.save <- y 112 y <- if (is(y, "rrvglm")) { 113 depvar(object.save) 114 } else { 115 as(as.matrix(y), "matrix") 116 } 117 if (length(dim(y)) != 2 || nrow(y) < 3 || ncol(y) < 3) 118 stop("argument 'y' must be a matrix with >= 3 rows & columns, or ", 119 "a rrvglm() object") 120 121 122 .rcim.df <- 123 if (!noroweffects) 124 data.frame("Row.2" = I.col(2, nrow(y))) else # See below 125 if (!nocoleffects) 126 data.frame("Col.2" = I.col(2, nrow(y))) else # See below 127 stop("at least one of 'noroweffects' and 'nocoleffects' ", 128 "must be FALSE") 129 130 131 min.row.val <- rindex[1] # == min(rindex) since its sorted # Usually 2 132 min.col.val <- cindex[1] # == min(cindex) since its sorted # Usually 2 133 if (!noroweffects) { 134 colnames( .rcim.df ) <- 135 paste(rprefix, as.character(min.row.val), # "2", 136 sep = "") # Overwrite "Row.2" 137 } else if (!nocoleffects) { 138 colnames( .rcim.df ) <- 139 paste(cprefix, as.character(min.col.val), # "2", 140 sep = "") # Overwrite "Col.2" 141 } 142 143 144 145 yn1 <- if (length(dimnames(y)[[1]])) dimnames(y)[[1]] else 146 param.names(iprefix, nrow(y)) 147 warn.save <- options()$warn 148 options(warn = -3) # Suppress the warnings (hopefully, temporarily) 149 if (any(!is.na(as.numeric(substring(yn1, 1, 1))))) 150 yn1 <- param.names(iprefix, nrow(y)) 151 options(warn = warn.save) 152 153 154 nrprefix <- as.name(rprefix) 155 ncprefix <- as.name(cprefix) 156 157 158 assign(rprefix, factor(1:nrow(y))) 159 modmat.row <- substitute( 160 model.matrix( ~ .rprefix ), list( .rprefix = nrprefix )) 161 162 LLL <- ifelse(special, M, ncol(y)) 163 assign(cprefix, factor(1:LLL)) 164 165 166 modmat.col <- substitute( 167 model.matrix( ~ .cprefix ), list( .cprefix = ncprefix )) 168 modmat.row <- eval( modmat.row ) 169 modmat.col <- eval( modmat.col ) 170 171 172 173 174 Hlist <- 175 if (has.intercept) { 176 list("(Intercept)" = matrix(1, LLL, 1)) 177 } else { 178 temp <- list("Row.2" = matrix(1, LLL, 1)) # Overwrite this name: 179 names(temp) <- paste(rprefix, as.character(min.row.val), sep = "") 180 temp 181 } 182 183 184 if (!noroweffects) 185 for (ii in rindex) { 186 Hlist[[paste(rprefix, ii, sep = "")]] <- matrix(1, LLL, 1) 187 .rcim.df[[paste(rprefix, ii, sep = "")]] <- modmat.row[, ii] 188 } 189 190 191 if (!nocoleffects) 192 for (ii in cindex) { 193 temp6.mat <- modmat.col[, ii, drop = FALSE] 194 Hlist[[paste(cprefix, ii, sep = "")]] <- temp6.mat 195 .rcim.df[[paste(cprefix, ii, sep = "")]] <- rep_len(1, nrow(y)) 196 } 197 198 199 if (Rank > 0) { 200 for (ii in iindex) { 201 202 Hlist[[yn1[ii]]] <- diag(LLL) 203 204 .rcim.df[[yn1[ii]]] <- I.col(ii, nrow(y)) 205 } 206 } 207 208 209 dimnames(.rcim.df) <- list(if (length(dimnames(y)[[1]])) 210 dimnames(y)[[1]] else 211 as.character(iindex), 212 dimnames( .rcim.df )[[2]]) 213 214 str1 <- paste(if (has.intercept) "~ 1 + " else "~ -1 + ", rprefix, 215 as.character(min.row.val), # "2", 216 sep = "") 217 218 219 if (nrow(y) > 2) 220 str1 <- paste(str1, 221 paste(rprefix, rindex[-1], sep = "", collapse = " + "), 222 sep = " + ") 223 224 225 str1 <- paste(str1, 226 paste(cprefix, cindex, sep = "", collapse = " + "), 227 sep = " + ") 228 229 230 231 232 str2 <- paste("y", str1) 233 if (Rank > 0) { 234 str2 <- paste(str2, 235 paste(yn1[iindex], sep = "", collapse = " + "), 236 sep = " + ") 237 } 238 239 240 241 242 controlfun <- if (Rank == 0) vglm.control else rrvglm.control # orig. 243 244 245 mycontrol <- controlfun(Rank = Rank, 246 Index.corner = Index.corner, 247 str0 = str0, ...) 248 249 if (mycontrol$trace) { 250 } 251 252 253 254 if ((mindim <- min(nrow(y), ncol(y))) <= Rank) { 255 stop("argument 'Rank' is too high. Must be a value from 0 ", 256 "to ", mindim - 1, " inclusive") 257 } 258 259 260 261 if (Rank > 0) 262 mycontrol$noRRR <- as.formula(str1) # Overwrite this 263 264 assign(".rcim.df", .rcim.df, envir = VGAM::VGAMenv) 265 266 warn.save <- options()$warn 267 options(warn = -3) # Suppress the warnings (hopefully, temporarily) 268 269 if (mycontrol$trace) { 270 } 271 272 273 if (M1 > 1) { 274 orig.Hlist <- Hlist 275 276 kmat1 <- matrix(0, nrow = M1, ncol = 1) 277 kmat1[which.linpred, 1] <- 1 278 kmat0 <- (diag(M1))[, -which.linpred, drop = FALSE] 279 280 for (ii in seq_along(Hlist)) { 281 Hlist[[ii]] <- kronecker(Hlist[[ii]], kmat1) 282 } 283 if (has.intercept) 284 Hlist[["(Intercept)"]] <- cbind(Hlist[["(Intercept)"]], 285 kronecker(matrix(1, ncoly, 1), 286 kmat0)) 287 288 289 if (mycontrol$trace) { 290 } 291 } 292 293 294 295 offset.matrix <- 296 matrix(offset, nrow = nrow(y), 297 ncol = M) # byrow = TRUE 298 299 300 301 answer <- if (Rank > 0) { 302 if (is(object.save, "rrvglm")) object.save else 303 rrvglm(as.formula(str2), 304 family = family, 305 constraints = Hlist, 306 offset = offset.matrix, 307 weights = if (length(weights)) 308 weights else rep_len(1, nrow(y)), 309 ..., 310 control = mycontrol, data = .rcim.df ) 311 } else { 312 if (is(object.save, "vglm")) object.save else 313 vglm(as.formula(str2), 314 family = family, 315 constraints = Hlist, 316 offset = offset.matrix, 317 weights = if (length(weights)) 318 weights else rep_len(1, nrow(y)), 319 ..., 320 control = mycontrol, data = .rcim.df ) 321 } 322 323 options(warn = warn.save) # Restore warnings back to prior state 324 325 326 answer <- if (summary.arg) { 327 if (Rank > 0) { 328 summary.rrvglm(as(answer, "rrvglm"), h.step = h.step) 329 } else { 330 summary(answer) 331 } 332 } else { 333 as(answer, ifelse(Rank > 0, "rcim", "rcim0")) 334 } 335 336 337 answer@misc$rbaseline <- rbaseline 338 answer@misc$cbaseline <- cbaseline 339 answer@misc$which.linpred <- which.linpred 340 answer@misc$offset <- offset.matrix 341 342 answer 343} 344 345 346 347 348 349 350 351 352summaryrcim <- function(object, ...) { 353 rcim(depvar(object), summary.arg = TRUE, ...) 354} 355 356 357 358 359 360 361 362 363 364 setClass("rcim0", representation(not.needed = "numeric"), 365 contains = "vglm") # Added 20110506 366 367 setClass("rcim", representation(not.needed = "numeric"), 368 contains = "rrvglm") 369 370 371setMethod("summary", "rcim0", 372 function(object, ...) 373 summaryrcim(object, ...)) 374 375 376setMethod("summary", "rcim", 377 function(object, ...) 378 summaryrcim(object, ...)) 379 380 381 382 383 384 385 386 387 388 389 Rcim <- function(mat, rbaseline = 1, cbaseline = 1) { 390 391 mat <- as.matrix(mat) 392 RRR <- dim(mat)[1] 393 CCC <- dim(mat)[2] 394 395 rnames <- if (is.null(rownames(mat))) { 396 param.names("X", RRR) 397 } else { 398 rownames(mat) 399 } 400 401 cnames <- if (is.null(colnames(mat))) { 402 param.names("Y", CCC) 403 } else { 404 colnames(mat) 405 } 406 407 r.index <- if (is.character(rbaseline)) 408 which(rownames(mat) == rbaseline) else 409 if (is.numeric(rbaseline)) rbaseline else 410 stop("argement 'rbaseline' must be numeric", 411 "or character of the level of row") 412 413 c.index <- if (is.character(cbaseline)) 414 which(colnames(mat) == cbaseline) else 415 if (is.numeric(cbaseline)) cbaseline else 416 stop("argement 'cbaseline' must be numeric", 417 "or character of the level of row") 418 419 if (length(r.index) != 1) 420 stop("Could not match with argument 'rbaseline'") 421 422 if (length(c.index) != 1) 423 stop("Could not match with argument 'cbaseline'") 424 425 426 yswap <- rbind(mat[r.index:RRR, ], 427 if (r.index > 1) mat[1:(r.index - 1),] else NULL) 428 yswap <- cbind(yswap[, c.index:CCC], 429 if (c.index > 1) yswap[, 1:(c.index - 1)] else NULL) 430 431 new.rnames <- rnames[c(r.index:RRR, 432 if (r.index > 1) 1:(r.index - 1) else NULL)] 433 new.cnames <- cnames[c(c.index:CCC, 434 if (c.index > 1) 1:(c.index - 1) else NULL)] 435 colnames(yswap) <- new.cnames 436 rownames(yswap) <- new.rnames 437 438 yswap 439} 440 441 442 443 444 445 446 447 448 449 450 451 452 453 plotrcim0 <- function(object, 454 centered = TRUE, which.plots = c(1, 2), 455 hline0 = TRUE, hlty = "dashed", hcol = par()$col, hlwd = par()$lwd, 456 rfirst = 1, cfirst = 1, 457 rtype = "h", ctype = "h", 458 rcex.lab = 1, rcex.axis = 1, # rlabels = FALSE, 459 rtick = FALSE, 460 ccex.lab = 1, ccex.axis = 1, # clabels = FALSE, 461 ctick = FALSE, 462 rmain = "Row effects", rsub = "", 463 rxlab = "", rylab = "Row effects", 464 cmain = "Column effects", csub = "", 465 cxlab = "", cylab = "Column effects", 466 rcol = par()$col, ccol = par()$col, 467 no.warning = FALSE, 468 ...) { 469 470 471 nparff <- if (is.numeric(object@family@infos()$M1)) { 472 object@family@infos()$M1 473 } else { 474 1 475 } 476 477 478 479 if (!no.warning && 480 is.numeric(object@control$Rank) && 481 object@control$Rank != 0) 482 warning("argument 'object' is not Rank-0") 483 484 485 n.lm <- nrow(object@y) 486 487 cobj <- coefficients(object) 488 489 upperbound <- if (!is.numeric(object@control$Rank) || 490 object@control$Rank == 0) length(cobj) else 491 length(object@control$colx1.index) 492 493 orig.roweff <- c("Row.1" = 0, cobj[(nparff + 1) : (nparff + n.lm - 1)]) 494 orig.coleff <- c("Col.1" = 0, cobj[(nparff + n.lm) : upperbound]) 495 last.r <- length(orig.roweff) 496 last.c <- length(orig.coleff) 497 498 499 orig.raxisl <- rownames(object@y) 500 orig.caxisl <- colnames(object@y) 501 if (is.null(orig.raxisl)) 502 orig.raxisl <- as.character(1:nrow(object@y)) 503 if (is.null(orig.caxisl)) 504 orig.caxisl <- as.character(1:ncol(object@y)) 505 506 roweff.orig <- 507 roweff <- orig.roweff[c(rfirst:last.r, 508 if (rfirst > 1) 1:(rfirst-1) else NULL)] 509 coleff.orig <- 510 coleff <- orig.coleff[c(cfirst:last.c, 511 if (cfirst > 1) 1:(cfirst-1) else NULL)] 512 513 if (centered) { 514 roweff <- scale(roweff, scale = FALSE) # Center it only 515 coleff <- scale(coleff, scale = FALSE) # Center it only 516 } 517 518 raxisl <- orig.raxisl[c(rfirst:last.r, 519 if (rfirst > 1) 1:(rfirst-1) else NULL)] 520 521 caxisl <- orig.caxisl[c(cfirst:last.c, 522 if (cfirst > 1) 1:(cfirst-1) else NULL)] 523 524 525 if (any(which.plots == 1, na.rm = TRUE)) { 526 plot(roweff, type = rtype, 527 axes = FALSE, col = rcol, main = rmain, 528 sub = rsub, xlab = rxlab, ylab = rylab, ...) 529 530 axis(1, at = seq_along(raxisl), 531 cex.lab = rcex.lab, 532 cex.axis = rcex.axis, 533 labels = raxisl) 534 axis(2, cex.lab = rcex.lab, ...) # las = rlas) 535 536 if (hline0) 537 abline(h = 0, lty = hlty, col = hcol, lwd = hlwd) 538 } 539 540 541 if (any(which.plots == 2, na.rm = TRUE)) { 542 plot(coleff, type = ctype, 543 axes = FALSE, col = ccol, main = cmain, # lwd = 2, xpd = FALSE, 544 sub = csub, xlab = cxlab, ylab = cylab, ...) 545 546 axis(1, at = seq_along(caxisl), 547 cex.lab = ccex.lab, 548 cex.axis = ccex.axis, 549 labels = caxisl) 550 axis(2, cex.lab = ccex.lab, ...) # las = clas) 551 552 if (hline0) 553 abline(h = 0, lty = hlty, col = hcol, lwd = hlwd) 554 } 555 556 557 558 559 object@post$row.effects = roweff 560 object@post$col.effects = coleff 561 object@post$raw.row.effects = roweff.orig 562 object@post$raw.col.effects = coleff.orig 563 564 invisible(object) 565} 566 567 568 569 570 571setMethod("plot", "rcim0", 572 function(x, y, ...) 573 plotrcim0(object = x, ...)) 574 575 576setMethod("plot", "rcim", 577 function(x, y, ...) 578 plotrcim0(object = x, ...)) 579 580 581 582 583 584 585 586 587 588 589 590 591 592 moffset <- 593 function(mat, roffset = 0, coffset = 0, postfix = "", 594 rprefix = "Row.", 595 cprefix = "Col." 596 ) { 597 598 599 600 601 602 if ((is.numeric(roffset) && (roffset == 0)) && 603 (is.numeric(coffset) && (coffset == 0))) 604 return(mat) 605 606 607 vecmat <- c(unlist(mat)) 608 ind1 <- if (is.character(roffset)) 609 which(rownames(mat) == roffset) else 610 if (is.numeric(roffset)) roffset + 1 else 611 stop("argument 'roffset' not matched (character).", 612 " It must be numeric, ", 613 "else character and match the ", 614 "row names of the response") 615 ind2 <- if (is.character(coffset)) 616 which(colnames(mat) == coffset) else 617 if (is.numeric(coffset)) coffset + 1 else 618 stop("argument 'coffset' not matched (character).", 619 " It must be numeric, ", 620 "else character and match the ", 621 "column names of the response") 622 623 if (!is.Numeric(ind1, positive = TRUE, 624 integer.valued = TRUE, length.arg = 1) || 625 !is.Numeric(ind2, positive = TRUE, 626 integer.valued = TRUE, length.arg = 1)) 627 stop("bad input for arguments 'roffset' and/or 'coffset'") 628 if (ind1 > nrow(mat)) 629 stop("too large a value for argument 'roffset'") 630 if (ind2 > ncol(mat)) 631 stop("too large a value for argument 'coffset'") 632 633 634 start.ind <- (ind2 - 1)* nrow(mat) + ind1 635 636 637 svecmat <- vecmat[c(start.ind:(nrow(mat) * ncol(mat)), 638 0:(start.ind - 1))] 639 640 rownames.mat <- rownames(mat) 641 if (length(rownames.mat) != nrow(mat)) 642 rownames.mat <- param.names(rprefix, nrow(mat)) 643 644 colnames.mat <- colnames(mat) 645 if (length(colnames.mat) != ncol(mat)) 646 colnames.mat <- param.names(cprefix, ncol(mat)) 647 648 649 newrn <- if (roffset > 0) 650 c(rownames.mat[c(ind1:nrow(mat))], 651 paste(rownames.mat[0:(ind1-1)], postfix, sep = "")) else 652 rownames.mat 653 654 newcn <- c(colnames.mat[c(ind2:ncol(mat), 0:(ind2 - 1))]) 655 if (roffset > 0) 656 newcn <- paste(newcn, postfix, sep = "") 657 658 newmat <- matrix(svecmat, nrow(mat), ncol(mat), 659 dimnames = list(newrn, newcn)) 660 newmat 661} 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681Confint.rrnb <- function(rrnb2, level = 0.95) { 682 683 if (class(rrnb2) != "rrvglm") 684 stop("argument 'rrnb2' does not appear to be a rrvglm() object") 685 686 if (!any(rrnb2@family@vfamily == "negbinomial")) 687 stop("argument 'rrnb2' does not appear to be a negbinomial() fit") 688 689 if (rrnb2@control$Rank != 1) 690 stop("argument 'rrnb2' is not Rank-1") 691 692 if (rrnb2@misc$M != 2) 693 stop("argument 'rrnb2' does not have M = 2") 694 695 if (!all(rrnb2@misc$link == "loglink")) 696 stop("argument 'rrnb2' does not have log links for both parameters") 697 698 a21.hat <- (Coef(rrnb2)@A)["loglink(size)", 1] 699 beta11.hat <- Coef(rrnb2)@B1["(Intercept)", "loglink(mu)"] 700 beta21.hat <- Coef(rrnb2)@B1["(Intercept)", "loglink(size)"] 701 delta1.hat <- exp(a21.hat * beta11.hat - beta21.hat) 702 delta2.hat <- 2 - a21.hat 703 704 SE.a21.hat <- sqrt(vcovrrvglm(rrnb2)["I(latvar.mat)", "I(latvar.mat)"]) 705 706 707 ci.a21 <- a21.hat + c(-1, 1) * qnorm(1 - (1 - level)/2) * SE.a21.hat 708 (ci.delta2 <- 2 - rev(ci.a21)) # e.g., the 95 percent CI 709 710 list(a21.hat = a21.hat, 711 beta11.hat = beta11.hat, 712 beta21.hat = beta21.hat, 713 CI.a21 = ci.a21, 714 CI.delta2 = ci.delta2, 715 delta1 = delta1.hat, 716 delta2 = delta2.hat, 717 SE.a21.hat = SE.a21.hat) 718} 719 720 721 722 723 724Confint.nb1 <- function(nb1, level = 0.95) { 725 726 727 728 729 if (class(nb1) != "vglm") 730 stop("argument 'nb1' does not appear to be a vglm() object") 731 732 if (!any(nb1@family@vfamily == "negbinomial")) 733 stop("argument 'nb1' does not appear to be a negbinomial() fit") 734 735 if (!all(unlist(constraints(nb1)[-1]) == 1)) 736 stop("argument 'nb1' does not appear to have 'parallel = TRUE'") 737 738 if (!all(unlist(constraints(nb1)[1]) == c(diag(nb1@misc$M)))) 739 stop("argument 'nb1' does not have 'parallel = FALSE' ", 740 "for the intercept") 741 742 if (nb1@misc$M != 2) 743 stop("argument 'nb1' does not have M = 2") 744 745 if (!all(nb1@misc$link == "loglink")) 746 stop("argument 'nb1' does not have log links for both parameters") 747 748 cnb1 <- coefficients(as(nb1, "vglm"), matrix = TRUE) 749 mydiff <- (cnb1["(Intercept)", "loglink(size)"] - 750 cnb1["(Intercept)", "loglink(mu)"]) 751 delta0.hat <- exp(mydiff) 752 (phi0.hat <- 1 + 1 / delta0.hat) # MLE of phi0 753 754 755 myvcov <- vcov(as(nb1, "vglm")) # Not great; improve this! 756 757 758 759 myvec <- cbind(c(-1, 1, rep_len(0, nrow(myvcov) - 2))) 760 se.mydiff <- c(sqrt(t(myvec) %*% myvcov %*% myvec)) # c() to drop the 2-D 761 762 ci.mydiff <- mydiff + c(-1, 1) * qnorm(1 - (1 - level)/2) * se.mydiff 763 764 ci.delta0 <- ci.exp.mydiff <- exp(ci.mydiff) 765 (ci.phi0 <- 1 + 1 / rev(ci.delta0)) # e.g., the 95 percent CI for phi0 766 767 list(CI.phi0 = ci.phi0, 768 CI.delta0 = ci.delta0, 769 delta0 = delta0.hat, 770 phi0 = phi0.hat) 771} 772 773 774 775 776 777 778plota21 <- function(rrvglm2, show.plot = TRUE, nseq.a21 = 31, 779 se.eachway = c(5, 5), # == c(LHS, RHS), 780 trace.arg = TRUE, 781 lwd = 2, ...) { 782 783 784 785 786 787 if (class(rrvglm2) != "rrvglm") 788 stop("argument 'rrvglm2' does not appear to be a rrvglm() object") 789 790 if (rrvglm2@control$Rank != 1) 791 stop("argument 'rrvglm2' is not Rank-1") 792 793 if (rrvglm2@misc$M != 2) 794 stop("argument 'rrvglm2' does not have M = 2") 795 796 797 loglik.orig <- logLik(rrvglm2) 798 temp1 <- Confint.rrnb(rrvglm2) # zz 799 800 a21.hat <- (Coef(rrvglm2)@A)[2, 1] 801 SE.a21.hat <- temp1$SE.a21.hat 802 803 804 SE.a21.hat <- sqrt(vcov(rrvglm2)["I(latvar.mat)", "I(latvar.mat)"]) 805 806 807 big.ci.a21 <- a21.hat + c(-1, 1) * se.eachway * SE.a21.hat 808 seq.a21 <- seq(big.ci.a21[1], big.ci.a21[2], length = nseq.a21) 809 Hlist.orig <- constraints.vlm(rrvglm2, type = "term") 810 811 812 alreadyComputed <- !is.null(rrvglm2@post$a21.matrix) 813 814 815 a21.matrix <- if (alreadyComputed) rrvglm2@post$a21.matrix else 816 cbind(a21 = seq.a21, loglikelihood = 0) 817 prev.etastart <- predict(rrvglm2) # Halves the computing time 818 funname <- "vglm" 819 listcall <- as.list(rrvglm2@call) 820 821 822 if (!alreadyComputed) 823 for (ii in 1:nseq.a21) { 824 if (trace.arg) 825 print(ii) 826 argslist <- vector("list", length(listcall) - 1) 827 for (kay in 2:(length(listcall))) 828 argslist[[kay - 1]] <- listcall[[kay]] 829 830 names(argslist) <- c(names(listcall)[-1]) 831 832 argslist$trace <- trace.arg 833 argslist$etastart <- prev.etastart 834 argslist$constraints <- Hlist.orig 835 836 837 for (kay in 2:length(argslist[["constraints"]])) { 838 argslist[["constraints"]][[kay]] <- rbind(1, a21.matrix[ii, 1]) 839 } 840 841 842 fitnew <- do.call(what = funname, args = argslist) 843 844 a21.matrix[ii, 2] <- logLik(fitnew) 845 846 prev.etastart <- predict(fitnew) 847 } 848 849 850 851 if (show.plot) { 852 plot(a21.matrix[ ,1], a21.matrix[ ,2], type = "l", 853 col = "blue", cex.lab = 1.1, 854 xlab = expression(a[21]), ylab = "Log-likelihood") # ... 855 856 abline(v = (Hlist.orig[[length(Hlist.orig)]])[2, 1], 857 col = "darkorange", lty = "dashed") 858 859 abline(h = loglik.orig, 860 col = "darkorange", lty = "dashed") 861 862 abline(h = loglik.orig - 863 qchisq(0.95, df = 1) / 2, 864 col = "darkorange", lty = "dashed") 865 866 abline(v = a21.hat + c(-1, 1) * 1.96 * SE.a21.hat, 867 col = "gray50", lty = "dashed", lwd = lwd) 868 869 } # End of (show.plot) 870 871 rrvglm2@post <- list(a21.matrix = a21.matrix) 872 invisible(rrvglm2) 873} 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 Qvar <- function(object, 889 factorname = NULL, 890 which.linpred = 1, 891 coef.indices = NULL, 892 labels = NULL, dispersion = NULL, 893 reference.name = "(reference)", 894 estimates = NULL 895 ) { 896 897 898 899 900 901 902 903 904 if (!is.Numeric(which.linpred, length.arg = 1, 905 integer.valued = TRUE, positive = TRUE)) 906 stop("argument 'which.linpred' must be a positive integer") 907 908 909 910 coef.indices.saved <- coef.indices 911 912 if (!is.matrix(object)) { 913 model <- object 914 if (is.null(factorname) && is.null(coef.indices)) { 915 stop("arguments 'factorname' and 'coef.indices' are ", 916 "both NULL") 917 } 918 919 920 if (is.null(coef.indices)) { 921 922 M <- npred(model) 923 if (M < which.linpred) 924 stop("argument 'which.linpred' must be a value from the set 1:", 925 M) 926 927 928 newfactorname <- if (M > 1) { 929 clist <- constraints(model, type = "term") 930 931 Hk <- clist[[factorname]] 932 Mdot <- ncol(Hk) 933 Hk.row <- Hk[which.linpred, ] 934 if (sum(Hk.row != 0) > 1) 935 stop("cannot handle rows of constraint matrices with more ", 936 "than one nonzero value") 937 938 foo <- function(ii) 939 switch(as.character(ii), '1'="1st", '2'="2nd", '3'="3rd", 940 paste(ii, "th", sep = "")) 941 if (sum(Hk.row != 0) == 0) 942 stop("factor '", factorname, "' is not used the ", 943 foo(which.linpred), " eta (linear predictor)") 944 945 row.index <- (1:Mdot)[Hk.row != 0] 946 947 all.labels <- vlabel(factorname, ncolHlist = Mdot, M = M) 948 all.labels[row.index] 949 } else { 950 factorname 951 } 952 953 colptr <- attr(model.matrix(object, type = "vlm"), "vassign") 954 colptr <- if (M > 1) { 955 colptr[newfactorname] 956 } else { 957 colptr[[newfactorname]] 958 } 959 960 coef.indices <- colptr 961 962 contmat <- if (length(model@xlevels[[factorname]]) == 963 length(coef.indices)) { 964 diag(length(coef.indices)) 965 } else { 966 eval(call(model@contrasts[[factorname]], 967 model@xlevels [[factorname]])) 968 } 969 rownames(contmat) <- model@xlevels[[factorname]] 970 971 if (is.null(estimates)) { 972 if (M > 1) { 973 estimates <- matrix(-1, nrow(contmat), 1) # Used to be nc = Mdot 974 ii <- 1 975 estimates[, ii] <- contmat %*% 976 (coefvlm(model)[(coef.indices[[ii]])]) 977 } else { 978 estimates <- contmat %*% (coefvlm(model)[coef.indices]) 979 } 980 } 981 982 983 984 985 Covmat <- vcov(model, dispersion = dispersion) 986 987 988 989 covmat <- Covmat[unlist(coef.indices), 990 unlist(coef.indices), drop = FALSE] 991 covmat <- if (M > 1) { 992 ii <- 1 993 ans <- contmat %*% Covmat[(colptr[[ii]]), 994 (colptr[[ii]])] %*% t(contmat) 995 ans 996 } else { 997 contmat %*% covmat %*% t(contmat) 998 } 999 } else { # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 1000 kk <- length(coef.indices) 1001 refPos <- numeric(0) 1002 if (0 %in% coef.indices) { 1003 refPos <- which(coef.indices == 0) 1004 coef.indices <- coef.indices[-refPos] 1005 } 1006 1007 1008 1009 1010 covmat <- vcov(model, dispersion = dispersion) 1011 1012 1013 1014 1015 covmat <- covmat[coef.indices, coef.indices, drop = FALSE] 1016 1017 if (is.null(estimates)) 1018 estimates <- coefvlm(model)[coef.indices] 1019 1020 if (length(refPos) == 1) { 1021 if (length(estimates) != kk) 1022 estimates <- c(0, estimates) 1023 covmat <- rbind(0, cbind(0, covmat)) 1024 names(estimates)[1] <- 1025 rownames(covmat)[1] <- 1026 colnames(covmat)[1] <- reference.name 1027 if (refPos != 1) { 1028 perm <- if (refPos == kk) c(2:kk, 1) else 1029 c(2:refPos, 1, (refPos + 1):kk) 1030 estimates <- estimates[perm] 1031 covmat <- covmat[perm, perm, drop = FALSE] 1032 } 1033 } 1034 } # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 1035 1036 1037 return(Recall(covmat, 1038 factorname = factorname, 1039 which.linpred = which.linpred, 1040 coef.indices = coef.indices.saved, 1041 labels = labels, 1042 dispersion = dispersion, 1043 estimates = estimates 1044 ) 1045 ) 1046 } else { # ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1047 1048 covmat <- object 1049 if (length(labels)) 1050 rownames(covmat) <- colnames(covmat) <- labels 1051 if ((LLL <- dim(covmat)[1]) <= 2) 1052 stop("This function works only for factors with 3 ", 1053 "or more levels") 1054 } # ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1055 1056 1057 1058 1059 allvcov <- covmat 1060 for (ilocal in 1:LLL) 1061 for (jlocal in ilocal:LLL) 1062 allvcov[ilocal, jlocal] <- 1063 allvcov[jlocal, ilocal] <- covmat[ilocal, ilocal] + 1064 covmat[jlocal, jlocal] - 1065 covmat[ilocal, jlocal] * 2 1066 1067 diag(allvcov) <- rep_len(1.0, LLL) # Any positive value should do 1068 1069 1070 wmat <- matrix(1.0, LLL, LLL) 1071 diag(wmat) <- sqrt( .Machine$double.eps ) 1072 1073 1074 logAllvcov <- log(allvcov) 1075 attr(logAllvcov, "Prior.Weights") <- wmat 1076 attr(logAllvcov, "estimates") <- estimates 1077 attr(logAllvcov, "coef.indices") <- coef.indices 1078 attr(logAllvcov, "factorname") <- factorname 1079 attr(logAllvcov, "regularVar") <- diag(covmat) 1080 attr(logAllvcov, "which.linpred") <- which.linpred 1081 1082 logAllvcov 1083} # End of Qvar() 1084 1085 1086 1087 1088 1089 1090 1091 1092WorstErrors <- function(qv.object) { 1093 stop("20110729; does not work") 1094 1095 reducedForm <- function(covmat, qvmat) { 1096 nlevels <- dim(covmat)[1] 1097 firstRow <- covmat[1, ] 1098 ones <- rep_len(1, nlevels) 1099 J <- outer(ones, ones) 1100 notzero <- 2:nlevels 1101 r.covmat <- covmat + (firstRow[1]*J) - 1102 outer(firstRow, ones) - 1103 outer(ones, firstRow) 1104 r.covmat <- r.covmat[notzero, notzero] 1105 qv1 <- qvmat[1, 1] 1106 r.qvmat <- (qvmat + qv1*J)[notzero, notzero] 1107 list(r.covmat, r.qvmat) 1108 } 1109 covmat <- qv.object$covmat 1110 qvmat <- diag(qv.object$qvframe$quasiVar) 1111 r.form <- reducedForm(covmat, qvmat) 1112 r.covmat <- r.form[[1]] 1113 r.qvmat <- r.form[[2]] 1114 inverse.sqrt <- solve(chol(r.covmat)) 1115 evalues <- eigen(t(inverse.sqrt) %*% r.qvmat %*% inverse.sqrt, 1116 symmetric = TRUE)$values 1117 sqrt(c(min(evalues), max(evalues))) - 1 1118} 1119 1120 1121 1122 1123IndentPrint <- function(object, indent = 4, ...) { 1124 stop("20110729; does not work") 1125 1126 zz <- "" 1127 tc <- textConnection("zz", "w", local = TRUE) 1128 sink(tc) 1129 try(print(object, ...)) 1130 sink() 1131 close(tc) 1132 indent <- paste(rep_len(" ", indent), sep = "", collapse = "") 1133 cat(paste(indent, zz, sep = ""), sep = "\n") 1134} 1135 1136 1137 1138Print.qv <- function(x, ...) { 1139 stop("20110729; does not work") 1140 1141} 1142 1143 1144 1145 1146summary.qvar <- function(object, ...) { 1147 1148 1149 relerrs <- 1 - sqrt(exp(residuals(object, type = "response"))) 1150 diag(relerrs) <- NA 1151 1152 minErrSimple <- round(100 * min(relerrs, na.rm = TRUE), 1) 1153 maxErrSimple <- round(100 * max(relerrs, na.rm = TRUE), 1) 1154 1155 1156 1157 estimates <- c(object@extra$attributes.y$estimates) 1158 if (!length(names(estimates)) && 1159 is.matrix(object@extra$attributes.y$estimates)) 1160 names( estimates) <- rownames(object@extra$attributes.y$estimates) 1161 if (!length(names(estimates))) 1162 names( estimates) <- param.names("Level", length(estimates)) 1163 1164 1165 regularVar <- c(object@extra$attributes.y$regularVar) 1166 QuasiVar <- exp(diag(fitted(object))) / 2 1167 QuasiSE <- sqrt(QuasiVar) 1168 1169 1170 structure(list(estimate = estimates, 1171 SE = sqrt(regularVar), 1172 minErrSimple = minErrSimple, 1173 maxErrSimple = maxErrSimple, 1174 quasiSE = QuasiSE, 1175 object = object, 1176 quasiVar = QuasiVar), 1177 class = "summary.qvar") 1178} 1179 1180 1181 1182 1183print.summary.qvar <- function(x, ...) { 1184 1185 object <- x$object 1186 minErrSimple <- x$minErrSimple 1187 maxErrSimple <- x$maxErrSimple 1188 1189 x$minErrSimple <- NULL 1190 x$maxErrSimple <- NULL 1191 x$object <- NULL 1192 1193 1194 if (length(cl <- object@call)) { 1195 cat("Call:\n") 1196 dput(cl) 1197 } 1198 1199 1200 facname <- c(object@extra$attributes.y$factorname) 1201 if (length(facname)) 1202 cat("Factor name: ", facname, "\n") 1203 1204 1205 if (length(object@dispersion)) 1206 cat("\nDispersion: ", object@dispersion, "\n\n") 1207 1208 x <- as.data.frame(c(x)) 1209 print.data.frame(x) 1210 1211 1212 cat("\nWorst relative errors in SEs of simple contrasts (%): ", 1213 minErrSimple, ", ", maxErrSimple, "\n") 1214 1215 invisible(x) 1216} 1217 1218 1219 1220 1221qvar <- function(object, se = FALSE, ...) { 1222 1223 1224 1225 if (!inherits(object, "rcim") && !inherits(object, "rcim0")) 1226 warning("argument 'object' should be an 'rcim' or 'rcim0' object. ", 1227 "This call may fail.") 1228 1229 if (!(object@family@vfamily %in% c("uninormal", "normal1"))) 1230 warning("argument 'object' does not seem to have used ", 1231 "a 'uninormal' family.") 1232 1233 if (!any(object@misc$link == "explink")) 1234 warning("argument 'object' does not seem to have used ", 1235 "a 'explink' link function.") 1236 1237 quasiVar <- diag(predict(object)[, c(TRUE, FALSE)]) / 2 1238 if (se) sqrt(quasiVar) else quasiVar 1239} 1240 1241 1242 1243 1244 1245 1246 1247 1248plotqvar <- 1249qvplot <- function(object, 1250 interval.width = 2, 1251 ylab = "Estimate", 1252 xlab = NULL, # x$factorname, 1253 ylim = NULL, 1254 main = "", 1255 level.names = NULL, 1256 conf.level = 0.95, 1257 warn.ratio = 10, 1258 border = "transparent", # None 1259 points.arg = TRUE, 1260 length.arrows = 0.25, angle = 30, 1261 lwd = par()$lwd, 1262 scol = par()$col, 1263 slwd = par()$lwd, 1264 slty = par()$lty, 1265 ...) { 1266 1267 1268 1269 1270 if (!is.numeric(interval.width) && 1271 !is.numeric(conf.level)) 1272 stop("at least one of arguments 'interval.width' and 'conf.level' ", 1273 "should be numeric") 1274 1275 1276 1277 1278 1279 if (!any("uninormal" %in% object@family@vfamily)) 1280 stop("argument 'object' dos not appear to be a ", 1281 "rcim(, uninormal) object") 1282 1283 estimates <- c(object@extra$attributes.y$estimates) 1284 if (!length(names(estimates)) && 1285 is.matrix(object@extra$attributes.y$estimates)) 1286 names(estimates) <- rownames(object@extra$attributes.y$estimates) 1287 1288 1289 if (length(level.names) == length(estimates)) { 1290 names(estimates) <- level.names 1291 } else if (!length(names(estimates))) 1292 names(estimates) <- param.names("Level", length(estimates)) 1293 1294 1295 1296 1297 QuasiVar <- exp(diag(fitted(object))) / 2 1298 QuasiSE <- sqrt(QuasiVar) 1299 1300 if (!is.numeric(estimates)) 1301 stop("Cannot plot, because there are no 'proper' ", 1302 "parameter estimates") 1303 if (!is.numeric(QuasiSE)) 1304 stop("Cannot plot, because there are no ", 1305 "quasi standard errors") 1306 1307 1308 1309 faclevels <- factor(names(estimates), levels = names(estimates)) 1310 1311 1312 xvalues <- seq(along = faclevels) 1313 tops <- estimates + interval.width * QuasiSE 1314 tails <- estimates - interval.width * QuasiSE 1315 1316 1317 1318 1319 if (is.numeric(conf.level)) { 1320 zedd <- abs(qnorm((1 - conf.level) / 2)) 1321 lsd.tops <- estimates + zedd * QuasiSE / sqrt(2) 1322 lsd.tails <- estimates - zedd * QuasiSE / sqrt(2) 1323 if (max(QuasiSE) / min(QuasiSE) > warn.ratio) 1324 warning("Quasi SEs appear to be quite different... the ", 1325 "LSD intervals may not be very accurate") 1326 } else { 1327 lsd.tops <- NULL 1328 lsd.tails <- NULL 1329 } 1330 1331 1332 1333 1334 if (is.null(ylim)) 1335 ylim <- range(c(tails, tops, lsd.tails, lsd.tops), 1336 na.rm = TRUE) 1337 1338 if (is.null(xlab)) 1339 xlab <- "Factor level" 1340 1341 plot(faclevels, estimates, 1342 border = border, 1343 ylim = ylim, xlab = xlab, ylab = ylab, 1344 lwd = lwd, 1345 main = main, ...) 1346 1347 1348 if (points.arg) 1349 points(estimates, ...) 1350 1351 1352 if (is.numeric(interval.width)) { 1353 segments(xvalues, tails, xvalues, tops, 1354 col = scol, lty = slty, lwd = slwd) 1355 } 1356 1357 1358 if (is.numeric(conf.level)) { 1359 arrows(xvalues, lsd.tails, xvalues, lsd.tops, 1360 col = scol, lty = slty, lwd = slwd, code = 3, 1361 length = length.arrows, angle = angle) 1362 1363 } 1364 1365 1366 1367 1368 if (any(slotNames(object) == "post")) { 1369 object@post$estimates <- estimates 1370 object@post$xvalues <- xvalues 1371 if (is.numeric(interval.width)) { 1372 object@post$tails <- tails 1373 object@post$tops <- tops 1374 } 1375 if (is.numeric(conf.level)) { 1376 object@post$lsd.tails <- lsd.tails 1377 object@post$lsd.tops <- lsd.tops 1378 } 1379 } 1380 1381 1382 invisible(object) 1383} 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398