1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.