Log-Likelihood Output in DeepPumas Model

[SaiAhil]

Hello Team,
Today I was working on covariate modelling using Neural networks in the DeepPumas model. So initially while we do the traditional NLME model, we get the marginal likelihood of the fit displayed in the results and it matches with the optimization as well

But when I do the DeepPumas based model, with NN in it, we still go the loglikelihood values during the fit but it resembled the NLME values. But to my understanding since we use regularization in the DeepPumas model, while the likelihood is calculated, it has to be a combination of Marginal and the penalty right. Thats when I saw that at the last optimization process, we can see a different value of Log likelihood and the final printed results gave a different one and resembled the Marginal lieklihood of the NLME approach

So here we can see the last optmization step ends at 4.88e02 but the results we got was 504 which matched the marginal Likelihood results without the prior influence. I just wanted to confirm whether this is the scenario and what is the inference and how both these values are useful?

Also have attached the jl file codes I used for the reference:

using Pumas, DeepPumas
using DataFrames, Random, StableRNGs, Printf

Random.seed!(20250715)

# ── constants  ────────────────────────────────────────
TV_CL, TV_VC = 10.0, 100.0     # CL (L/hr), Vc (L)
OM2          = 0.3^2           # IIV variance on CL & Vc
SIGMA        = 0.1             # residual SD (exponential error)
NSUBJ        = 50
TIMES        = [1.75, 7.0]     # 2 samples
DOSE         = 1.0
BETA_COV     = 0.5             # TRUE effect: +50% CL when c1 == 1

# ── models ───────────────────────────────────────────────────────────────────
# base: covariate-blind
basemodel = @model begin
    @param begin
        tvCL ∈ RealDomain(; lower = 0.0, init = TV_CL)
        tvVc ∈ RealDomain(; lower = 0.0, init = TV_VC)
        Ω    ∈ PDiagDomain(; init = [OM2, OM2])
        σ    ∈ RealDomain(; lower = 0.0, init = SIGMA)
    end
    @random η ~ MvNormal(Ω)
    @covariates c1 c2
    @pre begin
        CL = tvCL * exp(η[1])
        Vc = tvVc * exp(η[2])
    end
    @dynamics Central1
    @derived begin
        cp := @. Central / Vc
        dv ~ @. LogNormal(log(cp), σ)
    end
end

# classical covariate model: CL = tvCL*(1 + θCOV*c1)*exp(η)   (df = 1)
nlme_c1 = @model begin
    @param begin
        tvCL ∈ RealDomain(; lower = 0.0, init = TV_CL)
        θCOV ∈ RealDomain(; init = 0.0)
        tvVc ∈ RealDomain(; lower = 0.0, init = TV_VC)
        Ω    ∈ PDiagDomain(; init = [OM2, OM2])
        σ    ∈ RealDomain(; lower = 0.0, init = SIGMA)
    end
    @random η ~ MvNormal(Ω)
    @covariates c1 c2
    @pre begin
        CL = tvCL * (1 + θCOV * c1) * exp(η[1])
        Vc = tvVc * exp(η[2])
    end
    @dynamics Central1
    @derived begin
        cp := @. Central / Vc
        dv ~ @. LogNormal(log(cp), σ)
    end
end

# DeepNLME: NN(c1,c2) -> CL shift, weights co-estimated. L2(1.0) prior on the NN.
deepnlme = @model begin
    @param begin
        tvCL ∈ RealDomain(; lower = 0.0, init = TV_CL)
        tvVc ∈ RealDomain(; lower = 0.0, init = TV_VC)
        NN   ∈ MLPDomain(2, 3, 3, (1, identity); reg = L2(1.0))
        Ω    ∈ PDiagDomain(; init = [OM2, OM2])
        σ    ∈ RealDomain(; lower = 0.0, init = SIGMA)
    end
    @random η ~ MvNormal(Ω)
    @covariates c1 c2
    @pre begin
        CL = tvCL * exp(NN(c1, c2)[1] + η[1])
        Vc = tvVc * exp(η[2])
    end
    @dynamics Central1
    @derived begin
        cp := @. Central / Vc
        dv ~ @. LogNormal(log(cp), σ)
    end
end

# ── simulate from the TRUE model (nlme_c1 at θCOV = BETA_COV) ──────────────────
rng     = StableRNG(20250715)
c1v     = [rand(rng) < 0.5 ? 1.0 : 0.0 for _ in 1:NSUBJ]   # binary covariate (real)
c2v     = randn(rng, NSUBJ)                                 # continuous covariate (null decoy)
regimen = DosageRegimen(DOSE; cmt = 1, time = 0.0)
skel    = [Subject(; id = i, events = regimen, covariates = (; c1 = c1v[i], c2 = c2v[i]))
           for i in 1:NSUBJ]
truth   = merge(init_params(nlme_c1), (; θCOV = BETA_COV))
pop     = read_pumas(DataFrame(simobs(nlme_c1, skel, truth; obstimes = TIMES, rng = rng));
                     observations = [:dv], covariates = [:c1, :c2])

# ── fits ──────────────────────────────────────────────────────────────────────
fit_base = fit(basemodel, pop, init_params(basemodel), FOCE())
fit_nlme = fit(nlme_c1,   pop, init_params(nlme_c1),   FOCE())

# DeepNLME: random NN init, pin structural tv, MAP(FOCE)
init_deep = merge(sample_params(deepnlme; rng = StableRNG(20250715)),
                  (; tvCL = TV_CL, tvVc = TV_VC))
fit_deep  = fit(deepnlme, pop, init_deep, MAP(FOCE()); ignore_numerical_error = true)

# ── log-likelihoods ───────────────────────────────────────────────────────────
println("\n", "═"^60)
@printf("base      loglikelihood(fpm) = %.4f\n", loglikelihood(fit_base))
@printf("nlme_c1   loglikelihood(fpm) = %.4f   (θCOV = %.3f)\n",
        loglikelihood(fit_nlme), coef(fit_nlme).θCOV)
println("─"^60)
# The question: for the MAP fit, is this the marginal or the penalised value?
@printf("deepnlme  loglikelihood(fpm)                         = %.4f\n",
        loglikelihood(fit_deep))
@printf("deepnlme  loglikelihood(model, pop, coef, FOCE())    = %.4f  <- plain-FOCE marginal\n",
        loglikelihood(deepnlme, pop, coef(fit_deep), FOCE()))
println("═"^60)
```{julia}

Thank you

[korsbo]

Hi there,

That’s the expected behaviour and you’re all good. The parameters of the neural network are just fixed effects of the NLME model and the regularization is a prior over those fixed effects. The marginal likelihood does not factor such fixed-effect priors in - it only includes the prior over the things we marginalize out (the random effects).

For the priors over the fixed effects to have any effect, we need to use the maximum aposteriori (MAP) as a target for optimization and that’s what you follow in the optimization trace. So, in your deepnlme model, you fitted the MAP (not just the loglikelihood). But for model evaluation and comparison, the loglikelihood (on withheld data) is the way to go.

I hope this helps

[SaiAhil]

Hi Korsbo,
Thank you so much for the clarification!