Incorporating time lag in observation models

Hello. I am estimating parameters of a simple SIR epidemic model (shown below) by fitting the model to a known time series data:

Deterministic model:

function SIR_model(du, u, p, t)
    S, I, R = u
    β, γ = p

    du[1] = -β * S * I
    du[2] = β * S * I - γ * I
    du[3] = γ * I
end

Observation model:

@model function SIR_fitting(data)
    σ ~ Uniform(0, 50)
    β ~ truncated(Normal(5e-05, 1e-03), 0, 1)
    γ ~ truncated(Normal(0.6, 0.3), 0, 1)

    p = [β,γ]
    prob = ODEProblem(SIR_model, u0, (0, 50), p)
    predicted = solve(prob, Tsit5(), saveat=1.0)

    for i = 1:length(predicted)
        data[i] ~ Normal(predicted[i][2], σ)
    end
end

As you can see, I am sampling the observations of I(t) from a normal distribution and fitting them to data. Everything works perfectly with this setup.

Now, I assume there is a time lag τ between data and sampling, representing the delayed reporting of infections in reality. But modifying the last line in the observation model to data[i] ~ Normal(predicted[i-τ][2], σ) throws out Invalid indexing of solution error. I can understand this is due to non-integer indexing. Can someone please please help me and suggest how to avoid it?

Note: I have used the prior distribution of τ ~ Gamma(1,5)

Thanks in advance!

indexing the solution object accesses it at discrete timesteps.
I think what you want here is to interpolate your solution at arbitrary time steps, which can be done by calling (parantheses instead of square brackets) the solution object:

predicted(i - τ)

see here how to handle the solution object in various ways:
https://diffeq.sciml.ai/latest/tutorials/ode_example/#Handling-the-Solution-Type

however you’ll probably still have to make sure that i - τ is within bounds.

1 Like

I am no expert, but I think there are also other ways to model time lags.
This repository has a great overview of such models: https://github.com/epirecipes/sir-julia

Thank you @trahflow. This worked beautifully! :slight_smile:

1 Like