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
14endfpvgam <- function(object,
15                      nonlinear.edf = TRUE,
16                      diag.all = FALSE,
17                      return.endf = TRUE, ...) {
18
19
20
21
22
23  M <- npred(object)
24  n <- nobs(object, type = "lm")
25  wz <- weights(object, type = "working")
26  X.vlm.save <- model.matrix(object, type = "vlm")
27  U <- vchol(wz, M = M, n = n)
28  X.vlm <- mux111(U, X.vlm.save, M = M)
29  X.vlm.aug <- rbind(X.vlm,
30                 model.matrix(object, type = "penalty"))
31
32
33
34  if (!object@ospsslot$magicfit$gcv.info$fully.converged)
35    warning("fitted object has a GCV criterion that has not ",
36            "fully converged")
37
38
39
40
41
42  poststuff <-
43    mgcv::magic.post.proc(X.vlm.aug,
44                          object = object@ospsslot$magicfit, w = NULL)
45
46
47  if (!return.endf)
48    return(poststuff)
49
50
51  which.X.sm.osps <- object@ospsslot$sm.osps.list$which.X.sm.osps
52  all.ncol.Hk <- unlist(lapply(constraints(object, type = "term"), ncol))
53  names.which.X.sm.osps <- names(which.X.sm.osps)
54  endf <- rep_len(NA_real_, sum(all.ncol.Hk[names.which.X.sm.osps]))
55  names(endf) <- vlabel(names.which.X.sm.osps,
56                        all.ncol.Hk[names.which.X.sm.osps],
57                        M = npred(object))
58  use.index <- NULL
59
60
61  endf.all0 <-  diag(solve(crossprod(X.vlm.aug), crossprod(X.vlm)))
62
63
64  if (FALSE) {
65  qr1 <- qr(X.vlm.aug)
66  qr2 <- qr(X.vlm)
67  endf.all <-  diag(solve(crossprod(qr.R(qr1)), crossprod(qr.R(qr2))))
68  }
69
70
71
72  endf.all <- endf.all0
73
74
75
76
77  if (diag.all)
78    return(endf.all)
79
80
81
82
83  startstop <- startstoppvgam(object)
84  for (iterm in 1:length(startstop)) {
85    endf[iterm] <- sum(endf.all[(startstop[[iterm]])])
86  }
87  endf[endf < 1] <- 1  # Cannot be smoother than linear
88
89
90  if (nonlinear.edf) endf - 1 else endf
91}  # endfpvgam()
92
93
94
95
96
97show.pvgam <- function(object) {
98
99  digits <- 3
100
101
102  if (!is.null(cl <- object@call)) {
103    cat("\nCall:\n", paste(deparse(cl), sep = "\n", collapse = "\n"),
104        "\n\n", sep = "")
105  }
106
107
108
109  magicfit <- object@ospsslot$magicfit
110
111
112
113
114  if (FALSE) {
115  XX <- model.matrix(object, type = "vlm")
116  poststuff <-
117    mgcv::magic.post.proc(XX,
118                          object = object@ospsslot$magicfit, w = NULL)
119  }
120
121
122
123  if (FALSE) {
124    edf <- rep_len(NA_real_, n.smooth)
125    cat("\nEstimated degrees of freedom:\n")
126    for (i in 1:n.smooth)
127      edf[i] <- sum(x$edf[x$smooth[[i]]$first.para:x$smooth[[i]]$last.para])
128    edf.str <- format(round(edf, digits = 4), digits = 3, scientific = FALSE)
129    for (i in 1:n.smooth) {
130      cat(edf.str[i], " ", sep = "")
131      if (i%%7 == 0)
132        cat("\n")
133    }
134    cat(" total =", round(sum(poststuff$edf), digits = 2), "\n")
135  }
136
137
138  endf <- endfpvgam(object)
139  cat("\nEstimated nonlinear degrees of freedom:\n")  # based on endfpvgam()
140  print(round(endf, digits = digits + 2), digits = digits, scientific = FALSE)
141
142  if (length(endf) > 1)
143  cat("Total:",
144      format(sum(endf), digits = digits), "\n")
145
146  object@post$endf <- endf  # Good to save this on the object
147
148
149  if (FALSE)
150  cat("\nEstimated degrees of freedom based on poststuff:",
151      format(poststuff$edf, digits = digits),
152      "\nTotal:",
153      format(round(sum(poststuff$edf), digits = digits)), "\n")
154
155
156  cat("\nUBRE score:", format(magicfit$score, digits = digits + 1), "\n\n")
157
158
159  if (length(deviance(object)))
160    cat("Residual deviance:", format(deviance(object)), "\n")
161
162
163  llx <- logLik.vlm(object = object)
164  if (length(llx))
165    cat("Log-likelihood:", format(llx), "\n")
166
167
168
169
170  invisible(object)
171}
172
173
174
175setMethod("show", "pvgam", function(object) show.pvgam(object))
176
177
178
179
180
181
182
183if (!isGeneric("endf"))
184    setGeneric("endf", function(object, ...)
185    standardGeneric("endf"))
186
187
188setMethod("endf", "pvgam", function(object, ...)
189          endfpvgam(object, ...))
190
191setMethod("endf", "summary.pvgam", function(object, ...)
192          endfpvgam(object, ...))
193
194
195
196
197
198
199
200
201show.vglm <- function(object) {
202
203  if (!is.null(cl <- object@call)) {
204    cat("\nCall:\n", paste(deparse(cl), sep = "\n", collapse = "\n"),
205        "\n\n", sep = "")
206  }
207
208  coef <- object@coefficients
209  if (any(nas <- is.na(coef))) {
210    if (is.null(names(coef)))
211      names(coef) <- paste("b", seq_along(coef), sep = "")
212    cat("\nCoefficients: (", sum(nas),
213        " not defined because of singularities)\n", sep = "")
214  } else {
215    cat("\nCoefficients:\n")
216  }
217  print(coef)
218
219  rank <- object@rank
220  if (!length(rank))
221    rank <- sum(!nas)
222  nobs <- if (length(object@df.total)) object@df.total else
223          length(object@residuals)
224  rdf <- object@df.residual
225  if (!length(rdf))
226    rdf <- nobs - rank
227  cat("\nDegrees of Freedom:", nobs, "Total;", rdf, "Residual\n")
228
229  if (length(deviance(object)))
230    cat("Residual deviance:", format(deviance(object)), "\n")
231  llx <- logLik.vlm(object = object)
232
233  if (length(llx))
234    cat("Log-likelihood:", format(llx), "\n")
235
236  if (length(object@criterion)) {
237    ncrit <- names(object@criterion)
238    for (ii in ncrit)
239      if (ii != "loglikelihood" &&
240          ii != "deviance")
241          cat(paste(ii, ":", sep = ""),
242              format(object@criterion[[ii]]), "\n")
243  }
244
245
246
247
248
249  try.this <- findFirstMethod("showvglmS4VGAM", object@family@vfamily)
250  if (length(try.this)) {
251    showvglmS4VGAM(object = object,
252                   VGAMff = new(try.this))
253  } else {
254  }
255
256
257
258
259  invisible(object)
260}
261
262
263
264
265
266
267show.vgam <- function(object) {
268
269  digits <- 2
270
271
272  if (!is.null(cl <- object@call)) {
273    cat("\nCall:\n", paste(deparse(cl), sep = "\n", collapse = "\n"),
274        "\n\n", sep = "")
275  }
276
277  coef <- object@coefficients
278  nas <- is.na(coef)
279
280  rank <- object@rank
281  if (is.null(rank))
282      rank <- sum(!nas)
283  nobs <- if (length(object@df.total)) object@df.total else
284          length(object@residuals)
285  rdf <- object@df.residual
286  if (is.null(rdf))
287    rdf <- nobs - rank
288  cat("\nDegrees of Freedom:", nobs, "Total;",
289      format(round(rdf, digits = digits)), "Residual\n")
290
291  if (length(deviance(object)))
292    cat("Residual deviance:", format(deviance(object)), "\n")
293
294  llx <- logLik.vlm(object = object)
295
296  if (length(llx))
297    cat("Log-likelihood:", format(llx), "\n")
298
299  criterion <- attr(terms(object), "criterion")
300  if (!is.null(criterion) &&
301      criterion != "coefficients")
302    cat(paste(criterion, ":", sep = ""),
303        format(object[[criterion]]), "\n")
304
305
306
307
308  try.this <- findFirstMethod("showvgamS4VGAM", object@family@vfamily)
309  if (length(try.this)) {
310    showvgamS4VGAM(object = object,
311                   VGAMff = new(try.this))
312  } else {
313  }
314
315
316
317  invisible(object)
318}
319
320
321
322
323setMethod("show",  "vlm", function(object) show.vlm (object))
324setMethod("show", "vglm", function(object) show.vglm(object))
325setMethod("show", "vgam", function(object) show.vgam(object))
326
327
328
329
330
331
332
333
334
335