Simple SIR epidemiological model using pMCMC

[Bharadwaj Vemparala]

Hello,

I aim to use the pMCMC scheme for a complex epidemiology problem, but since I am a beginner with NIMBLE, I generated test data using the standard SIR (Susceptible-Infected-Recovered) model with known parameters and tried to estimate back those parameters. The R code and the model are given below.

Model:

dS/dt = -beta*S*I

dI/dt = beta*S*I - gamma*I

dR/dt = gamma*I

Code:

# STEP 1: PRELIMINARIES

library(nimble)
library(nimbleSMC)
library(ggplot2)
library(coda)
set.seed(1)

# STEP 2: DEFINING THE MODEL

SIRCode <- nimbleCode({
  # Known parameters
  N         <- 10000
  I0        <- 10
 
  # Initial conditions
  x[1,1] <- N-I0
  x[2,1] <- I0
  x[3,1] <- 0
  y[1]   <- N-I0
 
  # Discretized SIR model
  for (t in 2:totalTime){
    x[1,t] <- x[1,t-1] + deltaT*(-beta*x[1,t-1]*x[2,t-1])
    x[2,t] <- x[2,t-1] + deltaT*(beta*x[1,t-1]*x[2,t-1]
                                 - gamma*x[2,t-1])
    x[3,t] <- x[3,t-1] + deltaT*(gamma*x[2,t-1])
    y[t]   <- x[1,t]
  }
 
  # Prior distributions of unknown parameters
  # beta  ~ dunif(min=0, max=1e-03) # Actual value 1e-04
  # gamma ~ dunif(min=0.1, max=1)     # Actual value 0.75
  beta  ~ dbeta(shape1=2,shape2=2)
  gamma ~ dbeta(shape1=2,shape2=2)
})

# STEP 3: EXTRACTING THE REQUIRED DATA
# First column is time (days), second is the S population, and
# third column is the second column with added white noise

SIR_data <- read.table("new_SIR_data.txt",header=FALSE,sep=",")
y_data <- SIR_data[,2]

# STEP 4: CREATE AND COMPILE THE MODEL

SIRModel <- nimbleModel(code=SIRCode,
                        data=list(y=y_data),
                        constants=list(totalTime=50,deltaT=1),
                        inits=list(beta=1e-05,gamma=0.5))
compileNimble(SIRModel)

# STEP 5: CONFIGURE THE MCMC ALGORITHM

SIRModelMCMCSpec <- configureMCMC(SIRModel,nodes=NULL,
                                  monitors=c('beta','gamma',"x","y"))
SIRModelMCMCSpec$addSampler(
  target=c('beta','gamma'),
  type='RW_PF_block',
  control=list(
    adaptive=TRUE,
    latents="x",
    resample=TRUE,
    pfType="bootstrap",
    pfOptimizeNparticles=TRUE
  )
)

# STEP 6: BUILD, COMPILE, AND RUN

SIRModelMCMC <- buildMCMC(SIRModelMCMCSpec)
cMCMC <- compileNimble(SIRModelMCMC,
                       project=SIRModel,
                       resetFunctions=TRUE)
resultsMCMC <- runMCMC(cMCMC,niter=10000,
                       nchains=1,nburnin=200,
                       summary=TRUE)

Unfortunately, the whole program runs but the parameters are never sampled; traceplots show a straight line for both beta and gamma at the initial condition. I do not really understand where am I wrong! I tried with different prior distributions for the parameters, different data sets, and other conditions but it is the same situation with parameters just stuck at the initial guess.

Can someone please help me out as I am a novice in this field?

A ton of thanks in advance,

Bharadwaj

[Perry de Valpine]

Dear Bharadwaj,

Thanks for the message and reproducible example.

I think you need both state dynamics and observations to be stochastic.  That is, you need x[1, t], x[2, t] and x[3, t] each to follow a distribution.   Similarly, you need the observations, y[t], to follow a distribution.  (If one or the other is assumed deterministic, you wouldn't really need PMCMC.)

Let me know if that helps.

-Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/0258e914-3b7c-4120-a9f8-d834de750e25n%40googlegroups.com.