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)
Best Answer
It seems your
while
statement has the wrong inequality: the rhs should be larger thanepsilon
, not smaller. That is,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, namelyFALSE
.