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 18vcov.pvgam <- function(object, ...) { 19 vcovpvgam(object, ...) 20} 21 22 23 24 vcovpvgam <- 25 function(object, 26 special = FALSE, 27 frequentist = FALSE, dispersion = NULL, unconditional = FALSE, 28 ...) { 29 30 31 32 33 34 if (!special) { 35 return(vcovvlm(object, ...)) 36 } 37 38 39 40 warning("vcovpvgam() is only 50% finished") 41 42 43 print("in vcovpvgam; hi 2a") 44 print("class(object)") 45 print( class(object) ) 46 47 48 49 50 M <- npred(object) 51 n <- nobs(object, type = "lm") 52 wz <- weights(object, type = "working") 53 X.vlm.save <- model.matrix(object, type = "vlm") 54 U <- vchol(wz, M = M, n = n) 55 X.vlm <- mux111(U, X.vlm.save, M = M) 56 X.vlm.aug <- rbind(X.vlm, 57 model.matrix(object, type = "penalty")) 58 59 60 qr1 <- qr(X.vlm.aug) 61 qr2 <- qr(X.vlm) 62 poststuff <- 63 mgcv::magic.post.proc(X.vlm.aug, 64 object = object@ospsslot$magicfit, w = NULL) 65 magicfit <- object@ospsslot$magicfit 66 rV <- magicfit$rV 67 Vb <- poststuff$Vb 68 Ve <- poststuff$Ve 69 hhat <- poststuff$hat 70 eedf <- poststuff$edf 71 72 73 scale.param <- 1 # Assumed 74 75 76 vc <- if (frequentist) { 77 78 mat1 <- solve(crossprod(qr.R(qr1))) 79 scale.param * 80 (mat1 %*% crossprod(qr.R(qr2)) %*% mat1) 81 } else { 82 Vc <- NULL # Corrected ML or REML is not available. 83 84 Vp <- scale.param * tcrossprod(solve(qr.R(qr1))) 85 86 Vp2 <- rV %*% t(rV) # * sig2 # For checking 87 88 print("max(abs(Vp - Vp2)); should be 0") 89 print( max(abs(Vp - Vp2)) ) 90 91 92 if (FALSE) { 93 He <- SemiParFit$fit$hessian 94 He.eig <- eigen(He, symmetric=TRUE) 95 Vb <- He.eig$vectors %*% 96 tcrossprod(diag(1/He.eig$values), 97 He.eig$vectors) # this could be taken from magic as well 98 Vb <- (Vb + t(Vb) ) / 2 99 100 101 HeSh <- He - SemiParFit$fit$S.h 102 F <- Vb%*%HeSh # diag(SemiParFit$magpp$edf) 103 104 105 HeSh <- He 106 Ve <- Vb 107 F <- F1 <- diag(rep(1,dim(Vb)[1])) 108 R <- SemiParFit$bs.mgfit$R 109 } 110 111 112 113 if (unconditional && !is.null(Vc)) 114 Vc else Vp 115 } # Bayesian 116 if (is.null(dispersion)) { 117 sig2 <- 1 # zz 118 vc <- summary(object)@dispersion * vc / sig2 119 } else { 120 sig2 <- summary(object)@dispersion # 1 # zz 121 vc <- dispersion * vc / sig2 122 } 123 124 125 126 print("head(sort(diag(vc)))") 127 print( head(sort(diag(vc))) ) 128 print("head(sort(diag(Ve)))") 129 print( head(sort(diag(Ve))) ) 130 131 132 print("tail(sort(diag(vc)))") 133 print( tail(sort(diag(vc))) ) 134 print("tail(sort(diag(Ve)))") 135 print( tail(sort(diag(Ve))) ) 136 137 138 139 print("head(sort(diag(vc))) / head(sort(diag(Ve)))") 140 print( head(sort(diag(vc))) / head(sort(diag(Ve))) ) 141 142 143 144 145 print("max(abs(sort(diag(vc)) - sort(diag(Ve))))") 146 print( max(abs(sort(diag(vc)) - sort(diag(Ve)))) ) 147 148 149 150 151 152 vc 153} 154 155 156 157 158 159 160setMethod("vcov", "pvgam", 161 function(object, ...) 162 vcovpvgam(object, ...)) 163 164 165 166 167 168 169 170 171startstoppvgam <- 172 function(object, ...) { 173 174 which.X.sm.osps <- object@ospsslot$sm.osps.list$which.X.sm.osps 175 if (!length(which.X.sm.osps)) 176 stop("no 'sm.os()' or 'sm.ps()' term in 'object'") 177 all.ncol.Hk <- unlist(lapply(constraints(object, type = "term"), ncol)) 178 names.which.X.sm.osps <- names(which.X.sm.osps) 179 endf <- rep_len(NA_real_, sum(all.ncol.Hk[names.which.X.sm.osps])) 180 names(endf) <- vlabel(names.which.X.sm.osps, 181 all.ncol.Hk[names.which.X.sm.osps], 182 M = npred(object)) 183 stopstart <- NULL 184 185 186 iptr <- 1 187 iterm <- 1 188 for (ii in names(all.ncol.Hk)) { 189 if (length(which.X.sm.osps[[ii]])) { 190 temp3 <- -1 + iptr + all.ncol.Hk[ii] * length(which.X.sm.osps[[ii]]) 191 new.index <- iptr:temp3 # Includes all component functions wrt xk 192 iptr <- iptr + length(new.index) # temp3 193 mat.index <- matrix(new.index, ncol = all.ncol.Hk[ii], byrow = TRUE) 194 for (jay in 1:all.ncol.Hk[ii]) { 195 cf.index <- mat.index[, jay] 196 197 198 stopstart <- c(stopstart, list(cf.index)) 199 200 201 iterm <- iterm + 1 202 } # for 203 } else { 204 iptr <- iptr + all.ncol.Hk[ii] 205 } 206 } # ii 207 names(stopstart) <- names(endf) 208 stopstart 209} 210 211 212 213 214 215 216 217 218 219 220 221summarypvgam <- 222 function(object, dispersion = NULL, 223 digits = options()$digits-2, 224 presid = TRUE) { 225 226 stuff <- summaryvglm(object, dispersion = dispersion, 227 digits = digits, 228 presid = presid) 229 230 answer <- 231 new("summary.pvgam", 232 object, 233 call = stuff@call, 234 cov.unscaled = stuff@cov.unscaled, 235 correlation = stuff@correlation, 236 df = stuff@df, 237 sigma = stuff@sigma) 238 239 answer@misc$nopredictors <- stuff@misc$nopredictors 240 answer@ospsslot <- object@ospsslot 241 242 243 slot(answer, "coefficients") <- stuff@coefficients # Replace 244 245 coef3 <- stuff@coef3 246 aassign <- attr(model.matrix(object, type = "vlm"), "assign") 247 myterms <- names(object@ospsslot$sm.osps.list$which.X.sm.osps) 248 index.exclude <- NULL 249 for (ii in myterms) { 250 index.exclude <- c(index.exclude, unlist(aassign[[ii]])) 251 } 252 slot(answer, "coef3") <- coef3[-index.exclude, , drop = FALSE] 253 254 255 256 if (is.numeric(stuff@dispersion)) 257 slot(answer, "dispersion") <- stuff@dispersion 258 259 if (presid) { 260 Presid <- residuals(object, type = "pearson") 261 if (length(Presid)) 262 answer@pearson.resid <- as.matrix(Presid) 263 } 264 265 266 267 268 269 270 271 pinv <- function(V, M, rank.tol = 1e-6) { 272 273 274 275 276 D <- eigen(V, symmetric = TRUE) 277 M1 <- length(D$values[D$values > rank.tol * D$values[1]]) 278 if (M > M1) 279 M <- M1 # avoid problems with zero eigen-values 280 281 if (M+1 <= length(D$values)) 282 D$values[(M+1):length(D$values)] <- 1 283 D$values <- 1 / D$values 284 if (M+1 <= length(D$values)) 285 D$values[(M+1):length(D$values)] <- 0 286 res <- D$vectors %*% (D$values * t(D$vectors)) ##D$u%*%diag(D$d)%*%D$v 287 attr(res, "rank") <- M 288 res 289 } ## end of pinv 290 291 292 293 startstop <- startstoppvgam(object) 294 m <- length(startstop) 295 296 df <- edf1 <- edf <- s.pv <- chi.sq <- array(0, m) 297 names(chi.sq) <- names(startstop) 298 p.type <- 5 # Frequentist 299 est.disp <- if (is.logical(object@misc$estimated.dispersion)) 300 object@misc$estimated.dispersion else FALSE 301 302 pvgam.residual.df <- df.residual_pvgam(object) 303 304 305 306 for (i in 1:m) { 307 308 p <- coef(as(object, "pvgam"))[(startstop[[i]])] # params for smooth 309 310 endf <- endfpvgam(object, diag.all = TRUE) # This is ENDF+1 actually 311 312 313 edf1[i] <- edf[i] <- sum(endf[(startstop[[i]])]) 314 if (FALSE && !is.null(object$edf1)) 315 edf1[i] <- sum(object$edf1[(startstop[[i]])]) 316 317 V <- if (p.type == 5) { 318 Ve <- vcov(object, special = FALSE) 319 Ve[(startstop[[i]]), (startstop[[i]]), drop = FALSE] 320 } else { 321 Vp <- vcov(object, special = TRUE, frequentist = FALSE) 322 Vp[(startstop[[i]]), (startstop[[i]]), drop = FALSE] 323 } 324 325 if (p.type == 5) { 326 M1 <- length(startstop[[i]]) # zz 327 M <- min(M1, 328 ceiling(2*sum(endf[(startstop[[i]])])) 329 ) 330 V <- pinv(V, M) # , rank.tol = 1e-5 331 chi.sq[i] <- t(p) %*% V %*% p 332 df[i] <- attr(V, "rank") 333 } 334 335 336 337 if (p.type == 5) { 338 s.pv[i] <- if (est.disp) { 339 pf(chi.sq[i] / df[i], df1 = df[i], df2 = pvgam.residual.df, 340 lower.tail = FALSE) 341 } else { 342 pchisq(chi.sq[i], df = df[i], lower.tail = FALSE) 343 } 344 if (df[i] < 0.1) 345 s.pv[i] <- NA 346 } 347 348 349 if (est.disp) { 350 if (p.type == 5) { 351 s.table <- cbind(edf, df, chi.sq / df, s.pv) 352 dimnames(s.table) <- list(names(chi.sq), 353 c("edf", "Est.rank", "F", "p-value")) 354 } else { 355 s.table <- cbind(edf, df, chi.sq/df, s.pv) 356 dimnames(s.table) <- list(names(chi.sq), 357 c("edf", "Ref.df", "F", "p-value")) 358 } 359 } else { 360 if (p.type == 5) { 361 # This case is commonly executed 362 s.table <- cbind(edf, df, chi.sq, s.pv) 363 dimnames(s.table) <- list(names(chi.sq), 364 c("edf", "Est.rank", "Chi.sq", "p-value")) 365 } else { 366 s.table <- cbind(edf, df, chi.sq, s.pv) 367 dimnames(s.table) <- list(names(chi.sq), 368 c("edf", "Ref.df", "Chi.sq", "p-value")) 369 } 370 } # else 371 } # for (i) 372 answer@post$s.table <- s.table 373 374 375 376 377 378 aod <- data.frame(message = 'this does not work yet') 379 slot(answer, "anova") <- aod 380 381 answer 382} # summarypvgam() 383 384 385 386 387 388 389 390show.summary.pvgam <- 391 function(x, quote = TRUE, prefix = "", 392 digits = options()$digits-2, 393 signif.stars = getOption("show.signif.stars")) { 394 395 396 397 398 show.summary.vglm(x, quote = quote, prefix = prefix, 399 digits = digits, top.half.only = TRUE) 400 401 402 403 startstop <- startstoppvgam(x) 404 m <- length(startstop) 405 s.table <- x@post$s.table 406 if (0 < m && length(s.table)) { 407 cat("\nApproximate significance of smooth terms:\n") 408 printCoefmat(s.table, digits = digits, 409 signif.stars = signif.stars, has.Pvalue = TRUE, 410 na.print = "NA", cs.ind = 1) 411 } 412 413 414 415 416 417 418 M <- x@misc$M 419 420 421 Presid <- x@pearson.resid 422 rdf <- x@df[2] 423 424 425 cat("\nNumber of linear/additive predictors: ", M, "\n") 426 427 if (!is.null(x@misc$predictors.names)) 428 if (M == 1) 429 cat("\nName of linear/additive predictor:", 430 paste(x@misc$predictors.names, collapse = ", "), "\n") else 431 if (M <= 5) 432 cat("\nNames of linear/additive predictors:", 433 paste(x@misc$predictors.names, collapse = ", "), "\n") 434 435 prose <- "" 436 if (length(x@dispersion)) { 437 if (is.logical(x@misc$estimated.dispersion) && 438 x@misc$estimated.dispersion) { 439 prose <- "(Estimated) " 440 } else { 441 if (is.numeric(x@misc$default.dispersion) && 442 x@dispersion == x@misc$default.dispersion) 443 prose <- "(Default) " 444 445 if (is.numeric(x@misc$default.dispersion) && 446 x@dispersion != x@misc$default.dispersion) 447 prose <- "(Pre-specified) " 448 } 449 cat(paste("\n", prose, "Dispersion Parameter for ", 450 x@family@vfamily[1], 451 " family: ", 452 format(round(x@dispersion, digits)), "\n", sep = "")) 453 } 454 455 if (length(deviance(x))) 456 cat("\nResidual deviance: ", format(round(deviance(x), digits)), 457 "on", format(round(rdf, 3)), "degrees of freedom\n") 458 459 if (length(logLik.vlm(x))) 460 cat("\nLog-likelihood:", format(round(logLik.vlm(x), digits)), 461 "on", format(round(rdf, 3)), "degrees of freedom\n") 462 463 if (length(x@criterion)) { 464 ncrit <- names(x@criterion) 465 for (ii in ncrit) 466 if (ii != "loglikelihood" && ii != "deviance") 467 cat(paste(ii, ":", sep = ""), format(x@criterion[[ii]]), "\n") 468 } 469 470 471 472 473 if (is.Numeric(x@ospsslot$iter.outer)) { 474 cat("\nNumber of outer iterations: ", x@ospsslot$iter.outer, "\n") 475 cat("\nNumber of IRLS iterations at final outer iteration: ", x@iter, 476 "\n") 477 } else { 478 cat("\nNumber of IRLS iterations: ", x@iter, "\n") 479 } 480 481 482 if (FALSE && length(x@anova)) { 483 show.vanova(x@anova, digits = digits) # ".vanova" for Splus6 484 } 485 486 invisible(NULL) 487} # show.summary.pvgam() 488 489 490 491 492 493 494 495 496setMethod("summary", "pvgam", 497 function(object, ...) 498 summarypvgam(object, ...)) 499 500 501 502setMethod("show", "summary.pvgam", 503 function(object) 504 show.summary.pvgam(object)) 505 506 507 508 509 510 511psintpvgam <- function(object, ...) { 512 object@ospsslot$sm.osps.list$ps.int 513} 514 515 516if (!isGeneric("psint")) 517 setGeneric("psint", function(object, ...) 518 standardGeneric("psint"), 519 package = "VGAM") 520 521 522setMethod("psint", "pvgam", 523 function(object, ...) 524 psintpvgam(object, ...)) 525 526 527 528 529 530 531