Pourquoi SAS PROC GLIMMIX me donne-t-il des pentes aléatoires TRÈS différentes de glmer (lme4) pour un glmm binomial

12

Je connais mieux R et j'essaie d'estimer les pentes aléatoires (coefficients de sélection) pour environ 35 individus sur 5 ans pour quatre variables d'habitat. La variable de réponse est de savoir si un emplacement a été «utilisé» (1) ou «disponible» (0) de l'habitat («utilisation» ci-dessous).

J'utilise un ordinateur Windows 64 bits.

Dans R version 3.1.0, j'utilise les données et l'expression ci-dessous. PS, TH, RS et HW sont des effets fixes (distance mesurée normalisée aux types d'habitat). 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 donne des estimations de paramètres pour les effets fixes qui ont du sens pour moi, et les pentes aléatoires (que j'interprète comme des coefficients de sélection pour chaque type d'habitat) ont également un sens lorsque j'étudie les données qualitativement. La probabilité logarithmique du modèle est de -3050,8.

Cependant, la plupart des recherches en écologie animale n'utilisent pas R car avec les données de localisation des animaux, l'autocorrélation spatiale peut rendre les erreurs standard sujettes aux erreurs de type I. Alors que R utilise des erreurs standard basées sur un modèle, les erreurs standard empiriques (également Huber-white ou sandwich) sont préférées.

Bien que R n'offre pas actuellement cette option (à ma connaissance - S'IL VOUS PLAÎT, corrigez-moi si je me trompe), SAS le fait - bien que je n'ai pas accès à SAS, un collègue a accepté de me laisser emprunter son ordinateur pour déterminer si les erreurs standard changer de manière significative lorsque la méthode empirique est utilisée.

Premièrement, nous voulions nous assurer que lors de l'utilisation d'erreurs standard basées sur un modèle, SAS produirait des estimations similaires à R - pour être certain que le modèle est spécifié de la même manière dans les deux programmes. Je m'en fiche s'ils sont exactement les mêmes - juste similaires. J'ai essayé (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;

J'ai également essayé diverses autres formes, comme l'ajout de lignes

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

J'ai essayé sans préciser le

solution type = UN,

ou commentant

ddfm=betwithin;

Peu importe comment nous spécifions le modèle (et nous avons essayé de nombreuses façons), je ne peux pas obtenir les pentes aléatoires dans SAS pour ressembler à distance à celles de R - même si les effets fixes sont suffisamment similaires. Et quand je veux dire différent, je veux dire que même les signes ne sont pas les mêmes. La probabilité de -2 log dans SAS était de 71344,94.

Je ne peux pas télécharger mon ensemble de données complet; J'ai donc créé un jeu de données de jouets avec uniquement les enregistrements de trois personnes. SAS me donne une sortie en quelques minutes; dans R, cela prend plus d'une heure. Bizarre. Avec ce jeu de données de jouets, je reçois maintenant différentes estimations pour les effets fixes.

Ma question: quelqu'un peut-il nous expliquer pourquoi les estimations aléatoires des pentes peuvent être si différentes entre R et SAS? Puis-je faire quelque chose dans R ou SAS pour modifier mon code afin que les appels produisent des résultats similaires? Je préfère changer le code en SAS, car je "crois" que mes estimations R sont plus importantes.

Je suis vraiment préoccupé par ces différences et je veux aller au fond de ce problème!

Ma sortie d'un jeu de données jouet qui utilise seulement trois des 35 individus dans le jeu de données complet pour R et SAS est incluse en jpeg.

Sortie R Sortie SAS 1 Sortie SAS 2 Sortie SAS 3


MODIFIER ET METTRE À JOUR:

Comme @JakeWestfall a aidé à le découvrir, les pentes de SAS n'incluent pas les effets fixes. Lorsque j'ajoute les effets fixes, voici le résultat - en comparant les pentes R aux pentes SAS pour un effet fixe, "PS", entre les programmes: (coefficient de sélection = pente aléatoire). Notez la variation accrue de SAS.

R vs SAS pour PS

Nova
la source
Je remarque que ce IDn'est pas un facteur dans R; Vérifiez et voyez si cela change quelque chose.
Aaron a quitté Stack Overflow le
Je vois que vous ajustez les deux en utilisant l'approximation de Laplace pour la log-vraisemblance. Quels sont leurs scores respectifs de log-vraisemblance?
usεr11852
1
Avez-vous vérifié que vous modélisez la variable dépendante dans le même sens?
Peter Flom - Réintègre Monica
1
Soit dit en passant, ce à quoi Peter veut en venir, c'est que, par défaut avec les données binomiales étiquetées 0s et 1s, Rmodélisera la probabilité d'une réponse "1" tandis que SAS modélisera la probabilité d'une réponse "0". Pour rendre le modèle SAS la probabilité de "1", vous devez écrire votre variable de réponse sous la forme use(event='1'). Bien sûr, même sans cela, je pense que nous devrions toujours nous attendre aux mêmes estimations des variances à effet aléatoire, ainsi qu'aux mêmes estimations à effet fixe, bien que leurs signes soient inversés.
Jake Westfall du
1
@EricaN Une chose que vous venez de me rappeler est que vous devriez comparer les effets aléatoires de R à ceux de SAS en utilisant la ranef()fonction plutôt que coef(). Le premier donne les effets aléatoires réels, tandis que le second donne les effets aléatoires plus le vecteur à effets fixes. Cela explique donc en grande partie pourquoi les chiffres illustrés dans votre message diffèrent, mais il y a encore un écart substantiel que je ne peux pas totalement expliquer.
Jake Westfall

Réponses:

11

Il semble que je n'aurais pas dû m'attendre à ce que les pentes aléatoires soient similaires entre les paquets, selon Zhang et al 2011. Dans leur article sur l'ajustement des modèles à effets mixtes linéaires généralisés pour les réponses binaires utilisant différents packages statistiques , ils décrivent:

Abstrait:

Le modèle à effets mixtes linéaires généralisés (GLMM) est un paradigme populaire pour étendre les modèles de données transversales à un cadre longitudinal. Lorsqu'ils sont appliqués à la modélisation de réponses binaires, différents progiciels et même différentes procédures au sein d'un progiciel peuvent donner des résultats très différents. Dans ce rapport, nous décrivons les approches statistiques qui sous-tendent ces différentes procédures et discutons de leurs forces et faiblesses lorsqu'elles sont appliquées pour ajuster les réponses binaires corrélées. Nous illustrons ensuite ces considérations en appliquant ces procédures implémentées dans certains progiciels populaires à des données d'étude simulées et réelles. Nos résultats de simulation indiquent un manque de fiabilité pour la plupart des procédures considérées, ce qui entraîne des implications importantes pour l'application de tels progiciels populaires dans la pratique.

J'espère que @BenBolker et son équipe considéreront ma question comme un vote pour que R intègre les erreurs standard empiriques et la capacité de quadrature de Gauss-Hermite pour les modèles avec plusieurs termes de pente aléatoires à briller, car je préfère l'interface R et j'aimerais pouvoir appliquer quelques analyses supplémentaires dans ce programme. Heureusement, même si R et SAS n'ont pas de valeurs comparables pour les pentes aléatoires, les tendances générales sont les mêmes. Merci à tous pour votre contribution, j'apprécie vraiment le temps et la considération que vous y consacrez!

Nova
la source
désolé: qu'est-ce qu'une "erreur standard standard"? Voulez-vous dire les erreurs standard des composantes de la variance? Ou voulez-vous dire des erreurs standard sandwich?
Ben Bolker
désolé ... signifiait des SE empiriques / sandwich. J'ai modifié ma réponse.
Nova
@BenBolker Est-ce que cela a déjà été incorporé?
Lépidoptériste
Nan. Je continue d'essayer de comprendre comment je vais soutenir un développement comme celui-ci, car cela ne fait pas techniquement partie de mon programme de recherche ...
Ben Bolker
4

Un mélange de réponse et de commentaire / plus de questions:

J'ai équipé l'ensemble de données «jouet» de trois choix d'optimiseur différents. (* Note 1: il serait probablement plus utile à des fins de comparaison de faire un petit ensemble de données en sous-échantillonnant à l'intérieur de chaque année et id, plutôt qu'en sous-échantillonnant les variables de regroupement. En l'état, nous savons que le GLMM ne fonctionnera pas particulièrement bien avec un si petit nombre de niveaux variables de regroupement. Vous pouvez le faire via quelque chose comme:

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

Code de montage par lots:

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")
}

Ensuite, j'ai lu les résultats dans une nouvelle session:

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

Temps écoulé et déviance:

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

Ces écarts sont considérablement inférieurs à l'écart signalé par l'OP de R (6101.7), et légèrement inférieur à ceux signalés par l'OP de SAS (6078.9), bien que la comparaison des écarts entre les packages ne soit pas toujours judicieuse.

J'ai été en effet surpris que SAS n'ait convergé qu'en une centaine d'évaluations de fonctions!

Les durées varient de 17 minutes ( nloptwrap) à 80 minutes ( bobyqa) sur un Macbook Pro, conformément à l'expérience de l'OP. La déviance est un peu mieux pour 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

Les réponses semblent assez différentes avec nloptwrap- bien que les erreurs standard soient assez importantes ...

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

(le code ici donne des avertissements à year:idce sujet que je devrais retrouver)

À suivre ... ?

Ben Bolker
la source
Serait-il plus utile de vous envoyer l'ensemble de données complet? Le seul problème est que la convergence prend environ 9 heures avec l'ensemble de données complet, donc votre suggestion concernant l'échantillonnage est bonne. J'ai essayé de transformer les données avec une transformation logarithmique, mais le tracé résiduel groupé est toujours moche - pensez-vous que le tracé résiduel explique une partie du problème avec ces données? Enfin - vos résultats en SAS étaient-ils similaires à ceux de R?
Nova