Intervalle de prédiction pour une future proportion de succès dans le cadre binomial

9

Supposons que j'adapte une régression binomiale et que j'obtienne les estimations ponctuelles et la matrice de variance-covariance des coefficients de régression. Cela me permettra d'obtenir un IC pour la proportion attendue de succès dans une future expérience, , mais j'ai besoin d'un IC pour la proportion observée. Il y a eu quelques réponses connexes publiées, y compris la simulation (supposons que je ne veux pas faire ça) et un lien vers Krishnamoorthya et al (qui ne répond pas tout à fait à ma question).p

Mon raisonnement est le suivant: si nous utilisons uniquement le modèle binomial, nous sommes obligés de supposer que est échantillonné à partir de la distribution normale (avec le CI Wald correspondant) et qu'il est donc impossible d'obtenir un IC pour la proportion observée sous forme fermée. Si nous supposons que est échantillonné à partir de la distribution bêta, alors les choses sont beaucoup plus faciles car le nombre de succès suivra la distribution bêta-binomiale. Nous devrons supposer qu'il n'y a pas d'incertitude dans les paramètres bêta estimés, et .p αppαβ

Il y a trois questions:

1) Théorique: est-il acceptable d'utiliser uniquement les estimations ponctuelles des paramètres bêta? Je sais que pour construire un CI pour une future observation en régression linéaire multiple

Y=xβ+ϵ,ϵN(0,σ2)

ils font cette variance du terme d'erreur wrt, . Je suppose (corrigez-moi si je me trompe) que la justification est qu'en pratique est estimé avec une précision beaucoup plus grande que les coefficients de régression et nous ne gagnerons pas beaucoup en essayant d'incorporer l'incertitude de . Une justification similaire s'applique-t-elle aux paramètres bêta estimés, et ?σ 2 σ 2 α βσ2σ2σ2αβ

2) Quel est le meilleur package (R: gamlss-bb, betareg, aod?; J'ai également accès à SAS).

3) Compte tenu des paramètres bêta estimés, existe-t-il un raccourci (approximatif) pour obtenir les quantiles (2,5%, 97,5%) pour le nombre de succès futurs ou, mieux encore, pour la proportion de succès futurs sous distribution bêta-binomiale.

James
la source
À la première question, oui, c'est une chose valable que les gens font, cela s'appelle Empirical Bayes: en.wikipedia.org/wiki/Empirical_Bayes_method
Paul
1
Je ne pense pas que l'utilisation de la méthode XYZ pour estimer un paramètre de modèle puisse automatiquement impliquer qu'il est acceptable d'ignorer l'incertitude d'estimation lors de la production d'un IC pour une observation future. Par exemple, dans la régression linéaire multiple, ils utilisent OLS au lieu d'EB, et l'incertitude dans est également ignorée. Pourquoi donc? De plus, cet article du Wiki ne suggère jamais qu'en EB la précision d'estimation des hyperparamètres de haut niveau est généralement tellement plus élevée qu'il est acceptable de les considérer comme fixes pour des raisons pratiques. σ
James
1
«Lorsque la distribution réelle est fortement atteinte, l'intégrale déterminant peut ne pas être beaucoup modifiée en remplaçant la distribution de probabilité sur par une estimation ponctuelle représentant le pic de la distribution ». Que cela soit vrai dans votre cas dépend des spécificités de votre domaine problématique. p ( θ y ) η η p(ηy)p(θy)ηη
Paul
2
Bonne question! Vous ne pouvez pas obtenir de pivot, mais qu'en est-il de la probabilité de profil? Voir Quelles méthodes non bayésiennes existe-t-il pour l'inférence prédictive? .
Scortchi - Réintégrer Monica

Réponses:

1

Je vais aborder les 3 parties de la question.

Il y a deux problèmes confondus, le premier est la méthode que vous utilisez pour ajuster un modèle de régression dans ce cas. La seconde consiste à séparer les estimations de vos estimations pour prévoir une nouvelle estimation.

si vos variables de réponse sont distribuées binomialement, vous utiliserez généralement une régression logistique ou une régression probit (glm avec cdf normal comme fonction de lien).

Si vous effectuez une régression logistique, prenez la réponse comme étant le rapport des comptes observés divisé par la borne supérieure connue, c'est-à-dire . Prenez ensuite vos prédicteurs / covariables et mettez-les dans votre appel R à une fonction glm. L'objet retourné a tout ce dont vous avez besoin pour effectuer le reste de vos calculs. yi/ni

x<- rnorm(100, sd=2)
prob_true <- 1/(1+exp(-(1+5*x)))
counts <- rbinom(100, 50,prob_true)
print(d.AD <- data.frame(counts,x))
glm.D93 <- glm(counts/50 ~ x, family = binomial() )

Pour un modèle de régression linéaire , la formule d'un intervalle de prédiction est:

y^i±tnpsy1+1n+(xix¯)2(n1)sx2

Vous pouvez utiliser le modèle de régression linéaire comme approximation de la glm. Pour ce faire, vous devez utiliser une formule de régression linéaire pour la combinaison linéaire de prédicteurs avant d'effectuer la transformation de lien inverse pour récupérer les probabilités sur l'échelle 0-1. Le code pour ce faire est intégré dans la fonction Predict.glm () R. Voici un exemple de code qui fera également un joli tracé. ( EDIT : Ce code est pour l'intervalle de confiance, pas pour l'intervalle de prédiction)

y_hat <- predict(glm.D93, type="link", se.fit=TRUE)
t_np<- qt(.975, 100-2, ncp=0)

ub <- y_hat$fit + t_np * y_hat$se.fit
lb <- y_hat$fit - t_np * y_hat$se.fit

point <- y_hat$fit

p_hat <- glm.D93$family$linkinv(point)
p_hat_lb <- glm.D93$family$linkinv(lb)
p_hat_ub <- glm.D93$family$linkinv(ub)

plot(x,p_hat)
points(x, p_hat_ub, col='red')
points(x, p_hat_lb, col='blue')

Vous pouvez faire la même chose pour n'importe quel glm, par exemple Poisson, Gaussien inverse, gamma, etc. Dans chaque cas, faites l'intervalle de prédiction sur l'échelle de la combinaison linéaire des prédicteurs. Après avoir obtenu les deux points finaux de l'intervalle de prédiction, vous convertissez ces points finaux via le lien inverse. Pour chacun des glms que j'ai mentionnés, le lien inverse peut être différent du cas logit que j'ai écrit ici. J'espère que cela t'aides.

Lucas Roberts
la source