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.
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.
ID
n'est pas un facteur dans R; Vérifiez et voyez si cela change quelque chose.0
s et1
s,R
modé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 formeuse(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.ranef()
fonction plutôt quecoef()
. 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.Réponses:
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:
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!
la source
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:
Code de montage par lots:
Ensuite, j'ai lu les résultats dans une nouvelle session:
Temps écoulé et déviance:
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 pournloptwrap
.Les réponses semblent assez différentes avec
nloptwrap
- bien que les erreurs standard soient assez importantes ...(le code ici donne des avertissements à
year:id
ce sujet que je devrais retrouver)À suivre ... ?
la source