# Unbalanced mixed effect ANOVA for repeated measures

I have data from patients treated with 2 different kinds of treatments during surgery. I need to analyze its effect on heart rate. The heart rate measurement is taken every 15 minutes.

Given that the surgery length can be different for each patient, each patient can have between 7 and 10 heart rate measurements. So an unbalanced design should be used. I'm doing my analysis using R. And have been using the ez package to do repeated measure mixed effect ANOVA. But I do not know how to analyse unbalanced data. Can anyone help?

Suggestions on how to analyze the data are also welcomed.

Update:
As suggested, I fitted the data using the lmer function and found that the best model is:

heart.rate~ time + treatment + (1|id) + (0+time|id) + (0+treatment|time)


with the following result:

Random effects:
Groups   Name        Variance   Std.Dev. Corr
id       time        0.00037139 0.019271
id       (Intercept) 9.77814104 3.127002
time     treat0      0.09981062 0.315928
treat1      1.82667634 1.351546 -0.504
Residual             2.70163305 1.643665
Number of obs: 378, groups: subj, 60; time, 9

Fixed effects:
Estimate Std. Error t value
(Intercept) 72.786396   0.649285  112.10
time         0.040714   0.005378    7.57
treat1       2.209312   1.040471    2.12

Correlation of Fixed Effects:
(Intr) time
time   -0.302
treat1 -0.575 -0.121


Now I'm lost at interpreting the result. Am I right in concluding that the two treatments differed in affecting heart rate? What does the correlation of -504 between treat0 and treat1 means?

-
Before I update the answer, is treatment a repeated factor? i.e., does each subject receive both treatment "a" and treatment "b" or is this a between-subjects factor? –  Matt Albrecht Oct 28 '10 at 8:46
Treatment is a between-subjects factor. Each subject only receive 1 kind of treatment. I've coded the two treatments as 1 and 0 and set treatment as a factor variable. –  biostat_newbie Oct 28 '10 at 9:08

The lme/lmer functions from the nlme/lme4 packages are able to deal with unbalanced designs. You should make sure that time is a numeric variable. You would also probably want to test for different types of curves as well. The code will look something like this:

library(lme4)
#plot data with a plot per person including a regression line for each
xyplot(heart.rate ~ time|id, groups=treatment, type= c("p", "r"), data=heart)

#Mixed effects modelling
#variation in intercept by participant
lmera.1 <- lmer(heart.rate ~ treatment * time + (1|id), data=heart)
#variation in intercept and slope without correlation between the two
lmera.2 <- lmer(heart.rate ~ treatment * time + (1|id) + (0+time|id), data=heart)
#As lmera.1 but with correlation between slope and intercept
lmera.3 <- lmer(heart.rate ~ treatment * time + (1+time|id), data=heart)

#Determine which random effects structure fits the data best
anova(lmera.1, lmera.2, lmera.3)


To get quadratic models use the formula "heart.rate ~ treatment * time * I(time^2) + (random effects)".

Update:
In this case where treatment is a between-subjects factor, I would stick with the model specifications above. I don't think the term (0+treatment|time) is one that you want included in the model, to me it doesn't make sense in this instance to treat time as a random-effects grouping variable.

But to answer your question of "what does the correlation -0.504 mean between treat0 and treat1" this is the correlation coefficient between the two treatments where each time grouping is one pair of values. This makes more sense if id is the grouping factor and treatment is a within-subjects variable. Then you have an estimate of the correlation between the intercepts of the two conditions.

Before making any conclusions about the model, refit it with lmera.2 and include REML=F. Then load the "languageR" package and run:

plmera.2<-pvals.fnc(lmera.2)
plmera.2


Then you can get p-values, but by the looks of it, there is probably a significant effect of time and a significant effect of treatment.

-
Should one set lmer's REML argument to FALSE when generating those models since they'll eventually be compared using the anova() function? –  Mike Lawrence Oct 22 '10 at 5:47
When comparing models using likelihood-ratio tests, you can compare different random effects structures using REML (restricted/residual maximum likelihood, as above), but you must use ML (maximum likelihood) to compare different fixed effect models. –  onestop Oct 22 '10 at 7:36
@onestop: Ah! Good to know, thanks! –  Mike Lawrence Oct 22 '10 at 12:56
The term '(time|id)' on the random effects side tells the function to fit different (linear) slopes for each person. So you can have time as both a fixed effect and a random effect, but they mean different things. Have a look at the sleepstudy example in Douglas Bates' book: lme4.r-forge.r-project.org/book/Ch4.pdf –  Matt Albrecht Oct 26 '10 at 13:47
I don't understand how for repeated measures designs, lmer is suggested rather than the good old lme. In such desgins crossed random effects, the main strenght for lmer, are rare but quite often you want to model the correlation structure of the residuals. As far as I understand lmer does not support that but lme does. Am I wrong to assume in such cases lmer is an inferior tool compared to lme? –  AlefSin Apr 8 '11 at 12:32