Gamma回归迭代加权法NaN问题及glm迭代次数查询求助

1. 手动迭代加权法推导Gamma回归系数出现NaN #

问题描述 #

手动实现迭代加权法推导Gamma回归（inverse链接）系数时，无循环逻辑时运行正常，但加入循环后出现NaN值，原代码如下：

n<-10
y <- rgamma(n, 10, 0.1)
x1 <- rnorm(n, -1,1)
x2 <- rnorm(n, -1,1)
x3 <- rnorm(n, -1,1)

x<-as.matrix(cbind(1,x1,x2,x3))
reg <-glm(y~x1+x2+x3, family=Gamma(link = "inverse"))

### step 1
W<-G<-matrix(0,ncol=length(y),nrow=length(y))
b<-rep(0,4)

for(i in 1:50) {

### step 2
eta<-x%*%b
mu<-pnorm(eta)
diag(G)<-1/dnorm(eta)
z<-eta + G%*%(y - mu)
diag(W)<-(dnorm(eta)^2)/(mu*(1-mu))

### step 3
b <- solve(t(x)%*%W%*%x)%*%t(x)%*%W%*%z

}

问题原因与修正 #

原代码误用了逻辑回归的迭代公式（pnorm/dnorm/mu*(1-mu)），而非Gamma回归的正确推导逻辑。同时初始系数设为全0会导致eta=0，进而mu=1/eta无穷大，产生NaN。

Gamma回归（inverse链接）的正确迭代加权步骤：

• 链接函数：g(mu) = 1/mu = eta，即mu = 1/eta
• 工作变量z = eta + (y - mu)/g'(mu)，其中g'(mu) = -1/mu²
• 权重矩阵W的对角线元素：1/(Var(y|x)*(g'(mu))²)，Gamma分布的Var(y|x) = mu² * dispersion（dispersion为glm输出的离散参数）

修正后代码：

n<-10
y <- rgamma(n, 10, 0.1)
x1 <- rnorm(n, -1,1)
x2 <- rnorm(n, -1,1)
x3 <- rnorm(n, -1,1)

x<-as.matrix(cbind(1,x1,x2,x3))
reg <- glm(y~x1+x2+x3, family=Gamma(link = "inverse"))

### step 1
W <- matrix(0, ncol=n, nrow=n)
# 用glm的系数作为初始值，避免eta=0
b <- coef(reg) 

for(i in 1:50) {
  ### step 2
  eta <- x %*% b
  mu <- 1/eta  # inverse链接的均值计算
  # 计算工作变量z
  z <- eta + (y - mu) / (-1/mu^2) 
  # 计算权重矩阵对角线
  diag(W) <- mu^2 / reg$dispersion
  
  ### step 3
  b <- solve(t(x) %*% W %*% x) %*% t(x) %*% W %*% z
}

2. 查看glm()函数的迭代次数 #

直接提取glm模型对象的iter属性即可：

# 查看迭代次数
reg$iter

3. gnlr与glm拟合Gamma回归结果不一致 #

问题描述 #

使用gnlm包拟合Gamma回归时，无法得到与glm一致的结果，原代码如下：

library(gnlm)
# custom link / inverse 
inv <- function(eta)  -1/(eta)

n<-10
y <- rgamma(n, 10, 0.1)
x1 <- rnorm(n, -1,1)
x2 <- rnorm(n, -1,1)
x3 <- rnorm(n, -1,1)

x<-as.matrix(cbind(1,x1,x2,x3))
reg <-glm(y~x1+x2+x3, family=Gamma(link = "inverse"))
reg1<- gnlr(y=y,
     distribution = "gamma",
     mu = ~ inv(beta0 + beta1*x1 + beta2*x2 + beta3*x3),
     pmu = list(beta0=1, beta1=1, beta2=1, beta3=1),
     pshape=0.1
)

问题原因与修正 #

1. 链接函数错误：glm的inverse链接是mu = 1/eta，原代码的inv函数多了负号。
2. 参数化不一致：glm的Gamma分布离散参数dispersion = 1/shape，而gnlm需要传入的是shape参数，需转换。
3. 初始值不合理：初始值设为1可能导致收敛慢或结果偏差，建议用glm的系数作为初始值。

修正后代码：

library(gnlm)
# 正确的inverse链接函数
inv_link <- function(eta) 1/eta

n<-10
y <- rgamma(n, 10, 0.1)
x1 <- rnorm(n, -1,1)
x2 <- rnorm(n, -1,1)
x3 <- rnorm(n, -1,1)

reg <- glm(y~x1+x2+x3, family=Gamma(link = "inverse"))
# 转换glm的离散参数为gnlm需要的shape参数
shape_val <- 1/reg$dispersion

reg1<- gnlr(y=y,
     distribution = "gamma",
     mu = ~ inv_link(beta0 + beta1*x1 + beta2*x2 + beta3*x3),
     pmu = list(beta0=coef(reg)[1], beta1=coef(reg)[2], beta2=coef(reg)[3], beta3=coef(reg)[4]),
     pshape=shape_val
)

# 对比系数
coef(reg)
coef(reg1)

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