Gamma GLM - Derivar intervalos de predicción para nuevos x_i

Question

En un MLG Gamma, se supone que el modelo estadístico para cada observación es $Y_i \sim Gamma(shape, scale)$ où $E(Y_i) = \mu_i = f(X_i\beta)$ y $f$ es la función de enlace.

He utilizado MLE para estimar $\hat{\beta}$ y $\hat{scale}$ y desea obtener un intervalo de predicción del 90% sobre un nuevo punto $Y'$ dado $X'$ .

Puedo producir los intervalos de confianza en $E(Y|X') = \mu'$ mediante la función de enlace $f$ sobre los intervalos de confianza de distribución normal para $X\hat{\beta}$ . Digamos que $\hat{\mu'} = 10$ y los intervalos de confianza del 90% son [5, 30].

Sin embargo, queremos los intervalos de la distribución de $Y'$ no $\mu'$ . Intuitivamente, estos intervalos deberían ser mucho más amplios que los intervalos de confianza para $\mu'$ I piense en también deben ser más amplios que los percentiles 5 y 95 de una distribución Gamma única con $\mu=\hat{\mu'}$ ya que la incertidumbre en torno a $\hat{\mu'}$ debería traducirse en un aumento de la incertidumbre en torno a la distribución final, algo así como una vaga a priori sobre una distribución posterior bayesiana.

¿Cuál es la forma correcta de modelar los intervalos de predicción en el nuevo punto $Y'$ ?

El siguiente esquema muestra cómo la incertidumbre sobre $\mu'$ se traduce en muchas posibles distribuciones gamma y un amplio intervalo de predicción para $Y'$

Referencias:

https://www.rocscience.com/help/swedge/swedge/Gamma_Distribution.htm

https://www.statsmodels.org/stable/glm.html

Answer 1

El intervalo de predicción para una nueva observación depende tanto de la aleatoriedad inherente asumida en este caso dada por la distribución gamma, como de la incertidumbre procedente de los parámetros que se estiman y que no se asumen como conocidos.

En general, no existe una expresión analítica o de forma cerrada para la combinación de ambos efectos. Las dos opciones principales son ignorar la incertidumbre de los parámetros y utilizar métodos de simulación.

Ignorar la incertidumbre de los parámetros: Si ignoramos que los parámetros se estiman con cierta incertidumbre de muestreo, la distribución de una nueva observación viene dada simplemente por la distribución gamma supuesta. Podemos utilizar las estimaciones de la media y la escala para calcular los intervalos de predicción pertinentes utilizando, por ejemplo, los métodos de distribución de scipy.stats. La parametrización podría tener que ser transformada de la parametrización GLM a la parametrización scipy.stats.

Métodos de simulación: Una posibilidad es utilizar bootstrap completo en la muestra original para simular nuevas observaciones. El método más sencillo consiste en suponer que la distribución normal asintótica para las estimaciones de los parámetros es adecuada y simular los parámetros de la distribución gamma a partir de la media y la covarianza de las estimaciones de los parámetros. Para cada parámetro muestreado podemos muestrear una nueva observación y calcular un intervalo de confianza a partir de las observaciones simuladas. Un problema con esto es que GLM sólo estima la media de los parámetros directamente, la escala se estima usando la desviación o pearson chi2. Es decir, GLM en statsmodels en otros paquetes no proporciona una covarianza conjunta para la media y el parámetro de escala.

Debido a estos problemas, statsmodels proporciona actualmente intervalos de predicción para nuevas observaciones que tienen en cuenta la incertidumbre de los parámetros sólo para el caso normal lineal, es decir, OLS.

Answer 2

Esta publicación responde exactamente a su pregunta

"Una derivación de los intervalos de predicción para la regresión gamma"

https://www.tandfonline.com/doi/full/10.1080/00949655.2016.1169421

Editar

El documento no está disponible públicamente, pero si le sirve de ayuda, aquí tiene un código R para su implementación. Digamos que desea modelar el resultado y que se distribuye gamma en dos predictores x1 y x2. Se requiere una predicción y un intervalo de predicción un nuevo punto x0=(x01,x02)=(0.75,-1).

Generar datos falsos

set.seed(1)
x1   = runif(100,-1,1); x2 = rnorm(100)
Mean = exp(.1+x1+x2/10) 
y    = rgamma(100,5,5/Mean)

Ajustar un modelo gamma y extraer la matriz del modelo

mod   = glm(y~x1+x2,family=Gamma(link="log"))
alpha = (mod$deviance/mod$df.residual)^-1  # 
X     = model.matrix(mod)
n     = nrow(X); p = ncol(X) 

Supongamos que la predicción es necesaria en x1 =0,75 y x2 = -1

x0    = c(1,0.75,-1)
pred  = exp(x0%*%coef(mod)) # Prediction of y at x0 

Encuentre los cuantiles de G_hat como se propone en la ec(3) del documento

h     = as.vector(t(x0)%*%solve(t(X)%*%X)%*%x0)
a1    = rchisq(10^6,2*alpha) # 10^6 is arbitrary large.
a2    = rt(10^6,n-p)
G     = a1/(a2*sqrt(h/alpha) + 1)
G_25  = quantile(G,0.025); G_975 = quantile(G,0.975)

Un intervalo de predicción del 95% en x0

c(G_25*pred/(2*alpha),G_975*pred/(2*alpha))

[1] 0.5745241 4.8665244