Aidez-moi à comprendre les valeurs de

13

J'essaie d'exécuter un logit bayésien sur les données ici . J'utilise bayesglm()dans le armpackage en R. Le codage est assez simple:

df = read.csv("http://dl.dropbox.com/u/1791181/bayesglm.csv", header=T)
library(arm)
model = bayesglm(PASS ~ SEX + HIGH, family=binomial(link="logit"), data=df)

summary(model) donne la sortie suivante:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.10381    0.10240   1.014    0.311    
SEXMale      0.02408    0.09363   0.257    0.797    
HIGH        -0.27503    0.03562  -7.721 1.15e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2658.2  on 1999  degrees of freedom
Residual deviance: 2594.3  on 2000  degrees of freedom
AIC: 2600.3

glm()bayesglm()pz

r bayesian p-value user3671
la source
C'est un commentaire et non une réponse, mais c'est ce qui aurait du sens pour moi. Vous obtenez des estimations qui sont probablement les valeurs pour lesquelles la distribution postérieure est maximisée. Il pourrait également être possible qu'ils ne soient que les moyens du postérieur? Cela vaut la peine de vérifier si vous le pouvez. Mais quels que soient les détails exacts, une fois que vous avez des estimations, vous pouvez les tester par l'estimation / Std habituelle. Erreur -> procédure z-score qui fonctionne si la partie postérieure est suffisamment proche d'une normale (elle passe à une normale dans certaines conditions qui tiennent généralement).
Erik
Erik ... Vous avez raison: les coefficients sont en effet les moyennes des densités postérieures. Ma question concerne les valeurs p et z. Que représentent-ils ici?
user3671
D'accord. Si vous avez une densité qui est approximativement normalement distribuée, vous pouvez tester si sa moyenne est nulle en prenant le score z = moyenne / écart-type de la distribution et en la comparant avec la distribution normale standard. Ensuite, vous regardez dans quelle mesure votre valeur ou une valeur supérieure serait peu probable sous la distribution normale standard -> valeur p. Voir z-score sur Wikipedia pour plus de détails.
Erik
Hé bien oui. Mais pourquoi s'embêter à faire ça dans un cadre bayésien? Dans l'inférence bayésienne, l'estimation ponctuelle est ma meilleure estimation du paramètre aléatoire, il n'est donc pas nécessaire de le tester. Tout au plus, je peux inclure un "intervalle crédible" qui équivaut à un "intervalle de confiance" fréquentiste mais dont l'interprétation statistique est très différente. C'est la partie déroutante de la sortie summary (). L'esprit est bayésien, mais la sortie est fréquentiste?
user3671
Un point est que l'estimation que vous obtiendrez sera différente, car vous avez utilisé un préalable. Et tandis que l'estimation ponctuelle est la "meilleure estimation" si vous voulez montrer de manière bayésienne que quelque chose a un effet, vous essayez de montrer que l'intervalle crédible ne contient pas le zéro. Lorsque vous approximez le postérieur par une normale avec la même moyenne et sd (asymptotiquement correcte), l'intervalle de crédibilité (1-p / 2) est le plus grand intervalle de crédibilité symétrique contenant le zéro, donc votre réponse est fondamentalement la même. Le p est la valeur de p indiquée ci-dessus.
Erik

Réponses:

16

Grande question! Bien qu'il existe des valeurs p bayésiennes , et l'un des auteurs du paquet arm est un défenseur, ce que vous voyez dans votre sortie n'est pas une valeur p bayésienne. Vérifiez la classe demodel

class(model)
"bayesglm" "glm"      "lm"

et vous pouvez voir que la classe bayesglm hérite de glm. En outre, l'examen de l'ensemble de bras ne montre aucune méthode de résumé spécifique pour un objet bayesglm. Alors quand tu fais

summary(model)

vous faites réellement 

summary.glm(model)

et obtenir une interprétation fréquentiste des résultats. Si vous voulez une perspective plus bayésienne, la fonction dans arm estdisplay()

atiretoo - réintégrer monica
la source
+1 Excellente réponse! C'est le problème avec R, il y a tellement de statisticiens très intelligents qui écrivent un code horrible qui laisse ces types de mines terrestres traîner.
Bogdanovist
Cela semble être un choix délibéré de la part des concepteurs plutôt qu'un oubli.
atiretoo - réintègre monica le
Après avoir lu le lien, je suis d'accord avec l'intention, mais dans ce cas, summary () aurait dû être réimplémenté pour simplement appeler display () plutôt que de donner des résultats absurdes sans avertissement. La personne qui a posé cette question a été déclenchée par un code qui a brisé le modèle utilisateur de R qui a été établi par tous les autres objets qu'ils ont jamais utilisés. C'est une terrible pratique de programmation.
Bogdanovist
2
Merci beaucoup, atiretoo. Cela soulève une autre question. Quelle est la différence entre display () et summary ()? Il me semble que la sortie de la première n'est que la sortie de la seconde, moins deux colonnes, et arrondie à 2 chiffres. Il semblerait que oui, d'après le post de Gelman que vous avez lié ci-dessus.
user3671
Oui, et d'après la discussion sur le blog d'Andrew Gelman, il semblerait qu'ils corrigent ce problème dans les futures versions du package de bras.
atiretoo - réintègre monica le