With R, I want to get the multivariate normal marginal (unconditional) distribution of the response of a mixed model as a function of the parameters (fixed effects and variance components). For example for a one-way ANOVA model with random effects I want the variance-covariance matrix
$$
\begin{pmatrix}
\sigma^2_b + \sigma^2_w & \sigma^2_b & \sigma^2_b & \cdots \\
\sigma^2_b & \sigma^2_b + \sigma^2_w & \sigma^2_b & \cdots \\
\sigma^2_b & \sigma^2_b & \ddots & \ddots \\
\vdots & \vdots & \ddots & \ddots
\end{pmatrix}.
$$
I believe I already succeed to get it with lmer
and getME
but I lost my code. Moreover, if possible, I would like to get it without having to fit the model with lme4
.
1 Answer
I get it, I think. For the mean, just do model.matrix(fixed) %*% beta
. For the variance-covariance matrix:
data <- data.frame(
group = gl(2L, 4L),
treatment = gl(2L, 2L, 8L)
)
fixed <- ~ 1
random <- ~ group * treatment
library(lazyeval)
data <- droplevels(data)
tf <- terms.formula(random)
factors <- rownames(attr(tf, "factors"))
tvars <- attr(tf, "variables")
tlabs <- attr(tf, "term.labels")
tvars <- eval(tvars, envir = data)
rdat <- lapply(tvars, function(tvar) droplevels(as.factor(tvar)))
names(rdat) <- factors
RE <- lapply(tlabs, function(tlab){
droplevels(lazy_eval(as.lazy(tlab), data = rdat))
})
# just to show:
setNames(as.data.frame(RE), tlabs)
# group treatment group:treatment
# 1 1 1 1:1
# 2 1 2 1:2
# 3 1 3 1:3
# 4 2 1 2:1
# 5 2 2 2:2
# 6 2 3 2:3
n <- nrow(data)
RE2 <- c(RE, list(error = factor(1L:n))) # Adds the error effect
E <- c(vapply(RE, nlevels, integer(1L)), n)
Z <- NULL
for(i in seq_along(E)){ # Builds an indicator matrix for the effects
re_levels <- levels(RE2[[i]])
for(j in 1L:E[i]){
temp1 <- which(RE2[[i]] == re_levels[j])
temp2 <- integer(n)
temp2[temp1] <- 1L
Z <- cbind(Z, temp2)
}
}
variances <- c(1, 2, 3, 4)
Z %*% diag(rep(variances, times = E)) %*% t(Z)
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,] 10 6 1 1 2 2 0 0
# [2,] 6 10 1 1 2 2 0 0
# [3,] 1 1 10 6 0 0 2 2
# [4,] 1 1 6 10 0 0 2 2
# [5,] 2 2 0 0 10 6 1 1
# [6,] 2 2 0 0 6 10 1 1
# [7,] 0 0 2 2 1 1 10 6
# [8,] 0 0 2 2 1 1 6 10