Public Health Surveillance in High-Dimensions with Supervised Learning

CHAPTER 1

INTRODUCTION

1.1. Introduction

Public health surveillance is used to monitor health data for changes, which

is a special case of the more general problem where counts (or rates) of events

are monitored for changes. This is a problem that commonly occurs in dis-

ciplines such as manufacturing, economics, engineering, as well as in public

health. In manufacturing, the events might simply be surface flaws or incor-

rectly soldered components. In public health surveillance, the events might be

disease incidents, but other events such as accidents or preventive treatments

(or many others) might also be of interest. Even though this problem occurs

in many disciplines, the discrete nature of the health-related data as well as its

common reporting practices pose significant challenges to traditional monitoring

approaches.

The prevalent theme in public health surveillance has traditionally been out-

break detection, where health-related data is monitored for changes in disease

incidence rates. Increments in disease incidence rates present a threat to the

public (Joner Jr et al., 2008). Hence, it is important to have in place a surveil-

lance system that rapidly detects change increments and that allows one to con-

duct in-depth investigations to determine their root causes. Despite the focus in

outbreak detection (Buckeridge et al., 2005), public health surveillance is also

important for the following purposes: (1) estimating the magnitude of the public

health problem, (2) determining the geographic distribution of illness, (3) por-

traying the natural history of disease, (4) detecting epidemics, (5) generating

hypothesis and stimulating research, (6) evaluating control measures, (7) moni-

toring changes in infectious agents, (8) detecting changes in health practice, and

(9) facilitating planning.

1.2. Motivation

When there are few covariates, detecting a change in the overall rate of an

event is a relatively easy task. Consider the case of spatial surveillance where

a geographical code is the only covariate, and the goal is to detect a change

to the event rate in one or more geographic regions. Traditional approaches

such as the spatial scan statistic (Kulldorff, 1997) and Tango’s method (Tango,

1995) have been widely used for this problem. These methods use counts of

health-related events, which are aggregated over geographical regions, to trigger

a signal whenever one or more geographical regions have larger counts than ex-

pected. Nevertheless, an aggregated count of events is not sensitive to changes

in multiple dimensions. Such localized changes require a more subtle analysis,

particularly, when the number of covariates and/or levels associated with covari-

ates are more than a few. In addition, it is not easy to incorporate spatial and

other covariate information into traditional surveillance methods. Spatial infor-

mation would require counts in specific regions (such as counties, cites, or cen-

sus tracts), and these are usually pre-defined. Such region aggregates can greatly

increase the number of covariates that need to be monitored. Also, pre-defined

regions restrict the flexibility of the method to detect changes at higher and lower

resolutions. Furthermore, other covariates are difficult to add. Patient-level co-

variate information (such as age) would also need to be aggregated over similar

time/space intervals. For example, to monitor asthma events using a weekly

count, one might need to define covariates for each of several age groups. Then,

age covariates might need to be crossed with spatial regions to obtain counts

in every combination of age group and region. This generates a rapid expan-

sion to the number of covariates that limits such an approach. Separate monitors

for each covariate potentially increase false alarms and, even with multivariate

methods (such as multivariate control charts), the large number of covariates

increases the dimensionality of the problem and further decreases the sensitivity.

With modern public health data, it is also possible to have numerous demo-

graphic and geographic covariates associated with individual disease incidences.

The objective is to detect a change in any region of this high-dimensional space.

Frequently, surveillance methods are limited to monitoring disease incidence

rates in space or time. When they can directly incorporate more covariates into

the analysis, they still struggle to detect changes that occur in a relatively small

region or set of regions of the high-dimensional space. For instance, consider

a surveillance system that monitors individual diagnoses observed at the county

level. Few methods available to public health surveillance are able to detect an

emerging disease cluster in one county when there is a large number of coun-

ties simultaneously being monitored. Very frequently, the change gets diluted in

high dimensions.

On the other hand, monitoring disease incidence rates that are subject to sea-

sonality is another challenge that face public health surveillance methods. Many

traditional methods assume a constant disease incidence rate over time, which

can lead to a large number of false alarms when the disease incidence rate of in-

terest is indeed a function of time. It is important to have a method that is flexible

enough to allow for detecting changes in incidence rates only when the number

of incidences observed significantly exceeds the number of incidences expected

for that particular season. Equally important is the case where the population for

the study region is unevenly distributed. A large number of false alarms can arise

if one, for instance, compares the number of observed incidences within a geo-

graphical region with average population size against the number of incidences

observed in heavily populated geographical region. Clearly, this does not indi-

cate a change in the disease incidence rate, and the method needs to adequately

address this issue.

1.3. Original contributions

This research focuses on developing a series of methodologies for the early

detection of changes in disease incidence rates with a focus on the problem with

high-dimensional covariates, where traditional methods are challenged. It uses

all available covariates found in modern public health data, which provide nu-

merous demographic and geographic covariates that generate categories or re-

gions in which events may be organized. The proposed methods approach multi-

dimensional surveillance through a transform of public health surveillance to

supervised learning. The first method, labeled Feature Selection (FS) method,

defines a signal for changes in disease incidence rates using feature selection.

The second method, labeled Percentile method, proposes an alternative to the

FS method with simpler implementation requirements. The Percentile method

defines a signal for changes in disease incidence rates using a percentile of class

probability estimates from tree-based ensembles. Finally, the third method, la-

beled Rule-Based Feature Value Selection (RBFVS) method, can be used as a

follow-up step to either the FS method or the Percentile method. Once either

one of these methods identifies a change in the disease incidence rate of interest,

the RBFVS method can be used to define its specific location as defined by the

multi-dimensional interaction among the covariates involved in the cluster. This

third method uses a combination of tree-based classifiers with association rules

to extract the important values, categorical or numerical, which are responsible

for the observed change in the disease incidence rate.

1.4. Organization of the research

Chapter 2 presents both a formal definition of public health surveillance as a

point process and a summary of traditional public health surveillance methods.

Chapter 3 presents the Feature Selection (FS) method along with a variety of ap-

plications of the method on simulated data including cases with strong covariate

interactions, seasonality, and unevenly distributed populations. Chapter 3 also

provides the power curves for the FS method. Chapter 4 presents the descrip-

tion for the Percentile method. In addition, several applications of the Percentile

method on simulated data and a comparison to the FS method are shown in

Chapter 4. Chapter 5 presents the Rule-Based Feature Value Selection (RBFVS)

method, its results for a variety of simulated examples, and a comparison to Arti-

ficial Contrast Ensemble method. Lastly, Chapter 6 includes concluding remarks

and potential extensions to this work.

CHAPTER 2

BACKGROUND

2.1. Introduction

Public health data often exists in the form of counts (or rates) of events that

are monitored over time and compared to historical or baseline data. Less fre-

quently, public health data can be obtained in the case-event format, where each

row of data represents an individual disease incidence (or case) along with its

associated covariate information. The baseline is generally historical data of the

same type of event, but, in some cases, the baseline could be generated from a

different type of event.

Case-event data can be used to define the occurrence of disease incidences

in public health surveillance as a point process (O’Sullivan and Unwin, 2006).

Let Pbe the marked point process defined by a set of disease incidences that are

irregularly distributed within geographical region A. Note that the term “disease

incidences” refers to incidences associated with a particular disease of interest.

Let Zrepresent the binary response variable associated with marked point pro-

cess P. Since Pis the type of process where only the occurrence of events are

observed, let Z(x, y)=1for {x, y} ⊂ Arepresent the occurrence of a disease

incidence in geographical coordinates {x, y}. Coordinates {x, y}often represent

the latitude and longitude coordinates for the patient’s home address, although,

they could also represent the planar projection of these latitude and longitude

coordinates. If available, the latter would be preferred in order to remove the

effects of the curvature of the Earth in any distance calculations the surveillance

method of choice might require. The marks for point process Pbecome the

remaining set of patient-level covariates recorded with the disease incidence.

Again, potential patient-level covariates can represent the patient’s gender, age

group, race, etc.

A summary of the most commonly used public health surveillance methods

is discussed in the remainder of this chapter. Some of the methods discussed in-

clude: control charts, Tango’s tests for “generalized” and “focused” disease clus-

ters, scan statistics, generalized linear mixed models (GLMM), and wavelets.

Note that some of the information in this chapter might be summarized in the

following chapters.

2.2. Control Charts

Many methods that have proven successful in fields such as manufacturing

have been adapted to public health surveillance. For instance, control charting

tools from statistical process control (SPC) have proven effective in temporal

public health surveillance (Woodall, 2006). While the common objective in sta-

tistical process control and public health surveillance is to detect changes in the

process (or in disease incidence rates) as early and accurately as possible, there

are inherent differences between public health surveillance and statistical pro-

cess control. Tsui et al. (2008), Shmueli and Burkom (2007), and Woodall et al.

(2006) comment on some of the differences between control charting in public

health surveillance and standard SPC. Some of these differences include:

(a) the use of attribute data is more prevalent in health-related applications,

(b) public health surveillance baseline (Phase I) data is generally non-stationary

and correlated,

control charting,

(d) there are numerous types of outbreak patterns related to a disease as opposed

to the often clear definitions and patterns of shifts in an industrial process,

(e) the definition of a natural disease outbreak is not clear,

(f) risk-adjusted probability estimates are used in the construction of health-

related control charts, as opposed to assuming that there is a constant in-

control probability of failure in well-controlled conditions of industrial ap-

plications,

(g) public health surveillance does not emphasize the sampling of only a portion

of the output of a process at periodic intervals,

(h) it is not possible to adjust a health-related process to try to return it quickly

to in-control performance, and

(i) a health-related control chart might continue to provide alarms after its first

alarm.

Some control charts are suitable for case-event data, while others are suit-

able for aggregated counts. Control charts, nevertheless, are one of the public

health surveillance approaches that are most unreliable when monitoring for an

increase in disease incidence rates that occurs, for instance, within a specific ge-

ographical area and subpopulation. When using control charts for monitoring

such highly localized changes, one would have to monitor individually multiple

spatial-covariate regions where these regions would need to be pre-defined, and

each region would provide a disease incidence rate. Separate control charts for

each event type are well-known to increase false alarms and more importantly,

they do not incorporate the valuable information available in the correlations

among covariates.

2.2.1. CUSUM

The cumulative sum (CUSUM) chart is one of the most common tools from

statistical process control adapted to monitor adverse events in public health

surveillance (Morton et al., 2001; Woodall, 2006). The original CUSUM chart

was first proposed by Page (1954) as a visual technique for plotting the cumula-

tive sum of deviations from the mean, placing emphasis on keeping the process

on aim rather than allowing it to drift between the upper and lower control lim-

its (Lucas, 1985). An out-of-control signal in a CUSUM chart indicates that an

action should be taken to prevent the adverse effects from exceeding the target

(Limaye et al., 2008). In public health surveillance, this translates to monitoring

disease incidence rates in order to prevent a disease from spreading through the

population across a large region.

The cumulative sum control chart (CUSUM) directly incorporates all the in-

formation in the sequence of sample values by plotting the cumulative sums of

the deviations of the sample values from a target value (Montgomery, 2005). Let

µ0be the target for the process mean, then, the CUSUM control chart is formed

by plotting

Ci=

j=1

¯xj−µ0,(2.1)

where Ciis called the cumulative sum up to and including the ith sample and ¯xj

is the average measurement for the jth sample. Due to the combination of the in-

formation from several samples, cumulative sum charts are more effective than

Shewhart control charts for detecting small shifts in the process mean (or dis-

ease incidence rate) (Montgomery, 2005). Indeed, they are particularly effective

for monitoring individuals. Because tighter control is available with CUSUM

schemes, CUSUM charts also place more emphasis on keeping the process on

aim rather than allowing it to drift between the control chart limits (Lucas, 1985).

Note that if a process (or incidence rate) is under control at µ0, the cumulative

sum in Equation (2.1) is a random walk with mean equal to zero. However, if

the mean shifts upward to some value µ1> µ0, then, an upward or positive drift

will develop in the cumulative sum Ci. On the other hand, if the mean shifts

downward to µ1< µ0, then, a downward or negative drift in Ciwill develop.

Therefore, if a significant trend develops in the plotted points, whether upward or

downward, this trend should be considered as an indication that the process mean

has shifted and that a search for assignable causes is in order (Montgomery,

2005).

Health-related CUSUM charts are typically one-sided, with the part corre-

sponding to a decrease in the target many times not being considered. The upper

CUSUM is designed to detect worsening performance, while the lower CUSUM

is designed to detect improved performance (Woodall, 2006). By including both

upper and lower CUSUM limits, both an increase and decrease in disease inci-

dence rates can be detected (Limaye et al., 2008).

As a general rule, a CUSUM chart triggers an out-of-control signal when-

ever Ciexceeds a reference value kwhere kis determined by an acceptable

disease incidence rate. After selecting k, reference tables are used to choose the

parameter hbased on the acceptable frequency of false out-of-control signals

(Lucas, 1985). Evaluations on CUSUM control schemes are based on average

run lengths (ARL) performance. The in-control ARL is the average number of

samples taken before an out-of-control signal is obtained. Two-sided schemes to

detect either an increase or a decrease in disease incidence rates are obtained by

simultaneously implementing two one-sided schemes.

A standard CUSUM control chart will have a starting value S0= 0, while

a fast initial response (FIR) CUSUM will have a positive starting value. Lucas

(1985) recommends a head start (S0) value approximately equal to h/2to allow

the CUSUM to signal quickly if the process is off aim at start-up. For a process

running at the target level, the head start will soon go to zero. It is, nonetheless,

critical to view the choice of S0in terms of the αand βerrors of the equivalent

sequential probability ratio tests (SPRT) introduced by Wald (1947). As in stan-

dard hypothesis testing, αis the probability of declaring the process off target

when it is indeed running at the target level, and βis the probability of declar-

ing the process on target when it is not. Using a head-start value approximately

equal to h/2allows equating the αand βerrors for the initial SPRT.

2.2.1.1. Poisson CUSUM

Poisson CUSUM charts can be used in public healthcare whenever incidences

are reported as aggregated counts over a given sampling interval. For a Poisson

CUSUM used to detect an increase in incidences, the CUSUM statistic can be

written as

Si= max (0, Yi−k+Si−1),(2.2)

where kwill be chosen to be between the acceptable incidence rate (µa) and

the mean level of incidences (µd) that the CUSUM scheme is expected to de-

tect quickly. Note both µaand µdrepresent a mean number of incidences per

sampling interval. Even though the ideal value for µais zero, a µavalue of zero

is usually not used in the design formulas. If µais zero, then, the CUSUM is

designed with h= 1 and k= 0, which implies that any incidence will be con-

sidered an out-of-control signal. When incidence levels are low and a decrease

in the incidence level indicates a reduction in the incidence rate, every incidence

must be followed up to find and remove assignable causes. In practice, µais

often chosen near to the current mean level (Lucas, 1985). The reference value

for the Poisson CUSUM (kp) is typically selected to approximate the following

formula:

kp=µd−µa

ln(µd)−ln(µa).(2.3)

2.2.1.2. Binomial CUSUM

The Binomial CUSUM monitors a binomial random variable T, where Tis

the total number of defectives in a sample of size n. In public health surveillance,

the problem can be translated to nhospital visits in a given interval where only

Tincidences correspond to the disease of interest. To detect an increase in a

disease incidence rate (p), the binomial CUSUM control statistic can be written

as:

Sj= max (0, Sj−1+ (Tj−nγ)) , j = 1,2,..., (2.4)

where nγ is the reference value. To determine the value of γ, it is necessary

to specify a value p1> p0, which represents an out-of-control value of pthat

should be detected quickly (Reynolds and Stoumbos, 2000). For a given in-

control value p0and a given out-of-control value p1, the constants r1and r2can

be obtained by:

r1=−ln 1−p1

1−p0and r2= ln p1(1 −p0)

p0(1 −p1).(2.5)

Then, the appropriate choice for γbecomes:

γ=r1

.(2.6)

For simplicity in plotting and for ease in developing statistical properties, it will

usually be convenient to set γ= 1/m where mis an integer. In general, if γ=

1/m, where mis an integer, then the possible values for Tj−nγ will be integer

multiples of 1/m. Also, the possible values of Sjwill also be integer multiples

of 1/m. If p0and p1are small, then, a slight change in p1will be sufficient to

make γ= 1/m where mis an integer. In most applications, however, the precise

specification of p1will not be critical (Reynolds and Stoumbos, 2000).

Note that after the inspection of item j, the value of the CUSUM statistic

Sjis calculated by adding the increment (Tj−nγ)to Sj−1if Sj−1≥0, but if

Sj−1<0, then, there is a reset to zero before (Tj−nγ) is added. This def-

inition, thus, allows Sj−1to assume negative values (Reynolds and Stoumbos,

2000). Although this definition of the CUSUM statistic has been used previously

(Reynolds and Stoumbos, 1998, 1999), it is slightly different from the traditional

binomial CUSUM statistic, where negative values are not allowed. The reason

for using the definition in Equation (2.4) is that certain expressions for the sta-

tistical properties of the binomial CUSUM chart depend on knowing the value

of Sjwhen it drops below zero.

2.2.1.3. Bernoulli CUSUM

As opposed to the binomial CUSUM, the Bernoulli CUSUM does not require

to wait until a sample of nitems has been processed. The Bernoulli CUSUM

chart is based directly on the individual observations without using a summary

statistic based on any grouping of the items into samples (Reynolds and Stoum-

bos, 2000). Thus, Bernoulli CUSUM charts in public health surveillance can

be used to monitor individual incidences without the need to consider sampling

periods or incidences from diseases that are not of interest. The non-inspection

periods are ignored for purposes of defining the sequence of observations used

by the Bernoulli CUSUM chart. To detect an increase in a disease incidence rate

(p), the Bernoulli CUSUM control statistic can be expressed as:

Bk= max (0, Bk−1)+(Xk−γ), k = 1,2,..., (2.7)

where the reference value is γgiven by Equation (2.6). Note the definition in

Equation (2.7) for the Bernoulli CUSUM statistic (Bk) implies that certain ex-

pressions for the statistical properties of the Bernoulli CUSUM chart also depend

on knowing the value of Bkwhen it drops below zero.

Note that a binomial CUSUM with n= 1 is equivalent to a Bernoulli

CUSUM. Also, there exits almost no difference between the detection abili-

ties of the Bernoulli and binomial CUSUM charts, except for large increases

in p. Moreover, the binomial CUSUM chart will perform almost as well as the

Bernoulli CUSUM chart for small values of n. However, if implementing a bi-

nomial CUSUM chart is considered to be more convenient than implementing a

Bernoulli CUSUM, then, using the binomial CUSUM chart will only be reason-

able if nis small (Reynolds and Stoumbos, 2000).

2.2.1.4. MCUSUM

Multivariate control charts allow to simultaneously monitor for multiple co-

variates. At the same time, they allow to incorporate correlation information

among covariates. To detect a shift in the mean vector of an independent m-

variate normal distribution, Woodall and Ncube (1985) suggest operating m

one-sided (or two-sided) CUSUM schemes simultaneously and evaluate the per-

formance of multivariate scheme in terms of the average run length (ARL) of

the collection of schemes (Crosier, 1988). Suppose that for a case or group

of cases, mcharacteristics are observed and given by mindependent normal

random covariates Xn= (X1n , X2n , . . . , Xmn )0,n= 1,2,···. Each mean vec-

tor {Xn}=µn= (µ1n , . . . , µmn )0,n= 1,2, . . ., is assumed to have known

variance-covariance matrix Σ. The run length of a one-sided CUSUM scheme

for the ith covariate is given by

N(i) = min (n:Si,n ≥hi),(2.8)

where

Si= max (0, Si,n−1+Xin −Ki), n = 1,2,..., (2.9)

0≤Si,0< hi, and Ki≥0. For a two-sided CUSUM scheme, the run length for

ith covariate is given by

N(i) = min (n:Si,n ≥hior Ti,n ≤ −hi),(2.10)

where

Ti,n = min (0, Ti,n−1+Xin +Ki), n = 1,2, . . . (2.11)

and −hi< Ti,0≤0. The reference value Kidepends on the shift in the mean of

the ith covariate that needs to be detected quickly. For the MCUSUM scheme,

the run length depends on the run lengths of the individual covariates:

N= min {N(1), . . . , N(m)}.(2.12)

Thus, interpreting out-of-control signals using a MCUSUM is simple: any co-

variate corresponding to a signaling univariate CUSUM is considered to be out

of control. Nevertheless, a problem that arises with the MCUSUM and with al-

most any other any multivariate control scheme, is the fact that as the number

of covariates being monitored increases, the average run length performance to

detect a specified shift in the mean of these covariates also increases because the

shift is “diluted” in the multidimensional space of the covariates (Montgomery,

2005). Moreover, the assumption of a multivariate normal distribution in health-

related data sets where event counts may not adequately transform to normality,

particularly for events of low frequency, is typically inadequate.

2.2.2. Hotelling T2

The most well known multivariate control chart is based on Hotelling’s T2

statistic (Hotelling, 1947). The T2test statistic is optimal for testing the null

hypothesis that the mean vector of a multivariate normal distribution is equal to

a constant vector against the alternative hypothesis that the mean vector is not

equal to that constant vector (Montgomery, 2005). The Hotelling T2test statistic

is sensitive to shifts in the variance as well as to shifts in the mean.

Suppose a preliminary analysis of msamples of size nis taken when the

process is assumed to be in control (i.e., base disease incidence rate). Let xijk

be the ith observation on the jth characteristic in the kth sample. When mon-

itoring changes in disease incidence rates, xijk will often be a binary covariate

where xijk = 1 indicates an incidence and xijk = 0 indicates no incidence.

Now, suppose vector Sis an unbiased estimate of Σwhen there are no dis-

ease clusters, and vector ¯

¯x is the corresponding approximation of mean vector µ

(Montgomery, 2005). Then, the T2test statistic becomes

T2=n(¯x −¯

¯x)0S−1(¯x −¯

¯x).(2.13)

Phase I control limits are given by

UCL =p(m−1)(n−1)

mn −m−p+ 1Fα,p,mn−m−p+1

LCL = 0,

(2.14)

while Phase II control limits are given by

UCL =p(m+ 1)(n−1)

mn −m−p+ 1Fα,p,mn−m−p+1

LCL = 0.

(2.15)

When a large number of samples is used to estimate µand Σ, it is customary

to use UCL =χ2

α,p as the upper control limits in both Phase I and Phase II

monitoring. In the univariate setting, Phase I and Phase II limits will nearly

coincide with m≥20. However, a larger number of preliminary samples (m≥

50) is often required in the multivariate setting (Lowry and Montgomery, 1995).

Similarly to the MCUSUM, problems with the Hotelling T2are sometimes

rooted on the assumption of the normality of the underlying population, where

event counts may not adequately transform to normality. Other issues that com-

plicate the application of Hotelling’s T2to health-related data include: (a) an in-

sensitivity to small changes in the underlying population, (b) the departure from

the underlying assumptions, (c) the discrete nature of the data, (d) the structured

nature of shifts in disease incidence rates (only a few covariates are involved in

changes in disease incidence), and (e) the unavailability of time-homogeneous

samples (Sonesson and Bock, 2003).

Nevertheless, the main difficulty that arises with any multivariate control

chart is how to identify which covariate or set of covariates is responsible for

an out-of-control signal. Alternatives to address this issue with the T2chart

include:

(a) plotting univariate ¯xcharts on the individual covariates,

(b) plotting univariate ¯xcharts on the individual covariates using Bonferroni-

type control limits (Alt, 1985),

(d) using control charts based on the principal components of the original co-

variates (Jackson, 1980), and

(e) decomposing the T2statistic into components that reflect the contribution of

each individual covariate (Runger et al., 1996).

2.2.3. EWMA

The Exponentially Moving Average (EWMA) control chart was introduced

by Roberts (1959). This control chart gives the most recent samples the greatest

weight, and all previous samples weights decreasing in geometric progression

from the most recent back to the first. Let the test statistic Zidenote the geomet-

ric moving average at time i,i= 0,1,2,3, . . .. Frequently, Z0=µ0where µ0

is the ordinate of the central line of control chart (Roberts, 1959). Other times,

it makes sense to set Z0= ¯x0where ¯x0is the average of preliminary data. The

EWMA test statistic is, thus, given by

Z0=µ0or Z0= ¯x0

Zi= (1 −λ)Zi−1+λ¯xi, i > 0.

(2.16)

Notice Ziis located a fraction 0≤λ≤1of the distance from Zi−lto ¯xion their

connecting straight line. Also, notice Zican be rewritten as

Zi=λ

i−1

j=0

(1 −λ)jxi−j+ (1 −λ)iZ0,(2.17)

where the weights λ(1 −λ)jdecrease geometrically with the age of the sample

mean (Montgomery, 2005).

If the ¯x’s are independent with common standard deviation σ¯x, the standard

deviation of Zis

σZi=rλ

2−λh1−(1 −λ)2iiσ¯x,(2.18)

and for a two-sided EWMA scheme, an out-of control signal is generated when-

ever Ziexceeds the control limits given by

UCL =Z0+ 3σZ

LCL =Z0−3σZ.

(2.19)

2.2.3.1. MEWMA

In order to detect changes in disease incidence, a number of covariates need

to be simultaneously monitored. These covariates obviously increase the dimen-

sionality of the problem. Sonesson and Bock (2003) propose the multivariate

extension of the EWMA (MEWMA) as a more desirable alternative for health-

related data. The MEWMA uses a quadratic form structure of the Hotelling

T2, thus, it is sensitive to shifts in the mean as well as to shifts in the variance

(Montgomery, 2005). The MEWMA test statistic becomes:

Zi=λxi+ (1−λ)Zi−1,(2.20)

where 0≤λ≤1and Z0=0. On the other hand, the quantity plotted on the

MEWMA control chart is

i=Z0

iΣ−1

ZiZi,(2.21)

where the covariance matrix is given by

ΣZi=λ

2−λ1−(1 −λ)2iΣ.(2.22)

Notice Equation (2.22) is analogous to the variance of the univariate EWMA.

The shift in the mean vector is described by the noncentrality parameter δ:

δ=µ0Σ−1µ1/2.(2.23)

Except for very large values of δ(or large shifts in the mean), ARLs generally

tend to increase as λincreases (Montgomery, 2005). Smaller values of λare,

thus, preferred to detect small shifts in the mean. Stoumbos and Sullivan (2002)

and Testik et al. (2003) report that small values of λalso result in MEWMAs

that are very insensitive to the form of the underlying multivariate distribution

of the process data, making MEWMAs robust to the assumption of normality

(Montgomery, 2005).

Since the MEWMA with λ= 1 is equivalent to the T2control chart, the

MEWMA is more sensitive to smaller shifts (Montgomery, 2005). Also even

when the Hotelling’s T2statistic is optimal for detecting the hypothesis that

the mean vector of a multivariate normal distribution is equal to a constant,

the MEWMA can be designed to have faster detection capability (Montgomery,

2005). However, while the use of the MEWMA chart will be an improvement

over the T2chart in regards to its robustness to the normality assumption, as

the dimensionality of the problem increases, the shift in the mean vector is still

“diluted” in the high dimensional space of the covariates.

2.3. Tango’s Method

In addition to monitoring for out of control scenarios as in SPC, disease

counts can be spatially clustered through distance methods. A well known exam-

ple by Tango (1995) uses a class of tests for detecting “focused” and “general”

clusters on aggregated counts. “General” tests are concerned with the over-

all pattern of a disease over a study region whereas “focused” tests assess the

clustering of observed cases around pre-specified points (or foci) (Besag and

Newell, 1991). Tango (1995) performs the “general” test using an extension of

the method proposed by Tango (1984) that further allows to incorporate popu-

lation densities as well as a confounding factor. This method is performed in

two stages. First, the “general” test considers the null hypothesis of no disease

clusters. Then, if this hypothesis is rejected, “focused” tests are applied to ap-

propriately select candidates for disease clusters foci to estimate whether one or

more disease clusters exist.

Suppose a study area is divided into mgeographical regions. Let random

variable Nidenote the number of cases in region iwith observed value nifor i=

1, . . . , m. Under the null hypothesis of no clustering (H0), Niare independent

Poisson variables where the mean of Niis proportional to the population size ξ

in region i:

H0:E(Ni) = λξi, i = 1, . . . , m. (2.24)

Under H0, the total number N=N1+···+Nmof cases is a sufficient statistic

for the unknown parameter λ. A test statistic whose distribution is independent

of unknown λcan be obtained by conditioning on the total number of cases

N=n. Given n, the disease locations are the values of a random sample of

size nfrom a multinomial distribution with parameter p0= (p1,...,pm)where

pi=ξi/Pjξjfor i= 1, . . . , m (Tango, 1995).

On the other hand, a reasonable alternative hypothesis is to detect the cluster-

ing of a rare disease around a pre-specified focus i0:

H1:E(Ni) = λξi(1 + ai0iε), i = 1, . . . , m, (2.25)

where ε > 0and ai0iindicates the measure of closeness between the ith region

and the pre-specified focus. The measure of closeness aij is a monotonically

non-increasing function of distance dij and aii = 1. When the mode of the

disease cluster is ignored, Tango (1995) proposes using an exponential form of

aij:

aij = exp (−dij /τ).(2.26)

Parameter τin Equation (2.26) is a scale parameter. Large values of τwill lead

to a test sensitive to large clusters, while small values of τwill lead to a test

sensitive to small clusters. Nevertheless, the choice of τwill not drastically

change the test results.

The test statistic for the “focused” test is given by:

CF=w0A(r−p)(2.27)

with variance

var (CF) = w0AVpAw,(2.28)

where Vp=∆(p)−pp0,∆(p)is the m×mdiagonal matrix based on vector

p, and for a given set of foci Ω={i0, i1, . . . , iI},

wj=









1if j∈Ω

0otherwise.

The test statistic for the general test allows for non-uniform population size over

study area:

CG= (r−p)0A(r−p)(2.29)

where the elements of vector rare determined by ri=ni/n for i= 1, . . . , m.

The statistic CG, as opposed to CF, has a substantial amount of positive skew-

ness. As a result, a normal approximation is not valid for the standardized statis-

tic (Tango, 1995):

TG=CG−E(CG)

pvar (CG)(2.30)

with

E(CG) = 1

ntr(AVp)(2.31)

and

var (CG) = 2

n2tr(AVp)2.(2.32)

Nevertheless, TGcan be approximated with a chi-square distribution with de-

grees of freedom adjusted by skewness

ν=8

qβ1(CG)2.(2.33)

The significance of the test can be calculated using

Pr {CG> c}= 1 −I ν+TG√2ν

2,ν

2!,(2.34)

where

I(x, φ) = Zx

Γ(φ)e−ttφ−1dt. (2.35)

If the null hypothesis of no clustering is rejected, it is important to figure out

location of the disease cluster(s). This is done as follows (Tango, 1995):

1) Select the regions satisfying the criterion {ri−pi> θ}as candidates for po-

tential foci where θis a pre-specified threshold.

2) When considering only one disease cluster, do a focused test CFwhere each

of the selected regions is assumed to be a focus. When there might be more

than one disease cluster, apply the focused test CFusing several weight vec-

tors w.

3) Examine the location(s) of the significant cluster(s) based on the p-values

computed from CF.

This method has been widely applied to monitor for disease clusters. How-

ever, some of the difficulties that arise when implementing this approach in-

clude: (a) geographical regions are pre-defined, (b) the method does not handle

case-event data, (c) the method only allows for the inclusion of one confounding

factor, (d) the null hypothesis in the “general” test is tested around pre-specified

foci, and (e) regions in the “focused” test are defined by means of a pre-specified

threshold.

2.4. Scan Statistics

Similar to control charts, scan statistics were not originally developed as a

public health surveillance approach. The original scan statistic was developed

by Naus (1965) as a theoretical framework to obtain clusters of points in two di-

mensions using a rectangular scanning window of fixed size. Scan statistics were

first applied to public health surveillance nearly two decades later. In fact, Kull-

dorff (1997) made the method very popular when he introduced a spatial scan

statistic that uses circular windows of various sizes to monitor for disease clus-

ters in multiple geographical regions. Many variations of Kulldorff’s scan statis-

tic have been used in temporal (Naus and Wallenstein, 2006) and spatiotemporal

(Kulldorff, 2001) surveillance. Moreover, a multivariate extension of Kulldorff’s

scan statistic (Kulldorff et al., 2007) allows to simultaneously analyze multiple

data sets.

In general, scan statistics attempt to detect clusters in a multi-dimensional

point process which are not explained by the baseline process. The baseline

process may be any inhomogeneous Poisson or Bernoulli process with intensity

proportional to some known function (Kulldorff, 1997). In public health surveil-

lance, scan statistics can be used on both case-event data and aggregated counts.

They allow to detect a local excess of disease incidences (Lawson and Kleinman,

2005) by calculating the maximum of some function and comparing the observed

and expected numbers of events in a spatial, temporal, and/or spatiotemporal

window (Naus and Wallenstein, 2006). Public health officials test whether this

local excess of events may have reasonably occurred by chance or whether there

is evidence of a significant increase in the disease incidence rate being moni-

tored. The scanning area is allowed to vary in such a way that a moving window

can test for multiple clusters of various sizes. At different positions, the scanning

window includes different sets of neighboring areas (Kulldorff, 2001). The best

choice for the shape of the scanning window is application-dependent. Com-

mon choices for the shape of the window include: all circular subsets, all circles

centered at any of several foci on a fixed grid with a possible upper limit on

the circle size (Kulldorff and Nagarwalla, 1995), all circles centered at any of

several foci on a fixed grid with a fixed circle size (Turnbull et al., 1990), all

rectangles of a fixed size and shape (Naus, 1965), and even cylinders when si-

multaneously scanning space and time (Kulldorff, 2001). On the other hand,

although the size and the shape of the scan window can be varied, pre-defining a

shape for the scanning window imposes a restriction on the type of clusters that

can be detected.

Conceptually, The scan statistic could be applied to one (or several) other co-

variates rather than spatial and temporal data (Kulldorff and Information Man-

agement Services, Inc., 2009). Other covariates need to be integrated through

expected values which requires additional models, such as logistic regression).

However, as the number of covariates increases, computational demands also

increase significantly due to the complexity of analyzing multiple length inter-

vals for each covariate. Consequently, only spatial and temporal covariates are

typically incorporated into the scan statistic, making it difficult to allow for any

higher-order interactions in these models. Unfortunately, without the simultane-

ous study of multiple covariates, disease clusters defined by higher-order inter-

actions can be missed.

2.4.1. Temporal Scan Statistics

Temporal scan statistics resemble the well-known CUSUM control chart.

They are used to detect changes retrospectively and prospectively in disease

incidence rates (Joner Jr et al., 2008). In retrospective surveillance, a past se-

quence of data is examined to determine whether there has been an increase in

the incidence rate sometime in the study period. In prospective surveillance,

new information is used as it becomes available to help detect any rate increase

as soon as possible after it occurs (Joner Jr et al., 2008). Prospective methods

are intended to identify an increase in a disease incidence rate before the dis-

ease spreads to a larger sector of the population. The literature on prospective

methods, however, is not as mature as the literature on retrospective methods.

Public health officials and epidemiologists often focus on a large temporal

disease cluster, and evaluate whether this cluster is likely or unlikely to have

arisen by chance given random variations from the background incidence (Naus

and Wallenstein, 2006). It is important, nevertheless, to take into account the

multiple comparisons inherent in the scanning of a larger review period. Naus

and Wallenstein (2006) define three different scenarios in the application of tem-

poral scan statistics. In the first scenario, events are distributed over a continuous

time frame, and a window of predetermined length scans for time intervals with

an increased incidence rate. The scan statistic, then, becomes the maximum

number of events in a window of length wthat scans time period [0, T ]. Wallen-

stein (1980) assumes tis the proposed duration of the epidemic under study pe-

riod Tand r=t/T . Wallenstein (1980) provides some tables for p=P(n, N, r)

given N≤100, where p=P(n, N, r)represents the probability that a maxi-

mum number of points in any interval of length tis greater than or equal to n,

under the hypothesis that Nincidences are uniformly distributed over interval

T. If p=P(n, N, r)is the probability that if Npoints are randomly distributed

along a line of length one, then, there is an interval of length rthat contains at

least n<Ncases (Weinstock, 1981). Wallenstein and Naus (1973), Huntington

and Naus (1975), and Neff (1978) provide general formulas for p=P(n, N, r).

In the second scenario, time is divided into Tdisjoint intervals, so that the

scan statistic becomes the maximum number of cases within wconsecutive in-

tervals. The scan statistic searches for a predefined time of length window of

consecutive event occurrences. Finally, the third scenario deals with prospective

scanning for disease clusters using time windows that always cover the most re-

cent available data. The data is viewed as a sequence of Ttrials. Each trial is

recorded whether or not an event has occurred, and the scan statistic becomes the

maximum number of events in any wconsecutive trials. For all these scenarios,

simple approximations for p-values exist for the constant background case that

allow to test for the statistical significance of clusters (Glaz et al., 2001; Glaz

and Balakrishnan, 1999; Balakrishnan and Koutras, 2002) .

2.4.1.1. Bernoulli Scan Statistic

Generally, health-related applications monitor the rate of positive diagnoses

as Bernoulli trials. Disease diagnoses are repeatedly performed, and at the end

of ith repetition, an outcome Yiis recorded. If the outcome is an incidence, then

Yi= 1; otherwise, Yi= 0. It is further assumed that these Bernoulli trials

are independent with disease incidence rate (p) equal to a constant p0(Joner Jr

et al., 2008). However, there are situations where pis not constant. In these

situations, risk-adjusted methods are used to detect a shift in the constant odds

ratio assumed for the in-control case scenario.

The Bernoulli scan statistic has often been used in retrospective surveillance.

This scan statistic is defined as the maximum number of incidences occurring

within any mconsecutive Bernoulli trials. All of the N−m+ 1 possible subsets

of mconsecutive trials taken from a larger set of N > m consecutive opportuni-

ties need to be considered (Joner Jr et al., 2008). In this context, the scan statistic

is retrospective, as it is assumed that all of the necessary Nobservations have al-

ready been collected before the scan statistic is computed (Joner Jr et al., 2008).

Nevertheless, it is possible to create a control chart based on the Bernoulli scan

statistic to prospectively detect an increase in the incidence rate from p0to p1,

where p1is the smallest increase in the incidence rate that needs to be detected

quickly. The control chart statistic at the ith trial takes the form

Si=

j=max(1,i−m+1)

Yj.(2.36)

A signal or alarm is given if Si≥k,2≤k≤m. Unfortunately, it is not easy

to determine the most suitable set of values for kand m. Increasing kwhile

holding mconstant, implies that incidences will need to occur more frequently

in order to produce an out-of-control signal. On the other hand, increasing m

while holding kconstant means that an out-of-control signal is more likely to

occur, regardless of whether an actual rate increase has occurred. To obtain a

compromise between these contrasting goals, Joner Jr et al. (2008) developed an

algorithm to help target a specified average number of trials until a false alarm.

2.4.1.2. Generalized Scan Statistic

Temporal scan statistics are generally not valid in situations in which there

are changes in (a) the population at risk, (b) the detection rate of the disease,

and (e) the levels of known risk factors (Weinstock, 1981). The generalized scan

statistic proposed by Weinstock (1981) eliminates the assumption of constant

population at risk by changing the constant length intervals to intervals based on

a fixed number of person-years at risk (u). In other words, if at time xthe

population at risk is defined by p(x). Then, the intervals would become all

continuous intervals that satisfy

Zxf

p(x)dx =u. (2.37)

For changing p(x)functions, the intervals defined by Equation (2.37) will be of

different durations. Under the null hypothesis, all intervals are expected to have

approximately the same number of cases. As a result, if nbecomes the maximum

number of cases expected in any interval, RT

0p(x)dx =Uand r=u/U, then,

the tables in Naus (1965) and Wallenstein (1980) can be used to test the null

hypothesis of no disease clusters.

In order to overcome the restriction of constant detection rate, Weinstock

(1981) proposed defining intervals based on a(x)where a(x)is the probability

of detecting at time xa person has a given disease. At this point, intervals that

account for non-constant population at risk as well as non-constant detection

rates can be defined by

Zxf

a(x)p(x)dx =v, (2.38)

where vis fixed and pre-determined. Once again, the tables in Naus (1965) and

Wallenstein (1980) can be used if V=RT

0a(x)p(x)dx,r=v/V , and nbe

the maximum number of cases in any interval. Notice that for all intervals, the

expected number of detected cases is Nv/V regardless of interval length. Also,

parameter vis typically chosen by letting v=rV where 1/r is the number of

non-overlapping intervals which could fit in the time frame [0, T ](Weinstock,

1981).

When further accounting for changing levels in known risk factors, Weinstock

(1981) defines intervals using

Zxf

a(x)b(x)p(x)dx =v, (2.39)

where b(x)represents the a priori risk of developing the disease associated with a

person a time xby virtue of the history of the risk factor of the population at risk

at time x. Once more, vis predetermined, V=RT

0a(x)b(x)p(x)dx,r=v/V ,n

is the maximum number of cases in any interval, and the tables in Naus (1965)

and Wallenstein (1980) can be used to test for disease clusters.

2.4.2. Spatial Scan Statistic

Spatial scan statistics were developed to test for geographical clusters and to

identify their approximate location. This section presents a brief summary of

the spatial scan statistic for the Bernoulli model (or case-event data) using the

notation by Kulldorff (1997). The spatial scan statistic for the Poisson model (or

aggregated counts) is fully described in Kulldorff (1997).

One can define a process where individual cases define points on an two-

dimensional geographical area Awhere there is a window of size wthat is

moved along A. Typically, a circular window Zis used, and the size of the

window (w) becomes the circle’s radius. The center of circle is moved, so that

at different positions of the window, the window includes different neighboring

geographical locations. The size of the window (w) is also varied from 0 up to

a maximum radius, in such a manner that the window never includes more than

50% of the population at risk (Kulldorff, 2001). Notice that this method creates a

very large number of distinct windows where each window represents a potential

cluster of events.

Over all possible values of the radius, the maximum number of points in the

window is recorded and compared to its distribution under the null hypothesis of

a Bernoulli-based random process. Consider the hypothesis

H0:p=q

H1:p>q

(2.40)

where pis the probability a specific case is a point within window Zand qis

probability the case lies outside window Z. The probability of being a point

within window Zis independent of all the others. Equation (2.40) shows that a

window Zis a potential cluster when the probability of being a case within the

window exceeds the probability of being a case outside the window. With a fixed

scanning window, the maximum number of points in the window at any given

time becomes the test statistic for Equation (2.40). When the scanning window

is no longer fixed, a likelihood ratio test statistic (λ) is used instead (Loader,

1991).

Conditioning on the total number of cases (C), the test statistic λbecomes

the maximum likelihood over all possible circles Z,

λ=maxZL(Z)

= max

L(Z)

,(2.41)

where L(Z)is the maximum likelihood for an individual circular window Z,

and L0is the likelihood function under the null hypothesis. Let cdenote the

observed number of cases in window Z,nthe total number of cases and controls

within window Z, and Nthe total number of cases and controls observed. Then,

the likelihood ratio for an individual circle becomes

L(Z) =c

ncn−c

nn−cC−c

N−nC−c

∗. . .

. . . ∗(N−n)−(C−c)

N−n(N−n)−(C−c)

I(c>E(c)),

(2.42)

where I(c>E(c)) is an indicator function. When scanning for clusters with

high rates, I(c > E(c)) = 1 whenever the number of observed cases within

window Zis larger than the expected number of cases under the null hypothesis

and is zero otherwise. When scanning for clusters with low rates, the indicator

function in Equation (2.42) becomes I(c<E(c)), so that I(c<E(c)) = 1

whenever the window has less cases than expected under the null-hypothesis

and zero otherwise. Similarly, when scanning for either low or high incidence

rates, the indicator function becomes one for all possible values of cases within

window Z.

To detect the most likely cluster, there is the need to find the window ˆ

Zthat

maximizes the likelihood function in Equation (2.42) conditioned on Z. Next,

compare the likelihood of ˆ

Zagainst the likelihood of the null hypothesis using

λ. The distribution of λ, unfortunately, cannot be determined in closed analytical

form. Instead, Monte Carlo replicates of the data set under the null hypothesis

conditioned on the total number of cases and controls (N) within the study area

are used. These replications allow one to assess the likelihood of the most likely

cluster(s). The p-value for λis defined by

P−V alue =R

nrep + 1,(2.43)

where Rrefers to the rank of L(ˆ

Z)when compared with the maximum likeli-

hoods from the Monte Carlo replicates. The number of Monte Carlo replicates

(nrep) is typically restricted to 999 or another number ending in 999 to simplify

the task of comparing p-values against the most common choices for signifi-

cance levels–0.05, 0.01, and 0.001 (Kulldorff and Information Management Ser-

vices, Inc., 2009). Turnbull et al. (1990) discusses in detail how to apply the

Monte Carlo sampling technique proposed by Dwass (1957) in the context of

scan statistics.

When cases are classified by categorical covariates, such as different age and

ethnic groups, the spatial scan statistic allows one to adjust expected counts

based on the different levels of these covariates. Frequently, when adjusting

the Poisson-based scan statistic for covariates, the expected counts are calcu-

lated through indirect standardization. Let cist be the number of cases and pist be

the population in covariate group i, geographical area s, and time period t. The

covariate adjusted expected count becomes PipistCi/Pi, where Ciis the total

number of cases in level iand Piis the total population in level i(Kulldorff et al.,

2007). Note, however, that this type of adjustment provides no differentiation in

the informational content between cases in different covariate groups. For in-

stance, an influenza-related death in a young adult provides more evidence of an

outbreak with serious public health consequences than an additional influenza-

related death in a 90-year-old (Kulldorff et al., 2007).

2.4.3. Spatiotemporal Scan Statistic

Unlike the two-dimensional circular window used in the spatial scan statis-

tic, the spatiotemporal scan statistic proposed by Kulldorff (2001) uses a three-

dimensional cylindrical window Z. The circular base of the cylinder represents

a geographical area, in the same manner as the spatial scan statistic, whereas

height represents time (Kulldorff, 2001). The cylinder is flexible in its circular

geographical base as well as in its starting date, independently of each other.

For each possible location and size, the spatiotemporal scan statistic considers

only those cylinders that reach all the way to the end of the study period. In this

regard, the spatiotemporal scan statistic considers only disease clusters that are

that are still active, while ignoring those that may have existed historically, but

which no longer a threat to the public.

Let [Y1, Y2]be the study period, and let sand tbe the start and end dates,

respectively, of cylinder Z. The spatiotemporal scan statistic considers all cylin-

ders for which Y1≤s≤Y2and t≥Ymwhere Ymis the start of the study period.

For each cylinder Z, a log likelihood ratio LLR(Z)is calculated, and the test

statistic is defined as the maximum LLR(Z)over all Z(Kulldorff et al., 2007).

The likelihood ratio test statistic is constructed in the same way as for the purely

spatial scan statistic. At any given moment, an observed cluster is statistically

significantly at the α-level if the probability of having detected a cluster with

higher likelihood during any of the previous analyses or the present analysis is at

most α. Notice, however, that the spatiotemporal scan statistic algorithm is more

computer intensive, since the likelihood for each window has to be calculated in

three, rather than two, dimensions.

2.4.4. Multivariate Scan Statistics

When looking prospectively for disease clusters, monitoring for symptoms

(or syndromic surveillance) as opposed to waiting for diagnoses is becoming the

standard. Kulldorff et al. (2007) proposes an extension of the spatiotemporal

scan statistic that allows one to prospectively monitor for disease clusters while

simultaneously incorporating multiple data sets into a single likelihood function.

With the multivariate scan statistic, a signal is generated whether a signal occurs

in only one or in multiple data sets. Some major reasons for taking a multivariate

approach to syndromic surveillance include (Kulldorff et al., 2007):

(a) No single data source captures all individuals in the disease cluster. Depend-

ing on the disease, some patients will go to a pharmacy, call a physician or

a nurse hot-line, visit their regular physician, go to the emergency room, or

call the ambulance.

(b) Some diseases will typically manifest themselves with a single symptom,

while other diseases can cause a wide variety of symptoms in different indi-

viduals.

entire population or mainly the elderly. If a work-site exposure is causing

the outbreak, mainly working-age adults who are affected.

Furthermore, one not always know which disease or even which combination

of data sources, symptoms and/or age groups to monitor. The multivariate scan

statistic addresses this problem as follows (Kulldorff et al., 2007):

1) For each cylinder Z, the log likelihood ratios are calculated for each data set,

and it is noted whether the observed number of cases is larger or smaller than

expected.

2) For each cylinder Z, the log likelihood ratios for data sets with more than the

expected number of cases are added. The second likelihood is the sum of all

the log likelihood ratios for data sets with fewer cases than expected.

3) The multivariate scan statistic is defined as the maximum of all the added log

likelihood ratios taken over all the cylinders.

The mathematical notation for the Poisson model described in Steps 1-3 is

given by

LRi(High, Z) = cz

nzczC−cz

C−nzC−cz

I(cz> nz)(2.44)

LRi(Low, Z) = cz

nzczC−cz

C−nzC−cz

I(cz< nz),(2.45)

where LRi(High, Z)is the likelihood ratio for high clusters and LRi(Low, Z)

is the likelihood ratio for low clusters under cylinder Zin data set i. Note that

for each cylinder, only one of Equations (2.44)-(2.45) will be non-zero. The test

statistic can now be expressed as

T= max

Z max X

LLRi(High, Z),X

LLRi(Low, Z)!!,(2.46)

where LLRi(High, Z) = ln (LLRi(High, Z)) and LLRi(Low, Z) =

ln (LLRi(Low, Z)). Note that the most likely cluster may either be a “high”

disease cluster with more cases than expected or a “low” disease cluster with

less cases than expected. If only “high” disease clusters are of interest, only

Equation (2.44) should be considered.

To adjust for different informational content in different covariate groups,

one can treat each covariate as a separate data set with a separately calculated

likelihood function. These likelihoods are, then, combined as follows (Kulldorff

et al., 2007):

1) Divide the data into different data sets, one for each covariate group.

2) For each cylinder Z, calculate the log likelihood ratio for each data set. Also,

note whether the number of cases is higher or lower than expected.

3) For each cylinder Z, the log likelihood ratios for the data sets with more than

expected number of cases are added. Similarly, the log likelihood ratios for

the data sets with fewer number of cases than expected are added.

4) For each cylinder Z, subtract the second sum from the first, to create the

log likelihood for “high” disease clusters with an excess number of cases.

Subtract the first sum from the second, to create the log likelihood for “low”

disease clusters with fewer cases than expected. The larger of these is the log

likelihood for this particular cylinder.

5) The maximum of all the log likelihoods, taken over all the cylinders, is the

definition of the informational content-adjusted scan statistic.

When only clusters with high rates are of interest, the same procedure is used,

but the second subtraction in Step 4 is skipped. The mathematical notation for

the Poisson model is given by

T=X

LLRi(High, Z)−X

LLRi(Low, Z),(2.47)

where LLRi(High, Z)and LLRi(Low, Z)are defined as above. When moni-

toring for both high and low disease clusters,

T= max X

LLRi(High, Z)−X

LLRi(Low, Z), . . .

. . . X

LLRi(Low, Z)−X

LLRi(High, Z).

(2.48)

Note that different types of adjustment may be used simultaneously for different

covariates. Nevertheless, as is the case with other scan statistics, covariates can-

not be directly incorporated into the analysis, and still need to be accounted for

through adjustments.

2.5. Generalized Linear Mixed Models (GLMM)

Generalized linear mixed models (GLMM) are also suitable for health-related

data. Kleinman et al. (2004) propose a method based on logistic regression with

a random intercept to estimate the probability that a subject is a case in a given

area on a given day. When new data becomes available, observed cases are ag-

gregated and the probability of a subject being a case is used to estimate the

likelihood of such event. If it is highly unlikely that such a large number of

cases on any given day will be observed, then, one might suspect there is an

increase in the specific disease rate being monitored. This method can be used

on aggregated counts as well as on case-event data. Case-event data, nonethe-

less, makes the model fitting process highly computationally intensive. GLMMs

can manage spatial, temporal, and even patient-specific covariates in numerical

and/or categorical format. The model takes the form:

E(yijt |bi) = pijt and logit(pit ) = xijt β+bi.(2.49)

With case-event data, yijt indicates whether person jin area iis a case on

day t,pijt is the probability that the person is a case, and xijt is a set of covari-

ates measured for person jand/or area iand/or day t. On the other hand, β

is a vector of fixed effects, and biis random effect for area i. Covariates xijt

can include: measures on individuals, measures on geographical areas, day de-

scriptors, seasonal indicators, among others. Random effects (bi) are typically

assumed to have a multivariate distribution with mean zero. The estimated bi’s

are sometimes called shrinkage estimators (Kleinman et al., 2004). These shrink-

age estimators have smaller mean squared error than simple logistic regression

estimators at the expense of non-zero bias.

With case-event data, as the number of covariates on individuals increases,

so does the number of scenarios to monitor. With as little as seven three-level

categorical covariates, over 2000 scenarios must be monitored. Moreover, when

a cluster is spread out across all levels of a covariate, the signal can be eas-

ily missed when monitoring each level individually. For example, if there is

an increase in asthma that affects all age groups in a given region, monitoring

teenagers alone can fail to detect the cluster. Once again, there is a dimensional

explosion that tends to limit the sensitivity of the method.

2.6. Wavelets

Shmueli (2005) proposes detecting disease clusters of unknown form using

wavelet-based methods. The main idea is to develop a more flexible set of mon-

itoring tools that does not suffer from the restrictive assumptions of normality,

constant variance, stationary, and non-correlated data. Wavelets can operate on

multiple scales and, unlike univariate control charts, they are proficient at detect-

ing multiple types of abnormalities. Frequently used in signal processing and

compression, wavelets use a complete orthonormal basis that shrinks and selects

coefficients to create a sparse representation of the data (Hastie et al., 2001).

Wavelets allow the representation of smooth and/or locally bumpy functions,

such as time series, in terms of coefficients that are associated with a particular

time and frequency (Hastie et al., 2001; Percival and Walden, 2000). As op-

posed to Fourier transforms that maintain specificity only in frequency, wavelet

decompositions allow for a localized frequency analysis (Shmueli, 2005).

2.6.1. Discrete Wavelet Transform

The following description will be based on the notation by Percival and

Walden (2000). Let w=WXbe the discrete wavelet transform (DWT) of a

real-valued time series Xof length nwhere wis an n×1vector of coeffi-

cients and Wis an n×northonormal matrix. Then, the multiresolution analysis

(MRA) of Xis given by

Xt=

j=1

Dj,t +AJ,t,(2.50)

where Dj=W0

jwjis the detail associated with Xat scale j, and AJ,tis the

sample mean.

The first scale is obtained by passing Xthrough a high-pass filter. This yields

the first-scale detail coefficient cD1. Then, Xis passed through a low-pass filter,

which yields the first-scale approximation coefficient vector cA1. The same two

operations are performed on cA1to yield cA2and cD2. The process can be

repeated for as many as J= log2(n). The standard DWT, however, includes

a down-sampling step after each filtering, such that the filtered output is half

the input. Thus, a small number of scales is used in practice. For public health

surveillance, a maximum scale of 5 is common since it is reasonable to compare

daily counts with neighboring days up to a month (or the 32-day window in level

5) (Shmueli, 2005).

The vector of wavelet coefficients wcan be organized into J+ 1 vectors

w= [cD1,cD2,...,cDJ,cAJ]where cDjand cAjrepresent the detail and

approximation coefficients, respectively, at scale j. To reconstruct Xusing w,

apply upsampling and reconstruction filters starting on cDJand cAJto obtain

cDJ−1and cAJ−1. The procedure is repeated until Xis finally reconstructed

from cD1and cA1.

2.6.2. Wavelet Monitoring

The underlying idea of using wavelets to monitor for changes in disease inci-

dence rates is to decompose health-related data using a DWT and, then, to eval-

uate the individual detail and approximation levels. Depending on the method,

one can monitor wavelet coefficients (Aradhye et al., 2003) or reconstructed ap-

proximation and details (Renaud et al., 2005).

When retrospectively looking for changes in disease incidence rates, DWT

not only compares a case to its past, but also to its future. Wang (1995) suggests

that if there are patterns in a time series, they will manifest themselves as large

coefficients in high detail levels. In the absense of disease clusters (or an increase

in disease incidence), detail coefficients should be zero. Wang (1995) proposes

to use only wavelet coefficients that exceed Donoho’s universal threshold. Using

Donoho’s universal threshold, a coefficient will be zeroed out unless it exceeds

σ√2 log (n)/n, where σis estimated by the median absolute deviation of the

coefficients at the finest detail divided by 0.6745 (Wang, 1995). An alternate

method called Multiscale Statistical Process Control (MSSPC) is proposed by

Bakshi (1998). The method consists of the following steps (Shmueli, 2005):

1) Perform a DWT on data, known to be free of disease clusters, and use this

decomposition to estimate standard deviations and to compute control limits

from scale-specific coefficients.

2) Perform a DWT on test data.

3) Threshold the coefficients using the control charts from Step 1, and zero out

coefficients within the control limits.

4) Reconstruct the series using the thresholded coefficients. Control limits at

time tare based on estimated variances from scales where coefficients ex-

ceeded their limits at time t.

5) Trigger an alarm if the reconstructed series exceeds the control limits.

This method is computationally inexpensive and easy to interpret. It relies on

the decorrelation property of DWT, where coefficients within and across scales

are approximately uncorrelated. Aradhye et al. (2003) shows MSSPC outper-

forms traditional control charts in detecting a wide range of anomalies with and

without autocorrelation.

A disadvantage of MSSPC is that it introduces a multiple testing problem

when monitoring each level in the decomposition at every point in time. An

alternative to this problem is to skip the reconstruction step and account for mul-

tiple testing by integrating a False Discovery Rate (FDR) correction (Benjamini

and Hochberg, 1995; Storey, 2002) into Step 3 of MSSPC. This FDR correction

can be used to control the expected proportion of falsely rejected hypotheses

among all detections. Suppose there are mcoefficients at scale j, such that

P(1) ≤P(2) ≤. . . ≤P(m)are the ordered p-values for which H(i)is the null

hypothesis for ith coefficient corresponding to P(i)(Benjamini and Hochberg,

1995). For each scale j, let kis the largest ifor which

P(i)≤i

mq∗(2.51)

where q∗is the desired ratio of falsely rejected hypotheses over all detections

(Benjamini and Hochberg, 1995). Now, reject all H(i)for i= 1,2, . . . , k. In

other words, FDR thresholds coefficients based on the collection of the p-values

at each scale (Shmueli, 2005), such that coefficients with p-values larger than

P(k)are zeroed out.

As with many other public surveillance methods, the problem of prospec-

tively looking for disease clusters is more challenging. The downsampling in-

troduced by DWT introduces a time lag for calculating coefficients that becomes

more acute at coarser scales (Shmueli, 2005). Also, DWT is highly dependent

on the starting point of the time series. For these reasons, an alternative to DWT,

known as redundant or stationary DWT (SWT), that skips the downsampling

step and is independent of the starting point is commonly preferred. The trade-

off, however, is that the wavelet coefficients are no longer uncorrelated within

scales. Nevertheless, scale-level correlation is easier to handle than the overall

correlation in the original series.

Unfortunately, both DWT and SWT are dependent on the boundaries of the

data set. The most typical approach to reduce the impact of this problem is to use

wavelet families with short filters, such as the Haar filter of length two. Renaud

et al. (2005) and Xiao et al. (2005) propose doing prospective monitoring using

a “backward implementation” of SWT that applies a non-symmetric Haar filter

for forecasting a time series using SWT with the Haar filter. This method is

equivalent to ordinary SWT, except for the fact that it shifts in time coefficients

at scale jby 2jpoints. Coefficients in different scales, thus, start at different

points in time. However, this is not a major concern in prospective monitoring,

for our interest focuses on the most recent data. Some appealing features of this

method include (Shmueli, 2005): (a) New data points have no effect on past

coefficients. (b) Its approximation and detail coefficients are much smoother

than the ones obtained with the downsampled Haar-based DWT. (c) It does not

require the time series to be of special length.

Wavelet decompositions allow for a smoother representation of health-related

data. Nonetheless, as with many other multivariate methods, there are inter-

pretability issues. Even though DWT allows to define the specific point in time

where there is a change in the disease incidence rate, it is not trivial to identify

which covariate or set of covariates is responsible for the change.

CHAPTER 3

SURVEILLANCE AS SUPERVISED LEARNING

3.1. Introduction

Public health surveillance collects, analyzes, and interprets health-related

data in order to understand trends and detect changes in disease incidence rates

(Thacker, 2000). Frequently, interest is strictly focused on detecting increases in

disease rates. This problem is typically addressed by simply looking for disease

clusters within temporal and/or spatial regions. In order to do so, incidences

from different time periods are contrasted. The data for the first time period

often represents incidences that behave according to the base incidence rate of

the disease, and thus, is typically referred to as baseline data. Baseline data is

a point of reference, so it is usually “free from disease clusters”. Notice “free

from disease clusters” means the data behaves according to what is expected to

be the natural realization of the process for the disease of interest. Data for the

second time period is often associated with new test disease incidences that are

contrasted against the baseline data, in order to evaluate whether there have been

any changes in the disease incidence rate. Increases in incidence rates (or occur-

rences of disease clusters) are detected when there is at least one region with an

unusually large count of new test incidences.

This chapter proposes the Feature Selection (FS) method for detecting

changes in disease incidence rates through a transform of public health surveil-

lance to supervised learning. This approach can greatly increase the flexibil-

ity and robustness of numerous public health monitoring and informatics tasks.

The method uses computationally-intensive, but feasible, non-linear learners to

address public health surveillance with a large number of covariates. In ad-

dition, it uses case-event data from modern public health records, where indi-

vidual disease incidences are recorded along with a large number of additional

patient-level covariates. A change in any region (or multiple regions) of this

high-dimensional space is of interest, and the regions are not pre-specified. Fur-

thermore, a region can be defined from cross-products of all these patient-level

covariates which can be categorical, numerical, or even of mixed type. In this

work, the baseline represents historical incidences of the disease of interest; al-

though in some cases, the baseline can be generated from a different type of

event.

Section 3.2 defines the public health surveillance transform to supervised

learning. Section 3.2 also discusses the important topic of signal criteria. Sec-

tion 3.3 shows examples of the Feature Selection (FS) method on a variety of

simulated data sets and provides some comments on the spatial scan statistic.

Section 3.4 shows the power curves for the proposed method. Finally, Section

3.5 presents concluding remarks.

3.2. Transform to Supervised Learning

Supervised learning models an output or target (or multiple outputs) from in-

put measurements (e.g., regression and classification models). Supervised learn-

ing provides a defined target for the model (the output), so that models are easier

to score, evaluate, compare, and refine as needed. There is typically less subjec-

tivity in supervised learning models than in other traditional public health meth-

ods (such as density estimates). Based on an analogy with multivariate SPC,

one might start with count data and transform to supervised learning in a manner

similar to Hwang et al. (2007). However, counts would suffer from the need to

pre-specify multiple regions among spatial and other covariates (as previously

discussed in Chapter 1).

Using case-event data, public health surveillance can be transformed into a

supervised learning problem in a simple, but effective, manner. For each event

in a baseline (historical data), the class label becomes zero while for each event

in the current data the class label becomes one. A supervised learning model is,

then, used to predict the class of each case based on covariates. The basic idea

is that the model does not predict well if the new and baseline events have sim-

ilar covariate information. That is, if a model cannot separate the 0 - 1 classes;

then, no regions in covariate space are detected to have different rates. Con-

versely, a supervised learner that can detect different rates in any region of the

covariate space can be used to signal a change. Note that it is not required to

pre-specify the form of a region. One or multiple regions can be detected with

a suitable learner. Another advantage of the transform to supervised learning is

that a learner produces a number of quantitative results that can be used to detect

changes. Previous research in manufacturing shows how supervised learners are

used to identify unusual results (Hwang et al., 2007). In surveillance, a similar

learner can be used to detect differential rates between the classes in localized

regions of covariate space.

For example, refer to Figure 3.1. This figure illustrates a region of greater

intensity of disease incidences for a particular age group and geographical loca-

tion. This represents an interaction between the geographical location and age

covariates in the supervised learning model. Consequently, a suitable learner

needs to be sensitive to interactions between covariates. Also, this example illus-

trates three covariates, but higher order interactions are of interest. For example,

a cluster might comprise age, location, and poor air quality (PM10 measure-

ments). Moreover, covariates in public health data can be categorical (such as

gender, emergency department type, etc.), so that the supervised learner should

be able to integrate this type of information as well. A number of supervised

learners may be applied to this task. However, the ability to detect higher-order

interactions among the covariates is the key to identify clusters.

Fig. 3.1. Single disease cluster in 3-D point process.

The simplest example of a data table for the surveillance transform is shown

in Table 3.1. Notice that Table 3.1 is also an example of a case-event data set

where column vectors Age through Gender represent patient-level covariates.

Column vector Date stands for the date the incidence was observed. Column

vectors Lat and Long represent the latitude and longitude coordinates for the

patient’s address. In this example, it is assumed that there are no observed dis-

ease outbreaks during the month of October, and that disease incidences for the

Table 3.1. Example of supervised learning transform.

Event Date Age Lat Long Gender Age* Lat* Long* Gender* Class

1 11/5/2009 9 18 26 F 58 19 30 F 1

2 11/3/2009 16 16 22 M 16 14 26 F 1

3 11/2/2009 47 20 18 F 34 12 22 F 1

4 10/30/2009 63 14 22 F 47 20 18 M 0

5 10/28/2009 9 13 27 M 63 14 27 F 0

6 10/28/2009 11 12 27 F 9 18 22 M 0

7 10/26/2009 58 19 30 M 9 16 27 F 0

8 10/15/2009 34 14 27 F 11 13 27 M 0

month of November are currently being monitored. Notice no observed disease

outbreaks for the month of October imply that the incidences for the month of

October behave according to the base incidence rate of the disease or, in other

words, depict the in-control behavior of the disease. The historical records from

October become training records for Class 0, while the test records from Novem-

ber become training records for Class 1. A classifier is, then, used to discriminate

between the two classes in order to look for a deviation in the disease incidence

rate. In brief, transforming public health surveillance to supervised learning is as

simple as identifying a data set of historical records “free from disease clusters”,

another data set that consists of new records that can be used to assess the cur-

rent performance of the process, and assigning each set of records to a different

class label. In this manner, the incidences or records will be in a format that is

suitable for any classifier. The reasoning behind columns Age* -Gender* will

be discussed in Section 3.2.2.

3.2.1. Signal Criteria

A transform to supervised learning allows incidences (or cases) to be assigned

to classes, but an important issue is the criteria to generate a signal from a rate

change. Several approaches can be considered. In classification problems, error

rates are commonly used to evaluate performance. However, some disadvantages

of a class error rate are discussed as a signal and, instead, this work focuses on

feature selection. Of course, other choices are possible and more work is needed

to research alternatives. It is also important to consider diagnostics to identify the

region or regions that generate a signal. An approach to this problem is presented

here, but a more complete discussion is beyond the scope of this work.

Error rates are natural scores from a learner. If error rates are high, the learner

cannot distinguish between the new and baseline data. This indicates that no

changes are detected. Conversely, if the new data can be classified correctly,

then, covariate regions exist where the new data differs from the baseline. In

general, one expects baseline events to occur in many regions of covariate space.

It is, therefore, not surprising that some baseline events (Class 0) will be pre-

dicted to be Class 1 events. Of greater interest are the Class 1 events that are

predicted correctly. Consequently, a low error rate for Class 1 events suggests a

change in a specific region.

If cases in the new data fall primarily in a few regions, one expects a decrease

in Class 1 errors, so this error rate might be used to define a signal. However,

even with an increased rate in a localized region, one might also expect a number

of cases in the new data that are similar to the baseline. These more random cases

are not expected to be classified correctly by the learner, and this would attenuate

the signal for the Class 1 error rate for regions as the ratio of random cases to

cases in a localized region increases.

An advantage of the transform to supervised learning is that additional tools

are available for many learners. Specifically, feature selection of the most im-

portant covariates to a classification model has been widely studied (Guyon and

Elisseeff, 2003). The role of feature selection for surveillance is to identify co-

variates that distinguish between the baseline and new data and, thereby, con-

tribute to a region with rate changes. If no covariates are identified, no signal is

generated. On the contrary, selected covariates are those that define the region

with the rate change. Although a number of feature selection methods may be

used, it is important that the method is sensitive to interactions among the co-

variates that define the localized regions of interest for surveillance. The number

of cases in a region may also be small relative to the number of cases in the new

data.

This work adopts the feature selection method based on tree ensembles pro-

posed by Tuv et al. (2009). A decision tree is a supervised learner that succes-

sively partitions (or splits) rows of data, so that the class is as consistent (pure)

as possible within each partition. Trees have been widely used as learners and

were described in detail by Breiman et al. (1984). Nevertheless, a single tree is

generated by a greedy algorithm that can lead to an unstable model in which a

small change to the data can change the model substantially. Ensemble meth-

ods construct multiple models (called base learners) and, in a parallel ensemble,

the weighted outcome (or vote) is used to classify new data. An exemplar for

tree-based ensembles is a random forest model (Breiman, 2001).

3.2.2. Variable Importance

Tree ensembles are adopted based on their inherent ability to isolate regions

of different proportions of events. In addition, public health data is particularly

suited to tree-based ensembles because they can deal with the complexities asso-

ciated with such data. They can incorporate categorical, numerical, and missing

data; potentially large interactive effects, and even a large number of variables

and cases. Other advantageous features of tree-based ensembles that appeal to

public health surveillance include: ability to fit nonlinear models that can adapt

to the characteristics of the baseline environment, invariance to variable units,

insensitivity to outliers, overfitting avoidance, and solution accuracy assessment.

More importantly, tree-based ensembles provide embedded importance mea-

sures that can identify the feature relevant to a model. For public health surveil-

lance, this can simplify an analysis to a more detailed study of a few covariates

from a potentially extensive set. For a single decision tree, the measure of vari-

able importance was proposed by Breiman et al. (1984):

V I(xi, T ) = X

t∈T

∆I(xi, t),(3.1)

where ∆I(xi, t) = I(t)−pLI(tL)−pRI(tR)is the decrease in impurity due to

an actual (or potential) split on covariate xiat a node tof the optimally pruned

tree T. For classification I(t) = Gini(t), where Gini(t)is the Gini index of

node t:

Gini(t) = X

i6=j

ipt

j,(3.2)

where pt

iis the proportion of observations in twhose response label equals i

(y=i) and i, j run through all response class numbers. The Gini index measures

node impurity. It is zero when thas observations only from one class, and it is

maximized when classes are perfectly mixed.

An important question remains for tree based models: how to rank covariates

that were masked by others with slightly higher splitting scores, but could pro-

vide as accurate a model if used instead? This can be addressed by a key variable

of CART (Breiman et al., 1984) called surrogate splits. In CART methodology,

the surrogate splits are used to handle missing values, detect masking, and assess

variable-importance ranking.

The predictive association of a surrogate variable xsfor the best splitter x∗at

a tree node tis defined through the probability that xspredicts x∗action correctly

which could be estimated as:

p(xs, x∗) = pL(xs, x∗) + pR(xs, x∗)

where pL(xs, x∗)and pR(xs, x∗)define the estimated probabilities that both xs

and x∗send a case in tleft (right).

The predictive measure of association λ(x∗|xs)between xsand x∗is defined

λ(x∗|xs) = min(pL, pR)−(1 −p(xs, x∗))

min(pL, pR)(3.3)

where pL, pRare proportions of cases sent to the left (right) by x∗. It measures

the relative reduction in error due to using xsto predict x∗as compared with the

“naive” max(pL, pR)rule. If λ(x∗|xs)<0, then, xsis disregarded as a surrogate

for x∗.

The sum in Equation (3.1) is taken over all internal tree nodes where xiwas

a primary splitter or a surrogate variable (for any primary splitter CART keeps

a number of surrogate variables). Often, a covariate that does not appear in a

tree will be ranked high on the variable importance list constructed based on

surrogate scores.

For ensembles of Mtrees, this importance measure is easily generalized. It

is simply averaged over the trees

M(xi) = 1

m=1

V I(xi, Tm).(3.4)

The effect of averaging makes this measure more reliable, and because of the

stochastic nature of the ensembles, the masking issue is not a problem (especially

for the independently built random forest). Therefore, in Equation (3.4), the

sum is evaluated over internal nodes where the variable of interest is the primary

splitter. In this case, Equation (3.4) generalizes to

M(xi, k) = 1

m=1

V I(xi, Tkm)

and represents the relevance of xito separate class kfrom the rest of them. The

overall relevance of xican be obtained by averaging over all classes:

M(xi) = 1

k=1

V I(xi, k).

With only variable importance scores, there is no obvious threshold to se-

lect only important covariates from a potentially large subset of covariates. For

surveillance, this issue is important to control risks associated with false alarms

and missed signals. The feature selection method used by Tuv et al. (2009) in-

corporates a statistical criterion that is directly useful for surveillance and briefly

described. The approach generates artificial variables known to not affect the

target. That is, artificial variables are constructed from random permutations of

each of the actual covariates in the data set. Random permutations break rela-

tionships between the class and each covariate. Permutations still maintain the

same marginal distribution of the covariate. Note that these artificial variables

can be constructed for categorical or numerical covariates in exactly the same

manner. In Table 3.1, column vectors Age∗through Gender∗represent the arti-

ficial variables for covariates Age through Gender, where artificial variable xv∗i

is simply a column vector of the permuted values of covariate xi. Because these

artificial variables are by definition not related to the target, one would expect

them not to be important.

An alternative way to create artificial variables is to permute the row vectors

as opposed to permuting the column vectors as previously described. In this

regards, the joint distribution of the original covariates is maintained in the ar-

tificial variables. Nevertheless, several experiments were used to evaluate both

permutation schemes and the former lead to remarkably better results. Please re-

fer to Appendix A for a description of the experiments considered and the results

obtained for both permutation schemes.

Importance scores can be calculated for the actual and artificial variables.

The artificial scores provide a distribution for covariate scores unrelated to the

response. If this procedure is replicated (with random permutations), the dif-

ferences between the variable importance for covariate xiat replicate j(V Iij )

and the pth percentile of the variable importance for the artificial variables at

replicate j(p1−α

j) can be used to create a test for statistically significant mean

difference. Only a variable with a significant difference (as measured with a

p-value) from the percentile of artificial variables is selected as an important fea-

ture. For surveillance, this implies a covariate is detected as contributor to a rate

change. The p-values also allow for simple multiple comparison adjustments

(such as Bonferroni’s adjustment). This mechanism simultaneously provides an

objective measure to detect changes and identifies the covariates that define the

clusters. Percentiles are typically chosen to be between 0.8 and 0.95. Empirical

tests by Tuv et al. (2009) suggested using the 80th or 85th percentile from the

artificial scores.

The surveillance work here uses these significance tests as described by Tuv

et al. (2009). In addition to artificial covariates, further enhancements include the

split weight estimated from out-of-bag (OOB) samples. The OOB samples give

better estimates of variable importance, particularly for mixed categorical and

numerical data, and compensate for categorical covariates with a large number

of levels. Besides, the effect of identified important covariates are iteratively re-

moved with a serial ensemble model to allow detection of less important covari-

ates. Furthermore, in linear models it is important to detect collinearity among

predictors. More generally, covariates that mask each other can both be identi-

fied as important covariates and masked sets can be indicated. Refer to Tuv et al.

(2009) for further details and some examples of applications of this method on

real data sets. Again, there are many advantages with the chosen feature selec-

tion algorithm for surveillance, but other feature selection algorithms can also

be applied.

The proposed algorithm is summarized as follows.

Feature Selection (FS) Method:

1) Create a permuted covariate for each covariate in the original data set.

2) Augment original data set with permuted covariates.

3) Perform supervised learning transform.

4) Select learner.

5) Calculate importance scores for each covariate.

6) Select covariates with significant importance scores.

7) If at least one covariate is selected, conclude there has been a change in the disease incidence

rate.

3.3. Examples

3.3.1. High-dimensional data

Several examples are used to illustrate the effectiveness of the supervised

learning transform. The first set of examples consists of a fully simulated data

set that allows one to illustrate the ability of the method to detect changes even

among high-dimensional data sets with many covariates. The simulated data

set includes 10,000 baseline cases randomly generated among one hundred co-

variates. These covariates are uniformly distributed between 0 and 1, but any

distribution, or even categorical covariates could be used. New data with single

or multiple clusters were generated among multiple covariates with a different

number of cases. Clusters within covariates were created by defining smaller

ranges in the specific covariates for the cases involved in the clusters. For exam-

ple, one experiment formed a cluster with twenty cases within a region defined

by two covariates (labeled x1and x2). A second experiment formed three clus-

ters with twenty cases each. Covariates x1−x5were used to define the three

clusters, but only covariates x1and x2were involved in all three clusters. Figure

3.2 shows that the greatest variable importance scores clearly correspond to the

important covariates in these examples. In Figure 3.7(a), only covariates x1and

x2stand out as important; whereas in Figure 3.7(b), covariates x1−x5stand out

as important. Notice, still, that covariates x1and x2in Figure 3.7(b) are ranked

remarkably higher than the rest of the covariates in the cluster.

Permuted artificial covariates can be added to generate statistical significance

tests. Even so, there was no need to further look at those permutations to select

important covariates in this set of experiments. Separation was highly distinct

by looking only at the standard importance measures.

(a) Correctly shows greatest impor-

tance scores for two variables that

define a cluster.

(b) Correctly shows greatest impor-

tance scores for five variables that

define three clusters.

Fig. 3.2. Variable Importance for clusters among 100 variables (only the top 10

scores shown).

3.3.2. Data with multiple locations

A problem that commonly arises in public health surveillance is how to si-

multaneously monitor for disease clusters in a large number of cities or even

counties. This next example consists of fifty categorical geographical locations

(e.g., 50 counties) with 3,600 baseline cases and 1,200 test cases randomly dis-

tributed. The data set also includes fifteen numerical covariates uniformly dis-

tributed between 0 and 1, and ten extra categorical covariates. Each of the addi-

tional categorical covariates is discrete uniform with five to thirty levels.

A disease cluster was superimposed on geographical location by forcing this

categorical covariate to level 1 (County 1) in 100 test cases. Ten different repli-

cates for this example were considered, and the results are described in Table

3.2. Covariate selection criteria for this example was based on a 0.05 p-value

and a 0.9 percentile importance score threshold. Clearly, the proposed method

accurately detected county in all replicates. However, there were additional co-

variates detected (false alarms) in data sets 1, 4, 6, 7, and 8. Note that the random

sampling strategy used can create scenarios where covariates that were not in-

tended to be part of the superimposed clusters still turned out to be part of them,

simply because of the nature of the random sample. Nevertheless, none of these

additional covariates was selected more than twice in the ten replicates. On a

similar note, the system time for these experiments was, on average, 44.02 sec-

onds with a standard deviation of 10.01 seconds on an Intel Core i5 processor.

3.3.3. Disease Clusters with Covariate Interaction

The following set of examples are simulated based on case-event data sets.

There are six different experiments with 2,500 baseline cases and 625 test cases.

Table 3.2. Feature selection results for experiments with a cluster of 100 cases

in County 1.

Replicate Covariates

Detected

1county,x23

2county

3county

4county,x2

5county

6county,x12

7county,x23,x10

8county,x5,x15

9county

10 county

As one wants to be able to detect an increase in a disease incidence rate as soon

as possible, the baseline data set (Class 0) will often be much larger than the test

data set (Class 1). In this manner, one evaluates Class 1 cases as soon as they

become available. This is also the main reason why the number of cases in Class

1 was simulated to be only 25% of the number of cases in Class 0.

The patient-level covariates for these examples include: latitude,longitude,

gender,age group, and race.Latitude and longitude stand for the geographical

coordinates associated with the patient’s address and are, thus, numerical. The

remaining covariates in the data set are categorical. The levels for age group con-

sist of ten-year intervals between “0 - 10” and “91+” years, while the levels for

race include: “White”,“Native American”,“Hispanic”,“Black”,“Asian/Pacific

Islander”, and “Other”. The six different experiments described in Table 3.3

were considered. Experiment 0 represents the null experiment where all data

Table 3.3. Description of experiments with covariate interaction.

Experiment Number of

Clusters

Cases in

Cluster

Covariates

in Cluster 1

Covariates

in Cluster 2

0 0 0

1 1 50 x,y

2 1 100 x,y

3 1 200 x,y

4 1 100 x,y,age group

5 2 100/100 x,y,age group x,y,race

was simulated to be random within the state of Arizona and state equivalent

territories polygon (U.S. Census Bureau Geography Division, 2001). For Exper-

iment 1, all data was simulated randomly except for the geographical coordinates

of 50 cases in Class 1, which were simulated to be clustered around a random

point within the polygon. All clusters were defined to be circular and of radius

= 0.4 units based on the latitude and longitude coordinates for the Arizona state

polygon. Experiments 2 - 3 are defined in the same manner as Experiment 1, but

with a different number of cases within the cluster. Experiment 4 represents a

dual cluster of 100 cases in age and space. Besides modifying the geographical

coordinates for 100 cases in Class 1, the age group for these cases was fixed to

“11 - 20”. The first cluster in Experiment 5 was defined as in Experiment 4.

The second cluster in Experiment 5 was defined by setting race = “Hispanic”

for 100 cases in the cluster. Each experiment was replicated ten times, and a

planar projection of the geographical coordinates {x,y}was used for all further

analyses.

Figure 3.3 shows an example of the intensity functions associated with the

different experiments. The smoothing parameter for all these functions was set to

s=p|cov(x, y)|. Figure 3.3(a) shows the null experiment. On the other hand,

Figure 3.3(b) shows Experiment 3 with a single cluster of 200 cases. Notice in

terms of geographical location, Experiment 4 in Figure 3.3(c) looks somewhat

similar to Case 3, as the interaction between age and space cannot be captured

in this plot. However, the circular shape of the cluster with 200 cases in Figure

3.3(b) is more clearly defined. The intensity function for Experiment 5 is shown

in Figure 3.3(d). Now, there are two clearly defined clusters of 100 cases. In this

case, the intensity function is plotted individually for each level of age group and

race (Figure 3.4), it is clear that age group = “11 - 20” and race = “Hispanic”

are also interacting with space. Although, not shown here, an intensity plot for

Experiment 4 shows a very similar space-age interaction to the one depicted in

Figure 3.4(a).

Feature selection on the experiments described in Table 3.3 was performed

using the algorithm by Tuv et al. (2009) with the following parameters: 100

trees, 100 replicates, and the 0.9 percentile for the artificial importance scores.

Results are summarized in Table 3.4. All entries in Table 3.4 represent the per-

centage of times a given covariate had a p-value ≤0.01 in ten independent repli-

Table 3.4. Average detection rates (over ten replicates) for six experiments with

interactions between different covariates. Columns age group -yrepresent

patient-level covariates where xand ystand for the planar projection of the

patient’s address geographical coordinates.

Experiment age group gender race x y

0 0 0 0 0.1 0

1 0 0 0 1 0.9

2 0 0 0 0.9 0.9

3 0 0 0 1 0.9

410011

510111

cates. One can further assume that a given covariate was assumed to be identified

as important in discriminating between baseline and new test cases if its detec-

tion percentage in Table 3.4 was greater than 0.5. For instance, in Experiment 0;

where there are no clusters, no covariate was identified as important in discrim-

inating between the two classes. In Experiments 1 - 3, which refer to strictly

spatial clusters, covariates xand ywere identified. Similarly, in Experiment 4;

which is the cluster based on the space-age interaction, correctly identified x,

y, and age group. Finally, in Experiment 5, all covariates involved in the two

clusters were identified: x,y,race, and age group. Consequently, the algorithm

was successfully able to identify all covariates used to define the clusters with

zero false alarms. It is important to acknowledge that the average system time

for these runs is 29.55 seconds with a standard deviation of 17.18 seconds on an

Intel Core i5 processor.

Table 3.5 summarizes the error rates for Experiments 0 - 5. Error rates were

Table 3.5. Mean error rates for experiments with covariate interaction. Error rates

were defined as the summation over the number of covariates that were incorrectly

detected divided by the total number of covariates in the data set.

Experiment Replicate Mean Error Rate

1 2 3 4 5 6 7 8 9 10

0 0 0 0 0.2 0 0 0 0 0 0 0.02

1 0 0 0 0 0.2 0 0 0 0 0 0.02

2 0 0 0 0 0 0 0.4 0 0 0 0.04

3 0 0 0 0 0 0.2 0 0 0 0 0.02

4 0 0 0 0 0 0 0 0 0 0 0

5 0 0 0 0 0 0 0 0 0 0 0

defined as the summation over the number of covariates that were incorrectly

detected divided by the total number of covariates in the data set. A covariate

is considered to be incorrectly detected whenever it is failed to be identified as

important by the feature selection algorithm (with a p-value >0.01) when it

was used to define a cluster (false negative). Also, detecting a covariate that was

not used to define a cluster (with a p-value ≤0.01) is considered an incorrect

detection (false positive). Since there are only five covariates in the data set,

error rates will increase in 0.2 increments. For each experiment, the error rate

will be reported as the mean error rate across all replicates. Notice that for all

experiments the error rate reported was no more than 2%, with the exception of

Experiment 2, where the method failed to identify both xand yin one of the ten

independent replicates.

3.3.4. Disease Clusters with Seasonality

Let λ(t)be a continuous function of time tused to define the incidence rate

of a seasonal disease. The presence of seasonal factors that trigger the incidence

of the disease will result in larger values of λ(t), while the absence of these

will result in smaller values of λ(t). Ideally, any public health surveillance sys-

tem should account for seasonality. That is, a public health surveillance system

should not generate false alarms in time periods where historically the incidence

rate for the disease being monitored is higher than for any other periods. For

instance, when monitoring asthma, one would not expect to trigger a signal just

because there are more asthma attacks observed in the spring as opposed to win-

ter. The surveillance system should be able to trigger an alert when the asthma

attacks for the spring period being monitored significantly exceed the baseline

asthma incidence rate for the spring.

This section will show how the proposed method can successfully account

for seasonality. Suppose incidences are recorded for two consecutive years (see

Figure 3.5(a)). Year 1 is assumed to be free of disease clusters, however, the

observed daily counts vary substantially throughout the year. The daily counts

appear to be stable during spring showing the highest daily counts for the year.

They also seem to be somewhat stable at lower daily count level during the

summer. For the remaining seasons, the daily counts are either increasing or

decreasing steadily. On the other hand, Year 2 represents test data. The data for

Year 2 in Figure 3.5(a) follows the same seasonal behavior as Year 1. The data

for Year 2 in Figure 3.5(b) also follows the same overall behavior of Year 1 data.

Nevertheless, there is an increase in daily counts during days 152-181. The FS

method can be used to capture such an increase in a disease incidence rate, while

ignoring all other seasonal fluctuations throughout the year.

Table 3.6 describes a series of simulated scenarios where the disease of inter-

est can assume three different incidence rates depending the time of the year. In

these experiments, Class 0 refers to case-event baseline data for a previous year

(say Year 1), while Class 1 refers to case-event test data for the current year in the

study period (say Year 2). Notice the baseline period for Class 0 does not neces-

sarily need to be a year. It can consist of data for several years, months, or even

weeks. There is a total of 3,600 cases in Class 0 where 720 cases are distributed

uniformly between days 1 and 120; 2,160 cases are uniformly distributed be-

tween days 120 and 240, and 720 cases are uniformly distributed between days

241 and 365. On the other hand, there are only 1,200 cases in Class 1. However,

they are distributed in a like manner. Twenty percent of the cases are uniformly

distributed between days 1 - 120, 60% of the cases are uniformly distributed be-

Table 3.6. Case distribution in Experiments 1 - 8 with seasonality.

tCases in

Class 0

Cases in

Class 1

1 - 120 720 240

121 - 240 2160 720

241 - 365 720 240

tween days 121 - 240, and the remaining 20% are uniformly distributed between

days 241 - 365.

Tables 3.7 - 3.8 refer to a series of mixed scenarios. Class 0 cases are defined

as in Experiments 1 - 8. In contrast, Class 1 cases can assume six different

incidence rates depending on the time of the year. Nevertheless, each Class 1

period with a fixed incidence rate in the Experiments 1 - 8 is split into two periods

of equal length, and the same number of cases for the original longer period are

now distributed unevenly between the two shorter periods. For example, Case 9

splits the 240 observations in the first 120 days as 140 observations uniformly

distributed between days 1 and 60, and 100 observations uniformly distributed

between days 61 and 120. Cases 10 - 13 follow the same pattern, but these cases

divide the observations more unevenly. For example, the 240 observations in

period 1 are split into 160 observations between days 1 - 60 and 80 observations

between days 61 - 120.

Refer to Table 3.9 for a description of the disease clusters superimposed in

the different experiments in Tables 3.6 - 3.8. Each data set consisted of 4,800

Table 3.7. Case distribution in Ex-

periment 9 with seasonality.

tCases in

Class 0

Cases in

Class 1

1 - 60 720 140

61 - 120 100

121 - 180 2160 330

181 - 240 390

241 - 300 720 80

301 - 365 160

Table 3.8. Case distribution in Ex-

periments 10 - 13 with seasonality.

tCases in

Class 0

Cases in

Class 1

1 - 60 720 160

61 - 120 80

121 - 180 2160 420

181 - 240 300

241 - 300 720 170

301 - 365 70

cases, twenty five covariates uniformly distributed between 0 and 1, and a time

(t) covariate following a discrete uniform distribution between 1 and 365 days.

Whenever a cluster was superimposed on a data set, the range of the covariate(s)

involved in the disease cluster was limited to a subset of its original range. For

instance, Experiment 1 refers to a disease cluster of 100 cases in covariates: t,

x1,x2, and x3. In these 100 cases, tis discrete uniform between days 150 and

160, x1is uniformly distributed between 0.3 and 0.35, x2is uniformly distributed

between 0.6 and 0.7, and x3is uniformly distributed between 0.2 and 0.4. All

other covariates in the data set (x4-x25) are uniformly distributed between 0 and

1. For the remaining cases outside the disease cluster, all covariates in the data

set (x1-x25) are uniformly distributed between 0 and 1. Similarly; Experiments

7, 8, 11, and 12 consist of two disjoint disease clusters where cluster 1 is a cluster

on tand x1, and cluster 2 is a cluster on t,x2, and x3.

In this set of experiments, Year was used as the class covariate. Ten differ-

ent replicates of each experiment were considered, and ultimately, a covariate

was selected if it had a p-value ≤0.05 in at least five replicates. Note, in Table

3.9, Experiments 9 and 10 represent null cases (i.e., no clusters are superim-

posed). However, remember that for these experiments each period in Class 1

for Experiments 1 - 8 is split into two periods of essentially half the length of

the original period. Nevertheless, the number of cases are not distributed evenly

amongst them. For example, Experiments 10 - 13 split the first 120-day period

into two 60-day periods where the first 60-day period accounts for two thirds

of the cases and the second 60-day period accounts for one third of the cases in

the original 120-day period. Consequently, by making uneven the distribution of

cases, one is indirectly making ta relevant covariate in separating Class 0 and

Class 1 cases. For this reason, Table 3.10 shows these experiments do not have

any cases in a cluster, but tstill shows up under Covariates in Cluster. On the

other hand, Period(s) with Cluster(s) in Table 3.10 indicates the exact period(s)

where a cluster was superimposed. Covariates Detected and Covariates Missing

indicate which covariates were detected (missed) using the tree-based method

where the parenthesis refers to the number of times the covariate was selected in

ten replicates.

Results in Table 3.10 show the FS method is able to successfully account

for seasonality and correctly detect the superimposed clusters. Notice the Ex-

periments 9 and 10 accurately detected tas a driver in discriminating between

Class 0 and Class 1 cases. Nevertheless, given the distribution of cases was more

uneven in Experiment 10, twas detected in all data sets for Experiment 10 as

opposed to being detected in half of the data sets for Experiment 9. Furthermore,

all experiments with at least 85 cases in the cluster(s) successfully detected all

covariates involved in the superimposed clusters. The experiments with 75 and

80 cases, however, failed to detect some of the covariates involved in the disease

clusters. As the size of the disease cluster decreases relative to the size of the

test data set, the more challenging it is for the FS method to detect such change.

In fact, these experiments represent disease clusters of less than 7% of the test

data set. As a result, when monitoring for disease clusters one might be able to

define the frequency of monitoring (i.e., the size of the test data set) based on the

smallest cluster one is interesting in detecting.

Table 3.9. Description of experiments with seasonality.

Experiment Number of Time

Periods

Number of

Clusters Periods with Clusters Cases in

Cluster

1 3 0 – 0

2 3 1 2 100

3 3 1 2 90

4 3 1 2 85

5 3 1 2 80

6 3 1 2 75

7 3 2 1, 3 100, 100

8 3 2 1, 3 75, 75

9 6 0 – 0

10 6 0 – 0

11 6 1 2 100

12 6 1 2, 4 100, 100

13 6 1 2, 4, 6 100, 100, 100

Table 3.10. Feature selection results for experiments with seasonality.

Experiment Covariates in Clusters Covariates Detected Covariates Missed

1 – – –

2t,x1,x2,x3t(10), x1(10), x2(10), x3(9) –

3t,x1,x2,x3t(10), x1(10), x2(9), x3(7) –

4t,x1,x2,x3t(10), x1(10), x2(9), x3(7) –

5t,x1,x2,x3t(6), x1(7), x2(6) x3(4)

6t,x1,x2,x3t(6), x1(6) x2(3), x3(4)

7t,x1,x2,x3t(10), x1(8), x2(7), x3(8) –

8t,x1,x2,x3t(10) x1(3), x2(3), x3(2)

9t t (5) –

10 t t (10) –

11 t,x1,x2,x3t(10), x1(10), x2(10), x3(10) –

12 t,x1,x2,x3t(10), x1(10), x2(10), x3(9) –

13 t,x1,x2,x3t(10), x1(9), x2(9), x3(9) –

3.3.5. Disease Clusters with Uneven Population

An adequate surveillance scheme should also be able to account for an uneven

population distribution across the geographical region of interest. For example,

one would expect a larger portion of population to be clustered around metropoli-

tan areas. One would also expect that a larger territorial extension implies an

increased likelihood of incidences. The surveillance scheme, thus, should be

able to successfully break down the analysis at the subregion level. It should

only trigger a signal when the number of incidences within a given subregion

is remarkably larger than expected. There should be no signals triggered when

there is a larger number of incidences in one subregion as compared to other

subregions in the study area. Using geographical coordinates, such as latitude

and longitude, often alleviates the problem. A given combination of latitude and

longitude (assuming it is recorded with several decimal places) can be as specific

as indicating the location for a given household. As a result, one would not be

concerned with contrasting incidences for subregions that vary greatly in terms

of area. However, as the geographical region of interest increases, one might

also start to worry that the characteristics of the individual households can vary

greatly across the region. On the other hand, compensating for an uneven popu-

lation is key when the geographical location of disease incidences is recorded in

a categorical format. This generally implies that disease incidences are assigned

to geographical regions that are not necessarily of the same size or characteris-

tics. A typical example of such scenario would be monitoring disease incidences

in Arizona, where geographical location is recorded at the county level. Table

3.11 shows the 2009 estimate of the population in Arizona (AZ) by county. It

also shows the land area, as given by the 2000 Census, by county (U.S. Census

Bureau, 2010). Population roughly ranges between 8,000 and 4,000,000 per-

sons. Moreover, the second largest county in terms of population has only one

fourth of the population in the largest county. Land area shows a smaller degree

of variation. It ranges from roughly 1,200 square miles to 18,000 square miles.

Table 3.12 shows the description for a series of simulated experiments based

on the population distribution by counties in AZ. Six different experiments were

simulated with 500 cases in each class. The number of cases in Class 0 for each

county is proportional to the percentage of the AZ population in each county.

The same approach was used with Class 1 cases in Experiment 0, although, a

Gaussian error component was added to make sure the number of cases for each

county in each class was not identical. Class 1 cases for Experiments 1 - 4 were

simulated as in Experiment 0, with the exception that a fixed number of cases

was selected at random and the county covariate was set equal to “Greenlee”.

Column Cases in Cluster in Table 3.12 indicates the number of cases in each

class that was selected at random and whose county covariate was set to “Green-

lee”. Notice the cluster in Experiment 5 is simulated to be generated by the

two-way interaction between age group and county. As a result, the age group

covariate was set equal to “11 - 20” in the same 150 cases that were selected at

random and whose county covariate was set to “Greenlee”.

Average detection results based on ten replicates for each experiment are

shown in Table 3.13. Again, the algorithm parameters were set to: 100 trees,

100 replicates, and the 0.9 percentile for the artificial importance scores. From

now on, an accurate detection will be used to refer to covariates that were se-

lected as important in at least five replicates, whenever they were involved in a

disease cluster. Notice the FS method did not detect any covariates in all repli-

cates associated with the null experiment. It was also able to correctly detect

county in at least five out of ten replicates for Experiments 1 - 5. Moreover, the

FS method correctly detected age group in all ten replicates for Experiment 5.

3.3.6. Scan Statistics

An alternative to the FS method is the spatial scan statistic, which is designed

to detect disease clusters only in space. Additional covariates can be incorpo-

rated into the analysis through adjustments (Kulldorff, 1997), but the adjust-

Table 3.11. Arizona population

(2009 estimate) and land area (Cen-

sus 2000) by county.

County Population Land Area

(sq. miles)

Apache 70591 11204.88

Cochise 129518 6169.45

Coconino 129849 18617.42

Gila 52199 4767.7

Graham 37045 4629.32

Greenlee 8041 1847

La Paz 20012 4499.95

Maricopa 4023132 9203.14

Mohave 194825 13311.64

Navajo 112975 9953.18

Pima 1020200 9186.27

Pinal 340962 5369.59

Santa Cruz 43771 1237.63

Yavapai 215686 8123.3

Yuma 196972 5514.09

Table 3.12. Description for experiments

with uneven population based on 500

cases in each class.

Experiment Cases in

Cluster

Covariates in

Cluster

0 0 –

1 50 county

2 75 county

3 100 county

4 150 county

5 150 county, age group

Table 3.13. Average detection for exper-

iments with uneven population.

Experiment age group county gender race

0 0 0 0 0

1 0.1 0.9 0 0.2

2 0.2 0.8 0 0.1

3 0.4 0.5 0 0.2

4 0.2 0.8 0 0.1

5 1 0.5 0 0.1

ments require additional models that are difficult to generate in high dimensions

and the method itself does not detect clusters in the additional covariates. One

might interpret the tree ensembles as intelligent scan statistics that continue to di-

vide into regions in (both space and other covariates) that hold the most promise

for differences in baseline and new data. This allows for the higher-dimensional

covariates to be explored as there is a lack of methods to handle these covariates.

An alternative way to approach this comparison would be to use the Pois-

son model of the scan statistic. As opposed to the Bernoulli model (that needs

baseline data and new test observations), the Poisson model requires new test

observations in the form of counts for each zip code as well as a population

file. One way to generate a population file is to directly input the population in

each zip code. Another way would be to input the expected counts for each zip

code after indirectly adjusting for additional covariates. Typically, this indirect

standardization is done using regression modeling. One predicts the expected

counts for a given zip code for a given combination of the additional covari-

ates. For numerical covariates, one would have to strategically make predictions

at different points of the covariate range and figure out a reasonable approach

to combine these expected counts. For categorical covariates, one would have

an expected count for each possible combination of values and, then, aggregate

these predictions over the entire range. As the number of categorical covariates

and the number of levels for each categorical covariate increases, the number of

expected counts that need to be aggregated also increases. This aggregation of

counts makes it obvious that it would be relatively easy to miss an early signal

where new incidences increase in a very small number of these expected counts

for a zip code. The overall expected count for the zip code will be relatively

insensitive to this increase.

3.4. Power Curves

Figures 3.6 - 3.7 show the power curves for the proposed method. Both

graphs are based on data sets with 3,000 baseline cases and fifty uniformly dis-

tributed covariates. The series in Figure 3.6 refer to the number of cases in Class

1 and the range for the single covariate in the cluster (x50). For instance, the

first series with Class 1 = 250 and Range = 0.05 refers to an experiment where

there are 250 test cases and the range for the covariate in the cluster has been

reduced from 0 - 1 to 0 - 0.05. In general, the larger the number of cases in the

cluster, the easier it is for the proposed method to detect a change in the dis-

ease incidence rate. Also, the smaller the number of cases in Class 1 and the

smaller the range, the easier it is for the proposed method to detect any changes.

Notice each point in the graph represents the average number of times x50 was

detected in thirty independent replicates. Notice also false alarms are not a ma-

jor concern as, on average, x50 was detected 1.11% of the time when none of its

cases were in a cluster. Lastly, observation labels refer to the standard deviation

s=ppx50 (1 −px50 )/30 ∗100, where px50 refers to the proportion of replicates

in which covariate x50 was detected.

Figure 3.7 shows the power curves for a two-way interaction between covari-

ates x49 and x50. For simplicity, the percentage of detection for each covariate

over thirty independent replicates was shown in separate graphs. In general, one

can observe the same behavior as in the power curve for a main effect. Smaller

number of Class 1 cases, smaller ranges, and larger number of cases in cluster

lead to improved detection. In addition, very few false alarms are observed.

3.5. Summary

Point processes occur in numerous disciplines and change detection in such

data has been limited to low-dimensional solutions. Here, supervised learners

are used that can handle high-dimensional data of mixed numerical and cate-

gorical types, strong interaction effects, seasonality, and even with an unevenly

distributed population. Baseline observations (Class 0) are contrasted against

new data (Class 1) to look for those covariates that allow one to better discrimi-

nate between the classes and, in turn, generate the signal of a change.

A clear advantage of my method against alternatives, such as the spatial scan

statistic, is its ability to identify changes in all covariates. Alternatives have been

limited to a few covariates and, as data and computationally power increases,

there is a need to look throughout the data for signals. The proposed approach

to first identify covariates that contribute to a signal (if any) is expected to often

reduce the complexity of the problem to many fewer variables. Then, simpler

methods can isolate the locations of changes.

Although tree-ensembles have the advantage that they can explore high-

dimensional space and detect interactions, other supervised learners can be ap-

plied. The proposed learner provide a general approach to the problem, but other

learners and their parameter setting could be tuned to be most sensitive to par-

ticular problem characteristics. This generic approach is reasonable for many

applications, but tuning issues can be used to improve the performance, and is-

sues, such as sample size, are dependent on the type of learner and signal criteria

adopted. For example, if the new data comprises a much smaller proportion than

the baseline data, one might want to adjust prior probabilities in the classes. Al-

ternatives learners would use alternative feature selection methods (or alternative

criteria) to generate signals.

(a) (b)

Fig. 3.3. Kernel intensity functions for: (a) null experiment; (b) one cluster of

200 cases in space; (c) one cluster of 100 cases in space and age, and (d) two

clusters of 100 cases each in a space-age and space-race interaction.

(a) (b)

Fig. 3.4. Kernel intensity functions for Experiment 5 by: (a) age and (b) race.

(a) Seasonal daily counts without disease clus-

ters.

(b) Seasonal daily counts with a single disease

cluster (days 152-181).

Fig. 3.5. Example with seasonal disease incidences.

Fig. 3.6. Power curve for a main effect with fifty uniformly distributed covariates

between 0 and 1. Results on the y-axis indicate the percentage of times the

covariate in a cluster (x50) was detected in thirty replicates. Data labels indicate

the standard deviation for the detection percentage. Different series refer to the

total number of cases in Class 1 and the reduced range for the variable in a

cluster.

(a) First variable in two-way interaction.

(b) Second variable in two-way interaction.

Fig. 3.7. Power curve for a two-way interaction (between variables x49 and x50)

with 50 uniformly distributed variables. Results on the y-axis indicate the per-

centage of times the covariate in a cluster was detected in thirty replicates. Data

labels indicate the standard deviation for the detection percentage. Different se-

ries refer to the total number of cases in Class 1 and panel labels refer the reduced

range for the covariates in a cluster.

100

CHAPTER 4

ALTERNATIVE METHOD FOR SURVEILLANCE

4.1. Introduction

This chapter discusses the Percentile method as an alternative to the Feature

Selection (FS) method. Similar to the FS method, the Percentile method uses the

transform to supervised learning to detect changes in disease incidence rates. In

addition, it uses tree-based ensembles as the classifier of choice. The Percentile

method relies strictly on parallel tree-based ensembles, as opposed to the FS

method, which combines both parallel and tree-based ensembles. As it is the

case with the FS method, some of the advantages of using tree-based classifiers

include: (1) the form and size of a disease cluster(s) does not need to be pre-

specified, (2) one or multiple disease clusters and high order interactions can be

detected, and (3) missing values and mixed type covariates can be handled.

On the other hand, the Percentile method detects changes in disease inci-

dence rates using an alternative signal criterion based on the probability of clas-

sification into Class 1 for tree-based ensembles. The intent of developing the

Percentile method is to provide a less computationally intensive approach for de-

tecting changes in disease incidence rates, while still taking advantage of surveil-

lance’s transform to supervised learning.

101

Section 4.2 describes the new signal criterion for the Percentile method. Sec-

tion 4.3 describes a series of simulated experiments with a variety of superim-

posed disease clusters. Section 4.4 shows the results for the Percentile method

using different parameters and compares it to the FS method. Finally, Section

4.5 includes concluding remarks.

4.2. Methods

The Percentile method detects changes in disease incidence rates by identi-

fying the specific incidences that most significantly differ from baseline inci-

dences. Another way to see the problem is to identify all the incidences in those

subregions with an increased number of new test incidences. For example, re-

visit Figure 3.1. Most of the incidences from the two different time periods seem

to behave similarly. Nonetheless, there is a subregion dominated by new test in-

cidences. The Percentile method, thus, presents a non-computationally intensive

method for isolating the specific new test incidences that behave significantly

different from baseline incidences. In other words, it isolates the specific inci-

dences that lead to an increase in the disease incidence rate of interest.

Suppose there is a random forest of mtrees (Breiman, 2001), where each tree

was built on a different bootstrap sample. As a bootstrap sampling is done with

replacement, some incidences will be in the sample multiple times, while others

102

will be excluded. Incidence iis considered out-of-bag in tree t∈ {1, . . . , m}

when it is not in bootstrap sample t. Let n∗

ibe the number of times incidence i

was out-of-bag across the entire ensemble, then, the probability of classification

into Class 1 for incidence ibecomes

Pi(Y= 1) = 1

n∗

j=1

ˆyij.(4.1)

The probabilities of classification into Class 0, P(Y= 0), and Class 1,

P(Y= 1), are complimentary in binary classification problems. Only the prob-

ability of classification into Class 1, P(Y= 1), will be considered from now

on. Groups of test (Class 1) incidences that have remarkably large P(Y= 1)

indicate the subregions in the multidimensional space defined by patient-level

covariates where the new test incidences behave differently from baseline data.

The Percentile method, thus, uses the estimate of P(Y= 1) to trigger a sig-

nal when there is an indication that these unusual incidences indeed represent

disease clusters. Please note the terms incidences and cases will be used inter-

changeably throughout the remaining chapters.

Consider the algorithm for the Percentile method (summarized in page 105).

The first step is to assign all nBbaseline cases to Class 0 and all nTnew test

cases to Class 1. Next, draw Jreplicates of a random forest of mtrees. At each

replicate, record P(Y= 1) for each one of the (nB+nT)cases. Also, compare

103

the P(Y= 1) value for each one of the nTtest cases against the P(Y= 1)

values for all nBbaseline cases. Assign case i, i ∈ {1, . . . , nT}, at replicate

j, j ∈ {1, . . . , J}, a rank of rif it is among the rlargest P(Y= 1) values

observed in the baseline cases at replicate j. Transform the ranks into simulated

p-values using the following formula:

simulated p-value =Rankij

.(4.2)

Repeat the process at each replicate j, j ∈ {1, . . . , J}, until there are J

simulated p-values for each case. Then, calculate a percentile of the simulated

p-values for each case i, i ∈ {1, . . . , nT}, and compare them against a threshold.

If at least γ%of the nTtest cases have a median simulated p-value less than

or equal to the threshold, trigger a signal indicating there has been a change

in the disease incidence rate. Notice γis a function of the minimum size of

the cluster that one is interested in detecting and the number of cases in Class

1. For instance, suppose there are 500 cases in Class 1 and one is interested

in detecting emerging disease clusters that consist of at least 10 cases. In that

case, set γ= 10/500 = 0.02. In Section 4.4, different alternatives for potential

threshold values as well as the choice of the percentile for the simulated p-values

will be discussed.

104

Percentile Method:

1) Assign nBbaseline cases to Class 0.

2) Assign nTtest cases to Class 1.

3) FOR j= 1 TO Jreplicates

a. Build a random forest of mtrees.

b. Calculate P(Y = 1) for all (nB+nT)cases.

c. FOR i= 1 TO nT

i. Compare P(Y= 1) of test case iagainst P(Y= 1) for all nBbaseline cases at

replicate j.

ii. Calculate a simulated p-value for test case iusing Rankij / nBwhere Rankij =r

when P(Y= 1) of case iis among the rlargest P(Y= 1) values observed in the

nBbaseline cases at replicate j.

4) Signal when at least γ% of test cases have a median(simulated p-value)≤threshold.

4.3. Simulated Experiments

This section tests the effectiveness of the method on eleven different experi-

ments with mixed type covariates and clusters of different sizes. All experiments

have 500 cases in each class. Ideally, one has more baseline cases than testing

cases; and if that is the case, prior probabilities should be set to 0.5. That way

classification will not be biased toward the class with the larger number of inci-

dences.

Refer to the experiments described in Table 4.1. The numerical covariates in

these experiments are simulated as a planar projection of geographical coordi-

nates within the state of Arizona. The categorical covariates for each experiment

105

consist of a combination of covariates with five, seven, and ten levels. Moreover,

there are ten different replicates for each experiment.

Experiment 0 represents the null experiments. It has two numerical covari-

ates and twenty five categorical covariates, and none of the Class 1 cases are in a

cluster. The numerical covariates are used to represent a planar projection of the

geographical coordinates associated with the patient’s address. The twenty five

categorical covariates are simulated in the following manner: ten covariates with

five levels, ten covariates with seven levels, and five covariates with ten levels.

Experiment 1 is very similar to Experiment 0 with the exception that it has one

superimposed cluster. The cluster is simulated to be generated by the interaction

between geographical coordinates and three categorical covariates in 100 out of

the 500 cases in Class 1. For example, focus on one of the replicates for Ex-

periment 1. Let {x, y}represent projected geographical coordinates. The range

for xin cases outside the cluster is [-0.0487, 0.0487]; while it becomes narrower

for the 100 cases within the cluster, [-0.0477, 0.0453]. The same happens with

y. Its range for cases not in the cluster is [0.5473, 0.6457] and [0.5484, 0.6439]

for cases within the cluster. Similarly, values for categorical covariates in the

cluster are simulated to be a random sample from a subset of the levels for each

covariate. Covariates C1, C11, and C21 are involved in the simulated cluster.

106

Table 4.1. Description for Percentile method experiments.

Experiment Cases in

covariates

Numerical

covariates

Categorical

in Cluster

Numerical

covariates

Categorical

covariates

0 0 2 25 0 0

1 100 2 25 2 3

2 100 2 15 2 3

3 100 2 10 2 3

4 50 2 10 2 3

5 100 2 10 2 0

6 100 2 10 2 1

7 100 2 10 0 1

8 100 2 10 0 2

9 200 2 10 0 1

10 200 2 10 0 2

C1 originally has five different levels, but the cases in the cluster all belong to

the same level. In the same manner; the levels for covariates C11 and C21 are

reduced from seven and ten levels to two and three levels, respectively, for all

cases in the cluster.

4.4. Results

Figure 4.1 shows the kernel density estimates for the distribution of different

P(Y = 1) percentiles for the same Experiment 1 replicate discussed in detail

in Section 4.3. Clearly, the behavior in all these plots is almost the same. The

distribution for Class 0 cases and Class 1 cases not in a cluster is reasonably sim-

ilar, while the distribution for cases in a cluster is shifted toward the right. This

validates the hypothesis that cases in a disease cluster will have, on average, a

larger P(Y= 1) than cases that behave according to the base incidence rate. As

a result, either one of these plots could essentially provide the same information.

107

(a) median (b) third quartile

Fig. 4.1. Different P(Y = 1) percentile distributions for one replicate of Exper-

iment 1.

108

For the sake of simplicity, the signal for the Percentile method is based on the

median of the P(Y= 1) simulated p-values.

Tables 4.2 - 4.4 show the detection percentages based on the Percentile

method using 100 replicates of random forests with 100, 250, and 500 trees,

respectively. Six different thresholds were considered. For the hard thresholds;

0.03, 0.04, and 0.05 were used. On the other hand; remember in Step 4(b) of the

Percentile Method (summarized in page 105), the median simulated p-value for

each incidence in Class 1 was calculated. The 0.01, 0.015, and 0.02 percentiles

of this distribution were used as soft thresholds. Overall, the method triggered a

signal whenever there were at least 10 Class 1 cases (γ= 2%) with a simulated

p-value less than or equal to the threshold of choice. The detection percentage

is given as the percentage of signals triggered for each experiment in ten inde-

pendent replicates. Tables 4.2 - 4.4 also show the detection percentage for the

Feature Selection (FS) with parameters: 100 trees, 100 replicates, and the 0.9

percentile for the artificial importance scores. Using the FS method, a signal

was triggered whenever there was at least one covariate with a p-value less than

or equal to 0.05 divided by the total number of input covariates in the data set.

As in Chapter 3, an accurate detection is defined as a signal triggered in at least

five out of ten different replicates for those experiments with a superimposed

109

Table 4.2. Signal detection using 100 trees and 100 replicates.

Experiment Hard Thresholds Soft Thresholds FS Method

0.03 0.04 0.05 0.01 0.015 0.02 p-val ≤0.05/nVars

0 0 0 0 0 0 10.1

1 1 1 1 0.1 0.5 1 0.9

2 1 1 1 0.1 0.5 1 0.8

3 1 1 1 0.3 0.7 1 1

4 0.6 0.8 0.9 0 0.2 1 0.9

5 1 1 1 0 0.1 1 1

6 1 1 1 0.2 0.6 1 1

7 0 0.1 0.3 0 0 1 0.8

8 0.9 1 1 0 0.1 1 1

90.4 0.7 0.8 0 0.1 1 1

10 1 1 1 00.6 1 1

Table 4.3. Signal detection using 250 trees and 100 replicates.

Experiment Hard Thresholds Soft Thresholds FS Method

0.03 0.04 0.05 0.01 0.015 0.02 p-val ≤0.05/nVars

0 0 0.1 0.1 0 0 10.1

1 1 1 1 0.5 0.8 1 0.9

2 1 1 1 0.6 1 1 1

3 1 1 1 0.5 0.9 1 0.9

4 1 1 1 0.3 0.5 1 1

5 1 1 1 0.1 0.5 1 1

6 1 1 1 0.9 0.9 1 0.8

70.3 0.5 0.6 0 0 1 1

8 1 1 1 0.2 0.4 1 1

9 0.8 1 1 0 0 1 1

10 1 1 1 0.1 0.3 1 0.8

cluster. An accurate detection for the null case (Experiment 0), on the contrary,

is defined as less than five signals triggered in ten replicates.

Highlighted in bold face font are the inaccurate detections in Tables 4.2 - 4.4.

In general, the hard thresholds lead to the fewest number of inaccurate detec-

tions. They are also independent of the data and much easier to implement than

soft thresholds. Nevertheless, there appears to be issues with the sensitivity of

the proposed method when built around random forests of 100 trees. The results

110

Table 4.4. Signal detection using 500 trees and 100 replicates.

Experiment Hard Thresholds Soft Thresholds FS Method

0.03 0.04 0.05 0.01 0.015 0.02 p-val ≤0.05/nVars

0 0.1 0.3 0.8 0 0 10.1

1 1 1 1 0.4 0.8 1 0.9

2 1 1 1 0.4 0.5 1 0.8

3 1 1 1 0.7 0.9 1 1

4 1 1 1 0.2 0.4 1 0.9

5 1 1 1 0.3 0.6 1 1

6 1 1 1 0.9 0.9 1 1

7 0.5 0.7 1 0 0.1 1 0.8

8 1 1 1 0 0.4 1 1

9 1 1 1 0 0.2 1 1

10 1 1 1 0.2 0.5 1 1

for the 0.4 and 0.5 hard thresholds using random forests of 250 and 500 trees,

on the other hand, are highly comparable to those obtained with the more com-

plex FS method. Indeed, one can recommend the Percentile method using 100

replicates of random forests with 250 trees each and a 0.05 hard threshold. This

combination of parameters leads to an incorrect signal in 4 / 110 (11 experiments

with 10 replicates each) runs, which is comparable to the 6 / 110 incorrect signals

associated with the FS method. The system time for the proposed method based

on the combination of parameters suggested is, on average, ¯x= 6.63 minutes

with standard deviation s= 0.48 minutes when implemented in R (R Develop-

ment Core Team, 2009). Notice this is roughly half the computational time of

the model with random forests of 500 trees (¯x= 11.2min,0.9min), which leads

to highly similar sensitivity with a more complex model.

Table 4.5 considers different number of replicates (J)for random forests

111

of 250 trees. This table shows the complexity of method can be further de-

creased without decreasing the method’s sensitivity. The detection percentages

are roughly the same irrespective for all number of replicates (J). Nevertheless,

the system time can be considerably decreased when using fewer replicates. As

a consequence, one can further recommend using 30 replicates of random forests

with 250 trees each. This combination further reduces the average system time

to ¯x= 1.9minutes with standard deviation s= 0.48 minutes without reducing

the method’s sensitivity.

The system time for the experiments described in Table 4.1 using the FS

method is ¯x= 5.28 seconds with standard deviation s= 1.42 seconds. Clearly,

this appears counterintuitive as one is proposing a less computationally intensive

method. The issue is that the FS method was implemented in C++, while the

Percentile method was implemented in R (R Development Core Team, 2009).

Clearly, the system time for the Percentile method could be significantly im-

proved if implemented in C++. Irrespectively of which language is chosen, the

lead time for implementing the Percentile method is much shorter.

4.5. Summary

The Percentile method is a simpler alternative to the FS method. The idea

is that good feature selection algorithms might not be readily available, and the

112

Table 4.5. Signal detection for Percentile method using 250 trees and a 0.05 hard

threshold for different number of replicates (J).

Experiment Replicates (J)

30 50 75 100

0 0.1 0.1 0.2 0.1

1 1 1 1 1

2 1 1 1 1

3 1 1 1 1

4 1 1 1 1

5 1 1 1 1

6 1 1 1 1

7 0.6 0.6 0.7 0.6

8 1 1 1 1

9 1 1 1 1

10 1 1 1 1

Mean System Time (min) 1.90 3.22 4.88 6.63

Std System Time (min) 0.15 0.23 0.36 0.48

user might require extensive programming experience to implement a feature

selection algorithm such as the one used in the FS method (Tuv et al., 2009). The

Percentile method alleviates the problem by developing a non-computationally

intensive method that could be easily implemented even in free software such as

R. Both methods show highly satisfactory results for all simulated scenarios. In

fact, both methods lead to similar results. The system time requirements for the

FS is still much lower; yet, the methods are implemented in different languages.

It is expected that the system time requirements for the Percentile method can

be substantially reduced if implemented in the same language as the FS method.

Again, both the FS and the Percentile methods are not exclusive to public health

surveillance. In fact, they are general methodologies that can be applied to detect

changes in rates of events in many other disciplines. One of the simplest cases

113

being detecting out-of-control performance in a statistical process control.

114

CHAPTER 5

RULE-BASED FEATURE VALUE SELECTION

5.1. Introduction

Good feature selection methods are often successful at identifying the impor-

tant covariates in medium to large data sets. Identifying the important covariates,

however, is not always sufficient. Sometimes, it is equally important to fully un-

derstand which values or levels of a covariate contribute to its importance. For

instance, in public health surveillance, it is very useful to know whether there

are any disease clusters in the data. However, it is equally important to identify

the specific geographical regions and/or subpopulations where these clusters are

occurring. In statistical process control, it is also of interest to identify the spe-

cific combination of product characteristics that contribute to an out-of-control

signal. This problem can be referred to as feature value selection, where one

identifies the specific levels and/or ranges of the important covariates that are

leading to relevant changes in the data.

Suppose there is a binary classification problem where Class 0 refers to base-

line data. Again, baseline data refers to data that behaves according to the nat-

ural realization of a given process. In the problem of public health surveillance,

this means data that behaves according to the base incidence rate of the disease.

In statistical process control, baseline data can be related to non-defective pro-

115

duction. On the other hand, Class 1 refers to new test observations that are used

to monitor the current state of the process. One wants to know whether new dis-

ease clusters have appeared or if a manufacturing process became out-of-control,

and these new test observations are the key to unveiling this information. If data

from the different classes (different time periods) behaves similarly, it means

that no changes have occurred in the process. On the contrary, if a classifier can

successfully discriminate between Class 0 and Class 1, it means that there is a

specific combination of covariate values that leads to a change in the process,

and developing a method for identifying this specific combination is the aim of

this chapter.

In this chapter, rule-based classifiers are used to address feature value selec-

tion on the binary classification problem previously described. Rule-based clas-

sifiers based on mining algorithms, such as Apriori (Tan et al., 2006; Hahsler

et al., 2010), can only be used on categorical covariates. Section 5.2 addresses

this issue with a simple algorithm that can be used to automate the discretization

technique for continuous covariates proposed by Berrado and Runger (2009).

Section 5.3 briefly describes how to use rule-based classifiers in feature value

selection. Sections 5.4 - 5.5 cover some experiments and results that evaluate

the effectiveness of the rule-based feature value selection method on simulated

116

experiments. Finally, Section 5.6 summarizes the findings.

5.2. Random Forest Discretization

An important problem that arises when applying rule-based classifiers to

continuous attributes is how to determine an adequate discretization scheme that

splits the observed covariate range(s) into categories, at the same time, loss of

information is minimized. Berrado and Runger (2009) propose a discretization

method based on the covariate splits used by a random forest. In their paper,

however, they do not provide a means to automate the algorithm. They rely on

visual inspection of a histogram of the random forest splits for each continuous

covariate. In this work, a simple algorithm is used to automate the process

of splitting the range of a continuous covariate into intervals by means of

identifying the large frequency intervals or groups of intervals in a histogram.

Most important, this simple automation scheme greatly facilitates the use of

the method by Berrado and Runger (2009) in large-scale applications where

multiple runs of rule-based classifiers might be needed as inputs to a wide array

of additional methods. The proposed automation for the method by Berrado and

Runger (2009) is described next.

117

Random Forest Discretizing Algorithm:

I. Build random forest RF using mtrees.

II. SET numVars =number of numerical covariates in data set.

III. FOR (i= 1 TO numVars)

A. Create vector SPibased on all the split points for numerical covariate iin RF .

B. Create histogram Hiof SPiusing ni=dlog2|SPi|+ 1eequal-sized intervals.

C. SET interval j= 1.

D. WHILE (j <=ni)

1. FOR (k= 4 TO 1{BY −1})

a) Consider interval jand the following k−1intervals in histogram H.

b) IF ( percentage difference among the intervals within the group is less than γw

and greater than γbwith respect to the intervals outside the interval group and the

frequency of each interval in the group is greater than or equal to F)

(1) Add the midpoint of the kgroup interval to breaks.

(2) j=j+k

(3) EXIT FOR

c) ELSE k=k−1

E. IF (breaks ={})

1. Split the covariate into pequal-sized intervals.

2. NEXT j

F. Discretize numerical covariate iaccording to breaks.

Steps I−III(a)refer back to the method by Berrado and Runger (2009). Start

by building a random forest using m= 100 trees. For each numerical covari-

ate, create a vector that contains all the split points whenever the continuous

covariate of interest was used as the primary splitter in the forest. Next, plot

each vector of split points in a histogram with ni=dlog2|SPi|+ 1eequal-

sized intervals, where |SPi|refers to the length of the vector of split points (SPi)

118

for continuous covariate i. For each interval, consider every possible group of

k, k ∈ {1,...,4}, intervals. Starting with interval 1, consider every possible

group of k= 4 contiguous intervals. If none of these groups satisfies the con-

ditions for group importance, then, consider every possible group of k−1=3

contiguous intervals and so on. At any point, if a kgroup interval satisfies the

conditions for group importance, do not consider any additional group of inter-

vals for any of the intervals included in the detected group. Figure 5.1 clearly

shows what is meant by a group of k= 3 contiguous intervals including the

interval of interest. In this example, assume interval 3 is the interval of interest.

Hence, groups of k= 3 contiguous intervals for interval 3 can be formed by

including interval 3 along with: (1) intervals 4 - 5 (two intervals above it), (2) in-

tervals 1 - 2 (two intervals below it), and/or (3) intervals 2 and 4 (one interval

below and above it). The same approach is used for other kvalues. Of course,

there are more possible combinations for groups of k= 4 contiguous intervals,

and less possible combinations for groups of k≤2contiguous intervals as well

as for those intervals toward the lower end or upper end of the histogram. No-

tice Step D in Figure 5.2 discusses only forward comparisons for the sake of

simplicity in the summarized description.

Now, one can focus on the conditions that make a group of intervals impor-

119

(a) Group including interval 3 and two

intervals above it (intervals 4-5).

(b) Group including interval 3 and one

interval above and another interval be-

low it (intervals 2 and 4).

intervals below it (intervals 1-2).

Fig. 5.1. Interval importance evaluation considering k= 3 groups.

tant. Suppose there is a group of kcontiguous intervals including the interval of

interest. The idea is that a group of intervals is important if it satisfies three con-

ditions. First, the frequency in each interval must exceed a pre-defined threshold

F. Second, the percentage difference among frequencies of intervals within the

group should not exceed γw, where 0< γw<1. Finally, the percentage differ-

ence between the left most and right most intervals in the group with respect to

the intervals preceding and following the interval group, respectively, should be

greater than or equal to γb, where 0< γb<1and γb> γw.

Consider Figure 5.2. The group of k= 3 contiguous intervals being consid-

ered includes intervals 3 - 5. Suppose the selected threshold for this example

120

is set to F= 70. Clearly, all intervals in the group would immediately satisfy

Condition 1. The threshold Fcan be substantially smaller than 70 depending

on the application. Indeed, 5≤F≤10 is often reasonable in public health

applications. In order to verify whether Condition 2 is true, do k−1=2for-

ward pairwise comparisons. The percentage difference between the frequency

of interval 3 (f3)with respect to the frequency of interval 4 (f4)must satisfy

|(f3−f4)/f3| ≤ γw. Also, the percentage difference between the frequency

of interval 4 (f4)with respect to the frequency of interval 5 (f5)must satisfy

|(f4−f5)/f4| ≤ γw. Since intervals within a group are intended to be simi-

lar, values of γwbetween 0 and 0.1 are often reasonable. Finally, Condition 3

considers a comparison at each end of the interval group. If the group is indeed

important, there should be an increase when going from the interval preceding

the group to the first interval in the group. In addition, there should be a decrease

when going from the last interval in the group to the next interval. As a result,

(f3−f2)/f3≥γband (f5−f6)/f5≥γbmust be satisfied.

The value of γbcontrols the amount of split points one wants to detect. If γbis

close to 1, very few to no split points will be detected. On the other end, a value

of γbclose to zero implies almost every interval will be detected. In most of the

experiments presented in this chapter, values of γbbetween 0.10 and 0.15 lead

121

to satisfactory results. Note that the two-sided comparison is only possible for

intervals toward the middle of the histogram. The first and last intervals in the

histogram, for instance, can only be compared to intervals above it and below

it, respectively. When looking backward to compare groups of intervals where

the interval of interest is the last interval, the proposed method also takes into

consideration not to include those intervals which have already been detected as

part of another group of intervals.

In the example shown in Figure 5.3, two groups of intervals are detected

using F= 5,γw= 0.075, and γb= 0.10. One split point would be selected

for the group of intervals 3 - 5 and another one for the group of intervals 7 -

9. Since intervals 4 and 5 are detected when evaluating interval 3, there is no

need to perform further evaluations involving intervals 4 or 5. One would also

stop all comparisons after detecting the group of intervals 7 - 9 when evaluating

interval 7. Given the two detected groups of intervals, the continuous covariate

associated with the histogram in Figure 5.1(a) would be discretized using their

midpoints as follows: [0,0.39),[0.39,0.83), and [0.83,1). Each observed value

for the continuous covariate of interest is, now, replaced by the category defined

by its associated range. In the case the proposed method is not able to detect any

important intervals, one can simply divide the continuous covariate range into p

122

equal-sized intervals. The number of default intervals is typically set to p= 20.

Fig. 5.2. Comparison for intervals within and between kinterval group.

Fig. 5.3. Example of continuous covariate discretization.

5.3. Rule-Based Feature Value Selection

Once all covariates are in categorical format, one can proceed to perform

feature value selection using rule-based classifiers. In fact, interest is strictly

centered in mining the subset of rules with Class 1 as consequent. These rules

are the ones that allow identification of covariate ranges and/or levels that allow

one to discriminate between the two different classes. In other words, these rules

explain which covariate values (or combination of values) lead to a deviation

123

from the baseline. The Rule-Based Feature Value Selection (RBFVS) method is

briefly summarized next.

Rule-Based Feature Value Selection (RBFVS) Method:

1) Discretize each numerical covariate in data set Dusing the Random Forest Discretization Algo-

rithm.

2) Create data set D0by replacing the original continuous covariates with their discretized coun-

terparts.

3) Create a binary incidence matrix from data set D0.

4) Use the Apriori algorithm on the binary incidence matrix to mine all rules with support ≥Sand

confidence ≥C.

5) Extract all rules with Class 1 as its consequent.

6) Important covariate values appear in the antecedent of the rules found in Step 5.

First, discretize all continuous covariates according to the method described in

Section 5.2 (Steps 1 - 2). Next, create a binary incidence matrix from the updated

data set based on all categorical covariates. Each covariate range and/or level in

the updated data set becomes a column (or item) in the binary incidence matrix.

Values of one in the binary incidence matrix relate to items that were included

in the original row vector, while values of zero relate to the lack of items in

the original row vector. For example, the first row vector in Table 5.1 becomes

[X1== [0,0.3)] = 1,[X2== L1] = 1,[X3== C1] = 1,[Class == 0] =

1, and zero elsewhere in the binary incidence matrix shown in Table 5.2. The

remaining rows in the binary incidence matrix are calculated in the same manner.

124

Table 5.1. Updated data set with discretized numerical covariates.

x1x2x3Class

[0,0.3) L1 C1 0

[0.3,0.6) L2 C1 1

[0,0.3) L3 C3 1

[0.06,1) L1 C2 0

[0.3,0.6) L2 C4 0

Table 5.2. Binary incidence matrix.

Itemset Items

x1== [0,0.3) x1== [0.3,0.6) x1== [0.6,1) x2== L1x2== L2x2== L3

11 0 0 1 0 0

20 1 0 0 1 0

31 0 0 0 0 1

40 0 1 1 0 0

50 1 0 0 1 0

Itemset Items

x3== C1x3== C2x3== C3x3== C4Class == 0 Class == 1

11 0 0 0 1 0

21 0 0 0 0 1

30 0 1 0 0 1

40 1 0 0 1 0

50 0 0 1 1 0

The following step in the RBFVS method uses the Apriori algorithm to mine

all possible rules with support ≥Sand confidence ≥C. In the experiments

that will be discussed in the following sections, the number of row vectors that

contain a particular itemset is set to 5% (S≥0.05), and the reliability of the

inference made by any rule is set to 62.5% (C≥0.625). From this set of rules,

extract those rules with rule consequent (or right-hand side) equal to Class 1.

Identify all those categories or ranges that appear in the antecedent of the rules

extracted in Step 5 as the covariate values responsible for the observed change

in the process.

125

5.4. Experiments

Table 5.3 describes eleven simulated experiments used to evaluate the effec-

tiveness of the RBFVS method. For each experiment, there are 500 cases in

Class 0 and 500 cases in Class 1. Experiments 0 - 5 consist of ten mixed-type

covariates. The first five covariates in these data sets are continuous covariates.

Covariates 1 - 2 are simulated as geographical coordinates within the state of

Arizona (AZ). Covariates 3 - 5 are simulated as uniform covariates between 0

and 1. The remaining covariates are categorical with the number of levels vary-

ing between five and ten. Experiments 6 - 10 consist of twenty five numerical

covariates. The first ten covariates are continuous, where the first two are once

again simulated as geographical coordinates and the remaining eight continuous

covariates are uniformly distributed between 0 and 1. Similarly, the remaining

fifteen covariates are simulated as categorical with the number of levels varying

between five and ten.

Experiment 0 refers to the null experiment. There have been no changes in

the process, so the new test data follows the same distribution as the baseline

data. Experiment 1, on the other hand, has a cluster in one numerical covariate.

In this case, the range for a uniform covariate has been reduced from [0,1] to

[0,0.15] for 100 cases in Class 1. Experiment 2 has a cluster in a categorical co-

126

Table 5.3. Experiment Description.

Experiment Numerical

covariates

Categorical

covariates

Numerical covariates

in Cluster

Categorical

covariates in Cluster

0 5 5 0 0

1 5 5 1 0

2 5 5 0 1

3 5 5 1 1

4 5 5 1 2

5 5 5 2 1

6 10 15 1 0

7 10 15 0 1

8 10 15 1 1

9 10 15 1 2

10 10 15 2 1

variate where, for 100 cases in Class 1, the number of observed levels has been

reduced from five levels to two levels. Experiment 5 has a cluster in the two ge-

ographical coordinates where a geographical location is first selected at random,

and, then, a circular cluster of fixed radius is superimposed in the data. This

cluster also involves the interaction of these geographical coordinates with one

categorical covariate for which the cluster has been defined in the same manner

as in Experiment 2. The remaining experiments can be described accordingly.

Notice, also, Experiments 6 - 10 mirror Experiments 1 - 5 in the larger dimen-

sional space. Experiment 6 is the equivalent of Experiment 1 in 25 dimensions

and so on.

5.5. Results

In this section, the RBFVS method is compared against the ACE method.

An overview of the ACE method is described in page 130. For more details on

127

the t-test comparison described in page 130, refer to Tuv et al. (2009). Note

that the ACE feature value selection method can only be applied to categorical

covariates. As a result, the data sets for both RBFVS and ACE are discretized

using the Random Forest discretization method. For experiments with five and

ten continuous covariates, the system time for the Random Forest discretization

method was ¯x5= 79.45 (¯x10 = 207.14) seconds on average with a standard

deviation of ¯s5= 4.17 (¯s10 = 12.55) seconds on an Intel Core i5 processor.

Results for ten independent replicates of each experiment described in Table

5.3 are summarized in Table 5.4 where the terms sensitivity,specificity, and d

are defined as follows:

Sensitivity =Number of Correctly Detected Covariates

Total Number of Covariates in a Cluster

Specificity = 1 −Number of Covariates Detected in Excess

Total Number of covariates Not in a Cluster

d=p(1 −Specificity)2+ (1 −Sensitivity)2.

(5.1)

A successful method is expected to have large sensitivity and specificity. Nev-

ertheless, improving the sensitivity of a method often leads to a decrease in its

specificity. As a result, the combined criteria din Equation 5.1 is used to equally

weight the sensitivity and specificity criteria in order to compare the performance

128

of the RBFVS and ACE methods. The RBFVS method was implemented us-

ing the Apriori implementation by Borgelt (2004) and the following set of pa-

rameters: m= 100 trees, F= 5 observations, p= 20 equal-sized intervals,

γw= 0.075,γb= 0.1, support(S)=0.05, and confidence(C) = 0.0625. For the

ACE method, the parameters used are: number of replicates = 100, number of

trees = 100, artificial contrast quantile(q)= 0.9, and t-test p-value = 0.05. For

more information on the ACE parameters, refer to Tuv et al. (2009). Finally, the

level estimation mode in the ACE method was set to one vs. the rest as described

in the method’s summary in page 130.

129

Artificial Contrasts Ensemble (ACE) Feature Value Selection:

I. Build random forest RF using mtrees.

II. Create an artificial covariate for each covariate in the data set by permuting its observed values.

III. SET catV ars = total number of categorical covariates in data set.

IV. FOR i= 1 TO catV ars

A. SET C=total number of covariate icategories

B. FOR j= 1 TO C

1. Calculate category importance score using one of the following:(one vs. the rest,one

vs. the rest random, or influence estimation).

a) One vs. the Rest

(1) FOR t= 1 TO m

(a) Importance Scoreijt =PNt

n=1 {Split Weightnt ∗∆(Gini score)nt }

where:

Split Weightnt =





1if split made on var i, cat jat node n using tree t

0otherwise

(b) NEXT t

(2) Use a paired ttest to compare the Importance Score of covariate i, category j

against a large percentile from the artificial covariate importance scores across

RF .

b) One vs. the Rest Random

(1) When building RF , at each node, select a random sample of covariates.

(2) For each covariate in the sample, select a split a random.

(3) Select the random split that leads to the largest ∆(Gini score)nt.

(4) Calculate and compare importance scores as in One vs. the Rest.

c) Influence Estimation

(1) FOR t= 1 TO m

(a) Importance Scoreijt =PNt

n=1 Weight Scorent where:

Weight Scorent =max0,∆(Gini score, primary split)

−∆(Gini score, split with category j moved

to the opposite child node)

(2) Use a paired ttest to compare the Importance Score of covariate i, category j

against a large percentile from the artificial covariate importance scores across

RF .

2. NEXT j

C. NEXT i

130

Results in Table 5.4 show the sensitivity of the ACE feature value selection is

superior. Although the ACE method requires nearly 3 times the CPU processing

time of the rule-based method, both methods are reasonably fast as both require

on average less than one second of CPU time on an Intel Core i5 processor. Nev-

ertheless, the implementation of the ACE method is significantly more complex

than that of the of the rule-based method, which can even be very easily imple-

mented in open-source software such as R (R Development Core Team, 2010;

Hahsler et al., 2010). In addition, there are fewer false alarms associated with

the rule-based method (larger specificity). When both sensitivity and specificity

are simultaneously considered using d, results indicate that the RBFVS method

outperforms ACE in 7 / 11 experiments. The only two experiments where the

ACE method outperforms RBFVS involve those experiments with one or two

categorical (prior to discretization) covariates in a single disease cluster. Exper-

iments 2, 7, and 8 clearly represent the biggest challenge to the RBFVS method

as the method’s sensitivity drops from an average of 0.98 to an average of 0.86.

Even though this work does not intend to be an exhaustive evaluation of the

different parameters for each method, the sensitivity of the rule-based method

can be easily improved if one lowers the confidence threshold for the rules

(confidence(C)≤0.625). Of course, lowering the confidence threshold will

131

Table 5.4. Mean sensitivity and specificity for the RBFVS and ACE methods.

Experiment RBFVS ACE

Sensitivity Specificity dSensitivity Specificity d

0 – 0.98 0.02 – 0.98 0.02

1 1 0.99 0.01 1 0.94 0.06

2 0.8 0.99 0.20 1 0.92 0.08

3 1 0.99 0.01 1 0.9 0.10

4 1 0.99 0.01 1 0.93 0.07

5 1 1 0.00 1 0.97 0.03

6 0.9 0.98 0.10 0.9 0.93 0.12

7 0.8 0.99 0.20 0.9 0.94 0.12

8 0.85 0.99 0.15 0.9 0.97 0.10

9 0.97 0.99 0.03 0.9 0.95 0.11

10 1 0.99 0.01 1 0.96 0.04

increase the number of false alarms (lower specificity). Similarly, the specificity

of the ACE method can be improved by increasing the artificial contrast quantile

(q≥0.9), again, at the expense of sensitivity.

Tables 5.5 - 5.6 break down the information in Table 5.4 by the type of co-

variate. Table 5.5 summarizes the number of categories in a cluster along with

the number of categories detected correctly and in excess. As suggested by the

results in Table 5.4, the ACE method captures a larger number of correctly de-

tected categories in the majority of the experiments. Nevertheless, it also de-

tects a larger number of categories in excess. The number of false alarms with

the rule-based method is between 0.10 and 0.50 categories in excess, while for

the ACE method, there can be up to 2.2 categories in excess. Similar patterns

are also shown in Table 5.6. The ACE tends to capture a larger portion of the

range(s) of the covariates in a cluster, while it also captures a larger range out-

132

Table 5.5. Average number of categories detected correctly and in excess. Experi-

ments with zero categories in cluster indicate numerical clusters.

Experiment Categories

in Cluster

RBFVS ACE

Correctly Detected Detected in Excess Correctly Detected Detected in Excess

0 0 0 0 0 0

1 0 0 0 0 0.40

2 2 1.10 0.10 2 0.50

3 2 1.60 0 2 0.90

4 3 3 0 3 1.40

5 2 1.40 0 2 0.40

6 0 0 0.50 0 1.70

7 2 1 0.30 1.60 2.20

8 2 1.20 0.20 1.90 1.50

9 3 2.90 0.20 2.80 1.80

10 2 1.30 0.20 1.90 1.20

side the cluster. Nevertheless, both methods seem to work very efficiently. On

average, both methods detect at least 50% of the range(s) in the cluster. For

rule-based feature value selection (RBFVS), the largest range detected in excess

is only 0.002, while the largest range detected in excess for the ACE method is

0.026. On the other hand, it is worth mentioning that the higher the order of

the interaction associated with the cluster, the easier it is for both methods to

identify the cluster. In other words, clusters associated with a three covariate

interaction (i.e., Experiment 10) are easier to detect than those associated with

a single covariate in a cluster (i.e., Experiment 2). Clearly, this behavior is in-

tuitive given that clusters associated with higher order interactions define small

regions in space where it is easier to detect an increased number of new test

cases.

133

Table 5.6. Average range detected correctly and in excess. Experiments with zero

range in cluster indicate categorical clusters.

Experiment Range in

Cluster

RBFVS ACE

Correctly Detected Detected in Excess Correctly Detected Detected in Excess

0 0 0 0.002 0 0.004

1 0.15 0.088 0.001 0.125 0.003

2 0 0 0 0 0.008

3 0.15 0.111 0.001 0.121 0.022

4 0.15 0.097 0.001 0.137 0.017

5 0.02 0.02 0 0.02 0.006

6 0.15 0.100 0 0.106 0.026

7 0 0 0 0 0.006

8 0.15 0.082 0 0.072 0

9 0.15 0.102 0 0.103 0.011

10 0.02 0.02 0 0.02 0.014

5.6. Summary

Either the Feature Selection (FS) or the Percentile method can be used to ad-

dress the problem of “global clustering”. They assess whether the process of

interest has deviated from its baseline. In the case of public health surveillance,

“global clustering” relates to a change in the disease incidence rate of interest.

If no “global clustering” structure is detected, one can simply assess the perfor-

mance of the process at a later time. If there is “global clustering”, a follow up

analysis to determine the specific location of the “local clusters” in the data is

of interest. The RBFVS method is a computationally feasible and easy to im-

plement method that can be used to identify the location of these clusters in a

high-dimensional covariate space. First, it discretizes all covariates. Second,

it extracts the important covariate values and/or levels using a fast rule mining

134

classifier. Results in Section 5.5 indicate the method is successful at accurately

identifying the important covariate ranges and/or levels. The method’s sensitiv-

ity in these experiments ranges between 0.8 and 1. In addition, there are very

few false alarms associated with the method as its specificity is always greater

than or equal to 0.98. Clusters strictly involving categorical covariates represent

a challenge to the method, particularly, as the number of covariates in the data

set increases. Lastly, results for the RBFVS suggest the method is a competitor

to the more complex ACE method.

135

CHAPTER 6

CONCLUSIONS AND FUTURE WORK

6.1. Conclusions

There is a vast literature on public health surveillance methods that can only

detect changes in disease incidence rates that occur in geographical regions

(space) and/or time. As the complexities in modern public health data grow,

one needs to be able to develop or improve these methods in order to better

address the current needs of public health surveillance. Nowadays, one is often

faced with extensive databases that record disease incidences as transactions that

may include a large number of covariates, missing values, strong interactions ef-

fects, highly distinct geographical areas, and even disease incidence rates that

are subject to seasonality. This work proposes a supervised learning approach to

public health surveillance that can inherently account for the majority of these

complexities. Although many learners can be used to detect changes in disease

incidence rates, the use of tree-based ensembles is suggested because it is par-

ticularly effective at isolating changes in disease incidence rates that occur even

in a small region of a much larger high-dimensional space. Because tree-based

ensembles are very effective at isolating changes in high-dimensions, they can

also inherently account for seasonality effects and unevenly distributed popu-

lations. Lastly, tree-based ensembles have proven effective in accounting for

136

missing values, non-linear relationships, and even high-order interactions.

In the context of public health surveillance, there are two main problems.

First, there is the global clustering problem that intends to detect whether any

changes have occurred in a disease incidence rate. Two different approaches

to this problem were proposed. First, a signal using feature selection was de-

veloped, where a change in the disease incidence rate is defined as detecting at

least one important covariate using a measure of statistical significance. This

Feature Selection (FS) method is an efficient approach to the problem, although,

its implementation can be complex. Another simpler and easier to implement

method was developed using class probability estimates. Instead of identifying

the important covariates in the data set, the Percentile method identifies the test

incidences that most significantly differ from the baseline incidences and trig-

gers a signal when the number of these unusual incidences exceeds a pre-defined

threshold. The method has highly similar sensitivity to that of the FS method; at

the same time, it can be easily implemented even in open-source software such

as R (R Development Core Team, 2010).

The second main problem in health surveillance is a follow up step to global

clustering. Once a change in the disease incidence rate is detected, one wants to

be able to explore the high-dimensional space where the change is occurring and

137

identify the specific region or set of regions that are contributing to the change.

Rule-Based Feature Value Selection (RBFVS) is proposed as an approach this

problem. When compared to the more complex ACE feature value selection

method, the RBFVS shows smaller sensitivity in roughly 50% of the experi-

ments. Nevertheless, it shows better specificity than the ACE method in more

than 90% of the experiments. RBFVS is a simple methodology that can be eas-

ily explained to a non-technical audience. At the same time, it can also be very

easily implemented.

6.2. Future Work

A natural extension to this work would be to perform a more in-depth study

on the effects of an unevenly distributed population. Chapter 3 considers the

performance of the FS method on the problem where the size of the population

is unevenly distributed across the different geographical regions that conform a

larger study area. There are different approaches than can potentially improve

the sensitivity of tree-based classifiers in this problem. For example, one could

change the manner in which the bootstrap samples used to build the individual

decision trees are drawn. It is possible to perform sampling, such that replace-

ment is considered, but in a manner that the sample is balanced with respect

to the different geographical regions within the study area. Another alternative

138

would adjust the weights of the individual incidences going into the decision

trees. For example, the weight of individual incidences that belong to heavily

populated geographical regions could be decreased.

A similar problem arises when one wants to detect a change in disease inci-

dence rates, but there has been a change in the underlying population during the

study period. For instance, the population for the study area increased from one

year to the next. If the baseline incidences represent the underlying population

for the previous year, the method will detect a change in the disease incidence

rate even when the change being detected is an increase in the size of the popu-

lation. A potential approach to the problem is to adjust the baseline incidences.

One can either increase the number of incidences in the baseline data by sim-

ply creating replicates of several baseline incidences. Otherwise, if the observed

increase in population is known to have happened within a given geographical

region, then, create only replicates of the baseline incidences in that region.

In terms of the local clustering problem, the choice of support and confidence

parameters for the RBFVS method can be application-dependent. Rule pruning

methods can be used to make the method less sensitive to the choice of support

and confidence parameters. Similarly, partial dependence plots can be used as

an alternative approach to feature value selection. Preliminary results indicate

139

partial dependence plots are successful at identifying change points in numeri-

cal covariates. Nevertheless, a lack of sensitivity was observed for categorical

covariates. Further improvements on partial dependent plots for categorical co-

variates are needed.

140

REFERENCES

Alt, F. (1985). Encyclopedia of Statistical Sciences, Volume 6, Chapter Multi-

variate Quality Control. New York: Wiley.

Aradhye, H. B., B. R. Bakshi, R. A. Strauss, and J. F. Davis (2003). Multiscale

SPC using wavelets: Theoretical analysis and properties. AIChE Journal 49(4),

939–958.

Bakshi, B. (1998). Multiscale PCA with application to multivariate statistical

process monitoring. AIChE Journal 44, 1596–1610.

Balakrishnan, N. and M. Koutras (2002). Runs and Scans with Applications.

New York, NY: Wiley.

Benjamini, Y. and Y. Hochberg (1995). Controlling the false discovery rate:

A practical and powerful approach to multiple testing. Journal of the Royal

Statistical Society B 57, 289–300.

Berrado, A. and G. Runger (2009, may). Supervised multivariate discretization

in mixed data with random forests. In IEEE/ACS International Conference on

Computer Systems and Applications, pp. 211–217.

Besag, J. and J. Newell (1991). The detection of clusters in rare diseases.

Journal of the Royal Statistical Society. Series A (Statistics in Society) 154(1),

143–155.

Borgelt, C. (2004). Apriori – Finding Association Rules / Hyperedges with the

Apriori Algorithm. http://fuzzy.cs.uni-magdeburg.de/ borgelt/apriori.html.

Breiman, L. (2001). Random forests. Machine Learning 45(1), 5–32.

Breiman, L., J. Friedman, R. Olshen, and C. Stone (1984). Classification and

Regression Trees. Belmont, MA: Wadsworth.

141

Buckeridge, D. L., H. Burkom, M. Campbell, W. R. Hogan, and A. W. Moore

(2005). Algorithms for rapid outbreak detection: a research synthesis. Journal

of Biomedical Informatics 38(2), 99 – 113.

Crosier, R. B. (1988). Multivariate generalizations of cumulative sum quality-

control schemes. Technometrics 30(3), 291–303.

Dwass, M. (1957). Modified randomization tests for nonparametric hypotheses.

The Annals of Mathematical Statistics 28(1), 181–187.

Glaz, J. and N. Balakrishnan (1999). Scan Statistics and Applications. Boston,

MA: Birkhauser.

Glaz, J., J. Naus, and S. Wallenstein (2001). Scan Statistics. New York:

Springer.

Guyon, I. and A. Elisseeff (2003, Mar). An introduction to variable and feature

selection. Journal of Machine Learning Research 3, 1157–1182.

Hahsler, M., C. Buchta, B. Gruen, and K. Hornik (2010, March). Mining Asso-

ciation Rules and Frequent Itemsets.

Hastie, T., R. Tibshirani, and J. H. Friedman (2001). Chapter Boosting and

Additive Trees, pp. 333–334. Springer Series in Statistics. Springer-Verlag.

Hayter, A. and K. Tsui (1994). Identification and quantification in multivariate

quality control problems. Journal of Quality Technology 26(3), 197–208.

Hotelling, H. (1947). Multivariate Quality Control–illustrated by the Air Test-

ing of Sample Bombsights. New York: McGraw-Hill.

Huntington, R. J. and J. I. Naus (1975). A simpler expression for kth nearest

neighbor coincidence probabilities. The Annals of Probability 3(5), 894–896.

142

Hwang, W., G. Runger, and E. Tuv (2007). Multivariate statistical process

control with artificial contrasts. IIE Transactions: Special Issue on Data Min-

ing 36(6), 659–669.

Jackson, J. (1980). Principal components and factor analysis: Part I–Principal

components. Journal of Quality Technology 12(4), 201–213.

Joner Jr, M. D., W. H. Woodall, and M. R. Reynolds Jr (2008). Detecting a

rate increase using a Bernoulli scan statistic. Statistics in Medicine 27(14),

2555–2575.

Kleinman, K., R. Lazarus, and R. Platt (2004). A Generalized Linear Mixed

Models Approach for Detecting Incident Clusters of Disease in Small Areas,

with an Application to Biological Terrorism. American Journal of Epidemiol-

ogy 159(3), 217–224.

Kulldorff, M. (1997). A spatial scan statistic. Communications in Statistics–

Theory and Methods 26(6), 1481–1496.

Kulldorff, M. (2001). Prospective time periodic geographical disease surveil-

lance using a scan statistic. Journal of the Royal Statistical Society 164(1),

61–72.

Kulldorff, M. and Information Management Services, Inc. (2009). SaTScanTM

v8.0 User Guide.

Kulldorff, M., F. Mostashari, L. Duczmal, W. Yih, K. Kleinman, and R. Platt

(2007). Multivariate scan statistics for disease surveillance. Statistics in

Medicine 26(8), 1824–1833.

Kulldorff, M. and N. Nagarwalla (1995). Spatial disease clusters: Detection

and inference. Statistics in Medicine 14, 799–810.

Lawson, A. and K. Kleinman (2005). Spatial and Syndromic Surveillance.

Hoboken, NJ: Wiley.

143

Limaye, S. S., C. M. Mastrangelo, and D. M. Zerr (2008). A case study in

monitoring hospital-associated infections with count control charts. Quality

Engineering 20(4), 404–413.

Loader, C. R. (1991). Large-deviation approximations to the distribution of

scan statistics. Advances in Applied Probability 23(4), 751–771.

Lowry, C. and D. Montgomery (1995). A review of multivariate control charts.

IIE Transactions 27(6), 800–810.

Lucas, J. M. (1985). Counted data CUSUM’s. Technometrics 27(2), 129–144.

Montgomery, D. C. (2005). Introduction to Statistical Quality Control. USA:

Wiley.

Morton, A. P., M. Whitby, M.-L. McLaws, A. Dobson, S. McElwain, D. Looke,

J. Stackelroth, and A. Sartor (2001). The application of statistical process con-

trol charts to the detection and monitoring of hospital-acquired infections. Jour-

nal of Quality in Clinical Practice 21(4), 112–117.

Naus, J. and S. Wallenstein (2006). Temporal surveillance using scan statistics.

Statistics in Medicine 25(2).

Naus, J. I. (1965). Clustering of random points in two dimensions.

Biometrika 52(1/2), 263–267.

Neff, N. (1978). Piecewise polynomials for the probabilities of clustering on

the unit interval. Ph. D. thesis, Rutgers University, New Brunswick, NJ.

O’Sullivan, D. and D. J. Unwin (2006). Geographic Information Analysis.

Hoboken: Wiley.

Page, E. (1954). Continuous inspection schemes. Biometrics 41(1), 100–115.

144

Percival, D. B. and A. T. Walden (2000). Wavelet methods for time series

analysis. Cambridge, New York: Cambridge University Press.

R Development Core Team (2009). R: A Language and Environment for Sta-

tistical Computing. Vienna, Austria: R Foundation for Statistical Computing.

ISBN 3-900051-07-0.

R Development Core Team (2010). R: A Language and Environment for Sta-

tistical Computing. Vienna, Austria: R Foundation for Statistical Computing.

ISBN 3-900051-07-0.

Renaud, O., J. Starck, and F. Murtagh (2005, December). Wavelet-based com-

bined signal filtering and prediction. IEEE Transactions on Systems, Man, and

Cybernetics, Part B 35(6), 1241–1251.

Reynolds, M. and Z. Stoumbos (1998). The SPRT chart for monitoring a pro-

portion. IIE Transactions 30, 545–561.

Reynolds, M. and Z. Stoumbos (1999). A CUSUM chart for monitoring a

proportion when inspecting continuously. Journal of Quality Technology 31(1),

87–108.

Reynolds, M. R. and Z. G. Stoumbos (2000). A general approach to modeling

CUSUM charts for a proportion. IIE Transactions 32(6), 515–535.

Roberts, S. W. (1959). Control chart tests based on geometric moving averages.

Technometrics 1(3), 239–250.

Runger, G., F. Alt, and D. Montogmery (1996). Contributors to a multivariate

statistical process control signal. Communications in Statistics–Theory and

Methods 25(10), 2203–2213.

Shmueli, G. (2005). Wavelet-based monitoring for modern biosurveillance.

Technical Report RHS-06-002, Robert H. Smith School of Business, College

Park, Maryland.

145

Shmueli, G. and H. Burkom (2007). Statistical chellenges in modern biosurveil-

lance. Technometrics. (In Press).

Sonesson, C. and D. Bock (2003). A review and discussion of prospective sta-

tistical surveillance in public health. Journal of Royal Statistical Society 166(1),

5–21.

Storey, J. (2002). A direct approach to false discovery rates. Journal of the

Royal Statistical Society B. 64(3), 479498.

Stoumbos, Z. G. and J. H. Sullivan (2002). Robustness to non-normality of

the multivariate EWMA control chart. Journal of Quality Technology 34(3),

260–276.

Tan, P., M. Steinbach, and V. Kumar (2006). Introduction to Data Mining,

Chapter Association Analysis: Basic Concepts and Algorithms, pp. 327–314.

Boston, Massachussetts: Addison Wesley.

Tango, T. (1984). The detection of disease clustering in time. Biometrics 40(1),

15–26.

Tango, T. (1995). A class of tests for detecting ‘general ’and ‘focused ’cluster-

ing of rare diseases. Statistics in Medicine 14(21–22), 2323–2334.

Testik, M. C., G. C. Runger, and C. M. Borror (2003). Robustness properties

of multivariate EWMA control charts. Quality and Reliability Engineering

International 19(1), 31–38.

Thacker, S. (2000). Principles and Practice of Public Health Surveillance,

Chapter Historical development, pp. 1–16. New York: Oxford University Press.

Tsui, K.-L., W. Chiu, P. Gierlich, D. Goldsman, X. Liu, and T. Maschek (2008).

A review of healthcare, public health, and syndromic surveillance. Quality

Engineering 20(4), 435–450.

146

Turnbull, B. W., E. J. Iwano, W. S. Burnett, H. L. Howe, and L. C. Clark (1990).

Monitoring for clusters of disease: application to leukemia incidence in upstate

new york. American Journal of Epidemiology 132(supp1), 136–143.

Tuv, E., A. Borisov, G. Runger, and K. Torkkola (2009). Feature selection with

ensembles, artificial variables, and redundancy elimination. J. Mach. Learn.

Res. 10, 1341–1366.

U.S. Census Bureau (2010, August). State & County QuickFacts.

http://quickfacts.census.gov/qfd/states/04/04001.html. Last accessed: October

27, 2010.

U.S. Census Bureau Geography Division (2001, July). Census 2000:

State and State Equivalent Areas in ArcView Shapefile (.shp) format.

http://www.census.gov/geo/www/cob/st2000.html. Last accessed: October 27,

2010.

Wald, A. (1947). Sequential Analysis. New York, NY: Dover Publications, Inc.

Wallenstein, S. (1980). A Test for Detection of Clustering Over Time. Ameri-

can Journal of Epidemiology 111(3), 367–372.

Wallenstein, S. R. and J. I. Naus (1973). Probabilities for a kth nearest neighbor

problem on the line. The Annals of Probability 1(1), 188–190.

Wang, Y. (1995). Jump and sharp cusp detection by wavelets. Biometrika 82(2),

385–397.

Weinstock, M. A. (1981). A generalised scan statistic test for the detection of

clusters. International Journal of Epidemiology 10(3), 289–293.

Woodall, W. (2006). The use of control charts in health-care and public-health

surveillance. Journal of Quality Technology 38(2), 89–134.

147

Woodall, W. H., M. A. Mohammed, J. M. Lucas, and R. Watkins (2006). The

use of control charts in health-care and public-health surveillance. Journal of

Quality Technology 38(2), 89–134.

Woodall, W. H. and M. M. Ncube (1985). Multivariate CUSUM quality-control

procedures. Technometrics 27(3), 285–292.

Xiao, H., H. Sun, B. Ran, and Y. Oh (2005). Fuzzy-neural network traffic

prediction framework with wavelet decomposition. Transportation Research

Record: Journal of Transportation Research 1836(1), 16–36.

148

Public Health Surveillance in High-Dimensions with Supervised Learning

Students also viewed