Nested design with pairing

Question

I have been asked to help a colleague with analysing the following experimental design.

2 treatments (fixed), 5 time points (fixed) with 4 animals (random) at each time point. The animals within each time point are not the same.

Each animal had two samples taken, one from a treated limb and one from an untreated limb so within each time point the samples from each animal are paired.

I would like to analyse this data as a nested ANOVA in R. The formula I have currently is:

value ~ treatment + time/individual

Testing the measured value for the treatment at each time with individual nested within time. This doesn't account for the paired nature of the samples though. Also this is the first example of a nested design I've been asked to analyse and I wonder if someone could advise me as to whether I'm on the right path.

Answer 1

Based on what you have described, I would use the mixed model:

Value = intercept + treatment + time + treatment*time + individal + e

Where treatment, time and their interaction are fixed effects and individual is a random effect.

The way to formulate this model depends on the package you are using. If you are using lme() in package {nlme} then you can fit this model with:

lme(fixed = value ~ treatment * time, data = my.data ,  random = ~ 1 | individual)

Where my.data is the data.frame containing treatment, time and individual as factor variables and value as a numeric variable. This formulation is short form; it’s actually the same as:

lme(fixed = value ~ 1 + treatment + time + treatment:time, data = my.data, random = ~ 1 | individual)

which corresponds more directly to the model above. Make sure that your data.frame is ordered correctly and is univariate (vertical).

This model has two variance parameters: one for the variance explained by individual animals and one for the residual error. If you want to estimate individual animal variance at each time then you could use a model like:

lme(fixed = value ~ treatment * time, data = my.data, random = ~ time | individual)

but this might be over-parameterized for your data set and you might have convergence issues and the final hessian might not be positive definite.