R/model.matrix.vglm.q

In VGAM: Vector Generalized Linear and Additive Models

# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.













attrassigndefault <- function(mmat, tt) {
  if (!inherits(tt, "terms"))
    stop("need terms object")
  aa <- attr(mmat, "assign")
  if (is.null(aa))
    stop("argument is not really a model matrix")
  ll <- attr(tt, "term.labels")
  if (attr(tt, "intercept") > 0)
    ll <- c("(Intercept)", ll)
  aaa <- factor(aa, labels = ll)
  split(order(aa), aaa)
}


attrassignlm <- function(object, ...)
  attrassigndefault(model.matrix(object), object@terms)





 vlabel <-
  function(xn, ncolHlist, M, separator = ":", colon = FALSE) {

  if (length(xn) != length(ncolHlist))
    stop("length of first two arguments not equal")

  n1 <- rep(xn, ncolHlist)
  if (M == 1)
    return(n1)
  n2 <- as.list(ncolHlist)
  n2 <- lapply(n2, seq)
  n2 <- unlist(n2)
  n2 <- as.character(n2)
  n2 <- paste(separator, n2, sep = "")
  n3 <- rep(ncolHlist, ncolHlist)
  if (!colon)
    n2[n3 == 1] <- ""
  n1n2 <- paste(n1, n2, sep = "")
  n1n2
}  # vlabel






 vlm2lm.model.matrix <-
  function(x.vlm, Hlist = NULL,
           which.linpred = 1,
           M = NULL) {







  if (is.numeric(M)) {
    if (M != nrow(Hlist[[1]]))
      stop("argument 'M' does not match argument 'Hlist'")
  } else {
    M <- nrow(Hlist[[1]])
  }


  Hmatrices <- matrix(c(unlist(Hlist)), nrow = M)
  if (ncol(Hmatrices) != ncol(x.vlm))
    stop("ncol(Hmatrices) != ncol(x.vlm)")


  n.lm <- nrow(x.vlm) / M
  if (round(n.lm) != n.lm)
    stop("'n.lm' does not seem to be an integer")
  linpred.index <- which.linpred

  if (FALSE) {
    vecTF <- Hmatrices[linpred.index, ] != 0
    X.lm.jay <- x.vlm[(0:(n.lm - 1)) * M + linpred.index,
                      vecTF, drop = FALSE]
  }  # (FALSE)
  vecTF2 <- colSums(Hmatrices[linpred.index, ,
                              drop = FALSE]) != 0
  vecTF1 <- rep_len(FALSE, M)
  vecTF1[linpred.index] <- TRUE
  X.lm.jay <- x.vlm[vecTF1, vecTF2, drop = FALSE]  # Recycling

  X.lm.jay
}  # vlm2lm.model.matrix







 lm2vlm.model.matrix <-
  function(x, Hlist = NULL,
           assign.attributes = TRUE,
           M = NULL, xij = NULL,
           label.it = TRUE,  # Ignored for "lm".
           Xm2 = NULL) {




  if (length(Hlist) != ncol(x))
    stop("length(Hlist) != ncol(x)")

  if (length(xij)) {
    if (inherits(xij, "formula"))
      xij <- list(xij)
    if (!is.list(xij))
      stop("'xij' is not a list of formulae")
  }

  if (!is.numeric(M))
    M <- nrow(Hlist[[1]])

  nrow.X.lm <- nrow(x)
  if (all(trivial.constraints(Hlist) == 1)) {
    X.vlm <- if (M > 1) kronecker(x, diag(M)) else x
    ncolHlist <- rep(M, ncol(x))
  } else {
    allB <- matrix(unlist(Hlist), nrow = M)
    ncolHlist <- unlist(lapply(Hlist, ncol))
    Rsum <- sum(ncolHlist)

    X1 <- rep(c(t(x)), rep(ncolHlist, nrow.X.lm))
    dim(X1) <- c(Rsum, nrow.X.lm)
    X.vlm <- kronecker(t(X1), matrix(1, M, 1)) *
             kronecker(matrix(1, nrow.X.lm, 1),
                       allB)
    rm(X1)
  }



  if (label.it) {
    dn <- labels(x)
    yn <- dn[[1]]
    xn <- dn[[2]]
    dimnames(X.vlm) <- list(vlabel(yn, rep(M, nrow.X.lm), M),
                            vlabel(xn, ncolHlist, M))
  }  # label.it



  if (assign.attributes) {
    attr(X.vlm, "contrasts")   <- attr(x, "contrasts")
    attr(X.vlm, "factors")     <- attr(x, "factors")
    attr(X.vlm, "formula")     <- attr(x, "formula")
    attr(X.vlm, "class")       <- attr(x, "class")
    attr(X.vlm, "order")       <- attr(x, "order")
    attr(X.vlm, "term.labels") <- attr(x, "term.labels")

    orig.assign.lm <- attr(x, "orig.assign.lm")  # NULL if x = F
    nasgn <- oasgn <- attr(x, "assign")
    lowind <- 0
    for (ii in seq_along(oasgn)) {
      mylen <- length(oasgn[[ii]]) * ncolHlist[oasgn[[ii]][1]]
      nasgn[[ii]] <- (lowind+1):(lowind+mylen)
      lowind <- lowind + mylen
    } # End of ii
    if (lowind != ncol(X.vlm))
      stop("something gone wrong")
    attr(X.vlm, "assign") <- nasgn


    fred <- unlist(lapply(nasgn,
            length))/unlist(lapply(oasgn, length))
    vasgn <- vector("list", sum(fred))
    kk <- 0
    for (ii in seq_along(oasgn)) {
      temp <- matrix(nasgn[[ii]],
                     ncol = length(oasgn[[ii]]))
      for (jloc in 1:nrow(temp)) {
        kk <- kk + 1
        vasgn[[kk]] <- temp[jloc, ]
      }
    }
    names(vasgn) <- vlabel(names(oasgn), fred, M)
    attr(X.vlm, "vassign") <- vasgn

    attr(X.vlm, "constraints") <- Hlist


    attr(X.vlm, "orig.assign.lm") <- orig.assign.lm
  }  # End of if (assign.attributes)




  if (!length(xij))
    return(X.vlm)







  at.x <- attr(x, "assign")
  at.vlmx <- attr(X.vlm, "assign")
  at.Xm2 <- attr(Xm2, "assign")

  for (ii in seq_along(xij)) {
    form.xij <- xij[[ii]]
    if (length(form.xij) != 3)
      stop("xij[[", ii, "]] is not a formula with a response")
    tform.xij <- terms(form.xij)
    aterm.form <- attr(tform.xij, "term.labels")
    if (length(aterm.form) != M)
      stop("xij[[", ii, "]] does not contain ", M, " terms")

    name.term.y <- as.character(form.xij)[2]
    cols.X.vlm <- at.vlmx[[name.term.y]]  # May be > 1 in length.

    x.name.term.2 <- aterm.form[1]  # Choose the first one
    One.such.term <- at.Xm2[[x.name.term.2]]
    for (bbb in seq_along(One.such.term)) {
      use.cols.Xm2 <- NULL
      for (sss in 1:M) {
        x.name.term.2 <- aterm.form[sss]
        one.such.term <- at.Xm2[[x.name.term.2]]
        use.cols.Xm2 <- c(use.cols.Xm2, one.such.term[bbb])
      } # End of sss

      allXk <- Xm2[, use.cols.Xm2, drop = FALSE]
      cmat.no <- (at.x[[name.term.y]])[1]
      cmat <- Hlist[[cmat.no]]
      Rsum.k <- ncol(cmat)
      tmp44 <- kronecker(matrix(1, nrow.X.lm, 1), t(cmat)) *
       kronecker(allXk, matrix(1, ncol(cmat), 1))  # n*Rsum.k x M

      tmp44 <- array(t(tmp44), c(M, Rsum.k, nrow.X.lm))
      tmp44 <- aperm(tmp44, c(1, 3, 2))  # c(M, n, Rsum.k)
      rep.index <- cols.X.vlm[((bbb-1)*Rsum.k+1):(bbb*Rsum.k)]
      X.vlm[, rep.index] <- c(tmp44)
    }  # End of bbb
  }  # End of for (ii in seq_along(xij))

  if (assign.attributes) {
    attr(X.vlm, "vassign") <- vasgn
    attr(X.vlm, "assign") <- nasgn
    attr(X.vlm, "xij") <- xij
  }
  X.vlm
}  # lm2vlm.model.matrix









model.matrix.vlm <- function(object, ...)
  model.matrixvlm(object, ...)





 model.matrixvlm <-
  function(object,
           type = c("vlm", "lm", "lm2", "bothlmlm2"),
           linpred.index = NULL,
           label.it = TRUE,
           ...) {





  if (mode(type) != "character" && mode(type) != "name")
    type <- as.character(substitute(type))
  type <- match.arg(type, c("vlm", "lm", "lm2",
                            "bothlmlm2"))[1]

  linpred.index <- unique(sort(linpred.index))
  LLLL <- length(linpred.index)
  if (LLLL == 1 && type != "lm")
    stop("Must set 'type = \"lm\"' when 'linpred.index' is ",
         "assigned a single value")
  if (LLLL  > 1 && type != "vlm")
    stop("Must set 'type = \"vlm\"' when 'linpred.index' is ",
         "assigned more than a single value")

  if (length(linpred.index) &&
      length(object@control$xij))
    stop("Currently cannot handle 'xij' models when ",
         "'linpred.index' is assigned a value")


  x <- slot(object, "x")


  Xm2 <- if (any(slotNames(object) == "Xm2"))
         slot(object, "Xm2") else numeric(0)


 form2 <- if (any(slotNames(object) == "misc"))
          object@misc$form2 else NULL
  if (type == "lm2" && !length(form2))
    return(Xm2)


  if (!length(x)) {
    data <- model.frame(object, xlev = object@xlevels, ...)

    kill.con <- if (length(object@contrasts))
                object@contrasts else NULL

    x <- vmodel.matrix.default(object, data = data,
                               contrasts.arg = kill.con)
    tt <- terms(object)
    attr(x, "assign") <- attrassigndefault(x, tt)
  }




  if ((type == "lm2" || type == "bothlmlm2") &&
      !length(Xm2)) {
    object.copy2 <- object
    data <- model.frame(object.copy2,
                        xlev = object.copy2@xlevels, ...)

    kill.con <- if (length(object.copy2@contrasts))
                object.copy2@contrasts else NULL

    Xm2 <- vmodel.matrix.default(object.copy2, data = data,
                                 contrasts.arg = kill.con)


    if (length(form2)) {
      attr(Xm2, "assign") <- attrassigndefault(Xm2,
                                               terms(form2))
    }
  }





  if (type == "lm" && !length(linpred.index)) {
    return(x)
  } else if (type == "lm2") {
    return(Xm2)
  } else if (type == "bothlmlm2") {
    return(list(X = x, Xm2 = Xm2))
  }


  M <- object@misc$M  # Number of linear/additive predictors
  Hlist <- object@constraints  # == constraints(object,type="lm")
    X.vlm <-
      lm2vlm.model.matrix(x = x, Hlist = Hlist,
                          Xm2 = Xm2,
                          label.it = label.it,
                          xij = object@control$xij)

  if (type == "vlm" && !length(linpred.index))
    return(X.vlm)



  if (!is.Numeric(linpred.index,  # length.arg = 1,  20190625
                  integer.valued = TRUE, positive = TRUE))
    stop("bad input for argument 'linpred.index'")
  if (!length(intersect(linpred.index, 1:M)))
    stop("argument 'linpred.index' should have ",
         "values from the set 1:", M)

  Hlist <- Hlist
  n.lm <- nobs(object, type = "lm")  # nrow(the LM matrix)
  Hmatrices <- abs(constraints(object,
                               matrix = TRUE))  # by column
  jay <- linpred.index
  vecTF2 <- colSums(Hmatrices[jay, , drop = FALSE]) != 0
  index2 <- which(vecTF2)
  vecTF1 <- rep_len(FALSE, M)
  vecTF1[jay] <- TRUE
  X.lm.jay <- X.vlm[vecTF1, vecTF2, drop = FALSE]  # Recycling



  aasgn.copy <- aasgn <- attr(X.vlm,  'assign')
  vasgn.copy <- vasgn <- attr(X.vlm, 'vassign')
  for (iloc in length(vasgn.copy):1) {
    if (!any(is.element(index2, vasgn[[iloc]])))
      vasgn[[iloc]] <- NULL      
  }


  kasgn <- vasgn
  i.start <- 1
  for (jay in seq_len(length(vasgn))) {
    tmp4 <- kasgn[[jay]]
    kasgn[[jay]] <- i.start:(i.start + length(tmp4) - 1)
    i.start <- i.start + length(tmp4)
  }
    attr(X.lm.jay, "vassign") <- kasgn
    attr(X.lm.jay, "rm.vassign") <- vasgn  # Some elts removed




  nasgn <- names(aasgn)
  all.union <- NULL
  for (iloc in 1:length(aasgn.copy)) {
    if (length(intersect(aasgn[[iloc]], index2)))
      all.union <- c(all.union, nasgn[iloc])
  }
  all.union <- unique(all.union)
  ans4 <- aasgn[all.union]
  for (iloc in 1:length(ans4)) {
    tmp4 <- ans4[[iloc]]
    ans4[[iloc]] <- intersect(tmp4, index2)
  }
  ptr.start <- 1
  for (iloc in 1:length(ans4)) {
    tmp4 <- ans4[[iloc]]
    ans4[[iloc]] <- ptr.start:(ptr.start + length(tmp4) - 1)
    ptr.start <- ptr.start + length(tmp4)
  }




  attr(X.lm.jay, "assign") <- ans4
  attr(X.lm.jay, "rm.assign") <- aasgn[all.union]



  X.lm.jay
}  # model.matrixvlm











setMethod("model.matrix",  "vlm", function(object, ...)
           model.matrixvlm(object, ...))




 model.matrixvgam <-
  function(object,
           type = c("lm", "vlm", "lm", "lm2", "bothlmlm2"),
           linpred.index = NULL,
           ...) {
  model.matrixvlm(object = object,
                  type = type[1],
                  linpred.index = linpred.index, ...)
}
setMethod("model.matrix",  "vgam", function(object, ...)
           model.matrixvgam(object, ...))






 model.framevlm <- function(object,
                            setupsmart = TRUE,
                            wrapupsmart = TRUE, ...) {

  dots <- list(...)
  nargs <- dots[match(c("data", "na.action", "subset"), names(dots), 0)]
  if (length(nargs) || !length(object@model)) {
    fcall <- object@call
    fcall$method <- "model.frame"
    fcall[[1]] <- as.name("vlm")

    fcall$smart <- FALSE
    if (setupsmart && length(object@smart.prediction)) {
      setup.smart("read", smart.prediction=object@smart.prediction)
    }

    fcall[names(nargs)] <- nargs
    env <- environment(object@terms$terms)  # @terms or @terms$terms ??
    if (is.null(env))
      env <- parent.frame()
    ans <- eval(fcall, env, parent.frame())

    if (wrapupsmart && length(object@smart.prediction)) {
      wrapup.smart()
    }
    ans
  } else object@model
}  # model.framevlm


if (!isGeneric("model.frame"))
    setGeneric("model.frame", function(formula, ...)
        standardGeneric("model.frame"))

setMethod("model.frame",  "vlm", function(formula, ...)
           model.framevlm(object = formula, ...))





 vmodel.matrix.default <-
  function(object, data = environment(object),
           contrasts.arg = NULL, xlev = NULL, ...) {

  t <- if (missing(data)) terms(object) else
                          terms(object, data = data)
  if (is.null(attr(data, "terms")))
    data <- model.frame(object, data, xlev = xlev) else {
    reorder <- match(sapply(attr(t, "variables"), deparse,
                     width.cutoff = 500)[-1], names(data))
    if (anyNA(reorder))
      stop("model frame and formula mismatch in model.matrix()")
    if (!identical(reorder, seq_len(ncol(data))))
      data <- data[, reorder, drop = FALSE]
  }
  int <- attr(t, "response")
  if (length(data)) {
    contr.funs <- as.character(getOption("contrasts"))
    namD <- names(data)
    for (i in namD) if (is.character(data[[i]])) {
      data[[i]] <- factor(data[[i]])
      warning(gettextf("variable '%s' converted to a factor", i),
              domain = NA)
    }
    isF <- sapply(data, function(x) is.factor(x) || is.logical(x))
    isF[int] <- FALSE
    isOF <- sapply(data, is.ordered)
    for (nn in namD[isF])
      if (is.null(attr(data[[nn]], "contrasts")))
        contrasts(data[[nn]]) <- contr.funs[1 + isOF[nn]]
    if (!is.null(contrasts.arg) && is.list(contrasts.arg)) {
      if (is.null(namC <- names(contrasts.arg)))
        stop("invalid 'contrasts.arg' argument")
      for (nn in namC) {
        if (is.na(ni <- match(nn, namD)))
          warning(gettextf(
            "variable '%s' is absent, its contrast will be ignored",
            nn), domain = NA) else {
          ca <- contrasts.arg[[nn]]
          if (is.matrix(ca))
            contrasts(data[[ni]], ncol(ca)) <- ca else
            contrasts(data[[ni]]) <- contrasts.arg[[nn]]
        }
      }
    }
  } else {
    isF <- FALSE
    data <- list(x = rep(0, nrow(data)))
  }


  ans  <-          (model.matrix(t, data))




  cons <- if (any(isF))
    lapply(data[isF], function(x) attr(x, "contrasts")) else NULL
  attr(ans, "contrasts") <- cons
  ans
}  # vmodel.matrix.default





depvar.vlm <-
  function(object,
           type = c("lm", "lm2"),
           drop = FALSE,
           ...) {
  type <- match.arg(type, c("lm", "lm2"))[1]
  ans <- if (type == "lm") {
    object@y
  } else {
    object@Ym2
  }
  ans[, , drop = drop]
}



if (!isGeneric("depvar"))
    setGeneric("depvar",
               function(object, ...)
                 standardGeneric("depvar"),
               package = "VGAM")


setMethod("depvar",  "vlm", function(object, ...)
           depvar.vlm(object, ...))
setMethod("depvar",  "rrvglm", function(object, ...)
           depvar.vlm(object, ...))
setMethod("depvar",  "qrrvglm", function(object, ...)
           depvar.vlm(object, ...))
setMethod("depvar",  "rrvgam", function(object, ...)
           depvar.vlm(object, ...))
setMethod("depvar",  "rcim", function(object, ...)
           depvar.vlm(object, ...))





npred.vlm <- function(object,
                      type = c("total", "one.response"),
                      ...) {
  if (!missing(type))
    type <- as.character(substitute(type))
  type.arg <- match.arg(type, c("total", "one.response"))[1]


  MM <-
    if (length(object@misc$M))
      object@misc$M else
    if (NCOL(predict(object)) > 0)
      NCOL(predict(object)) else
    stop("cannot seem to obtain 'M'")


  if (type.arg == "one.response") {
    M1.infos <- NULL
    infos.fun <- object@family@infos
    Ans.infos <- infos.fun()
    if (is.list(Ans.infos) && length(Ans.infos$M1))
      M1.infos <- Ans.infos$M1

    Q1 <- Ans.infos$Q1
    if (is.numeric(Q1)) {
      S <- ncol(depvar(object)) / Q1  # No. of (multiple) responses
      if (is.numeric(M1.infos) && M1.infos * S != MM)
        warning("contradiction in values after computing it two ways")
    }


    M1 <- if (is.numeric(M1.infos)) M1.infos else
          if (is.numeric(MM      )) MM       else
          stop("failed to compute 'M'")
    M1
  } else {  # One response is assumed, by default
    MM
  }
}  # npred.vlm


if (!isGeneric("npred"))
    setGeneric("npred", function(object, ...) standardGeneric("npred"),
               package = "VGAM")


setMethod("npred",  "vlm", function(object, ...)
           npred.vlm(object, ...))
setMethod("npred",  "rrvglm", function(object, ...)
           npred.vlm(object, ...))
setMethod("npred",  "qrrvglm", function(object, ...)
           npred.vlm(object, ...))
setMethod("npred",  "rrvgam", function(object, ...)
           npred.vlm(object, ...))
setMethod("npred",  "rcim", function(object, ...)
           npred.vlm(object, ...))






hatvaluesvlm <-
  function(model,
           type = c("diagonal", "matrix", "centralBlocks"), ...) {


  if (!missing(type))
    type <- as.character(substitute(type))
  type.arg <- match.arg(type,
                        c("diagonal", "matrix", "centralBlocks"))[1]


  qrSlot <- model@qr

  if (!is.list(qrSlot) && !inherits(qrSlot, "qr"))
    stop("slot 'qr' should be a list")

  M  <- npred(model)
  nn <- nobs(model, type = "lm")

  if (is.empty.list(qrSlot)) {

    wzedd <- weights(model, type = "working")
  UU <- vchol(wzedd, M = M, n = nn, silent = TRUE)
    X.vlm <- model.matrix(model, type = "vlm")
    UU.X.vlm <- mux111(cc = UU, xmat = X.vlm, M = M)
    qrSlot <- qr(UU.X.vlm)
  } else {
    X.vlm <- NULL
    class(qrSlot) <- "qr"  # S3 class
  }
  Q.S3 <- qr.Q(qrSlot)



  if (type.arg == "diagonal") {
    Diag.Hat <- rowSums(Q.S3^2)
    Diag.Elts <- matrix(Diag.Hat, nn, M, byrow = TRUE)

    if (length(model@misc$predictors.names) == M)
      colnames(Diag.Elts) <- model@misc$predictors.names
    if (length(rownames(model.matrix(model, type = "lm"))))
      rownames(Diag.Elts) <- rownames(model.matrix(model,
                                                   type = "lm"))

    attr(Diag.Elts, "predictors.names") <- model@misc$predictors.names
    attr(Diag.Elts, "ncol.X.vlm") <- model@misc$ncol.X.vlm

    Diag.Elts
  } else if (type.arg == "matrix") {
    all.mat <- Q.S3 %*% t(Q.S3)
    if (!length(X.vlm))
      X.vlm <- model.matrix(model, type = "vlm")
    dimnames(all.mat) <- list(rownames(X.vlm), rownames(X.vlm))

    attr(all.mat, "M") <- M
    attr(all.mat, "predictors.names") <- model@misc$predictors.names
    attr(all.mat, "ncol.X.vlm") <- model@misc$ncol.X.vlm

    all.mat
  } else {
    ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
    MMp1d2 <- M * (M + 1) / 2
    all.rows.index <- rep((0:(nn-1))*M, rep(MMp1d2, nn))+ind1$row.index
    all.cols.index <- rep((0:(nn-1))*M, rep(MMp1d2, nn))+ind1$col.index

    H.ss <- rowSums(Q.S3[all.rows.index, ] *
                    Q.S3[all.cols.index, ])

    H.ss <- matrix(H.ss, nn, MMp1d2, byrow = TRUE)
    H.ss
  }
}  # hatvaluesvlm





setMethod("hatvalues",  "vlm", function(model, ...)
           hatvaluesvlm(model, ...))
setMethod("hatvalues",  "vglm", function(model, ...)
           hatvaluesvlm(model, ...))
setMethod("hatvalues",  "rrvglm", function(model, ...)
           hatvaluesvlm(model, ...))
setMethod("hatvalues",  "qrrvglm", function(model, ...)
           hatvaluesvlm(model, ...))
setMethod("hatvalues",  "rrvgam", function(model, ...)
           hatvaluesvlm(model, ...))
setMethod("hatvalues",  "rcim", function(model, ...)
           hatvaluesvlm(model, ...))








hatplot.vlm <-
  function(model, multiplier = c(2, 3),
           lty = "dashed",
           xlab = "Observation",
           ylab = "Hat values",
           ylim = NULL, ...) {

  if (is(model, "vlm")) {
    hatval <- hatvalues(model, diag = TRUE)
  } else {
    hatval <- model
  }

  if (!is.matrix(hatval))
    stop("argument 'model' seems not a vglm() object or a matrix")

  ncol.X.vlm <- attr(hatval, "ncol.X.vlm")
  M <- attr(hatval, "M")
  predictors.names <- attr(hatval, "predictors.names")
  if (!length(predictors.names)) {
    predictors.names <- param.names("Linear/additive predictor ", M)
  }

  if (length(M)) {
    N <- nrow(hatval) / M
    hatval <- matrix(hatval, N, M, byrow = TRUE)
  } else {
    M <- ncol(hatval)
    N <- nrow(hatval)
  }

  if (is.null(ylim))
    ylim <- c(0, max(hatval))
  for (jay in 1:M) {
    plot(hatval[, jay], type = "n", main = predictors.names[jay],
         ylim = ylim, xlab = xlab, ylab = ylab,
         ...)
    points(1:N, hatval[, jay], ...)
    abline(h = multiplier * ncol.X.vlm / (N * M), lty = lty, ...)
  }
}  # hatplot.vlm




if (!isGeneric("hatplot"))
    setGeneric("hatplot", function(model, ...)
      standardGeneric("hatplot"), package = "VGAM")


setMethod("hatplot",  "matrix", function(model, ...)
           hatplot.vlm(model, ...))

setMethod("hatplot",  "vlm", function(model, ...)
           hatplot.vlm(model, ...))
setMethod("hatplot",  "vglm", function(model, ...)
           hatplot.vlm(model, ...))

setMethod("hatplot",  "rrvglm", function(model, ...)
           hatplot.vlm(model, ...))
setMethod("hatplot",  "qrrvglm", function(model, ...)
           hatplot.vlm(model, ...))
setMethod("hatplot",  "rrvgam", function(model, ...)
           hatplot.vlm(model, ...))
setMethod("hatplot",  "rcim", function(model, ...)
           hatplot.vlm(model, ...))








dfbetavlm <-
  function(model,
           maxit.new = 1,
           trace.new = FALSE,
           smallno = 1.0e-8,
           ...) {

  if (!is(model, "vlm"))
    stop("argument 'model' does not seem to be a vglm() object")

  n.lm <- nobs(model, type = "lm")
  X.lm <- model.matrix(model, type = "lm")
  X.vlm <- model.matrix(model, type = "vlm")
  p.vlm <- ncol(X.vlm)  # nvar(model, type = "vlm")
  M    <- npred(model)
  etastart <- predict(model)
  offset <- matrix(model@offset, n.lm, M)
  new.control <- model@control
  pweights <- weights(model, type = "prior")
  orig.w <- if (is.numeric(model@extra$orig.w))
              model@extra$orig.w else 1
  y.integer <- if (is.logical(model@extra$y.integer))
                 model@extra$y.integer else FALSE
  coef.model <- coef(model)


  new.control$trace <- trace.new
  new.control$maxit <- maxit.new

  dfbeta <- matrix(0, n.lm, p.vlm)

  Terms.zz <- NULL





  for (ii in 1:n.lm) {
    if (trace.new) {
      cat("\n", "Observation ", ii, "\n")
      flush.console()
    }

    w.orig <- if (length(orig.w) != n.lm)
                rep_len(orig.w, n.lm) else
                orig.w
    w.orig[ii] <- w.orig[ii] * smallno  # Relative

    fit <- vglm.fit(x = X.lm,
                X.vlm.arg = X.vlm,  # Should be more efficient
                y = if (y.integer)
                round(depvar(model) * c(pweights) / c(orig.w)) else
                     (depvar(model) * c(pweights) / c(orig.w)),
                w = w.orig,  # Set to zero so that it is 'deleted'.
                Xm2 = NULL, Ym2 = NULL,
                etastart = etastart,  # coefstart = NULL,
                offset = offset,
                family = model@family,
                control = new.control,
             criterion =  new.control$criterion,  # "coefficients",
                qr.arg = FALSE,
                constraints = constraints(model, type = "term"),
                extra = model@extra,
                Terms = Terms.zz,
                function.name = "vglm")

    dfbeta[ii, ] <- coef.model - fit$coeff
  }


  dimnames(dfbeta) <- list(rownames(X.lm), names(coef.model))
  dfbeta
}  # dfbetavlm






setMethod("dfbeta",  "matrix", function(model, ...)
           dfbetavlm(model, ...))


setMethod("dfbeta",  "vlm", function(model, ...)
           dfbetavlm(model, ...))
setMethod("dfbeta",  "vglm", function(model, ...)
           dfbetavlm(model, ...))


setMethod("dfbeta",  "rrvglm", function(model, ...)
           dfbetavlm(model, ...))
setMethod("dfbeta",  "qrrvglm", function(model, ...)
           dfbetavlm(model, ...))
setMethod("dfbeta",  "rrvgam", function(model, ...)
           dfbetavlm(model, ...))
setMethod("dfbeta",  "rcim", function(model, ...)
           dfbetavlm(model, ...))





hatvaluesbasic <- function(X.vlm,
                           diagWm,
                           M = 1) {


  if (M  > 1)
    stop("currently argument 'M' must be 1")

  nn <- nrow(X.vlm)
  ncol.X.vlm <- ncol(X.vlm)

  XtW <- t(c(diagWm) * X.vlm)


  UU <- sqrt(diagWm)  # Only for M == 1
  UU.X.vlm <- c(UU) * X.vlm # c(UU) okay for M==1

  qrSlot <- qr(UU.X.vlm)
  Rmat <- qr.R(qrSlot)

  rinv <- diag(ncol.X.vlm)
  rinv <- backsolve(Rmat, rinv)


  Diag.Hat <- if (FALSE) {
    covun <- rinv %*% t(rinv)
    rhs.mat <- covun %*% XtW
    colSums(t(X.vlm) * rhs.mat)
  } else {
    mymat <- X.vlm %*% rinv
    rowSums(diagWm * mymat^2)
  }
  Diag.Hat
}  # hatvaluesbasic






 model.matrixpvgam <-
  function(object,
           type = c("vlm", "lm", "lm2", "bothlmlm2",
                    "augmentedvlm", "penalty"),  # This line is new
           linpred.index = NULL,
           ...) {
  type <- match.arg(type, c("vlm", "lm", "lm2", "bothlmlm2",
                            "augmentedvlm", "penalty"))[1]

  if (type == "augmentedvlm" ||
      type == "penalty") {
    rbind(if (type == "penalty") NULL else
          model.matrixvlm(object, type = "vlm",
                          linpred.index = linpred.index, ...),
          get.X.VLM.aug(constraints  = constraints(object,
                                                   type = "term"),
                        sm.osps.list = object@ospsslot$sm.osps.list))
  } else {
    model.matrixvlm(object, type = type,
                    linpred.index = linpred.index, ...)
  }
}  # model.matrixpvgam



setMethod("model.matrix",  "pvgam", function(object, ...)
           model.matrixpvgam(object, ...))