# Iteratively Reweighted Least Squares Poisson Regression with Non-Canonical Link, eventually in R

I am attempting to implement iteratively reweighted least squares (IRLS) for a poisson regression with a non-canonical link function. My understanding of IRLS is only basic at this point, so please bear with me. For the typical poisson IRLS model, I would implement the following:

$$\beta_{old} = 1$$

For i in maxit:


$$Z = X\beta_{old} + \dfrac{y - \exp(X\beta_{old})}{\exp(X\beta_{old})} \implies X\beta_{old} + \dfrac{y}{\exp(X\beta_{old})}-1$$

$$W = \text{diag} [\exp(X\beta_{old})]$$

$$\beta_{new} = (X'WX)^{-1}X'WZ$$

STOP IF $$\sqrt{(\beta_{old}-\beta_{new})'(\beta_{old}-\beta_{new})} <$$ tolerance

with the following R-code implementing this procedure:

pois_irls = function(y, X, bstarts=NULL, maxit=1000, tol=1e-12){

if(is.null(bstarts)) b_old = rep(1, ncol(X)) else b_old=bstarts

cat("Starting iterations at:", Sys.time(), "\n")

for(i in seq_len(maxit)){

z = as.vector(X%*%b_old + (y - exp(X%*%b_old))/exp(X%*%b_old))

W = diag(as.vector(exp(X%*%b_old)))

b_new = as.vector(solve(crossprod(X, W%*%X), crossprod(X, W%*%z)))

if(sqrt(crossprod(b_new - b_old)) < tol) break

b_old = as.vector(b_new)
}

cat("Iterations completed at:", Sys.time(), "\n")
cat("Total Iterations:", i, "\n")

return(list(params = b_new, infomat = crossprod(X, W%*%X)))
}

# RUN MODEL
N=100

X = cbind(int=rep(1, N), x1 = rnorm(N), x2 = runif(N, -1,1))
B = c(1.2,2,-3.375)
y = rpois(N, exp(X%*%B))

myirls = pois_irls(y, X)
myirls


Okay, so, this works as expected. However, I am wondering how I can implement the same algorithm but with the following link function: (is it even possible?)

$$\lambda = \exp(X\beta) + \rho y_{t-1}$$

where $$y_{t-1}$$ is the first lag of the outcome. Thus, my understanding is that the link is a mix between identity and log. Is this possible to do, both mathematically speaking, and in code? I have looked into the math of IRLS, but it seems that the derivatives in the Taylor Expansion are always written with respect to a single parameter, so I haven't been able to figure this out independently. In any case, below is some R-code that I've used in an attempt to do so:

linfeed_sim=function(n,X,B,rho){ #function to simulate data from DGP

lambda=matrix(NA, nrow=n)
y=matrix(NA, nrow=n)
p=length(rho)
eXB=exp(X%*%B)

for(i in 1:n){

if(i<(p+1)) {
lambda[i]=eXB[i]
y[i]=rpois(1,lambda[i])
}
else {
lambda[i]=y[(i-p):(i-1)]%*%rev(rho)+eXB[i]
y[i]=rpois(1, lambda[i])
}

}
return(y)
}

N=100

X = cbind(int=rep(1, N), x1 = rnorm(N), x2 = runif(N, -1,1))
B = c(1.2,2,-3.375)
rho = .35

y = linfeed_sim(N, X, B, rho)
lagy = dplyr::lag(y, 1)

df = cbind.data.frame(X, lagy, y)
df = na.omit(df)

X = as.matrix(df[,1:ncol(X)])
lagy = as.vector(df[,(ncol(df)-1)])
y = as.vector(df[,ncol(df)])

#ATTEMPTED FUNCTION
pois_irls_lf = function(y, X, lagy, bstarts=NULL, maxit=1000, tol=1e-8){

if(is.null(bstarts)) b_old = rep(1, ncol(X))
rho_old = 1

cat("Starting iterations at:", Sys.time(), "\n")
for(i in seq_len(maxit)){
cat(i, "\n")
z = as.vector(X%*%b_old + rho_old*lagy + (y - (exp(X%*%b_old)+rho_old*lagy))/(exp(X%*%b_old)+rho_old*lagy))

W = diag(as.vector(exp(X%*%b_old)+rho_old*lagy))

Xlagy = cbind(X, lagy)

b_new = as.vector(solve(crossprod(Xlagy, W%*%Xlagy), crossprod(Xlagy, W%*%z))) # W is diagonal so multiplication in this way is okay

if(sqrt(crossprod(b_new - c(b_old,rho_old))) < tol) break

b_old = as.vector(b_new[-length(b_old)])
rho_old = as.vector(b_new[length(b_old)])

}
cat("Iterations completed at:", Sys.time(), "\n")
#cat("Total Iterations:", i, "\n")
return(list(params = b_new, infomat = crossprod(Xlagy, W%*%Xlagy)) )
}

pois_irls_lf(y, X, lagy) # ends up singular before itmax


but this runs into singularity issues after about 600 iterations. Any help is appreciated. I wasn't sure if here or Overflow was better for this post, as I am wondering about both the math and the code; let me know if this should be migrated or if you need more info from me. Thanks in advance.

• As for on-topicness, if it requires statistical expertise to answer it should be on-topic; so at least those aspects of the question that require such expertise should be fine, but please focus the question on what you're trying to do rather than on the specifics of your implementation (e.g. if you make a coding error in implementing what you're trying to do that's not a statistical issue and fixing that would be off topic) – Glen_b Mar 25 at 4:42