¿Cómo puedo ajustar un modelo lineal de efectos mixtos para medidas repetidas de datos utilizando nlmer()?

Question

Estoy tratando de analizar medidas repetidas de datos y estoy luchando para hacer que funcione en R. Mis datos es esencialmente la siguiente, tengo dos grupos de tratamiento. Cada sujeto en cada grupo se ponen a prueba cada día y le asigna un puntaje (el porcentaje correcto en una prueba). Los datos en el 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 asemeja a una curva logística, los sujetos que hacerlo muy mal para un par de días, seguido por una rápida mejoría, seguido de una meseta. Me gustaría saber si el tratamiento tiene un efecto sobre el rendimiento de la prueba de la curva. Mi idea fue usar nlmer() en la lme4 paquete en R. Puedo líneas de ajuste para cada grupo utilizando el siguiente:

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

Yo puede que comparar los grupos por mirar las estimaciones para los diferentes parámetros y desviaciones estándar de la estimación de líneas, pero no estoy seguro de que esta es la manera correcta de hacerlo. Cualquier ayuda sería muy apreciada.

Answer 1

Usted puede utilizar normal de probabilidad de la relación de las pruebas. He aquí un ejemplo simple. En primer lugar, vamos a crear las observaciones a partir de 10 personas con base 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 vamos a la mitad de ellos tienen diferentes asíntotas y el punto medio de los parámetros:

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 los valores de las respuestas para todas las personas, basada en 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, las diferencias que los modelos deben ser capaces de recoger. Ahora vamos a probar primero para adaptarse a un simple modelo, haciendo caso omiso de los 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 era de esperar, las estimaciones para Asym y xmid están en algún lugar entre los verdaderos valores de los parámetros para los dos grupos. (Este sería el caso no es obvio, sin embargo, ya que el parámetro de escala es también cambió, para ajustar el modelo misspecification.) Ahora vamos a ajustar un completo modelo, con parámetros diferentes 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

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

> 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 extremadamente pequeño p-valor muestra claramente que el modelo simple era demasiado simple; los dos grupos hacen difieren en sus parámetros.

Sin embargo, los dos de la escala de los parámetros estimados son casi idénticos, con una diferencia de sólo .1. Tal vez necesitamos sólo necesita un parámetro de escala? (Por supuesto, sabemos que la respuesta es sí, ya que contamos con datos simulados.)

(La diferencia entre los dos asíntota parámetros también está a sólo .1, pero eso es una gran diferencia cuando tomamos los errores estándar en cuenta – ver summary(fm2).)

Así que el ajuste de un nuevo modelo, con un común scale de parámetro para los dos grupos, pero con diferentes Asym y xmid pará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 puesto que el modelo reducido es anidada en el modelo completo, se puede volver a hacer una prueba de razón de verosimilitud:

> 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 gran p-valor indica que el modelo reducido se adapta así como el modelo completo, como se esperaba.

Por supuesto podemos hacer pruebas similares para comprobar si los diferentes valores de los parámetros son necesarios para que sólo Asym, sólo xmid o ambos. Dicho esto, yo no recomiendo hacer de regresión paso a paso como este para eliminar parámetros. En lugar de eso, prueba el modelo completo (fm2) contra el modelo simple (fm1), y ser feliz con los resultados. Para cuantificar las diferencias, las parcelas será de ayuda.