Unbalanced mixed effect ANOVA for repeated measures

Question

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?

Answer 1

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.