对于第一个代码块，该算法用于概率单位回归，而不是γ。要使用glm的默认值（family = Gamma(link = "inverse")没有权重和偏移）手动执行迭代，请按如下方式更新代码。

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("(Intercept)" = 1,x1,x2,x3))
reg <- glm(y~x1+x2+x3, family = Gamma(link = "inverse"))

### step 1
eta <- 1/y

for(i in 1:reg$iter) {
  tX <- t(X <- x/eta)
  b <- drop(solve(tX%*%X)%*%tX%*%(2 - y*eta))
  eta <- drop(x %*% b)
}

reg$iter是glm函数执行的迭代次数。检查b是否等于glm给出的系数：

all.equal(reg$coefficients, b)
#> [1] TRUE

• 你的反函数是负的。去掉负号。
• 另外，将pshape更改为1.0。
• 我是在为可复制性埋下种子。
• 小数据集的初始值是关键。如果你能得到足够相似的链接，使用glm结果设置它们是一种常见的方法。另一种方法是@jblood94的答案。还有一种方法是使用nls()进行（粗略的）初始估计。
• glm()中的参数trace=TRUE将显示迭代次数

set.seed(111)
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"), trace=TRUE)
library(gnlm)
reg1<- gnlr(y=y,
     distribution = "gamma",
     mu = ~ inv(beta0 + beta1*x1 + beta2*x2 + beta3*x3),
     pmu = c(0.002, -0.002, -0.001, -0.001), ## or set to reg$coeff,
     pshape=1
)

cbind(c(reg$coeff,NA), reg1$coeff)

其中：

> cbind(c(reg$coeff,NA), reg1$coeff)
                     [,1]          [,2]
(Intercept)  0.0033899338  0.0033914440
x1          -0.0037481699 -0.0037476263
x2          -0.0007462714 -0.0007463346
x3          -0.0014941431 -0.0014936034
                       NA  2.8592334563

不同链接和使用nls获取起始值的示例：

nls.init3 <-
nls(y ~  beta0 + 1/(beta1+1)*x1 + sqrt(beta2)*x2 + beta3^2*x3,
    data=data.frame(y=y, x1=x1, x2=x2, x3=x3),
    start=list(beta0=1,beta1=.1,beta2=.1,beta3=.1)
    )
summary(nls.init3)$coefficients[,1]

reg3<- gnlr(y=y,
     distribution = "gamma",
     mu = ~  beta0 + 1/(beta1+1)*x1 + sqrt(beta2)*x2 + beta3^2*x3,
     pmu = summary(nls.init3)$coefficients[,1],
     pshape=1
)

reg3$coeff

还有一个

nls.init4 <-
nls(y ~  exp(beta0 + 1/(beta1+1)*x1),
    data=data.frame(y=y, x1=x1),
    start=list(beta0=0, beta1=0)
    )
summary(nls.init4)$coefficients[,1]

reg4<- gnlr(y=y,
     distribution = "gamma",
     mu = ~  exp(beta0 + 1/(beta1+1)*x1),
     pmu = summary(nls.init4)$coefficients[,1],
     pshape=1
)

reg4$coeff