# Confidence intervals in a mixed model

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.

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).

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)[]["position:tube","Mean Sq"] > ag <- aggregate(y~position, data=dat, FUN=mean) > estimate <- ag$y-ag$y > > ### 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))) )  -0.317764 > ( upp.bound <- estimate + sqrt(2*S2/J/K*qf(0.95,1,(J-1)*(I-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
 -0.317764  4.734348
attr(,"conf.level")
 0.95