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.

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[1] ~ dpois(lambda = beta[1] * S[1] * I[1] / population)
deltaR[1] ~ dpois(lambda = gamma[1] * I[1])
beta[1] ~ dnorm(mean = 0, sd = sqrt(5 * sigma2beta))
gamma[1] ~ 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 `SIRInits`

contains 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)
```

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[1] ~ dpois(lambda = beta[1] * S[1] * I[1] / population)
deltaR[1] ~ dpois(lambda = gamma[1] * I[1])
beta[1] ~ dnorm(mean = 0, sd = sqrt(5 * sigma2beta))
gamma[1] ~ 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)
}
})
```

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[1] ~ dpois(lambda = beta[1] * S[1] * I[1] / population)
deltaR[1] ~ dpois(lambda = gamma[1] * I[1])
beta[1] ~ dnorm(mean = 0, sd = sqrt(5 * sigma2beta))
gamma[1] ~ 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)
}
})
```

Gilks, W. R., A. Thomas, and D. J. Spiegelhalter. 1994. “A Language and Program for Complex Bayesian Modelling.” *Journal of the Royal Statistical Society. Series D (the Statistician)* 43 (1): 169–77.

Valpine, Perry de, Daniel Turek, Christopher J Paciorek, Clifford Anderson-Bergman, Duncan Temple Lang, and Rastislav Bodik. 2017. “Programming with Models: Writing Statistical Algorithms for General Model Structures with NIMBLE.” *Journal of Computational and Graphical Statistics* 26 (2): 403–13.