如何手动计算Gamma族恒等链接GLM的回归系数？

你之前用的$\beta=(X^TX){-1}X^TY$是高斯族（正态）GLM的解析解，但Gamma族恒等链接的GLM没有解析解，必须通过**迭代加权最小二乘法（IRLS）**迭代求解，这也是R中glm()函数的底层逻辑。以下是具体实现步骤和代码：

核心原理 #

Gamma族恒等链接的模型中，因变量服从Gamma分布，其方差与均值的平方成正比（$\text{Var}(y)=\mu^2/\alpha$，$\alpha$为形状参数）。参数估计需要通过迭代构造加权伪数据，反复用加权最小二乘法更新系数，直到收敛：

1. 用初始值计算当前拟合值$\hat{\mu}=X\hat{\beta}$
2. 构造权重矩阵$W=\text{diag}(1/\hat{\mu}^2)$（权重与拟合值平方成反比）
3. 用加权最小二乘法更新系数：$\hat{\beta}=(X^TWX){-1}X^TWy$
4. 重复上述步骤，直到系数变化小于收敛阈值

手动实现R代码 #

# 假设x是自变量矩阵（无截距，对应公式~.+0），y是因变量向量
# 初始化系数（对应你设置的starter1，每列初始值100）
beta <- rep(100, ncol(x))
# 迭代参数设置
max_iter <- 100
tol <- 1e-8
iter <- 0
diff <- Inf

# IRLS迭代过程
while (iter < max_iter && diff > tol) {
  beta_old <- beta
  # 计算当前拟合值
  mu <- x %*% beta
  # 构造权重矩阵
  W <- diag(as.vector(1/(mu^2)))
  # 加权最小二乘更新系数
  xtw <- t(x) %*% W
  xtwx <- xtw %*% x
  xtwy <- xtw %*% y
  beta <- solve(xtwx, xtwy)
  # 计算系数变化量，判断是否收敛
  diff <- max(abs(beta - beta_old))
  iter <- iter + 1
}

# 输出最终系数，与glm(family=Gamma(link="identity"))结果一致
beta

关键注意事项 #

• Gamma族GLM的系数估计不存在解析解，必须依赖迭代，这是你之前用OLS公式得到高斯族结果的核心原因
• 初始值starter1会影响迭代的起始点，但只要模型收敛，最终结果会稳定到全局最优解
• 分散参数$\alpha$的估计不影响系数结果，若需要可以通过轮廓似然后续计算

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