5.3 Example: Poisson Regression

Here is an example of the IRWLS algorithm applied to a Poisson-family GLM using the aids data set. Recall that this data set contains the number of deaths from AIDS in Australia for thee-month periods from 1983 to 1896, demonstrating an increase in the number of cases over time. The underlying GLM for this data set is:

$${\eta }_i={\beta }_1+{\beta }_2{x}_i\mathit{\hspace{1em}}\mathrm{and}\mathit{\hspace{1em}}{\eta }_i=\log ({\mu }_i).$$

The derivative of the link function is:

$$\frac{d{\eta }_i}{d{\mu }_i}=\frac{1}{{\mu }_i}$$

With this information, the following R function performs IRWLS for a Poisson regression model.

irwls <- function(param0, y, X, I){
  #param0 -- Initial parameter vector
  #y -- Realisations
  #X -- Design Matrix
  #I -- Number of iterations to perform
  beta <- param0
  for(i in 1:I){
    eta <- as.vector(X %*% beta)    #Linear Predictor
    mu <- exp(eta)                  #Mean
    var <- diag(mu)                 #Variance
    linkderiv <- 1 / muΨΨΨΨΨΨΨ#Derivative of Link Function
    z <- as.matrix(eta + (y - mu) * linkderiv)  #Adjusted Responses
    W <- solve(var * linkderiv^2)   #Weights Matrix
    betanew <- solve( t(X) %*% W %*% X) %*% t(X) %*% W %*% z
    cat("Params(", i, "):", betanew, "\n")
    beta <- betanew
  }
}

Using starting values ${\beta }_1={\beta }_2=1$, we use the irwls() function to find the MLEs for the AIDS data set:

aids <- read.table("aids.dat")
Responses <- aids$number
DesignMatrix <- cbind(1, aids$time)
irwls(param0 = c(1, 1), y = Responses, X = DesignMatrix , I = 15)

Params( 1 ): 0.001549063 0.9998877
Params( 2 ): -0.9942415 0.9995826
Params( 3 ): -1.982809 0.9987539
Params( 4 ): -2.951808 0.9965066
Params( 5 ): -3.868091 0.9904381
Params( 6 ): -4.644543 0.9742338
Params( 7 ): -5.06551 0.9322646
Params( 8 ): -4.683145 0.8321057
Params( 9 ): -3.013692 0.6391485
Params( 10 ): -0.8658216 0.4160803
Params( 11 ): 0.05815249 0.3004687
Params( 12 ): 0.3091476 0.2617296
Params( 13 ): 0.3391299 0.256613
Params( 14 ): 0.3396338 0.2565236
Params( 15 ): 0.3396339 0.2565236

Verify convergence by using different starting values to obtain the same result. The estimates using the glm() are:

coef( glm(number ~ 1 + time, family=poisson, data=aids) )

(Intercept)         time
     0.3396       0.2565

Summary

• •

  The IRWLS algorithm is Fisher’s scoring method and the resulting solution correspond to maximum likelihood estimations.

• •

  The algorithm consists of a number of simple steps using weighted least squares which are iterated to produce a sequence of parameter values.

• •

  The iterations terminate when this sequence converges to a specified accuracy.