Issues related to Gamma Regression

Question

I want to derive coefficients of Gamma regression by iterated reweighted method (manually). When I run this code with out for{} loop it works properly but with loop it produce NaN. My code is:

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

}

Kindly help. My 2nd question is related to glm(). Is there any way which describe that how many iterations has glm() used? Regards.
Updates
with help of this I update this code but its not working.

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"))
library(gnlm)
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
)

I want to derive reg and reg1 same results. Kindly help.

Answer 1

For the first code chunk, the algorithm is for probit regression, not gamma. To perform the iterations manually using glm's default of no weights and no offset for family = Gamma(link = "inverse"), update the code as follows.

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 is the number of iterations performed by the glm function. Check that b is equal to the coefficients given by glm:

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

Answer 2

• Your inverse function is negative. Take away the minus sign.
• Also, change pshape to 1.0.
• I'm setting a seed for reproducibility.
• Initial values for small datasets is key. Setting them using glm results is a common approach if you can get a similar enough link. Another approach would be that in the answer by @jblood94. Yet another one would be to use nls() for (rough) initial estimates.
• argument trace=TRUE in glm() will show how many iterations

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)

Which gives:

> 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

An example of different link and using nls to get starting values:

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

And another

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