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 16summaryvglm <- 17 function(object, correlation = FALSE, 18 dispersion = NULL, digits = NULL, 19 presid = FALSE, # TRUE, 20 HDEtest = TRUE, # Added 20180203 21 hde.NA = TRUE, 22 threshold.hde = 0.001, 23 signif.stars = getOption("show.signif.stars"), 24 nopredictors = FALSE, 25 lrt0.arg = FALSE, 26 score0.arg = FALSE, 27 wald0.arg = FALSE, 28 values0 = 0, 29 subset = NULL, 30 omit1s = TRUE, 31 ... # Added 20151211 32 ) { 33 34 missing.HDEtest <- missing(HDEtest) 35 36 37 38 39 40 41 42 43 44 if (length(dispersion) && 45 dispersion == 0 && 46 length(object@family@summary.dispersion) && 47 !object@family@summary.dispersion) { 48 stop("cannot use the general VGLM formula (based on a residual ", 49 "sum of squares) for computing the dispersion parameter") 50 } 51 52 53 54 stuff <- summaryvlm( 55 object, 56 57 presid = FALSE, 58 59 correlation = correlation, 60 dispersion = dispersion, 61 62 lrt0.arg = lrt0.arg, 63 score0.arg = score0.arg, 64 wald0.arg = wald0.arg, 65 values0 = values0, 66 subset = subset, 67 omit1s = omit1s) 68 69 70 71 72 73 infos.fun <- object@family@infos 74 infos.list <- infos.fun() 75 summary.pvalues <- if (is.logical(infos.list$summary.pvalues)) 76 infos.list$summary.pvalues else TRUE 77 if (!summary.pvalues && ncol(stuff@coef3) == 4) 78 stuff@coef3 <- stuff@coef3[, -4] # Delete the pvalues column 79 80 81 82 83 84 85 86 87 88 answer <- 89 new("summary.vglm", 90 object, 91 coef3 = stuff@coef3, 92 coef4lrt0 = stuff@coef4lrt0, # Might be an empty "matrix" 93 coef4score0 = stuff@coef4score0, # Might be an empty "matrix" 94 coef4wald0 = stuff@coef4wald0, # Might be an empty "matrix" 95 cov.unscaled = stuff@cov.unscaled, 96 correlation = stuff@correlation, 97 df = stuff@df, 98 sigma = stuff@sigma) 99 100 101 if (presid) { 102 Presid <- resid(object, type = "pearson") 103 if (length(Presid)) 104 answer@pearson.resid <- as.matrix(Presid) 105 } 106 107 slot(answer, "misc") <- stuff@misc # Replace 108 109 110 answer@misc$signif.stars <- signif.stars # 20140728 111 answer@misc$nopredictors <- nopredictors # 20150831 112 113 114 if (is.numeric(stuff@dispersion)) 115 slot(answer, "dispersion") <- stuff@dispersion 116 117 118 119 120 121 122 try.this <- findFirstMethod("summaryvglmS4VGAM", object@family@vfamily) 123 if (length(try.this)) { 124 new.postslot <- 125 summaryvglmS4VGAM(object = object, 126 VGAMff = new(try.this), 127 ...) 128 answer@post <- new.postslot 129 } else { 130 } 131 132 133 134 control <- object@control 135 if (missing.HDEtest && 136 length(temp <- object@control$summary.HDEtest)) { 137 HDEtest <- temp 138 } 139 140 if (HDEtest) { 141 answer@post$hdeff <- hdeff(object, derivative = 1, se.arg = TRUE) 142 answer@post$hde.NA <- hde.NA 143 answer@post$threshold.hde <- threshold.hde 144 } 145 146 147 148 answer 149} # summary.vglm 150 151 152 153 154 155 156 157 158setMethod("summaryvglmS4VGAM", signature(VGAMff = "cumulative"), 159 function(object, 160 VGAMff, 161 ...) { 162 object@post <- 163 callNextMethod(VGAMff = VGAMff, 164 object = object, 165 ...) 166 object@post$reverse <- object@misc$reverse 167 168 169 cfit <- coef(object, matrix = TRUE) 170 M <- ncol(cfit) 171 if (rownames(cfit)[1] == "(Intercept)") 172 object@post$expcoeffs <- exp(coef(object)[-(1:M)]) 173 174 175 object@post 176}) 177 178 179 180setMethod("showsummaryvglmS4VGAM", signature(VGAMff = "cumulative"), 181 function(object, 182 VGAMff, 183 ...) { 184 185 if (length(object@post$expcoeffs)) { 186 cat("\nExponentiated coefficients:\n") 187 print(object@post$expcoeffs) 188 } 189 if (FALSE) { 190 if (object@post$reverse) 191 cat("Reversed\n\n") else 192 cat("Not reversed\n\n") 193 } 194}) 195 196 197 198 199 200 201setMethod("summaryvglmS4VGAM", signature(VGAMff = "multinomial"), 202 function(object, 203 VGAMff, 204 ...) { 205 object@post <- 206 callNextMethod(VGAMff = VGAMff, 207 object = object, 208 ...) 209 object@post$refLevel <- object@misc$refLevel 210 object@post 211}) 212 213 214 215setMethod("showsummaryvglmS4VGAM", signature(VGAMff = "multinomial"), 216 function(object, 217 VGAMff, 218 ...) { 219 cat("\nReference group is level ", object@post$refLevel, 220 " of the response\n") 221 callNextMethod(VGAMff = VGAMff, 222 object = object, 223 ...) 224}) 225 226 227 228setMethod("summaryvglmS4VGAM", signature(VGAMff = "VGAMcategorical"), 229 function(object, 230 VGAMff, 231 ...) { 232 object@post 233}) 234 235 236setMethod("showsummaryvglmS4VGAM", signature(VGAMff = "VGAMcategorical"), 237 function(object, 238 VGAMff, 239 ...) { 240}) 241 242 243 244 245 246 247 248 249 250setMethod("logLik", "summary.vglm", function(object, ...) 251 logLik.vlm(object, ...)) 252 253 254 255show.summary.vglm <- 256 function(x, 257 digits = max(3L, getOption("digits") - 3L), # Same as glm() 258 quote = TRUE, 259 prefix = "", 260 presid = length(x@pearson.resid) > 0, # FALSE, # TRUE, 261 HDEtest = TRUE, 262 hde.NA = TRUE, 263 threshold.hde = 0.001, 264 signif.stars = NULL, # Use this if logical; 20140728 265 nopredictors = NULL, # Use this if logical; 20150831 266 top.half.only = FALSE, # Added 20160803 267 ... # Added 20151214 268 ) { 269 270 271 272 M <- x@misc$M 273 coef3 <- x@coef3 # icients 274 correl <- x@correlation 275 276 digits <- if (is.null(digits)) options()$digits - 2 else digits 277 278 cat("\nCall:\n", paste(deparse(x@call), sep = "\n", collapse = "\n"), 279 "\n", sep = "") 280 281 282 283 284 if (HDEtest) 285 if (is.logical(x@post$hde.NA) && x@post$hde.NA) { 286 if (length(hado <- x@post$hdeff)) { 287 HDE <- is.Numeric(hado[, "deriv1"]) & # Could be all NAs 288 hado[, "deriv1"] < 0 289 if (any(HDE) && ncol(coef3) == 4) { 290 HDE <- HDE & (x@post$threshold.hde < coef3[, 4]) 291 coef3[HDE, 3:4] <- NA # 3:4 means WaldStat and p-value 292 } 293 } 294 } # is.logical(x@post$hde.NA) && x@post$hde.NA 295 296 297 298 299 Presid <- x@pearson.resid 300 rdf <- x@df[2] 301 if (presid && 302 length(Presid) && 303 all(!is.na(Presid)) && 304 is.finite(rdf)) { 305 306 if (rdf/M > 5) { 307 rq <- apply(as.matrix(Presid), 2, quantile) # 5 x M 308 dimnames(rq) <- list(c("Min", "1Q", "Median", "3Q", "Max"), 309 x@misc$predictors.names) 310 cat("\nPearson residuals:\n") 311 print(t(rq), digits = digits) 312 } else 313 if (rdf > 0) { 314 cat("\nPearson residuals:\n") 315 print(Presid, digits = digits) 316 } 317 } 318 319 320 use.signif.stars <- if (is.logical(signif.stars)) 321 signif.stars else x@misc$signif.stars # 20140728 322 if (!is.logical(use.signif.stars)) 323 use.signif.stars <- getOption("show.signif.stars") 324 325 326 use.nopredictors <- if (is.logical(nopredictors)) 327 nopredictors else x@misc$nopredictors # 20140728 328 if (!is.logical(use.nopredictors)) { 329 warning("cannot determine 'nopredictors'; choosing FALSE") 330 use.nopredictors <- FALSE 331 } 332 333 334 Situation <- -1 335 how.many <- c(length(x@coef4lrt0), 336 length(x@coef4score0), 337 length(x@coef4wald0)) 338 339 if (length(x@coef4lrt0)) { # && wald0.arg 340 cat(if (top.half.only) "\nParametric coefficients:" else 341 "\nLikelihood ratio test coefficients:", "\n") 342 printCoefmat(x@coef4lrt0, digits = digits, 343 signif.stars = use.signif.stars, 344 na.print = "NA", 345 P.values = TRUE, has.Pvalue = TRUE, 346 signif.legend = sum(how.many[-1]) == 0) # Last one 347 Situation <- 3 348 } 349 350 351 352 if (length(x@coef4score0)) { # && wald0.arg 353 cat(if (top.half.only) "\nParametric coefficients:" else 354 "\nRao score test coefficients:", "\n") 355 printCoefmat(x@coef4score0, digits = digits, 356 signif.stars = use.signif.stars, 357 na.print = "NA", 358 signif.legend = sum(how.many[3]) == 0) # Last one 359 Situation <- 4 360 } 361 362 363 364 if (length(x@coef4wald0)) { # && wald0.arg 365 cat(if (top.half.only) "\nParametric coefficients:" else 366 "\nWald (modified by IRLS iterations) coefficients:", "\n") 367 printCoefmat(x@coef4wald0, digits = digits, 368 signif.stars = use.signif.stars, 369 na.print = "NA") 370 Situation <- 1 371 } else 372 if (length(coef3) && Situation < 0) { 373 cat(if (top.half.only) "\nParametric coefficients:" else 374 "\nCoefficients:", "\n") 375 printCoefmat(coef3, digits = digits, 376 signif.stars = use.signif.stars, 377 na.print = "NA") 378 Situation <- 2 379 } # length(coef3) 380 381 382 383 384 385 386 387 if (top.half.only) 388 return(invisible(NULL)) 389 390 391 392 if (M >= 5) 393 cat("\nNumber of linear predictors: ", M, "\n") 394 395 if (!is.null(x@misc$predictors.names) && !use.nopredictors) { 396 if (M == 1) { 397 cat("\nName of linear predictor:", 398 paste(x@misc$predictors.names, collapse = ", "), "\n") 399 } else 400 if (M <= 12) { 401 LLL <- length(x@misc$predictors.names) 402 cat("\nNames of linear predictors:", 403 if (LLL == 1) 404 x@misc$predictors.names else 405 c(paste0(x@misc$predictors.names[-LLL], sep = ","), 406 x@misc$predictors.names[LLL]), fill = TRUE) 407 } 408 } 409 410 prose <- "" 411 if (length(x@dispersion)) { 412 if (is.logical(x@misc$estimated.dispersion) && 413 x@misc$estimated.dispersion) { 414 prose <- "(Estimated) " 415 } else { 416 if (is.numeric(x@misc$default.dispersion) && 417 x@dispersion == x@misc$default.dispersion) 418 prose <- "(Default) " 419 420 if (is.numeric(x@misc$default.dispersion) && 421 x@dispersion != x@misc$default.dispersion) 422 prose <- "(Pre-specified) " 423 } 424 if (any(x@dispersion != 1)) 425 cat(paste("\n", prose, "Dispersion Parameter for ", 426 x@family@vfamily[1], 427 " family: ", yformat(x@dispersion, digits), "\n", 428 sep = "")) 429 } 430 431 432 if (length(deviance(x))) { 433 cat("\nResidual deviance:", yformat(deviance(x), digits)) 434 if (is.finite(rdf)) 435 cat(" on", round(rdf, digits), "degrees of freedom\n") else 436 cat("\n") 437 } 438 439 440 if (length(vll <- logLik.vlm(x))) { 441 cat("\nLog-likelihood:", yformat(vll, digits)) 442 if (is.finite(rdf)) 443 cat(" on", round(rdf, digits), "degrees of freedom\n") else 444 cat("\n") 445 } 446 447 448 if (length(x@criterion)) { 449 ncrit <- names(x@criterion) 450 for (ii in ncrit) 451 if (ii != "loglikelihood" && ii != "deviance") 452 cat(paste(ii, ":", sep = ""), 453 yformat(x@criterion[[ii]], digits), "\n") 454 } 455 456 457 cat("\nNumber of Fisher scoring iterations:", 458 format(trunc(x@iter)), "\n\n") 459 460 461 if (!is.null(correl)) { 462 ncol.X.vlm <- dim(correl)[2] 463 if (ncol.X.vlm > 1) { 464 cat("Correlation of Coefficients:\n\n") 465 ll <- lower.tri(correl) 466 correl[ll] <- format(round(correl[ll], digits)) 467 correl[!ll] <- "" 468 print(correl[-1, -ncol.X.vlm, drop = FALSE], quote = FALSE, 469 digits = digits) 470 cat("\n") 471 } 472 } 473 474 475 476 477 478 if (HDEtest) 479 if (Situation == 2 && 480 length(hado <- x@post$hdeff)) { 481 if (is.Numeric(hado[, "deriv1"]) & # Could be all NAs 482 all(hado[, "deriv1"] > 0)) 483 cat("No Hauck-Donner effect found in any of the estimates\n\n") 484 if (is.Numeric(hado[, "deriv1"]) & # Could be all NAs 485 any(hado[, "deriv1"] < 0)) { 486 cat("Warning: Hauck-Donner effect detected in the", 487 "following estimate(s):\n") 488 cat(paste("'", rownames(hado)[hado[, "deriv1"] < 0], 489 "'", collapse = ", ", sep = "")) 490 cat("\n\n") 491 } 492 } # Situation == 2 && length(hado) 493 494 495 496 497 try.this <- findFirstMethod("showsummaryvglmS4VGAM", x@family@vfamily) 498 if (length(try.this)) { 499 showsummaryvglmS4VGAM(object = x, 500 VGAMff = new(try.this), 501 ...) 502 } else { 503 } 504 505 506 507 invisible(NULL) 508} # show.summary.vglm 509 510 511 512 513setMethod("summary", "vglm", 514 function(object, ...) 515 summaryvglm(object, ...)) 516 517 518 519 520 521 522setMethod("show", "summary.vglm", 523 function(object) 524 show.summary.vglm(object)) 525 526 527 528 529 530 531 532 533 534 535if (FALSE) 536show.summary.binom2.or <- 537 function(x, 538 digits = max(3L, getOption("digits") - 3L) # Same as glm() 539 ) { 540 541 if (length(x@post$oratio) == 1 && 542 is.numeric(x@post$oratio)) { 543 cat("\nOdds ratio: ", round(x@post$oratio, digits), "\n") 544 } 545} 546 547 548 549 550if (FALSE) 551setMethod("show", "summary.binom2.or", 552 function(object) 553 show.summary.vglm(object)) 554 555 556 557 558 559 560 561 562vcovdefault <- function(object, ...) { 563 if (is.null(object@vcov)) 564 stop("no default") 565 object@vcov 566} 567 568 569 570 571 572 573 574 575vcov.vlm <- function(object, ...) { 576 577 vcovvlm(object, ...) 578} # vcov.vlm 579 580 581 582 vcovvlm <- 583function(object, dispersion = NULL, untransform = FALSE, 584 complete = TRUE) { 585 586 587 588 so <- summaryvlm(object, correlation = FALSE, 589 dispersion = dispersion) 590 d <- if (any(slotNames(so) == "dispersion") && 591 is.Numeric(so@dispersion)) 592 so@dispersion else 1 593 answer <- d * so@cov.unscaled 594 595 if (is.logical(OKRC <- object@misc$RegCondOK) && !OKRC) 596 warning("MLE regularity conditions were violated ", 597 "at the final iteration of the fitted object") 598 599 if (!untransform) 600 return(answer) 601 602 603 604 605 606 new.way <- TRUE 607 608 609 610 if (!is.logical(object@misc$intercept.only)) 611 stop("cannot determine whether the object is", 612 "an intercept-only fit, i.e., 'y ~ 1' is the response") 613 if (!object@misc$intercept.only) 614 stop("object must be an intercept-only fit, i.e., ", 615 "y ~ 1 is the response") 616 617 if (!all(trivial.constraints(constraints(object)) == 1)) 618 stop("object must have trivial constraints") 619 620 M <- object@misc$M 621 622 623 624 625 tvector <- numeric(M) 626 etavector <- predict(object)[1, ] # Contains transformed parameters 627 LINK <- object@misc$link 628 EARG <- object@misc$earg # This could be a NULL 629 if (is.null(EARG)) 630 EARG <- list(theta = NULL) 631 if (!is.list(EARG)) 632 stop("the 'earg' component of 'object@misc' must be a list") 633 634 635 636 637 if (length(LINK) != M && 638 length(LINK) != 1) 639 stop("cannot obtain the link functions to untransform 'object'") 640 641 642 643 if (!is.character(LINK)) 644 stop("the 'link' component of 'object@misc' should ", 645 "be a character vector") 646 647 learg <- length(EARG) 648 llink <- length(LINK) 649 if (llink != learg) 650 stop("the 'earg' component of 'object@misc' should ", 651 "be a list of length ", learg) 652 653 654 level1 <- length(EARG) > 3 && 655 length(intersect(names(EARG), 656 c("theta", "inverse", "deriv", "short", "tag"))) > 3 657 if (level1) 658 EARG <- list(oneOnly = EARG) 659 660 661 662 learg <- length(EARG) 663 for (ii in 1:M) { 664 TTheta <- etavector[ii] # Transformed theta 665 666 use.earg <- 667 if (llink == 1) EARG[[1]] else EARG[[ii]] 668 function.name <- 669 if (llink == 1) LINK else LINK[ii] 670 671 672 if (new.way) { 673 use.earg[["inverse"]] <- TRUE # New 674 use.earg[["theta"]] <- TTheta # New 675 Theta <- do.call(function.name, use.earg) 676 677 678 use.earg[["inverse"]] <- TRUE # Reset this 679 680 681 use.earg[["deriv"]] <- 1 # New 682 use.earg[["theta"]] <- Theta # Renew this 683 tvector[ii] <- do.call(function.name, use.earg) 684 } else { 685 stop("link functions handled in the new way now") 686 687 } 688 } # of for (ii in 1:M) 689 690 tvector <- abs(tvector) 691 answer <- (cbind(tvector) %*% rbind(tvector)) * answer 692 if (length(dmn2 <- names(object@misc$link)) == M) 693 dimnames(answer) <- list(dmn2, dmn2) 694 answer 695} # vcovvlm 696 697 698 699 700 701 702 703setMethod("vcov", "vlm", 704 function(object, ...) 705 vcovvlm(object, ...)) 706 707 708setMethod("vcov", "vglm", 709 function(object, ...) 710 vcovvlm(object, ...)) 711 712 713 714 715 716 717 718 719yformat <- function(x, digits = options()$digits) { 720 format(ifelse(abs(x) < 0.001, signif(x, digits), round(x, digits))) 721} 722 723 724 725 726 727 728 729