This webpage is created as an online supplementary material for the manuscript Time Fused Coefficient SIR Model with Application to COVID-19 Epidemic in the United States. MCMC sampling in the manuscript is implemented in nimble, an integrated system that enables one to program for general model structures in R. It facilitates flexible model specification, enjoys nice high-level programmibility, and generates c++ code for fast computation. The syntax is an extension of the BUGS language (Gilks, Thomas, and Spiegelhalter 1994). For interested readers, details about the nimble package can be found in Valpine et al. (2017). We present the set up of the three models in nimble here.

## 1. The $$t$$-Shrinkage Prior

We start with the $$t$$-shrinkage prior. The model is wrapped in the nimbleCode function. We start the model by declaring the conditional distributions and the priors for first difference at $$t=1$$ (lines 2–7). Note that the code is written based on simulated data with time ranging from 0 to 80, so that all successive differences are computed for $$t=1,\ldots,80$$. The other three variables are the number of observations within each of the S, I, and R compartments. Parameter deltaM denotes the quantity $$\Delta M(t)$$, and deltaR denotes the quantity $$\Delta R(t)$$ in Equation (4) in the manuscript. Both deltaM and deltaR are vectors of length 80, corresponding to $$\Delta M(t)$$ and $$\Delta R(t)$$ for $$t=1,\ldots,80$$. Variables beta and gamma are also vectors of length 80, which represent $$\beta(t)$$ and $$\gamma(t)$$ for $$t=1,\ldots,80$$. Variables sigma2beta and sigma2gamma correspond to components in the normal distributions for $$\beta(t) - \beta(t-1)$$ and $$\gamma(t) - \gamma(t-1)$$, and they are both assigned non-informative priors of $$\mbox{IG}(0.5,0.5)$$.

Next we iterate through the time after $$t=1$$ (lines 9–17). Variables deltaM and deltaR still follow the specified Poisson distribution, and the incremental in $$\beta(t)$$, $$\Delta \beta(t)$$, denoted in the code as theta_beta[t-1] for convenience of indexing, follows $$\mbox{N}(0,\lambda_{t-1}\sigma_\beta^2)$$ with the prior for $$\lambda_{t-1}$$ specified to be $$\mbox{IG}(2, 0.0001)$$. A similar setting is applied on incremental in $$\gamma(t)$$ as well, denoted by theta_gamma.

SIRCode <- nimbleCode({
deltaM ~ dpois(lambda = beta * S * I / population)
deltaR ~ dpois(lambda = gamma * I)
beta ~ dnorm(mean = 0, sd = sqrt(5 * sigma2beta))
gamma ~ dnorm(mean = 0, sd = sqrt(5 * sigma2gamma))
sigma2beta ~ dinvgamma(0.5, 0.5)
sigma2gamma ~ dinvgamma(0.5, 0.5)

for (t in 2:T) {
deltaM[t] ~ dpois(lambda = beta[t] * S[t] * I[t] / population)
deltaR[t] ~ dpois(lambda = gamma[t] * I[t])
beta[t] <- beta[t - 1] + theta_beta[t - 1]
theta_beta[t - 1] ~ dnorm(0, sd = sqrt(sigma2beta * lambda[t - 1]))
lambda[t - 1] ~ dinvgamma(2, 0.0001)
gamma[t] <- gamma[t - 1] + theta_gamma[t - 1]
theta_gamma[t - 1] ~ dnorm(0, sd = sqrt(sigma2gamma * eta[t - 1]))
eta[t - 1] ~ dinvgamma(2, 0.0001)
}
})

After the model is defined, model constants and initial values are supplied. The variable population is the total number of observations, which is set to 1,000,000, and remains constant. The list SIRInitscontains initial values used for simulation, and finally the MCMC is compiled and invoked by the function nimbleMCMC(). With niter = 50, thin = 10, nchains = 1, the final output, mcmc.out, contains 5,000 rows of results, which after the burn-in component is discarded, can be used for Bayesian estimation and inference.

T <- 80
population <- 10 ^ 6

SIRdata <- list(population = population, deltaN = deltaN,
deltaR = deltaR, S = S, I = I)
SIRConsts <- list(T = T)

SIRInits <- list(beta = rep(0.3, T), gamma = rep(0.01, T),
eta = rep(1, T - 1), lambda = rep(5, T - 1),
theta_beta = rep(0, T - 1), theta_gamma = rep(0, T - 1),
sigma2beta = 1,sigma2gamma = 1)

mcmc.out <- nimbleMCMC(code = SIRCode, data = SIRdata,
constants = SIRConsts,
inits = SIRInits,
monitors = c("beta", "gamma", "lambda", "eta"),
niter = 50000,
thin = 10,
nchains = 1)

## 2. The Horseshoe Prior

For the horseshoe prior, the initial states are set similar to above. Therefore we place the emphasis on the iteration part to update the day-to-day increments. Here the variable theta_beta[t-1] follows $$N(0,c^2 \lambda_t^2)$$, where the half Cauchy-distributed $$\lambda_t$$ is constructed as the absolute value of the ratio between two standard normally distributed variables a_lam[t-1] and b_lam[t-1] as shown in line 14 in the code chunk below. Similarly, $$\eta_t$$ is constructed using a_eta[t-1] and b_eta[t-1] for theta_gamma[t-1].

SIRCode <- nimbleCode({
deltaN ~ dpois(lambda = beta * S * I / population)
deltaR ~ dpois(lambda = gamma * I)
beta ~ dnorm(mean = 0, sd = sqrt(5 * sigma2beta))
gamma ~ dnorm(mean = 0, sd = sqrt(5 * sigma2gamma))
sigma2beta ~ dinvgamma(0.5, 0.5)
sigma2gamma ~ dinvgamma(0.5, 0.5)

for (t in 2:T) {
deltaN[t] ~ dpois(lambda = beta[t] * S[t] * I[t] / population)
deltaR[t] ~ dpois(lambda = gamma[t] * I[t])
beta[t] <- beta[t - 1] + theta_beta[t - 1]
theta_beta[t - 1] ~ dnorm(0, sd = sqrt(sigma2beta * lambda[t - 1] ^ 2))
lambda[t - 1] <- abs(a_lam[t - 1] / b_lam[t - 1])
a_lam[t - 1] ~ dnorm(0, 1)
b_lam[t - 1] ~ dnorm(0, 1)
gamma[t] <- gamma[t - 1] + theta_gamma[t - 1]
theta_gamma[t - 1] ~ dnorm(0, sd = sqrt(sigma2gamma * eta[t - 1] ^ 2))
eta[t - 1] <- abs(a_eta[t - 1] / b_eta[t - 1])
a_eta[t - 1] ~ dnorm(0, 1)
b_eta[t - 1] ~ dnorm(0, 1)
}
})

## 3. The Spike-and-slab Prior

Similar to Section 2, we focus on the iteration part. The spike and slab for $$\Delta \beta(t)$$ are, respectively, denoted by theta_beta2[t-1] and theta_beta1[t-1]. The variable lambda[t-1] represents the parameter $$\lambda_{t-1}$$ in Equation (6), and is assumed to follow Bernoulli(0.5). Similar constructions are applied to the spike and slab for variables related to $$\Delta \gamma(t)$$, including theta_gamma2[t-1], theta_gamma1[t-1], and eta[t-1].

SIRCode <- nimbleCode({
deltaN ~ dpois(lambda = beta * S * I / population)
deltaR ~ dpois(lambda = gamma * I)
beta ~ dnorm(mean = 0, sd = sqrt(5 * sigma2beta))
gamma ~ dnorm(mean = 0, sd = sqrt(5 * sigma2gamma))
sigma2beta ~ dinvgamma(0.5, 0.5)
sigma2gamma ~ dinvgamma(0.5, 0.5)

for (t in 2:T) {
deltaN[t] ~ dpois(lambda = beta[t] * S[t] * I[t] / population)
deltaR[t] ~ dpois(lambda = gamma[t] * I[t])
beta[t] <- beta[t - 1] + lambda[t - 1] * theta_beta1[t - 1] +
(1 - lambda[t - 1]) * theta_beta2[t - 1]
theta_beta1[t - 1] ~ dnorm(0, sd = sqrt(sigma2beta))
theta_beta2[t - 1] ~ dnorm(0, sd = sqrt(sigma2beta / 100))
lambda[t - 1] ~ dbern(0.5)
gamma[t] <- gamma[t - 1] + eta[t - 1] * theta_gamma1[t - 1] +
(1 - eta[t - 1]) * theta_gamma2[t - 1]
theta_gamma1[t - 1] ~ dnorm(0, sd = sqrt(sigma2gamma))
theta_gamma2[t - 1] ~ dnorm(0, sd = sqrt(sigma2gamma / 100))
eta[t - 1] ~ dbern(0.5)
}
})