For simple glm
objects, I can use predict(fit, type = "terms")
to retrieve a matrix with fitted values for each term.
What is the equivalent for lmer
resp. glmer
fitted models? As far as I can see, the predict.merMod
function does not support type = terms
.
lme4: Linear Mixed-Effects Models using 'Eigen' and S4 The models and their components are represented using S4 classes and methods. The core computational algorithms are implemented using the 'Eigen' C++ library for numerical linear algebra and 'RcppEigen' "glue". Version: 1.1-30.
Mixed-model formulas. Like most model-fitting functions in R, lmer takes as its first two arguments a formula spec- ifying the model and the data with which to evaluate the formula. This second argument, data, is optional but recommended and is usually the name of an R data frame.
What is the equivalent for
lmer
resp.glmer
fitted models?
I do not think there is one. Though, you can easily make one as follows
#####
# fit model with one terms which is a matrix
library(lme4)
fit <- lmer(Reaction ~ cbind(Days, (Days > 3) * Days) + (Days | Subject),
sleepstudy)
#####
# very similar code to `predict.lm`
pred_terms_merMod <- function(fit, newdata){
tt <- terms(fit)
beta <- fixef(fit)
mm <- model.matrix(tt, newdata)
aa <- attr(mm, "assign")
ll <- attr(tt, "term.labels")
hasintercept <- attr(tt, "intercept") > 0L
if (hasintercept)
ll <- c("(Intercept)", ll)
aaa <- factor(aa, labels = ll)
asgn <- split(order(aa), aaa)
if (hasintercept) {
asgn$"(Intercept)" <- NULL
avx <- colMeans(mm)
termsconst <- sum(avx * beta)
}
nterms <- length(asgn)
if (nterms > 0) {
predictor <- matrix(ncol = nterms, nrow = NROW(mm))
dimnames(predictor) <- list(rownames(mm), names(asgn))
if (hasintercept)
mm <- sweep(mm, 2L, avx, check.margin = FALSE)
for (i in seq.int(1L, nterms, length.out = nterms)) {
idx <- asgn[[i]]
predictor[, i] <- mm[, idx, drop = FALSE] %*% beta[idx]
}
} else {
predictor <- ip <- matrix(0, n, 0L)
}
attr(predictor, "constant") <- if (hasintercept) termsconst else 0
predictor
}
# use the function
newdata <- data.frame(Days = c(1, 5), Reaction = c(0, 0))
(out <- pred_terms_merMod(fit, newdata))
#R> cbind(Days, (Days > 3) * Days)
#R> 1 -21.173
#R> 2 21.173
#R> attr(,"constant")
#R> [1] 283.24
#####
# confirm results
beta. <- fixef(fit)
beta.[1] + beta.[2]
#R> (Intercept)
#R> 262.07
out[1] + attr(out, "constant")
#R> [1] 262.07
beta.[1] + (beta.[2] + beta.[3]) * 5
#R> (Intercept)
#R> 304.41
out[2] + attr(out, "constant")
#R> [1] 304.41
Extending the above to also include standard errors should be straightforward as far as I gather.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With