2
$\begingroup$
  • All infected individuals go through an incubation period, which lasts on average for 4 days. During this time, individuals are not infectious, nor do they have any excess mortality risk.

  • All infected individuals eventually develop symptoms, and the mean duration of symptoms before either recovery or death is 5 days. Symptomatic individuals are infectious, as well as having a 3% case fatality rate. Those who survive the infection are thought to have long-term immunity.

  • In the source country, the peak prevalence (i.e. maximum number of symptomatic people during the epidemic) was observed to be 8% of the population.

The total population of this country is N=11.000.000

\begin{align*} \text{Case Fatality Rate (CFR)} &= 3\% = 0.03 \\ \text{Mean duration of symptoms (}D\text{)} &= 5 \text{ days} \end{align*}

The rate of recovery ($\gamma$) can be calculated using the formula:

$$ \gamma = \frac{1 - \text{CFR}}{D} $$

$$\gamma = \frac{0.97}{5} = 0.194 \text{ per day}$$

Then ($\mu$) is $$CFR = \frac{\mu}{\mu + \gamma}$$ So $$\mu=0.006$$ The estimated value of $(I_{peak})$ can be calculated as: $$I_{\text{peak}} = \text{Peak Prevalence} \times N = 0.08 \times 11,000,000 = 880,000$$

The differential equations for this model is \begin{align*} \frac{dS}{dt}&= -\beta \cdot S(t)\cdot \frac{I(t)}{N} - \mu \cdot S(t)\\ \frac{dI}{dt} &= \beta \cdot S(t) \cdot \frac{I(t)}{N} - \gamma \cdot I(t) - \mu \cdot I(t)\\ \frac{dR}{dt} &= \gamma \cdot I(t) - \mu \cdot R(t) \end{align*} What estimate do get for $\beta$ which corresponds to a peak prevalence (of symptomatic infection) of 0.08(i.e. 8%)

(*total size of population*)n = 11000000;
times = 100;
reportedData = n*0.08
sireqns = {s'[t] == -\[Beta]*s[t]*i[t]/n - \[Mu]*s[t], 
   i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t] - \[Mu]*i[t], 
   r'[t] == \[Gamma]*i[t] - \[Mu]*r[t]};

initialConditions = {s[0] == n, i[0] == 1, r[0] == 0};

sol = ParametricNDSolve[{sireqns, initialConditions}, {s, i, r}, {t, 
   0, times}, {\[Beta], \[Gamma], \[Mu]}]



(*Extract the infected values from the solution*)
infectedValues = i /. sol;

I have tried the following code to Define a function to calculate the Poisson log-likelihood but it does not work.

loglik[\[Beta]_, \[Gamma]_, \[Mu]_, data_] := 
 Module[{output, modelI, \[Lambda], loglik}, 
  output = {i} /. sol[\[Beta], \[Gamma], \[Mu]];
  modelI = 
   Interpolation[Transpose[{Range[0, times], output[[1]]}], 
    InterpolationOrder -> 1];
  \[Lambda] = 0.6*modelI[data[[All, 1]]];
  loglik = 
   Total[Log[PoissonDistribution[\[Lambda]] /@ data[[All, 2]]]];
  loglik]

result = 
 NMaximize[{loglik[\[Beta], \[Gamma], \[Mu], reportedData], 
   0 <= \[Beta] <= 1, 0 <= \[Gamma] <= 1, 
   0 <= \[Mu] <= 1}, {\[Beta], \[Gamma], \[Mu]}]

estimatedParameters = {\[Beta], \[Gamma], \[Mu]} /. result[[2]]

enter image description here I want to define a function that simulates the model for a given combination of parameters and calculates the Poisson log-likelihood for the epidemic curve of reported cases as R code shows us:

# DISTANCE FUNCTION

loglik_function <- function(parameters, dat) {   # takes as inputs the parameter values and dataset

  beta <- parameters[1]    # extract and save the first value in the "parameters" input argument as beta
  gamma <- parameters[2]   # extract and save the second value in the "parameters" input argument as gamma
    mu <- parameters[2]   # extract and save the second value in the "parameters" input argument as mu
 
  # Simulate the model with initial conditions and timesteps defined above, and parameter values from function call
  output <- as.data.frame(ode(y = initial_state_values, 
                              times = times, 
                              func = sir_model,
                              parms = c(beta = beta,       # ode() takes the values for beta and gamma extracted from
                                        gamma = gamma,
mu=mu)))   # the "parameters" input argument of the loglik_function()

So far I have figured out how to code it. My question is how do we generate this code in Mathematica, do I need Poisson Distribution? Could some one explain to me how could i code that?

  LL <- sum(dpois(x = dat$number_reported, lambda = 0.6 * output$I[output$time %in% dat$time], log = TRUE))
   
  return(LL) 
}


# OPTIMISATION

optim(par = c(1.7, 0.1),           # starting values for beta and gamma - you should get the same result no matter 
                                  # which values you choose here
     fn = loglik_function,        # the distance function to optimise
     dat = reported_data,         # the dataset to fit to ("dat" argument is passed to the function specified in fn)
     control = list(fnscale=-1))  # tells optim() to look for the maximum number instead of the minimum (the default)

Also, I know that the result of $\beta=0.38$

If someone needs it this the whole code I used in R

# Load required packages
library(deSolve)

# Define initial conditions
initial_state_values <- c(S = 11000000,  
                        I = 1,       
                        R = 0)

# Define time points
times <- seq(from = 0, to = 100, by = 0.1)

# Define SIR model function
sir_model <- function(time, state, parameters) {  
  with(as.list(c(state, parameters)), {
      N <- S + I + R
      lambda <- beta * I / N
      dS <- -lambda * S - mu * S              
      dI <- lambda * S - gamma * I - mu * I
      dR <- gamma * I - mu * R             
      return(list(c(dS, dI, dR))) 
  })
}

# Known values
gamma <- 0.194
mu <- 0.006
Ipeak <- 880000

# Define a range of possible beta values
possible_beta_values <- seq(0, 1, by = 0.01)

# Find beta that results in Ipeak
for (beta in possible_beta_values) {
  parameters <- c(beta = beta, gamma = gamma, mu = mu)
  output <- as.data.frame(ode(y = initial_state_values, 
                               times = times, 
                               func = sir_model,
                               parms = parameters))
  
  peak_time <- times[which.max(output$I)]
  
  if (output$I[peak_time / 0.1 + 1] >= Ipeak) {
      cat("Found beta:", beta, "\n")
      break
  }
}
$\endgroup$
4
  • 2
    $\begingroup$ "... but it does not work" - how does it fail exactly? Where is the failure? $\endgroup$
    – MarcoB
    Sep 2, 2023 at 1:54
  • $\begingroup$ Parameter data is not defined. Is it the same as in your previous question? $\endgroup$ Sep 2, 2023 at 11:44
  • $\begingroup$ You need something like Total[Map[Log[PDF[PoissonDistribution[\Lambda], #]&, data[[All,2]] ].]. This is equivalent to using the density dpois. The PDF is missing in what you did. $\endgroup$
    – Asim
    Sep 2, 2023 at 13:10
  • $\begingroup$ @MarcoB I get the following message when I run the code, you will see it in my updated post $\endgroup$ Sep 2, 2023 at 14:14

1 Answer 1

2
$\begingroup$

It looks like we try to estimate parameter $\beta$ with known $\gamma, \mu$ using one point only reportedData = n*0.08. In this case we can compute $\beta$ as follows

n = 11000000;
times = 100;
reportedData = n*0.08; 
model[\[Beta]_?NumberQ] := (model[\[Beta]] = 
   First[i /. 
     NDSolve[{s'[t] == -\[Beta]*s[t]*i[t]/n - \[Mu]*s[t], 
        i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t] - \[Mu]*i[t], 
        r'[t] == \[Gamma]*i[t] - \[Mu]*r[t], s[0] == n, i[0] == 1, 
        r[0] == 0} /. {\[Gamma] -> 0.194, \[Mu] -> 0.006}, 
      i, {t, times}]]); sol = 
 NMinimize[{(Max[model[b][#] & /@ Range[0, times]] - 
      reportedData)^2, {b > 0}}, {b}, MaxIterations -> 1000]

(*Out[]= {1.95698*10^-17, {b -> 0.510274}}*)

Visualization

Plot[model[b][t] /. sol[[2]], {t, 0, 100}]

Figure 1

We also can try to estimate $\beta, \gamma, \mu$ using one point only

model1[\[Beta]_?NumberQ, \[Gamma]_?NumberQ, \[Mu]_?
   NumberQ] := (model1[\[Beta], \[Gamma], \[Mu]] = 
   First[i /. 
     NDSolve[{s'[t] == -\[Beta]*s[t]*i[t]/n - \[Mu]*s[t], 
       i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t] - \[Mu]*i[t], 
       r'[t] == \[Gamma]*i[t] - \[Mu]*r[t], s[0] == n, i[0] == 1, 
       r[0] == 0}, i, {t, times}]]); sol1 = 
 NMinimize[{(Max[
       model1[\[Beta], \[Gamma], \[Mu]][#] & /@ Range[0, times]] - 
      reportedData)^2, {0 <= \[Beta] <= 1, 0 <= \[Gamma] <= 1, 
    0 <= \[Mu] <= 1}}, {\[Beta], \[Gamma], \[Mu]}, 
  MaxIterations -> 1000]

(*Out[]= {4.87078*10^-11, {\[Beta] -> 0.88748, \[Gamma] -> 
   0.035164, \[Mu] -> 0.0405108}}*)

Visualization

Plot[
 model1[\[Beta], \[Gamma], \[Mu]][t] /. sol1[[2]], {t, 0, 100}]

Figure 2

What is the advance to use LogLikelihood[] or loglik?

$\endgroup$
3
  • $\begingroup$ The professor had said that the poisson log-likelihood gives better accuracy $\endgroup$ Sep 2, 2023 at 13:20
  • $\begingroup$ Because the estimate value of $\beta$=0.38 $\endgroup$ Sep 2, 2023 at 14:29
  • 1
    $\begingroup$ Try model with $\beta=0.38$ at $\gamma= 0.194, \mu -> 0.006$, You will get Max[model[0.38][#] & /@ Range[0, times]] of about 5178.92, or 0.047%. No, I don't believe your professor :) $\endgroup$ Sep 2, 2023 at 14:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.