Jak dopasować nieliniowy model efektów mieszanych dla danych z powtarzanymi pomiarami za pomocą nlmer ()?

Próbuję analizować dane z powtarzanych pomiarów i staram się, aby to zadziałało R. Moje dane są zasadniczo następujące, mam dwie grupy leczenia. Każdy przedmiot w każdej grupie jest codziennie testowany i otrzymuje wynik (procent poprawny na teście). Dane są w długim formacie:

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

Dane przypominają krzywą logistyczną, badani radzą sobie bardzo słabo przez kilka dni, po czym następuje szybka poprawa, a następnie plateau. Chciałbym wiedzieć, czy leczenie ma wpływ na krzywą wydajności testu. Myślałem o użyciu nlmer()w lme4pakiecie w R. Mogę dopasować linie dla każdej grupy, używając:

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

Mogę porównywać grupy, patrząc na szacunki dla różnych parametrów i odchylenia standardowe szacowanych linii, ale nie jestem pewien, czy jest to właściwy sposób, aby to zrobić. Każda pomoc byłaby bardzo mile widziana.

Możesz użyć normalnych testów współczynnika wiarygodności. Oto prosty przykład. Najpierw utwórzmy obserwacje od 10 osób na podstawie twoich parametrów:

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)))

Teraz połowa z nich ma różne asymptoty i parametry punktu środkowego:

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)

Możemy symulować wartości odpowiedzi dla wszystkich osób w oparciu o model:

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")

Widzimy wyraźne różnice między dwiema grupami, różnice, które modele powinny być w stanie wychwycić. Teraz spróbujmy najpierw dopasować prosty model, ignorując grupy:

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

Być może zgodnie z oczekiwaniami, szacunki Asymi xmidsą gdzieś pomiędzy rzeczywistymi wartościami parametrów dla dwóch grup. (To nie byłoby oczywiste, ponieważ zmieniono również parametr skali, aby dostosować się do błędnej specyfikacji modelu.) Teraz dopasujmy pełny model z różnymi parametrami dla dwóch grup:

> 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

Ponieważ oba modele są zagnieżdżone, możemy wykonać test współczynnika wiarygodności:

> 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

Niezwykle mała wartość p wyraźnie pokazuje, że prosty model był zbyt prosty; dwie grupy nie różnią się pod względem parametrów.

Jednak szacunki dwóch parametrów skali są prawie identyczne, z różnicą zaledwie .1. Być może potrzebujemy tylko jednego parametru skali? (Oczywiście wiemy, że odpowiedź brzmi tak, ponieważ symulowaliśmy dane.)

(Różnica między dwoma asymptotycznymi parametrami jest również tylko .1, ale to duża różnica, jeśli weźmiemy pod uwagę standardowe błędy - patrz summary(fm2).)

Więc dopasować nowy model, z wspólny scaleparametr dla dwóch grup, ale inny Asymi xmidparametrów, jak poprzednio:

> 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

A ponieważ model zredukowany jest zagnieżdżony w pełnym modelu, możemy ponownie wykonać test współczynnika wiarygodności:

> 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

Duża wartość p wskazuje, że zredukowany model pasuje, podobnie jak pełny model, zgodnie z oczekiwaniami.

Możemy oczywiście wykonać podobne testy, aby sprawdzić, czy potrzebne są różne wartości parametrów dla just Asym, just xmidczy oba. To powiedziawszy, nie zalecałbym wykonywania regresji krokowej w ten sposób, aby wyeliminować parametry. Zamiast tego po prostu przetestuj pełny model ( fm2) na prostym modelu ( fm1) i ciesz się z wyników. Aby oszacować różnice, pomocne będą wykresy.