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