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