# How to calculate confidence intervals for linear mixed effects models when default methods with default settings fail?

I have a simple linear model describing a set of straight lines and would like to estimate confidence intervals for the parameters and the covariance matrix describing the hidden parameters. First attempts using linear mixed effects models as implemented in R using default options failed.

The confidence intervals coud be calculated but the coverage of the 90% confidence intervals was, in general, too small. Instead of having 10% of the confidence intervals that did not include the correct value of the parameters, about 20% did not include the correct values. This was despite quite lengthy calculation times that were required to determine the confidence intervals.

How can one improve upon this???

The remainder of the question desribe the model and attempts at getting the confidence intervals using bootstrap.

### The model

A model is considered with response variable $y_{i,j}$. The model has three parameters: $\theta_{0}$, an offset, $\theta_{1}$, a slope, and $\theta_{2}$ a factor for the explanatory variable $X_i$. There are two hidden variables $\epsilon_{0,i}$ and $\epsilon_{1,i}$ on the offset and the slope, and there is a residual error, $\epsilon_{2,i,j}$. The model can be written as

$$y_{i,j} = \theta_{0} + \theta_{1} * t_j + \theta_{2} * x_i + \epsilon_{0,i} + \epsilon_{1,i} * t_j + \epsilon_{2,i,j}$$

with the residual error defined by

$$\epsilon_{2,i,j} \sim \mathrm{Norm}(0,vr)$$

and the hidden parameters described by

$$\left( \begin{array}{c} \epsilon_{0,i}\\ \epsilon_{1,i}\\ \end{array} \right) \sim \mathrm{Norm}(0,\mathbf{V})$$.

The aim is to estimate with confidence intervals the parameters, $\theta_{0}$, $\theta_{1}$, $\theta_{2}$, the variance of the resifual error, $vr$, and the 2-dimensional covariance matrix, $\mathbf{V}$, of the hidden effects.

### Example

The model was evaluated setting the parameters to 1, the covariance matrix to a unit matrix, and the standard deviation of the residual error to 0.1. The explanatory variable was set to $x = \{2,4,6,8,10\}$, and the time variable to $t = \{1,2,3,4,5\}$.
Data was simulated repeatedly and then analyzed using the lmer and confint commands of the lme4 package with bootstrap to determine the confidence intervals. The commands used were

m1 <- lmer(V3~1+V1+V4+(1+V1|V2),data=dat2)

ci <- confint(m1,method="boot",level=0.9)


Full code is provided below. Using these settings the confidence intervals coud be calculated but the coverage was, in general, too small. Instead of having 10% of the confidence intervals that did not include the correct value of the parameters, about 20% did not include the correct values. This was despite quite lengthy calculation times that were required to determine the confidence intervals.

How can one improve upon this???

### Full code

require(lme4)

set.seed(123)

# set 'true' values of parameters
vp<-c(.sig01=1,.sig02=0,.sig03=1,.sigma=0.1,"(Intercept)"=0,V1=0,V4=1)

# setup counts to evaluate error rates
cnts<-rep(0,length(vp))

# simulate data, determine CI, compare to 'true' values
ni<-200
it<-0
for(i in 1:ni) {
cat(i)
ns<-5 # number of subjects
np<-3 # parameters per subject
subj <- matrix(NA,ns,np)
subj[,1] <- rnorm(ns,vp["(Intercept)"],vp[".sig01"])
subj[,2] <- rnorm(ns,vp["V1"],vp[".sig03"])
subj[,3] <- (1:5)*2

dx<-1:5
nm<-length(dx)

dat1<-matrix(NA,ns*nm,4)
dat1[,1]<-dx
dat1[,2]<-sort(rep(1:ns,nm))
dat1[,4]<-subj[dat1[,2],3]
dat1[,3]<-subj[dat1[,2],1] + dat1[,1]*subj[dat1[,2],2] + vp["V4"]*dat1[,4] + rnorm(ns*nm,0,vp[".sigma"])

dat2 <- as.data.frame(dat1)

try({
m1 <- lmer(V3~1+V1+V4+(1+V1|V2),data=dat2)
ci <- confint(m1,method="boot",level=0.9)
cnts<-cnts+(ci[,1]<=vp&ci[,2]>=vp)
it <- it+1
})

}


### Results of executing the code

print(cnts/it)


Except for the residual error .sigma, the coverage is too low, with some of the parameters, e.g., V4, having a coverage of 70% instead of 90%.

##      .sig01      .sig02      .sig03      .sigma (Intercept)          V1
##   0.8131313   0.7424242   0.7929293   0.8939394   0.7323232   0.8383838
##          V4
##   0.7070707


This seems to be a challenging question. Having asked it some time ago, it remained unanswered as tumbleweed and got finally removed. The manuscript cited below makes progress on the question, the proposed method could deal with the question; an efficient implementation is still missing.

Key elements are

• An efficient integration method for repeated calculation of statistical integrals for a set of hypotheses (e.g., p-values, confidence intervals) using importance sampling
• Introduction of pointwise mutual information as an efficient test statistics that was shown to be optimal under certain conditions.
• Eliminating the need for complex minimax optimizations by demonstrating that priors may be used to derive loss functions, rather than trying to determine optimal minimax priors given the loss function.

The proposed approach has a few fringe benefits:

• The proposed test statistics (point wise mutual information) can be used for frequentist and Bayesian inference and results in consistent confidence intervals and credible intervals.
• The proposed approach is entirely based on the Likelihood function
• It includes a generic proposal on how to use prior information for frequentist tests.

Still missing is an efficient implementation of marginalization.

Ref: Using prior knowledge in frequentist tests. figshare. https://doi.org/10.6084/m9.figshare.4819597.v3