# Compute variance-covariance matrix of adjusted means in R lme4

I am trying to fit some mixed models for unbalanced data as follows.

library(easyanova)
data(data13)


### 1) genotypes as random

frmla <- "yield ~ 1 + (1|blocks) + (1|genotypes)"

model1 <- lmer(formula(frmla), data = data13)

# adjusted means - BLUPs for genotypes
newdata13 <- expand.grid(genotypes = levels(data13$genotypes), blocks = levels(data13$blocks))
newdata13$pred <- predict(model1, newdata=newdata13) tapply(newdata13$pred, newdata13$genotypes, mean)  ### 1) genotypes as fixed frmla <- "yield ~ 1 + (1|blocks) + genotypes" model2 <- lmer(formula(frmla), data = data13) # adjusted means - BLUEs for genotypes newdata13 <- expand.grid(genotypes = levels(data13$genotypes), blocks = levels(data13$blocks)) newdata13$pred <- predict(model2, newdata=newdata13)
tapply(newdata13$pred, newdata13$genotypes, mean)


For further calculations I need to compute mean variance of difference of adjusted means (BLUPs or BLUEs). A method is given (https://static-content.springer.com/esm/art%3A10.1186%2F1471-2164-14-860/MediaObjects/12864_2013_5591_MOESM1_ESM.doc) to compute it from variance-covariance matrix of adjusted means.

How to get the variance-covariance matrix of adjusted means for model1 and model2 ?

When genotypes are fixed in model2 does vcov(model2) give the variance-covariance matrix of adjusted means ?