# \[ANN\] ProfileLikelihood.jl: Computing profile likelihoods

**URL:** <https://discourse.julialang.org/t/ann-profilelikelihood-jl-computing-profile-likelihoods/91823>\
**Category:** Package Announcements\
**Tags:** package, statistics, inference\
**Created:** [December 19, 2022, 12:00am UTC](https://discourse.julialang.org/t/ann-profilelikelihood-jl-computing-profile-likelihoods/91823 "2022-12-19T00:00:34Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [December 19, 2022, 12:00am UTC](https://discourse.julialang.org/t/ann-profilelikelihood-jl-computing-profile-likelihoods/91823/1 "2022-12-19T00:00:34Z")

</div>

Hi all,

I have just released my package [ProfileLikelihood.jl](https://github.com/DanielVandH/ProfileLikelihood.jl) that will be fully registered in a [few days](https://github.com/JuliaRegistries/General/pull/74323). The aim of this package is to provide an easy interface for computing maximum likelihood estimates (MLEs), profile likelihoods, and confidence and prediction intervals. By building on top of the great interface from [Optimization.jl](https://github.com/SciML/Optimization.jl), it is very easy to switch between different optimisers, use different types of automatic differentiation, to implement constraints.

I define what profile likelihoods, etc. are in the README of the package, but to quickly summarise: Suppose we have a likelihood function \mathcal L(\boldsymbol \theta) (basically a joint probability density function for some sample), and assume that we have some interest parameter \psi so that \boldsymbol \theta is decomposed as \boldsymbol \theta = (\psi, \boldsymbol \omega), where \boldsymbol \omega are the remaining ‘nuisance’ parameters. Currently I only support scalar \psi (edit: **as of v0.2, bivariate profiles are now also supported** ). The _profile likelihood_ of \psi is defined by \mathcal L\_p(\psi\_0) = \sup\_{\boldsymbol \omega \in \Omega \mid \psi = \psi\_0} \mathcal L(\psi, \boldsymbol \omega), i.e. the maximum of the likelihood when we fix \psi at a given value. This profile likelihood makes it easy to compute confidence intervals and approximate prediction intervals. Again, more detail is in the README, and a nice very recent resource is [Simpson and Maclaren (2022)](https://doi.org/10.1101/2022.12.14.520367).

To illustrate the package, I give one of the five examples from the docs below, where I estimate parameters for a logistic ODE, reproducing the first case study of [Simpson and Maclaren (2022)](https://doi.org/10.1101/2022.12.14.520367). The main idea is to first define the likelihood function so that a `LikelihoodProblem` can then be built. Computing the results is as easy as then using `mle`, `profile`, and then `compute_prediction_intervals` (if you have a prediction function you are interested in). Other examples for multiple linear regression, a simpler ODE \mathrm du/\mathrm dt = \lambda u, and estimating parameters for a PDE (the heat equation) are given in the README.

Suppose we have \mathrm du/\mathrm dt = \lambda u(1-u/K), u(0)=u\_0, and our interest is in estimating (\lambda, K, u\_0), fixing the standard deviation of the noise, \sigma, at \sigma=10 for this example. First, we need some data for this synthetic study:

```julia
using OrdinaryDiffEq, Random
λ = 0.01
K = 100.0
u₀ = 10.0
t = 0:100:1000
σ = 10.0
@inline function ode_fnc(u, p, t)
    λ, K = p
    du = λ * u * (1 - u / K)
    return du
end
# Initial data is obtained by solving the ODE 
tspan = extrema(t)
p = (λ, K)
prob = ODEProblem(ode_fnc, u₀, tspan, p)
sol = solve(prob, Rosenbrock23(), saveat=t)
Random.seed!(2828)
uᵒ = sol.u + σ * randn(length(t))

```

Next, a likelihood function is defined. I define the log-likelihood function, and use a normal likelihood since I added Gaussian noise to the data:

```julia
using ProfileLikelihood
function loglik_fnc2(θ, data, integrator)
    λ, K, u₀ = θ
    uᵒ, σ = data
    integrator.p[1] = λ
    integrator.p[2] = K
    reinit!(integrator, u₀)
    solve!(integrator)
    return gaussian_loglikelihood(uᵒ, integrator.sol.u, σ, length(uᵒ))
end

```

When solving ODEs, it is assumed that the third argument is an `integrator`, coming from the integrator interface from DifferentialEquations.jl. In other cases, there are only the first two arguments (see the regression example in the README). You could re-define an `ODEProblem` many times, but this is slow. To now build the `LikelihoodProblem` from all this, and constraining 0 \leq \lambda \leq 0.05, 50 \leq K \leq 150, and 0 \leq u\_0 \leq 50:

```julia
using Optimization
lb = [0.0, 50.0, 0.0] # λ, K, u₀
ub = [0.05, 150.0, 50.0]
θ₀ = [λ, K, u₀]
syms = [:λ, :K, :u₀]
prob = LikelihoodProblem(
    loglik_fnc2, θ₀, ode_fnc, u₀, maximum(t); 
    syms=syms,
    data=(uᵒ, σ),
    ode_parameters=[1.0, 1.0], # temp values for [λ, K]
    ode_kwargs=(verbose=false, saveat=t),
    f_kwargs=(adtype=Optimization.AutoFiniteDiff(),),
    prob_kwargs=(lb=lb, ub=ub),
    ode_alg=Rosenbrock23()
)

```

This function will construct the `LikelihoodProblem`, a function of two arguments that in this case is a closure that contains the integrator for the problem. From this point, it is very easy to find the MLEs and the profile likelihoods:

```julia
using OptimizationNLopt
sol = mle(prob, NLopt.LD_LBFGS)

```

```julia
LikelihoodSolution. retcode: Failure
Maximum likelihood: -38.99053694428977
Maximum likelihood estimates: 3-element Vector{Float64}
     λ: 0.010438031266786045
     K: 99.59921873132551
     u₀: 8.098422110755225

```

```julia
prof = profile(prob, sol; alg=NLopt.LN_NELDERMEAD, parallel=false)

```

```julia
ProfileLikelihoodSolution. MLE retcode: Failure
Confidence intervals: 
     95.0% CI for λ: (0.006400992274213644, 0.01786032876226762)
     95.0% CI for K: (90.81154862835605, 109.54214763511888)
     95.0% CI for u₀: (1.5919805025139593, 19.070831536649305)

```

We now have estimates for our parameters, along with quantified uncertainty from the confidence intervals. Using the `plot_profiles` function, you can also visualise the profile likelihoods to assess e.g. smoothness and identifiability of the individual parameters:

```julia
using CairoMakie, LaTeXStrings
fig = plot_profiles(prof; latex_names=[L"\lambda", L"K", L"u_0"],
    show_mles=true, shade_ci=true, nrow=1, ncol=3,
    true_vals=[λ, K, u₀], # optional, gives the black lines below
    fig_kwargs=(fontsize=30, resolution=(2109.644f0, 444.242f0)),
    axis_kwargs=(width=600, height=300))

```

 ![logistic_example](https://global.discourse-cdn.com/julialang/original/3X/2/9/2942e845e9693a7b450d7688a4e70d27aae9a283.png)

A last feature to show, again following [Simpson and Maclaren (2022)](https://doi.org/10.1101/2022.12.14.520367), is the construction of profile-wise profile likelihoods, which give profile likelihoods for some prediction function q(\boldsymbol \theta) conditional on an interest parameter \psi. In this case, suppose we are interested in quantifying the uncertainty in the function u(t) over many times. First, we define a prediction function that will give us our vector of interest:

```julia
function prediction_function(θ, data)
    λ, K, u₀ = θ
    t = data
    prob = ODEProblem(ode_fnc, u₀, extrema(t), (λ, K))
    sol = solve(prob, Rosenbrock23(), saveat=t)
    return sol.u
end

```

Computing prediction intervals from this is simple using `get_prediction_intervals`:

```julia
t_many_pts = LinRange(extrema(t)..., 1000)
parameter_wise, union_intervals, all_curves, param_range =
    get_prediction_intervals(prediction_function, prof,
        t_many_pts; parallel = true)

```

The `parameter_wise` output gives prediction intervals based on varying one parameter at a time, and `union_intervals` gives the union of the prediction intervals at each time; see the docs for more detail. Plotting these prediction intervals together with the true curve in red + the curve from the MLE in blue, we obtain the figure below (code for making this plot is in the docs):

 ![logistic_example_prediction](https://global.discourse-cdn.com/julialang/original/3X/6/d/6d40d3f6076f27bb08a8cc718007768fcce1e915.png)

In (a), (b), and (c), the black curves show 95% profile-wise prediction intervals for the respective parameters. These views allow us to identify how each parameter contributes to the uncertainty in our prediction, e.g. the uncertainty from \lambda is largest for early time, which makes sense since only K matters for large time, and this is what we see in (b) where the interval becomes large only for long time. In (d) I show the union of all these intervals in the black curve, compared to the interval coming from the full likelihood (which is more accurate, but takes a significantly longer time to compute) in magenta. We see that we have fully enclosed the true solution, and this sensitivity analysis of the individual contributions to the total uncertainty is sensible.

All feedback welcome!  
Thanks,  
Daniel.

---

<div class="post-metadata">

**Author:** ![Julius\_Martensen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/julius_martensen/32/17077_2.png) [@Julius\_Martensen](https://discourse.julialang.org/u/Julius_Martensen)\
**Post date:** [December 19, 2022, 5:48am UTC](https://discourse.julialang.org/t/ann-profilelikelihood-jl-computing-profile-likelihoods/91823/2 "2022-12-19T05:48:39Z")

</div>

Thanks! I literally stumbled upon the mentioned paper two days ago and wanted to implement this 🫶 This looks awesome and well done 👍

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [December 19, 2022, 6:31am UTC](https://discourse.julialang.org/t/ann-profilelikelihood-jl-computing-profile-likelihoods/91823/3 "2022-12-19T06:31:49Z")

</div>

The package README.md is very impressive. The documentation is great.

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [December 28, 2022, 3:00am UTC](https://discourse.julialang.org/t/ann-profilelikelihood-jl-computing-profile-likelihoods/91823/4 "2022-12-28T03:00:00Z")

</div>

An update as of v0.2.0: Bivariate profiles are now supported, with a complete example given in [Example V: Lotka-Volterra ODE, GeneralLazyBufferCache, and computing bivarate profile likelihoods · ProfileLikelihood.jl](https://danielvandh.github.io/ProfileLikelihood.jl/dev/lotka/), where I demonstrate how to compute bivariate profile likelihoods for a Lotka-Volterra ODE problem, how to support automatic differentiation, and how to compute prediction intervals from a bivariate profile (this is the second case study of [Simpson and Maclaren (2022)](https://doi.org/10.1101/2022.12.14.520367)).