Issues related to Gamma Regression-CodePudding

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.

CodePudding user response：

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

CodePudding user response：

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