Hosmer-Leme 展示 Python 中拟合检验的优点

Question

我估计了 python 中的 glm。我怎样才能表现出 Hosmer-Lemeshow 的优点

在 python 中对该模型进行拟合测试？

Answer 1

import matplotlib.pyplot as plt

import pandas as pd 
import numpy as np
from scipy.stats import chi2

# This could be made into a neat function of Hosmer-Lemeshow. 
def HosmerLemeshow (model,Y):
    pihat=model.predict()
    pihatcat=pd.cut(pihat, np.percentile(pihat,[0,25,50,75,100]),labels=False,include_lowest=True) #here we've chosen only 4 groups


    meanprobs =[0]*4 
    expevents =[0]*4
    obsevents =[0]*4 
    meanprobs2=[0]*4 
    expevents2=[0]*4
    obsevents2=[0]*4 

    for i in range(4):
       meanprobs[i]=np.mean(pihat[pihatcat==i])
       expevents[i]=np.sum(pihatcat==i)*np.array(meanprobs[i])
       obsevents[i]=np.sum(Y[pihatcat==i])
       meanprobs2[i]=np.mean(1-pihat[pihatcat==i])
       expevents2[i]=np.sum(pihatcat==i)*np.array(meanprobs2[i])
       obsevents2[i]=np.sum(1-Y[pihatcat==i]) 


    data1={'meanprobs':meanprobs,'meanprobs2':meanprobs2}
    data2={'expevents':expevents,'expevents2':expevents2}
    data3={'obsevents':obsevents,'obsevents2':obsevents2}
    m=pd.DataFrame(data1)
    e=pd.DataFrame(data2)
    o=pd.DataFrame(data3)
    
    # The statistic for the test, which follows, under the null hypothesis,
    # The chi-squared distribution with degrees of freedom equal to amount of groups - 2. Thus 4 - 2 = 2
    tt=sum(sum((np.array(o)-np.array(e))**2/np.array(e))) 
    pvalue=1-chi2.cdf(tt,2)

    return pd.DataFrame([[chi2.cdf(tt,2).round(2), pvalue.round(2)]],columns = ["Chi2", "p - value"])
    
HosmerLemeshow(glm_full,Y)

Answer 2

我找到了一种方法，代码质量不是最好的但是它有效：

import pandas as pd
import numpy as np
from scipy.stats import chi2
pihat=model.predict()
pihatcat=pd.cut(pihat, np.percentile(pihat,[0,25,50,75,100]),labels=False,include_lowest=True) #here I've chosen only 4 groups


meanprobs =[0]*4 
expevents =[0]*4
obsevents =[0]*4 
meanprobs2=[0]*4 
expevents2=[0]*4
obsevents2=[0]*4 

for i in range(4):
   meanprobs[i]=np.mean(pihat[pihatcat==i])
   expevents[i]=np.sum(pihatcat==i)*np.array(meanprobs[i])
   obsevents[i]=np.sum(data.r[pihatcat==i])
   meanprobs2[i]=np.mean(1-pihat[pihatcat==i])
   expevents2[i]=np.sum(pihatcat==i)*np.array(meanprobs2[i])
   obsevents2[i]=np.sum(1-data.r[pihatcat==i]) 


data1={'meanprobs':meanprobs,'meanprobs2':meanprobs2}
data2={'expevents':expevents,'expevents2':expevents2}
data3={'obsevents':obsevents,'obsevents2':obsevents2}
m=pd.DataFrame(data1)
e=pd.DataFrame(data2)
o=pd.DataFrame(data3)

tt=sum(sum((np.array(o)-np.array(e))**2/np.array(e))) #the statistic for the test, which follows,under the null hypothesis, the chi-squared distribution with degrees of freedom equal to amount of groups - 2 
pvalue=1-chi2.cdf(tt,2)
pvalue  

Answer 3

Hosmer-Lemeshow 是一个测试，仅当您的响应变量是二进制时才应用。我这么说是因为简单线性回归和泊松回归属于广义线性模型，但它们中的响应变量不是二元的。

我还没有在 Python 中找到任何函数来应用 Hosmer-Lemeshow 检验，所以我将使用传统方式进行此检验（计算检验统计量的值，然后计算 pvalue）。

作为测试统计我们将使用这个函数

地点：

O_{1g} 是第 g 个十分位数组中观察到的 Y = 1 个事件

O_{0g} 是第 g 个十分位数组中观察到的 Y = 0 事件

E_{1g} 是第 g 个十分位数组中预期的 Y = 1 个事件

E_{0g} 是第 g 个十分位数组中预期的 Y = 0 事件

另外，我们应该记住，这个检验统计量渐近地遵循 X^2_{(G - 2)}

来源

为了做这个测试我需要这些导入

import numpy as np

import pandas as pd

from scipy.stats import chi2 #In order to be able to calculate the pvalue

假设您的代码中有一个因变量 Y，它是一个 Series 对象，它包含我们观察到的二进制值（0 或 1）。另外，我们假设您的代码中有一个变量 X，它是一个数据框，包含模型的所有自变量。

#We will calculate the the predicted probabilities of success (Y = 1) for every data via
#the model we have created

#Let's suppose that you have named your model log_reg

predictions = log_reg(X)

#Now we will create a dataframe with two columns. The first column will represent the
#predicted probabilies and the second column will represent the binary data we observed

hl_df = pd.DataFrame({

"P_i": predictions,

"Outcome": Y

})

为了进行 Hosmer-Lemeshow 检验，我们必须根据样本分位数将保存预测概率的变量（在我们的例子中为 P_i）离散化为大小相等的桶。在 Hosmer-Lemeshow 测试中，我们通常创建 10 个分位数，但您可以使用例如 4 个分位数。这取决于你。

为此，我们将使用 pd.qcut() 函数。

hl_df["decile"] = pd.qcut(hl_df["P_i"],10)

现在，我们将计算$O_{0g}$和$O_{1g}$

#We will calculate all the observed ones in every decile

obsevents_1 = hl_df["Outcome"].groupby(hl_df.decile).sum()

#We will find all the observed zeroes of every decile if we substract the obsevents_1 from the
#number of elements in every decile

obsevents_0 = hl_df["Outcome"].groupby(hl_df.decile).count() - obsevents_1

现在我们来计算E_{0g} E_{1g}

我们知道变量 Y 是一个二进制变量，这意味着它服从 Bernoulli(p_i)。因此，变量 Y 的期望值为 p_i。此外，我们应该记住，Y 的每个值都相互独立。因此，在第 g 个十分位数中 Y = 1 的预期结果数为 p_1 + p_2 + p_3 + ... + p_n，其中 n 是属于第 g 个十分位数的 p_i 的数量。

expevents_1 = hl_df["P_i"].groupby(hl_df.decile).sum()

#We will find the expected number of events Y = 0 for every decile by substracting the
#expevents_1 from the number of elements in every decile

expevents_0 = hl_df["P_i"].groupby(hl_df.decile).count() - expevents_1

现在，我们准备计算检验统计量的值

hl = (((obsevents_0 - expevents_0)**2)/(expevents_0)).sum() + (((obsevents_1 - expevents_1)**2)/(expevents_1)).sum()

现在，我们将计算pvalue

pvalue = 1 - chi2.cdf(hl , 10 - 2)

#if you choose 4 quartiles instead of 10 deciles you should write pvalue = 1 - chi2.cdf(hl,4 - 2)

现在，根据您要应用的 p 值和显着性水平，您可以决定是否拒绝原假设。

我想提一下，对于上面的代码，我使用了this源。