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.


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)


data = load(url("https://www.math.ntnu.no/emner/TMA4315/2020h/hoge-veluwe.Rdata"))
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.

  • $\begingroup$ 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. $\endgroup$
    – Pame
    Sep 28 '20 at 19:52
  • 1
    $\begingroup$ 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. $\endgroup$
    – baruuum
    Sep 28 '20 at 20:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.