---
title: Reducing a model to get confidence intervals
date : 2014-04-12
--- &lead
```{r setup0, echo=FALSE}
opts_chunk$set(fig.path="assets/fig/ModelReduction", tidy=FALSE)
```
```{r options, echo=FALSE, message=FALSE}
library(pander)
```
## Reducing the one-way ANOVA model with random effects
Consider the balanced one-way ANOVA model with random effects, call $\cal M$
this model:
$$({\cal M})\colon \qquad \begin{cases}
(y_{ij} \mid \mu_i) \sim_{\text{iid}} {\cal N}(\mu_i, \sigma^2_w) & j=1,\ldots,J \\
\mu_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2_b) & i=1,\ldots,I
\end{cases}.$$
Let us derive the distribution of the observed group means $\bar{y}_{i\bullet}$:
$$(\bar{y}_{i\bullet} \mid \mu_i) \sim {\cal N}\left(\mu_i, \frac{\sigma^2_w}{J}\right)$$
for every $i=1,\ldots,I$, and finally
$$\bar{y}_{i\bullet}\sim_{\text{iid}} {\cal N}\left(\mu, \sigma^2_b+\frac{\sigma^2_w}{J}\right) \qquad i=1,\ldots,I.$$
Denote $z_i=\bar{y}_{i\bullet}$ for more clarity. Thus ${(z_i)}_{i=1}^I$ is a sample
of *iid* Gaussian distributions:
$$({\cal M}')\colon \qquad
z_i \sim_{\text{iid}} {\cal N}\left(\mu, \delta^2\right), \quad i=1,\ldots,I$$
with $\delta^2=\sigma^2_b+\frac{\sigma^2_w}{J}$.
Thus, we have a new model ${\cal M}'$ derived from $\cal M$ and for which there is
a well-known $100(1-\alpha)\%$-confidence interval around $\mu$:
$$\bar{z}_{\bullet} \pm t^*_{I-1}(\alpha/2)\frac{SS(z)}{(I-1)I}$$
where $SS(z)=\sum_{i=1}^I{(z_i-\bar{z}_{\bullet})}^2$ is the total variation
of the $z_i$.
Compare this confidence interval to the one derived at
the end of my article about
[the balanced one-way ANOVA model with random effects](http://stla.github.io/stlapblog/posts/Anova1random.html): *they are exactly the same !*
Thus, there's no lost of information about $\mu$ when we observe only
the group means $\bar{y}_{i\bullet}$.
## General principle of model reduction
Was the previous remark an expected one ?
Given the original model $\cal M$ with observations $y=(y_{ij})$, the reduced
model ${\cal M}'$ has been derived by considering the $z_i:=f_i(y)$ as the
observations for some functions $f_i$.
The parameter $\mu$ appears in both models.
The obvious thing is that any $100(1-\alpha)\%$-confidence interval
about $\mu$ in the reduced model ${\cal M}'$ also is a
$100(1-\alpha)\%$-confidence interval
about $\mu$ in the original model $\cal M$.
But generally such a method will not provide an *optimal* confidence interval
(in a sense to be made precise), as shown by the following other example
## Reducing the two-ways ANOVA model without interaction
Consider the fixed effects two-ways ANOVA model
with replications but without interaction :
$$({\cal M})\colon \qquad
y_{ijk} \sim {\cal N}(\mu+\alpha_i+\beta_j, \sigma^2),
\qquad i=1,\ldots,I, \quad j=1,\ldots,J, \quad k=1,\ldots,K,$$
with $\sum \alpha_i=\sum \beta_j=0$ and
assuming as usual independent observations $y_{ijk}$.
Now consider the reduced model obtained by averaging the observations
in the same "cell" $(i,j)$:
$$({\cal M}')\colon \qquad
\bar y_{ij\bullet} \sim {\cal N}\left(\mu+\alpha_i+\beta_j, \frac{\sigma^2}{K}\right), \qquad i=1,\ldots,I, \quad j=1,\ldots,J.$$
The reduced model ${\cal M}'$ is a fixed effects two-ways ANOVA model
without interaction and without replications.
We will consider for our illustration that $I=2$, hence $\alpha_2=-\alpha_1$,
and we will compare the confidence intervals about $\alpha_1$ obtained by
considering either model $\cal M$ or model ${\cal M}'$.
Let us simulate some data
```{r}
I <- 2
J <- 3
K <- 4
effects1 <- c(1,-1)
effects2 <- 1:J - mean(1:J)
dat <- data.frame(factor1=rep(1:I,each=J*K), factor2=rep(1:J,times=I*K))
set.seed(666)
dat <- transform(dat,
y=rnorm(I*J*K, effects1[factor1]+effects2[factor2]),
factor1=factor(factor1),
factor2=factor(factor2)
)
head(dat)
```
Fitting the model as follows, the parameter of interest $\alpha_1=1$ is named `factor11` in R:
```{r, results='asis'}
fit1 <- lm(y~factor1+factor2,data=dat,contrasts=list(factor1="contr.sum",factor2="contr.sum"))
pandoc.table(confint(fit1), style="rmarkdown", digits=3, emphasize.rows=2)
```
Now we look at the confidence interval in the reducel model ${\cal M}'$:
```{r, results='asis'}
mdat <- aggregate(y~factor1+factor2, data=dat, FUN=mean)
fit2 <- lm(y~factor1+factor2,data=mdat,contrasts=list(factor1="contr.sum",factor2="contr.sum"))
pandoc.table(confint(fit2), style="rmarkdown", digits=3, emphasize.rows=2)
```
The confidence interval in ${\cal M}'$ is wider than the one in $\cal M$.
Let us use simulations to see this more precisely.
Below we use simulations to derive the distributions of the upper bound
of the confidence interval in $\cal M$ and ${\cal M}'$.
We use the `local` function to preserve our preliminary defined
objects `dat`, `fit1`, and `fit2`.
We also evaluate, for each model the power of the test for
$H_0\colon\{\alpha_1=0\}$ derived from the confidence interval.
```{r, fig.width=7, fig.height=5}
nsims <- 5000
confint1 <- confint2 <- NULL
power1 <- power2 <- rep(NA,nsims)
local({
for(i in 1:nsims){
# new data
dat <- within(dat, y <- rnorm(I*J*K, effects1[factor1]+effects2[factor2]))
# model M
fit1 <- lm(y~factor1+factor2,data=dat,contrasts=list(factor1="contr.sum",factor2="contr.sum"))
ci <- confint(fit1)[2,]
power1[i] <<- ci[1]>0 || ci[2]<0
confint1 <<- rbind(confint1,ci)
# model M'
mdat <- aggregate(y~factor1+factor2, data=dat, FUN=mean)
fit2 <- lm(y~factor1+factor2,data=mdat,contrasts=list(factor1="contr.sum",factor2="contr.sum"))
ci <- confint(fit2)[2,]
power2[i] <<- ci[1]>0 || ci[2]<0
confint2 <<- rbind(confint2,ci)
}
})
plot(density(confint1[,2]),
main="Distribution of the upper bound of the confidence interval",
ylab=NA, xlab=NA,
axes=FALSE, col="blue",
xlim=range(confint2[,2]))
lines(density(confint2[,2]), col="red")
alpha1 <- effects1[1]
axis(1, at=alpha1, labels=expression(alpha[1]))
abline(v=alpha1, lty="dashed")
```
As we said in the previous section, the confidence interval derived from
the reduced model ${\cal M}'$ is correct: for both curves, the area at the left
of $\alpha_1$ is $2.5\%$.
But the upper bound is considerably more dispersed for model ${\cal M}'$,
thereby showing that the confidence interval is more likely wider
for model ${\cal M}'$.
As a consequence, the power of the test of $H_0\colon\{\alpha_1=0\}$ is higher
for $\cal M$:
```{r}
mean(power1)
mean(power2)
```
Note that $H_0$ means that the first factor has no effect, and actually it can
be shown that the test based on the confidence interval is exactly the same
as the classical Fisher test provided by the `anova` function:
```{r}
anova(fit1)
anova(fit2)
```
We can understand why the power is suboptimal in ${\cal M}'$.
As you know, the $F$ statistic is the ratio of the mean square of
factor $1$ over the residual mean square.
The residual sum of squares in the
decomposition of the variation provided by `anova(fit1)`:
$$
SS(y) = SS_A(y) + SS_B(y) + SS_{ {(1+A+B)}^\perp}(y)
$$
can be itself decomposed by including the interaction term, thereby
yielding a new decomposition :
$$
SS(y) = SS_A(y) + SS_B(y) + SS_{AB} + SS_{ {(1+A+B+AB)}^\perp}(y)
$$
which is here:
```{r}
anova(update(fit1, y~factor1*factor2))
```
The first three mean squares are exactly $K=4$ times the sum of squares in
${\cal M}'$:
```{r}
K*anova(fit2)[,2:3]
```
Thus, we have the following equality about the classical Fisher statistic $F'$
for $H_0\colon\{\alpha_1=0\}$ in the reduced model ${\cal M}'$:
$$
F' := \frac{MS_A(z)}{MS_{ {(1+A+B)}^\perp}(z)}
= \frac{MS_A(y)}{MS_{AB}(y)}
\sim F_{I-1, (I-1)(J-1)}
$$
whereas the classical Fisher statistic $F$
for $H_0\colon\{\alpha_1=0\}$ in the unreduced model ${\cal M}$ is:
$$
F = \frac{MS_A(y)}{MS_{ {(1+A+B)}^\perp}(y)}
\sim F_{I-1, (I-1)(J-1) + IJ(K-1)}
$$
thereby showing that the reducel model ${\cal M}'$ yields a less
powerful test.
## Why does it work for the random effects one-way ANOVA model ?
Using the hierarchical writing of the model
$$\begin{cases}
(y_{ij} \mid \mu_i) \sim_{\text{iid}} {\cal N}(\mu_i, \sigma^2_w) & j=1,\ldots,J \\
\mu_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2_b) & i=1,\ldots,I
\end{cases},$$
it is quite easy to write down the joint distribution of the $y_{ij}$ and the $\mu_i$.
Setting $\sigma^2=J\sigma^2_b+\sigma^2_w$ and integrating over the $\mu_i$ yields the likelihood
$$
L\bigl(\mu, \sigma_b, \sigma_w \mid \{y_{ij}\}\bigr) \propto
{(\sigma_w^2)}^{-\frac{I(J-1)}{2}}{(\sigma^2)}^{-\frac{I}{2}}
\exp\left[-\frac{1}{2}\left(
\frac{IJ{(\bar{y}_{\bullet\bullet}-\mu)}^2}{\sigma^2} + \frac{SS_b(y)}{\sigma^2} + \frac{SS_w(y)}{\sigma^2_w}
\right)\right],
$$
thereby showing that the triple of the three summary statistics
$\bigl\{\bar{y}_{\bullet\bullet}, SS_b(y), SS_w(y)\bigr\}$ is sufficient
for $(\mu, \sigma_b, \sigma_w)$.
Moreover, because of the equality
$$IJ{(\bar{y}_{\bullet\bullet}-\mu)}^2 + SS_b(y)
= J\sum_{i=1}^I {(\bar{y}_{i\bullet}-\mu)}^2,$$
the likelihood can also be written
$$
\begin{align}
L\bigl(\mu, \sigma_b, \sigma_w \mid \{y_{ij}\}\bigr) & \propto
{(\sigma_w^2)}^{-\frac{I(J-1)}{2}}{(\sigma^2)}^{-\frac{I}{2}}
\exp\left[-\frac{1}{2}\left(
\frac{J\sum_{i=1}^I {(\bar{y}_{i\bullet}-\mu)}^2}{\sigma^2} + \frac{SS_w(y)}{\sigma^2_w}
\right)\right] \\
& \propto
\left({(\sigma_w^2)}^{-\frac{I(J-1)}{2}}
\exp\left[-\frac{1}{2}
\frac{SS_w(y)}{\sigma^2_w}
\right]\right)
\left(
{(\sigma^2)}^{-\frac{I}{2}}
\exp\left[-\frac{1}{2}
\frac{J\sum_{i=1}^I {(\bar{y}_{i\bullet}-\mu)}^2}{\sigma^2}
\right]
\right).
\end{align}
$$
That shows that the group means $\bar{y}_{i\bullet}$ are sufficient for
$(\mu, \sigma)$.