¿Por qué SAS PROC GLIMMIX me da pendientes aleatorias MUY diferentes que glmer (lme4) para un binomial glmm

Soy un usuario más familiarizado con R, y he estado tratando de estimar pendientes aleatorias (coeficientes de selección) para aproximadamente 35 individuos durante 5 años para cuatro variables de hábitat. La variable de respuesta es si una ubicación fue "utilizada" (1) o "disponible" (0) hábitat ("uso" a continuación).

Estoy usando una computadora con Windows de 64 bits.

En R versión 3.1.0, uso los datos y la expresión a continuación. PS, TH, RS y HW son efectos fijos (distancia estandarizada y medida a los tipos de hábitat). lme4 V 1.1-7.

str(dat)
'data.frame':   359756 obs. of  7 variables:
 $ use     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Year    : Factor w/ 5 levels "1","2","3","4",..: 4 4 4 4 4 4 4 4 3 4 ...
 $ ID      : num  306 306 306 306 306 306 306 306 162 306 ...
 $ PS: num  -0.32 -0.317 -0.317 -0.318 -0.317 ...
 $ TH: num  -0.211 -0.211 -0.211 -0.213 -0.22 ...
 $ RS: num  -0.337 -0.337 -0.337 -0.337 -0.337 ...
 $ HW: num  -0.0258 -0.19 -0.19 -0.19 -0.4561 ...

glmer(use ~  PS + TH + RS + HW +
     (1 + PS + TH + RS + HW |ID/Year),
     family = binomial, data = dat, control=glmerControl(optimizer="bobyqa"))

glmer me da estimaciones de parámetros para los efectos fijos que tienen sentido para mí, y las pendientes aleatorias (que interpreto como coeficientes de selección para cada tipo de hábitat) también tienen sentido cuando investigo los datos cualitativamente. La probabilidad de registro para el modelo es -3050.8.

Sin embargo, la mayoría de las investigaciones en ecología animal no usan R porque con los datos de ubicación de los animales, la autocorrelación espacial puede hacer que los errores estándar sean propensos al error tipo I. Mientras que R usa errores estándar basados ​​en el modelo, se prefieren los errores estándar empíricos (también Huber-white o sandwich).

Si bien R no ofrece actualmente esta opción (que yo sepa, POR FAVOR, corríjame si estoy equivocado), SAS sí, aunque no tengo acceso a SAS, un colega accedió a prestarme su computadora para determinar si los errores estándar cambiar significativamente cuando se utiliza el método empírico.

Primero, deseamos asegurarnos de que cuando se usan errores estándar basados ​​en modelos, SAS produciría estimaciones similares a R, para estar seguros de que el modelo se especifica de la misma manera en ambos programas. No me importa si son exactamente iguales, simplemente similares. Intenté (SAS V 9.2):

proc glimmix data=dat method=laplace;
   class year id;
   model use =  PS TH RS HW / dist=bin solution ddfm=betwithin;
   random intercept PS TH RS HW / subject = year(id) solution type=UN;
run;title;

También probé varias otras formas, como agregar líneas

random intercept / subject = year(id) solution type=UN;
random intercept PS TH RS HW / subject = id solution type=UN;

Lo intenté sin especificar el

solution type = UN,

o comentando

ddfm=betwithin;

No importa cómo especifiquemos el modelo (y lo hemos intentado de muchas maneras), no puedo lograr que las pendientes aleatorias en SAS se parezcan remotamente a los resultados de R, a pesar de que los efectos fijos son lo suficientemente similares. Y cuando quiero decir diferente, quiero decir que ni siquiera los signos son iguales. La probabilidad de registro -2 en SAS fue 71344.94.

No puedo cargar mi conjunto de datos completo; así que hice un conjunto de datos de juguete con solo los registros de tres personas. SAS me da salida en unos minutos; en R lleva más de una hora. Extraño. Con este conjunto de datos de juguetes, ahora obtengo diferentes estimaciones para los efectos fijos.

Mi pregunta: ¿Alguien puede arrojar luz sobre por qué las estimaciones de pendientes aleatorias pueden ser tan diferentes entre R y SAS? ¿Hay algo que pueda hacer en R o SAS para modificar mi código para que las llamadas produzcan resultados similares? Prefiero cambiar el código en SAS, ya que "creo" que mi R estima más.

¡Estoy realmente preocupado por estas diferencias y quiero llegar al fondo de este problema!

Mi salida de un conjunto de datos de juguete que usa solo tres de los 35 individuos en el conjunto de datos completo para R y SAS se incluye como jpegs.

EDITAR Y ACTUALIZAR:

Como @JakeWestfall ayudó a descubrir, las pendientes en SAS no incluyen los efectos fijos. Cuando agrego los efectos fijos, aquí está el resultado: comparar las pendientes R con las pendientes SAS para un efecto fijo, "PS", entre programas: (Coeficiente de selección = pendiente aleatoria). Tenga en cuenta la mayor variación en SAS.

Parece que no debería haber esperado que las pendientes aleatorias fueran similares entre los paquetes, de acuerdo con Zhang et al 2011. En su artículo Sobre el ajuste de modelos lineales generalizados de efectos mixtos para respuestas binarias usando diferentes paquetes estadísticos , describen:

Resumen:

El modelo de efectos mixtos lineales generalizados (GLMM) es un paradigma popular para extender modelos de datos de sección transversal a un entorno longitudinal. Cuando se aplica al modelado de respuestas binarias, diferentes paquetes de software e incluso diferentes procedimientos dentro de un paquete pueden dar resultados bastante diferentes. En este informe, describimos los enfoques estadísticos que subyacen a estos diferentes procedimientos y discutimos sus fortalezas y debilidades cuando se aplican para ajustar las respuestas binarias correlacionadas. Luego ilustramos estas consideraciones aplicando estos procedimientos implementados en algunos paquetes de software populares a datos de estudio simulados y reales. Nuestros resultados de simulación indican una falta de confiabilidad para la mayoría de los procedimientos considerados, lo que conlleva implicaciones significativas para la aplicación de paquetes de software populares en la práctica.

Espero que @BenBolker y el equipo consideren mi pregunta como un voto para que R incorpore errores estándar empíricos y la capacidad de Cuadratura de Gauss-Hermite para modelos con varios términos de pendiente aleatorios para deslumbrar, ya que prefiero la interfaz R y me encantaría poder aplicar algunos análisis adicionales en ese programa. Afortunadamente, aunque R y SAS no tienen valores comparables para pendientes aleatorias, las tendencias generales son las mismas. Gracias a todos por su aporte, ¡realmente aprecio el tiempo y la consideración que pusieron en esto!

Una mezcla de una respuesta y comentario / más preguntas:

Puse el conjunto de datos 'juguete' con tres opciones diferentes de optimizador. (* Nota 1: probablemente sería más útil para fines comparativos hacer un pequeño conjunto de datos submuestreando dentro de cada año e id, en lugar de submuestrear las variables de agrupación. Tal como está, sabemos que el GLMM no funcionará particularmente bien con un número tan pequeño de niveles variables de agrupación. Puede hacerlo a través de algo como:

library(plyr)
subdata <- ddply(fulldata,c("year","id"),
    function(x) x[sample(nrow(x),size=round(nrow(x)*0.1)),])

Código de ajuste de lote:

Ntoy <- readRDS("Newton_toy.RDS")
library(lme4)
fitfun <- function(opt) {
    tt <- system.time(fit1 <- glmer(use ~  ps + th + rs + hw +
                                    (1 + ps + th + rs + hw |id/year),
                                    family = binomial, data = Ntoy,
                                    control=glmerControl(optimizer=opt),
                                    verbose=100))
    return(list(time=tt,fit=fit1))
}

opts <- c("nloptwrap","nlminbwrap","bobyqa")
## use for() instead of lapply so we can checkpoint more easily
res <- setNames(vector("list",length(opts)),opts)
for (i in opts) {
    res[[i]] <- fitfun(i)
    save("res",file="Newton_batch.RData")
}

Luego leí los resultados en una nueva sesión:

load("Newton_batch.RData")
library(lme4)

Tiempo transcurrido y desviación:

cbind(time=unname(sapply(res,function(x) x$time["elapsed"])),
          dev=sapply(res,function(x) deviance(x$fit)))
##                time      dev
## nloptwrap  1001.824 6067.706
## nlminbwrap 3495.671 6068.730
## bobyqa     4945.332 6068.731

Estas desviaciones están considerablemente por debajo de la desviación informada por el OP de R (6101.7), y ligeramente por debajo de las reportadas por el OP de SAS (6078.9), aunque comparar las desviaciones entre paquetes no siempre es sensato.

¡De hecho, me sorprendió que SAS convergiera en solo 100 evaluaciones de funciones!

Los tiempos varían de 17 minutos ( nloptwrap) a 80 minutos ( bobyqa) en un Macbook Pro, de acuerdo con la experiencia del OP. La desviación es un poco mejor para nloptwrap.

round(cbind(sapply(res,function(x) fixef(x$fit))),3)
##             nloptwrap nlminbwrap bobyqa
## (Intercept)    -5.815     -5.322 -5.322
## ps             -0.989      0.171  0.171
## th             -0.033     -1.342 -1.341
## rs              1.361     -0.140 -0.139
## hw             -2.100     -2.082 -2.082

Las respuestas parecen bastante diferentes nloptwrap, aunque los errores estándar son bastante grandes ...

round(coef(summary(res[[1]]$fit)),3)
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -5.815      0.750  -7.750    0.000
## ps            -0.989      1.275  -0.776    0.438
## th            -0.033      2.482  -0.013    0.989
## rs             1.361      2.799   0.486    0.627
## hw            -2.100      0.490  -4.283    0.000

(el código aquí da algunas advertencias sobre year:ideso que debería rastrear)

Continuará ... ?