0

I have the following data set:

y <- c(5,8,6,2,3,1,2,4,5)
x <- c(-1,-1,-1,0,0,0,1,1,1)
d1 <- as.data.frame(cbind(y=y,x=x))

I when I fit a model to this data set with glm(), using a Poisson distribution with a log link:

model <- glm(y~x, data=d1, family = poisson(link="log"))
summary(model)

I get the following output:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.3948     0.1671   8.345   <2e-16 ***
x            -0.3038     0.2250  -1.350    0.177

I want to write a function for the iterative re-weighted least squares regression that will obtain the same estimates. So far I have been able to do this using an identity link, but not a log link, as I do in the glm.

X <- cbind(1,x)

#write an interatively reweighted least squares function with log link
glmfunc.log <- function(d,betas,iterations=1)
{
X <- cbind(1,d[,"x"])
z <- as.matrix(betas[1]+betas[2]*d[,"x"]+((d[,"y"]-exp(betas[1]+betas[2]*d[,"x"]))/exp(betas[1]+betas[2]*d[,"x"])))


for(i in 1:iterations) {
W <- diag(exp(betas[1]+betas[2]*d[,"x"]))
betas <- solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%z
}
return(list(betas=betas,Information=t(X)%*%W%*%X))
}

#run the function
model <- glmfunc.log(d=d1,betas=c(1,0),iterations=1000)

Which gives output:

#print betas
model$betas
           [,1]
[1,]  1.5042000
[2,] -0.6851218

Does anyone know where I am going wrong in writing the custom function and how I would correct this to replicate the output from the glm() function

Cettt
  • 11,460
  • 7
  • 35
  • 58
solventarmadillo
  • 3
  • 2
  • there's an example of a simple IRLS (but using `lm.wfit` rather than crude linear algebra) at https://ms.mcmaster.ca/~bolker/classes/s4c03/notes/week4A.pdf – Ben Bolker Sep 20 '17 at 00:23

1 Answers1

1

It appears your 'z' needs to be inside of your loop, as your 'betas' get updated each iteration, thus so should your 'z' as it is based on those values.

The implementation looks right otherwise.

EVB
  • 26
  • 2