Intervalo de predicción para una proporción futura de éxitos bajo configuración binomial

9

Supongamos que ajusto una regresión binomial y obtengo las estimaciones puntuales y la matriz de varianza-covarianza de los coeficientes de regresión. Eso me permitirá obtener un IC para la proporción esperada de éxitos en un experimento futuro, , pero necesito un IC para la proporción observada. Se han publicado algunas respuestas relacionadas, incluida la simulación (supongamos que no quiero hacer eso) y un enlace a Krishnamoorthya et al (que no responde a mi pregunta).p

Mi razonamiento es el siguiente: si usamos solo el modelo Binomial, nos vemos obligados a suponer que se muestrea de la distribución Normal (con el IC de Wald correspondiente) y, por lo tanto, es imposible obtener CI para la proporción observada en forma cerrada. Si suponemos que se muestrea a partir de la distribución beta, entonces las cosas son mucho más fáciles porque la cuenta de éxitos seguirá a la distribución beta-binomial. Tendremos que suponer que no hay incertidumbre en los parámetros beta estimados, y .ppαβ

Hay tres preguntas:

1) Teórico: ¿está bien usar solo las estimaciones puntuales de los parámetros beta? Sé que para construir un IC para la observación futura en regresión lineal múltiple

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

hacen esa varianza de término de error de wrt, . Supongo (corríjame si me equivoco) que la justificación es que en la práctica se estima con una precisión mucho mayor que los coeficientes de regresión y no ganaremos mucho al tratar de incorporar la incertidumbre de . ¿Se aplica una justificación similar a los parámetros beta estimados, y ?σ2σ2σ2αβ

2) ¿Qué paquete es mejor (R: gamlss-bb, betareg, aod ?; También tengo acceso a SAS).

3) Dados los parámetros beta estimados, ¿existe un atajo (aproximado) para obtener los cuantiles (2.5%, 97.5%) para el conteo de éxitos futuros o, mejor aún, para la proporción de éxitos futuros bajo la distribución Beta-Binomial.

James
fuente
En la pregunta uno, sí, esto es algo válido que la gente hace, se llama Empirical Bayes: en.wikipedia.org/wiki/Empirical_Bayes_method
Paul
1
No creo que usar el método XYZ para estimar un parámetro del modelo pueda implicar automáticamente que está bien ignorar la incertidumbre de la estimación al producir un IC para una observación futura. Por ejemplo, en la regresión lineal múltiple, utilizan OLS en lugar de EB, y la incertidumbre en se ignora. ¿Porqué es eso? Además, ese artículo de Wiki nunca sugiere que en EB la precisión de la estimación de los hiperparámetros de nivel superior es típicamente mucho mayor que está bien considerarlos fijos para fines prácticos. σ
James
1
“Cuando la distribución verdadera tiene un pico máximo, la determinación integral puede no cambiar mucho al reemplazar la distribución de probabilidad sobre con una estimación puntual representa el pico de la distribución ". Si eso es cierto en su caso depende de los detalles de su dominio problemático. p(ηy)p(θy)ηη
Paul
2
¡Buena pregunta! No puede obtener un pivote, pero ¿qué pasa con el uso de probabilidad de perfil? Vea ¿Qué métodos no bayesianos existen para la inferencia predictiva? .
Scortchi - Restablece a Monica

Respuestas:

1

Dirigiré las 3 partes a la pregunta.

Hay dos problemas combinados, primero es el método que utiliza para ajustar un modelo de regresión en este caso. El segundo es cómo hacer intervalos de estimaciones de sus estimaciones para predecir una nueva estimación.

si sus variables de respuesta están distribuidas binomialmente, normalmente usaría una regresión logística o una regresión probit (glm con cdf normal como función de enlace).

Si realiza una regresión logística, tome la respuesta como la razón de los recuentos observados dividida por el límite superior conocido, es decir, . Luego tome sus predictores / covariables y colóquelos en su llamada R a una función glm. El objeto devuelto tiene todo lo que necesita para hacer el resto de sus cálculos. 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() )

Para un modelo de regresión lineal , la fórmula para un intervalo de predicción es:

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

Puede usar el modelo de regresión lineal como una aproximación para el glm. Para hacer esto, utilizaría la fórmula de regresión lineal para la combinación lineal de predictores antes de realizar la transformación del enlace inverso para recuperar las probabilidades en la escala 0-1. El código para hacer esto está integrado en la función predict.glm () R. Aquí hay un código de ejemplo que también será un buen argumento. ( EDITAR : este código es para el intervalo de confianza, no para el intervalo de predicción)

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')

Puede hacer lo mismo para cualquier glm, por ejemplo, Poisson, Gaussiana inversa, gamma, etc. En cada caso, haga el intervalo de predicción en la escala de la combinación lineal de los predictores. Después de obtener los dos puntos finales del intervalo de predicción, convierte estos puntos finales a través del enlace inverso. Para cada uno de los glms que mencioné, el enlace inverso podría ser diferente al caso logit que escribí aquí. Espero que esto ayude.

Lucas Roberts
fuente