Solved – Maximum Likelihood Estimate of Infection Model Parameters

compartmental-modelsdifferential equationsepidemiologymaximum likelihood

I'm using the standard infection model on some data I am working with.

$ dS = -\beta SI $

$ dI = \beta SI – \gamma I $

$ dR = \gamma I $

Where $S$ is the number of susceptible subjects, $I$ is the infected, and $R$ is the recovered. I'm trying various methods for estimating the parameters $\beta$ and $\gamma$.

For any given discrete, fixed-width time period, I know the number of infected and the total population, which is fixed. One of the methods I have used to estimate the parameters is to feed the initial state into a differential equation solver in R and loop through several values for $\beta$ and $\gamma$ until they minimized the Mean Square Error.

I've been told it's possible to do a maximum likelihood estimate of the parameters, but am at a total loss as to how to begin doing this.

One idea that was presented to me involved using a normal curve and estimating the parameters of the distribution using the well-known maximum likelihood estimates for the parameters of normal distributions. My problem with this is that I'm dealing with the number of infected (or even proportion of infected) in a population, not with anything that follows the necessary assumptions of a probability distribution.

If I were to do this, I would need to introduce another parameter to shift the normal curve upward towards my data. By that I mean, if $f(t;\mu,\sigma)$ is the normal distribution, I would need another parameter $k>1$ such that

$$ I_t\approx kf(t;\mu,\sigma)$$

where $I_t$ is the number of infected (or the proportion, doesn't matter) at time $t$.

If I were to do this method, I believe I would estimate the parameters $\beta$ and $\gamma$ by:

  1. Use the maximum likelihood estimate of the normal(ish) curve using the above method.
  2. Use the least-squares estimate as in this question to fit the infection model against the normal curve instead of the actual data.

I'm not really sure what else to do, so I greatly appreciate any insight you can provide.

Best Answer

Here's a possibility in which the model is modified (1) to be explicitly probabilistic, and (2) to take place in discrete time.

The code below explains the modified model, simulates it, and then uses MLE to recover the parameters (whose true value is known in this toy example, since we simulated the data). Careful: my beta will not be exactly equivalent to your beta -- see "story" in the comments below.

library(ggplot2)
library(reshape2)

## S(t) susceptible, I(t) infected, R(t) recovered at time t
## Probabilistic model in discrete time:
## S(t+1) = S(t) - DeltaS(t)
## I(t+1) = I(t) + DeltaS(t) - DeltaR(t)
## R(t+1) = R(t) + DeltaR(t)
## DeltaR(t) ~ Binomial(I(t), gamma) >= 0
## DeltaS(t) ~ Binomial(S(t), 1 - (1 - beta)^I(t)) >= 0
## Story: each infected has probability gamma of recovering during the period;
## before recoveries are realized, each susceptible interacts with each infected;
## each interaction leads to infection with probability beta;
## susceptible becomes infected if >= 1 of her interactions leads to infection

simulate <- function(T=100, S1=100, I1=10, R1=0, beta=0.005, gamma=0.10) {
    stopifnot(T > 0)
    stopifnot(beta >= 0 && beta <= 1)
    stopifnot(gamma >= 0 && gamma <= 1)
    total_pop <- S1 + I1 + R1
    df <- data.frame(t=seq_len(T))
    df[, c("S", "I", "R")] <- NA
    for(t in seq_len(T)) {
        if(t == 1) {
            df$S[t] <- S1
                df$I[t] <- I1
            df$R[t] <- R1
                next
            }
            DeltaS <- rbinom(n=1, size=df$S[t-1], prob=1 - (1-beta)^df$I[t-1])
            DeltaR <- rbinom(n=1, size=df$I[t-1], prob=gamma)
        df$S[t] <- df$S[t-1] - DeltaS
        df$I[t] <- df$I[t-1] + DeltaS - DeltaR
        df$R[t] <- df$R[t-1] + DeltaR
        stopifnot(df$S[t] + df$I[t] + df$R[t] == total_pop)  # Sanity check
    }
    return(df)
}

inverse_logit <- function(x) {
    p <- exp(x) / (1 + exp(x))  # Maps R to [0, 1]
    return(p)
}
curve(inverse_logit, -10, 10)  # Sanity check

loglik <- function(logit_beta_gamma, df) {
    stopifnot(length(logit_beta_gamma) == 2)
    beta <- inverse_logit(logit_beta_gamma[1])
    gamma <- inverse_logit(logit_beta_gamma[2])
    dS <- -diff(df$S)
        dR <- diff(df$R)
    n <- nrow(df)
    pr_dS <- 1 - (1-beta)^df$I[seq_len(n-1)]  # Careful, problematic if 1 or 0
        return(sum(dbinom(dS, size=df$S[seq_len(n-1)], prob=pr_dS, log=TRUE) +
               dbinom(dR, size=df$I[seq_len(n-1)], prob=gamma, log=TRUE)))
}

get_estimates <- function() {
    df <- simulate()
    mle <- optim(par=c(-4, 0), fn=loglik, control=list(fnscale=-1), df=df)
    beta_gamma_hat <- inverse_logit(mle$par)
    names(beta_gamma_hat) <- c("beta", "gamma")
    return(beta_gamma_hat)
}

set.seed(54321999)

df <- simulate()
df_melted <- melt(df, id.vars="t")
p <- (ggplot(df_melted, aes(x=t, y=value, color=variable)) +
      geom_line(size=1.1) + theme_bw() +
      xlab("time") +
      theme(legend.key=element_blank()) +
      theme(panel.border=element_blank()))
p

## Sampling distribution of beta_gamma_hat
estimates <- replicate(100, get_estimates())
df_estimates <- as.data.frame(t(estimates))
summary(df_estimates)  # Looks reasonable given true values of (0.005, 0.10)

Let me know if anything is not self-explanatory.

Disclaimer: I haven't studied the SIR model except once very briefly in a college class, several years ago. The model I simulate and estimate above is not exactly the classic differential equation SIR model you stated in your question. Also I'm feeling a bit feverish today, so check the code for mistakes!

Related Question