Iterative re-weighted least squares(IRLS) for Poisson regression with identity link and canonical log link

 

enwul

2020-10-25

y<-c(2,3,6,7,8,9,10,12,15)
x<-matrix(c(1,1,1,1,1,1,1,1,1,-1, -1, 0, 0, 0, 0, 1, 1, 1), nrow=9,ncol=2)
b.init<-c(7,5)
x1<-c(-1, -1, 0, 0, 0, 0, 1, 1, 1)
data_4.3<-data.frame(y,x1)

poisson_IRLS<-function(x,y,b.init,tol=1e-8){
  change <- Inf
  b.old <- b.init
  while(change > tol) {
    eta <- x %*% b.old  # linear predictor
    w<-diag(1/as.vector(eta),nrow=nrow(x),ncol=nrow(x))
    z<-y
    b.new<-solve(t(x)%*%w%*%x)%*%t(x)%*%w%*%z
    change <- sqrt(sum((b.new - b.old)^2))
    b.old <- b.new
  }
  b.new
}

poisson_IRLS(x,y,b.init)

##          [,1]
## [1,] 7.451633
## [2,] 4.935300

m1<-glm(y ~ x1, family=poisson(link="identity"),data=data_4.3)
summary(m1)

## 
## Call:
## glm(formula = y ~ x1, family = poisson(link = "identity"), data = data_4.3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7019  -0.3377  -0.1105   0.2958   0.7184  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.4516     0.8841   8.428  < 2e-16 ***
## x1            4.9353     1.0892   4.531 5.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.4206  on 8  degrees of freedom
## Residual deviance:  1.8947  on 7  degrees of freedom
## AIC: 40.008
## 
## Number of Fisher Scoring iterations: 3

inv_link<-function(x){exp(x)}

poisson_IRLS_log_link<-function(x,y,b.init,inv_link,tol=1e-8){
  change <- Inf
  b.old <- b.init
  while(change > tol){
    eta <- x %*% b.old  # linear predictor
    y.hat<-inv_link(eta)  #y.hat is mu in fact
    w<-diag(as.vector(y.hat),nrow=nrow(x),ncol=nrow(x))
    z<-solve(w)%*%(y-y.hat)
    b.new<-b.old+solve(t(x)%*%w%*%x)%*%t(x)%*%w%*%z
    change <- sqrt(sum((b.new - b.old)^2))
    b.old <- b.new
  }
  b.new
}

poisson_IRLS_log_link(x,y,rep(1,2),inv_link)

##           [,1]
## [1,] 1.8892720
## [2,] 0.6697856

m2<-glm(y ~ x1, family=poisson(link="identity"),data=data_4.3)
summary(m2)

## 
## Call:
## glm(formula = y ~ x1, family = poisson(link = "identity"), data = data_4.3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7019  -0.3377  -0.1105   0.2958   0.7184  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.4516     0.8841   8.428  < 2e-16 ***
## x1            4.9353     1.0892   4.531 5.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.4206  on 8  degrees of freedom
## Residual deviance:  1.8947  on 7  degrees of freedom
## AIC: 40.008
## 
## Number of Fisher Scoring iterations: 3