¿Cómo ajusto un modelo de efectos mixtos no lineales para datos de medidas repetidas usando nlmer ()?

Estoy tratando de analizar datos de medidas repetidas y estoy luchando para que funcione R. Mis datos son esencialmente los siguientes, tengo dos grupos de tratamiento. Todos los sujetos de cada grupo se evalúan todos los días y se les asigna una puntuación (el porcentaje correcto en una prueba). Los datos están en formato largo:

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

Los datos se asemejan a una curva logística, los sujetos lo hacen muy mal durante unos días seguido de una mejora rápida, seguida de una meseta. Me gustaría saber si el tratamiento tiene un efecto en la curva de rendimiento de la prueba. Mi pensamiento era usarlo nlmer()en el lme4paquete R. Puedo ajustar líneas para cada grupo usando lo siguiente:

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

Puedo comparar grupos mirando las estimaciones de los diferentes parámetros y las desviaciones estándar de las líneas estimadas, pero no estoy seguro de que esta sea la forma correcta de hacerlo. Cualquier ayuda sería muy apreciada.

Puede usar pruebas de razón de probabilidad normal. Aquí hay un ejemplo simple. Primero, creemos observaciones de 10 personas basadas en sus parámetros:

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

Ahora deje que la mitad de ellos tengan diferentes asíntotas y parámetros de punto medio:

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)

Podemos simular valores de respuesta para todos los individuos, según el modelo:

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

Podemos ver claras diferencias entre los dos grupos, diferencias que los modelos deberían poder detectar. Ahora intentemos primero ajustar un modelo simple , ignorando grupos:

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

Tal vez como se esperaba, las estimaciones Asymy xmidestán en algún lugar entre los valores de parámetros reales para los dos grupos. (Sin embargo, este sería el caso no es obvio, ya que el parámetro de escala también se cambia para ajustar la especificación errónea del modelo). Ahora ajustemos un modelo completo , con diferentes parámetros para los dos grupos:

> 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

Como los dos modelos están anidados, podemos hacer una prueba de razón de probabilidad:

> 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

El valor p extremadamente pequeño muestra claramente que el modelo simple era demasiado simple; los dos grupos no difieren en sus parámetros.

Sin embargo, las estimaciones de los dos parámetros de escala son casi idénticas, con una diferencia de solo .1. ¿Quizás solo necesitamos un parámetro de escala? (Por supuesto, sabemos que la respuesta es sí, ya que tenemos datos simulados).

(La diferencia entre los dos parámetros de asíntota también es solo .1, pero esa es una gran diferencia cuando tomamos en cuenta los errores estándar, ver summary(fm2)).

Entonces ajustamos un nuevo modelo, con un scaleparámetro común para los dos grupos, pero diferentes Asymy xmidparámetros, como antes:

> 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

Y dado que el modelo reducido está anidado en el modelo completo, nuevamente podemos hacer una prueba de razón de probabilidad:

> 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

El valor p grande indica que el modelo reducido se ajusta tan bien como el modelo completo, como se esperaba.

Por supuesto, podemos hacer pruebas similares para verificar si se necesitan valores de parámetros diferentes para just Asym, just xmido both. Dicho esto, no recomendaría hacer una regresión gradual como esta para eliminar los parámetros. En cambio, solo pruebe el modelo completo ( fm2) contra el modelo simple ( fm1), y esté satisfecho con los resultados. Para cuantificar las diferencias, las gráficas serán útiles.