I am trying to extract the loglikelihood function that maps the parameter vector to each observation. The goal is to use the likelihood function to generate the score and hessian for each observation to use for inference.

The ideal would be that there is a method of `GLM.loglik_obs`

spit out after the estimation that maps the model parameter vector to a vector equal to the sample size.

Is there a way to do that?

The alternative I am using is to code the loglikelihood myself. See below a working example for a logit MLE.

Thanks!

```
#Julia Version 1.9.4
#Commit 8e5136fa297 (2023-11-14 08:46 UTC)
#Build Info:
# Official https://julialang.org/ release
# [38e38edf] GLM v1.8.3
import GLM
import StatsFuns as StatF
import DataFrames as DF
import StatsBase
import CovarianceMatrices as COV
function main()
data = DF.DataFrame(Y=[1,1,1,1,0,0,0,0], X=[1,2,4,2,11,7,9,12])
estimates = GLM.glm(GLM.@formula(Y ~ 0+ X), data, GLM.Bernoulli(), GLM.LogitLink())
θhat , ll0 = GLM.coef(estimates) , GLM.loglikelihood(estimates)
function loglike_obs(θ)
let Y=GLM.response(estimates), XX=estimates.mm.m
Z = XX * θ
ll = map( (y,z) -> y == 1 ? -StatF.log1pexp(-z) : -StatF.log1pexp(-z) - z, Y,Z )
return ll
end
end
loglike(θ)= sum(loglike_obs(θ))
ϵ = 1e-7
fscore(θ) = ( loglike_obs(θ .+ ϵ) .- loglike_obs(θ) ) ./ ϵ
score = fscore(θhat)
hess = ( fscore(θhat .+ ϵ) .- fscore(θhat) ) ./ ϵ
SSp = StatsBase.mean( x->x^2, score)
Vh = inv( - StatsBase.mean(hess) )
Vs = Vh*SSp*Vh
println(" GLM Log likelihood: ", ll0, " , manual ll: ", loglike(θhat), " , θ: ", θhat )
println(" Manual SE from hessian: ", sqrt( Vh / GLM.nobs(estimates) ) , ", SE GLM: ", GLM.stderror(estimates) |> only )
println(" Manual SE from sandwich: ", sqrt( Vs / GLM.nobs(estimates) ) , ", SE CovarianceMatrices HC0: ", COV.stderror(COV.HC0(), estimates) |> only )
return estimates
end #main
main()
```