Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Bei Verwendung von SLS mit REML, auch als gemischte Effekte bezeichnet, erhält ein Kunde unterschiedliche Ergebnisse in verschiedenen Paketen (R / Python / JMP) für den Standardfehler für das Konfidenzintervall. (R / Python / JMP: 0,05 / 0,08 / 0,002) Dies ist ein Faktor von 40, was ziemlich hoch ist. Haben wir eine Dokumentation darüber, wie es genau berechnet wird (Formel)? Die BLUPS selbst sowie die Schätzungen für FE und ME sind das Sema in den verschiedenen Paketen, aber das Vorhersageintervall und das Konfidenzintervall.
Ich dachte, es könnte etwas mit einer anderen zugrunde liegenden wiederholten Kovarianzstruktur als Standard oder der Kacker-Harville-Korrektur zu tun haben?
Ein Beispiel hinzufügen: https://www.rdocumentation.org/packages/doBy/versions/4.6-3/topics/data-dietox
- Gebäudemodell Gewicht = Zeit + Futter + Schwein (Zufallseffekt)
@ Jiancao @laural_jmp
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Hallo Martin,
Können Sie die Daten und das Skript für das Modell freigeben?Das wäre super hilfreich, wenn du kannst.Welche Standardfehler und Konfidenzintervalle sind unterschiedlich?Ich versuche auch, Elizabeth zu diesem Track hinzuzufügen.
Vielen Dank.
Laura
@eclaassen
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Ah.Jetzt sehe ich die JMP-Datentabelle.Es kam nicht in der E-Mail durch.Vielen Dank.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Martin,
Wenn R-Code bereitgestellt werden kann, der die Konfidenzintervalle und Standardfehler erzeugt, wäre dies hilfreich. Ich bin in meinem R-Wissen zurückgeblieben und konnte das Modell nur erfolgreich erstellen (unter Verwendung der Dietox-Daten und lmer ()) und die marginalen und bedingten vorhergesagten Werte erhalten. Diese Werte stimmen mit JMP überein.
Unter Verwendung der Dietox-Daten habe ich das in JMP und SAS (PROC MIXED) angegebene Modell ausgeführt, und die vorhergesagten Grenz- und Bedingungswerte und Konfidenzintervalle stimmen überein (wenn die entsprechende DDFM-Methode in SAS - KR2 für Rand und Satterthwaite für Bedingt angegeben ist). . SAS hat keine individuellen CIs. Es gab ein Problem mit den Freiheitsgraden für die einzelnen CIs, das in 16.0 behoben wurde.
Ich weiß, dass die Kackar-Harville-korrigierte Kovarianzmatrix nur für Tests / Intervalle für Schätzungen mit festen Effekten verwendet wird. Die Satterthwaite-Methode wird verwendet, wenn ein zufälliger Effekt vorliegt.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Der Standardfehler von 0,002 wird vom Benutzer aus den gespeicherten Conditional Mean CI-Spalten berechnet.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Ich hatte die Screeshots ursprünglich hinzugefügt, aber sie erschienen nicht im Beitrag. Hier versuche ich es noch einmal:
2.125 als Stderr für die Varianzkomponenten
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Ich frage nach dem R-Code.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?In Bezug auf die Freiheitsgrade hatten wir gerade diesen Gedanken, und wenn sie JMP 15 verwenden, könnte dies der Grund für das unterschiedliche Ergebnis sein? Wenn ich dieses Beispiel in 16 und 15 ausführe, sehe ich genau die gleichen Ergebnisse.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Hallo Martin,
Ich glaube, die Methode zur Berechnung der Kovarianzmatrix der Varianzkomponenten befindet sich auf Seite 1306 des folgenden Dokuments:
Wolfinger, R. D., Tobias, R. D. und Sall, J. (1994). "Berechnung der Gaußschen Wahrscheinlichkeiten und ihrer Ableitungen für allgemeine lineare gemischte Modelle." SIAM Journal on Scientific Computing 15: 1294–1310.
Die SAS-Dokumentation für Proc Mixed enthält einige weitere Details, wenn dies hilfreich ist.
SAS-Hilfe: Standardausgabe
Elizabeth könnte eine andere gute Referenz kennen?
Laura
@eclaassen
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?@ Jiancao ,
Ich sehe diese Formelspalte nicht? Ich wäre sehr daran interessiert, es zu sehen, weil ich nicht weiß, wie Sie es ohne die Freiheitsgrade berechnen würden, die nicht einfach sind.
Laut Martin verglichen sie die Standardfehler ...Laut Martin verglichen sie die Standardfehler für das Konfidenzintervall. Nachdem ich die bedingten mittleren CI-Spalten gespeichert habe, nehme ich einfach die Differenz und verwende dann die Tabelle, um den Standfehler zu erhalten. Da es verschiedene Arten von CI gibt, habe ich versucht zu sehen, welche CI sie verglichen haben.
Ich habe auch SAS Proc Mixed ausgeführt, aber die Ergebnisse unterscheiden sich geringfügig von denen von JMP. Beispielsweise sind die Schätzungen der Varianzkomponenten nicht genau gleich.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Hallo Elisabeth,
Im Anhang befindet sich jetzt auch der R-Code. Sie benutzten
checkintervals_dietox.R
* R 3.6.3
Benötigen Sie noch etwas? Es scheint, dass sie auch lmer benutzt haben.
Re: Laut Martin haben sie den Standardfehler verglichen ...Pfui. Ich hatte gerade eine lange Antwort geschrieben, und als ich auf "Post" drückte, wurde ein Authentifizierungsfehler angezeigt, und die Antwort ging verloren. Ich werde versuchen, neu zu erstellen ... und hoffe, dass wir keinen doppelten Beitrag erhalten.
Diese Menge ist nicht die SE? Die Differenz = obs + krit * se - (obs - kritisch * se) = 2 * kritisch * se. Die Varianz davon zu nehmen, wird nicht der Standardfehler sein, an dem wir interessiert sind. Es wäre die Varianz der Breite der Konfidenzintervalle, die meines Wissens nichts 'Interessantes' ist. Wenn Sie den Standardfehler der vorhergesagten Werte oder der Konfidenzgrenzen erhalten, erhalten Sie eine SE ~ 0,834, was eher einer durchschnittlichen SE entspricht, da jede Beobachtung ihre eigene SE hat.
Die Diskrepanz zwischen SAS und JMP ist gering. Die (-2 log) Wahrscheinlichkeiten sind gleich, daher stoßen wir entweder auf ein bisschen lokales Optima oder auf einen numerischen Fuzz, da die Werte in der Datentabelle nicht gut gerundet sind. Sogar die vorhergesagten Werte variieren nach der 4. Dezimalstelle, obwohl die festen Effekte gleich sind.
/* outpm gives the marginal prediction information (values, se, CI) that corresponds to JMPs Mean Confidence Interval */
proc mixed data=sasuser.dietox nobound cl covtest;
class pig;
model weight= time feed / ddfm=kr2 s outpm=krpredmean;
random pig ;
run;
/* outp gives the conditional prediction information that corresponds to JMPs Conditional Mean CI */
proc mixed data=sasuser.dietox nobound;
class pig;
model weight= time feed / ddfm=satterth s outp=spred;
random pig;
run;
Tatsächlich ist jede Diskrepanz in den vorhergesagten SE-Werten zwischen JMP und Python / R auf den Unterschied in der SE für die Gruppen- / Schweinevarianz zurückzuführen, wie in Martins Screenshots angegeben. Die Varianz der Varianz kommt ins Spiel. Ich weiß nicht, wie Python es berechnet, und es scheint seltsam, dass alle anderen SEs zustimmen und einer nicht. Auch hier stimmen SAS und JMP im Wesentlichen überein, und die JMPs SLS-REML und Mixed stimmen überein, während leicht unterschiedliche Routinen verwendet werden.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Vielen Dank, @martindemel . Ich werde einen Blick darauf werfen.
Re: Laut Martin haben sie den Standardfehler verglichen ...Ah, die Informationen, die beim Repost verloren gegangen sind:
Die bedingten Vorhersage-SEs sind in JMPPro's Mixed verfügbar. Ich weiß nicht, warum sie in SLS-REML nicht verfügbar sind, wenn beide Plattformen die marginale SE (Std Error of Predicted) bieten.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Es sieht so aus, als ob eine Diskrepanz in den Konfidenzintervallen von den berechneten Standardfehlern der Varianzkomponenten herrührt.
Die Python-Ausgabe hat die Gruppen- / Schweinevarianz 31.836 mit einer SE von 2.125. (Aus Martins Screenshot)
Der R-Ausgang ist:
Random effects:
Groups Name Variance Std.Dev.
Pig (Intercept) 31.836 5.642
Residual 7.364 2.714
Number of obs: 789, groups: Pig, 72
JMP hat:
SAS hat:
| Schätzungen der Kovarianzparameter |
|---|
| Cov Parm |
Schätzen |
Standard
Error |
Z-Wert |
Pr Z. |
Alpha |
Niedriger |
Oberer, höher |
|---|
| Schwein |
31.8276 |
5.4936 |
5.79 |
<.0001 |
0,05 |
21.0604 |
42,5949 |
|---|
| Restwert |
7.3645 |
0,3898 |
18.89 |
<.0001 |
0,05 |
6,6565 |
8.1924 |
|---|
Alle Pakete erhalten die gleiche Loglikelihood.
Die Parameterschätzungen sind im Wesentlichen für alle Pakete gleich.
Die SE für die Schweine- / Gruppenvarianzkomponente ist für SAS, JMP und R vergleichbar, in Python jedoch sehr unterschiedlich.
Die SE für die Restvarianz ist für SAS und JMP vergleichbar, aber R ist sehr unterschiedlich und ich sehe keine SE für diese Schätzung im Python-Ausgabe-Screenshot.
Da diese SEs der Varianzkomponenten für die weiteren Berechnungen benötigt werden, erhalten wir Unterschiede in den Vorhersagekonfidenzintervallen. Ich sehe keinen Weg, um die Freiheitsgrade und / oder SEs zu erhalten, die R in den Konfidenzintervallen verwendet, um diese direkt zu vergleichen. Wir sehen die resultierenden Abweichungen nur in den Intervallen.
Um die Freiheitsgrade für verschiedene Intervalle anzuzeigen, können benutzerdefinierte Tests in JMP hinzugefügt werden. Dieses Skript enthält benutzerdefinierte Tests, die dem mittleren Konfidenzintervall und dem bedingten mittleren CI entsprechen. Sie können sehen, dass der vorhergesagte Wert und SE mit denen übereinstimmen, die beim Speichern dieser Spalten erhalten werden. Die zusätzlichen Informationen sind die DFDen, die für den Test / CI verwendet werden. Mit diesem DFDen kann der kritische Wert ermittelt werden, der mit der SE für das CI verwendet werden soll.
Fit Model(
Y( :Weight ),
Effects( :Feed, :Time ),
Random Effects( :Pig ),
NoBounds( 1 ),
Personality( "Standard Least Squares" ),
Method( "REML" ),
Emphasis( "Minimal Report" ),
Run(
:Weight << {Summary of Fit( 1 ), Analysis of Variance( 0 ),
Parameter Estimates( 1 ), Scaled Estimates( 0 ),
Plot Actual by Predicted( 0 ), Plot Regression( 0 ),
Plot Residual by Predicted( 0 ), Plot Studentized Residuals( 0 ),
Plot Effect Leverage( 1 ), Plot Residual by Normal Quantiles( 0 ),
Custom Test(
[1 5.200005 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0],
Label( "obs2 marginal" )
), Custom Test(
[1 5.200005 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0],
Label( "obs2 conditional" )
)}
)
);
Ich denke, dies "erklärt", woher das Problem kommt, erklärt aber nicht warum. Ich weiß nicht, wie R oder Python die SEs der Varianzkomponenten erhalten, aber da gibt es eindeutig einen Unterschied. JMP- und SAS-Prozeduren sind vollständig dokumentiert (ich denke, das JMP-Dokument verweist auf das SAS-Dokument für alle tiefen statistischen Inferenzinformationen).Relevantes Dokument für das Viya LMIXED-Verfahren (wie beim regulären PROC MIXED)
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Ich habe festgestellt, dass JMP und SAS es als Standardfehler kennzeichnen, während das R-Paket Std.Dev aufruft. Ist da ein Unterschied? In seinem Bericht mit festem Effekt bezeichnet R R als Std. Error.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Wahrscheinlich nicht. Obwohl dies technisch nicht dasselbe ist, werden die Begriffe normalerweise synonym verwendet.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Danke Elisabeth!
In Bezug auf den verwendeten Optimierungsalgorithmus wird im Dokument Newton-Raphson erwähnt, jedoch mit einer Änderung, die in LMIXED geändert werden muss. Ich denke in JMP verwenden wir nur NR und haben keine Änderung zu wechseln, richtig?
Und es wird auch keine Möglichkeit geben, die hessische Matrix auszuschreiben. Mit dieser Ausgabe können wir möglicherweise vergleichen, ob es einen Unterschied zwischen Itereation und Iteration gibt, der sich vom R-Algorithmus unterscheidet, obwohl es mathematisch keinen geben sollte. Die Vergleichsarbeit würde ich dem Kunden überlassen).
Nochmals vielen Dank für Ihre Untersuchungen! Das war sehr hilfreich.
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Es gibt keine Möglichkeit, die Schätzmethode in JMP zu ändern.
JMP hat nicht die Option, die hessische Matrix auszugeben (wir speichern sie nicht direkt, weil sie normalerweise groß ist), aber in JMP 16 in Fit Mixed gibt es jetzt Optionen, um die endgültigen asymptotischen Kovarianzmatrizen in den Bericht aufzunehmen.
Ich entschuldige mich für die langsame Antwort. Obwohl ich dieses Thema abonniert habe, werde ich nicht sehr schnell über neue Antworten informiert. Wenn Sie auf irgendetwas schneller eine Antwort benötigen, können Sie mich gerne markieren. Ich glaube, ich werde viel schneller über Tags benachrichtigt!
Betreff: Gemischte Modelle / REML-Ergebnisse: Haben wir eine Dokumentation, wie der Standardfehler berechnet wird?Elisabeth, keine Sorge, du hast schnell genug geantwortet:)
Vielen Dank für Ihre Untersuchungen. Dies hilft sehr, die Fragen zu beantworten und zu verstehen, was passiert ist. Obwohl wir einige Unterschiede in den Ergebnissen zwischen R, Python und JMP möglicherweise nicht vollständig aufdecken können, ist es jetzt viel klarer.
Alles Gute,
Martin
