Como calcular o intervalo de predição em GLM (Gamma) / TweedieRegression em Python?
Eu verifiquei muitas fontes de sites sobre como conduzir o intervalo de previsão, especialmente na função GLM. Uma das abordagens é sobre intervalos de previsão para aprendizado de máquinahttps://machinelearningmastery.com/prediction-intervals-for-machine-learning/de Jason Brownlee. No entanto, seu método visa a regressão linear e pode não ser apropriado para o GLM (Gama) em alguns graus. Outra abordagem que descobri é usar o método de bootstrapping para conduzir o intervalo de predição. No entanto, o cálculo era muito demorado e a memória do meu computador foi interrompida ao executar a função do artigohttps://saattrupdan.github.io/2020-03-01-bootstrap-prediction/. Estou confuso sobre como conduzir o intervalo de predição de maneira apropriada em GLM (provavelmente Gamma) em Python em vez de em R. Eu encontrei um pacote relacionado em R, mas não quero usar R para conduzir o intervalo. Outra informação relacionada que encontrei na web é Gamma GLM - Derivar intervalos de predição para novos x_i: Gamma GLM - Derivar intervalos de predição para novos x_i .
Respostas
É um pouco complicado, mas deve ser viável.
Como diz aquele post, para obter um intervalo de predição, você deve integrar a incerteza nos coeficientes. Isso é difícil de fazer analiticamente, mas podemos simular isso. Aqui estão alguns dados de regressão gama
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()
Agora, vou ajustar a regressão gama a esses dados
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
Contanto que eu tenha dados suficientes, podemos fazer uma aproximação normal para a distribuição amostral dos coeficientes.
A média e a covariância podem ser obtidas no resumo do modelo.
beta_samp_mean = gamma_results.params
beta_samp_cov = gamma_results.cov_params()
dispersion = gamma_results.scale
Agora, é apenas uma questão de amostrar dados falsos usando essas estimativas e tomando quantis.
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)
É fácil traçar o intervalo de previsão
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ática do que está acontecendo
Nossos dados $y$ são distribuídos de acordo com
$$ y\vert X \sim \mbox{Gamma}(\phi, \mu(x)/\phi) $$
Pelo menos acho que é a parametrização correta do Gama, nunca consigo acertar. Em qualquer caso, presumindo que usamos um link de registro para o modelo, isso significa
$$ \mu(x) = \exp(X\beta)$$
A questão é que nunca sabemos $\beta$, nós só temos $\hat{\beta}$porque temos que estimar os parâmetros do modelo. Os parâmetros são, portanto, uma variável aleatória (porque dados diferentes podem produzir parâmetros diferentes). A teoria diz que com dados suficientes, podemos considerar
$$ \hat{\beta} \sim \mbox{Normal}(\beta, \Sigma) $$
e um pouco mais de teoria diz que ligando nossa estimativa para $\beta$ e $\Sigma$deve ser bom o suficiente. Deixei$\tilde{y}\vert X$ ser dados que posso ver para observações com covariáveis $X$. Se eu pudesse, eu realmente computaria
$$ \tilde{y} \vert X \sim \int p(y\vert X,\beta)p (\beta) \, d \beta $$
e então pegue quantis dessa distribuição. Mas essa integral é muito difícil, então, em vez disso, apenas a aproximamos simulando$p(\beta)$ (a distribuição normal) e passando tudo o que simulamos para $p(y\vert X, \beta)$ (neste caso, a distribuição gama).
Agora, eu sei que fui muito rápido e solto aqui, então se algum leitor quiser colocar um pouco mais de rigor na minha explicação, por favor, me avise em um comentário e eu irei limpar. Acho que isso deve ser bom o suficiente para dar ao OP uma ideia de como isso funciona.