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
15attrassigndefault <- function(mmat, tt) {
16  if (!inherits(tt, "terms"))
17    stop("need terms object")
18  aa <- attr(mmat, "assign")
19  if (is.null(aa))
20    stop("argument is not really a model matrix")
21  ll <- attr(tt, "term.labels")
22  if (attr(tt, "intercept") > 0)
23    ll <- c("(Intercept)", ll)
24  aaa <- factor(aa, labels = ll)
25  split(order(aa), aaa)
26}
27
28
29attrassignlm <- function(object, ...)
30  attrassigndefault(model.matrix(object), object@terms)
31
32
33
34
35
36 vlabel <-
37  function(xn, ncolHlist, M, separator = ":", colon = FALSE) {
38
39  if (length(xn) != length(ncolHlist))
40    stop("length of first two arguments not equal")
41
42  n1 <- rep(xn, ncolHlist)
43  if (M == 1)
44    return(n1)
45  n2 <- as.list(ncolHlist)
46  n2 <- lapply(n2, seq)
47  n2 <- unlist(n2)
48  n2 <- as.character(n2)
49  n2 <- paste(separator, n2, sep = "")
50  n3 <- rep(ncolHlist, ncolHlist)
51  if (!colon)
52    n2[n3 == 1] <- ""
53  n1n2 <- paste(n1, n2, sep = "")
54  n1n2
55}
56
57
58
59
60
61
62 vlm2lm.model.matrix <-
63  function(x.vlm, Hlist = NULL,
64           which.linpred = 1,
65           M = NULL) {
66
67
68
69
70
71
72
73  if (is.numeric(M)) {
74    if (M != nrow(Hlist[[1]]))
75      stop("argument 'M' does not match argument 'Hlist'")
76  } else {
77    M <- nrow(Hlist[[1]])
78  }
79
80
81  Hmatrices <- matrix(c(unlist(Hlist)), nrow = M)
82  if (ncol(Hmatrices) != ncol(x.vlm))
83    stop("ncol(Hmatrices) != ncol(x.vlm)")
84
85
86  n.lm <- nrow(x.vlm) / M
87  if (round(n.lm) != n.lm)
88    stop("'n.lm' does not seem to be an integer")
89  linpred.index <- which.linpred
90
91  if (FALSE) {
92    vecTF <- Hmatrices[linpred.index, ] != 0
93    X.lm.jay <- x.vlm[(0:(n.lm - 1)) * M + linpred.index, vecTF,
94                      drop = FALSE]
95  }
96  vecTF2 <- colSums(Hmatrices[linpred.index, , drop = FALSE]) != 0
97  vecTF1 <- rep_len(FALSE, M)
98  vecTF1[linpred.index] <- TRUE
99  X.lm.jay <- x.vlm[vecTF1, vecTF2, drop = FALSE]  # Recycling
100
101  X.lm.jay
102}  # vlm2lm.model.matrix
103
104
105
106
107
108
109 lm2vlm.model.matrix <-
110  function(x, Hlist = NULL, assign.attributes = TRUE,
111           M = NULL, xij = NULL, Xm2 = NULL) {
112
113
114
115
116  if (length(Hlist) != ncol(x))
117    stop("length(Hlist) != ncol(x)")
118
119  if (length(xij)) {
120    if (inherits(xij, "formula"))
121      xij <- list(xij)
122    if (!is.list(xij))
123      stop("'xij' is not a list of formulae")
124  }
125
126  if (!is.numeric(M))
127    M <- nrow(Hlist[[1]])
128
129  nrow.X.lm <- nrow(x)
130  if (all(trivial.constraints(Hlist) == 1)) {
131    X.vlm <- if (M > 1) kronecker(x, diag(M)) else x
132    ncolHlist <- rep(M, ncol(x))
133  } else {
134    allB <- matrix(unlist(Hlist), nrow = M)
135    ncolHlist <- unlist(lapply(Hlist, ncol))
136    Rsum <- sum(ncolHlist)
137
138    X1 <- rep(c(t(x)), rep(ncolHlist, nrow.X.lm))
139    dim(X1) <- c(Rsum, nrow.X.lm)
140    X.vlm <- kronecker(t(X1), matrix(1, M, 1)) *
141             kronecker(matrix(1, nrow.X.lm, 1), allB)
142    rm(X1)
143  }
144
145  dn <- labels(x)
146  yn <- dn[[1]]
147  xn <- dn[[2]]
148  dimnames(X.vlm) <- list(vlabel(yn, rep(M, nrow.X.lm), M),
149                          vlabel(xn, ncolHlist, M))
150
151  if (assign.attributes) {
152    attr(X.vlm, "contrasts")   <- attr(x, "contrasts")
153    attr(X.vlm, "factors")     <- attr(x, "factors")
154    attr(X.vlm, "formula")     <- attr(x, "formula")
155    attr(X.vlm, "class")       <- attr(x, "class")
156    attr(X.vlm, "order")       <- attr(x, "order")
157    attr(X.vlm, "term.labels") <- attr(x, "term.labels")
158
159    orig.assign.lm <- attr(x, "orig.assign.lm")  # NULL if x = F
160    nasgn <- oasgn <- attr(x, "assign")
161    lowind <- 0
162    for (ii in seq_along(oasgn)) {
163      mylen <- length(oasgn[[ii]]) * ncolHlist[oasgn[[ii]][1]]
164      nasgn[[ii]] <- (lowind+1):(lowind+mylen)
165      lowind <- lowind + mylen
166    } # End of ii
167    if (lowind != ncol(X.vlm))
168      stop("something gone wrong")
169    attr(X.vlm, "assign") <- nasgn
170
171
172    fred <- unlist(lapply(nasgn,
173                          length))/unlist(lapply(oasgn, length))
174    vasgn <- vector("list", sum(fred))
175    kk <- 0
176    for (ii in seq_along(oasgn)) {
177      temp <- matrix(nasgn[[ii]], ncol = length(oasgn[[ii]]))
178      for (jloc in 1:nrow(temp)) {
179        kk <- kk + 1
180        vasgn[[kk]] <- temp[jloc, ]
181      }
182    }
183    names(vasgn) <- vlabel(names(oasgn), fred, M)
184    attr(X.vlm, "vassign") <- vasgn
185
186    attr(X.vlm, "constraints") <- Hlist
187
188
189    attr(X.vlm, "orig.assign.lm") <- orig.assign.lm
190  }  # End of if (assign.attributes)
191
192
193
194
195  if (!length(xij))
196    return(X.vlm)
197
198
199
200
201
202
203
204  at.x <- attr(x, "assign")
205  at.vlmx <- attr(X.vlm, "assign")
206  at.Xm2 <- attr(Xm2, "assign")
207
208  for (ii in seq_along(xij)) {
209    form.xij <- xij[[ii]]
210    if (length(form.xij) != 3)
211      stop("xij[[", ii, "]] is not a formula with a response")
212    tform.xij <- terms(form.xij)
213    aterm.form <- attr(tform.xij, "term.labels")
214    if (length(aterm.form) != M)
215      stop("xij[[", ii, "]] does not contain ", M, " terms")
216
217    name.term.y <- as.character(form.xij)[2]
218    cols.X.vlm <- at.vlmx[[name.term.y]]  # May be > 1 in length.
219
220    x.name.term.2 <- aterm.form[1]  # Choose the first one
221    One.such.term <- at.Xm2[[x.name.term.2]]
222    for (bbb in seq_along(One.such.term)) {
223      use.cols.Xm2 <- NULL
224      for (sss in 1:M) {
225        x.name.term.2 <- aterm.form[sss]
226        one.such.term <- at.Xm2[[x.name.term.2]]
227        use.cols.Xm2 <- c(use.cols.Xm2, one.such.term[bbb])
228      } # End of sss
229
230      allXk <- Xm2[, use.cols.Xm2, drop = FALSE]
231      cmat.no <- (at.x[[name.term.y]])[1]
232      cmat <- Hlist[[cmat.no]]
233      Rsum.k <- ncol(cmat)
234      tmp44 <- kronecker(matrix(1, nrow.X.lm, 1), t(cmat)) *
235       kronecker(allXk, matrix(1, ncol(cmat), 1))  # n*Rsum.k x M
236
237      tmp44 <- array(t(tmp44), c(M, Rsum.k, nrow.X.lm))
238      tmp44 <- aperm(tmp44, c(1, 3, 2))  # c(M, n, Rsum.k)
239      rep.index <- cols.X.vlm[((bbb-1)*Rsum.k+1):(bbb*Rsum.k)]
240      X.vlm[, rep.index] <- c(tmp44)
241    }  # End of bbb
242  }  # End of for (ii in seq_along(xij))
243
244  if (assign.attributes) {
245    attr(X.vlm, "vassign") <- vasgn
246    attr(X.vlm, "assign") <- nasgn
247    attr(X.vlm, "xij") <- xij
248  }
249  X.vlm
250}  # lm2vlm.model.matrix
251
252
253
254
255
256
257
258
259
260model.matrix.vlm <- function(object, ...)
261  model.matrixvlm(object, ...)
262
263
264
265
266
267 model.matrixvlm <-
268  function(object,
269           type = c("vlm", "lm", "lm2", "bothlmlm2"),
270           linpred.index = NULL,
271           ...) {
272
273
274
275
276
277  if (mode(type) != "character" && mode(type) != "name")
278    type <- as.character(substitute(type))
279  type <- match.arg(type, c("vlm", "lm", "lm2", "bothlmlm2"))[1]
280
281  linpred.index <- unique(sort(linpred.index))
282  LLLL <- length(linpred.index)
283  if (LLLL == 1 && type != "lm")
284    stop("Must set 'type = \"lm\"' when 'linpred.index' is ",
285         "assigned a single value")
286  if (LLLL  > 1 && type != "vlm")
287    stop("Must set 'type = \"vlm\"' when 'linpred.index' is ",
288         "assigned more than a single value")
289
290  if (length(linpred.index) &&
291      length(object@control$xij))
292    stop("Currently cannot handle 'xij' models when ",
293         "'linpred.index' is assigned a value")
294
295
296  x   <- slot(object, "x")
297
298
299  Xm2 <- if (any(slotNames(object) == "Xm2"))
300         slot(object, "Xm2") else numeric(0)
301
302
303 form2 <- if (any(slotNames(object) == "misc"))
304          object@misc$form2 else NULL
305  if (type == "lm2" && !length(form2))
306    return(Xm2)
307
308
309  if (!length(x)) {
310    data <- model.frame(object, xlev = object@xlevels, ...)
311
312    kill.con <- if (length(object@contrasts))
313                object@contrasts else NULL
314
315    x <- vmodel.matrix.default(object, data = data,
316                               contrasts.arg = kill.con)
317    tt <- terms(object)
318    attr(x, "assign") <- attrassigndefault(x, tt)
319  }
320
321
322
323
324  if ((type == "lm2" || type == "bothlmlm2") &&
325      !length(Xm2)) {
326    object.copy2 <- object
327    data <- model.frame(object.copy2,
328                        xlev = object.copy2@xlevels, ...)
329
330    kill.con <- if (length(object.copy2@contrasts))
331                object.copy2@contrasts else NULL
332
333    Xm2 <- vmodel.matrix.default(object.copy2, data = data,
334                                 contrasts.arg = kill.con)
335
336
337    if (length(form2)) {
338      attr(Xm2, "assign") <- attrassigndefault(Xm2, terms(form2))
339    }
340  }
341
342
343
344
345
346  if (type == "lm" && !length(linpred.index)) {
347    return(x)
348  } else if (type == "lm2") {
349    return(Xm2)
350  } else if (type == "bothlmlm2") {
351    return(list(X = x, Xm2 = Xm2))
352  }
353
354
355  M <- object@misc$M  # Number of linear/additive predictors
356  Hlist <- object@constraints  # == constraints(object, type = "lm")
357  X.vlm <- lm2vlm.model.matrix(x = x, Hlist = Hlist, Xm2 = Xm2,
358                               xij = object@control$xij)
359
360  if (type == "vlm" && !length(linpred.index))
361    return(X.vlm)
362
363
364
365  if (!is.Numeric(linpred.index,  # length.arg = 1,  20190625
366                  integer.valued = TRUE, positive = TRUE))
367    stop("bad input for argument 'linpred.index'")
368  if (!length(intersect(linpred.index, 1:M)))
369    stop("argument 'linpred.index' should have ",
370         "values from the set 1:", M)
371
372  Hlist <- Hlist
373  n.lm <- nobs(object, type = "lm")  # Number of rows of the LM matrix
374  Hmatrices <- abs(constraints(object, matrix = TRUE))  # by column
375  jay <- linpred.index
376  vecTF2 <- colSums(Hmatrices[jay, , drop = FALSE]) != 0
377  index2 <- which(vecTF2)
378  vecTF1 <- rep_len(FALSE, M)
379  vecTF1[jay] <- TRUE
380  X.lm.jay <- X.vlm[vecTF1, vecTF2, drop = FALSE]  # Recycling
381
382
383
384  aasgn.copy <- aasgn <- attr(X.vlm,  'assign')
385  vasgn.copy <- vasgn <- attr(X.vlm, 'vassign')
386  for (iloc in length(vasgn.copy):1) {
387    if (!any(is.element(index2, vasgn[[iloc]])))
388      vasgn[[iloc]] <- NULL
389  }
390
391
392  kasgn <- vasgn
393  i.start <- 1
394  for (jay in seq_len(length(vasgn))) {
395    tmp4 <- kasgn[[jay]]
396    kasgn[[jay]] <- i.start:(i.start + length(tmp4) - 1)
397    i.start <- i.start + length(tmp4)
398  }
399    attr(X.lm.jay, "vassign") <- kasgn
400    attr(X.lm.jay, "rm.vassign") <- vasgn  # Some elts removed
401
402
403
404
405  nasgn <- names(aasgn)
406  all.union <- NULL
407  for (iloc in 1:length(aasgn.copy)) {
408    if (length(intersect(aasgn[[iloc]], index2)))
409      all.union <- c(all.union, nasgn[iloc])
410  }
411  all.union <- unique(all.union)
412  ans4 <- aasgn[all.union]
413  for (iloc in 1:length(ans4)) {
414    tmp4 <- ans4[[iloc]]
415    ans4[[iloc]] <- intersect(tmp4, index2)
416  }
417  ptr.start <- 1
418  for (iloc in 1:length(ans4)) {
419    tmp4 <- ans4[[iloc]]
420    ans4[[iloc]] <- ptr.start:(ptr.start + length(tmp4) - 1)
421    ptr.start <- ptr.start + length(tmp4)
422  }
423
424
425
426
427  attr(X.lm.jay, "assign") <- ans4
428  attr(X.lm.jay, "rm.assign") <- aasgn[all.union]  # Some elts gone
429
430
431
432
433  X.lm.jay
434}  # model.matrixvlm
435
436
437
438
439
440
441
442
443
444
445
446setMethod("model.matrix",  "vlm", function(object, ...)
447           model.matrixvlm(object, ...))
448
449
450
451
452 model.matrixvgam <-
453  function(object,
454           type = c("lm", "vlm", "lm", "lm2", "bothlmlm2"),
455           linpred.index = NULL,
456           ...) {
457  model.matrixvlm(object = object,
458                  type = type[1],
459                  linpred.index = linpred.index, ...)
460}
461setMethod("model.matrix",  "vgam", function(object, ...)
462           model.matrixvgam(object, ...))
463
464
465
466
467
468
469 model.framevlm <- function(object,
470                            setupsmart = TRUE,
471                            wrapupsmart = TRUE, ...) {
472
473  dots <- list(...)
474  nargs <- dots[match(c("data", "na.action", "subset"), names(dots), 0)]
475  if (length(nargs) || !length(object@model)) {
476    fcall <- object@call
477    fcall$method <- "model.frame"
478    fcall[[1]] <- as.name("vlm")
479
480    fcall$smart <- FALSE
481    if (setupsmart && length(object@smart.prediction)) {
482      setup.smart("read", smart.prediction=object@smart.prediction)
483    }
484
485    fcall[names(nargs)] <- nargs
486    env <- environment(object@terms$terms)  # @terms or @terms$terms ??
487    if (is.null(env))
488      env <- parent.frame()
489    ans <- eval(fcall, env, parent.frame())
490
491    if (wrapupsmart && length(object@smart.prediction)) {
492      wrapup.smart()
493    }
494    ans
495  } else object@model
496}  # model.framevlm
497
498
499if (!isGeneric("model.frame"))
500    setGeneric("model.frame", function(formula, ...)
501        standardGeneric("model.frame"))
502
503setMethod("model.frame",  "vlm", function(formula, ...)
504           model.framevlm(object = formula, ...))
505
506
507
508
509
510 vmodel.matrix.default <-
511  function(object, data = environment(object),
512           contrasts.arg = NULL, xlev = NULL, ...) {
513
514  t <- if (missing(data)) terms(object) else
515                          terms(object, data = data)
516  if (is.null(attr(data, "terms")))
517    data <- model.frame(object, data, xlev = xlev) else {
518    reorder <- match(sapply(attr(t, "variables"), deparse,
519                     width.cutoff = 500)[-1], names(data))
520    if (anyNA(reorder))
521      stop("model frame and formula mismatch in model.matrix()")
522    if (!identical(reorder, seq_len(ncol(data))))
523      data <- data[, reorder, drop = FALSE]
524  }
525  int <- attr(t, "response")
526  if (length(data)) {
527    contr.funs <- as.character(getOption("contrasts"))
528    namD <- names(data)
529    for (i in namD) if (is.character(data[[i]])) {
530      data[[i]] <- factor(data[[i]])
531      warning(gettextf("variable '%s' converted to a factor", i),
532              domain = NA)
533    }
534    isF <- sapply(data, function(x) is.factor(x) || is.logical(x))
535    isF[int] <- FALSE
536    isOF <- sapply(data, is.ordered)
537    for (nn in namD[isF])
538      if (is.null(attr(data[[nn]], "contrasts")))
539        contrasts(data[[nn]]) <- contr.funs[1 + isOF[nn]]
540    if (!is.null(contrasts.arg) && is.list(contrasts.arg)) {
541      if (is.null(namC <- names(contrasts.arg)))
542        stop("invalid 'contrasts.arg' argument")
543      for (nn in namC) {
544        if (is.na(ni <- match(nn, namD)))
545          warning(gettextf(
546            "variable '%s' is absent, its contrast will be ignored",
547            nn), domain = NA) else {
548          ca <- contrasts.arg[[nn]]
549          if (is.matrix(ca))
550            contrasts(data[[ni]], ncol(ca)) <- ca else
551            contrasts(data[[ni]]) <- contrasts.arg[[nn]]
552        }
553      }
554    }
555  } else {
556    isF <- FALSE
557    data <- list(x = rep(0, nrow(data)))
558  }
559
560
561  ans  <-          (model.matrix(t, data))
562
563
564
565
566  cons <- if (any(isF))
567    lapply(data[isF], function(x) attr(x, "contrasts")) else NULL
568  attr(ans, "contrasts") <- cons
569  ans
570}  # vmodel.matrix.default
571
572
573
574
575
576depvar.vlm <-
577  function(object,
578           type = c("lm", "lm2"),
579           drop = FALSE,
580           ...) {
581  type <- match.arg(type, c("lm", "lm2"))[1]
582  ans <- if (type == "lm") {
583    object@y
584  } else {
585    object@Ym2
586  }
587  ans[, , drop = drop]
588}
589
590
591
592if (!isGeneric("depvar"))
593    setGeneric("depvar",
594               function(object, ...)
595                 standardGeneric("depvar"),
596               package = "VGAM")
597
598
599setMethod("depvar",  "vlm", function(object, ...)
600           depvar.vlm(object, ...))
601setMethod("depvar",  "rrvglm", function(object, ...)
602           depvar.vlm(object, ...))
603setMethod("depvar",  "qrrvglm", function(object, ...)
604           depvar.vlm(object, ...))
605setMethod("depvar",  "rrvgam", function(object, ...)
606           depvar.vlm(object, ...))
607setMethod("depvar",  "rcim", function(object, ...)
608           depvar.vlm(object, ...))
609
610
611
612
613
614npred.vlm <- function(object,
615                      type = c("total", "one.response"),
616                      ...) {
617  if (!missing(type))
618    type <- as.character(substitute(type))
619  type.arg <- match.arg(type, c("total", "one.response"))[1]
620
621
622  MM <-
623    if (length(object@misc$M))
624      object@misc$M else
625    if (NCOL(predict(object)) > 0)
626      NCOL(predict(object)) else
627    stop("cannot seem to obtain 'M'")
628
629
630  if (type.arg == "one.response") {
631    M1.infos <- NULL
632    infos.fun <- object@family@infos
633    Ans.infos <- infos.fun()
634    if (is.list(Ans.infos) && length(Ans.infos$M1))
635      M1.infos <- Ans.infos$M1
636
637    Q1 <- Ans.infos$Q1
638    if (is.numeric(Q1)) {
639      S <- ncol(depvar(object)) / Q1  # No. of (multiple) responses
640      if (is.numeric(M1.infos) && M1.infos * S != MM)
641        warning("contradiction in values after computing it two ways")
642    }
643
644
645    M1 <- if (is.numeric(M1.infos)) M1.infos else
646          if (is.numeric(MM      )) MM       else
647          stop("failed to compute 'M'")
648    M1
649  } else {  # One response is assumed, by default
650    MM
651  }
652}  # npred.vlm
653
654
655if (!isGeneric("npred"))
656    setGeneric("npred", function(object, ...) standardGeneric("npred"),
657               package = "VGAM")
658
659
660setMethod("npred",  "vlm", function(object, ...)
661           npred.vlm(object, ...))
662setMethod("npred",  "rrvglm", function(object, ...)
663           npred.vlm(object, ...))
664setMethod("npred",  "qrrvglm", function(object, ...)
665           npred.vlm(object, ...))
666setMethod("npred",  "rrvgam", function(object, ...)
667           npred.vlm(object, ...))
668setMethod("npred",  "rcim", function(object, ...)
669           npred.vlm(object, ...))
670
671
672
673
674
675
676hatvaluesvlm <-
677  function(model,
678           type = c("diagonal", "matrix", "centralBlocks"), ...) {
679
680
681  if (!missing(type))
682    type <- as.character(substitute(type))
683  type.arg <- match.arg(type,
684                        c("diagonal", "matrix", "centralBlocks"))[1]
685
686
687  qrSlot <- model@qr
688
689  if (!is.list(qrSlot) && class(qrSlot) != "qr")
690    stop("slot 'qr' should be a list")
691
692  M  <- npred(model)
693  nn <- nobs(model, type = "lm")
694
695  if (is.empty.list(qrSlot)) {
696
697    wzedd <- weights(model, type = "working")
698  UU <- vchol(wzedd, M = M, n = nn, silent = TRUE)
699    X.vlm <- model.matrix(model, type = "vlm")
700    UU.X.vlm <- mux111(cc = UU, xmat = X.vlm, M = M)
701    qrSlot <- qr(UU.X.vlm)
702  } else {
703    X.vlm <- NULL
704    class(qrSlot) <- "qr" # S3 class
705  }
706  Q.S3 <- qr.Q(qrSlot)
707
708
709
710  if (type.arg == "diagonal") {
711    Diag.Hat <- rowSums(Q.S3^2)
712    Diag.Elts <- matrix(Diag.Hat, nn, M, byrow = TRUE)
713
714    if (length(model@misc$predictors.names) == M)
715      colnames(Diag.Elts) <- model@misc$predictors.names
716    if (length(rownames(model.matrix(model, type = "lm"))))
717      rownames(Diag.Elts) <- rownames(model.matrix(model,
718                                                   type = "lm"))
719
720    attr(Diag.Elts, "predictors.names") <- model@misc$predictors.names
721    attr(Diag.Elts, "ncol.X.vlm") <- model@misc$ncol.X.vlm
722
723    Diag.Elts
724  } else if (type.arg == "matrix") {
725    all.mat <- Q.S3 %*% t(Q.S3)
726    if (!length(X.vlm))
727      X.vlm <- model.matrix(model, type = "vlm")
728    dimnames(all.mat) <- list(rownames(X.vlm), rownames(X.vlm))
729
730    attr(all.mat, "M") <- M
731    attr(all.mat, "predictors.names") <- model@misc$predictors.names
732    attr(all.mat, "ncol.X.vlm") <- model@misc$ncol.X.vlm
733
734    all.mat
735  } else {
736    ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
737    MMp1d2 <- M * (M + 1) / 2
738    all.rows.index <- rep((0:(nn-1))*M, rep(MMp1d2, nn))+ind1$row.index
739    all.cols.index <- rep((0:(nn-1))*M, rep(MMp1d2, nn))+ind1$col.index
740
741    H.ss <- rowSums(Q.S3[all.rows.index, ] *
742                    Q.S3[all.cols.index, ])
743
744    H.ss <- matrix(H.ss, nn, MMp1d2, byrow = TRUE)
745    H.ss
746  }
747}  # hatvaluesvlm
748
749
750
751
752
753setMethod("hatvalues",  "vlm", function(model, ...)
754           hatvaluesvlm(model, ...))
755setMethod("hatvalues",  "vglm", function(model, ...)
756           hatvaluesvlm(model, ...))
757setMethod("hatvalues",  "rrvglm", function(model, ...)
758           hatvaluesvlm(model, ...))
759setMethod("hatvalues",  "qrrvglm", function(model, ...)
760           hatvaluesvlm(model, ...))
761setMethod("hatvalues",  "rrvgam", function(model, ...)
762           hatvaluesvlm(model, ...))
763setMethod("hatvalues",  "rcim", function(model, ...)
764           hatvaluesvlm(model, ...))
765
766
767
768
769
770
771
772
773hatplot.vlm <-
774  function(model, multiplier = c(2, 3),
775           lty = "dashed",
776           xlab = "Observation",
777           ylab = "Hat values",
778           ylim = NULL, ...) {
779
780  if (is(model, "vlm")) {
781    hatval <- hatvalues(model, diag = TRUE)
782  } else {
783    hatval <- model
784  }
785
786  if (!is.matrix(hatval))
787    stop("argument 'model' seems not a vglm() object or a matrix")
788
789  ncol.X.vlm <- attr(hatval, "ncol.X.vlm")
790  M <- attr(hatval, "M")
791  predictors.names <- attr(hatval, "predictors.names")
792  if (!length(predictors.names)) {
793    predictors.names <- param.names("Linear/additive predictor ", M)
794  }
795
796  if (length(M)) {
797    N <- nrow(hatval) / M
798    hatval <- matrix(hatval, N, M, byrow = TRUE)
799  } else {
800    M <- ncol(hatval)
801    N <- nrow(hatval)
802  }
803
804  if (is.null(ylim))
805    ylim <- c(0, max(hatval))
806  for (jay in 1:M) {
807    plot(hatval[, jay], type = "n", main = predictors.names[jay],
808         ylim = ylim, xlab = xlab, ylab = ylab,
809         ...)
810    points(1:N, hatval[, jay], ...)
811    abline(h = multiplier * ncol.X.vlm / (N * M), lty = lty, ...)
812  }
813}  # hatplot.vlm
814
815
816
817
818if (!isGeneric("hatplot"))
819    setGeneric("hatplot", function(model, ...)
820      standardGeneric("hatplot"), package = "VGAM")
821
822
823setMethod("hatplot",  "matrix", function(model, ...)
824           hatplot.vlm(model, ...))
825
826setMethod("hatplot",  "vlm", function(model, ...)
827           hatplot.vlm(model, ...))
828setMethod("hatplot",  "vglm", function(model, ...)
829           hatplot.vlm(model, ...))
830
831setMethod("hatplot",  "rrvglm", function(model, ...)
832           hatplot.vlm(model, ...))
833setMethod("hatplot",  "qrrvglm", function(model, ...)
834           hatplot.vlm(model, ...))
835setMethod("hatplot",  "rrvgam", function(model, ...)
836           hatplot.vlm(model, ...))
837setMethod("hatplot",  "rcim", function(model, ...)
838           hatplot.vlm(model, ...))
839
840
841
842
843
844
845
846
847dfbetavlm <-
848  function(model,
849           maxit.new = 1,
850           trace.new = FALSE,
851           smallno = 1.0e-8,
852           ...) {
853
854  if (!is(model, "vlm"))
855    stop("argument 'model' does not seem to be a vglm() object")
856
857  n.lm <- nobs(model, type = "lm")
858  X.lm <- model.matrix(model, type = "lm")
859  X.vlm <- model.matrix(model, type = "vlm")
860  p.vlm <- ncol(X.vlm)  # nvar(model, type = "vlm")
861  M    <- npred(model)
862  etastart <- predict(model)
863  offset <- matrix(model@offset, n.lm, M)
864  new.control <- model@control
865  pweights <- weights(model, type = "prior")
866  orig.w <- if (is.numeric(model@extra$orig.w))
867              model@extra$orig.w else 1
868  y.integer <- if (is.logical(model@extra$y.integer))
869                 model@extra$y.integer else FALSE
870  coef.model <- coef(model)
871
872
873  new.control$trace <- trace.new
874  new.control$maxit <- maxit.new
875
876  dfbeta <- matrix(0, n.lm, p.vlm)
877
878  Terms.zz <- NULL
879
880
881
882
883
884  for (ii in 1:n.lm) {
885    if (trace.new) {
886      cat("\n", "Observation ", ii, "\n")
887      flush.console()
888    }
889
890    w.orig <- if (length(orig.w) != n.lm)
891                rep_len(orig.w, n.lm) else
892                orig.w
893    w.orig[ii] <- w.orig[ii] * smallno  # Relative
894
895    fit <- vglm.fit(x = X.lm,
896                X.vlm.arg = X.vlm,  # Should be more efficient
897                y = if (y.integer)
898                round(depvar(model) * c(pweights) / c(orig.w)) else
899                     (depvar(model) * c(pweights) / c(orig.w)),
900                w = w.orig,  # Set to zero so that it is 'deleted'.
901                Xm2 = NULL, Ym2 = NULL,
902                etastart = etastart,  # coefstart = NULL,
903                offset = offset,
904                family = model@family,
905                control = new.control,
906             criterion =  new.control$criterion,  # "coefficients",
907                qr.arg = FALSE,
908                constraints = constraints(model, type = "term"),
909                extra = model@extra,
910                Terms = Terms.zz,
911                function.name = "vglm")
912
913    dfbeta[ii, ] <- coef.model - fit$coeff
914  }
915
916
917  dimnames(dfbeta) <- list(rownames(X.lm), names(coef.model))
918  dfbeta
919}  # dfbetavlm
920
921
922
923
924
925
926setMethod("dfbeta",  "matrix", function(model, ...)
927           dfbetavlm(model, ...))
928
929
930setMethod("dfbeta",  "vlm", function(model, ...)
931           dfbetavlm(model, ...))
932setMethod("dfbeta",  "vglm", function(model, ...)
933           dfbetavlm(model, ...))
934
935
936setMethod("dfbeta",  "rrvglm", function(model, ...)
937           dfbetavlm(model, ...))
938setMethod("dfbeta",  "qrrvglm", function(model, ...)
939           dfbetavlm(model, ...))
940setMethod("dfbeta",  "rrvgam", function(model, ...)
941           dfbetavlm(model, ...))
942setMethod("dfbeta",  "rcim", function(model, ...)
943           dfbetavlm(model, ...))
944
945
946
947
948
949hatvaluesbasic <- function(X.vlm,
950                           diagWm,
951                           M = 1) {
952
953
954  if (M  > 1)
955    stop("currently argument 'M' must be 1")
956
957  nn <- nrow(X.vlm)
958  ncol.X.vlm <- ncol(X.vlm)
959
960  XtW <- t(c(diagWm) * X.vlm)
961
962
963  UU <- sqrt(diagWm)  # Only for M == 1
964  UU.X.vlm <- c(UU) * X.vlm # c(UU) okay for M==1
965
966  qrSlot <- qr(UU.X.vlm)
967  Rmat <- qr.R(qrSlot)
968
969  rinv <- diag(ncol.X.vlm)
970  rinv <- backsolve(Rmat, rinv)
971
972
973  Diag.Hat <- if (FALSE) {
974    covun <- rinv %*% t(rinv)
975    rhs.mat <- covun %*% XtW
976    colSums(t(X.vlm) * rhs.mat)
977  } else {
978    mymat <- X.vlm %*% rinv
979    rowSums(diagWm * mymat^2)
980  }
981  Diag.Hat
982}  # hatvaluesbasic
983
984
985
986
987
988
989 model.matrixpvgam <-
990  function(object,
991           type = c("vlm", "lm", "lm2", "bothlmlm2",
992                    "augmentedvlm", "penalty"),  # This line is new
993           linpred.index = NULL,
994           ...) {
995  type <- match.arg(type, c("vlm", "lm", "lm2", "bothlmlm2",
996                            "augmentedvlm", "penalty"))[1]
997
998  if (type == "augmentedvlm" ||
999      type == "penalty") {
1000    rbind(if (type == "penalty") NULL else
1001          model.matrixvlm(object, type = "vlm",
1002                          linpred.index = linpred.index, ...),
1003          get.X.VLM.aug(constraints  = constraints(object,
1004                                                   type = "term"),
1005                        sm.osps.list = object@ospsslot$sm.osps.list))
1006  } else {
1007    model.matrixvlm(object, type = type,
1008                    linpred.index = linpred.index, ...)
1009  }
1010}  # model.matrixpvgam
1011
1012
1013
1014setMethod("model.matrix",  "pvgam", function(object, ...)
1015           model.matrixpvgam(object, ...))
1016
1017
1018
1019
1020
1021
1022