Comment ajuster un modèle d'effets mixtes non linéaires pour les données de mesures répétées à l'aide de nlmer ()?

J'essaie d'analyser les données de mesures répétées et j'ai du mal à les faire fonctionner R. Mes données sont essentiellement les suivantes, j'ai deux groupes de traitement. Chaque sujet de chaque groupe est testé tous les jours et obtient un score (le pourcentage correct sur un test). Les données sont au format long:

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

Les données ressemblent à une courbe logistique, les sujets se portent très mal pendant quelques jours suivis d'une amélioration rapide, suivie d'un plateau. Je voudrais savoir si le traitement a un effet sur la courbe de performance du test. Ma pensée était d'utiliser nlmer()dans le lme4package en R. Je peux ajuster des lignes pour chaque groupe en utilisant les éléments suivants:

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

Je peux comparer les groupes en examinant les estimations des différents paramètres et écarts-types des lignes estimées, mais je ne suis pas sûr que ce soit la bonne façon de procéder. Toute aide serait grandement appréciée.

Vous pouvez utiliser des tests de rapport de vraisemblance normaux. Voici un exemple simple. Commençons par créer des observations à partir de 10 individus en fonction de vos paramètres:

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

Maintenant, la moitié d'entre eux ont des asymptotes et des paramètres médians différents:

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)

Nous pouvons simuler des valeurs de réponse pour tous les individus, sur la base du modèle:

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

Nous pouvons voir des différences claires entre les deux groupes, différences que les modèles devraient être capables de détecter. Essayons maintenant d'adapter un modèle simple , en ignorant les groupes:

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

Peut-être comme prévu, les estimations pour Asymet se xmidsituent quelque part entre les valeurs réelles des paramètres pour les deux groupes. (Ce ne serait pas le cas, cependant, puisque le paramètre d'échelle est également modifié, pour ajuster la spécification incorrecte du modèle.) A présent, ajustons un modèle complet , avec des paramètres différents pour les deux groupes:

> 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

Puisque les deux modèles sont imbriqués, nous pouvons faire un test de rapport de vraisemblance:

> 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

La valeur p extrêmement faible montre clairement que le modèle simple était trop simple; les deux groupes ne diffèrent dans leurs paramètres.

Cependant, les deux estimations des paramètres d'échelle sont presque identiques, avec une différence de seulement 0,1. Peut-être avons-nous besoin d'un seul paramètre d'échelle? (Bien sûr, nous savons que la réponse est oui, car nous avons des données simulées.)

(La différence entre les deux paramètres asymptotes n'est également que de 0,1, mais c'est une grande différence lorsque nous tenons compte des erreurs standard - voir summary(fm2).)

Donc , nous nous situons un nouveau modèle, avec un commun scaleparamètre pour les deux groupes, mais différents Asymet xmidparamètres, comme avant:

> 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

Et puisque le modèle réduit est imbriqué dans le modèle complet, nous pouvons à nouveau faire un test de rapport de vraisemblance:

> 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

La grande valeur p indique que le modèle réduit convient aussi bien que le modèle complet, comme prévu.

Nous pouvons bien sûr faire des tests similaires pour vérifier si différentes valeurs de paramètres sont nécessaires pour juste Asym, juste xmidou les deux. Cela dit, je ne recommanderais faire une régression comme celle-ci pour éliminer les paramètres. Au lieu de cela, testez simplement le modèle complet ( fm2) par rapport au modèle simple ( fm1) et soyez satisfait des résultats. Pour quantifier les différences, les graphiques seront utiles.