Modèles MIxed / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeEn utilisant SLS avec REML, aka effets mixtes, un client obtient des résultats différents dans différents packages (R / Python / JMP) pour l'erreur standard de l'intervalle de confiance. (R / Python / JMP: 0,05 / 0,08 / 0,002) c'est un facteur de 40, ce qui est assez élevé. Avons-nous une documentation sur la façon dont il est calculé exactement (formule). Le BLUPS lui-même ainsi que les estimations pour FE et ME sont le sema dans les différents packages, mais l'intervalle de prédiction et l'intervalle de confiance.
Je pensais que cela pourrait avoir quelque chose à voir avec une structure de covariance répétée sous-jacente différente par défaut ou avec la correction de Kacker-Harville?
Ajout d'un exemple: https://www.rdocumentation.org/packages/doBy/versions/4.6-3/topics/data-dietox
- modèle de construction Poids = Temps + Alimentation + Cochon (effet aléatoire)
@jiancao @laural_jmp
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeSalut Martin,
Êtes-vous en mesure de partager les données et le script du modèle?Ce serait super utile si vous le pouvez.De plus, quelles erreurs standard et quels intervalles de confiance sont différents?J'essaye également d'ajouter Elizabeth à cette piste.
Merci.
Laura
@eclaassen
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeAh.Maintenant, je vois la table de données JMP.Il n'est pas arrivé sur l'e-mail.Merci.
Re: Modèles mixtes / résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeMartin,
Si le code R produisant les intervalles de confiance et les erreurs standard peut être fourni, cela serait utile. Je suis à la traîne dans mes connaissances R et je n'ai réussi qu'à construire le modèle (en utilisant les données Dietox et lmer ()) et à obtenir les valeurs prédites marginales et conditionnelles. Ces valeurs correspondent à JMP.
En utilisant les données de Dietox, j'ai exécuté le modèle spécifié à la fois dans JMP et SAS (PROC MIXED), et les valeurs prédites marginales et conditionnelles et les intervalles de confiance correspondent (lorsque la méthode DDFM correspondante est spécifiée dans SAS - KR2 pour marginal et Satterthwaite pour conditionnel) . SAS n'a pas de CI individuels. Il y avait un problème avec les degrés de liberté pour les CI individuels qui a été corrigé dans 16.0.
Je sais que la matrice de covariance corrigée de Kackar-Harville est utilisée pour les tests / intervalles pour les estimations à effets fixes uniquement. La méthode Satterthwaite est utilisée si un effet aléatoire est impliqué.
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeL'erreur standard de 0,002 est calculée par l'utilisateur à partir des colonnes CI moyenne conditionnelle enregistrées.
Re: Modèles mixtes / résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeJ'avais ajouté les captures d'écran à l'origine, mais cela n'apparaissait pas dans le message. Ici, je réessaye:
2.125 comme Stderr pour les composantes de variance
Re: Modèles mixtes / résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeJe demande le code R.
Re: Modèles mixtes / résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeConcernant les degrés de liberté, nous venons d'avoir cette pensée, et s'ils utilisent JMP 15, est-ce que cela pourrait être la raison du résultat différent? Lorsque j'exécute cet exemple en 16 et 15, je vois exactement les mêmes résultats.
Re: Modèles mixtes / résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeSalut Martin,
Je crois que la méthode utilisée pour calculer la matrice de covariance des composantes de la variance se trouve à la page 1306 de l'article suivant:
Wolfinger, R. D., Tobias, R. D. et Sall, J. (1994). "Calcul des probabilités gaussiennes et de leurs dérivés pour les modèles mixtes linéaires généraux." Journal du SIAM sur le calcul scientifique 15: 1294–1310.
La documentation SAS pour Proc Mixed donne un peu plus de détails si cela est utile?
Centre d'aide SAS: sortie par défaut
Elizabeth pourrait connaître une autre bonne référence?
Laura
@eclaassen
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculée@jiancao ,
Je ne vois pas cette colonne de formules? Je serais très intéressé à le voir, car je ne sais pas comment vous le calculeriez sans les degrés de liberté, qui ne sont pas simples.
Selon Martin, ils comparaient l'erro standard ...Selon Martin, ils comparaient les erreurs standard pour l'intervalle de confiance. Ainsi, après avoir enregistré les colonnes CI de moyenne conditionnelle, je prends simplement la différence, puis j'utilise le tableau pour obtenir l'erreur de support. Puisqu'il existe différents types de CI, j'essayais de voir quel CI ils comparaient.
J'ai également exécuté le SAS Proc Mixed mais les résultats sont légèrement différents de ceux de JMP. Par exemple, les estimations des composantes de la variance ne sont pas exactement les mêmes.
Re: Modèles mixtes / résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeSalut Elisabeth,
ci-joint est maintenant le code R. Ils ont utilisé
checkintervals_dietox.R
* R 3.6.3
Avez-vous besoin d'autre chose? Il semble qu'ils en aient aussi utilisé.
Re: Selon Martin, ils comparaient l'erro standard ...Pouah. Je venais juste d'écrire une longue réponse, et quand j'ai frappé post, j'ai eu une erreur d'authentification et la réponse a été perdue. J'essaierai de recréer ... et j'espère que nous n'obtiendrons pas de publication en double.
Cette quantité n'est pas la SE? La différence = obs + crit * se - (obs - crit * se) = 2 * crit * se. Prendre la variance de cela ne sera pas l'erreur standard qui nous intéresse. Ce serait la variance de la largeur des intervalles de confiance, ce qui, pour autant que je sache, n'a rien d'intéressant. Si vous obtenez l'erreur standard des valeurs prédites ou des limites de confiance, vous obtenez un SE ~ 0,834, ce qui ressemble plus à un SE moyen car chaque observation a son propre SE.
L'écart entre SAS et JMP est mineur. Les probabilités (-2 log) sont les mêmes, donc nous sommes confrontés à un peu d'optima local ou à un flou numérique en raison du fait que les valeurs de la table de données ne sont pas bien arrondies. Même les valeurs prédites varient après la 4e décimale, bien que les effets fixes soient les mêmes.
/* 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;
En réalité, tout écart dans la valeur prédite des SEs entre JMP et Python / R provient de la différence de SE pour la variance groupe / porc, comme indiqué dans les captures d'écran de Martin. La variance de la variance entre en jeu. Je ne sais pas comment Python le calcule, et il semble étrange que tous les autres SE soient d'accord et que ce ne soit pas le cas. Encore une fois, SAS et JMP sont essentiellement d'accord et les JMP SLS-REML et Mixed sont d'accord tout en utilisant des routines légèrement différentes.
Re: Modèles mixtes / résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeMerci, @martindemel . Je regarderai.
Re: Selon Martin, ils comparaient l'erro standard ...Ah, le peu d'informations perdues dans le repost:
Les SE de prédiction conditionnelle sont disponibles dans Mixte de JMPPro. Je ne sais pas pourquoi ils ne sont pas disponibles dans SLS-REML alors que les deux plates-formes offrent le SE marginal (Std Error of Predicted).
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeIl semble que tout écart dans les intervalles de confiance provienne des erreurs-types calculées des composantes de la variance.
La sortie Python a la variance groupe / pig comme 31.836 avec SE de 2.125. (D'après la capture d'écran de Martin)
La sortie R est:
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 a:
SAS a:
| Estimations des paramètres de covariance |
|---|
| Cov Parm |
Estimation |
la norme
Erreur |
Valeur Z |
Pr Z |
Alpha |
Inférieur |
Plus haut |
|---|
| Porc |
31.8276 |
5,4936 |
5,79 |
<.0001 |
0,05 |
21,0604 |
42,5949 |
|---|
| Résiduel |
7,3645 |
0,3898 |
18,89 |
<.0001 |
0,05 |
6,6565 |
8.1924 |
|---|
Tous les packages obtiennent la même probabilité logicielle.
Les estimations des paramètres sont essentiellement les mêmes dans tous les packages.
Le SE pour la composante de variance porc / groupe est comparable pour SAS, JMP et R, mais très différent en Python.
Le SE pour la variance résiduelle est comparable pour SAS et JMP, mais R est très différent et je ne vois pas de SE pour cette estimation dans la capture d'écran de sortie Python.
Puisque ces SE des composantes de la variance sont nécessaires dans les calculs ultérieurs, nous obtenons des différences dans les intervalles de confiance de prédiction. Je ne vois pas de moyen d'obtenir les degrés de liberté et / ou les SE que R utilise dans les intervalles de confiance pour les comparer directement. Nous ne voyons que les écarts résultants dans les intervalles.
Pour voir les degrés de liberté pour différents intervalles, des tests personnalisés peuvent être ajoutés dans JMP. Ce script comprend des tests personnalisés équivalents à l'intervalle de confiance moyen et à l'IC moyen conditionnel. Vous pouvez voir la valeur prédite et SE correspondent à ce qui est obtenu lors de l'enregistrement de ces colonnes. Les informations supplémentaires sont le DFDen utilisé pour le test / IC. Ce DFDen peut être utilisé pour trouver la valeur critique à utiliser avec le SE pour le 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" )
)}
)
);
Je suppose que cela «explique» d'où vient le problème, mais n'explique pas pourquoi. Je ne sais pas comment R ou Python obtiennent les SE des composants de variance, mais il y a clairement une différence. Les procédures JMP et SAS sont entièrement documentées (je pense que le document JMP pointe vers le document SAS pour toutes les informations d'inférence statistique approfondie).Document pertinent pour la procédure Viya LMIXED (identique à PROC MIXED classique)
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeJ'ai remarqué que JMP et SAS le qualifient d'erreur standard, tandis que le package R appelle Std.Dev. Y a-t-il une différence? Dans son rapport à effet fixe, R est étiqueté Std. Erreur.
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeProbablement pas. Bien que techniquement ce ne soit pas la même chose, les termes sont généralement utilisés de manière interchangeable.
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeMerci Elisabeth!
En ce qui concerne l'algorithme d'optimisation utilisé, il est mentionné dans la doc Newton-Raphson, mais avec une otion à changer dans LMIXED. Je suppose que dans JMP, nous utilisons simplement NR et n'avons pas de changement pour changer, correct?
Et il n'y aura pas non plus d'option pour écrire la matrice de Hesse. Ayant cette sortie, nous pourrions être en mesure de comparer s'il y a eu une différence entre l'itération et l'itération qui diffère de l'algorithme R, bien qu'il ne devrait pas y en avoir mathématiquement. Le travail de comparaison que je laisserais faire au client).
Merci encore pour vos investigations! Cela a été très utile.
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeIl n'y a pas d'option pour changer la méthode d'estimation dans JMP.
JMP n'a pas la possibilité de sortir la matrice de Hesse (nous ne la stockons pas directement car elle est grande, généralement), mais dans JMP 16 dans Fit Mixed, il existe maintenant des options pour inclure les matrices de covariance asymptotiques finales dans le rapport.
Je m'excuse pour la lenteur de la réponse. Malgré mon 'abonnement' à ce sujet, je ne suis pas informé très rapidement des nouvelles réponses. Si vous avez besoin d'une réponse plus rapide sur quelque chose, n'hésitez pas à me taguer. Je pense que je suis notifié des balises beaucoup plus rapidement!
Re: Modèles MIxés / Résultats REML: Avons-nous une documentation sur la façon dont l'erreur standard sera calculéeElisabeth, pas de soucis, tu as répondu assez vite:)
Merci beaucoup pour vos investigations. Cela aide beaucoup à répondre aux questions et à comprendre ce qui s'est passé. Bien que nous ne soyons pas en mesure de découvrir complètement certaines différences dans les résultats entre R, Python et JMP, c'est beaucoup plus clair maintenant.
Bonne chance,
Martin
