Как мне подогнать нелинейную модель смешанных эффектов для данных повторных измерений, используя nlmer ()?

Я пытаюсь проанализировать данные повторных измерений и пытаюсь заставить их работать R. Мои данные по существу следующие, у меня есть две группы лечения. Каждый субъект в каждой группе тестируется каждый день и получает оценку (процент правильных на тесте). Данные в длинном формате:

Time Percent Subject   Group
   1       0    GK11 Ethanol
   2       0    GK11 Ethanol
   3       0    GK11 Ethanol
   4       0    GK11 Ethanol
   5       0    GK11 Ethanol
   6       0    GK11 Ethanol

Данные напоминают логистическую кривую, субъекты в течение нескольких дней чувствуют себя очень плохо, за ними следует быстрое улучшение, за которым следует плато. Я хотел бы знать, влияет ли лечение на кривую эффективности теста. Моя мысль была использовать nlmer()в lme4пакете R. Я могу подобрать строки для каждой группы, используя следующее:

print(nm1 <- nlmer(Percent ~ SSlogis(Time,Asym, xmid, scal) ~ Asym | Subject,
salinedata, start = c(Asym =.60,  xmid = 23, scal = 5)), corr = FALSE)

Я могу сравнить группы, посмотрев оценки для различных параметров и стандартных отклонений оценочных линий, но я не уверен, что это правильный способ сделать это. Любая помощь будет принята с благодарностью.

Вы можете использовать нормальные тесты отношения правдоподобия. Вот простой пример. Во-первых, давайте создадим наблюдения от 10 человек на основе ваших параметров:

Asym = .6
xmid = 23
scal = 5

n = 10
time = seq(1,60,5)

d = data.frame(time=rep(time,10),
               Asym, xmid, scal, group=0)
d$subj = factor(rep(1:n, each=length(time)))

Теперь пусть половина из них имеет разные асимптоты и параметры средней точки:

ind = (nrow(d)/2):nrow(d)
d$Asym[ind] = d$Asym[ind] + .1
d$xmid[ind] = d$xmid[ind] + 10
d$group[ind] = 1
d$group=factor(d$group)

Мы можем смоделировать значения ответов для всех людей на основе модели:

set.seed(1)
d = transform(d, y = Asym/(1+exp((xmid-time)/scal)) +
                     rnorm(nrow(d), sd=.04))
library(lattice)
xyplot(y~time | group, group=subj,
       data=d, type=c("g","l"), col="black")

Мы можем видеть четкие различия между двумя группами, различия, которые модели должны быть в состоянии уловить. Теперь давайте сначала попробуем подобрать простую модель, игнорируя группы:

> fm1 = nls(y ~ SSlogis(time, Asym, xmid, scal), data=d)
> coef(fm1)
      Asym       xmid       scal 
 0.6633042 28.5219166  5.8286082

Возможно, как и ожидалось, оценки Asymи xmidнаходятся где-то между реальными значениями параметров для двух групп. (То, что это имело бы место, не очевидно, хотя, поскольку параметр масштаба также изменен, чтобы приспособиться к неправильной спецификации модели.) Теперь давайте подберем полную модель с различными параметрами для двух групп:

> fm2 = nls(y ~ SSlogis(time, Asym[group], xmid[group], scal[group]),
          data=d,
          start=list(Asym=rep(.6,2), xmid=rep(23,2), scal=rep(5,2)))
> coef(fm2)
    Asym1     Asym2     xmid1     xmid2     scal1     scal2 
 0.602768  0.714199 22.769315 33.331976  4.629332  4.749555

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

> anova(fm1, fm2)
Analysis of Variance Table

Model 1: y ~ SSlogis(time, Asym, xmid, scal)
Model 2: y ~ SSlogis(time, Asym[group], xmid[group], scal[group])
  Res.Df Res.Sum Sq Df  Sum Sq F value    Pr(>F)    
1    117    0.70968                                 
2    114    0.13934  3 0.57034  155.54 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Чрезвычайно малое значение p ясно показывает, что простая модель была слишком простой; две группы действительно отличаются по своим параметрам.

Однако две оценки параметров шкалы практически идентичны, с разницей всего в 0,1. Возможно, нам нужен только один масштабный параметр? (Конечно, мы знаем, что ответ - да, так как мы смоделировали данные.)

(Разница между двумя параметрами асимптоты также составляет всего 0,1, но это большая разница, когда мы учитываем стандартные ошибки - см summary(fm2).)

Итак, мы подходим к новой модели, с общими scaleпараметрами для двух групп, но разными Asymи xmidпараметрами, как и раньше:

> fm3 = nls(y ~ SSlogis(time, Asym[group], xmid[group], scal),
          data=d,
          start=list(Asym=rep(.6,2), xmid=rep(23,2), scal=5))
> coef(fm3)
     Asym1      Asym2      xmid1      xmid2       scal 
 0.6035251  0.7129002 22.7821155 33.3080264  4.6928316

А поскольку уменьшенная модель вложена в полную модель, мы снова можем выполнить тест отношения правдоподобия:

> anova(fm3, fm2)
Analysis of Variance Table

Model 1: y ~ SSlogis(time, Asym[group], xmid[group], scal)
Model 2: y ~ SSlogis(time, Asym[group], xmid[group], scal[group])
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1    115    0.13945                             
2    114    0.13934  1 0.00010637   0.087 0.7685

Большое значение p указывает, что уменьшенная модель подходит так же, как и полная модель, как и ожидалось.

Конечно, мы можем сделать аналогичные тесты, чтобы проверить, нужны ли разные значения параметров для just Asym, just xmidили обоих. Тем не менее, я бы не рекомендовал делать пошаговую регрессию таким образом, чтобы исключить параметры. Вместо этого просто протестируйте полную модель ( fm2) с простой моделью ( fm1), и будьте довольны результатами. Для количественной оценки любых различий, графики будут полезны.