R/olmm-methods.R

In vcrpart: Tree-Based Varying Coefficient Regression for Generalized Linear and Ordinal Mixed Models

## --------------------------------------------------------- #
## Author:       Reto Burgin
## E-Mail:       rbuergin@gmx.ch
## Date:         2019-12-15
##
## Description:
## methods for olmm objects.
##
## anova:       Likelihood-ratio tests for the comparison of
##              models
## coef, coefficients: Extract model coefficients
## deviance:    -2*Log-likelihood at the ML estimator
## drop1:       Drop single fixed effects
## estfun:      Negative scores
## extractAIC:  Extract the AIC
## fitted:      Extract fitted values from the model
## fixef:       Extract fixed effect parameters
## formula:     Extracts 'formula'
## gefp:        Extract cumulated decorrelated score process
## getCall:     Extracts 'call'
## logLik:      Log-likelihood at the ML estimator
## model.frame: Model frame (all needed variables)
## model.matrix: Model matrix (for the fixed effects)
## predict:     Predict from the fitted model
## print:       Print summary output (method for olmm and
##              olmm.summary objects)
## ranef:       Extract predicted random effects
## ranefCov:    Covariance-matrix of random effect terms
## resid, residuals Extract different types of residuals from the
##              fitted model
## show:        Print summary output (method for olmm and
##              olmm.summary objects)
## simulate:    Simulate responses based on a fitted model
## summary:     Extract summary information
## terms, terms.olmm: Extracting the terms of the model frame
##              for fixed effects
## update:      Refits a model
## VarCorr, print.VarCorr.olmm: Extract variance and standard
##              deviation of random effects and their correlation
## vcov:        Variance-covariance matrix for fixed effect parameters
## weights:     Weights
##
## Modifications:
## 2025-08-11: Implement parametric bootstrap for anova
## 2024-05-03: rename 'gefp.olmm' to 'olmm_gefp' and 'olmm_estfun' to 'olmm_estfun'
## 2019-12-15: modify checks for classes (newly use function 'inherits')
## 2017-08-19: prettified code
## 2016-10-31: checked new implementation of C-code
## 2016-02-22: removed 'rdig' argument from 'VarCorr' method
## 2014-01-16: - improve 'predict.olmm' function
## 2014-12-07: - add argument 'center' to 'predecor_control'
## 2014-10-24: - improve simulate.olmm
##             - improved 'olmm_estfun' call in 'gefp.olmm'
## 2014-10-23: - fix bug in predict.olmm
## 2014-09-22: - (internal) change 'Ninpute' to 'Nimpute' in olmm_estfun
## 2014-09-20: - use tile case in titles
## 2014-09-08: - partial substitution of 'rep' by 'rep.int'
##             - replace 'do.call' by 'call' in 'resid.olmm'
## 2013-03-17: changed many methods to S3 methods (as in lme4)
## 2013-09-06: modify formula() method. Now the formula slot
##             is called
## 2013-09-06: add S3 for terms() method
## 2013-09-06: add drop1() method
##
## To do:
## - improve update method
## - plot methods
## - olmm_estfun: handle equal zero random effects
## - anova with a single model
## --------------------------------------------------------- #

anova.olmm <- function(object,  ..., boot = FALSE, boot.nsim = 199, boot.type = c("parametric"),
                       boot.mc.cores = 1) {

  if (!is.logical(boot) && length(boot) == 1)
    stop("'boot' must be a logical of length 1.")
  if (boot) {
    if (!(is.numeric(boot.nsim) && length(boot.nsim) == 1 && as.integer(boot.nsim) >= 1))
      stop("'boot.nsim' must be an integer >= 1 of length 1.")
    boot.nsim <- as.integer(boot.nsim)
    boot.type <- match.arg(boot.type)
    if (!(is.numeric(boot.mc.cores) && length(boot.mc.cores) == 1 && as.integer(boot.mc.cores) >= 1))
      stop("'boot.mc.cores' must be an integer >= 1 of length 1.")
}
  mCall <- match.call(expand.dots = TRUE)
  dots <- list(...)
  modp <- if (length(dots) > 0) {
    sapply(dots, is, "olmm")
  } else {
    logical(0)
  }

  if (any(modp)) {

    ## multiple models
    opts <- dots[!modp]
    mods <- c(list(object), dots[modp])
    names(mods) <- sapply(
      X = as.list(mCall)[c(2, (2 + seq.int(dots))[modp])],
      as.character)
    mods <- mods[order(sapply(lapply(mods, logLik),
                              attr, "df"))]
    calls <- lapply(mods, getCall)
    data <- lapply(calls, "[[", "data")
    if (any(data != data[[1]]))
      stop("all models must be fit to the same data object")
    header <- paste("Data:", data[[1]])
    subset <- lapply(calls, "[[", "subset")
    if (any(subset != subset[[1]]))
      stop("all models must use the same subset")
    if (!is.null(subset[[1]]))
      header <-
      c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
    llks <- lapply(mods, logLik)
    Df <- sapply(llks, attr, "df")
    llk <- unlist(llks)
    chisq <- 2 * pmax(0, c(NA, diff(llk)))
    dfChisq <- c(NA, diff(Df))
    val <- data.frame(
      Df = Df,
      AIC = sapply(llks, AIC),
      BIC = sapply(llks, BIC),
      logLik = llk,
      "Chisq" = chisq,
      "Chi Df" = dfChisq,
      "Pr(>Chisq)" = pchisq(
        chisq,
        dfChisq,
        lower.tail = FALSE),
      row.names = names(mods), check.names = FALSE)
    class(val) <- c("anova", class(val))
    attr(val, "heading") <-
      c(header, "Models:", paste(rep(names(mods), times = unlist(lapply(lapply(lapply(calls, "[[", "formula"), deparse), length))), unlist(lapply(lapply(calls, "[[", "formula"), deparse)), sep = ": "))

    if (boot) {

      chisq.boot <- mclapply(
        X = seq.int(boot.nsim),
        FUN = function(i) {
          data <- get(x = mods[[1]]$call$data, envir = parent.frame())
          data[, all.vars(mods[[1]]$formula)[1L]] <- simulate(
            object = mods[[1]], nsim = 1L)[, 1L, drop = TRUE]
          llk <- sapply(
            X = mods,
            FUN = function(m) {
              mNew <- try(update(m, data = data), silent = TRUE)
              if (!inherits(mNew, "try-error")) {
                return(as.numeric(logLik(mNew)))
              } else {
                return(as.numeric(NA))
              }
            })
          return(pmax(0, -2 * (llk[1] - llk[-1])))
        },
        mc.cores = boot.mc.cores)
      chisq.boot <- do.call("rbind", chisq.boot)
      if (all(is.na(chisq.boot)))
        stop("bootstrap computations failed.")
      if (any(is.na(chisq.boot)))
        warning(paste0(
          "During the bootstrap calculations, there were some errors in the ",
          "model estimates. These are ignored in the calculations."))
      val[2:nrow(val), "Pr(>Chisq)"] <- sapply(
        X = seq.int(ncol(chisq.boot)),
        FUN = function(j) {
          if (all(is.na(chisq.boot[, j]))) return(as.numeric(NA))
          return(
            sum(chisq.boot[!is.na(chisq.boot[, j]), j] >= chisq[j + 1]) /
              sum(!is.na(chisq.boot[, j])))
        }
      )
      attr(val, "chisq.boot") <- chisq.boot

    }


    return(val)

  } else {
    ## single model
    stop("single argument anova for 'olmm' objects not yet implemented")
  }
}


coef.olmm <- function(object, which = c("all", "fe"), ...) {

  which <- match.arg(which)
  if (which == "fe") return(fixef(object))

  dims <- object$dims
  if (object$family$family == "adjacent") {
    T <- diag(dims["nPar"])
    subsRows <- seq(1, dims["pCe"] * (dims["nEta"] - 1L), 1L)
    subsCols <- seq(dims["pCe"] + 1L,
                    dims["pCe"] * dims["nEta"], 1L)
    if (length(subsRows) == 1L) {
      T[subsRows, subsCols] <- c(-1.0)
    } else {
      diag(T[subsRows, subsCols]) <- c(-1.0)
    }
    rval <- c(T %*% object$coefficients)
    names(rval) <- names(object$coefficients)
  } else {
    rval <- object$coefficients
  }
  if (dims["hasRanef"] == 0L)
    rval <- rval[!grepl("ranefCholFac", names(rval))]

  return(rval)
}


coefficients.olmm <- coef.olmm


deviance.olmm <- function(object, ...)
  return(-as.numeric(2.0 * logLik(object)))


predecor_control <- function(impute = TRUE, seed = NULL,
                             symmetric = TRUE, center = FALSE,
                             reltol = 1e-6, maxit = 250L, minsize = 1L,
                             include = c("observed", "all"),
                             verbose = FALSE, silent = FALSE) {
  stopifnot(is.logical(impute) && length(impute) == 1L)
  stopifnot(is.null(seed) | is.numeric(seed) && length(seed) == 1L)
  stopifnot(is.logical(symmetric) && length(symmetric) == 1L)
  stopifnot(is.logical(center) && length(center) == 1L)
  stopifnot(is.numeric(reltol) && reltol > 0 && length(reltol) == 1L)
  stopifnot(is.numeric(maxit) && maxit > 0 && length(maxit) == 1L)
  stopifnot(is.numeric(minsize) && minsize > 0 && length(minsize) == 1L)
  include <- match.arg(include)
  stopifnot(is.logical(verbose) && length(verbose) == 1L)
  stopifnot(is.logical(silent) && length(silent) == 1L)
  return(structure(list(impute = impute, seed = seed,
                        symmetric = symmetric, center = center,
                        reltol = reltol, maxit = maxit,
                        minsize = minsize, include = include,
                        verbose = verbose, silent = silent),
                   class = "predecor_control"))
}


olmm_estfun <- function(x, predecor = FALSE, control = predecor_control(),
                        nuisance = NULL, ...) {

  ## append '...' arguments to control
  cArgs <- list(...)
  cArgs <- cArgs[intersect(names(cArgs), names(formals(olmm_control)))]
  cArgsNames <- names(cArgs)
  cArgs <- do.call("predecor_control", cArgs)
  control[cArgsNames] <- cArgs[cArgsNames]

  if (control$verbose) cat("* extract original scores ... ")

  ## check 'x'
  stopifnot(inherits(x, "olmm"))
  Ni <- table(x$subject)
  Nmax <- as.integer(max(Ni))
  xOld <- x # store the original model for restoring at the end

  ## check 'predecor'
  if (x$dims["hasRanef"] == 0L) predecor <- FALSE
  stopifnot(is.logical(predecor))
  predecor <- predecor[1L]

  ## check 'control'
  stopifnot(inherits(control,"predecor_control"))

  ## check 'nuisance'
  if (is.character(nuisance))
    nuisance <- which(colnames(x$score_obs) %in% nuisance)
  stopifnot(all(nuisance %in% seq_along(colnames(x$score_obs))))
  parm <- seq_along(x$coefficients) # internal variable
  if (!is.null(nuisance) & is.character(nuisance))
    nuisance <- which(names(coef(x)) %in% nuisance)
  nuisance <- sort(union(nuisance, which(x$restricted)))
  parm <- setdiff(parm, nuisance)
  attr <- list() # default attributes

  scores <- x$score_obs

  if (control$verbose) cat("OK")

  ## impute data

  subsImp <- rep.int(FALSE, nrow(scores))
  if (predecor && any(Ni != Nmax)) {

    Nimpute <- Nmax - Ni
    subsImp <- c(rep.int(FALSE, x$dims["n"]), rep.int(TRUE, sum(Nimpute)))
    sbjImp <- factor(rep.int(names(Ni), Nimpute), names(Ni))
    ranef <- ranef(x)
    ranefImp <- ranef[rownames(ranef) %in% unique(sbjImp),,drop = FALSE]

    ## get predictors from empirical distribution
    yName <- deparse(lhs(formula(x)))
    yLevs <- levels(x$y)
    newFrame <- x$frame[rep.int(1L, sum(Nimpute)),,drop=FALSE]
    newFrame[, x$subjectName] <- rep.int(names(Nimpute), Nimpute)
    newX <- x$X[rep.int(1L, sum(Nimpute)),,drop=FALSE]
    newW <- x$W[rep.int(1L, sum(Nimpute)),,drop=FALSE]

    ## add imputations to model
    x$frame <- rbind(x$frame, newFrame)
    x$y <- ordered(c(as.character(x$y),
                     as.character(newFrame[, yName])), yLevs)
    x$X <- rbind(x$X, newX)
    x$W <- rbind(x$W, newW)
    x$subject <-
      factor(c(as.character(x$subject), newFrame[, x$subjectName]),
             levels = names(Ni))
    x$weights <- x$weights_sbj[as.integer(x$subject)]
    x$offset <- rbind(x$offset, matrix(0.0, sum(Nimpute), x$dims["nEta"]))
    x$dims["n"] <- nrow(x$frame)
    x$eta <- rbind(x$eta, matrix(0.0, sum(Nimpute), x$dims["nEta"]))
    x$score_obs <- rbind(
      x$score_obs, matrix(0.0, sum(Nimpute), x$dims["nPar"]))

    ## simulate responses
    if (control$verbose) cat("\n* impute scores ... ")

    ## set seed
    if (!is.null(control$seed)) set.seed(control$seed)

    if (control$impute) {

      ## impute predictors
      times <- Nimpute[x$subject[!subsImp]]
      rows <- unlist(tapply(1:sum(Ni), x$subject[!subsImp], function(x) sample(x, times[x[1L]], replace = TRUE)))
      x$frame[subsImp,] <- x$frame[rows,,drop=FALSE]
      x$X[subsImp, ] <- x$X[rows,,drop=FALSE]
      x$W[subsImp, ] <- x$W[rows,,drop=FALSE]

      ## draw responses
      subsW <- c(rep(which(attr(xOld$W, "merge") == 1L),
                     x$dims["nEta"]),
                 which(attr(xOld$W, "merge") == 2L))
      tmatW <- rbind(kronecker(diag(x$dims["nEta"]),
                               rep(1,x$dims["qCe"])),
                     matrix(1, x$dims["qGe"], x$dims["nEta"]))
      etaFixef <- x$X[subsImp, ] %*% x$fixef
      etaRanef <- (x$W[subsImp, subsW,drop = FALSE] *
                     ranef[as.integer(x$subject[subsImp])]) %*% tmatW
      eta <- etaFixef + etaRanef
      probs <- x$family$linkinv(eta)
      x$y[subsImp] <- # simulate responses
        ordered(apply(probs, 1L, function(x) sample(yLevs, 1L, prob = x)), yLevs)

      ## recompute scores
      new <-
        .Call("olmm_update_marg", x,
              x$coefficients, PACKAGE = "vcrpart")
      x <- modifyList(x, new)
    }

    scores <- x$score_obs

    if (control$center && max(abs(cSums <- colSums(scores))) > 1e-6)
      scores <- scores -
      matrix(cSums / nrow(scores),
             nrow(scores), ncol(scores), byrow = TRUE)
  }

  ## drop the nuisance coefficients
  scores <- scores[, parm, drop = FALSE]

  if (predecor) {

    ## compute transformation matrix
    if (control$verbose) cat("\n* compute transformation matrix ...")

    subsT <- if (control$include == "observed")
      !subsImp else rep(TRUE, nrow(scores))
    T <- olmm_decormat(scores = scores[subsT,,drop = FALSE],
                       subject = x$subject[subsT],
                       control = control)

    ## if transformation failed, return raw scores
    if (attr(T, "conv") > 0L) {

      if (control$verbose) cat("\n* transforming scores ... ")

      ## transformation matrix for one subject
      Ti <- kronecker(matrix(1, Nmax, Nmax) - diag(Nmax), T)
      diag(Ti) <- 1
      subsOrd <- order(x$subject)
      sTmp <- matrix(c(t(scores[subsOrd,,drop=FALSE])),
                     Nmax * ncol(scores),length(Ni))
      sTmp <- matrix(c(Ti %*% sTmp),
                     nrow(scores), ncol(scores), byrow = TRUE)
      sTmp <- sTmp[order(subsOrd),,drop = FALSE]
      scores[] <- sTmp

      if (control$verbose) cat("OK")
    }
    scores <- subset(scores, !subsImp)
    attr(scores, "T") <- T
  }

  ## ^hack: recompute old model
  x <- xOld
  new <- .Call("olmm_update_marg", x, x$coefficients, PACKAGE = "vcrpart")
  x <- modifyList(x, new)

  if (control$verbose) cat("\n* return negative scores\n")

  ## return scores
  return(-scores)
}


## thanks lme4
extractAIC.olmm <- function(fit, scale, k = 2, ...) {
  L <- logLik(fit)
  edf <- attr(L,"df")
  c(edf,-2*L + k*edf)
}


fitted.olmm <- function(object, ...) vcrpart_fitted(object, ...)


fixef.olmm <- function(object, which = c("all", "ce", "ge"), ...) {

  which <- match.arg(which)
  dims <- object$dims
  coef <- coef(object)
  rval <- c()

  ## category-specific coefficients
  if (which %in% c("all", "ce") && dims["pCe"] > 0) {
    subs <- seq(from  = 1, to = dims["pCe"] * dims["nEta"], by = 1)
    rval <- c(rval, coef[subs])
  }

  ## global coefficients
  if (which %in% c("all", "ge") && dims["pGe"] > 0) {
    subs <- seq(from = dims["pCe"] * dims["nEta"] + 1,
                to = dims["pCe"] * dims["nEta"] + dims["pGe"],
                by = 1)
    rval <- c(rval, coef[subs])
  }

  return(rval)
}


formula.olmm <- function(x, ...) as.formula(x$formula, env = parent.frame())


olmm_gefp <- function(object, scores = NULL, order.by = NULL, subset = NULL,
                      predecor = TRUE, parm = NULL, center = TRUE,
                      drop = TRUE, silent = FALSE, ...) {

  ## extract scores (if scores is not a matrix)
  if (is.null(scores)) {
    estfunCall <- list(name = as.name("olmm_estfun"),
                       x = quote(object),
                       predecor = quote(predecor))
    dotargs <- list(...)
    dotargs <- dotargs[intersect(names(formals(olmm_estfun)), names(dotargs))]
    estfunCall[names(dotargs)] <- dotargs
    mode(estfunCall) <- "call"
    scores <- try(eval(estfunCall))
  } else if (is.function(scores)) {
    scores <- scores(object)
  } else if (is.matrix(scores)) {
    if (!silent && !is.null(attr(scores, "predecor")) &&
        as.integer(predecor) != attr(scores, "predecor"))
      warning("'scores' is not pre-decorrelated")
  }

  if (!is.matrix(scores)) stop("extracting the score function failed.")

  ## set 'order.by'
  if (is.null(order.by)) order.by <- 1:nobs(object)
  if (is.factor(order.by)) order.by <- droplevels(order.by)

  ## set subset
  if (!is.null(subset)) {
    if (is.character(subset)) subset <- rownames(scores) %in% subset
    if (is.numeric(subset)) subset <- (1:nobs(object)) %in% subset
  } else {
    subset <- rep.int(TRUE, nobs(object))
  }
  subsScores <- rownames(model.frame(object)) %in% rownames(scores)

  ## create process
  process <- scores[subset[subsScores],,drop=FALSE]
  cn <- colnames(process)
  order.by <- order.by[subset & subsScores]

  ## get dimensions
  n <- nrow(scores)
  k <- ncol(scores)

  ## if necessary, subtract the column means
  if (center & max(abs(cMeans <- colMeans(process))) > 1e-6)
    process <- process - matrix(cMeans, nrow(process), ncol(process), byrow = TRUE)

  ## scale scores by the number of observations
  process <- process / sqrt(n)

  ## multiply scores with the inverse of the square root of their crossproduct
  subs <- rep.int(TRUE, k)
  J12 <- crossprod(process)
  J12Inv <- matrix(0, k, k)
  if (drop)
    subs <- subs & apply(process, 2, function(x) max(abs(x)) > .Machine$double.eps)
  mat <- try(chol2inv(chol(root.matrix(J12[subs, subs, drop = FALSE]))), TRUE)
  if (inherits(mat, "try-error") && drop) {
    subs <- subs & (diag(J12) / max(diag(J12)) > 1e-6)
    mat <- try(chol2inv(chol(root.matrix(J12[subs, subs, drop = FALSE]))), TRUE)
  }
  if (inherits(mat, "try-error") && drop && length(parm) > 0L) {
    subs <- subs & colnames(process) %in% parm
    mat <- try(chol2inv(chol(root.matrix(J12[subs, subs, drop = FALSE]))), TRUE)
  }

  if (inherits(mat, "try-error")) return(mat)

  if (!silent && !all(subs))
    warning("covariance matrix is not positive semidefinite. Omit terms: ",
            paste(cn[!subs], collapse = ", "))

  J12Inv[subs, subs] <- mat
  process <- t(J12Inv %*% t(process))

  ## order and cumulate the process
  index <- order(order.by)
  process <- apply(process[index, , drop = FALSE], 2, cumsum)
  process <- rbind(0, process)
  colnames(process) <- cn

  ## transform process to a multivariate time series
  time <- order.by[index]
  if (is.factor(time)) time <- as.numeric(droplevels(time))
  time <- suppressWarnings(c(time[1] - as.numeric(diff(time[1:2])), time))

  ## extract terms
  if (!is.null(parm)) process <- process[, parm, drop = FALSE]

  ## return a list of class "gefp"
  rval <- list(process = suppressWarnings(zoo(process, time)),
               nreg = k,
               nobs = n,
               call = match.call(),
               fit = NULL,
               scores = NULL,
               fitted.model = getCall(object),
               par = NULL,
               lim.process = "Brownian bridge",
               type.name = "M-fluctuation test",
               order.name = deparse(substitute(order.by)),
               subset = rownames(model.frame(object))[subset & subsScores],
               J12 = NULL)
  class(rval) <- "gefp"
  return(rval)
}


getCall.olmm <- function(x, ...) return(x$call)


logLik.olmm <- function(object, ...) {
  dims <- object$dims
  rval <- object$logLik
  attr(rval, "nall") <- attr(rval, "nobs") <- dims[["n"]]
  nPar <- dims[["nPar"]] - (1 - dims[["hasRanef"]])
  attr(rval, "df") <- nPar
  class(rval) <- "logLik"
  return(rval)
}


model.frame.olmm <- function(formula, ...) formula$frame


model.matrix.olmm <- function(object, which = c("fe", "fe-ce", "fe-ge",
                                                "re", "re-ce", "re-ge"), ...) {
  which <- match.arg(which)
  rval <- switch(substr(which, 1L, 2L),
                 fe = object$X,
                 re = object$W)
  if (!which %in% c("fe", "re")) {
    attr <- attributes(rval)
    subs <- attr(rval, "merge") == switch(substr(which, 4, 5), ce = 1, ge = 2)
    attr$assign <- attr$assign[subs]
    attr$merge <- attr$merge[subs]
    attr$orig.colnames <- attr$orig.colnames[subs]
    rval <- rval[, subs, drop = FALSE]
    attributes(rval) <- appendDefArgs(attributes(rval), attr)
  }
  return(rval)
}

neglogLik2.olmm <- function(object, ...)
  return(-2 * as.numeric(logLik(object, ...)))


nobs.olmm <- function(object, ...) object$dims[["n"]]


predict.olmm <- function(object, newdata = NULL,
                         type = c("link", "response", "prob",
                                  "class", "ranef"),
                         ranef = FALSE, na.action = na.pass, ...) {

  ## extract data
  type <- match.arg(type) # retrieve type

  ## resolve conflicts with the 'ranef' argument
  if (type == "ranef" & !is.null(newdata))
    stop("prediction for random effects for 'newdata' is not ",
         "implemented.")
  if (type == "ranef") return(ranef(object, ...))

  if (type == "prob") type <- "response"
  formList <-
    vcrpart_formula(formula(object), object$family) # extract formulas

  offset <- list(...)$offset
  subset <- list(...)$subset
  dims <- object$dims

  ## checks
  if (!is.null(newdata) && !inherits(newdata, "data.frame"))
    stop("'newdata' must be a 'data.frame'.")
  if (!(inherits(ranef, "logical") | inherits(ranef, "matrix")))
    stop("'ranef' must be a 'logical' or a 'matrix'.")

  if (dims["hasRanef"] < 1L) {
    ranef <- FALSE
    formList$re$eta$ge <- as.formula(~ 1)
    formList$re$eta$ce <- as.formula(~ -1)
  }

  if (is.null(newdata)) {

    ## extract data from the model fit
    X <- model.matrix(object, "fe")
    W <- model.matrix(object, "re")
    subject <- object$subject
    offset <- object$offset

    ## check and extract random effects
    if (is.logical(ranef)) {
      if (ranef) ranef <- ranef(object)
    } else {
      if (any(dim(ranef) != c(dims["N"], dims["qEta"])))
        stop("'ranef' matrix has wrong dimensions")
    }

  } else {

    ## data preparation
    mfForm <- formList$allTerms
    if (!object$subjectName %in% colnames(newdata))
      mfForm <- update(mfForm, paste(". ~ . -", object$subjectName))
    mf <- model.frame(object)
    Terms <- delete.response(terms(mfForm))
    xlevels <- .getXlevels(attr(mf, "terms"), mf)
    xlevels <- xlevels[names(xlevels) %in%  all.vars(Terms)]
    xlevels <- xlevels[names(xlevels) != object$subjectName]

    newdata <- as.data.frame(model.frame(Terms, newdata,
                                         na.action = na.action,
                                         xlev = xlevels))

    if (!is.null(cl <- attr(Terms, "dataClasses")))
      .checkMFClasses(cl, mf)

    ## extract fixed effect model matrix from newdata
    X <- olmm_merge_mm(model.matrix(terms(formList$fe$eta$ce, keep.order = TRUE), newdata, attr(object$X, "contrasts")[intersect(all.vars(formList$fe$eta$ce), names(attr(object$X, "contrasts")))]), model.matrix(terms(formList$fe$eta$ge, keep.order = TRUE), newdata, attr(object$X, "contrasts")[intersect(all.vars(formList$fe$eta$ge), names(attr(object$X, "contrasts")))]), TRUE)
    rownames(X) <- rownames(newdata)

    ## delete columns of dropped terms
    X <- X[, colnames(object$X), drop = FALSE]

    if (is.null(offset))
      offset <- matrix(0.0, nrow(X), dims["nEta"])

    if (is.logical(ranef) && ranef) {

      if (object$subjectName %in% colnames(newdata)) {
        subjLevs <- unique(newdata[, object$subjectName])
        ranef <- matrix(0.0, length(subjLevs), ncol(object$u),
                        dimnames =
                          list(subjLevs, colnames(object$u)))
        if (!all(subjLevs %in% levels(object$subject)))
          message(paste("set random effects of new subjects",
                        paste(setdiff(subjLevs,
                                      levels(object$subject)),
                              collapse = ", "), "to 0."))
        ranef[intersect(rownames(ranef), rownames(object$u)), ] <-
          ranef(object)[intersect(rownames(ranef),
                                  rownames(object$u)),,
                        drop = FALSE]
      } else {
        message("set random effects of new subjects to 0.")
        ranef <- matrix(0.0, nrow(X), ncol(object$u),
                        dimnames = list(
                          paste("New", 1:nrow(X), sep = ""),
                          colnames(object$u)))
      }
    }

    if (is.matrix(ranef) && ncol(ranef) != ncol(object$u))
      stop("random effects matrix has wrong dimensions")

    ## random effects
    if (object$subjectName %in% colnames(newdata)) {

      subject <- factor(newdata[, object$subjectName])

      if (is.matrix(ranef)) {
        if (any(!levels(subject) %in% rownames(ranef))) {
          stop(paste("random effects missing for subjects",
                     paste(setdiff(levels(subject), rownames(ranef)),
                           collapse = ", ")))
        }
        ## !!! has problems with ids = ""
        ranef <- ranef[rownames(ranef) %in% levels(subject),,drop = FALSE]
      }

      subject <- factor(subject,
                        levels = unique(c(levels(object$subject),
                                          levels(subject))))

    } else {

      subject <- factor(paste("New", 1:nrow(X), sep = ""))

    }

    ## extract model formulas
    W <- olmm_merge_mm(model.matrix(terms(
      formList$re$eta$ce, keep.order = TRUE),
      newdata, attr(object$W, "contrasts")),
      model.matrix(terms(formList$re$eta$ge, keep.order = TRUE),
                   newdata, attr(object$W, "contrasts")),
      FALSE)
    rownames(W) <- rownames(newdata)
  }

  if (!is.null(subset)) {
    X <- X[subset, , drop = FALSE]
    W <- W[subset, , drop = FALSE]
    if (!is.null(subject)) subject <- subject[subset]
    offset <- offset[subset,,drop = FALSE]
    if (is.matrix(ranef))
      ranef <-
      ranef[rownames(ranef) %in% levels(subject), , drop = FALSE]
  }

  yLevs <- levels(object$y)
  if (nrow(X) > 0) {

    ## compute linear predictor
    eta <- offset + X %*% object$fixef

    ## predict marginal ...
    if (is.logical(ranef) && !ranef &
        type %in% c("response", "class")) {

      if (any(subject %in% levels(object$subject))) {

        ## cluster-averaged expectation (works also if person is
        ## not present)
        probs <- matrix(0, nrow(X), ncol(eta) + 1L)
        colnames(probs) <- levels(object$y)
        rownames(probs) <- rownames(X)
        probs <-
          .Call("olmm_pred_margNew", object, eta, W, subject,
                nrow(X), probs, PACKAGE = "vcrpart")

      } else {

        ## population-averaged expectation
        probs <- matrix(0, nrow(X), ncol(eta) + 1L)
        colnames(probs) <- levels(object$y)
        rownames(probs) <- rownames(X)
        probs <-
          .Call("olmm_pred_marg", object, eta, W, nrow(X), probs,
                PACKAGE = "vcrpart")

      }


    } else {
      ## or conditional expected probabilities (or linear predictor)

      subsW <- c(rep(which(attr(W, "merge") == 1L), dims["nEta"]),
                 which(attr(W, "merge") == 2L))
      tmatW <- rbind(kronecker(diag(dims["nEta"]),
                               rep.int(1,dims["qCe"])),
                     matrix(1, dims["qGe"], dims["nEta"]))

      ## extend linear predictor
      if (is.matrix(ranef))
        eta <- eta +
        (W[, subsW,drop = FALSE] *
           ranef[as.character(subject),,drop=FALSE]) %*% tmatW
      rownames(eta) <- rownames(X)
      colnames(eta) <- colnames(object$eta)

      if (type == "link") return(eta)

      probs <- object$family$linkinv(eta)
      colnames(probs) <- yLevs
      rownames(probs) <- rownames(X)
    }

    if (type == "response") {

      ## return probabilites
      rval <- probs
    } else if (type == "class") {

      ## return most probable category
      rval <- apply(probs, 1, which.max)
      rval <- ordered(yLevs[rval], levels = yLevs)
      names(rval) <- rownames(X)
    }

  } else { # useful for predict.tvcm
    if (type == "response") {
      rval <- matrix(, 0, dims["J"], dimnames = list(NULL, yLevs))
    } else if (type == "class") {
      rval <- NULL
    }
  }
  return(rval)
}


print.olmm <- function(x, etalab = c("int", "char", "eta"), ...) {

  etalab <- match.arg(etalab)

  so <- summary.olmm(x, silent = TRUE)

  if (length(so$methTitle) > 0) cat(so$methTitle, "\n\n")
  if (length(so$family) > 0)
    cat(" Family:", so$family$family, so$family$link, "\n")
  if (length(so$formula) > 0) cat("Formula:", so$formula, "\n")
  if (length(so$data) > 0) cat("   Data:", so$data, "\n")
  if (length(so$subset) > 0) cat(" Subset:", so$subset ,"\n")

  if (length(so$na.action) > 0L) { cat("\n"); cat(so$na.action, "\n"); }

  if (length(so$AICtab) > 0) {
    cat("\nGoodness of fit:\n")
    print(so$AICtab, ...)
  }

  if (length(so$REmat) > 0) {
    cat("\nRandom effects:\n")
    print.VarCorr.olmm(olmm_rename(
      so$REmat, so$yLevs, so$family, etalab), ...)
    cat(sprintf("Number of obs: %d, subjects: %d\n",
                so$dims["n"], so$dims["N"]))
  }

  if (length(so$feMatGe) > 0 && nrow(so$feMatGe) > 0) {
    cat("\nGlobal fixed effects:\n")
    text <- so$feMatGe[, 1]; names(text) <- rownames(so$feMatGe);
    print(text, ...)
  }

  if (length(so$feMatCe) > 0 && nrow(so$feMatCe) > 0) {
    cat("\nCategory-specific fixed effects:\n")
    text <- so$feMatCe[, 1]; names(text) <- rownames(so$feMatCe);
    print(olmm_rename(text, so$yLevs, so$family, etalab), ...)
  }
}


ranef.olmm <- function(object, norm = FALSE, ...) {
  if (!norm) {
    rval <- object$u %*% t(object$ranefCholFac)
  } else {
    rval <- object$u
  }
  return(rval)
}


ranefCov.olmm <- function(object, ...) {

  dims <- object$dims

  ## !!! 2025-07-24: needs to be checked !!!
  ## transformation for the adjacent-categories model
  if (object$family$family == "adjacent" & dims["qCe"] > 0) {

    ranefCholFac <- rval <- object$ranefCholFac

    ## row-wise subtraction of category-specific effects
    for (i in 1L:dims["qCe"]) {
      subs <- seq(i, dims["qCe"] * dims["nEta"], dims["qCe"])
      for (j in 1L:(dims["nEta"] - 1L)) {
        rval[subs[j], ] <- ranefCholFac[subs[j], ] -
          ranefCholFac[subs[j + 1L], ]
      }
    }

    for (i in 1L:(dims["nEta"] - 1L)) {
      rval[dims["qCe"]*(i-1)+1:dims["qCe"], ] <-
        ranefCholFac[dims["qCe"]*(i-1L)+1:dims["qCe"], ] -
        ranefCholFac[dims["qCe"]*i+1L:dims["qCe"], ]
    }
    return(rval %*% t(rval))

  } else {

    ## cumulative-link model or baseline-category model
    return(object$ranefCholFac %*% t(object$ranefCholFac))
  }
}


resid.olmm <- function(object, norm = FALSE, ...) {
  call <- list(name = as.name("predict"),
               object = quote(object), type = "response")
  call <- appendDefArgs(call, list(...))
  mode(call) <- "call"
  fitted <- eval(call)
  y <- as.integer(model.response(model.frame(object)))
  J <- object$dims["J"]
  n <- length(y)
  rval <- sapply(1:n, function(i) {
    sum(fitted[i, 1L:J > y[i]]) - sum(fitted[i, 1L:J < y[i]])
  })
  if (norm) {
    var <- (1.0 - apply(fitted^3, 1, sum)) / 3.0
    rval <- rval / sqrt(var)
  }
  return(rval)
}


residuals.olmm <- resid.olmm


simulate.olmm <- function(object, nsim = 1, seed = NULL,
                          newdata = NULL, ranef = TRUE,
                          ranef.simulate = TRUE, ...) {

  dotArgs <- list(...)
  if (!is.null(seed)) set.seed(seed)
  if (!exists(".Random.seed", envir = .GlobalEnv)) runif(1)
  if (!is.logical(ranef.simulate)) stop("'ranef.simulate' must be a single logical.")
  RNGstate <- .Random.seed
  dotArgs$type <- "response"

  ## simulate random effects
  if (object$dims["hasRanef"] && is.logical(ranef) && ranef && ranef.simulate) {
    if (is.null(newdata)) {
      ranef <- matrix(
        data = rnorm(nrow(object$u) * ncol(object$u)),
        nrow = nrow(object$u), ncol = ncol(object$u),
        dimnames = list(rownames(object$u), colnames(object$u)))
      ranef <- ranef %*% t(object$ranefCholFac)

    } else {
      if (!object$subjectName %in% colnames(newdata))
        stop(paste0(
          "'newdata' must contain a column with name '",
          object$subjectName, "'", sep = ""))
      if (!is.factor(newdata[, object$subjectName]))
        stop(paste0(
          "column '", object$subjectName, "' must be a factor."))
      subject <- droplevels(newdata[, object$subjectName])
      # simulate random effects
      ranef <- matrix(
        data = rnorm(nlevels(subject) * ncol(object$u)),
        nrow = nlevels(subject), ncol = ncol(object$u),
        dimnames = list(levels(subject), colnames(object$u)))
      ranef <- ranef %*% t(object$ranefCholFac)
    }
  }
  pred <- predict(object, newdata = newdata, type = "prob", ranef = ranef, ...)
  FUN <- function(x) sample(levels(object$y), 1, prob = x)
  rval <- as.data.frame(replicate(nsim, apply(pred, 1L, FUN)))
  for (i in 1:nsim)
    rval[,i] <-
    factor(rval[, i], levels = levels(object$y), ordered = TRUE)
  if (nsim == 1) {
    colnames(rval) <-
      colnames(model.frame(object))[1]
  } else {
    colnames(rval) <-
      paste(colnames(model.frame(object))[1], 1:nsim, sep = ".")
  }
  attr(rval, "seed") <- RNGstate
  return(rval)
}


summary.olmm <- function(object, etalab = c("int", "char", "eta"),
                         silent = FALSE, ...) {

  etalab <- match.arg(etalab)
  dims <- object$dims

  ## goodness of fit measures
  lLik <- logLik(object)
  AICtab <- data.frame(AIC = AIC(lLik),
                       BIC = BIC(lLik),
                       logLik = as.vector(lLik),
                       row.names = "")

  ## fixed-effect coefficients
  fixef <- fixef(object)

  ## fixed-effect coefficient-covariance matrix
  vcov <- try(vcov(object), silent = TRUE)
  validVcov <- !inherits(vcov, "try-error") && min(diag(vcov)) > 0
  if (!silent && !validVcov)
    warning("computation of variance-covariance matrix failed")

  ## global fixed effects
  if (dims["pGe"] > 0) {
    subs <- seq(dims["pCe"] * dims["nEta"] + 1L,
                dims["pCe"] * dims["nEta"] + dims["pGe"], 1L)
    feMatGe <-
      cbind("Estimate" = fixef[subs],
            "Std. Error" = rep.int(NaN, length(subs)),
            "z value" = rep.int(NaN, length(subs)))
    if (validVcov) {
      feMatGe[, 2L] <- sqrt(diag(vcov)[subs])
      feMatGe[, 3L] <- feMatGe[, 1L] / feMatGe[, 2L]
    }

  } else { # empty matrix
    feMatGe <- matrix(, 0L, 3L, dimnames = list(c(), c("Estimate", "Std. Error", "z value")))
  }

  ## category-specific fixed effects
  if (dims["pCe"] > 0L) {
    subs <- seq(1L, dims["pCe"] * dims["nEta"], 1L)
    feMatCe <-
      cbind("Estimate" = fixef[subs],
            "Std. Error" = rep.int(NaN, length(subs)),
            "z value" = rep.int(NaN, length(subs)))
    if (validVcov) {
      feMatCe[, 2L] <- sqrt(diag(vcov)[subs])
      feMatCe[, 3L] <- feMatCe[, 1L] / feMatCe[, 2L]
    }
  } else { # empty matrix
    feMatCe <- matrix(, 0L, 3L, dimnames = list(c(), c("Estimate", "Std. Error", "z value")))
  }

  ## random effects
  if (dims["hasRanef"] > 0L) {
    VarCorr <- VarCorr(object)
  } else {
    VarCorr <-
      matrix(, 0L, 3L,
             dimnames = list(c(), c("Variance", "StdDev", "")))
  }

  ## title
  methTitle <- "Linear"
  if (dims["hasRanef"] > 0L) methTitle <- paste(methTitle, "Mixed")
  methTitle <- paste(methTitle, "Model")
  if (dims["hasRanef"] > 0L)
    methTitle <- paste(methTitle, " fit by Marginal Maximum\n",
                       "Likelihood with Gauss-Hermite Quadrature",
                       sep = "")

  na.action <- naprint(attr(model.frame(object), "na.action"))
  na.action <- if (na.action == "") character() else paste("(", na.action, ")", sep = "")
  call <- getCall(object)

  yLevs <- if (is.factor(object$y)) levels(object$y) else 1L

  ## return a 'summary.olmm' object
  return(structure(
    list(methTitle = methTitle,
         family = object$family,
         formula = paste(deparse(formula(object)), collapse = "\n"),
         data = deparseCall(call$data),
         subset = deparseCall(call$subset),
         AICtab = AICtab,
         feMatCe = feMatCe,
         feMatGe = feMatGe,
         REmat = VarCorr,
         na.action = na.action,
         dims = dims,
         yLevs = levels(object$y),
         etalab = etalab,
         dotargs = list(...)), class = "summary.olmm"))
}


print.summary.olmm <- function(x, ...) {

  args <- appendDefArgs(list(...), x$dotargs)

  if (length(x$methTitle) > 0L) cat(x$methTitle, "\n\n")
  if (length(x$family) > 0L) cat(" Family:",
                                 x$family$family, x$family$link, "\n")
  if (length(x$formula) > 0L) cat("Formula:", x$formula, "\n")
  if (length(x$data) > 0L) cat("   Data:", x$data, "\n")
  if (length(x$subset) > 0L) cat(" Subset:", x$subset ,"\n")

  if (length(x$AICtab) > 0L) {
    cat("\nGoodness of fit:\n")
    args$x <- x$AICtab
    do.call("print", args)
  }

  if (length(x$REmat) > 0L) {
    cat("\nRandom effects:\n")
    args$x <- olmm_rename(x$REmat, x$yLevs, x$family, x$etalab)
    do.call("print", args)
    cat(sprintf("Number of obs: %d, subjects: %d\n",
                x$dims["n"], x$dims["N"]))
  }

  if (length(x$na.action) > 0L) cat(x$na.action, "\n")

  if (length(x$feMatGe) > 0 && nrow(x$feMatGe) > 0L) {
    cat("\nGlobal fixed effects:\n")
    args$x <- x$feMatGe
    do.call("printCoefmat", args)
  }

  if (length(x$feMatCe) > 0 && nrow(x$feMatCe) > 0L) {
    cat("\nCategory-specific fixed effects:\n")
    args$x <- olmm_rename(x$feMatCe, x$yLevs, x$family, x$etalab)
    do.call("printCoefmat", args)
  }
}


terms.olmm <- function(x, which = c("fe-ce", "fe-ge",
                                    "re-ce", "re-ge"), ...) {
  which <- match.arg(which)
  which <- switch(which,
                  "fe-ge" = "feGe",
                  "fe-ce" = "feCe",
                  "re-ge" = "reGe",
                  "re-ce" = "reCe")
  return(x$terms[[which]])
}


update.olmm <- function(object, formula., evaluate = TRUE, ...) {

  call <- getCall(object)
  if (is.null(call))
    stop("need an object with call slot")
  extras <- match.call(expand.dots = FALSE)$...
  if (!missing(formula.))
    call$formula <- update.formula(formula(object), formula.)
  if (length(extras) > 0) {
    existing <- !is.na(match(names(extras), names(call)))
    for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
    if (any(!existing)) {
      call <- c(as.list(call), extras[!existing])
      call <- as.call(call)
    }
  }
  if (evaluate)
    eval(call, parent.frame())
  else call
}


VarCorr.olmm <- function(x, sigma = 1., ...) {

  ## create formatted output according to VarCorr
  RECovMat <- ranefCov(x)
  REmat <- cbind(Variance = diag(RECovMat),
                 StdDev = sqrt(diag(RECovMat)))
  rval <- cbind(REmat, cov2cor(RECovMat))
  attr(rval, "title") <- paste("Subject:", x$subjectName)
  class(rval) <- "VarCorr.olmm"
  return(rval)
}


print.VarCorr.olmm <- function(x, ...) { # S3 method

  rval <- format(x, ...)

  if (nrow(rval) > 1L) {

    ## prettify Corr matrix
    Corr <- rval[, 3L:ncol(rval)]
    Corr[upper.tri(Corr)] <- ""
    diag(Corr) <- colnames(Corr)
    Corr <- Corr[, -nrow(Corr), drop = FALSE]
    dimnames(Corr) <- NULL
    colnames(Corr) <- c("Corr", rep.int("", nrow(Corr) - 2L))
    rval <- cbind(rval[, 1L:2L, drop = FALSE], Corr)
  } else {

    ## random intercept models need no correlation terms
    rval <- rval[, 1L:2L, drop = FALSE]
  }

  ## print the output
  if (!is.null(attr(x, "title"))) {
    cat(attr( x, "title" ), "\n")
    attr(x, "title") <- NULL
  }
  print(rval, quote = FALSE)
  invisible(x)
}


vcov.olmm <- function(object, ...) {

  dims <- object$dims
  info <- object$info
  if (dims["hasRanef"] == 0L)
    info <- info[1L:dims["p"], 1:dims["p"]]

  ## extract inverse of negative info-matrix
  rval <- chol2inv(chol(-info))
  dimnames(rval) <- dimnames(info)

  ## parameter transformation for adjacent-category models
  if (object$family$family == "adjacent") {

    ## matrix T with partial derivates of transformation
    T <- diag(nrow(rval))
    subsRows <- seq(1, dims["pCe"] * (dims["nEta"] - 1L), 1L)
    subsCols <- seq(dims["pCe"] + 1L,
                    dims["pCe"] * dims["nEta"], 1L)
    if (length(subsRows) == 1L) {
      T[subsRows, subsCols] <- c(-1.0)
    } else {
      diag(T[subsRows, subsCols]) <- c(-1.0)
    }
    # subsRows <- seq(dims["pCe"] + 1L,
    #                 dims["pCe"] * dims["nEta"], 1L)
    # subsCols <- seq(1,dims["pCe"] * dims["nEta"], 1L)
    # if (length(subsRows) == 1L) {
    #     T[subsRows, subsCols] <- c(-1.0)
    # } else {
    #     diag(T[subsRows, subsCols]) <- c(-1.0)
    # }

    ## transform covariance-matrix
    rval <- T %*% rval %*% t(T)
    dimnames(rval) <- dimnames(info)
  }

  return(rval)
}


weights.olmm <- function(object, level = c("observation", "subject"), ...) {
  return(switch(match.arg(level),
                observation = object$weights,
                subject = object$weights_sbj))
}