Nothing
## --------------------------------------------------------- #
## Author: Reto Burgin
## E-Mail: rbuergin@gmx.ch
## Date: 2021-04-15
##
## References:
## ordinal: http://cran.r-project.org/web/packages/ordinal/index.html
## lme4: http://cran.r-project.org/web/packages/lme4/index.html
## matrixcalc: http://cran.r-project.org/web/packages/matrixcalc/index.html
## statmod: http://cran.r-project.org/web/packages/statmod/index.html
##
## Dependencies:
## ucminf: http://cran.r-project.org/web/packages/ucminf/index.html
##
##
## Modifications:
## 2018-02-15: Fixed bug for 'numHess = TRUE'
## 2018-02-13: Few minor modification when exploring for the memory
## overflow problem
## 2017-08-19: - Add check for the left-hand of the formula.
## - Remove data from 'olmm' output
## 2017-08-13: Also return the data in the 'olmm' output object.
## 2017-08-11: In 'olmm' replace 'inVar <- all.vars(...)' by
## attr(terms(), "term.labels")
## 2016-11-03: modification for new C-code implementation
## 2016-04-12: modified 'olmm_control' for the 'optim' optimizier. Now
## 'BFGS' is set as the default and a warning is shown
## if a non-gradient based method is chosen while
## 'numGrad = FALSE'.
## 2015-10-30: set default 'na.action = na.omit' on 'olmm'
## 2015-09-02: started with integration auf gaussian mixed model
## 2015-01-15: improved predict.olmm
## 2014-09-25: - removed bug for numeric estimation of covariance of
## 'olmm' objects
## - define 'score_sbj' and 'score_obs' slot even if
## 'numGrad = FALSE' (otherwise olmm_update_marg gives error)
## 2014-09-19: allow 'family' to be of class 'function'
## 2014-09-08: partial substitution of 'rep' by 'rep.int'
## 2014-06-17: convert routine to S3 class
## 2014-05-03: moved several control parameters to 'control' argument
## 2014-05-02: added 'linkinv' function to families 'cumulative' etc.
## 2014-05-01: change offset argument: not it must be a 'matrix'
## 2014-04-22: implement change in 'form' object
## 2013-09-15: Free() commands were added in olmm.c
## 2013-09-07: C implementation for updating the marginal Likelihood
## and predicting random-effects was stabilized by
## replacing many alloca's by the R built-in function
## Calloc, which may slow the estimation
## 2013-07-27: change 'start' handling and add 'restricted'
## argument
## 2013-07-19: correct use of numGrad argument (from now the slots
## score_sbj and score_obs remain empty)
## 2013-07-12: improve use of contrasts. There were irritating
## warnings under correct use and now the slot
## 'contrasts' also contains contrasts from the
## model matrix for random effects
## 2021-04-15: olmm_control, when length(fit) > 1 we set fit as fit[1]
## (length(fit) > 1 caused an error in R 4.1). Fixed by GR.
##
## To do:
## - check 'nlopr' package
## - add to family 'link' and 'linkinv' and incorporate
## that in 'predict'
## - implement further family options
## - find better initial parameter values (see polr.R)
## - extract covariance matrix directly from optimizer
## - standardized coefficients
## - unconstrained covariance-matrix for random-effects
## --------------------------------------------------------- #
dev.resids <- function(y, mu, wt) {
sapply(1:nrow(y), function(i) - 2 * log(mu[i, which(y[i, ] > 0)]))
}
cumulative <- function(link = c("logit", "probit", "cauchy")) {
link <- match.arg(link)
linkinv <- function(eta) {
cumProbs <- plogis(eta)
cumProbs <- cbind(cumProbs, rep.int(1.0, nrow(eta)))
tmp <- apply(cumProbs, 1, diff)
if (ncol(cumProbs) > 2) tmp <- t(tmp)
probs <- cbind(cumProbs[, 1], tmp)
return(probs)
}
rval <- structure(list(family = "cumulative",
link = link,
linkinv = linkinv,
dev.resids = dev.resids),
class = "family.olmm")
return(rval)
}
baseline <- function(link = "logit") {
link <- match.arg(link)
linkinv <- function(eta) {
probs <- exp(eta) / (1 + matrix(rowSums(exp(eta)), nrow(eta), ncol(eta)))
probs <- cbind(probs, 1 - rowSums(probs))
return(probs)
}
rval <- structure(list(family = "baseline",
link = link,
linkinv = linkinv,
dev.resids = dev.resids),
class = "family.olmm")
return(rval)
}
adjacent <- function(link = "logit") {
link <- match.arg(link)
linkinv <- function(eta) {
probs <- exp(eta) / (1 + matrix(rowSums(exp(eta)), nrow(eta), ncol(eta)))
probs <- cbind(probs, 1 - rowSums(probs))
return(probs)
}
rval <- structure(list(family = "adjacent",
link = link,
linkinv = linkinv,
dev.resids = dev.resids),
class = "family.olmm")
return(rval)
}
print.family.olmm <- function(x, ...) {
cat("\nFamily:", x$family, "\n")
cat("Link function:", x$link, "\n\n")
invisible(x)
}
olmm_control <- function(fit = c("nlminb", "ucminf", "optim"), doFit = TRUE,
numGrad = FALSE, numHess = numGrad, nGHQ = 7L,
start = NULL, restricted = NULL, verbose = FALSE,
...) {
## intial fit vector fit = c("nlminb", "ucminf", "optim") caused an error
## if (length(fit) > 1) fit <- fit[1]
fit <- match.arg(fit)
dArgs <- list(...)
dArgs <- dArgs[intersect(names(dArgs), names(formals(fit)))]
dArgs <- dArgs[!names(dArgs) %in% c("par", "fn", "gr")]
if (fit == "optim" && is.null(dArgs$method)) {
dArgs$method <- "BFGS"
}
if (fit == "optim" &&
!dArgs$method %in% c("BFGS", "CG", "L-BFGS-B") &&
!numGrad) {
numGrad <- TRUE
warning("'numGrad' is set to TRUE since chosen optimization method ",
"is not based on gradient search")
}
stopifnot(is.logical(doFit) && length(doFit) == 1)
stopifnot(is.logical(numGrad) && length(numGrad) == 1)
stopifnot(is.logical(numHess) && length(numHess) == 1)
if (!numHess & numGrad)
stop("'numHess' must be TRUE if numGrad is 'TRUE'")
stopifnot(is.numeric(nGHQ) && length(nGHQ) == 1)
if (nGHQ != round(nGHQ))
warning("'nGHQ' is set to ", nGHQ, ".")
nGHQ <- as.integer(round(nGHQ))
stopifnot(is.null(start) | is.numeric(start))
stopifnot(is.null(restricted) | is.character(restricted))
stopifnot(is.logical(verbose) && length(verbose) == 1)
rval <- append(list(fit = fit,
doFit = doFit,
numGrad = numGrad,
numHess = numHess,
nGHQ = nGHQ,
start = start,
restricted = restricted,
verbose = verbose), dArgs)
class(rval) <- "olmm_control"
return(rval)
}
olmm <- function(formula, data, family = cumulative(),
weights, subset, na.action = na.omit,
offset, contrasts, control = olmm_control(), ...) {
## check arguments
## append '...' arguments to control
cArgs <- list(...)
## cArgs <- cArgs[intersect(names(cArgs), names(formals(olmm_control)))]
cArgsNames <- names(cArgs)
cArgs <- do.call("olmm_control", cArgs)
control[cArgsNames] <- cArgs[cArgsNames]
if (control$verbose) cat("* checking arguments ... ")
mc <- match.call(expand.dots = FALSE)
stopifnot(inherits(formula, "formula"))
stopifnot(inherits(control, "olmm_control"))
## formula
if (length(attr(terms(
as.formula(paste("~", deparse(lhs(formula))))), "term.labels")) > 1L)
stop("the left.hand of the formula must consist of a single term.")
## link and family
if (is.character(family)) {
family <- get(family, mode = "function", envir = parent.frame())
} else if (is.function(family)) {
family <- family()
}
if (!inherits(family, c("family", "family.olmm"))) {
print(family)
stop("'family' not recognized")
}
if (inherits(family, "family")) {
stop("'family' not recognized")
## 2015-10-05: work on gaussian model stopped for a moment
## if (family$family %in% "gaussian") {
## class(family) <- c("family.olmm", "olmm")
## } else {
## stop("'family' not recognized")
## }
## if (!family$link %in% "identity") stop("'link' not recognized")
}
linkNum <- switch(
family$link,
logit = 1L, probit = 2L, cauchy = 3L, identity = 11L)
famNum <- switch(
family$family,
cumulative = 1L, baseline = 2L, adjacent = 3L, gaussian = 11L)
## evaluate contrasts
con <- lapply(1:ncol(data), function(i) attr(data[, i], "contrasts"))
names(con) <- colnames(data)
con <- con[which(!unlist(lapply(con, is.null)))]
if (missing(contrasts)) contrasts <- NULL
contrasts <- appendDefArgs(contrasts, con)
## optimizer control option
optim <- olmm_optim_setup(x = control, env = environment())
## control$numGrad <- control$numHess <- is.null(optim$gr)
## set environment
env <- if (!is.null(list(...)$env))
list(...)$env else parent.frame(n = 1L)
## extract model frames
if (control$verbose) cat("OK\n* extracting model frames ... ")
## decompose model formula
if (any(substr(all.vars(formula), 1, 3) == "Eta"))
stop("'Eta' is a reserved label and cannot be used as ",
"variable name nor as prefix of a variable name.")
formList <- vcrpart_formula(eval.parent(mc$formula),
family = family,
env = env)
if (!is.null(formList$vc)) stop("'vc' terms are not allowed in 'olmm'.")
## set full model frame
m <- match(c("data", "subset", "weights", "na.action"), names(mc), 0L)
mf <- mc[c(1L, m)]
mf$formula <- formList$allTerms
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
fixefmf <- ranefmf <- mf
fullmf <- eval.parent(mf)
## extract responses
y <- model.response(fullmf)
if (famNum < 10) {
if (!is.factor(y)) stop("response must be a 'factor'.")
if (nlevels(y) < 2L)
stop("response variable has less than 2 categories.")
}
## extract fixed effect model matrix
fixefmf$formula <- terms(formList$fe$eta$ce, keep.order = TRUE)
fixefmfCe <- eval.parent(fixefmf)
fixefmf$formula <- terms(formList$fe$eta$ge, keep.order = TRUE)
fixefmfGe <- eval.parent(fixefmf)
conCe <-
contrasts[intersect(names(contrasts), all.vars(formList$fe$eta$ce))]
conGe <-
contrasts[intersect(names(contrasts), all.vars(formList$fe$eta$ge))]
X <- olmm_merge_mm(x = model.matrix(terms(fixefmfCe), fullmf, conCe),
y = model.matrix(terms(fixefmfGe), fullmf, conGe),
TRUE)
rownames(X) <- rownames(fullmf)
X <- olmm_check_mm(X)
## intercept term
if (attr(terms(fixefmfCe), "intercept") == 1L) {
intVar <- colnames(X)[1L]
intTerms <- colnames(X)[1L]
} else {
intVar <- attr(terms(fixefmfCe, keep.order = TRUE), "term.labels")
intVar <- if (length(intVar) > 0L) intVar[1L] else NULL
intTerms <- colnames(X)[attr(X, "assign") == 1L &
attr(X, "merge") == 1L]
if (!is.null(intVar) && !is.factor(fullmf[, intVar])) {
intVar <- intTerms <- NULL
}
}
## extract random effect grouping factor 'subject'
hasRanef <- TRUE
subjectName <- all.vars(formList$re$cond)
if (length(subjectName) == 0L) { # hack to permit models without
## random effects
subjectName <- as.character("id")
formList$re$eta$ge <- as.formula(~ 1)
formList$re$eta$ce <- as.formula(~ -1)
formList$re$eta$cond <- as.formula(~id)
fullmf$id <- factor(1:nrow(fullmf))
control$start["ranefCholFac1"] <- as.numeric(0)
control$restricted <- unique(c(control$restricted, "ranefCholFac1"))
control$nGHQ <- as.integer(1L)
hasRanef <- FALSE
}
subject <- fullmf[, subjectName, drop = TRUE]
if (!is.factor(subject)) stop("subject variable must be a factor")
## extract random effect model matrix W
ranefmf$formula <- terms(formList$re$eta$ce, keep.order = TRUE)
ranefmfCe <- eval.parent(ranefmf)
ranefmf$formula <- terms(formList$re$eta$ge, keep.order = TRUE)
ranefmfGe <- eval.parent(ranefmf)
conCe <- contrasts[intersect(names(contrasts),
all.vars(formList$re$eta$ce))]
conGe <- contrasts[intersect(names(contrasts),
all.vars(formList$re$eta$ge))]
W <- olmm_merge_mm(x = model.matrix(terms(ranefmfCe), fullmf, conCe),
y = model.matrix(terms(ranefmfGe), fullmf, conGe),
FALSE)
rownames(W) <- rownames(fullmf)
W <- olmm_check_mm(W)
## contrasts
cons <- append(attr(X, "contrasts"), attr(W, "contrasts"))
if (!is.null(cons)) cons <- cons[!duplicated(names(cons))]
if (is.null(cons)) storage.mode(cons) <- "list"
## vector for dimensions etc.
nEta <- if (is.factor(y)) nlevels(y) - 1 else 1
dims <- as.integer(c(n = nrow(X), N = nlevels(subject), p = nEta * sum(attr(X, "merge") == 1L) + sum(attr(X, "merge") == 2L), pEta = ncol(X), pInt = length(intTerms), pCe = sum(attr(X, "merge") == 1L), pGe = sum(attr(X, "merge") == 2L), q = nEta * sum(attr(W, "merge") == 1L) + sum(attr(W, "merge") == 2L), qEta = ncol(W), qCe = sum(attr(W, "merge") == 1L), qGe = sum(attr(W, "merge") == 2L), J = nEta + 1, nEta = nEta, nPar = nEta * sum(attr(X, "merge") == 1L) + sum(attr(X, "merge") == 2L) + (nEta * sum(attr(W, "merge") == 1L) + sum(attr(W, "merge") == 2L)) * (1L + nEta * sum(attr(W, "merge") == 1L) + sum(attr(W, "merge") == 2L)) / 2L, nGHQ = control$nGHQ, nQP = control$nGHQ^(nEta * sum(attr(W, "merge") == 1L) + sum(attr(W, "merge") == 2L)), family = famNum, link = linkNum, verb = control$verbose, numGrad = control$numGrad, numHess = control$numHess, doFit = control$doFit, hasRanef = hasRanef))
names(dims) <- c("n", "N", "p", "pEta", "pInt", "pCe", "pGe", "q", "qEta", "qCe", "qGe", "J", "nEta", "nPar", "nGHQ", "nQP", "family", "link", "verb", "numGrad", "numHess", "doFit", "hasRanef")
## parameter names
parNames <- list(fixef = c(if (dims["pCe"] > 0) paste("Eta", rep(seq(1L, dims["nEta"], 1), each = dims["pCe"]), ":", rep.int(colnames(X)[attr(X, "merge") == 1L], dims["nEta"]), sep = ""), if (dims["pGe"] > 0) colnames(X)[attr(X, "merge") == 2L]), ranefCholFac = paste("ranefCholFac", 1L:(dims["q"] * (dims["q"] + 1L) / 2L ), sep = ""))
## set the weights
if (is.null(model.weights(fullmf))) {
weights <- as.double(rep.int(1.0, dims["n"]))
weights_sbj <- as.double(rep.int(1.0, dims["N"]))
} else {
weights <- model.weights(fullmf)
weights_sbj <- tapply(weights, subject, unique)
if (is.list(weights_sbj)) {
stop("'weights' must be constant for subjects")
} else {
weights_sbj <- as.double(weights_sbj)
}
if (length(weights_sbj) != dims["N"]) {
stop("'weights' must be constant for subjects")
}
if (any(weights < 0.0)) stop("negative 'weights' are not allowed")
}
## set the offset
if (missing(offset)) offset <- NULL
if (!is.null(offset) & !is.null(model.offset(fullmf)))
stop("duplicated specification of 'offset'.")
if (is.null(offset)) {
offset <- matrix(0.0, dims["n"], dims["nEta"],
dimnames = list(rownames(fullmf),
paste("Eta", 1L:dims["nEta"], sep = "")))
} else {
if (!is.null(model.offset(fullmf))) offset <- model.offset(fullmf)
if (NCOL(offset) == 1L) offset <- matrix(offset, dims["n"], dims["nEta"])
if (!is.matrix(offset)) stop("'offset' must be a 'matrix'")
if (ncol(offset) != dims["nEta"])
stop("'offset should be a 'matrix' with ", dims["nEta"], " columns")
if (nrow(offset) != nrow(fullmf))
offset <- offset[-attr(fullmf, "na.action"), , drop = FALSE]
if (any(is.na(offset))) stop("'offset' contains NA's.")
if (nrow(offset) != dims["n"]) stop("'offset' has wrong dimensions.")
}
## weights and nodes for the Gauss-Hermite quadrature integration
if (hasRanef) {
gh <- gauss.quad(dims["nGHQ"], "hermite")
ghx <- olmm_expandQP(gh$nodes, dims["q"]) # with correction
ghw <- olmm_expandQP(gh$weights * 1 / sqrt(2 * pi) * exp((gh$nodes^2) / 2),
dims["q"])
} else {
ghx <- matrix(0.0, 1L, 1L)
ghw <- matrix(1.0, 1L, 1L)
}
## elimination matrix for lower triangular matrices
ranefElMat <- L.matrix(n = dims["q"])
## Likelihood function
ll_sbj <- rep.int(0.0, dims["N"])
names(ll_sbj) <- levels(subject)
ll <- c(0.0)
## score function
score_obs <- matrix(0, dims["n"], dims["nPar"])
rownames(score_obs) <- rownames(X)
colnames(score_obs) <- unlist(parNames)
score_sbj <- matrix(0, dims["N"], dims["nPar"],
dimnames = list(levels(subject), unlist(parNames)))
score <- rep.int(0, dims["nPar"])
names(score) <- unlist(parNames)
## info matrix
info <- matrix(0, dims["nPar"], dims["nPar"],
dimnames = list(unlist(parNames), unlist(parNames)))
## linear predictor (without contributions of random effects)
eta <- matrix(0, dims["n"], dims["nEta"],
dimnames = list(rownames(fullmf),
paste("Eta", 1L:dims["nEta"], sep = "")))
## inital values
if (control$verbose) cat("OK\n* setting inital values ... ")
start <- olmm_start(control$start, dims, parNames, X, W, eta, ranefElMat)
## restricted
restr <- rep.int(FALSE, dims["nPar"])
names(restr) <- unlist(parNames)
if (dims["family"] == 3L) control$restricted <- NULL
## check and set 'restr'
if (!is.null(control$restricted)) {
stopifnot(is.character(control$restricted))
if (!all(control$restricted %in% names(restr)))
stop(paste("the coefficient(s) ", paste("'", control$restricted[!control$restricted %in% names(restr)], "'", sep = "", collapse = ", "), " in 'restricted' were not found. The coefficient names are ", paste("'", names(restr), "'", sep = "", collapse = ", "), ".", sep = ""))
restr[control$restricted] <- TRUE
}
## xlevels
xlevels <- .getXlevels(attr(fullmf, "terms"), fullmf)
if (is.null(xlevels)) storage.mode(xlevels) <- "list"
## set the transformed random effect matrix
u <- matrix(0, dims["N"], dims["q"],
dimnames = list(levels(subject),
rownames(start$ranefCholFac)))
## get terms and delete environments
formList <- vcrpart_formula_delEnv(formList)
terms <- list(feCe = terms(formList$fe$eta$ce, keep.order = TRUE),
feGe = terms(formList$fe$eta$ge, keep.order = TRUE),
reCe = terms(formList$re$eta$ce, keep.order = TRUE),
reGe = terms(formList$re$eta$ge, keep.order = TRUE))
environment(formula) <- NULL
attr(attr(fullmf, "terms"), ".Environment") <- NULL
## define fit object
if (control$verbose) cat("OK\n* building the model object ... ")
object <- structure(
list(call = mc,
frame = fullmf,
formula = formula,
terms = terms,
family = family,
y = y,
X = X,
W = W,
subject = subject,
subjectName = subjectName,
weights = weights,
weights_sbj = weights_sbj,
offset = offset,
xlevels = xlevels,
contrasts = cons,
dims = dims,
fixef = start$fixef,
ranefCholFac = start$ranefCholFac,
coefficients = start$coefficients,
restricted = restr,
eta = eta,
u = u,
logLik_sbj = ll_sbj,
logLik = ll,
score_obs = score_obs,
score_sbj = score_sbj,
score = score,
info = info,
ghx = ghx,
ghw = ghw,
ranefElMat = ranefElMat,
control = control,
optim = optim,
output = list(),
converged = FALSE),
class = "olmm")
## delete big data blocks
# rm(list = ls()[!ls() %in% c("dims", "object")]) # commented out this after mail of konstanze.lauseker@wu.ac.at at 2024-05-03
if (dims["doFit"] > 0L) {
## set fitting evironment
if (dims["verb"] > 0L)
cat("OK\n* setting up the fitting environment ... ")
## set start parameters
object$optim[[1L]] <- object$coefficients
object$optim[[4L]] <- object$restricted
## fit the model
if (dims["verb"] > 0L) cat("OK\n* fitting the model ... ")
if (!is.null(object$optim$control$trace) &&
object$optim$control$trace > 0L)
cat("\n")
## extract the function for fitting the model
FUN <- object$optim$fit
subs <- which(names(object$optim) == "fit")
object$optim <- object$optim[-subs]
object$optim$env <- environment()
systemTime <- system.time(
object$output <-
suppressWarnings(do.call(FUN, object$optim)))
object$optim$fit <- FUN
## overwrite slots
new <- .Call("olmm_update_marg", object, object$output$par,
PACKAGE = "vcrpart")
object <- modifyList(object, new)
## print messages for opimization
if (dims["verb"] > 0L) {
cat(paste("OK\n\toptimization time:",
signif(systemTime[3L], 3L),
"seconds", sep = " "))
if (is.null(object$output$message)) {
cat("\n\tno message returned by the optimizer")
} else {
cat(paste("\n\tmessage: ", object$output$message, sep = ""))
}
}
## warnings from optimization
olmm_optim_warnings(object$output, FUN)
object$converged <- switch(
object$optim$fit,
optim = object$output$convergence == 0,
nlminb = object$output$convergence == 0,
ucminf = object$output$convergence %in% c(1, 2, 4))
## numeric estimate of fisher information
if (dims["numHess"] == 1L) {
if (dims["verb"] > 0L)
cat("\n* computing the approximative hessian matrix ... ")
object$info[] <- # replace the info slot
- hessian(
func = object$optim[[2L]],
x = object$coefficients,
method.args = list(func = if (dims["numGrad"]) object$optim[[3L]] else NULL),
restricted = rep.int(FALSE, dims["nPar"]),
env = object$optim$env)
if (dims["verb"] > 0L) cat("OK")
}
if (dims["verb"] > 0L) {
eigenHess <- eigen(object$info, only.values = TRUE)$values
condHess <- abs(max(eigenHess) / min(eigenHess))
cat("\n\tcondition number of Hessian matrix:",
format(condHess, digits = 2L, scientific = TRUE))
}
## fit / predict random effects
## compute expected standardized random effects
if (dims["hasRanef"] > 0) {
if (dims["verb"] > 0L) cat("\n* predicting random effects ... ")
object$u <- .Call("olmm_update_u", object, PACKAGE = "vcrpart")
if (dims["verb"] > 0L) cat("OK")
}
## reset environment of estimation equations
environment(object$optim[[2L]]) <- baseenv()
if (dims["numGrad"] < 1L)
environment(object$optim[[3]]) <- baseenv()
if (dims["verb"] > 0L)
cat("\n* computations finished, return model object\n")
} else {
## update the object with the current estimates
new <- .Call("olmm_update_marg", object, object$coefficients,
PACKAGE = "vcrpart")
object <- modifyList(object, new)
if (dims["hasRanef"] > 0L)
object$u <- .Call("olmm_update_u", object, PACKAGE = "vcrpart")
if (dims["verb"] > 0L)
cat("\n* no computations processed, return model object\n")
}
return(object)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.