Regresión lineal con medidas repetidas en R

12

No pude descubrir cómo realizar una regresión lineal en R para un diseño de medida repetida. En una pregunta anterior (aún sin respuesta) me sugirieron que no usara lmsino que usara modelos mixtos. Utilicé lmde la siguiente manera:

lm.velocity_vs_Velocity_response <- lm(Velocity_response~Velocity*Subject, data=mydata)

(se pueden encontrar más detalles sobre el conjunto de datos en el enlace de arriba)

Sin embargo, no pude encontrar en Internet ningún ejemplo con código R que mostrara cómo realizar un análisis de regresión lineal.

Lo que quiero es, por un lado, un gráfico de los datos con la línea que se ajuste a los datos y, por otro lado, el valor junto con el valor p para la prueba de significación para el modelo.R2

¿Hay alguien que pueda proporcionar algunas sugerencias? Cualquier ejemplo de código R podría ser de gran ayuda.


Editar
De acuerdo con la sugerencia que recibí hasta ahora, la solución para analizar mis datos para comprender si existe una relación lineal entre las dos variables Velocity_response (derivada del cuestionario) y Velocity (derivada del rendimiento) debería ser esta:

library(nlme)
summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))

El resultado del resumen da esto:

    > summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))
    Linear mixed-effects model fit by REML
     Data: scrd 
           AIC      BIC   logLik
      104.2542 126.1603 -30.1271

    Random effects:
     Formula: ~1 | Subject
            (Intercept) Residual
    StdDev:    2.833804 2.125353

Fixed effects: Velocity_response ~ Velocity * Subject 
                              Value Std.Error DF    t-value p-value
(Intercept)               -26.99558  25.82249 20 -1.0454288  0.3083
Velocity                   24.52675  19.28159 20  1.2720292  0.2180
SubjectSubject10           21.69377  27.18904  0  0.7978865     NaN
SubjectSubject11           11.31468  33.51749  0  0.3375754     NaN
SubjectSubject13           52.45966  53.96342  0  0.9721337     NaN
SubjectSubject2           -14.90571  34.16940  0 -0.4362299     NaN
SubjectSubject3            26.65853  29.41574  0  0.9062674     NaN
SubjectSubject6            37.28252  50.06033  0  0.7447517     NaN
SubjectSubject7            12.66581  26.58159  0  0.4764880     NaN
SubjectSubject8            14.28029  31.88142  0  0.4479188     NaN
SubjectSubject9             5.65504  34.54357  0  0.1637076     NaN
Velocity:SubjectSubject10 -11.89464  21.07070 20 -0.5645111  0.5787
Velocity:SubjectSubject11  -5.22544  27.68192 20 -0.1887672  0.8522
Velocity:SubjectSubject13 -41.06777  44.43318 20 -0.9242591  0.3664
Velocity:SubjectSubject2   11.53397  25.41780 20  0.4537754  0.6549
Velocity:SubjectSubject3  -19.47392  23.26966 20 -0.8368804  0.4125
Velocity:SubjectSubject6  -29.60138  41.47500 20 -0.7137162  0.4836
Velocity:SubjectSubject7   -6.85539  19.92271 20 -0.3440992  0.7344
Velocity:SubjectSubject8  -12.51390  22.54724 20 -0.5550080  0.5850
Velocity:SubjectSubject9   -2.22888  27.49938 20 -0.0810519  0.9362
 Correlation: 
                          (Intr) Velcty SbjS10 SbjS11 SbjS13 SbjcS2 SbjcS3 SbjcS6 SbjcS7 SbjcS8 SbjcS9 V:SS10 V:SS11 V:SS13 Vl:SS2 Vl:SS3
Velocity                  -0.993                                                                                                         
SubjectSubject10          -0.950  0.943                                                                                                  
SubjectSubject11          -0.770  0.765  0.732                                                                                           
SubjectSubject13          -0.479  0.475  0.454  0.369                                                                                    
SubjectSubject2           -0.756  0.751  0.718  0.582  0.362                                                                             
SubjectSubject3           -0.878  0.872  0.834  0.676  0.420  0.663                                                                      
SubjectSubject6           -0.516  0.512  0.490  0.397  0.247  0.390  0.453                                                               
SubjectSubject7           -0.971  0.965  0.923  0.748  0.465  0.734  0.853  0.501                                                        
SubjectSubject8           -0.810  0.804  0.769  0.624  0.388  0.612  0.711  0.418  0.787                                                 
SubjectSubject9           -0.748  0.742  0.710  0.576  0.358  0.565  0.656  0.386  0.726  0.605                                          
Velocity:SubjectSubject10  0.909 -0.915 -0.981 -0.700 -0.435 -0.687 -0.798 -0.469 -0.883 -0.736 -0.679                                   
Velocity:SubjectSubject11  0.692 -0.697 -0.657 -0.986 -0.331 -0.523 -0.607 -0.357 -0.672 -0.560 -0.517  0.637                            
Velocity:SubjectSubject13  0.431 -0.434 -0.409 -0.332 -0.996 -0.326 -0.378 -0.222 -0.419 -0.349 -0.322  0.397  0.302                     
Velocity:SubjectSubject2   0.753 -0.759 -0.715 -0.580 -0.360 -0.992 -0.661 -0.389 -0.732 -0.610 -0.563  0.694  0.528  0.329              
Velocity:SubjectSubject3   0.823 -0.829 -0.782 -0.634 -0.394 -0.622 -0.984 -0.424 -0.799 -0.667 -0.615  0.758  0.577  0.360  0.629       
Velocity:SubjectSubject6   0.462 -0.465 -0.438 -0.356 -0.221 -0.349 -0.405 -0.995 -0.449 -0.374 -0.345  0.425  0.324  0.202  0.353  0.385
Velocity:SubjectSubject7   0.961 -0.968 -0.913 -0.740 -0.460 -0.726 -0.844 -0.496 -0.986 -0.778 -0.718  0.886  0.674  0.420  0.734  0.802
Velocity:SubjectSubject8   0.849 -0.855 -0.807 -0.654 -0.406 -0.642 -0.746 -0.438 -0.825 -0.988 -0.635  0.783  0.596  0.371  0.649  0.709
Velocity:SubjectSubject9   0.696 -0.701 -0.661 -0.536 -0.333 -0.526 -0.611 -0.359 -0.676 -0.564 -0.990  0.642  0.488  0.304  0.532  0.581
                          Vl:SS6 Vl:SS7 Vl:SS8
Velocity                                      
SubjectSubject10                              
SubjectSubject11                              
SubjectSubject13                              
SubjectSubject2                               
SubjectSubject3                               
SubjectSubject6                               
SubjectSubject7                               
SubjectSubject8                               
SubjectSubject9                               
Velocity:SubjectSubject10                     
Velocity:SubjectSubject11                     
Velocity:SubjectSubject13                     
Velocity:SubjectSubject2                      
Velocity:SubjectSubject3                      
Velocity:SubjectSubject6                      
Velocity:SubjectSubject7   0.450              
Velocity:SubjectSubject8   0.398  0.828       
Velocity:SubjectSubject9   0.326  0.679  0.600

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.47194581 -0.46509026 -0.05537193  0.39069634  1.89436646 

Number of Observations: 40
Number of Groups: 10 
Warning message:
In pt(q, df, lower.tail, log.p) : NaNs produced
> 

Ahora, no entiendo de dónde puedo obtener el R ^ 2 y los valores p correspondientes que me indican si existe una relación lineal entre las dos variables o no, ni he entendido cómo mis datos se pueden trazar con la línea que se ajusta regresión.

¿Alguien puede ser tan amable de iluminarme? Realmente necesito su ayuda chicos ...

L_T
fuente
"Modelos de efectos mixtos y extensiones en ecología con R" por Zuur et al. es una buena introducción a los modelos lineales de efectos mixtos, que se centra menos en la teoría y más en la aplicación de la metodología.
Roland
Estimado Roland, creo que ese libro es útil, pero estoy buscando algo en línea ... ¿tiene alguna página web que sugerir?
L_T
1
Como dije en tu publicación anterior, lm () tiene una trama asociada. Entonces, si su modelo es M1, puede usar plot (M1).
Peter Flom - Restablece a Monica
querido @PeterFlom sí, pero también me dijiste que evitara usar lm para el diseño de medidas repetidas. Entonces, mi pregunta es si tengo que usar lm para analizar mis datos u otra función. ¿Cualquier sugerencia?
L_T
1
Como dije, busca modelos de varios niveles. En R, puedes mirar el nlmepaquete. Además, busque el tema en este sitio, aquí hay mucho escrito al respecto.
Peter Flom - Restablece a Monica

Respuestas:

17

Lo que haces realmente depende de los objetivos del análisis. No estoy seguro de cuáles son exactamente los objetivos de su análisis, pero analizaré varios ejemplos y espero que uno de ellos sea aplicable a su situación.

Caso 1 : una variable cuantitativa medida dos veces

Supongamos que realizó un estudio de sujetos humanos en el que hizo que los participantes realizaran una prueba de estadísticas dos veces y quisiera averiguar si los puntajes promedio en la segunda medición eran diferentes de la primera medición (para determinar si se produjo el aprendizaje). Si los puntajes test1 y test2 se almacenan en el marco de datos d, puede hacerlo completamente usando la función lm (), como en:

mod <- lm(test2 - test1 ~ 1, data = d)
summary(mod)

La prueba de la intersección es la prueba de la diferencia entre test1 y test2. Tenga en cuenta que no tendrá delta-R ^ 2 para la diferencia entre test1 y test2; en cambio, su medida del tamaño del efecto debería ser algo así como la d de Cohen.

Caso 2a : una variable cuantitativa medida dos veces, una variable dicotómica, medida totalmente entre sujetos

Digamos que tenemos el mismo diseño de estudio, pero queremos saber si ocurrieron diferentes tasas de aprendizaje para hombres y mujeres. Entonces, tenemos una variable cuantitativa (rendimiento de la prueba) que se mide dos veces, y una variable dicotómica, medida una vez. Suponiendo que test1, test2 y género están todos contenidos en el marco de datos d, también podríamos probar este modelo solo usando lm (), como en:

mod <- lm(test2 - test1 ~ gender, data = d)
summary(mod)
lm.sumSquares(mod) # lm.sumSquares() is located in the lmSupport package, and gives the change in R^2 due to the between-subjects part of the model

Suponiendo que el género está centrado (es decir, codificado, por ejemplo, masculino = -.5 y femenino = +.5), la intersección en este modelo es la prueba de la diferencia entre la prueba 1 y la prueba 2, promediada entre hombres y mujeres. El coeficiente de género es la interacción entre tiempo y género. Para obtener el efecto de género, promediado en el tiempo, tendría que hacer:

mod <- lm(rowMeans(cbind(test2, test1)) ~ gender, data = d)
summary(mod)

Caso 2b : una variable cuantitativa medida dos veces, una variable cuantitativa, solo medida una vez

Supongamos que nuevamente tenemos una variable cuantitativa medida dos veces y una variable cuantitativa medida una vez. Entonces, por ejemplo, digamos que teníamos una medida de interés de línea de base en las estadísticas y queríamos determinar si las personas que tenían niveles más altos de interés de línea de base aprendieron más de vez en cuando 2. Primero tendríamos que centrar el interés, como en :

d$interestc <- d$interest - mean(d$interest)

Suponiendo que test1, test2 e interésc están todos en el marco de datos d, esta pregunta podría entonces probarse de manera muy similar al Caso 1a:

mod <- lm(test2 - test1 ~ interestc, data = d)
summary(mod)
lm.sumSquares(mod)

Una vez más, la intercepción en este modelo evalúa si, en promedio a través del interés, los puntajes de las pruebas cambiaron del tiempo 1 al tiempo 2. Sin embargo, esta interpretación solo se cumple cuando el interés está centrado. El coeficiente de interés es si el efecto del tiempo depende del interés inicial. Podríamos obtener el efecto de interés, promediado en el tiempo, promediando juntos la prueba 1 y la prueba 2, como se indicó anteriormente, y probando el efecto de interés en esta variable compuesta.

Caso 2c : una variable cuantitativa medida dos veces, una variable categórica, solo medida una vez

Supongamos que su variable entre sujetos era una categoría, medida solo una vez. Entonces, por ejemplo, supongamos que estaba interesado en saber si las personas de diferentes razas (blancas vs asiáticas vs negras vs hispanas) tenían diferentes cantidades de aprendizaje de vez en cuando 2. Suponiendo que test1, test2 y raza están en el marco de datos d , primero deberías contrastar el código de carrera. Esto podría hacerse usando contrastes ortogonales planificados, códigos ficticios o códigos de efectos, dependiendo de las hipótesis / preguntas específicas que desee probar (recomiendo mirar lm.setContrasts () si está buscando una función auxiliar para hacer esto) . Suponiendo que la variable carrera ya esté codificada en contraste, usaría lm () de manera muy similar a los dos casos anteriores, como en:

mod <- lm(test2 - test1 ~ race, data = d)
summary(mod)
lm.sumSquares(mod)

Suponiendo que los contrastes de carrera están centrados, la intercepción en este modelo es, una vez más, el "efecto principal" del tiempo. Los coeficientes para los contrastes raciales son las interacciones entre esos contrastes y el tiempo. Para obtener los efectos generales de la raza, use el siguiente código:

Anova(mod, type = 3)

Caso 3 : una variable cuantitativa medida 3 veces (es decir, una manipulación dentro de los sujetos de tres niveles)

Supongamos que agregó un tercer punto de medición al diseño del caso uno. Entonces, sus participantes tomaron una prueba de estadísticas tres veces en lugar de dos. Aquí tiene un par de opciones, dependiendo de si desea una prueba general de las diferencias entre los puntos de tiempo (a veces no).

Por ejemplo, supongamos que su hipótesis principal es que los puntajes de las pruebas aumentarán linealmente de tiempo 1 a tiempo 3. Suponiendo que test1, test2 y test3 están en el marco de datos d, esta hipótesis podría probarse creando primero el siguiente compuesto:

d$lin <- d[, paste("test", sep = "", 1:3)] %*% c(-1, 0, 1)

Luego probaría si un modelo de solo intercepción que usa lin como variable dependiente tiene una intersección que es diferente de 0, como en:

mod <- lm(lin ~ 1, data = d)
summary(mod)

Esto le dará su prueba de si las puntuaciones de las estadísticas aumentaron con el tiempo. Por supuesto, puede crear otros tipos de puntajes de diferencia personalizados, según sus hipótesis particulares.

Si le interesan las pruebas de importancia omnibus, debe usar la función Anova () del paquete del automóvil. La implementación específica es un poco complicada. Básicamente, usted especifica qué variables están dentro de los sujetos y cuáles están entre sujetos usando lm (). Luego crea la parte del modelo dentro de los sujetos (es decir, especifica cuál de test1, test2 y test3 se midieron primero, segundo y tercero) y luego pasa ese modelo a Anova () creando un marco de datos llamado idata. Usando mi ejemplo hipotético:

mod <- lm(cbind(test1, test2, test3) ~ 1, data = d) # No between-subjects portion of the model
idata <- data.frame(time = c("time1", "time2", "time3")) # Specify the within-subjects portion of the model
mod.A <- Anova(mod, idata = idata, idesign = ~time, type = 3) # Gives multivariate tests.  For univariate tests, add multivariate = FALSE
summary(mod.A)

La declaración idesign le dice a Anova que incluya la variable de tiempo (compuesta de test1, test2 y test3) en el modelo. Este código le dará sus pruebas generales de los efectos del tiempo en los puntajes de las pruebas.

Caso 4 : una variable cuantitativa medida 3 veces, una variable cuantitativa entre sujetos

Este caso es una extensión simple del Caso 3. Como se indicó anteriormente, si simplemente le interesan las pruebas de 1 grado de libertad, simplemente puede crear un puntaje de diferencia personalizado con su variable dentro de las materias. Entonces, suponiendo que test1, test2, test3 e interés estén todos en el marco de datos d, y suponiendo que estamos interesados ​​en los efectos lineales del tiempo en los puntajes de las pruebas (y cómo esos efectos del tiempo varían con el interés de referencia), usted haría el seguimiento:

d$lin <- d[, paste("test", sep = "", 1:3)] %*% c(-1, 0, 1)

Luego, haga lo siguiente (con interés centrado en la media):

mod <- lm(lin ~ interestc, data = d)
summary(mod)
lm.sumSquares(mod)

Si desea pruebas omnibus, haga lo siguiente:

mod <- lm(cbind(test1, test2, test3) ~ interest, data = d) # We now have a between-subjects portion of the model
idata <- data.frame(time = c("time1", "time2", "time3"))
mod.A <- Anova(mod, idata = idata, idesign = ~time * interest, type = 3) # The idesign statement assumes that we're interested in the interaction between time and interest
summary(mod.A)

Otros casos: los omitiré por brevedad, pero son extensiones simples de lo que ya he descrito.

Tenga en cuenta que las pruebas omnibus (univariadas) de tiempo en las que el tiempo tiene más de 2 niveles, todas suponen esfericidad. Esta suposición se vuelve bastante insostenible a medida que aumenta el número de niveles. Si tiene bastantes puntos de medición en su diseño (por ejemplo, 4+), le recomiendo encarecidamente que use algo como el modelado multinivel y se mueva a un paquete especializado para esta técnica (como nlme o lme4 .

¡Espero que esto ayude!

Patrick S. Forscher
fuente
Estimado Patrick @ user1188407, muchas gracias por haber sido muy amable al dar esa respuesta. Lamentablemente, mi caso probablemente se ajusta a lo que escribió en las últimas oraciones ... por lo que necesitaría un ejemplo de código R para comprender cómo tratar mis datos. De hecho, si miras el diseño de mi experimento descrito en una publicación anterior stackoverflow.com/questions/12182373/… puedes ver que tengo una variable medida 4 veces (es decir, la velocidad medida en 4 condiciones)
L_T
y quiero encontrar si hay una relación lineal con una variable (respuesta_velocidad) que expresa la velocidad percibida en las cuatro condiciones. Entonces cada participante se sometió a 4 condiciones y luego evaluó la percepción de esas condiciones. Quiero saber si el rendimiento está relacionado con la percepción ...
L_T
Bueno, si desea utilizar una solución de modelado multinivel, puede recurrir a muchos recursos en línea diferentes. Para empezar, debería echar un vistazo al paquete nlme y esta viñeta . La viñeta está un poco desactualizada (2002), me pareció útil cuando estaba aprendiendo sobre el modelado multinivel. Finalmente, puede consultar el libro publicado por los creadores del paquete nlme.
Patrick S. Forscher
Estimado Patrick @ user1188407 gracias. Estudié los modelos multinivel y llegué a esta fórmula para analizar mis datos: lme (Velocity_response ~ Velocity * Subject, data = scrd, random = ~ 1 | Subject) ¿Pueden confirmarme que esta fórmula es correcta para el análisis? desea realizar en mis datos? Sin embargo, no entiendo cómo puedo obtener R ^ 2 y valores p, ni cómo trazar el gráfico con los puntos y la línea que se ajusta a la regresión. ¿Me podría ayudar? No soy un estático ...
L_T
La fórmula me parece correcta según mi comprensión de su estudio (y suponiendo que haya formateado sus datos en formato de período personal). Obtendría sus valores p guardando los resultados de su análisis en un objeto (como lo hago en mis ejemplos) y obteniendo un resumen de ese objeto. Sin embargo, debido a las diferencias entre los modelos multinivel y la regresión tradicional (p. Ej., En las métricas de tamaño del efecto, no existe una métrica estándar en los modelos multinivel), le sugiero que lea más sobre esta técnica antes de usarla. Parece que los otros usuarios han recomendado varias buenas opciones.
Patrick S. Forscher