R语言GLM与Python minimize方法的SE估计差异探究

问题核心：SE计算方法的差异 #

你在Python中用np.sqrt(np.diag(opt.hess_inv))计算的SE并不完全正确，原因在于：

• R的glm()拟合二项模型时，默认使用观测信息矩阵（Observed Information Matrix）的逆来计算标准误。观测信息矩阵是对数似然函数关于参数的二阶导数的负矩阵，在最优参数值处求值。
• 而scipy.optimize.minimize默认使用拟牛顿法（如BFGS）时，opt.hess_inv是算法迭代过程中对Hessian矩阵逆的近似值，并非严格的观测/期望信息矩阵逆，这会导致不同数据集下SE估计出现差异。

让R与Python SE结果一致的方法 #

要对齐两者的结果，需要在Python中严格计算与R一致的信息矩阵逆，具体步骤如下：

1. 确保对数似然函数与R的logit链接二项模型完全匹配 #

R的binomial()默认使用logit链接，对应的对数似然函数（负对数似然，用于minimize求最小值）应为：

import numpy as np

def neg_log_likelihood(params, X, y):
    beta0, beta1 = params
    eta = beta0 + beta1 * X  # 线性预测器，对应R的Age1模型（含截距）
    p = 1 / (1 + np.exp(-eta))  # logit链接转换为概率
    # 0/1二项变量的负对数似然
    neg_ll = -np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))
    return neg_ll

2. 手动计算观测信息矩阵并求逆 #

观测信息矩阵是负对数似然函数的二阶导数矩阵，直接对应R中用于计算SE的矩阵：

def observed_information_matrix(params, X, y):
    beta0, beta1 = params
    eta = beta0 + beta1 * X
    p = 1 / (1 + np.exp(-eta))
    w = p * (1 - p)  # 权重项
    # 构建信息矩阵
    I00 = np.sum(w)
    I01 = np.sum(w * X)
    I11 = np.sum(w * X**2)
    return np.array([[I00, I01], [I01, I11]])

拟合模型后用该矩阵计算SE：

from scipy.optimize import minimize

# 假设X是数据集的Age1列，y是PurchasedProb列（0/1变量）
X = dat['Age1'].values
y = dat['PurchasedProb'].values

# 初始参数猜测
initial_guess = [0, 0]

# 拟合模型
opt_result = minimize(neg_log_likelihood, initial_guess, args=(X, y), method='BFGS')

# 计算观测信息矩阵并求逆得到协方差矩阵
info_mat = observed_information_matrix(opt_result.x, X, y)
cov_matrix = np.linalg.inv(info_mat)
# 提取标准误
se = np.sqrt(np.diag(cov_matrix))

3. 让scipy直接计算精确Hessian #

如果不想手动编写信息矩阵，可以在minimize中指定使用精确Hessian计算，比如结合Newton-CG方法：

# 定义梯度（雅可比矩阵）
def gradient(params, X, y):
    beta0, beta1 = params
    eta = beta0 + beta1 * X
    p = 1 / (1 + np.exp(-eta))
    grad0 = np.sum(p - y)
    grad1 = np.sum((p - y) * X)
    return np.array([grad0, grad1])

# 用Newton-CG拟合，指定Hessian计算函数
opt_result = minimize(neg_log_likelihood, initial_guess, args=(X, y), 
                      method='Newton-CG', jac=gradient, 
                      hess=observed_information_matrix)

# 此时opt_result.hess_inv就是观测信息矩阵的逆，直接计算SE
se = np.sqrt(np.diag(opt_result.hess_inv))

4. 直接使用statsmodels匹配R的IRLS逻辑 #

R的glm()使用迭代重加权最小二乘（IRLS）拟合，Python的statsmodels库实现了完全一致的逻辑，无需手动编写似然函数：

import statsmodels.api as sm

# 添加截距项（R的~Age1默认包含截距）
X = sm.add_constant(dat['Age1'])
y = dat['PurchasedProb']

# 拟合二项GLM
model = sm.GLM(y, X, family=sm.families.Binomial())
result = model.fit()

# 直接获取与R一致的SE
se = result.bse

关键注意事项 #

• 确保数据集一致性：比如R和Python中变量类型、是否包含截距项（R默认添加截距，Python需手动用sm.add_constant处理）。
• 检查收敛性：若拟合不收敛，参数估计的细微差异会导致SE不同，需调整初始猜测或拟合方法。

内容的提问来源于stack exchange，提问作者Bogaso