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
17summaryvlm <-
18  function(object, correlation = FALSE, dispersion = NULL,
19           Colnames = c("Estimate", "Std. Error", "z value", "Pr(>|z|)"),
20           presid = TRUE,  # FALSE
21           lrt0.arg = FALSE,
22           score0.arg = FALSE,
23           wald0.arg = FALSE,
24           values0 = 0,
25           subset = NULL,
26           omit1s = TRUE) {
27
28
29
30  if (is.logical(object@misc$BFGS) && object@misc$BFGS)
31    warning("the estimated variance-covariance matrix is ",
32            "usually inaccurate because the working weight matrices ",
33            "are obtained by a crude BFGS quasi-Newton approximation")
34
35  M <- object@misc$M
36  n <- object@misc$n
37  nrow.X.vlm <- object@misc$nrow.X.vlm
38  ncol.X.vlm <- object@misc$ncol.X.vlm  # May be NULL for CQO objects
39
40  Coefs <- object@coefficients
41  cnames <- names(Coefs)
42
43  Presid <- if (presid) {
44    Presid <- residualsvlm(object, type = "pearson")
45    Presid
46  } else {
47    NULL
48  }
49
50  if (anyNA(Coefs)) {
51    warning("Some NAs in the coefficients---no summary is ",
52            " provided; returning 'object'")
53    return(object)
54  }
55  rdf <- object@df.residual
56
57  if (!length(dispersion)) {
58    if (is.numeric(object@misc$dispersion)) {
59      dispersion <- object@misc$dispersion
60      if (all(dispersion == 0))
61        stop("dispersion shouldn't be zero here!")
62    } else {
63      dispersion <- 1
64      object@misc$estimated.dispersion <- FALSE
65    }
66  } else if (dispersion == 0) {
67      dispersion <-
68        if (!length(object@ResSS)) {
69          stop("object@ResSS is empty")
70      } else {
71        object@ResSS / object@df.residual
72      }
73      object@misc$estimated.dispersion <- TRUE
74  } else {
75    if (is.numeric(object@misc$dispersion) &&
76        object@misc$dispersion != dispersion)
77      warning("overriding the value of object@misc$dispersion")
78    object@misc$estimated.dispersion <- FALSE
79  }
80  sigma <- sqrt(dispersion)  # Can be a vector
81
82  if (is.Numeric(ncol.X.vlm)) {
83    R <- object@R
84
85    if (ncol.X.vlm < max(dim(R)))
86      stop("'R' is rank deficient")
87
88
89
90
91
92    covun <- chol2inv(R)
93
94    dimnames(covun) <- list(cnames, cnames)
95  }
96  coef3 <- matrix(rep(Coefs, 4), ncol = 4)
97  dimnames(coef3) <- list(cnames, Colnames)
98  SEs <- sqrt(diag(covun))
99  if (length(sigma) == 1 && is.Numeric(ncol.X.vlm)) {
100    coef3[, 2] <- SEs %o% sigma  # Fails here when sigma is a vector
101    coef3[, 3] <- coef3[, 1] / coef3[, 2]
102    pvalue <- 2 * pnorm(-abs(coef3[, 3]))
103    coef3[, 4] <- pvalue
104
105    if (is.logical(object@misc$estimated.dispersion) &&
106        object@misc$estimated.dispersion)
107      coef3 <- coef3[, -4]  # Delete the pvalues column
108  } else {
109    coef3[, 1] <- coef3[, 2] <- coef3[, 3] <- coef3[, 4] <- NA
110    coef3 <- coef3[, -4]  # Delete the pvalues column
111  }
112
113
114
115
116
117
118  if (lrt0.arg) {
119    coef4lrt0 <- coef3[, -2, drop = FALSE]  # Omit SEs
120    lrt.list <- lrt.stat(object, all.out = TRUE,
121                    values0 = values0, subset = subset,
122                    trace = FALSE,
123                    omit1s = omit1s)  # Intercept-only model: NULL
124    lrt.list.values0 <- lrt.list$values0
125    SEs <- NA  # For correlation = TRUE
126
127
128    if (length(lrt.list)) {  # Usually omit intercepts:
129      coef4lrt0 <- coef4lrt0[names(lrt.list.values0), , drop = FALSE]
130
131
132      if (length(sigma) == 1 && is.Numeric(ncol.X.vlm)) {
133        coef4lrt0[, 'z value'] <- lrt.list$lrt.stat
134        coef4lrt0[, 'Pr(>|z|)'] <- lrt.list$pvalues
135
136        if (is.logical(object@misc$estimated.dispersion) &&
137            object@misc$estimated.dispersion)
138          coef4lrt0 <- coef4lrt0[, -3]  # Delete the pvalues column
139      } else {
140        coef4lrt0[, 1] <- coef4lrt0[, 2] <-
141        coef4lrt0[, 3] <- NA
142        coef4lrt0 <- coef4lrt0[, -3]  # Delete the pvalues column
143      }
144    } else {
145      coef4lrt0 <- new("matrix")  # Empty matrix, of length 0
146    }
147  } else {
148    coef4lrt0 <- new("matrix")  # Empty matrix, of length 0
149  }
150
151
152
153
154
155
156  if (score0.arg) {
157    coef4score0 <- coef3  # Overwrite some columns
158    score.list <- score.stat(object, all.out = TRUE,
159                             values0 = values0, subset = subset,
160                             trace = FALSE,
161                           omit1s = omit1s)  # Intercept-only model: NULL
162    SEs <- score.list$SE0
163    if (length(score.list)) {  # Usually omit intercepts:
164      coef4score0 <- coef4score0[names(SEs), , drop = FALSE]
165      if (length(sigma) == 1 && is.Numeric(ncol.X.vlm)) {
166        coef4score0[, 2] <- SEs %o% sigma  # Fails if sigma is a vector
167        coef4score0[, 3] <- score.list$score.stat
168        pvalue <- 2 * pnorm(-abs(coef4score0[, 3]))
169        coef4score0[, 4] <- pvalue
170
171        if (is.logical(object@misc$estimated.dispersion) &&
172            object@misc$estimated.dispersion)
173          coef4score0 <- coef4score0[, -4]  # Delete the pvalues column
174      } else {
175        coef4score0[, 1] <- coef4score0[, 2] <-
176        coef4score0[, 3] <- coef4score0[, 4] <- NA
177        coef4score0 <- coef4score0[, -4]  # Delete the pvalues column
178      }
179    } else {
180      coef4score0 <- new("matrix")  # Empty matrix, of length 0
181    }
182  } else {
183    coef4score0 <- new("matrix")  # Empty matrix, of length 0
184  }
185
186
187
188
189
190
191
192  if (wald0.arg) {
193    coef4wald0 <- coef3  # Overwrite some columns
194    SEs <- wald.stat(object, all.out = TRUE,
195                     values0 = values0, subset = subset,
196                     trace = FALSE,
197                     omit1s = omit1s)$SE0  # Intercept-only model: NULL
198    if (length(SEs)) {  # Usually omit intercepts:
199      coef4wald0 <- coef4wald0[names(SEs), , drop = FALSE]
200      if (length(sigma) == 1 && is.Numeric(ncol.X.vlm)) {
201        coef4wald0[, 2] <- SEs %o% sigma  # Fails if sigma is a vector
202        coef4wald0[, 3] <- coef4wald0[, 1] / coef4wald0[, 2]
203        pvalue <- 2 * pnorm(-abs(coef4wald0[, 3]))
204        coef4wald0[, 4] <- pvalue
205
206        if (is.logical(object@misc$estimated.dispersion) &&
207            object@misc$estimated.dispersion)
208          coef4wald0 <- coef4wald0[, -4]  # Delete the pvalues column
209      } else {
210        coef4wald0[, 1] <- coef4wald0[, 2] <-
211        coef4wald0[, 3] <- coef4wald0[, 4] <- NA
212        coef4wald0 <- coef4wald0[, -4]  # Delete the pvalues column
213      }
214    } else {
215      coef4wald0 <- new("matrix")  # Empty matrix, of length 0
216    }
217  } else {
218    coef4wald0 <- new("matrix")  # Empty matrix, of length 0
219  }
220
221
222
223  if (correlation) {
224    correl <- covun * outer(1 / SEs, 1 / SEs)
225
226    diag(correl) <- 1.0
227
228    dimnames(correl) <- list(cnames, cnames)
229  } else {
230    correl <- matrix(0, 0, 0)  # was NULL, but now a special matrix
231  }
232
233
234
235
236  answer <-
237  new("summary.vlm",
238      object,
239      coef3       = coef3,
240      coef4lrt0   = coef4lrt0,
241      coef4score0 = coef4score0,
242      coef4wald0  = coef4wald0,
243      correlation = correl,
244      df          = c(ncol.X.vlm, rdf),
245      sigma       = sigma)
246
247  if (is.Numeric(ncol.X.vlm))
248    answer@cov.unscaled <- covun
249  answer@dispersion <- dispersion  # Overwrite this
250
251  if (length(Presid))
252    answer@pearson.resid <- as.matrix(Presid)
253
254
255  answer
256}
257
258
259
260
261
262
263
264
265show.summary.vlm <- function(x, digits = NULL, quote = TRUE,
266                             prefix = "") {
267
268
269  M <- x@misc$M
270  coef3 <- x@coef3 # ficients
271  correl <- x@correlation
272
273  if (is.null(digits)) {
274    digits <- options()$digits
275  } else {
276    old.digits <- options(digits = digits)
277    on.exit(options(old.digits))
278  }
279
280  cat("\nCall:\n")
281  dput(x@call)
282
283  Presid <- x@pearson.resid
284  rdf <- x@df[2]
285  if (length(Presid) && all(!is.na(Presid))) {
286    if (rdf/M > 5) {
287      rq <-  apply(as.matrix(Presid), 2, quantile)  # 5 x M
288      dimnames(rq) <- list(c("Min", "1Q", "Median", "3Q", "Max"),
289                           x@misc$predictors.names)
290      cat("\nPearson residuals:\n")
291      print(t(rq), digits = digits)
292    } else
293    if (rdf > 0) {
294      cat("\nPearson residuals:\n")
295      print(Presid, digits = digits)
296    }
297  }
298
299  if (!all(is.na(coef3))) {
300    cat("\nCoefficients:\n")
301    print(coef3, digits = digits)
302  }
303
304  cat("\nNumber of responses: ", M, "\n")
305
306
307  if (length(x@misc$predictors.names))
308  if (M == 1) {
309    cat("\nName of response:",
310        paste(x@misc$predictors.names, collapse = ", "), "\n")
311  } else {
312    UUU <- paste(x@misc$predictors.names, collapse = ", ")
313    UUU <- x@misc$predictors.names
314    cat("\nNames of responses:\n")
315    cat(UUU, fill = TRUE, sep = ", ")
316  }
317
318
319  if (!is.null(x@ResSS))
320    cat("\nResidual Sum of Squares:", format(round(x@ResSS, digits)),
321        "on", round(rdf, digits), "degrees of freedom\n")
322
323
324  if (length(correl)) {
325    ncol.X.vlm <- dim(correl)[2]
326    if (ncol.X.vlm > 1) {
327      cat("\nCorrelation of Coefficients:\n")
328      ll <- lower.tri(correl)
329      correl[ll] <- format(round(correl[ll], digits))
330      correl[!ll] <- ""
331      print(correl[-1, -ncol.X.vlm, drop = FALSE],
332            quote = FALSE, digits = digits)
333    }
334  }
335
336  invisible(NULL)
337}
338
339
340setMethod("summary", "vlm",
341          function(object, ...)
342          summaryvlm(object, ...))
343
344
345
346
347setMethod("show", "summary.vlm",
348          function(object)
349          show.summary.vlm(object))
350
351
352
353