如何在statsmodels中使用gamma GLM的尺度和形状参数

Question

任务

我有这样的数据：

我想使用 statsmodels 将广义线性模型 (glm) 拟合到伽玛族中。使用此模型，对于我的每个观察结果，我想计算观察到小于（或等于）该值的值的概率。换句话说我要计算：

P(y <= y_i | x_i)

我的问题

如何从 statsmodels 中的拟合 glm 获取形状和比例参数？根据，statsmodels 中的比例参数未以正常方式进行参数化。我可以直接将它用作 scipy 中伽马分布的输入吗？还是需要先转型？
如何使用这些参数（形状和比例）来获得概率？目前我正在使用 scipy 为每个 x_i 生成分布并从中获取概率。请参阅下面的实施。

我目前的实现

import scipy.stats as stat
import patsy
import statsmodels.api as sm

# Generate data in correct form
y, X = patsy.dmatrices('y ~ x', data=myData, return_type='dataframe')

# Fit model with gamma family and log link
mod = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.log())).fit()

# Predict mean
myData['mu'] = mod.predict(exog=X) 

# Predict probabilities (note that for a gamma distribution mean = shape * scale)
probabilities = np.array(
    [stat.gamma(m_i/mod.scale, scale=mod.scale).cdf(y_i) for m_i, y_i in zip(myData['mu'], myData['y'])]
)

但是，当我执行此过程时，我得到以下结果：

目前看来预测的概率都很高。图中的红线是预测平均值。但即使对于低于这条线的点，预测的累积概率也约为 80%。这让我想知道我使用的比例参数是否确实是正确的。

Answer 1

在 R 中，您可以使用 1/dispersion 获得形状估计（检查此 post）。不幸的是，statsmodels 中的离散估计的命名是 scale。所以你确实取了它的倒数来得到形状估计。我用下面的例子来展示它：

values = gamma.rvs(2,scale=5,size=500)
fit = sm.GLM(values, np.repeat(1,500), family=sm.families.Gamma(sm.families.links.log())).fit()

这是一个仅截距模型，我们检查截距和色散（命名比例）：

[fit.params,fit.scale]
[array([2.27875973]), 0.563667465203953]

所以平均值是 exp(2.2599) = 9.582131 如果我们使用形状作为 1/dispersion ，shape = 1/0.563667465203953 = 1.774096 这就是我们模拟的。

如果我使用模拟数据集，它工作得很好。这是它的样子，形状为 10:

from scipy.stats import gamma
import numpy as np
import matplotlib.pyplot as plt
import patsy
import statsmodels.api as sm
import pandas as pd

_shape = 10
myData = pd.DataFrame({'x':np.random.uniform(0,10,size=500)})
myData['y'] = gamma.rvs(_shape,scale=np.exp(-myData['x']/3 + 0.5)/_shape,size=500)

myData.plot("x","y",kind="scatter")

然后我们像您一样拟合模型：

y, X = patsy.dmatrices('y ~ x', data=myData, return_type='dataframe')
mod = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.log())).fit()
mu = mod.predict(exog=X) 

shape_from_model = 1/mod.scale

probabilities = [gamma(shape_from_model, scale=m_i/shape_from_model).cdf(y_i) for m_i, y_i in zip(mu,myData['y'])]

和情节：

fig, ax = plt.subplots()
im = ax.scatter(myData["x"],myData["y"],c=probabilities)
im = ax.scatter(myData['x'],mu,c="r",s=1)
fig.colorbar(im, ax=ax)

如何在statsmodels中使用gamma GLM的尺度和形状参数

How to use scale and shape parameters of gamma GLM in statsmodels

python

statistics

regression

statsmodels