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

Question

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)[1]
p = dim(X)[2]
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))
  }

data = load(url("https://www.math.ntnu.no/emner/TMA4315/2020h/hoge-veluwe.Rdata"))
myglm(y~t, data)

Answer 1

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.