¿Cómo calcular el intervalo de predicción en GLM (Gamma) / TweedieRegression en Python?
He comprobado muchas fuentes de las web sobre cómo realizar el intervalo de predicción, especialmente en la función GLM. Uno de los enfoques es sobre intervalos de predicción para el aprendizaje automático.https://machinelearningmastery.com/prediction-intervals-for-machine-learning/de Jason Brownlee. Sin embargo, su método apunta a la regresión lineal, y podría no ser apropiado para GLM (Gamma) en algunos grados. Otro enfoque que encontré es utilizar el método de arranque para realizar el intervalo de predicción. Sin embargo, el cálculo llevó mucho tiempo y la memoria de mi computadora se mató al ejecutar la función del artículo.https://saattrupdan.github.io/2020-03-01-bootstrap-prediction/. Estoy confundido sobre cómo realizar el intervalo de predicción de una manera apropiada en GLM (Gamma muy probablemente) en Python en lugar de R.He encontrado un paquete relacionado en R, pero no quiero usar R para realizar el intervalo. Otra información relacionada que encontré en la web es Gamma GLM - Derivar intervalos de predicción para nuevos x_i: Gamma GLM - Derivar intervalos de predicción para nuevos x_i .
Respuestas
Es un poco complicado, pero debería ser factible.
Como dice esa publicación, para obtener un intervalo de predicción, debe integrar la incertidumbre en los coeficientes. Eso es difícil de hacer analíticamente, pero en cambio podemos simularlo. Aquí hay algunos datos de regresión gamma
N = 100
x = np.random.normal(size = N)
true_beta = np.array([0.3])
eta = 0.8 + x*true_beta
mu = np.exp(eta)
shape = 10
#parameterize gamma in terms of shaope and scale
y = gamma(a=shape, scale=mu/shape).rvs()
Ahora, ajustaré la regresión gamma a estos datos.
X = sm.tools.add_constant(x)
gamma_model = sm.GLM(y, X, family=sm.families.Gamma(link = sm.families.links.log()))
gamma_results = gamma_model.fit()
gamma_results.summary()
Generalized Linear Model Regression Results
Dep. Variable: ,y , No. Observations: , 100
Model: ,GLM , Df Residuals: , 98
Model Family: ,Gamma , Df Model: , 1
Link Function: ,log , Scale: ,0.075594
Method: ,IRLS , Log-Likelihood: , -96.426
Date: ,Mon, 30 Nov 2020, Deviance: , 7.7252
Time: ,22:45:07 , Pearson chi2: , 7.41
No. Iterations: ,7 , ,
Covariance Type:,nonrobust , ,
, coef , std err , z ,P>|z| , [0.025 , 0.975]
const, 0.8172, 0.028, 29.264, 0.000, 0.762, 0.872
x1 , 0.2392, 0.029, 8.333, 0.000, 0.183, 0.296
Siempre que tenga suficientes datos, podemos hacer una aproximación normal a la distribución muestral de los coeficientes.
La media y la covarianza se pueden obtener del resumen del modelo.
beta_samp_mean = gamma_results.params
beta_samp_cov = gamma_results.cov_params()
dispersion = gamma_results.scale
Ahora, es solo cuestión de muestrear datos falsos usando estas estimaciones y tomando cuantiles.
X_pred = np.linspace(-2, 2)
X_pred = sm.tools.add_constant(X_pred)
num_samps = 100_000
possible_coefficients = np.random.multivariate_normal(mean = beta_samp_mean, cov = beta_samp_cov, size = num_samps)
linear_predictions = [X_pred@b for b in possible_coefficients]
y_hyp = gamma(a=1/dispersion, scale = np.exp(linear_predictions)*dispersion).rvs()
# Here is the prediction interval
l, u = np.quantile(y_hyp, q=[0.025, 0.975], axis = 0)
Es fácil luego trazar el intervalo de predicción
yhat = gamma_results.predict(X_pred)
fig, ax = plt.subplots(dpi = 120)
plt.plot(X_pred[:,1], yhat, color = 'red', label = 'Estimated')
plt.plot(X_pred[:, 1], np.exp(0.8 + X_pred[:, 1]*true_beta), label = 'Truth')
plt.fill_between(X_pred[:, 1], l, u, color = 'red', alpha = 0.1, label = 'Prediction Interval')
for i in range(10):
y_tilde = gamma(a=shape, scale=np.exp(0.8 + X_pred[:, 1]*true_beta)/shape).rvs()
plt.scatter(X_pred[:, 1], y_tilde, s = 1, color = 'k')
plt.scatter(X_pred[:, 1], y_tilde, s = 1, color = 'k', label = 'New Data')
plt.legend()
Matemáticas de lo que está pasando
Nuestros datos $y$ se distribuyen según
$$ y\vert X \sim \mbox{Gamma}(\phi, \mu(x)/\phi) $$
Al menos creo que esa es la parametrización correcta del Gamma, nunca puedo hacerlo bien. En cualquier caso, asumiendo que usamos un enlace de registro para el modelo, esto significa
$$ \mu(x) = \exp(X\beta)$$
La cosa es que nunca lo sabemos $\beta$, solo obtenemos $\hat{\beta}$porque tenemos que estimar los parámetros del modelo. Por tanto, los parámetros son una variable aleatoria (porque diferentes datos pueden producir diferentes parámetros). La teoría dice que con suficientes datos, podemos considerar
$$ \hat{\beta} \sim \mbox{Normal}(\beta, \Sigma) $$
y algo más de teoría dice que, al incorporar nuestra estimación de $\beta$ y $\Sigma$debería ser lo suficientemente bueno. Dejar$\tilde{y}\vert X$ ser datos que podría ver para observaciones con covariables $X$. Si pudiera, realmente calcularía
$$ \tilde{y} \vert X \sim \int p(y\vert X,\beta)p (\beta) \, d \beta $$
y luego tomar cuantiles de esta distribución. Pero esta integral es realmente difícil, así que en su lugar simplemente la aproximamos simulando desde$p(\beta)$ (la distribución normal) y pasar todo lo que simulamos a $p(y\vert X, \beta)$ (en este caso, la distribución gamma).
Ahora, me doy cuenta de que he sido bastante rápido y relajado aquí, así que si algún lector quiere poner un poco más de rigor en mi explicación, hágamelo saber en un comentario y lo limpiaré. Creo que esto debería ser lo suficientemente bueno como para darle a OP una idea de cómo funciona.