Extracting the loglikelihood function from GLM?

Sorry, missed that you want to pass own parameters.
In that case, I guess that GLM does not support it and your manual approach is fine. I would probably just go ahead and code the full GLM manually (evtl. even using Flux or something for parameter handling):

import Distributions as PD
import GLM
import StatsModels

function myglmloglik(fm, data, params)
    y, X = StatsModels.modelcols(StatsModels.apply_schema(fm, StatsModels.schema(fm, data)), data)
    p = GLM.linkinv.(GLM.LogitLink(), X * params)
    PD.logpdf.(PD.Bernoulli.(p), y)
end

loss(θ) = sum(myglmloglik(fm, data, θ))

At least, this would work with several AD backends and the model could easily be generalized by passing a Flux layer instead of params.

2 Likes