# Implementing The Fisher Scoring Algorithm in R for a Poisson GLM, Fisher Matrix not Square

I've been trying to implement the Fisher Scoring Algorithm in R for a Poisson GLM with the canonical choice of link function. I want the function to return a list containing:

coef: a matrix with the beta parameter estimates in the first column, and their standard errors in the second.

deviance

vcov: Covariance matrix of the coefficients.

The Fisher matrix comes out as a $$135X1$$ matrix which means I cant invert it later to find the covariance matrix. In summation form I found the score function and Fisher matrix to be: $$s(\beta) = \sum_{i=1}^nx_i(y_i-\lambda_i)$$ and $$F = \sum_{i=1}^n x_i x_i^T \lambda_i$$ And $$\eta = ln(\lambda) = X \beta$$

$$\lambda = exp(X \beta)$$

myglm <-  function(formula,data,start = 0) {

X = model.matrix(formula,data)
Y = data[,1]
n = dim(X)
p = dim(X)
beta_0 = rep(1,p)
beta = solve(t(X) %*%X) %*% t(X) %*% Y #Least Squares Estimate
epsilon = 0.01

#Run Fisher Iterations
while (norm(beta-beta_0,type = "2")/norm(beta_0, type = "2") > epsilon) {

beta_0 = beta
eta = X %*% beta
lambda = exp(eta)
F = X %*% t(X) %*% lambda #Fisher information matrix
s = X %*% (Y - lambda) #Score function

beta = beta + solve(F) %*% s

}

vcov = solve(F)
coef[,1] = beta
coef[,2] = sqrt(diag(vcov)) #Standard errors of the coefficients

#Calculate Deviance
mod_sat = glm(formula, family = poisson(link = "log"))
log_likelihood = Y %*% eta - exp(eta)
deviance = 2*(LogLik(mod_sat) - log_likelihood)

return(list(coef,deviance,vcov))
}

myglm(y~t, data)



It seems your while statement has the wrong inequality: the rhs should be larger than epsilon, not smaller. That is,

while (norm(beta-beta_0,type = "2")/norm(beta_0, type = "2") > epsilon)


is probably what you want.

With the wrong inequality, it is highly likely that your program will finish without even starting the Fisher iterations. This leads to the problem that the object F is never assigned and will evaluate as its default value, namely FALSE.

• Thank you for the tip. Having the inequality to be $> \epsilon$ was definitely what I intended, so this helped me move in the right direction. Now I am rid of the $F = False$ problem, however $F$ now comes out as a $135X1$ matrix(not square) which means I can't invert it later to find vcov. I don't see any obvious errors or changes I should do so that it will turn out square.
– Pame
Sep 28 '20 at 19:52
• I think you should check this line: F = X %*% t(X) %*% lambda. Try to calculate the components in small steps and see where it goes wrong. I would guess that X %*% t(X) seems to give the wrong dimensions. Sep 28 '20 at 20:02