CODEWITHSUNDEEP
rregressionglmlog-likelihood
solventarmadillo
solventarmadillo

Reputation: 3

Obtain the MLE of betas through iterative re-weighted least squares regression

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

Upvotes: 0

Views: 1244

Answers (1)

EVB
EVB

Reputation: 26

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.

Upvotes: 1

Related Questions