Confidence intervals in a mixed model

Question

Consider the following mixed model: $$(y_{ijk} \mid \mu_{ij}) \sim_{\text{i.i.d}} {\cal N}(\mu_{ij}, \sigma^2), \quad k=1, \ldots K $$ $$ \begin{pmatrix} \mu_{i1} \\ \vdots \\ \mu_{iJ} \end{pmatrix} \sim_{\text{i.i.d}} {\cal N}\left(\begin{pmatrix} \mu_{1} \\ \vdots \\ \mu_{J} \end{pmatrix}, \Sigma\right), \quad i=1, \ldots, I $$ with unknown fixed parameters $\sigma^2$, $\Sigma$, $\mu_i$. Moreover I assume $\Sigma$ has a "compound symmetry" structure. Therefore the above model is marginally equivalent to the 2-way ANOVA model with mixed effects $$y_{ijk}=\mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ijk}$$ with $\beta_j$ and $(\alpha\beta)_{ij}$ the random effects. I want some confidence intervals on the fixed parameters $\alpha_i$. I'm using PROC MIXED in SAS but that does not work for some datasets (that does not work with lme/lmer in R too). Do you know some confidence intervals based on least-squares instead of likelihood ? I am interested in the formulas giving these intervals and/or the way to calculate them with a software.

Answer 1

I have just found an answer for the ANOVA model in the book "Design and Analysis of Gauge R&R Studies" (Richard K. Burdick, Connie M. Borror, Douglas C. Montgomery).

Answer 2

Very amazing !!!! The confidence interval is exactly the same as the one obtained with the paired t-test on the means :

> library(mvtnorm)
> 
> ### SIMULATES DATA ###
> I <- 2 # positions  
> J <- 6 # tubes  
> K <- 5 # repeats 
> n <- I*J*K
> tube <- rep(1:J, each=I)
> position <- rep(LETTERS[1:I], times=J)
> Mu_i <- 3*(1:I)
> tube <- rep(tube, each=K)
> position <- rep(position, each=K)
> sigmaw <- 2 
> Mu_ij <- c(t(rmvnorm(J, mean=Mu_i, sigma=diag(I)+2)) )  
> Mu_ij <- rep(Mu_ij, each=K)
> dat <- data.frame(tube, position)
> dat$y <- rnorm(n, Mu_ij, sigmaw)
> dat$tube <- factor(dat$tube)
> an <- aov(y~position*tube, data=dat)
> S2 <- summary(an)[[1]]["position:tube","Mean Sq"]
> ag <- aggregate(y~position, data=dat, FUN=mean)
> estimate <- ag$y[2]-ag$y[1]
> 
> ### 95% CONFIDENCE INTERVAL BASED ON THE MIXED ANOVA MODEL ###
> ( low.bound <- estimate - sqrt(2*S2/J/K*qf(0.95,1,(J-1)*(I-1))) )
[1] -0.317764
> ( upp.bound <- estimate + sqrt(2*S2/J/K*qf(0.95,1,(J-1)*(I-1))) )
[1] 4.734348
> 
> ### 95% CONFIDENCE INTERVAL BASED ON THE PAIRED T-TEST ###
> ag <- aggregate(y~position+tube, data=dat, FUN=mean)
> posA <- subset(ag, subset= position=="A")$y
> posB <- subset(ag, subset= position=="B")$y
> t.test(x=posB, y=posA, paired=TRUE)$conf.int
[1] -0.317764  4.734348
attr(,"conf.level")
[1] 0.95