# How to describe variation on multiple levels?

This is continuation for a series of questions (1, 2). I have a data set from an experiment with 2x2 design. A replicate consists of 20-50 normally distributed replicate measurements. Each treatment combination has 3 of these replicates. Measured average response for each replicate seems not to be normally distributed.

Figure 1. Example of the problem. Error bars show standard error of mean.

Since this is an experiment and I am supposed to talk about an effect of a treatment, I would like to combine the variation within measurements (error bars) and between replicates (dots in the figure) into one number / treatment combination and use the error bars with one mean value to show the overall effect of a treatment (a dot instead of three). Ideally I would like to use confidence intervals of some sort (asked here). I have understood that since this is an unbalanced nested design, I should use mixed models or a simple 2-way ANOVA for mean values, but how to compress the variation on two levels into one figure?

Would the delta method be something for me?

R code with an example of the data:

library(ggplot2)

x1 <- c(rnorm(50, 14,3),rnorm(35, 7,1),rnorm(40, 15,9))
x2 <- c(rnorm(43, 6,3),rnorm(32, 7,1),rnorm(40, 8,4))
x3 <- c(rnorm(50, 15,5), rnorm(50, 10,7), rnorm(50, 13,9))
x4 <- c(rnorm(26, 14,2), rnorm(43, 25,10), rnorm(45, 15,9))

dx1 <- data.frame(Treatment = rep("T1", length(x1)), Temp = rep(10, length(x1)), Rep = rep(c(1,2,3), times = c(50,35,40)), Meas = x1)
dx2 <- data.frame(Treatment = rep("T1", length(x2)), Temp = rep(20, length(x2)), Rep = rep(c(4,5,6), times = c(43,32,40)), Meas = x2)
dx3 <- data.frame(Treatment = rep("T2", length(x3)), Temp = rep(10, length(x3)), Rep = rep(c(7,8,9), times = c(50,50,50)), Meas = x3)
dx4 <- data.frame(Treatment = rep("T2", length(x4)), Temp = rep(20, length(x4)), Rep = rep(c(10,11,12), times = c(26,43,45)), Meas = x4)

# Entire data set

dat <- rbind(dx1,dx2,dx3,dx4)

# Plot overview

w <- ggplot(dat, aes(x = factor(Rep), y = Meas))
w + geom_boxplot(aes(fill = factor(Temp)))

# Averages (with se)

aveg <- aggregate(Meas ~ Treatment + Temp + Rep, data = dat, FUN = mean)
se <- function(x) sd(x)/sqrt(length(x))
SE <- aggregate(Meas ~ Treatment + Temp + Rep, data = dat, FUN = se)

dat2 <- merge(aveg, SE, by = c("Treatment", "Temp", "Rep"), sort = F)

colnames(dat2)[colnames(dat2) %in% grep("\\.x", colnames(dat2), value = T)] <- "mean"
colnames(dat2)[colnames(dat2) %in% grep("\\.y", colnames(dat2), value = T)] <- "se"

# Plot entire data set

p <- ggplot(dat2, aes(x = Treatment, y = mean, ymax = mean + se/2,
ymin = mean - se/2))

p + geom_pointrange(aes(color = factor(Temp)),
position=position_dodge(width=0.50), size = 1)

For simplicity, first consider only one group. Your data can be modeled as $$Y_{ij} = \mu + \delta_i + \epsilon_{ij},$$ where $\delta_i \sim N(0, \tau^2), i=1,\ldots,k$ is the random day effect, and $\epsilon_{ij}\sim N(0, \sigma^2), j=1,\ldots, n_i$ is the random within-day replicate effect. Then you can obtain an unbiased estimate of the group average $$\hat\mu = \frac 1k \sum_{j=1}^j \bar{Y}_{i.} \sim N(\mu, \frac{\tau^2}{k} + \frac 1k \sum_{i=1}^k \frac{\sigma^2}{n_i}).$$ It is fairly simple to obtain estimates of the variance term with $\tau^2$ being the variance of the day means, and $\sigma^2$ the pooled within-day variance.

Start edit

The variance term in this formula gives you the square of the standard error for your plot. The point to show would be the mean of the within-day means. If you don't want to assume that the within-day variability is constant, just replace $\sigma$ with $\sigma_i$ - the within-day standard deviation on day 1, so the summation would be $\sum_iSE_i^2$, where $SE_i$ is the standard error of the within-day mean.

End edit

If for a moment you assume $n_i=n$, then the variance is $\frac{\tau^2}{k} + \frac{\sigma^2}{n}$. In your application the key observation is that $n >> k$, so the first term will dominate. In essence, you can ignore the within-day replicates, and just use them as an expensive way to obtain one (day-specific) observation.

I have run a "correct" mixed ANOVA analysis and a simple ANOVA on group means, and they give essentially the same result:

library(lme4)
m1 <- lmer(Meas ~ Treatment*Temp + (1|Rep), data=dat)
m1

Output:

Linear mixed model fit by REML
Formula: Meas ~ Treatment * Temp + (1 | Rep)
Data: dat
AIC  BIC logLik deviance REMLdev
3351 3376  -1669     3343    3339
Random effects:
Groups   Name        Variance Std.Dev.
Rep      (Intercept) 15.414   3.9260
Residual             42.000   6.4808
Number of obs: 504, groups: Rep, 12

Fixed effects:
Estimate Std. Error t value
(Intercept)       16.1579     5.2369   3.085
TreatmentT2       -9.3838     7.3903  -1.270
Temp              -0.4723     0.3315  -1.425
TreatmentT2:Temp   1.0377     0.4683   2.216

And

m2 <- lm(mean ~ Treatment*Temp, data=dat2)
summary(m2)

with output

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)       16.1037     5.2451   3.070   0.0153 *
TreatmentT2       -9.2846     7.4177  -1.252   0.2460
Temp              -0.4694     0.3317  -1.415   0.1948
TreatmentT2:Temp   1.0302     0.4691   2.196   0.0594 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.063 on 8 degrees of freedom
Multiple R-squared: 0.5955,     Adjusted R-squared: 0.4438
F-statistic: 3.926 on 3 and 8 DF,  p-value: 0.05412
• Thanks Aniko! This confirms that 2-way ANOVA is not much worse an approach than a mixed model in this case. I'll need a couple of days to digest your excellent answer, but as far as I understand you don't give a suggestion how to combine the variation for a plot? I do remember you telling just to show the points, but maybe one point would be easier to read, only if it was possible to compress the variation on two levels so that it shows the same than ANOVA / mixed model results? – Mikko May 2 '12 at 10:07
• I have added an explanation on how to use the formula to get a standard error, which in fact would be very close to just the standard error of the three within-day means. – Aniko May 2 '12 at 12:35