2
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.