Несбалансированный смешанный эффект ANOVA для повторных измерений

У меня есть данные от пациентов, получавших 2 разных вида лечения во время операции. Мне нужно проанализировать его влияние на частоту сердечных сокращений. Измерение частоты сердечных сокращений проводится каждые 15 минут.

Учитывая, что продолжительность операции может быть разной для каждого пациента, у каждого пациента может быть от 7 до 10 измерений сердечного ритма. Поэтому следует использовать несбалансированный дизайн. Я делаю свой анализ с использованием R. И я использую пакет ez для повторного измерения смешанного эффекта ANOVA. Но я не знаю, как анализировать несбалансированные данные. Кто-нибудь может помочь?

Предложения о том, как анализировать данные, также приветствуются.

Обновление:
Как предложено, я подгонял данные, используя lmerфункцию, и обнаружил, что лучшая модель:

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

со следующим результатом:

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

Теперь я потерял в интерпретации результата. Правильно ли я пришел к выводу, что два метода лечения влияли на частоту сердечных сокращений? Что означает корреляция -504 между Treat0 и Treat1?

Функции lme / lmer из пакетов nlme / lme4 способны работать с несбалансированными конструкциями. Вы должны убедиться, что время является числовой переменной. Вы также, вероятно, захотите проверить различные типы кривых. Код будет выглядеть примерно так:

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)

Для получения квадратичных моделей используйте формулу «сердечный ритм ~ лечение * время * I (время ^ 2) + (случайные эффекты)».

Обновление:
в этом случае, когда лечение является фактором между субъектами, я бы придерживался приведенных выше характеристик модели. Я не думаю, что термин (0 + обработка | время) является тем, который вы хотите включить в модель, для меня не имеет смысла рассматривать время как переменную группировки случайных эффектов.

Но, чтобы ответить на ваш вопрос о том, «что означает корреляция -0,504 между Treat0 и Treat1 », это коэффициент корреляции между двумя обработками, где каждая временная группировка представляет собой одну пару значений. Это имеет больше смысла, если id является фактором группировки, а лечение является переменной внутри субъекта. Затем у вас есть оценка корреляции между перехватами двух условий.

Прежде чем делать какие-либо выводы о модели, замените ее на lmera.2 и включите REML = F. Затем загрузите пакет languageR и запустите:

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

Тогда вы можете получить p-значения, но, судя по всему, есть значительный эффект времени и значительный эффект лечения.