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.

enter image description here 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:


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:

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


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


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

with output

                 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 
  • $\begingroup$ 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? $\endgroup$
    – Mikko
    May 2 '12 at 10:07
  • $\begingroup$ 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. $\endgroup$
    – Aniko
    May 2 '12 at 12:35

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.