Estimación de máxima verosimilitud de los parámetros del modelo de infección

Question

Estoy usando el modelo estándar de infección en algunos datos con los que estoy trabajando.

$ dS = -\beta SI $

$ dI = \beta SI - \gamma I $

$ dR = \gamma I $

Donde $S$ es el número de sujetos susceptibles, $I$ es el infectado y $R$ es el recuperado. Estoy probando varios métodos para estimar los parámetros $\beta$ y $\gamma$.

---

Para cualquier período de tiempo discreto y de ancho fijo, conozco el número de infectados y la población total, la cual es fija. Uno de los métodos que he utilizado para estimar los parámetros es alimentar el estado inicial en un solucionador de ecuaciones diferenciales en R y recorrer varios valores para $\beta$ y $\gamma$ hasta que minimicen el Error Cuadrático Medio.

Me han dicho que es posible hacer una estimación de máxima verosimilitud de los parámetros, pero no tengo idea de cómo empezar a hacer esto.

---

Una idea que me presentaron implicaba usar una curva normal y estimar los parámetros de la distribución utilizando las conocidas estimaciones de máxima verosimilitud para los parámetros de las distribuciones normales. Mi problema con esto es que estoy tratando con el número de infectados (o incluso la proporción de infectados) en una población, no con algo que siga las suposiciones necesarias de una distribución de probabilidad.

Si quisiera hacer esto, tendría que introducir otro parámetro para desplazar la curva normal hacia arriba hacia mis datos. Es decir, si $f(t;\mu,\sigma)$ es la distribución normal, necesitaría otro parámetro $k>1$ tal que

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

donde $I_t$ es el número de infectados (o la proporción, no importa) en el tiempo $t$.

Si usara este método, creo que estimaría los parámetros $\beta$ y $\gamma$ de la siguiente manera:

1. Usar la estimación de máxima verosimilitud de la curva normal usando el método anterior.
2. Usar la estimación de mínimos cuadrados como en esta pregunta para ajustar el modelo de infección a la curva normal en lugar de a los datos reales.

---

Realmente no estoy seguro de qué más hacer, así que agradezco mucho cualquier información que puedas proporcionar.

Answer 1

Aquí hay una posibilidad en la que el modelo se modifica (1) para ser explícitamente probabilístico y (2) para tener lugar en tiempo discreto.

El código a continuación explica el modelo modificado, lo simula y luego utiliza MLE para recuperar los parámetros (cuyo valor real se conoce en este ejemplo ilustrativo, ya que simulamos los datos). Cuidado: mi beta no será exactamente equivalente a tu beta, consulta la "historia" en los comentarios a continuación.

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
## Historia: cada infectado tiene probabilidad gamma de recuperarse durante el período; antes de que las recuperaciones se realicen, cada susceptible interactúa con cada infectado; cada interacción conduce a una infección con probabilidad beta; susceptible se contagia si >= 1 de sus interacciones conduce a una infección

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)  # Comprobación de cordura
    }
    return(df)
}

inverse_logit <- function(x) {
    p <- exp(x) / (1 + exp(x))  # Mapea R a [0, 1]
    return(p)
}
curve(inverse_logit, -10, 10)  # Comprobación de cordura

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)]  # Cuidado, problemático si es 1 o 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

## Distribución de muestreo de beta_gamma_hat
estimates <- replicate(100, get_estimates())
df_estimates <- as.data.frame(t(estimates))
summary(df_estimates)  # Parece razonable dado los valores reales de (0.005, 0.10)

Avísame si algo no es autoexplicativo.

Aviso legal: No he estudiado mucho el modelo SIR excepto una vez muy brevemente en una clase universitaria, hace varios años. El modelo que simulo y estimo arriba no es exactamente el clásico modelo de ecuaciones diferenciales SIR que mencionaste en tu pregunta. Además, me siento un poco febril hoy, ¡así que revisa el código por si hay errores!

Answer 2

Estas diapositivas proporcionan un buen comienzo inicial al hacer ajustes iterativos de mínimos cuadrados, sin embargo, esto es solo para el modelo SI.

Esta entrada de blog da a este problema un tratamiento más extenso. Esto parece ser similar a tu enfoque. Ajustar un núcleo gaussiano a través de mínimos cuadrados no lineales podría ser una buena forma de obtener estimaciones iniciales de parámetros. Luego debes identificar algún tipo de relación entre los parámetros del núcleo gaussiano y los parámetros en tu modelo.