Comment calculer l'intervalle de prédiction dans GLM (Gamma) / TweedieRegression en Python?
J'ai vérifié de nombreuses sources sur les sites Web sur la conduite de l'intervalle de prédiction, en particulier dans la fonction GLM. L'une des approches concerne les intervalles de prédiction pour l'apprentissage automatiquehttps://machinelearningmastery.com/prediction-intervals-for-machine-learning/de Jason Brownlee. Cependant, sa méthode cible la régression linéaire, et elle pourrait ne pas convenir à certains degrés au GLM (Gamma). Une autre approche que j'ai trouvée consiste à utiliser la méthode de bootstrap pour effectuer l'intervalle de prédiction. Cependant, le calcul prenait tellement de temps et la mémoire de mon ordinateur a été tuée lors de l'exécution de la fonction à partir de l'articlehttps://saattrupdan.github.io/2020-03-01-bootstrap-prediction/. Je ne sais pas comment conduire l'intervalle de prédiction de manière appropriée dans GLM (Gamma très probablement) en Python au lieu de R. J'ai trouvé un package associé dans R, mais je ne veux pas utiliser R pour conduire l'intervalle. Une autre information connexe que j'ai trouvée sur le Web est Gamma GLM - Dérivez les intervalles de prédiction pour new x_i: Gamma GLM - Dérivez les intervalles de prédiction pour le nouveau x_i .
Réponses
C'est un peu compliqué, mais cela devrait être faisable.
Comme le dit cet article, pour obtenir un intervalle de prédiction, vous devez intégrer l'incertitude dans les coefficients. C'est difficile à faire analytiquement, mais nous pouvons plutôt le simuler. Voici quelques données de régression 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()
Maintenant, je vais adapter la régression gamma à ces données
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
Tant que j'ai suffisamment de données, nous pouvons faire une approximation normale de la distribution d'échantillonnage des coefficients.
La moyenne et la covariance peuvent être obtenues à partir du résumé du modèle.
beta_samp_mean = gamma_results.params
beta_samp_cov = gamma_results.cov_params()
dispersion = gamma_results.scale
Maintenant, il s'agit simplement d'échantillonner de fausses données à l'aide de ces estimations et de prendre des quantiles.
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)
Il est alors facile de tracer l'intervalle de prédiction
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()
Mathématiques de ce qui se passe
Nos données $y$ sont distribués selon
$$ y\vert X \sim \mbox{Gamma}(\phi, \mu(x)/\phi) $$
Au moins je pense que c'est le paramétrage correct du Gamma, je ne peux jamais faire les choses correctement. Dans tous les cas, en supposant que nous utilisons un lien de journal pour le modèle, cela signifie
$$ \mu(x) = \exp(X\beta)$$
Le truc c'est qu'on ne sait jamais $\beta$, on obtient seulement $\hat{\beta}$car il faut estimer les paramètres du modèle. Les paramètres sont donc une variable aléatoire (car des données différentes peuvent produire des paramètres différents). La théorie dit qu'avec suffisamment de données, nous pouvons envisager
$$ \hat{\beta} \sim \mbox{Normal}(\beta, \Sigma) $$
et un peu plus de théorie dit que brancher notre estimation pour $\beta$ et $\Sigma$devrait être assez bon. Laisser$\tilde{y}\vert X$ être des données que je pourrais voir pour les observations avec des covariables $X$. Si je pouvais, je calculerais vraiment
$$ \tilde{y} \vert X \sim \int p(y\vert X,\beta)p (\beta) \, d \beta $$
puis prenez des quantiles de cette distribution. Mais cette intégrale est vraiment difficile, donc à la place, nous l'avons juste approximée en simulant à partir de$p(\beta)$ (la distribution normale) et en passant tout ce que nous avons simulé $p(y\vert X, \beta)$ (dans ce cas, la distribution gamma).
Maintenant, je me rends compte que j'ai été assez rapide et lâche ici, donc si des lecteurs veulent mettre un peu plus de rigueur dans mon explication, faites-le moi savoir dans un commentaire et je vais le nettoyer. Je pense que cela devrait être assez bon pour donner à OP une idée de la façon dont cela fonctionne.