At the optimum the covariance is available as vcov(lmer-model). This function
computes cov(beta) at non (RE)ML estimates of varpar.
get_covbeta(varpar, devfun)The covariances matrix of the fixed effects at supplied varpar values.
variance parameters; varpar = c(theta, sigma).
deviance function as a function of theta only.
Rune Haubo B. Christensen. Adapted to pbkrtest by Søren Højsgaard.