This vignette demonstrates basic usage of **surveil** for public health surveillance. **surveil** leverages the principles of Bayesian inference and Markov chain Monte Carlo (MCMC) (MacKay 2003; Gelman et al. 2014) to infer population risk of disease or death given time series data consisting of case counts and population at risk. Models were built using the Stan modeling language, but users only need to be familiar with the R language.

The package also contains special methods for age-standardization, printing and plotting model results, and for measuring and visualizing health inequalities. For age-standardization see `vignette("age-standardization")`

. For discussion and demonstration analysis see Donegan, Hughes, and Lee (2022).

Surveillance data minimally contain case counts, reliable population at risk estimates, and a discrete time period variable. They also may include one or more grouping variables, such as race-ethnicity.Time periods should consist of equally spaced intervals.

This vignette analyzes age-specific (ages 50-79) colorectal cancer incidence data by race-ethnicity, year, and Texas MSA, obtained through CDC Wonder. The race-ethnicity grouping includes (non-Hispanic) Black, (non-Hispanic) White, and Hispanic, and the MSAs include those centered on the cities of Austin, Dallas, Houston, and San Antonio.

```
head(msa) %>%
kable(booktabs = TRUE,
caption = "Glimpse of colorectal cancer incidence data (CDC Wonder)")
```

Year | Race | MSA | Count | Population |
---|---|---|---|---|

1999 | Black or African American | Austin-Round Rock, TX | 28 | 14421 |

2000 | Black or African American | Austin-Round Rock, TX | 16 | 15215 |

2001 | Black or African American | Austin-Round Rock, TX | 22 | 16000 |

2002 | Black or African American | Austin-Round Rock, TX | 24 | 16694 |

2003 | Black or African American | Austin-Round Rock, TX | 34 | 17513 |

2004 | Black or African American | Austin-Round Rock, TX | 26 | 18429 |

**surveil**’s model fitting function, `stan_rw`

, requires that the user provide a `data.frame`

with specific column names. There must be one column named `Count`

containing case counts, and another column named `Population`

, containing the sizes of the populations at risk. The user must provide the name of the column containing the time period, and may optionally provide a grouping factor. For the MSA data printed above, the grouping column is Race and the time column is Year.

We will demonstrate using aggregated CRC cases across Texas’s top four MSAs. The `msa`

data from CDC Wonder already has the necessary format (column names and contents), but these data are dis-aggregated by MSA. So for this analysis, we first group the data by year and race, and then combine cases across MSAs.

The following code chunk aggregates the data using the `dplyr`

package:

```
tx.msa <- msa %>%
group_by(Year, Race) %>%
summarise(Count = sum(Count),
Population = sum(Population))
```

The following code provides a glimpse of the aggregated data (Table 2):

```
head(tx.msa) %>%
kable(booktabs = TRUE,
caption = "Glimpse of aggregated Texas metropolitan CRC cases, by race and year")
```

Year | Race | Count | Population |
---|---|---|---|

1999 | Black or African American | 471 | 270430 |

1999 | Hispanic | 361 | 377471 |

1999 | White | 2231 | 1654251 |

2000 | Black or African American | 455 | 283280 |

2000 | Hispanic | 460 | 405562 |

2000 | White | 2343 | 1704425 |

The base **surveil** model is specified as follows. The Poisson model is used as the likelihood: the probability of observing a given number of cases, \(y_t\), conditional on a given level of risk, \(e^{\phi_t}\), and known population at risk, \(p_t\), is: \[y_t \sim \text{Pois}(p_t \cdot e^{\phi_t})\] where \(t\) indexes the time period.

Alternatively, the binomial model is available: \[y_t \sim \text{Binom}(p_t \cdot g^{-1}(\phi_t))\] where \(g\) is the logit function and \(g^{-1}(x) = \frac{exp(x)}{1 + exp(x)}\) (the inverse-logit function). The Poisson model is often preferred for ‘rare’ events (such as rates below .01), otherwise the binomial model is generally more appropriate. The remainder of this vignette will proceed using the Poisson model only.

Next, we build a model for the log-rates, \({\phi_t}\) (for the binomial model, the rates are logit-transformed, rather than log-transformed). The first-difference prior states that our expectation for the log-rate at any time is its previous value, and we assign a Gaussian probability distribution to deviations from the previous value (Clayton 1996). This is also known as the random-walk prior: \[\phi_t \sim \text{Gau}(\phi_{t-1}, \tau^2)\] This places higher probability on a smooth trend through time, specifically implying that underlying disease risk tends to have less variation than crude incidence.

The log-risk for time \(t=1\) has no previous value to anchor its expectation; thus, we assign a prior probability distribution directly to \(\phi_1\). For this prior, **surveil** uses a Gaussian distribution. The scale parameter, \(\tau\), also requires a prior distribution, and again **surveil** uses a Gaussian model.

For multiple time series, **surveil** allows users to add a correlation structure to the model. This allows our inferences about each population to be mutually informed by inferences about all other observed populations. Note that in many cases, the base model specification described above will be entirely sufficient.

The log-rates for \(k\) populations, \(\boldsymbol \phi_t\), are assigned a multivariate Gaussian model (Brandt and Williams 2007): \[\boldsymbol \phi_t \sim \text{Gau}(\boldsymbol \phi_{t-1}, \boldsymbol \Sigma),\] where \(\boldsymbol \Sigma\) is a \(k \times k\) covariance matrix.

The covariance matrix can be decomposed into a diagonal matrix containing scale parameters for each variable, \(\boldsymbol \Delta = diag(\tau_1,\dots \tau_k)\), and a symmetric correlation matrix, \(\boldsymbol \Omega\) (Stan Development Team 2021): \[\boldsymbol \Sigma = \boldsymbol \Delta \boldsymbol \Omega \boldsymbol \Delta\] When the correlation structure is added to the model, then a prior distribution is also required for the correlation matrix. **surveil** uses the LKJ model, which has a single shape parameter, \(\eta\) (Stan Development Team 2021). If \(\eta=1\), the LKJ model will place uniform prior probability on any \(k \times k\) correlation matrix; as \(\eta\) increases from one, it expresses ever greater skepticism towards large correlations. When \(\eta <1\), the LKJ model becomes ‘concave’—expressing skepticism towards correlations of zero.

The time series model is fit by passing surveillance data to the `stan_rw`

function. Here, `Year`

and `Race`

indicate the appropriate time and grouping columns in the `tx.msa`

data frame.

```
fit <- stan_rw(tx.msa,
time = Year,
group = Race,
iter = 1500,
chains = 2 #, for speed only; use default chains=4
)
#> Distribution: normal
#> Distribution: normal
#> [1] "Setting normal prior(s) for eta_1: "
#> location scale
#> -6 5
#> [1] "\nSetting half-normal prior for sigma: "
#> location scale
#> 0 1
#>
#> SAMPLING FOR MODEL 'RW' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 1500 [ 0%] (Warmup)
#> Chain 1: Iteration: 751 / 1500 [ 50%] (Sampling)
#> Chain 1: Iteration: 1500 / 1500 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.258363 seconds (Warm-up)
#> Chain 1: 0.231846 seconds (Sampling)
#> Chain 1: 0.490209 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'RW' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 1.3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 1500 [ 0%] (Warmup)
#> Chain 2: Iteration: 751 / 1500 [ 50%] (Sampling)
#> Chain 2: Iteration: 1500 / 1500 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.316586 seconds (Warm-up)
#> Chain 2: 0.245818 seconds (Sampling)
#> Chain 2: 0.562404 seconds (Total)
#> Chain 2:
```

If we wanted to add a correlation structure to the model, we would add `cor = TRUE`

(as opposed to the default, `cor = FALSE`

). To speed things up, we could take advantage of parallel processing using the `cores`

argument (e.g., add `cores = 4`

) to run on 4 cores simultaneously. For age-standardization, see `vignette("age-standardization")`

.

In this case, Stan did not issue any warnings; this means that all the R-hat values are near to 1 and the Monte Carlo standard errors are reasonably small. We can verify that chains have converged (and R-hats equal about 1):

```
rstan::stan_rhat(fit$samples)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
```

Other useful functions for examining MCMC diagnostics are `rstan::stan_ess(fit$samples)`

and `rstan::stan_mcse(fit$samples)`

. For an introduction to MCMC analysis with **surveil**, including MCMC diagnostics, see the vignette `vignette("surveil-mcmc")`

.

If we call `plot`

on a **surveil** model, we get a `ggplot`

object depicting risk estimates with 95% credible intervals:

```
plot(fit, scale = 100e3, base_size = 11)
#> Plotted rates are per 100,000
```

Crude incidence rates are also plotted as points.

The `plot`

method has a number of options to control its appearance. For example, the `base_size`

argument controls the size of labels. The size of the points for the crude rates can be adjusted using `size`

, and `size = 0`

removes them altogether. We can also use `ggplot`

to add custom modifications:

```
fig <- plot(fit, scale = 100e3, base_size = 11, size = 0)
#> Plotted rates are per 100,000
fig +
theme(legend.position = "right") +
labs(title = "CRC incidence per 100,000",
subtitle = "Texas MSAs, 50-79 y.o.")
```

The plot method has a `style`

argument that controls how the probability distribution is represented. The default, `style = "mean_qi"`

, shows the mean of the posterior distribution of the risk at each time period with a shaded 95% credible interval (as above). The alternative, `style = "lines"`

, plots MCMC samples from the joint probability distribution across all time periods:

```
plot(fit, scale = 100e3, base_size = 11, style = "lines")
#> Plotted rates are per 100,000
```

By default, `M = 250`

samples are plotted. The `style`

option is available for all of the **surveil** plot methods.

The `print`

method will print the estimates with 95% credible intervals to the console:

```
print(fit, scale = 100e3)
#> Summary of surveil model results
#> Time periods: 19
#> Grouping variable: Race
#> Correlation matrix: FALSE
#> time Race mean lwr_2.5 upr_97.5
#> 1 1999 Black or African American 170.01246 157.98969 181.81119
#> 2 2000 Black or African American 166.28707 156.30463 176.64155
#> 3 2001 Black or African American 168.21979 158.49241 177.76281
#> 4 2002 Black or African American 169.06382 159.68532 178.70871
#> 5 2003 Black or African American 166.99554 157.25642 176.69423
#> 6 2004 Black or African American 166.64963 157.32718 177.46076
#> 7 2005 Black or African American 158.89944 150.02064 167.73513
#> 8 2006 Black or African American 154.73518 146.33568 163.52624
#> 9 2007 Black or African American 152.62056 144.37622 161.03935
#> 10 2008 Black or African American 149.77311 141.70215 157.95968
#> 11 2009 Black or African American 143.57548 135.83268 151.64576
#> 12 2010 Black or African American 138.89103 131.98652 146.06689
#> 13 2011 Black or African American 131.01169 123.97772 138.05325
#> 14 2012 Black or African American 125.46654 118.15800 132.54584
#> 15 2013 Black or African American 125.11527 118.26142 131.73641
#> 16 2014 Black or African American 126.24511 120.02506 133.05571
#> 17 2015 Black or African American 122.46925 116.03108 128.70314
#> 18 2016 Black or African American 121.86017 115.45556 128.75496
#> 19 2017 Black or African American 122.36429 115.73665 129.49943
#> 20 1999 Hispanic 101.53401 94.78429 108.00777
#> 21 2000 Hispanic 104.12644 98.14130 110.67825
#> 22 2001 Hispanic 102.51482 96.98401 108.18937
#> 23 2002 Hispanic 101.46026 96.32235 107.24877
#> 24 2003 Hispanic 100.49400 95.17999 105.66433
#> 25 2004 Hispanic 100.86395 96.07851 106.57438
#> 26 2005 Hispanic 98.57103 93.91013 103.23694
#> 27 2006 Hispanic 96.87815 92.21519 102.07269
#> 28 2007 Hispanic 94.18909 89.90937 98.89390
#> 29 2008 Hispanic 90.49879 85.47552 94.94830
#> 30 2009 Hispanic 88.70993 84.27980 92.86947
#> 31 2010 Hispanic 87.91644 83.67085 91.84159
#> 32 2011 Hispanic 87.33348 83.09833 91.42645
#> 33 2012 Hispanic 87.87272 83.73342 92.03119
#> 34 2013 Hispanic 87.30706 83.66046 91.21094
#> 35 2014 Hispanic 85.84189 82.05655 89.54929
#> 36 2015 Hispanic 85.83022 82.11391 89.58177
#> 37 2016 Hispanic 84.79289 80.84719 88.74490
#> 38 2017 Hispanic 85.89517 81.82401 90.05726
#> 39 1999 White 135.16203 130.18088 140.06740
#> 40 2000 White 136.50468 131.70298 141.59965
#> 41 2001 White 134.14095 129.62203 138.95257
#> 42 2002 White 130.17990 125.97038 134.31126
#> 43 2003 White 127.41164 123.17908 131.84235
#> 44 2004 White 120.03485 115.90386 124.17778
#> 45 2005 White 115.06502 111.06505 119.02083
#> 46 2006 White 109.69405 105.90416 113.68644
#> 47 2007 White 108.20948 104.66659 112.09869
#> 48 2008 White 103.39353 99.75091 107.18740
#> 49 2009 White 98.86180 95.36861 102.41334
#> 50 2010 White 96.13789 92.89184 99.39938
#> 51 2011 White 93.65297 90.52846 96.79312
#> 52 2012 White 91.05922 87.90101 94.23520
#> 53 2013 White 90.09742 87.09160 93.24469
#> 54 2014 White 91.41577 88.27352 94.64249
#> 55 2015 White 93.06601 89.97668 96.37288
#> 56 2016 White 88.99003 86.04799 91.98494
#> 57 2017 White 92.25455 89.05752 95.52087
```

This information is also stored in a data frame, `fit$summary`

:

```
head(fit$summary)
#> time mean lwr_2.5 upr_97.5 Race Year Count
#> 1 1999 0.001700125 0.001579897 0.001818112 Black or African American 1999 471
#> 2 2000 0.001662871 0.001563046 0.001766415 Black or African American 2000 455
#> 3 2001 0.001682198 0.001584924 0.001777628 Black or African American 2001 505
#> 4 2002 0.001690638 0.001596853 0.001787087 Black or African American 2002 539
#> 5 2003 0.001669955 0.001572564 0.001766942 Black or African American 2003 546
#> 6 2004 0.001666496 0.001573272 0.001774608 Black or African American 2004 602
#> Population Crude
#> 1 270430 0.001741671
#> 2 283280 0.001606185
#> 3 298287 0.001693000
#> 4 313133 0.001721313
#> 5 329481 0.001657152
#> 6 346886 0.001735440
```

**surveil** provides a number of functions and methods for measuring health inequalities.

A selection of complementary pairwise inequality measures can be calculated using the `group_diff`

function. The function requires a fitted **surveil** model and character strings corresponding, respectively, to the target population (indicating which group is the target of our inference, typically the overburdened or disadvantaged group), and the reference population. You can also use `group_diff`

to compare two age-stratified populations with age-standardized rates (for details, see `vignette("age-standardization")`

and `?group_diff`

).

It returns probability distributions and summary statements for the following quantities, where \(A\) is the incidence rate for the advantaged group, \(D\) is the incidence rate of the disadvantaged group, and \(P_d\) is the size of the population at risk for the disadvantaged group.

Concept | Formula |
---|---|

Rate Ratio (RR) | \(\frac{D}{A}\) |

Rate Difference (RD) | \(D - A\) |

Excess Cases (EC) | \((D-A) \times P_d\) |

Proportion Attributable Risk (PAR) | \(\frac{D-A}{D}\) |

Notice that the PAR is simply the rate difference expressed as a fraction of total risk; it indicates the fraction of risk in the target population that would have been removed had the target rate been equal to the reference rate (Menvielle, Kulhánová, and Machenbach 2019). Each of these quantities (RD, RR, PAR, EC) provide a different perspective on the magnitude of a health disparity.

To calculate all of these measures for two groups in our data, we call `group_diff`

on our fitted model:

```
gd <- group_diff(fit, target = "Black or African American", reference = "White")
print(gd, scale = 100e3)
#> Summary of Pairwise Inequality
#> Target group: Black or African American
#> Reference group: White
#> Time periods observed: 19
#> Rate scale: per 100,000
#> Cumulative excess cases (EC): 3,209 [2971, 3441]
#> Cumulative EC as a fraction of group risk (PAR): 0.27 [0.25, 0.28]
#> time Rate RD PAR RR EC
#> 1999 170 35 0.20 1.3 94
#> 2000 166 30 0.18 1.2 84
#> 2001 168 34 0.20 1.3 102
#> 2002 169 39 0.23 1.3 122
#> 2003 167 40 0.24 1.3 130
#> 2004 167 47 0.28 1.4 162
#> 2005 159 44 0.28 1.4 160
#> 2006 155 45 0.29 1.4 180
#> 2007 153 44 0.29 1.4 186
#> 2008 150 46 0.31 1.4 204
#> 2009 144 45 0.31 1.5 208
#> 2010 139 43 0.31 1.4 209
#> 2011 131 37 0.28 1.4 192
#> 2012 125 34 0.27 1.4 185
#> 2013 125 35 0.28 1.4 197
#> 2014 126 35 0.28 1.4 205
#> 2015 122 29 0.24 1.3 180
#> 2016 122 33 0.27 1.4 209
#> 2017 122 30 0.25 1.3 198
```

All of the **surveil** plotting and printing methods provide an option to scale rates by a custom value. By setting `scale = 100e3`

(100,000), the RD is printed as cases per 100,000. Note that none of the other inequality measures (PAR, RR, EC) are ever impacted by this choice.

The plot method for `surveil_diff`

produces one time series `ggplot`

each for RD, PAR, and EC. The means of the probability distributions for each measure are plotted as lines, while the shading indicates a 95% credible interval:

```
plot(gd, scale = 100e3)
#> Rate differences (RD) are per 100,000 at risk
```

If we wanted to replace the plot of the PAR with one of the RR, we would set the `PAR`

option to `FALSE`

:

```
plot(gd, scale = 100e3, PAR = FALSE)
#> Rate differences (RD) are per 100,000 at risk
```

Pairwise measures are important, but they cannot provide a summary of inequality across multiple socially situated groups. Theil’s T is an entropy-based inequality index with many favorable qualities, including that it naturally accommodates complex grouping structures (Theil 1972; Conceição and Galbraith 2000; Conceição and Ferreira 2000).

Theil’s T measures the extent to which certain populations are overburdened by disease, meaning precisely that the proportion of cases accounted for by a particular group, \(\omega_j\), is higher than the proportion of the population constituted by that same group, \(\eta_j\). With \(k\) groups, Theil’s index is \[T = \sum_{j=1}^k \omega_j \big[ log(\omega_j / \eta_j) \big].\] This is zero when case shares equal population shares and it increases monotonically as the two diverge for any group. Theil’s T is thus a weighted mean of log-ratios of case shares to population shares, where each log-ratio (which we may describe as a raw inequality score) is weighted by its share of total cases.

Theil’s T can be computed from a fitted **surveil** model, the only requirement is that the model includes multiple groups (through the `group`

argument):

```
Ts <- theil(fit)
print(Ts)
#> Summary of Theil's Inequality Index
#> Groups:
#> Time periods observed: 19
#> Theil's T (times 100) with 95% credible intervals
#> time Theil .lower .upper
#> 1999 0.938 0.615 1.279
#> 2000 0.799 0.540 1.090
#> 2001 0.902 0.625 1.206
#> 2002 0.980 0.703 1.295
#> 2003 0.994 0.707 1.327
#> 2004 1.067 0.756 1.454
#> 2005 0.996 0.701 1.335
#> 2006 1.049 0.735 1.414
#> 2007 1.101 0.775 1.469
#> 2008 1.250 0.884 1.674
#> 2009 1.209 0.833 1.622
#> 2010 1.137 0.814 1.495
#> 2011 0.917 0.611 1.286
#> 2012 0.778 0.485 1.128
#> 2013 0.815 0.542 1.131
#> 2014 0.858 0.590 1.184
#> 2015 0.680 0.439 0.951
#> 2016 0.802 0.523 1.163
#> 2017 0.711 0.457 1.049
```

The probability distribution for Theil’s T can be summarized visualy using the `"lines"`

style plot or by plotting estimates with shaded 95% credible intervals:

`plot(Ts)`

While the minimum of Theil’s index is always zero, the maximum value varies with the structure of the population under observation. The index is useful for comparisons such as monitoring change over time, and should generally not be used as a indication of the absolute level of inequality.

The index also has interesting extensions; for example, given disease data for a nested population structure—such as racial-ethnic groups within states—Theil’s index can provide a measure of geographic inequality across states (between-state inequality), and social inequality within states (within-state inequality) (Conceição, Galbraith, and Bradford 2001). For details, see `?theil`

.

Brandt, P, and JT Williams. 2007. *Multiple Time Series Models*. Sage.

Clayton, DG. 1996. “Generalized Linear Mixed Models.” In *Markov Chain Monte Carlo in Practice: Interdisciplinary Statistics*, edited by WR Gilks, S Richardson, and DJ Spiegelhalter, 275–302. CRC Press.

Conceição, Pedro, and Pedro Ferreira. 2000. “The Young Person’s Guide to the Theil Index: Suggesting Intuitive Interpretations and Exploring Analytical Applications.” *University of Texas Inequality Project (UTIP)*. https://utip.gov.utexas.edu/papers.html.

Conceição, Pedro, and James K. Galbraith. 2000. “Constructing Long and Dense Time Series of Inequality Using the Theil Index.” *Eastern Economic Journal* 26 (1): 61–74.

Conceição, Pedro, James K Galbraith, and Peter Bradford. 2001. “The Theil Index in Sequences of Nested and Hierarchic Grouping Structures: Implications for the Measurement of Inequality Through Time, with Data Aggregated at Different Levels of Industrial Classification.” *Eastern Economic Journal* 27 (4): 491–514.

Donegan, Connor, Amy E. Hughes, and Simon J. Craddock Lee. 2022. “Colorectal Cancer Incidence, Inequality, and Prevention Priorities in Urban Texas: Surveillance Study with the ‘Surveil’ Software Pakcage.” *JMIR Public Health & Surveillance* 8 (8): e34589. https://doi.org/10.2196/34589.

Gelman, Andrew, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin. 2014. *Bayesian Data Analysis*. Third. CRC Press.

MacKay, David M. 2003. *Information Theory, Inference, and Learning Algorithms*. Cambridge University Press.

Menvielle, Gwenn, Kulhánová, and Johan P. Machenbach. 2019. “Assessing the Impact of a Public Health Intervention to Reduce Social Inequalities in Cancer.” In *Reducing Social Inequalities in Cancer: Evidence and Priorities for Research*, edited by Salvatore Vaccarella, Joannie Lortet-Tieulent, Rodolfo Saracci, David I. Conway, Kurt Straif, and Christopher P. Wild, 185–92. Geneva, Switzerland: WHO Press.

Stan Development Team. 2021. *Stan Modeling Language Users Guide and Reference Manual, 2.28*. https://mc-stan.org.

Theil, Henry. 1972. *Statistical Decomposition Analysis*. Amsterdam, The Netherlands; London, UK: North-Holland Publishing Company.