la conversión de medidas repetidas mezclados modelo fórmula de SAS a R

Question

Hay varias preguntas y puestos sobre los modelos mixtos para diseños experimentales más complejos, por lo que considera que la más sencillo modelo ayudaría a otros principiantes en este proceso, así como I.

por lo tanto, mi pregunta es me gustaría formular un ANCOVA de medidas repetidas en I sas procedimiento mixto 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 utilizando los datos creados en I (continuación):

.   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 la forma en que formulé el modelo de R utilizando el paquete 'nlme', pero no estoy haciendo estimaciones de los 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 se dan 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 más cerca a la especificación, pero no está seguro de qué pieza que me falta para que los resultados coinciden (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) 

de desplazamiento hacia abajo de las estimaciones de los parámetros - mira! 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 

Answer 1

Inténtelo más adelante:

model1 <- lme(
    Y ~ GROUP + X1, 
    random = ~ GROUP | person, 
    correlation = corCompSymm(form = ~ day | person), 
    na.action = na.exclude, data = df1, method = "REML" 
) 
summary(model1) 

creo random = ~ groupvar | subjvar opción con Rlme proporciona un resultado similar de repeated/subject = subjvar group = groupvar opción con SAS/MIXED en este caso.

Editar:

SAS/MIXTA

R (a model2 revisada)

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) 

Por lo tanto, creo que estos COVAR estructuras iance son muy similares (σ _g1 = τ _g + σ).

Editar 2: estimaciones

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 

Así

TEST group diagonal element 
    = 125.79 + 8789.23 
    = 8915.02 
CONTROL group diagonal element 
    = 33297 + 82775 
    = 116072 

donde elemento diagonal = σ _k1 + σ _k .

estimaciones 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 

Así

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 + σ + σ _g .

Answer 2

Oooh, esto va a ser complicado, y si es incluso posible con las funciones nlme estándar, va a tomar un estudio serio de Pinheiro/Bates.

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.

primer lugar, aquí es mi opinión sobre lo que está haciendo en el SAS con esta línea:

repeated/type=cs subject=person group=GROUP; 

Este type=cs subject=person está induciendo correlación entre todas las mediciones en la misma persona, y que la correlación es el mismo para todos pares de días El permite que la correlación para cada grupo sea diferente.

Por el contrario, aquí es mi opinión sobre lo que su código está haciendo R:

random = ~ +1 | person, 
correlation=corCompSymm(form=~day|person) 

Este código está añadiendo casi el mismo efecto de dos maneras diferentes; la línea random está agregando un efecto aleatorio para cada persona, y la línea correlation está induciendo una correlación entre todas las mediciones 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 debe 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 la 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.