---
title: "Scan Statistics for Space-Time Cluster Detection"
subtitle: "Licentiate presentation"
author: "Benjamin Allévius"
date: "2018-10-25"
output:
  xaringan::moon_reader:
    css: ["default", "su", "tamu-fonts", "custom.css"]
    chakra: libs/remark-latest.min.js
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
bibliography: bibliography.bib
---
background-image: url(./figures/piechart.png)
background-position: 40% 90%
background-size: 40% 40%;

```{r setup, include=FALSE, echo=FALSE, cache=FALSE}
options(htmltools.dir.version = FALSE)

library(RefManageR)
BibOptions(check.entries = FALSE, 
           bib.style = "authoryear", 
           cite.style = 'alphabetic', 
           style = "markdown",
           hyperlink = FALSE, 
           dashed = FALSE)
myBib <- ReadBib("./bibliography.bib", check = FALSE)
```


.pull-left[
## Outline

1. Context: disease surveillance

1. Problem: local outbreaks

1. Solution: scan statistics
]
.pull-right[
## Papers

.small[
1. "Prospective Detection of Outbreaks" (with Michael Höhle),
   _Handbook of Infectious Disease Data Analysis_, 2018.
   
1. "An unconditional space–time scan statistic for ZIP-distributed data"  (with Michael Höhle),
   _Scandinavian Journal of Statistics_, in press.
   
1. "scanstatistics: space-time anomaly detection using scan statistics",
   _The Journal of Open Source Software_, 2018.
   
1. "On the precision matrix of an irregularly sampled AR(1) process", technical
   report, 2018.
]]

???

Since I'm assuming most of you are unfamiliar with scan statistics

* I will frame this cluster detection problem in the context of disease surveillance

Where the cluster detection part corresponds to the detection of local disease outbreaks

To get there, I will first give you a more general introduction to disease surveillance
and some of the statistical methods used in this field

* This part of the talk will culminate in an introduction to scan statistics,
  which then leads me to talk about a particular scan statistic that I 
  have worked on

* Following that, I will speak a bit about some software I have written which
  provides a number of scan statistics in what I hope is a user-friendly format

* Lastly, I will speak a bit about some technical and computational results for 
  a particular time series processes

Disease surveillance is just one aspect of public health


---
background-image: url(https://foodsafetynewsfullservice.marlersites.com/files/2014/03/public-health406x250.jpg)
background-position: 95% 15%
background-size: cover;

## Public health

.large[

Broadly, deals with 

* Health _promotion_

* Disease _prevention_

Governmental agencies and NGOs, like:

* Folkhälsomyndigheten (Sweden)

* Centers for Disease Control and Prevention (US)

Here, focus is on _infectious_ diseases. One part of disease _prevention_ is 
then the detection and thwarting of disease __outbreaks__

> A disease outbreak is the occurrence of cases of disease in excess of what would normally be expected in a defined community, geographical area or season.

]

???

* Prevention one aspect, promoting health another.

* Carried out by governmental agencies and NGOs
    * __End users of the methods we discuss__

* Epidemiology, biostatistics and health services 

* Public Health Agency of Sweden, Centers for Disease Control and Prevention

To detect outbreaks, we must look for them, actively or passively. This 
"looking" is what we call disease surveillance

[Outbreak definition from WHO](http://www.searo.who.int/topics/disease_outbreaks/en/)


[Image credit: "Local Health Departments Gain Accreditation From New Independent Board" by Food Safety News](https://www.foodsafetynews.com/2014/03/local-health-departments-gain-accreditation-from-new-independent-board/)

---
background-image: url(./figures/CASE_small.png)
background-position: 80% 85%
background-size: 50% 50%;

## Disease surveillance


.large[
* _Systematic_ and _continuous_ collection of data

* Surveillance _systems_

* _Support_ decision-making by:
    * Automating
    
    * Summarizing
    
    * Visualizing

Example: [CASE](https://case.folkhalsomyndigheten.se/)
]



???

* WHO: Surveillance can be defined as 'ongoing systematic collection, collation, analysis and interpretation of data and the dissemination of information to those who __need to know__ in order that action may be taken' 

* Automate:
    * Data collection/centralization
    
    * Detection algorithms: indicate whether outbreak is likely
    
    * Report-making: __visual guides__ to help non-statisticians to make sense
      of results
    
* Examples: CASE (Computer Assisted Search for Epidemics)
    * Image from user interface
    * Select diseases to monitor, algorithms to use, other parameters
    
What are these algorithms?

<!-- * Collection: laboratories/physicians report cases to local health authorities -->
<!--     * In Sweden, each county (21 of them) has a medical officer who is in  -->
<!--       charge of regional infectious disease prevention and control. -->
<!--     * Officers + hospitals + laboratories connected to SmiNet database. -->
        
    
[Quote credit: "Integrated disease surveillance and response" by World Health Organization](https://www.who.int/countries/eth/areas/surveillance/en/)

---

# Univariate surveillance methods

.large[
Data:

* $y_1, y_2, \ldots, y_t$, where $y_t$ is the most recent observation.

Idea:

* Compare observed $y_t$ to quantiles from suitable distribution

* Which data?
    
* Which distribution?

Features:

* Low integer number of cases

* Possible tends and seasonality

]

???

Forget about cluster detection for a minute and focus on just detecting an outbreak

* Time index: daily, weekly, monthly

How to determine if $y_t$ is indicative of an outbreak?

* Compare most recent observation to a __threshold__, perhaps based on some 
  assumptions about how the data is distributed

Features making surveillance a challenge:

* Low integer number of reported cases

* Possible long-term trends and seasonality

* Under-reporting and reporting delays
    * Not seeing the full extent of number of cases
    
    
Still you have to start somewhere, and perhaps start with simple methods

<!-- * Potential past outbreaks -->

<!-- * Unclear case definition (IMD -> microbial investigation) -->

<!-- * Spikes in incidence could be due to increased testing -->

<!-- * Underestimation: under-ascertainment and under-reporting -->

<!-- * Reporting delays -->

---

## Univariate methods

.left-code[
Examples:

* EARS C1 (Early Aberration Reporting System, CDC)

* Farrington algorithm

Limitations:

* Aggregation hides outbreak signals?

* Follow multiple time series at once?



]
.right-plot[
![](./figures/farrington.png)
]


???

Data: Invasive meningococcal disease (IMD), 2002-2008, Germany
    
* EARS: 
    * Assume normality
    * Use 7 previous data points to calculate sample mean and sd
    * Less worry about seasonality and trends
    * Simplest case: Use mean + 3 sd as threshold
    * Better: use upper limit of prediction interval
    * Continuous distribution for discrete data - can do better

* Farrington: 
    * Use count data distribution (overdispersed Poisson GLM)
    * Handles trend and seasonality through regression
       * Use fixed window of historic data

* Univariate data typically aggregated
    * For example looking at total sum of cases in a country
    * Aggregation could hide outbreak signals in individual time series
    * Maybe want to look at multiple geographic regions instead
    * Or monitor secondary sources (e.g. OTC medicine sales)

---

# Multivariate surveillance methods

.left-code[

Observations $\{y_{it}\}$, where $i=1,\ldots,N$ and $t=1,\ldots,T$. 

Could apply univariate methods in parallel, perhaps hierarchically at 
different levels of aggregation.
* Could then use count distributions.
* The pragmatic choice?

Problems:

* Multiple testing 
* Miss dependencies

Look instead at joint distribution of $\boldsymbol{y}_{t} = (y_{1,t},\ldots,y_{N,t})^T$

* Can we model these dependencies properly?

]
.right-plot[
![](./figures/multivariate.png)
]


???

* $i=1,\ldots,N$ (regions, secondary data sources, etc)
* $t=1,\ldots,T$ (days, weeks, etc)
    
* Hierarchical example (Höhle et al 2009): 
    * Pragmatic choice: multivariate changepoint detection methods assume 
      continuous distributions; not realistic for low count time series.
    * Keep in-control alarm rates comparable between levels by having higher
      thresholds (lower significance levels) at lower levels
    
* Problem with univariate methods applied in parallell:
    * More false alarms for same significance level $\alpha$
    
    * Lower $\alpha$ $\Rightarrow$ outbreaks harder to detect 
        * Sample size per statistic is unchanged
        
        * So bigger effect size needed


---

## Multivariate methods

.large[

* Consider vector of counts/statistics to be multivariate normal. Examples:
    * Hotelling's $T^2$ statistic: 
    $$T_t^2 = (\boldsymbol{y}_t - \hat{\boldsymbol{\mu}})' \hat{\boldsymbol{\Sigma}}^{-1}
    (\boldsymbol{y}_t - \hat{\boldsymbol{\mu}})$$.
    
    * CUSUM applied to $T^2$: $S_0 = 0$ and $$S_{t+1} = \max\{0, S_t + T_t - \omega_t\}.$$
    
    * $T^2$ applied to univariate statistics ("vector accumulation" as opposed to
    "scalar reduction")
]
    
???

* Hotelling (1947): estimate $\mu$ and $\Sigma$ and calculate quadratic form with
  inverse covariance matrix
    * Squared Mahalanobis distance: unitless and scale-free squared distance of 
      observed vector to its mean when correlations are accounted for.
      
* CUSUM (cumulative sum control chart): Take square-root of $T^2$ and accumulate over time if sufficiently
  large value (calibrate threshold to FPR).
  
* E.g. univariate EMWA statistics

* Working with counts: normal dist not good
      
* Problem with Hotelling and similar methods is that they don't specify 
  _which_ variables are responsible for exceeding the threshold
    * That is, _where_ is the potential outbreak located, 
      if data is spatio-temporal
      
    * Thus more investigation is needed
      
* This is where _scan statistics_ come into the picture


---
background-image: url(./figures/us_map.jpg)
background-position: 50% 100%
background-size: 75%
class: center, top

# Scan Statistics 
#### for
## Space-Time Cluster Detection


???

They help us not only detect _that_ an outbreak may be ongoing, but also
  _where_ it may be located
    
__Image__ alludes to procedure: 

* I __envision__ it as __scanning__ over a map with a __magnifying glass__

* For captured regions, we calculate a __statistic__
    
* In a minute I will get into what we do with the statistics so calculated

[Image credit: Jin Chen, Robert E Roth, Adam T Naito, Eugene J Lengerich & Alan M MacEachren. 2008. "Geovisual analytics to enhance spatial scan statistic interpretation: an analysis of U.S. cervical cancer mortality". International Journal of Health Geographics 7:57 ](https://doi.org/10.1186/1476-072X-7-57)

---

## Scan statistics: setup

.large[
Data: Area-referenced discrete time series $\{y_{it}\}$

Cluster detection is the goal:

* Find a subset $Z$ of regions in which an outbreak is suspected to have emerged 
  during the last $D \in \{1,2,\ldots\}$ time periods. 
  
* The spatial _zone_ $Z$ and the _duration_ $D$ together define a 
  space-time window $W$, which stretches backwards in time from the present
  (in prospective surveillance).
  
* Zones $Z$ pre-specified or data-determined
  
_Many_ possible windows $W$ to consider. Test each one separately?
]

???

* Counts from a fixed region like a county, municipaly or zip code,
  rather than the coordinates of each event

* Window stretches backwards in time from most recent
  time period
    * Same duration in all regions - for computational feasibility
    
* We typically want the zones - which are what the __magnifying glass__ I 
  mentioned earlier captures at a given instance - to be geographically compact.
    * Increasing KNN or connected subsets thereof
    
    * Or more data-driven approaches
    
---


## Scan statistics: definition, computation, testing

.large[
No. Instead, consider only _most anomalous_ cluster.

* Formulate hypotheses: $H_0$ says there is no outbreak anywhere, 
  $H_1$ says there is some (unknown) window $W$ where outbreak is emerging.

* Calculate a statistic ([loglihood](https://goo.gl/vDET5K) ratio) for each 
  potential space-time outbreak window.

* Look at the maximum—this is the scan statistic. The corresponding window
  is the _most likely cluster_.
  
* How to test?
    * Distribution of scan statistic is generally unknown; exact and asymptotic 
      $p$-values hard to come by.
    * But, can simulate new data from null and redo calculations to obtain critical values.
    * Permutation tests also sometimes possible.

]



???

* __Calculate statistic__ under some assumption of how the 
  __counts are distributed__

* How do we know if the value we got is indicative of an outbreak?


---

## Example: Kulldorff's prospective scan statistic

.left-code[
Assumptions:

* Counts distributed as $Y_{it} \sim \text{Poisson}(q_{it} b_{it})$.

* Null: $q_{it} = q$ for all $i$ and $t$.

* Alternative: for $(i,t) \in W$, $q_{it} = q_{\text{in}} > q_{\text{out}}$.

Calculation:


$$\lambda_W = 
\log [
\left( \frac{Y_{\text{in}}}{B_{\text{in}}} \right)^{\!\! Y_{\text{in}}}
\!\!
\left( \frac{Y - Y_{\text{in}}}{Y - B_{\text{in}}} \right)^{\!\! Y - Y_{\text{in}}}
\!\!\!\!\!\!\!\!\! ]
\boldsymbol{1}_\left\{ Y_{\text{in}} > B_{\text{in}} \right\}$$

Do this calculation for each $W$ under consideration.

Does the Poisson assumption hold for disease case data?

]
.right-plot[
![](./figures/first_score_plot.png)
]

???

* __Canonical__ example dating back to two papers written in 1997 and 2001 by
  Martin Kulldorff which really set the modern field of scan statistics going.

* $b_{it}$ is an _expected count_ or _baseline_, proportional to the population
  at risk in Kulldorff's model
  
* Baselines sum to total observed count in original model

* $q$ is the _relative risk_
    * This type of multiplicative increase in the mean is how outbreaks are 
      typically modelled for scan statistics
      
* $Y_{\text{in}}$ is the sum of counts inside $W$

* $B_{\text{in}}$ is the sum of baseline parameters inside $W$
      
* $Y$ is the total count for the entire study area


---

## Example: Kulldorff's prospective scan statistic

.left-code[
Assumptions:

* Counts distributed as $Y_{it} \sim \text{Poisson}(q_{it} b_{it})$.

* Null: $q_{it} = q$ for all $i$ and $t$.

* Alternative: for $(i,t) \in W$, $q_{it} = q_{\text{in}} > q_{\text{out}}$.

Calculation:


$$\lambda_W = 
\log 
\left( \frac{Y_{\text{in}}}{B_{\text{in}}} \right)^{\!\! Y_{\text{in}}}
\!\!
\left( \frac{Y - Y_{\text{in}}}{Y - B_{\text{in}}} \right)^{\!\! Y - Y_{\text{in}}}
\!\!\!\!\!\!\!\!\!
\boldsymbol{1}_\left\{ Y_{\text{in}} > B_{\text{in}} \right\}$$

Do this calculation for each $W$ under consideration.

Does the Poisson assumption hold for disease case data?

]
.right-plot[
![](./figures/scoreplot.gif)
]

???

* Good: 
    * Theoretical guarantees regarding test power 
    * That has been elusive under other model assumptions
    * "Individually most powerful test"
    

* Problem with Kulldorff's Poisson statistic:
    * Empirically, disease case data often shows greater variation than Poisson 
      distribution allows
    
    * Another issue is that there are often more zero counts in real-world disease
      data than suggested by a Poisson distribution
    
* This last points brings us to the work done by me and Michael on a scan statistic
  based on the ZIP distribution

`r NoCite(myBib, c("Kulldorff1997", "Kulldorff2001"))`
---

## Unconditional ZIP scan statistic (Allévius & Höhle, 2018)
.small[

Model case counts by zero-inflated Poisson distribution: 
  $$\mathbb{P}(Y = 0) = p + (1-p) e^{-\mu}$$
  $$~~~~~~~~~~~~~~~~~~~\mathbb{P}(Y = y) = (1-p) e^{-\mu} \mu^y / y!, \quad y > 0$$

* Allow these _baseline parameters_ $~p$ and $\mu$ to vary spatially and
  temporally (and possibly by other covariates), and fit to historical data 
  using regression.

For outbreaks, assume multiplicative increase in Poisson mean: $\mu_{it} \to q_W \mu_{it}$
for $(i,t) \in W$ and $q_W > 1$.

* Hotspot outbreak model: same $q_W$ for all counts in outbreak window $W$
* Estimate $q_W$ by EM algorithm, by first initiating $\hat{q}_W = \hat{q}_W^{(0)}$ and then 
  for $k \geq 1$ alternating between 
    * __Expectation step__: for all $(i,t) \in W$ for which $y_{it} = 0$, set
      $$\hat{\delta}_{it}^{(k)} = \frac{p_{it}}{p_{it} + (1-p_{it}) \exp\left( -\hat{q}_W^{(k-1)} \mu_{it} \right)}$$
    * __Maximization step__: set 
      $$\hat{q}_W^{(k)} = \max\left\{ 1, \frac{\sum_{(i,t) \in W} y_{it}}{\sum_{(i,t) \in W} (1-\hat{\delta}_it^{(k)} \mu_{it}} \right\}$$
      
Increase $k$ by 1 at the end of each iteration and continue until convergence
of the likelihood.
]


???

* Mix of point mass at zero and Poisson distribution.

* Allow baseline parameters to vary; fit by regression.

* Separates count modelling from cluster detection task


* Hotspot: 
    * Same increase in all regions, at all timepoints
    * Good model for food-borne outbreaks
    * For emerging person-to-person outbreaks, increasing $q$ would be better
      (higher detection power, but more computational effort)
    * Side with simplicity and computational speed; hotspot model still has some
      power to detect gradually increasing outbreaks.
      
* Derive EM algorithm by first looking at the __complete data likelihood__,
  i.e. when we __know which zeros are from the point mass__.
    * $\delta_{it} = 1$ if obs. at $(i,t)$ is a structural zero, else 0.
    * Complete data likelihood has nice expression if all structural zero 
      indicators are known, and $q_W$ can be easily estimated.
    * For the incomplete data likelihood, when indicators are unknown, we can
      replace indicators by their expectations given a value for $q_W$.
    * Given these expecations, we can estimate $q_W$ anew.
      
    
`r NoCite(myBib, c("Allevius2018a"))`
      
---

## Results

.left-code[

#### Evaluation

* Simulated data

* Detection timeliness

* Spatial accuracy

#### Contribution

* Excess zeros + model flexibility

* Unconditional analysis

* Fast open-source implementation available on CRAN.

]
.right-plot[
![](./figures/zip_res.png)
]

???

* Two main components to performance:
    * How soon do we get an alarm signal?
    * How accurate is the detected cluster?

* Precision: what fraction of the detected regions truly are affected by an outbreak?
  $TP / (TP + FP)$

* Recall: what fraction of true outbreak regions did we capture?
  $TP / (TP + FN)$
  
* Separate baseline parameter estimation from outbreak parameter estimation

* Previous work: 
    * Only spatial
    * Uses population at risk and same $p$ for all regions
    
* Unconditional: don't condition on total count (Tango)

[Image credit: Allévius, Benjamin, and Michael Höhle. 2018, in press. “An unconditional space-time
scan statistic for ZIP-distributed data.” Scandinavian Journal of Statistics.](https://doi.org/10.1111/sjos.12341)

---

## scanstatistics R package (Allévius, 2018)

.large[

* Open-source R package currently providing 6 scan statistics plus related
  methods
  
* Partially overlapping complement to [SaTScan™](https://www.satscan.org/) `r NoCite(myBib, "satscan")`

* Upcoming: multivariate methods, new interface
]

.center[
<img src="./figures/nm_map.png" width=70% align="middle">
]

???

* Statistics:
    * Conditional Poisson and space-time permutation ScanStat by Kulldorff (in SaTScan)
    * Unconditional Poisson, ZIP, Negative Binomial (with increasing $q$)
    * Bayesian SS
* Methods:
    * KNN clusters, connected subsets
    * Heuristic score for regions for visualization
* Upcoming:
    * Fast (Penalized) Subset Scan
        * Score aggregation
        * Subset aggregation
    * Multivariate Bayesian Scan Statistic
    * Simplified interface a lá `caret::train`
    
`r NoCite(myBib, c("Allevius2018b"))`
    
---
background-image: url(./figures/ar1sim.png)
background-position: 50% 100%
background-size: 25% 25%;

## Irregularly sampled AR(1) processes

.large[

* Consider autoregressive process $X_t = \rho X_{t-1} + \epsilon_t$, $t=1,2,\ldots$

* Under assumptions on $\rho$, $\{\epsilon_t\}$ and distribution of starting value $X_0$, consecutively sampled process (i.e. at times $j,j+1,j+2,\ldots,j+m-1$) is stationary with easily expressed inverse 
  covariance matrix. 

    * Sample this process at times $t_1 < t_2 < \ldots < t_m$ instead, where $t_{\ell+1} - t_{\ell} \geq 1$, . What does 
      this matrix look like?

* All implemented in the R package `irregulAR1`.
]


???

Motivating example:

* Sampled calendar years (integers) according to some distribution, and needed
  to evaluate a density for a process observed at those years, and also simulate
  values of the process at itermediate years. These computations needed to be 
  fast.
  
* The results I'm about to describe are usable in missing data scenario as well.

Now, if the process had been sampled at consecutive time points, then my needs
would be met by a stationary Gaussian AR(1) process since its inverse covariance matrix is 
tridiagonal with easily expressible elements

* I was unable to find a corresponding expression of this matrix when the process
  was sampled at irregular frequencies, which I then set out to prove for myself

It turned out that the irregularly sampled process also had a nice tridiagnal
shape with elements that could easily be expressed using the distance in time
between one observation and the next

* This then led to some fast algorithms for unconditional and conditional 
  density evaluation and random number generation, which I implement in the
  irregulAR1 package.

---

# Future

* __Multivariate__ methods: space + time + multiple variables
    * Subset scan framework: ability to prioritize or exclude regions before
      finding highest-scoring subset `r Citep(myBib, c("Neill2012", "McFowland2013", "Speakman2016"))`
    
    * Bayesian multivariate methods: incorporate prior information and distinguish
      between multiple event types (e.g. diseases) `r Citep(myBib, c("Neill2010", "Neill2011"))`

* Integration of scan statistics with __machine learning__ methods
    * Support Vector Machines and Penalized Fast Subset Scan `r Citep(myBib, "Fitzpatrick2017")`
    
    * Gaussian Process Subset Scan `r Citep(myBib, "GPSS2018")`

* User-friendly __software__: keep end users in mind
    * [SaTScan™](https://www.satscan.org/) `r Citep(myBib, "satscan")`
    
    * [SpatialEpiApp](https://paula-moraga.github.io/software/) `r Citep(myBib, "SpatialEpiApp")`
    
    * scanstatistics?

???

Machine learning integration:

* Combine scan statistics with e.g. SVM:
    * Alternate between PFSS scan statistic that finds non-compact clusters and 
      Support Vector Machine, modify point-wise penalty by distance to boundary

* Also some promising Gaussian Process approaches that can be extended to a 
  __multivariate__ setting, another path forward I think is promising, seeing
  as public health authorities monitor many diseases at once.
  
    
* Software:
    * Usable by non-statisticians (automate as much as possible)
    * Visualizations: maps, heat-maps,  zoom in on flagged time series, etc
    * SpatialEpiApp

---

## References
.smallest[
```{r, results='asis', echo=FALSE}
RefManageR::PrintBibliography(myBib)
```
]
???
