overview.qmd

---
title: " "
bibliography: references.bib
format:
  html: 
    df-print: kable
    code-fold: show
    code-summary: "Hide code"
    code-overflow: wrap
    toc-title: Page Contents
    toc: true
    toc-depth: 2
    toc-location: right
    number-sections: false
    html-math-method: katex
    smooth-scroll: true
editor: source
editor_options: 
  chunk_output_type: console
---
```{=html}
<style type="text/css">

body, td {
   font-size: 13pt;
}
code.r{
  font-size: 9pt;
}
pre {
  font-size: 11pt
}
</style>
```

::: {.callout-note icon=false}
#### Supporting Information   
This site provides a demonstrative workflow and supporting information for the publication:  
**Interrogating Genomes and Geography to Unravel Multiyear Vesicular Stomatitis Epizootics**  
  
Humphreys, Shults, Velazquez Salinas, Bertram, Pelzel-McCluskey, Peters, and  Rodriguez  
  
**Viruses**. 2024; 16(7):1118. [https://doi.org/10.3390/v16071118](https://doi.org/10.3390/v16071118)  
  
::: 
  
::: {.callout-note icon=false} 
#### External Resources     
All code and data are available for download:    
[GitHub site: https://github.com/geoepi/vs-epizootics](https://github.com/geoepi/vs-epizootics)    

Additional supporting information and documents available on the Open Science Framework:  
[Project OSF site: https://osf.io/ghzfq/](https://osf.io/ghzfq/)    
  
:::

  
### Navigation  
The menu at right can be used to navigate through major sections on a webpage.  Links across top of page will direct to genomic data ([GENOMES](https://geoepi.github.io/vs-epizootics/sequences.html)), the GeoEpi home page ([GEOEPI](https://geoepi.github.io/)), or a listing of other GeoEpi projects ([RESEARCH](https://geoepi.github.io/pindex.html)).
   


## Load Libraries  
Load needed R-packages to execute code.  Note that the r-INLA package (@Rue2009, @Lindgren_2015) is not available on CRAN and must be installed from the INLA team website.
```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(tidyterra)
library(raster)
library(ggtree)
library(terra)
library(raptr)
library(rgeos)
library(pals)
library(here)

select <- dplyr::select

library(INLA)
#install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
```

## Custom Functions  
Custom functions to perform data process and create figure plots.  These are available from the linked repository at the top of the script.  
```{r message=FALSE, warning=FALSE}
set.seed(1976)

source(here("./R/utilities.R"))
source_dir("./R")
```


## Get Boundaries  
The *assets* folder includes raster and other data files to perform the analysis
```{r message=FALSE, warning=FALSE}
domain <- vect(here("assets/domain")) # geographic boundary polygons

obs_stack.r <- rast( # Rasters depicting Spatial Random Fields estimated from true data
  list(
    rast(here("assets/rf_2014.tif")),
    rast(here("assets/rf_2015.tif"))
       )
)
```

## Simulate Outbreak  
Virus and disease detection data presented in the manuscript includes the geographic coordinates recording locations for private residences and businesses.  This data is confidential and cannot be shared, therefore, this demonstration first simulates outbreak count data from an SEIR-vector model and then randomly assigns geographic locations based on the Spatial Random Fields estimated from the true data.  This information is included in the raster files loaded in the previous code chunk.  
   
### SEIR-Vector Dynamics  
Simulate Vesicular stomatitis disease detections for 2014.  
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=6}
#| label: fig-seir14
#| fig-cap: "Simulated VSV outbreak for 2014 with the number of susceptible (Sh), exposed (Eh), infectious (Ih), and removed (Rh) individuals."
#|  
seir_vector_14 <- seirsei_ode(Nh = 4100, # number of livestock hosts
                           vect_mult = 250, # number of vectors per host
                           obs_bias = 0.10, # proportion observed/detected
                           targ_year = "2014" # to approximate dates 
                           )

sum(seir_vector_14$adj_inc)

plot_seirsei(seir_vector_14)
```


Repeat for the 2015 outbreak year.
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=6}
#| label: fig-seir15
#| fig-cap: "Simulated VSV outbreak for 2015 with the number of susceptible (Sh), exposed (Eh), infectious (Ih), and removed (Rh) individuals."
#|
seir_vector_15 <- seirsei_ode(Nh = 5000, # number of livestock hosts
                           vect_mult = 250, # number of vectors per host
                           obs_bias = 0.10, # proportion observed/detected
                           targ_year = "2015" # to approximate dates 
                           )

sum(seir_vector_15$adj_inc)

plot_seirsei(seir_vector_15)
```

### Spatial Randomization   
**NOTE:** The data produced through simulation appears similar to that presented in the publication, but is different in several regards.  Because of these differences, results presented here will not exactly match those presented in the paper.  
```{r message=FALSE, warning=FALSE}
sample_2014 <- spatial_sample(inc_df = seir_vector_14, obs_rast = obs_stack.r[[1]])

sample_2015 <- spatial_sample(inc_df = seir_vector_15, obs_rast = obs_stack.r[[2]])
```

View the simulated VS cases in relation to  
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=7}
#| label: fig-simulation
#| fig-cap: "Simulated disease detections overlying Spatial Random Fields estimated from true data as presented in the publication.  Note that although clustering is evident, it is less pronounced as that in true data."  
#| 
join_sims <- rbind(sample_2014, sample_2015)

plot_rast_panels(obs_stack.r, domain, join_sims, brewer.ylorbr(4))
```

## Spatial Triangulation   
Construct a triangulated mesh to apply Stochastic Partial Differential Equations (see, @SPDEbook).  This mesh is at more coarse resolution than that used for the actual analysis in effort to reduce computational demand. 
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=6}
#| label: fig-mesh
#| fig-cap: "Trianglulated mesh with vvertices aligned to simulated disease locations.  Note that vertex intersections represent nodes."
#|
hull <- st_convex_hull(
  st_buffer(
    st_union(
      st_as_sf(domain)),
    100))

dom_bnds <- inla.sp2segment(as(hull, "Spatial"))

set.seed(1976)
mesh.dom <- inla.mesh.2d(boundary = dom_bnds, 
                        loc = join_sims[,c("x","y")],
                        cutoff = 20, 
                        max.edge = c(60, 300),
                        offset = c(100,250),
                        min.angle = 30) 


plot_mesh(mesh.dom)
```

### Node Coordinates  
Extracting the geographic coordinates for node locations. These are used for model integration and represent background characteristics with respect to environmental conditions.
```{r message=FALSE, warning=FALSE}
dd = as.data.frame(cbind(mesh.dom$loc[,1], 
                         mesh.dom$loc[,2]))

names(dd) = c("x", "y")

dd$set <- "node"
```

```{r}
#| label: tbl-node_df
#| tbl-cap: Quick look at node data attributes

head(dd)
```

### Copy to Year  
Create a separate copy of the nodes for each year of interest, join the nodes with simulated outbreaks, and assign a random calendar date to each.  The nodes are joined to the observation data where a fictitious predictor variable is created (*pred_var*) and a column obs is added to code a 1 for disease observations and a 0 for mesh nodes.  Note that the *pred_var* are simulated using different mean values so that years show a different average.  This process is performed for each year separately.   
```{r message=FALSE, warning=FALSE}
seq_2014_dates <- seq(as.Date("2014-01-01"), as.Date("2014-12-31"), by = 1)

sample_2014 <- rbind(sample_2014, mutate(dd, 
                                         year = 2014,
                                         date = sample(seq_2014_dates, n(), replace = TRUE))) %>%
  mutate(
    pred_var = rnorm(n = n(), mean = 0, sd = 1),
    obs = ifelse(set == "node", 0, 1)
  )


seq_2015_dates <- seq(as.Date("2015-01-01"), as.Date("2015-12-31"), by = 1)

sample_2015 <- rbind(sample_2015, mutate(dd, 
                                         year = 2015,
                                         date = sample(seq_2015_dates, n(), replace = TRUE))) %>%
  mutate(
    pred_var = rnorm(n = n(), mean = 0.5, sd = 1),
    obs = ifelse(set == "node", 0, 1)
  )
```
```{r}
#| label: tbl-samp_df
#| tbl-cap: Quick look at resulting data attributes

head(sample_2015)
```


## Cluster Distance  
A constructed covariate, or *cluster* variable is made to account for spatial correlation below the resolution of the mesh. This is done by calculating the distance to nearest locations with a simulated disease detection. (see, @Illian_2012 for discussion of this approach)
```{r warning=FALSE}
# 2014
split_factor <- as.factor(sample_2014$set)

NN <- as.data.frame(
  spatstat.geom::nndist(sample_2014[,c("x","y")], by = split_factor, k = 1) 
)["sim"]

sample_2014$NN <- round(NN$sim, 0)



# 2015
split_factor <- as.factor(sample_2015$set)

NN <- as.data.frame(
  spatstat.geom::nndist(sample_2015[,c("x","y")], by = split_factor, k = 1)
)["sim"]

sample_2015$NN <- round(NN$sim, 0)
```


### Combine Data  
Joining data from both 2014 and 2015 to a common data set.  
```{r message=FALSE, warning=FALSE}
combined_data = rbind(sample_2014, sample_2015)
```


## Phylodynamics   
A maximum clade credibility tree is loaded.  This tree includes actual disease genomic data but does not have location information.  After reading the phylogeny, phylodynamic analysis is performed using the **phylodyn** package (@Karcher_2017).  The resulting effective population size estimate (**Ne**) is matched to the data frame by date.  
```{r message=FALSE, warning=FALSE}
vsv_tree = read.tree(here("assets/mcc.tre"))

phylodynamics <- phylodynamic_process(vsv_tree)

head(phylodynamics$dyn)

combined_data <- append_nearest_date(combined_data, phylodynamics$dyn, "Ne")
```

```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=8}
#| label: fig-tree
#| fig-cap: "VSV Phylogeny"

plot(phylodynamics$tree)
```


```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=8}
#| label: fig-phylodyn
#| fig-cap: "VSV Effective Population size (Ne)."

plot(phylodynamics$Ne)
```


## Host Exposure  
Because the goal is to estimate disease intensity for a geographic area, not just individual point locations, it is necessary to identify the exposure rate of livestock (disease hosts) to VSV. That is, estimates for host population must be assigned to point locations.  Host population, represented here by the **host_pop** raster file, must be added to the data frame.  However, coordinates need to represent a geographic area, not just a point, thus a simple point extraction from the raster would be naive.  Instead, *natural neighborhoods* are defined based on node locations, then total host population is summed across this areal extent before adding the value to the node point associated with that neighborhood.  This is process is carried out using the *mesh_tessellation()* function below.      
```{r message=FALSE, warning=FALSE}
host_pop.r <- rast(here("assets/host_pop.tif"))

tessel_mesh <- vect(mesh_tessellation(mesh.dom))

crs(tessel_mesh) <- crs(domain)

# scale the population
tessel_mesh$host_pop <- extract(host_pop.r, tessel_mesh, fun="sum", na.rm=TRUE)[,"host_pop"]/10^3

host_pop <- as.data.frame(
  extract(tessel_mesh, combined_data[,c("x","y")])) %>% 
  distinct(id.y, .keep_all = TRUE) 

combined_data$host_pop <- host_pop$host_pop
```

```{r message=FALSE, warning=FALSE, warning=FALSE, fig.height=8, fig.width=8}
#| label: fig-exposure
#| fig-cap: "Plots depict the exposure estimation process with color indicating host abundance(darker tones indicate more hosts). Mesh overlying a continuous host abudnance raster grid (top left), removing mesh edges from plot and calculating natural neighborhoods around nodes using tesselation (top right), removing nodes to show total host abundance summed to neighborhood areas (bottom left), and showing assignment of neighborhood sums back to nodes (bottom right). "

plot_exposure_seq(mesh = mesh.dom, dd_nodes = dd, rast_data = host_pop.r, 
                  rast_scale = 10^3, xmin = -11500, xmax = -10500, ymin = 3500, ymax = 4500)
```

## Time Index  
Recoding year and date variables as integer based indices.
```{r message=FALSE, warning=FALSE}
combined_data$year_step <- as.integer(as.factor(combined_data$year))
combined_data$day_step <- as.integer(as.factor(combined_data$date))
```

## Split Predictor Variable  
Ensuring a difference is evident between years within the fictitious prediction variables, which are intended to represent environmental conditions.  
```{r message=FALSE, warning=FALSE}
combined_data <- combined_data %>%
  mutate(pred_var_14 = if_else(year == 2014, pred_var, NA),
         pred_var_15 = if_else(year == 2015, pred_var, NA))
```



## Projection  

### Random Field  
Coordinate locations in the data frame are matched or *projected* to the spatial mesh and a corresponding index is created to identify these locations during  spatial field estimation.
```{r message=FALSE, warning=FALSE}
k = length(unique(combined_data$year_step)) # number of years to model, each with a separate spatial field.

locs = cbind(combined_data$x, combined_data$y)

# Match locations in data frame to locations in mesh
A.pf = inla.spde.make.A(mesh.dom, 
                          alpha = 2,
                          loc=locs,
                          group = combined_data$year_step)

spde0 = inla.spde2.pcmatern(mesh.dom, alpha = 2,
                            prior.range=c(250, 0.01), # a 0.01 probability that auto correlation falls to zero at approximately 250km
                            prior.sigma=c(1, 0.01),
                            constr = TRUE)

# index spatial field to have k replicates (2 replicates, one for each year)
Field.pf = inla.spde.make.index("Field.pf", 
                               spde0$n.spde,
                               n.group=k)
```


### Organize Data
First tier of the spatial model is organized as a *list()* object.  This includes an intercept ($\beta_{\Lambda}$) and index for the spatial field ($\text{W}_{st}\}$).
```{r}
pf.lst = list(c(Field.pf, # index for spatial field
                list(intercept1 = 1)), # intercept, 
                list(year_step1 = combined_data[,"year_step"])) # time index for year (t)

pf.stk = inla.stack(data = list(Y = cbind(combined_data$obs, NA)), # creating a two bivariate matrix, 
                                A = list(A.pf, 1), # projection matrix
                          effects = pf.lst, # data frame 
                              tag = "pf.0") # label to pull data later

```

### Estimate Stack  
As in prior chunk, data for the second tier of the model is organized as a *list()*.  This includes an index (*Field.st*) to copy the spatial field from tier 1 ($\text{W}_{st}\}$) to tier 2, an intercept ($\beta_0$), the fictitious prediction variables (, *pred_var*, $\beta_bX_{bst}$), cluster variable ($\gamma_{clust}$), effective population size ($\gamma_{ne}$), and needed time indices (t).  
  
**NOTE:** Spatial fields are estimated for each year (two time steps), where as effective population size is daily.    
```{r message=FALSE, warning=FALSE}
A.est = inla.spde.make.A(mesh.dom, # projection as previously done
                         alpha = 2,
                         loc=locs,
                         group = combined_data$year_step)

Field.est = inla.spde.make.index("Field.est", # index spatial filed to copy W_st from tier 1
                                 spde0$n.spde,
                                 n.group=k)

est.lst = list(c(Field.est, # spatial field
                 list(intercept2 = 1)), # level specific intercept
                 list(pred_var_14 = combined_data[,"pred_var_14"], # fictitious fixed/linear effect
                      pred_var_15 = combined_data[,"pred_var_15"],
                      NN = combined_data[,"NN"], # cluster distances
                      Ne = combined_data[,"Ne"], # effective population size
                      year_step1 = combined_data[,"year_step"], # time indices
                      day_step = combined_data[,"day_step"]))

est.stk = inla.stack(data = list(Y = cbind(NA, combined_data$obs), # creating a two bivariate matrix
                                 e = combined_data$host_pop), # estimated host exposure
                                 A = list(A.est, 1), 
                           effects = est.lst,   
                               tag = "est.0")



joint.stk = inla.stack(pf.stk, est.stk) # join Tier 1 and Tier 2 to a common list object
```

## Model Formula  
The model is described fully in the manuscript but is provided below as a reference for comparison to code.

$$\begin{align}
   \text{log}(\Lambda_{st}) &= W_{st} \nonumber \\ \nonumber
    \quad \quad \Lambda_{st} &= \text{exp}\{\beta_{\Lambda} + \text{W}_{st}\}\\ \nonumber
    \text{W}_{st} &\overset{\textit{iid}}{\sim} \textit{N}(0, \text{Q}(\textit{r},\sigma)) \\ \nonumber
    \textit{r} &\sim \textit{Pr}(250, 0.01) \\ \nonumber
    \sigma &\sim \textit{Pr}(1, 0.01) \\ \nonumber
   \text{Y}_{st}|\lambda_{st}  &\sim \text{Poisson}(\mu_{st})  \nonumber \\ \nonumber
   \mu_{st} &= \textit{E}_{\textit{st}}\lambda_{st} \\ \nonumber
   log(\mu_{\textit{st}}) &= log(\textit{E}_{\textit{st}}) + log(\lambda_{\textit{st}}) \\ \nonumber
   log(\lambda_{\textit{st}}) &=  \beta_0 + \sum_{b=1}^{B} \beta_bX_{bst} + \gamma_{clust} + \gamma_{river} + \gamma_{ne} + \alpha \text{W}_{st}  \\  \nonumber
   \beta_b &= (\beta_{eddi}, \beta_{wet}, \beta_{temp}, \beta_{seas}, \beta_{carb}, \beta_{shrb}, \beta_{elev}, \beta_{ndvi}) \\ \nonumber
   \beta_b &\sim \textit{N}(1, 0.001) \\ \nonumber
   \gamma_{clust}, \gamma_{river} &\sim \textit{Pr}(1, 0.01) \\ \nonumber
   \gamma_{ne} &\sim \textit{Pr}(3, 0.01) \\ \nonumber
\end{align}$$

```{r message=FALSE, warning=FALSE}
pcprior1 = list(prec = list(prior="pc.prec", param = c(1, 0.01))) 
hc1 = list(theta = list(prior = 'normal', param = c(0, 10)))
ctr.g = list(model = 'iid', hyper = hc1)

Formula.1 = Y ~ -1 + intercept1 + # intercept tier 1
                     intercept2 + # intercept tier 2
          				f(Field.pf, # spatial index for tier 1 spatial field
          				  model=spde0,
          				  group = Field.pf.group, # group by year, k=2
                            control.group=ctr.g) +
                  f(Field.est, # create copy of tier 1 spatial field
                    copy = "Field.pf", 
          				  group = Field.est.group, 
                            fixed = FALSE,
                            hyper = hc1) + # IID between years
         				  f(NN,     # Cluster distance
                    model="rw1", # order-1 random walk
                    constr=TRUE, # center values
         				    replicate = year_step1, # estimate independently (IID) for each year
                    scale.model=TRUE, # internal scaling
         				    hyper=pcprior1) + #Penalizing complexity priors
                  f(day_step, Ne,   # Similar to above, but for effective population size
                    model="rw1", 
                    constr=TRUE,
                    scale.model=TRUE,
         				    hyper=pcprior1) +
                  pred_var_14 + pred_var_15 # linear effects 
```

## Run Model  
This light, demonstrative model took approximately 1hr to run.  The saved model is loaded below to examine outputs and results.
```{r eval=FALSE, message=FALSE, warning=FALSE}
theta1 = c(3.091600, 6.136307, -2.297550, -1.228109, -2.166505, -1.035271) # mean hyperparameters from initial run (speeds up runs)  

Model.demo <- inla(Formula.1, # modle formula
      				 num.threads = 12, # threads
      				 data = inla.stack.data(joint.stk),  # data organized above as list
      				 family = c("gaussian", "poisson"), # tier 1 and tier 2 distribution families
      				 verbose = TRUE,
      				 E = inla.stack.data(joint.stk)$e, # exposure
      				 control.fixed = list(prec = 1, prec.intercept=1), # proper intercept
      				 control.predictor = list(
      										A = inla.stack.A(joint.stk), 
      										compute = TRUE, # calculate fitted values
      										link = 1), # default link functions
      				 control.mode = list(restart = TRUE, theta = theta1), # hyperparameters to speed up run
      				 control.inla = list(strategy="adaptive", # fast efficient stratigy
      									 int.strategy = "eb"),
      				 control.compute=list(dic = FALSE, cpo = FALSE, waic = FALSE)) 

save(list=c("Model.demo"), file=here("assets/demo.RData"), version = 2) # save model
```

Load previously run model  
```{r message=FALSE, warning=FALSE}
load(here("assets/demo.RData"))
```

## Results

### Create Grid  
Creatin a point grid and raster to visulaize results.
```{r message=FALSE, warning=FALSE}
blank.r <- obs_stack.r[[1]]
blank.r[!is.na(blank.r)] <- 0

grid_pnts <- as.points(blank.r)
names(grid_pnts) <- "cell_value"


grid_coords <- grid_pnts %>%
  geom() %>%
  as.data.frame() 

Ap = inla.spde.make.A(mesh.dom, 
                      loc = cbind(grid_coords[,"x"], 
                                  grid_coords[,"y"]))
```

### Cluster Distance
```{r message=FALSE, warning=FALSE, warning=FALSE, fig.height=6, fig.width=8}
#| label: fig-clust
#| fig-cap: "Spatial clustering of simulated disease detections"

plot_cluster(Model.demo)
```

### Linear Coefficients  
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=6}
#| label: fig-fixed
#| fig-cap: "Fictitious linear effects constructed with different mean averages."
#| 
plot_fixed_marginals(Model.demo, c("pred_var_14", "pred_var_15"))
```

### Phylodynamic Effect
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=8}
#| label: fig-phylo_out
#| fig-cap: "Effective population size effect."
#| 
plot_phylodynamic(Model.demo, parameter = "day_step")
```

### Spatial Random Fields  
Extract estimates from model and rasterize them for visualization. 
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=5}
mrf_pf <- cbind(Model.demo$summary.random$Field.pf$mean, 
                Field.pf$Field.pf.group)

mrf_pf_v <- list()
mrf_pf_v  <- split(mrf_pf[,1], mrf_pf[,2])

grid_pnts$pf_2014 <- drop(Ap %*% mrf_pf_v[[1]]) 
grid_pnts$pf_2015 <- drop(Ap %*% mrf_pf_v[[2]]) 

pf_14.r <- rasterize(grid_pnts, 
                    blank.r, 
                    "pf_2014",
                    background = NA)

pf_15.r <- rasterize(grid_pnts, 
                    blank.r, 
                    "pf_2015",
                    background = NA)

pf_stack.r <-list(pf_14.r, pf_15.r)
```
```{r message=FALSE, warning=FALSE, fig.height=6, fig.width=7}
#| label: fig-rfs
#| fig-cap: "Year specific Spatial Random Fields estimated by demonstarted model.  Similiar to those used to used for simulation."

plot_rast_panels(pf_stack.r, domain, join_sims, rev(cubehelix(50)), plot_pnts = FALSE, center=TRUE)
```