Model Inference
Constructing formula and running the model.
Overview
Having completed data preprocessing and data organization, a formula is constructed and the model run.
Note: This script requires that both the data Preprocessing and the Data Organization have already been completed with all resulting objects available in the working environment.
Model
\begin{align} \textit{y}_{1\textit{st}} &= \beta_0 + \beta_1{Pop} + \varphi_{\textit{t}} + \textit{W}_1(\textit{st}) \nonumber \\ \textit{y}_{2\textit{st}}|\textit{y}_{1\textit{st}} &= \beta_0 + \beta_2{Pop} + \gamma_{seas} + \varphi_{\textit{t}} +\delta_{\textit{st}} + \textit{W}_2(\textit{st}) + \alpha_1 \cdot \textit{W}_1(\textit{st}) \nonumber \\ \textit{y}_{3\textit{st}}|\textit{y}_{1\textit{st}}, \textit{y}_{2\textit{st}} &= \beta_0 + \beta_3{Pop} + \varphi_{\textit{t}} + \alpha_2 \cdot \gamma_{seas} + \alpha_3 \cdot \textit{W}_1(\textit{st}) + \alpha_4 \cdot \delta_{\textit{st}} + \alpha_5 \cdot \textit{W}_2(\textit{st}) + \textit{W}_3(\textit{st}) \nonumber \\ \nonumber \end{align}
Where,
- \textit{s} (\textit{s} = 1, 2, 3, \ldots,56) are U.S. States and Territories
- \textit{t} are time steps (yearly for spatial effects (\textit{W}_(\textit{st}))
- \textit{y}_{1\textit{st}} is the percent infected from the ILI data set (Level 1)
- \textit{y}_{2\textit{st}} is the percent confirmed positive Influenza A from NREVSS (Level 2)
- \textit{y}_{3\textit{st}} is hospital incidence from HHS (Level 3)
- \beta_0’s are level-specific intercepts
- \beta_2{Pop} are level-specific population estimates for each state
- \textit{W}_(\textit{st}) are BYM space-time effects specific to each level (1,2,3)
- \gamma_{seas} is a week-based seasonal effect estimated from NREVSS
- \delta_{\textit{st}} is a space-time interaction (STI) term estimated from NREVSS
- \varphi_{\textit{t}} is a second-order random walk estimated from FluSurv data
Prior Specifications
Hide code
inla.setOption(inla.mode= "experimental")
#Besag-Mollie-York
bym_hyper <- list(
prec = list(
initial = 2,
prior = "pc.prec",
param = c(0.6, 0.01)),
phi = list(
initial = 2,
prior = "pc",
param = c(0.5, 0.5)))
#iid prior
pc_prec_iid <- list(theta = list(prior="pc.prec",
param=c(1, 0.01)))
#ar1 prior
pc_cor_ar1 <- list(theta = list(prior = 'pccor1',
param = c(0.7, 0.5)))
#rw prior
pcprior1 = list(prec = list(prior="pc.prec",
param = c(1, 0.01)))Joint Model Formula
Hide code
Frm.1 = Y ~ -1 + intercept1 + #intercept for ILI submodel
intercept2 + #intercept for NREVSS submodel
intercept3 + #intercept for HHS submodel
l.pop.ili + #log-population: could have been combined, but
l.pop.nv + #wanted to see if effect differd
l.pop.hhs + #by submodel
f(Region.1.ili, # location index specific to ILI submodel
model="bym2", #Besag-Mollie-York model (spatial + region IID)
graph=J, #neighborhood adjacency matrix
constr=TRUE, #enforce a zero-mean latent field
group = year.ili, #Annual time steps,year-specific spatial fields
control.group=list(model="ar1"), #order-1 autoregressive between annual fields
hyper=bym_hyper) + #prior
f(Region.1.nv, # As above but specific to NREVSS submodel
model="bym2",
graph=J,
constr=TRUE,
group = year.nv.1,
control.group=list(model="ar1"),
hyper=bym_hyper) +
f(Region.2.nv, # Create a copy of ILI BYM model, and share it with NREVSS submodel
copy="Region.1.ili", #what to copy
group = year.nv.2,
fixed=FALSE) +
f(Region.1.hhs, # Besag-Mollie-York model specific to HHS submodel
model="bym2", #as above
graph=J,
constr=TRUE,
group = year.hhs.1,
control.group=list(model="ar1"),
hyper=bym_hyper) +
f(Region.2.hhs, # Create a copy of ILI BYM model, and share it with HHS submodel
copy="Region.1.ili",
group = year.hhs.2,
fixed=FALSE) +
f(Region.3.hhs, # Create a copy of NREVSS BYM model, and share it with HHS submodel
copy="Region.1.nv",
group = year.hhs.3,
fixed=FALSE) +
f(ts_week.1.nv, # Use model etimates from NREVSS submodel to identify seasonality
constr=TRUE,
scale.model = TRUE,
model="seasonal",
season.length = 52, #guess at how many weeks between seasonal cycles
hyper=pcprior1) +
f(ts_week.2.nv, # IID term to account for variation outside of
constr=TRUE, # seasonal effect above
model="iid",
hyper=pc_prec_iid) +
f(ts_week.c, # copy seasonal effect to HHS submodel
copy = "ts_week.1.nv",
fixed=FALSE) +
f(ts_week.iid.c, # copy IID extra-seasonal effect to HHS submodel
copy = "ts_week.2.nv",
fixed=FALSE) +
f(ts_week.fs, en_est.s, # perform order-2 random walk across all levels based
constr=TRUE, # on the temporal pattern exhibited by the FluSurv data
scale.model = TRUE,
model="rw2",
hyper=pcprior1) +
f(st_int.nv, #NREVSS space-time interaction (location-week)
constr=TRUE,
model="iid",
hyper=pc_prec_iid) +
f(st_int.c, #copy NREVSS space-time interaction to HHS level
copy = "st_int.nv",
fixed=FALSE)Hyperparameter Estimates
Mean estimates for all 22 model hyperparameters. These were estimated in advance by running a smaller, 5yr model.
Hide code
theta8 = c(0.681697728, -0.492510866, -2.058167639, 0.006865191, -1.002600890, 2.926536468, -0.321834109, 0.153080825, 1.959042282,
-1.832427192, -1.118539999, 1.942998892, 3.403507238, -1.016877442, -3.634452225, -0.407023837, -0.023686173,
0.529309680, 0.332056118, 0.999033065, 0.857719181, 0.528751714)Run Model
Hide code
Joint.mod = inla(Frm.1, #formula from above
num.threads = 8, #cores to use
data = inla.stack.data(joint_stack), # combined DataStack with substacks for all submodels
family = c("gaussian", "gaussian","gaussian"), #responses are normal
verbose = TRUE, # print run process to screen
control.fixed = list(prec = 1, # proper intercept
prec.intercept = 1),
control.predictor = list(
A = inla.stack.A(joint_stack),
compute = TRUE, # include fitted values
link = 1), #default link function settings
control.mode = list(restart = TRUE, theta = theta8), #initial values above
control.inla = list(strategy="gaussian", #strategy to speed up
int.strategy = "eb"), #empirical Bayes (based on modes)
control.compute=list(dic = F, cpo = F, waic = F)) # these just slow the run down
#save resulting model and inputs
#save(list=c("Joint.mod", "J", "joint_stack"), file="./data/runs/model-2023-06-09.RData", version = 2)