Ich bin ein mit R vertrauterer Benutzer und habe versucht, zufällige Steigungen (Auswahlkoeffizienten) für ungefähr 35 Personen über 5 Jahre für vier Lebensraumvariablen zu schätzen. Die Antwortvariable gibt an, ob ein Standort als Lebensraum "verwendet" (1) oder "verfügbar" (0) war ("Verwendung" unten).
Ich verwende einen Windows 64-Bit-Computer.
In R-Version 3.1.0 verwende ich die folgenden Daten und Ausdrücke. PS, TH, RS und HW sind feste Effekte (standardisierte, gemessene Entfernung zu Lebensraumtypen). Ime4 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 gibt mir Parameterschätzungen für die für mich sinnvollen Fixeffekte, und die Zufallssteigungen (die ich als Auswahlkoeffizienten für jeden Lebensraumtyp interpretiere) sind auch sinnvoll, wenn ich die Daten qualitativ untersuche. Die Log-Wahrscheinlichkeit für das Modell beträgt -3050,8.
In den meisten Untersuchungen zur Tierökologie wird R jedoch nicht verwendet, da die räumliche Autokorrelation bei Tierstandortdaten die Standardfehler anfällig für Fehler vom Typ I machen kann. Während R modellbasierte Standardfehler verwendet, werden empirische (auch Huber-White- oder Sandwich-) Standardfehler bevorzugt.
Während R diese Option derzeit nicht anbietet (meines Wissens nach - bitte korrigieren Sie mich, wenn ich falsch liege), stimmt SAS zu - obwohl ich keinen Zugriff auf SAS habe, hat ein Kollege zugestimmt, dass ich mir seinen Computer ausleihen kann, um festzustellen, ob die Standardfehler vorliegen ändern sich signifikant, wenn die empirische Methode verwendet wird.
Zunächst wollten wir sicherstellen, dass SAS bei Verwendung modellbasierter Standardfehler ähnliche Schätzungen wie R erstellt, um sicherzustellen, dass das Modell in beiden Programmen auf dieselbe Weise angegeben wird. Es ist mir egal, ob sie genau gleich sind - nur ähnlich. Ich habe versucht (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;
Ich habe auch verschiedene andere Formen ausprobiert, beispielsweise das Hinzufügen von Linien
random intercept / subject = year(id) solution type=UN;
random intercept PS TH RS HW / subject = id solution type=UN;
Ich habe versucht, ohne die Angabe
solution type = UN,
oder auskommentieren
ddfm=betwithin;
Egal wie wir das Modell spezifizieren (und wir haben viele Möglichkeiten ausprobiert), ich kann nicht erreichen, dass die zufälligen Steigungen in SAS den Ausgaben von R aus der Ferne ähneln - obwohl die festen Effekte ähnlich genug sind. Und wenn ich anders meine, meine ich, dass nicht einmal die Zeichen gleich sind. Die -2-Log-Wahrscheinlichkeit in SAS betrug 71344,94.
Ich kann meinen vollständigen Datensatz nicht hochladen. Also habe ich einen Spielzeugdatensatz mit nur den Aufzeichnungen von drei Personen erstellt. SAS gibt mir die Ausgabe in wenigen Minuten; in R dauert es über eine Stunde. Seltsam. Mit diesem Spielzeug-Datensatz erhalte ich jetzt verschiedene Schätzungen für die festen Effekte.
Meine Frage: Kann jemand Aufschluss darüber geben, warum die Schätzungen der zufälligen Steigungen zwischen R und SAS so unterschiedlich sein können? Kann ich irgendetwas in R oder SAS tun, um meinen Code so zu ändern, dass die Aufrufe zu ähnlichen Ergebnissen führen? Ich möchte lieber den Code in SAS ändern, da ich glaube, dass mein R mehr schätzt.
Ich bin sehr besorgt über diese Unterschiede und möchte diesem Problem auf den Grund gehen!
Meine Ausgabe eines Spielzeugdatensatzes, der nur drei der 35 Personen im vollständigen Datensatz für R und SAS verwendet, ist als JPEG enthalten.
BEARBEITEN UND AKTUALISIEREN:
Wie @JakeWestfall herausgefunden hat, enthalten die Steigungen in SAS keine festen Effekte. Wenn ich die festen Effekte hinzufüge, ist hier das Ergebnis - Vergleichen von R-Steigungen mit SAS-Steigungen für einen festen Effekt, "PS", zwischen Programmen: (Auswahlkoeffizient = zufällige Steigung). Beachten Sie die erhöhte Variation in SAS.
ID
kein Faktor in R ist; Überprüfen Sie, ob sich dadurch etwas ändert.0
s und1
s gekennzeichnete BinomialdatenR
die Wahrscheinlichkeit einer "1" -Reaktion modellieren, während SAS die Wahrscheinlichkeit einer "0" -Reaktion modelliert. Um das SAS-Modell auf die Wahrscheinlichkeit "1" einzustellen, müssen Sie Ihre Antwortvariable als schreibenuse(event='1')
. Selbstverständlich sollten wir auch ohne dies meiner Meinung nach dieselben Schätzungen der zufälligen Effektvarianzen sowie dieselben Schätzungen der festen Effekte erwarten, wenn auch mit umgekehrten Vorzeichen.ranef()
Funktion verwenden und nichtcoef()
. Ersteres gibt die tatsächlichen Zufallseffekte an, während letzteres die Zufallseffekte plus den Vektor mit festen Effekten angibt. Das erklärt, warum die in Ihrem Beitrag abgebildeten Zahlen sehr unterschiedlich sind, aber es verbleibt immer noch eine erhebliche Diskrepanz, die ich nicht vollständig erklären kann.Antworten:
Es scheint , dass ich nicht die zufälligen Hängen zwischen Paketen ähnlich zu erwarten haben, nach Zhang et al 2011 In ihrem Papier auf Montage Generalized Linear Mixed-Effects - Modelle für binäre Antworten mit verschiedenen statistischen Pakete , die sie beschreiben:
Abstrakt:
Ich hoffe, @BenBolker und das Team betrachten meine Frage als eine Abstimmung für R, um empirische Standardfehler und die Gauß-Hermite-Quadratur-Fähigkeit für Modelle mit mehreren zufälligen Steigungsbegriffen zu integrieren, da ich die R-Schnittstelle bevorzuge und gerne anwenden könnte Einige weitere Analysen in diesem Programm. Erfreulicherweise sind die allgemeinen Trends gleich, auch wenn R und SAS keine vergleichbaren Werte für zufällige Steigungen aufweisen. Vielen Dank für Ihren Beitrag. Ich schätze die Zeit und die Überlegung, die Sie in diese Sache gesteckt haben, sehr!
quelle
Eine Mischung aus Antwort und Kommentar / mehr Fragen:
Ich habe den 'Spielzeug'-Datensatz mit drei verschiedenen Optimierungsoptionen ausgestattet. (* Anmerkung 1: Zu Vergleichszwecken wäre es wahrscheinlich sinnvoller, einen kleinen Datensatz durch Unterabtastung aus jedem Jahr und jeder ID zu erstellen, als durch Unterabtastung der Gruppierungsvariablen. Wir wissen, dass der GLMM keine Leistung erbringt besonders gut mit einer so kleinen Anzahl von Gruppen variabler Ebenen. Sie können dies über etwas tun wie:
Seriennummer:
Dann habe ich in einer neuen Sitzung die Ergebnisse eingelesen:
Verstrichene Zeit und Abweichung:
Diese Abweichungen liegen erheblich unter der Abweichung, die vom OP von R (6101,7) gemeldet wurde, und geringfügig unter der Abweichung, die vom OP von SAS (6078,9) gemeldet wurde, obwohl der Vergleich von Abweichungen zwischen Paketen nicht immer sinnvoll ist.
Ich war in der Tat überrascht, dass SAS in nur etwa 100 Funktionsbewertungen konvergierte!
Die Zeiten reichen von 17 Minuten (
nloptwrap
) bis 80 Minuten (bobyqa
) auf einem Macbook Pro, entsprechend der Erfahrung des OP. Abweichung ist ein bisschen besser fürnloptwrap
.Bei den Antworten sieht es ganz anders aus
nloptwrap
- obwohl die Standardfehler recht groß sind ...(Code hier gibt einige Warnungen über
year:id
dass ich aufspüren sollte)Fortsetzung folgt ... ?
quelle