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
