Lme4 calculates covariance confidence intervals

Question

Please see the answer from Ben Bolker 05/16/2016 for an appropriate solution. OP below.

---

I install several tiered models with lme4. I would like to report deviations and covariances of random effects and automate this process.

I know that I can get deviations with as.data.frame(VarCorr(mymodel)) , and I know that I can get confidence intervals using confint(mymodel) . It’s clear that I can join / relate the two tables and put the confidence intervals around the variances by simply squaring the confint() output square in the corresponding rows and columns, but I could not find a convincing method for calculating covariances, if not manually.

Let's say the result is confint :

 conf <- NULL a <- c(6.2,-0.4,2.2,1.5,-0.4,-0.5,2.8,-0.9,1.3,3.9) b <- c(6.8,-0.2,2.5,2.5,0.1,0.2,4.8,-0.7,2.3,5) conf <- data.frame(a,b,row.names = c("sd_(Intercept)|ID","cor_Time.(Intercept)|ID","sd_Time|ID","sd_(Intercept)|Group","cor_Time.(Intercept)|Group","cor_I(Time^2).(Intercept)|Group","sd_Time|Group","cor_I(Time^2).Time|Group","sd_I(Time^2)|Group","sigma")) colnames(conf) <- c("2.5%","97.5%") conf 

How can I automate various multiplications to obtain covariances such as

 cov.time.intercept <- conf[1,2]*conf[1,1]*conf[1,3] 

?

I tried to separate the standard deviations and correlations, creating the variables "ID", "Time", "I (Time ^ 2)" and "(Intercept)", and then juxtaposing two columns, but I will not go anywhere. The problem is that each time the model changes, you may have a different number of variances and covariances and different triangular matrices.

Thanks for any help

to.

Answer 1

Your calculations seem to give plausible answers, but it doesn't make sense (for me, I'm ready to be corrected / enlightened ...). Suppose cov = corr*var1*var2 . Assume that ci(.) Is the (lower or upper) confidence limit for the quantity. It is by no means true that ci(cov) = ci(corr)*ci(var1)*ci(var2) (it is interesting that you get reasonable answers, I think this is most likely to happen when the amount is approximately uncorrelated. ..) If you had deviations from each component and the covariance among them (I do not mean the variance of random effects and the covariances themselves, but their sampling / covariance sampling), you can distribute them approximately using the delta method, but there are enough of them hard to get (see here ).

The “right” way to do this, as far as I can tell, is to calculate the likelihood profile on the variance-covariance scale instead of the standard deviation scale. It was not possible, but now (with the development version on Github).

Install the latest version:

 library(ghit) ## for install_github (or library(devtools)) install_github("lme4/lme4") 

Qualifying:

 chap12 <- foreign::read.dta(file = "ch12.dta") library(lme4) snijders <- lmer(prox_pup ~ 1 + prox_sel + (1 + occ|teacher), data = chap12) as.data.frame(VarCorr(snijders)) ## grp var1 var2 vcov sdcor ## 1 teacher (Intercept) <NA> 0.15617962 0.3951957 ## 2 teacher occ <NA> 0.01205317 0.1097869 ## 3 teacher (Intercept) occ -0.03883458 -0.8950676 ## 4 Residual <NA> <NA> 0.04979762 0.2231538 

We have to be a little careful when comparing the results, because the profile.merMod , which we will use in the near future, will automatically (and silently!) Convert the settings from REML from default to maximum likelihood (since a profile based on REML can be statistically risky ); however, it does not look like it makes a huge difference.

 s2 <- refitML(snijders) as.data.frame(VarCorr(s2)) ## grp var1 var2 vcov sdcor ## 1 teacher (Intercept) <NA> 0.15426049 0.3927601 ## 2 teacher occ <NA> 0.01202631 0.1096645 ## 3 teacher (Intercept) occ -0.03884427 -0.9018483 ## 4 Residual <NA> <NA> 0.04955549 0.2226106 p.sd <- profile(s2,which="theta_", signames=FALSE) p.vcov <- profile(s2,which="theta_",prof.scale="varcov", signames=FALSE) 

We get some warnings about nonmonotonic profiles ...

 confint(p.vcov) ## 2.5 % 97.5 % ## var_(Intercept)|teacher 0.08888931 0.26131067 ## cov_occ.(Intercept)|teacher -0.07553263 -0.01589043 ## var_occ|teacher 0.00000000 0.02783863 ## sigma 0.03463184 0.07258777 

What if we check the square of the corresponding (sd / variance) elements?

 confint(p.sd)[c(1,3,4),]^2 ## 2.5 % 97.5 % ## sd_(Intercept)|teacher 0.089089363 0.26130970 ## sd_occ|teacher 0.002467408 0.02779329 ## sigma 0.034631759 0.07263869 

They combine well, with the exception of the lower boundary of the variance occ ; they also match your results above. However, the result of the covariance (which, it seems to me, is difficult) gives me (-0.0755, -0.0159) for me, against (-0.0588, -0.0148) for you, about a 20% difference. This may not be very important, depending on what you are trying to do.

Try brute force as well:

 sumfun <- function(x) { vv <- as.data.frame(VarCorr(x),order="lower.tri")[,"vcov"] ## cheating a bit here, using internal lme4 naming functions ... return(setNames(vv, c(lme4:::tnames(x,old=FALSE,prefix=c("var","cov")), "sigmasq"))) } cc <- confint(s2,method="boot",nsim=1000,FUN=sumfun,seed=101, .progress="txt", PBargs=list(style=3)) ## .progress/PBargs just cosmetic ... ## 2.5 % 97.5 % ## var_(Intercept)|teacher 0.079429623 0.24053633 ## cov_occ.(Intercept)|teacher -0.067063911 -0.01479572 ## var_occ|teacher 0.002733402 0.02378310 ## sigmasq 0.031952508 0.06736664 

The "gold standard" here seems to be partially between the result of my profile and your result: the lower limit of the covariance is -0.067 here versus -0.0755 (profile) or -0.0588.

Answer 2

Solved, thanks for the contribution. I will update the initial entry. The result can be checked using the dataset from Snijders and Bosker, available here .

Import from

 library(foreign) chap12 <- read.dta(file = "<your path>/ch12.dta") 

Improvised Model:

 snijders <- lmer(prox_pup ~ 1 + prox_sel + (1 + occ|teacher), data = chap12) 

Enter function:

 ExtractVarCovCI <- function(Model) { v <- NULL v <- as.data.frame(VarCorr(Model),order = "lower.tri") #Extract variances and covariances conf <- confint(Model, parm ="theta_", oldNames = F) #extract CIs v.conf <- cbind(v,conf) #bind confidence intervals covs <- as.data.frame(v.conf[!is.na(v[,3]),]) #separate variance from covariance components vars <- as.data.frame(v.conf[is.na(v[,3]),]) #separate variance from covariance components vars.sq <- vars[,6:7]^2 #calculate square of variance components colnames(vars.sq) <- sub("[%]", "% sq.", colnames(vars.sq)) vars2 <- cbind(vars,vars.sq) #bind squares of variance components covs$`2.5 % sq.` <- c(rep(NA,nrow(covs))) #create empty columns for later covs$`97.5 % sq.` <- c(rep(NA,nrow(covs))) #create empty columns for later lcovs <- length(row.names(covs)) #now we re-organise the table so that each covariance is below the variance of its variables k <- NULL for (i in seq(1:lcovs)) { k <- rbind(k,vars2[vars2$grp %in% covs[i,1] & vars2$var1 %in% covs[i,2],],vars2[vars2$grp %in% covs[i,1] & vars2$var1 %in% covs[i,3],],covs[i,]) } k2 <- rbind(k,vars2["sigma",]) #bind the level-1 residuals at the end k2.covrow <- grep("^cor",rownames(k2)) # isolate covariance row position k2[k2.covrow,8] <- k2[k2.covrow,6]*k2[k2.covrow-1,6]*k2[k2.covrow-2,6] #calculate covariance 2.5% k2[k2.covrow,9] <- k2[k2.covrow,7]*k2[k2.covrow-1,7]*k2[k2.covrow-2,7] #calculate covariance 97.5% p <- NULL p <- k2[,c(4,8:9)] #retain only the estimates and the confidence intervals rownames(p) <- sub("^sd","var",rownames(p)) #now it clear that we have proper variances and covariances rownames(p) <- sub("^cor","cov",rownames(p)) #now it clear that we have proper variances and covariances colnames(p) <- c("Estimate", "2.5%", "97.5%") return(p) } 

Run the function:

 ExtractVarCovCI(snijders) 

My result:

  Estimate 2.5% 97.5% var_(Intercept)|teacher 0.15617962 0.089020350 0.26130969 var_occ|teacher 0.01205317 0.002467408 0.02779329 cov_occ.(Intercept)|teacher -0.03883458 -0.014820577 -0.05887660 sigma 0.04979762 0.034631759 0.07263837 

Now we have a variance-covariance table that uses non-standardized random effects with their upper and lower confidence limits. I am sure there is a better way to do this, but this is the beginning ...

to.

Answer 3

Note that the standard deviation of random effects in lme4 NOT a standard variance error! This is just the square root of the variance!

If you need confidence intervals for a variance of random effects, you need profile() probability. See ?lme4::profile .