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