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