2
$\begingroup$

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:

# INPUT
initial_state_values <- c(S = 762,  
                          I = 1,       
                          R = 0)

times <- seq(from = 0, to = 14, by = 0.1)

# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {  
  
  with(as.list(c(state, parameters)), {
    
    N <- S+I+R
    
    lambda <- beta * I/N
    
    # The differential equations
    dS <- -lambda * S               
    dI <- lambda * S - gamma * I
    dR <- gamma * I             
    
    # Output
    return(list(c(dS, dI, dR))) 
  })
}

The translated Mathematica code

data = {{10, 19}, {17, 4}, {24, 5}, {32, 47}, {36, 38}, {38, 54}, {41,
     84}, {42, 55}, {47, 116}, {50, 102}, {52, 126}, {58, 153}, {59, 
    164}, {62, 172}, {63, 177}, {67, 142}, {68, 158}, {69, 151}, {75, 
    155}, {80, 125}, {81, 116}, {84, 101}, {87, 106}, {88, 116}, {94, 
    85}, {98, 59}, {102, 67}, {105, 57}, {114, 59}, {117, 24}, {118, 
    39}, {120, 31}, {128, 36}, {134, 8}, {140, 7}, {142, 3}, {147, 
    20}, {150, 7}, {155, 2}, {159, 25}, {163, 2}, {186, 10}, {191, 
    17}, {200, 7}};

dataTable = 
 TableForm[data, 
  TableHeadings -> {None, {"Time", "Number Infected"}}]

(*total size of population*)
 n = 763;

times = data[[All, 1]];

(*Flu dataset*)
numberinfected = data[[All, 2]]

ClearAll[\[Beta], \[Gamma]]
sireqns = {s'[t] == -\[Beta]*s[t]*i[t]/n, 
   i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t], 
   r'[t] == \[Gamma]*i[t]};

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

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

Now I am trying to translate the following code from R to Mathematica, so I could calculate the parameters. But I face difficulties how to do it.

My main difficulty how can I create this code that calculates the Poisson log-likelihood

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

So then I will try to find through this method the parameters with greater accuracy than with the (sum-of-squares) method.

# 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
    
   # 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)))   # the "parameters" input argument of the loglik_function()
    
  
   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)

Could anyone help me, please?

$\endgroup$
2
  • 2
    $\begingroup$ Since not all of use know R, we should post the ideas and the process of R instead of only post the code of R. $\endgroup$
    – cvgmt
    Sep 1, 2023 at 8:47
  • 1
    $\begingroup$ "But I face difficulties..." Specifying the difficulties would likely get you more and specific help. $\endgroup$
    – JimB
    Sep 1, 2023 at 14:21

1 Answer 1

3
$\begingroup$

It looks like we need to compute $\beta, \gamma$ using optimization with NMinimize in the form

data = {{10, 19}, {17, 4}, {24, 5}, {32, 47}, {36, 38}, {38, 54}, {41,
     84}, {42, 55}, {47, 116}, {50, 102}, {52, 126}, {58, 153}, {59, 
    164}, {62, 172}, {63, 177}, {67, 142}, {68, 158}, {69, 151}, {75, 
    155}, {80, 125}, {81, 116}, {84, 101}, {87, 106}, {88, 116}, {94, 
    85}, {98, 59}, {102, 67}, {105, 57}, {114, 59}, {117, 24}, {118, 
    39}, {120, 31}, {128, 36}, {134, 8}, {140, 7}, {142, 3}, {147, 
    20}, {150, 7}, {155, 2}, {159, 25}, {163, 2}, {186, 10}, {191, 
    17}, {200, 7}};

dataTable = 
  TableForm[data, 
   TableHeadings -> {None, {"Time", "Number Infected"}}];

(*total size of population*)
n = 763;

times = data[[All, 1]];

(*Flu dataset*)
numberinfected = data[[All, 2]];

ClearAll[\[Beta], \[Gamma]]
sireqns = {s'[t] == -\[Beta]*s[t]*i[t]/n, 
   i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t], 
   r'[t] == \[Gamma]*i[t]};

initialConditions = {s[0] == n, i[0] == 1, r[0] == 0};
infVal = 
  ParametricNDSolveValue[{sireqns, initialConditions}, 
   i[#] & /@ times, {t, 0, Last[times]}, {\[Beta], \[Gamma]}];



func[b_?NumericQ, g_?NumericQ] := Norm[infVal[b, g] - numberinfected];


sol1 = NMinimize[{func[b, g], {b > 0, g > 0}}, {b, g}]

(*Out[]= {84.3981, {b -> 0.180905, g -> 0.0748557}}*)

Visualization

Show[ListLinePlot[
  Transpose[{times, Evaluate[infVal[b, g] /. sol1[[2]]]}]], 
 ListPlot[data, PlotStyle -> Red]]

Figure 1

$\endgroup$
2
  • $\begingroup$ First of all, thank you for your answer. Is this way like the Poisson log-likelihood? $\endgroup$ Sep 1, 2023 at 20:10
  • 1
    $\begingroup$ @AthanasiosParaskevopoulos If you need exactly the Poisson log-likelihood, then we have function LogLikelihood[] and PoissonDistribution[\[Mu]] to optimize parameters. $\endgroup$ Sep 2, 2023 at 6:02

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.