Conversión de fórmula de modelo mixto de medidas repetidas de SAS a R
mixed-models (2)
Por favor intente abajo:
model1 <- lme(
Y ~ GROUP + X1,
random = ~ GROUP | person,
correlation = corCompSymm(form = ~ day | person),
na.action = na.exclude, data = df1, method = "REML"
)
summary(model1)
Creo que al random = ~ groupvar | subjvar random = ~ groupvar | subjvar opción random = ~ groupvar | subjvar con R lme proporciona un resultado similar de la repeated / subject = subjvar group = groupvar con SAS/MIXED en este caso.
Editar:
SAS / MIXED
R (un modelo revisado2)
model2 <- lme(
Y ~ GROUP + X1,
random = list(person = pdDiag(form = ~ GROUP - 1)),
correlation = corCompSymm(form = ~ day | person),
weights = varIdent(form = ~ 1 | GROUP),
na.action = na.exclude, data = df1, method = "REML"
)
summary(model2)
Entonces, creo que estas estructuras de covarianza son muy similares (σ g1 = τ g 2 + σ 1 ).
Editar 2:
Estimaciones de covariables (SAS / MIXED):
Variance person GROUP TEST 8789.23
CS person GROUP TEST 125.79
Variance person GROUP CONTROL 82775
CS person GROUP CONTROL 33297
Asi que
TEST group diagonal element
= 125.79 + 8789.23
= 8915.02
CONTROL group diagonal element
= 33297 + 82775
= 116072
donde elemento diagonal = σ k1 + σ k 2 .
Estimaciones de covariables (R lme):
Random effects:
Formula: ~GROUP - 1 | person
Structure: Diagonal
GROUP1TEST GROUP2CONTROL Residual
StdDev: 14.56864 184.692 93.28885
Correlation Structure: Compound symmetry
Formula: ~day | person
Parameter estimate(s):
Rho
-0.009929987
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | GROUP
Parameter estimates:
1TEST 2CONTROL
1.000000 3.068837
Asi que
TEST group diagonal element
= 14.56864^2 + (3.068837^0.5 * 93.28885 * -0.009929987) + 93.28885^2
= 8913.432
CONTROL group diagonal element
= 184.692^2 + (3.068837^0.5 * 93.28885 * -0.009929987) + (3.068837 * 93.28885)^2
= 116070.5
donde elemento diagonal = τ g 2 + σ 1 + σ g 2 .
Hay varias preguntas y publicaciones sobre modelos mixtos para diseños experimentales más complejos, así que pensé que este modelo más simple ayudaría a otros principiantes en este proceso, así como a mí mismo.
Entonces, mi pregunta es: me gustaría formular un ancova de medidas repetidas en R a partir del procedimiento mixto sas proc:
proc mixed data=df1;
FitStatistics=akaike
class GROUP person day;
model Y = GROUP X1 / solution alpha=.1 cl;
repeated / type=cs subject=person group=GROUP;
lsmeans GROUP;
run;
Aquí está la salida SAS usando los datos creados en R (abajo):
. Effect panel Estimate Error DF t Value Pr > |t| Alpha Lower Upper
Intercept -9.8693 251.04 7 -0.04 0.9697 0.1 -485.49 465.75
panel 1 -247.17 112.86 7 -2.19 0.0647 0.1 -460.99 -33.3510
panel 2 0 . . . . . . .
X1 20.4125 10.0228 7 2.04 0.0811 0.1 1.4235 39.4016
A continuación se muestra cómo formulé el modelo en R usando el paquete ''nlme'', pero no obtengo estimaciones de coeficientes similares:
## create reproducible example fake panel data set:
set.seed(94); subject.id = abs(round(rnorm(10)*10000,0))
set.seed(99); sds = rnorm(10,15,5);means = 1:10*runif(10,7,13);trends = runif(10,0.5,2.5)
this = NULL; set.seed(98)
for(i in 1:10) { this = c(this,rnorm(6, mean = means[i], sd = sds[i])*trends[i]*1:6)}
set.seed(97)
that = sort(rep(rnorm(10,mean = 20, sd = 3),6))
df1 = data.frame(day = rep(1:6,10), GROUP = c(rep(''TEST'',30),rep(''CONTROL'',30)),
Y = this,
X1 = that,
person = sort(rep(subject.id,6)))
## use package nlme
require(nlme)
## run repeated measures mixed model using compound symmetry covariance structure:
summary(lme(Y ~ GROUP + X1, random = ~ +1 | person,
correlation=corCompSymm(form=~day|person), na.action = na.exclude,
data = df1,method=''REML''))
Ahora, la salida de R, que ahora me doy cuenta es similar a la salida de lm() :
Value Std.Error DF t-value p-value
(Intercept) -626.1622 527.9890 50 -1.1859379 0.2413
GROUPTEST -101.3647 156.2940 7 -0.6485518 0.5373
X1 47.0919 22.6698 7 2.0772934 0.0764
Creo que estoy cerca de la especificación, pero no estoy seguro de qué pieza me falta para que los resultados coincidan (dentro de lo razonable ...). ¡Cualquier ayuda sería apreciada!
ACTUALIZACIÓN: Usando el código en la respuesta a continuación, la salida R se convierte en:
> summary(model2)
Desplácese hacia abajo para ver las estimaciones de los parámetros: ¡mire! idéntico a SAS.
Linear mixed-effects model fit by REML
Data: df1
AIC BIC logLik
776.942 793.2864 -380.471
Random effects:
Formula: ~GROUP - 1 | person
Structure: Diagonal
GROUPCONTROL GROUPTEST Residual
StdDev: 184.692 14.56864 93.28885
Correlation Structure: Compound symmetry
Formula: ~day | person
Parameter estimate(s):
Rho
-0.009929987
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | GROUP
Parameter estimates:
TEST CONTROL
1.000000 3.068837
Fixed effects: Y ~ GROUP + X1
Value Std.Error DF t-value p-value
(Intercept) -9.8706 251.04678 50 -0.0393178 0.9688
GROUPTEST -247.1712 112.85945 7 -2.1900795 0.0647
X1 20.4126 10.02292 7 2.0365914 0.0811
Oooh, esto va a ser complicado, y si es posible con las funciones estándar de nlme, va a tomar un estudio serio de Pinheiro / Bates.
Sin embargo, antes de dedicar el tiempo a hacerlo, debe asegurarse de que este sea el modelo exacto que necesita. Quizás haya algo más que encaje mejor en la historia de tus datos. O tal vez haya algo que R pueda hacer más fácilmente que sea igual de bueno, pero no exactamente igual.
Primero, aquí está mi opinión sobre lo que estás haciendo en SAS con esta línea:
repeated / type=cs subject=person group=GROUP;
Este type=cs subject=person está induciendo la correlación entre todas las medidas en la misma persona, y esa correlación es la misma para todos los pares de días. El group=GROUP permite que la correlación para cada grupo sea diferente.
Por el contrario, aquí está mi opinión sobre lo que está haciendo su código R:
random = ~ +1 | person,
correlation=corCompSymm(form=~day|person)
Este código en realidad está agregando casi el mismo efecto de dos maneras diferentes; la línea random está agregando un efecto aleatorio para cada persona, y la línea de correlation está induciendo la correlación entre todas las medidas en la misma persona. Sin embargo, estas dos cosas son casi idénticas; si la correlación es positiva, obtienes exactamente el mismo resultado al incluir cualquiera de ellos. No estoy seguro de qué sucede cuando incluye ambos, pero sí sé que solo uno es necesario. De todos modos, este código tiene la misma correlación para todas las personas, no permite que cada grupo tenga su propia correlación.
Para que cada grupo tenga su propia correlación, creo que debes construir una estructura de correlación más complicada a partir de dos piezas diferentes; Nunca he hecho esto, pero estoy bastante seguro de recordar que Pinheiro / Bates lo está haciendo.
En su lugar, podría considerar agregar un efecto aleatorio para una persona y luego dejar que la varianza sea diferente para los diferentes grupos con weights=varIdent(form=~1|group) (desde la memoria, verifique mi sintaxis, por favor). Esto no será lo mismo pero cuenta una historia similar. La historia en SAS es que las mediciones en algunos individuos están más correlacionadas que las mediciones en otros individuos. Pensando en lo que eso significa, las mediciones para individuos con mayor correlación estarán más juntas que las mediciones para individuos con una menor correlación. En contraste, la historia en R es que la variabilidad de las mediciones dentro de los individuos varía; pensando en eso, las mediciones con mayor variabilidad tienen una correlación menor. Entonces sí cuentan historias similares, pero llegan desde lados opuestos.
Incluso es posible (pero me sorprendería) que estos dos modelos terminen siendo parametrizaciones diferentes de la misma cosa. Mi intuición es que la variabilidad general de la medición será diferente de alguna manera. Pero incluso si no son lo mismo, valdría la pena escribir las parametrizaciones solo para asegurarse de que las comprende y para asegurarse de que estén describiendo adecuadamente la historia de sus datos.