Modelli misti / Risultati REML: disponiamo di documentazione su come verrà calcolato l'errore standardUtilizzando SLS con REML, ovvero effetti misti, un cliente ottiene risultati diversi in pacchetti diversi (R / Python / JMP) per l'errore standard per l'intervallo di confidenza. (R / Python / JMP: 0,05 / 0,08 / 0,002) questo è un fattore di 40, che è piuttosto alto. Abbiamo documentazione su come viene calcolato esattamente (formula). Lo stesso BLUPS e le stime sia per FE che per ME sono il sema nei diversi pacchetti, ma l'intervallo di previsione e l'intervallo di confidenza.
Ho pensato che potesse avere qualcosa a che fare con una diversa struttura di covarianza ripetuta sottostante come predefinita o con la correzione di Kacker-Harville?
Aggiungendo un esempio: https://www.rdocumentation.org/packages/doBy/versions/4.6-3/topics/data-dietox
- Peso del modello di costruzione = Tempo + Mangime + Maiale (effetto casuale)
@jiancao @laural_jmp
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardCiao Martin,
Sei in grado di condividere i dati e lo script per il modello?Sarebbe di grande aiuto se tu potessi.Inoltre, quali errori standard e intervalli di confidenza sono diversi?Sto anche tentando di aggiungere Elizabeth a questa traccia.
Grazie.
Laura
@eclaassen
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardAh.Ora vedo la tabella dei dati di JMP.Non è arrivato nell'e-mail.Grazie.
Oggetto: Modelli misti / Risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardMartin,
Se è possibile fornire il codice R che produce gli intervalli di confidenza e gli errori standard, sarebbe utile. Sono in ritardo nella mia conoscenza di R e sono stato solo in grado di costruire il modello con successo (usando i dati dietox e lmer ()) e ottenere i valori previsti marginali e condizionali. Questi valori corrispondono a JMP.
Utilizzando i dati dietox, ho eseguito il modello specificato sia in JMP che in SAS (PROC MIXED), e i valori previsti marginali e condizionali e gli intervalli di confidenza corrispondono (quando il metodo DDFM corrispondente è specificato in SAS - KR2 per marginale e Satterthwaite per condizionale) . SAS non dispone di elementi della configurazione individuali. Si è verificato un problema con i gradi di libertà per gli elementi della configurazione individuali che è stato risolto nella versione 16.0.
So che la matrice di covarianza corretta di Kackar-Harville viene utilizzata per test / intervalli solo per stime di effetti fissi. Il metodo Satterthwaite viene utilizzato se è coinvolto un effetto casuale.
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardL'errore standard di 0,002 viene calcolato dall'utente dalle colonne dell'IC medio condizionale salvate.
Oggetto: Modelli misti / Risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardIn origine avevo aggiunto gli screeshot, ma non appariva nel post. Qui lo riprovo:
2.125 come Stderr per le componenti della varianza
Oggetto: Modelli misti / Risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardChiedo il codice R.
Oggetto: Modelli misti / Risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardPer quanto riguarda i gradi di libertà, abbiamo pensato proprio a questo, e se usano JMP 15, potrebbe essere questo il motivo del diverso risultato? Quando eseguo questo esempio sia in 16 che in 15 vedo esattamente gli stessi risultati.
Oggetto: Modelli misti / Risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardCiao Martin,
Credo che il metodo utilizzato per calcolare la matrice di covarianza delle componenti della varianza sia a pagina 1306 del documento seguente:
Wolfinger, R. D., Tobias, R. D. e Sall, J. (1994). "Calcolo delle verosimiglianze gaussiane e dei loro derivati per modelli misti lineari generali". SIAM Journal on Scientific Computing 15: 1294–1310.
La documentazione SAS per Proc Mixed fornisce qualche dettaglio in più se è utile?
Centro assistenza SAS: output predefinito
Elizabeth potrebbe sapere di un altro buon riferimento?
Laura
@eclaassen
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standard@jiancao ,
Non vedo questa colonna formula? Sarei molto interessato a vederlo, perché non so come lo potresti calcolare senza i gradi di libertà, che non sono semplici.
Secondo Martin stavano confrontando l'erro standard ...Secondo Martin stavano confrontando il file errori standard per l'intervallo di confidenza. Quindi, dopo aver salvato le colonne CI della media condizionale, prendo semplicemente la differenza e quindi uso la tabella per ottenere l'errore di supporto. Poiché ci sono diversi tipi di CI, stavo cercando di vedere quale CI stavano confrontando.
Ho anche eseguito SAS Proc Mixed ma i risultati sono leggermente diversi da quelli di JMP. Ad esempio, le stime dei componenti della varianza non sono esattamente le stesse.
Oggetto: Modelli misti / Risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardCiao Elisabeth,
allegato è ora anche il codice R. Loro hanno usato
checkintervals_dietox.R
* R 3.6.3
Qualcos'altro di cui hai bisogno? Sembra che abbiano usato anche lmer.
Ri: Secondo Martin stavano confrontando lo standard erro ...Uffa. Avevo appena scritto una lunga risposta e quando ho cliccato su post, ho ricevuto un errore di autenticazione e la risposta è andata persa. Proverò a ricreare ... e spero di non ricevere un post duplicato.
Quella quantità non è la SE? La differenza = obs + crit * se - (obs - crit * se) = 2 * crit * se. Prendere la varianza di questo non sarà l'errore standard a cui siamo interessati. Sarebbe la varianza dell'ampiezza degli intervalli di confidenza, che per quanto ne so non è niente di "interessante". Se ottieni l'errore standard dei valori previsti o dei limiti di confidenza, ottieni un SE ~ 0,834, che è più simile a un SE medio perché ogni osservazione ha il suo SE.
La discrepanza tra SAS e JMP è minore. Le probabilità (-2 log) sono le stesse, quindi ci imbattiamo in un po 'di ottimali locali o in qualche fuzz numerico a causa dei valori nel datatable non essendo ben arrotondati. Anche i valori previsti variano dopo la quarta cifra decimale nonostante gli effetti fissi siano gli stessi.
/* 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;
In realtà, qualsiasi discrepanza nel valore SE previsto tra JMP e Python / R deriva dalla differenza nell'SE per la varianza gruppo / maiale come indicato negli screenshot di Martin. La varianza della varianza entra in gioco. Non so come Python lo stia calcolando, e sembra strano che tutti gli altri SE siano d'accordo e quello no. Ancora una volta, SAS e JMP sostanzialmente concordano e SLS-REML e Mixed di JMP concordano mentre utilizzano routine leggermente diverse.
Oggetto: Modelli misti / Risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardGrazie, @martindemel . Darò un'occhiata.
Ri: Secondo Martin stavano confrontando lo standard erro ...Ah, il bit di informazione persa nel repost:
I SE di previsione condizionale sono disponibili in Mixed di JMPPro. Non so perché non siano disponibili in SLS-REML quando entrambe le piattaforme offrono l'SE marginale (errore standard di Predicted).
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardSembra che qualsiasi discrepanza negli intervalli di confidenza provenga dagli errori standard calcolati dei componenti della varianza.
L'output di Python ha la varianza group / pig come 31.836 con SE di 2.125. (Dallo screenshot di Martin)
L'output R è:
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 ha:
SAS ha:
| Stime dei parametri di covarianza |
|---|
| Cov Parm |
Stima |
Standard
Errore |
Valore Z |
Pr Z |
Alfa |
Inferiore |
Superiore |
|---|
| Maiale |
31.8276 |
5.4936 |
5.79 |
<.0001 |
0,05 |
21.0604 |
42.5949 |
|---|
| Residuo |
7.3645 |
0.3898 |
18.89 |
<.0001 |
0,05 |
6.6565 |
8.1924 |
|---|
Tutti i pacchetti ottengono la stessa verosimiglianza.
Le stime dei parametri sono essenzialmente le stesse in tutti i pacchetti.
L'SE per il componente varianza maiale / gruppo è paragonabile per SAS, JMP e R, ma molto diverso in Python.
L'SE per la varianza residua è paragonabile per SAS e JMP, ma R è molto diverso e non vedo alcun SE per quella stima nello screenshot dell'output di Python.
Poiché questi SE delle componenti della varianza sono necessari nei calcoli successivi, otteniamo differenze negli intervalli di confidenza della previsione. Non vedo un modo per ottenere i gradi di libertà e / o SE che R sta usando negli intervalli di confidenza per confrontarli direttamente. Vediamo solo le discrepanze risultanti negli intervalli.
Per vedere i gradi di libertà per vari intervalli, è possibile aggiungere test personalizzati in JMP. Questo script include test personalizzati equivalenti all'intervallo di confidenza medio e alla media condizionale CI. È possibile vedere il valore previsto e SE corrispondere a ciò che si ottiene salvando quelle colonne. Le informazioni aggiuntive sono il DFDen utilizzato per il test / CI. Quel DFDen può essere utilizzato per trovare il valore critico da utilizzare con SE per CI.
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" )
)}
)
);
Immagino che questo "spieghi" da dove proviene il problema, ma non spiega perché. Non so come R o Python stiano ottenendo gli SE dei componenti della varianza, ma c'è chiaramente una differenza lì. Le procedure JMP e SAS sono completamente documentate (credo che il documento JMP punti al documento SAS per tutte le informazioni di inferenza statistica profonda).Documento pertinente per la procedura Viya LMIXED (uguale al normale PROC MIXED)
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardHo notato che JMP e SAS lo etichettano come Standard Error, mentre il pacchetto R chiama Std.Dev. C'è una differenza? Nel suo rapporto ad effetto fisso, R etichetta come Std. Errore.
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardProbabilmente no. Sebbene tecnicamente non sia la stessa cosa, di solito i termini vengono usati in modo intercambiabile.
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardGrazie Elisabeth!
Per quanto riguarda l'algoritmo di ottimizzazione utilizzato, nel documento viene citato Newton-Raphson, ma con un'ozione da modificare in LMIXED. Immagino che in JMP usiamo solo NR e non abbiamo una modifica per cambiare, giusto?
Inoltre, non ci sarà alcuna opzione per scrivere la matrice dell'Assia. Avendo questo output potremmo essere in grado di confrontare se c'è stata una differenza tra l'iterazione e l'iterazione che differisce dall'algoritmo R, anche se non ce ne dovrebbe essere uno matematicamente. Il lavoro di confronto lo lascerei fare al cliente).
Grazie ancora per le tue indagini! È stato molto utile.
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardNon è possibile modificare il metodo di stima in JMP.
JMP non ha l'opzione per l'output della matrice hessiana (non la memorizziamo direttamente perché è grande, in genere), ma in JMP 16 in Fit Mixed ora ci sono opzioni per includere le matrici di covarianza asintotiche finali nel report.
Mi scuso per la risposta lenta. Nonostante l''iscrizione 'a questo argomento, non vengo informato molto rapidamente di nuove risposte. Se hai bisogno di una risposta più rapida su qualsiasi cosa, sentiti libero di taggarmi. Penso di essere avvisato dei tag molto più velocemente!
Oggetto: modelli misti / risultati REML: abbiamo la documentazione su come verrà calcolato l'errore standardElisabeth, non preoccuparti, hai risposto abbastanza velocemente:)
Grazie mille per le tue indagini. Questo aiuta molto a rispondere alle domande ea capire cosa è successo. Sebbene potremmo non essere in grado di scoprire completamente alcune differenze nei risultati tra R, Python e JMP, ora è molto più chiaro.
Ti auguro il meglio,
Martin
