With an increased focus on incorporating fairness in machine learning models, it becomes imperative not only to assess and mitigate bias at each stage of the machine learning pipeline but also to understand the downstream impacts of bias across stages. Here we consider a general, but realistic, scenario in which a predictive model is learned from (potentially biased) training data, and model predictions are assessed post-hoc for fairness by some auditing method. We provide a theoretical analysis of how a specific form of data bias, differential sampling bias, propagates from the data stage to the prediction stage. Unlike prior work, we evaluate the downstream impacts of data biases quantitatively rather than qualitatively and prove theoretical guarantees for detection. Under reasonable assumptions, we quantify how the amount of bias in the model predictions varies as a function of the amount of differential sampling bias in the data, and at what point this bias becomes provably detectable by the auditor. Through experiments on two criminal justice datasets-- the well-known COMPAS dataset and historical data from NYPD's stop and frisk policy-- we demonstrate that the theoretical results hold in practice even when our assumptions are relaxed.

The opioid overdose epidemic represents a serious public health crisis, with fatality rates rising considerably over the past several years. To help address the abuse of prescription opioids, state governments collect data on dispensed prescriptions, yet the use of these data is typically limited to manual searches. In this paper, we propose a novel graph-based framework for detecting anomalous opioid prescribing patterns in state Prescription Drug Monitoring Program (PDMP) data, which could aid governments in deterring opioid diversion and abuse. Specifically, we seek to identify connected networks of opioid prescribers and dispensers who engage in high-risk and possibly illicit activity. We develop and apply a novel extension of the Non-Parametric Heterogeneous Graph Scan (NPHGS) to two years of de-identified PDMP data from the state of Kansas, and find that NPHGS identifies subgraphs that are significantly more anomalous than those detected by other graph-based methods. NPHGS also reveals clusters of potentially illicit activity, which may strengthen state law enforcement and regulatory capabilities. Our paper is the first to demonstrate how prescription data can systematically identify anomalous opioid prescribers and dispensers, as well as illustrating the efficacy of a network-based approach. Additionally, our technical extensions to NPHGS offer both improved flexibility and graph density reduction, enabling the framework to be replicated across jurisdictions and extended to other problem domains.

Graph neural networks (GNNs) provide a powerful and scalable solution for modeling continuous spatial data. However, they often rely on Euclidean distances to construct the input graphs. This assumption can be improbable in many real-world settings, where the spatial structure is more complex and explicitly non-Euclidean (e.g., road networks). Here, we propose PE-GNN, a new framework that incorporates spatial context and correlation explicitly into the models. Building on recent advances in geospatial auxiliary task learning and semantic spatial embeddings, our proposed method (1) learns a context-aware vector encoding of the geographic coordinates and (2) predicts spatial autocorrelation in the data in parallel with the main task. On spatial interpolation and regression tasks, we show the effectiveness of our approach, improving performance over different state-of-the-art GNN approaches. We observe that our approach not only vastly improves over the GNN baselines, but can match Gaussian processes, the most commonly utilized method for spatial interpolation problems.

Regression discontinuity (RD) designs are widely used to estimate causal effects in the absence of a randomized experiment. However, standard approaches to RD analysis face two significant limitations. First, they require a priori knowledge of discontinuities in treatment. Second, they yield doubly-local treatment effect estimates, and fail to provide more general causal effect estimates away from the discontinuity. To address these limitations, we introduce a novel method for automatically detecting RDs at scale, integrating information from multiple discovered discontinuities with an observational estimator, and extrapolating away from discovered, local RDs. We demonstrate the performance of our method on two synthetic datasets, showing improved performance compared to direct use of an observational estimator, direct extrapolation of RD estimates, and existing methods for combining multiple causal effect estimates. Finally, we apply our novel method to estimate spatially heterogeneous treatment effects in the context of a recent economic development problem.

From ecology to atmospheric sciences, many academic disciplines deal with data characterized by intricate spatio-temporal complexities, the modeling of which often requires specialized approaches. Generative models of these data are of particular interest, as they enable a range of impactful downstream applications like simulation or creating synthetic training data. Recently, COT-GAN, a new GAN algorithm inspired by the theory of causal optimal transport (COT), was proposed in an attempt to improve generation of sequential data. However, the task of learning complex patterns over time and space requires additional knowledge of the specific data structures. In this study, we propose a novel loss objective combined with COT-GAN based on an autoregressive embedding to reinforce the learning of spatio-temporal dynamics. We devise SPATE (spatio-temporal association), a new metric measuring spatio-temporal autocorrelation. We compute SPATE for real and synthetic data samples and use it to compute an embedding loss that considers space-time interactions, nudging the GAN to learn outputs that are faithful to the observed dynamics. We test our new SPATE-GAN on a diverse set of spatio-temporal patterns: turbulent flows, log-Gaussian Cox processes and global weather data. We show that our novel embedding loss improves performance without any changes to the architecture of the GAN backbone, highlighting our model's increased capacity for capturing autoregressive structures.

We propose a new approach, the calibrated nonparametric scan statistic (CNSS), for more accurate detection of anomalous patterns in large-scale, real-world graphs. Scan statistics identify connected subgraphs that are interesting or unexpected through maximization of a likelihood ratio statistic; in particular, nonparametric scan statistics (NPSSs) identify subgraphs with a higher than expected proportion of individually significant nodes. However, we show that recently proposed NPSS methods are miscalibrated, failing to account for the maximization of the statistic over the multiplicity of subgraphs. This results in both reduced detection power for subtle signals, and low precision of the detected subgraph even for stronger signals. Thus we develop a new statistical approach to recalibrate NPSSs, correctly adjusting for multiple hypothesis testing and taking the underlying graph structure into account. While the recalibration, based on randomization testing, is computationally expensive, we propose both an efficient (approximate) algorithm and new, closed-form lower bounds (on the expected maximum proportion of significant nodes for subgraphs of a given size, under the null hypothesis of no anomalous patterns). These advances, along with the integration of recent core-tree decomposition methods, enable CNSS to scale to large real-world graphs, with substantial improvement in the accuracy of detected subgraphs. Extensive experiments on both semi-synthetic and real-world datasets are demonstrated to validate the effectiveness of our proposed methods, in comparison with state-of-the-art counterparts.

Machine learning is gaining popularity in a broad range of areas working with geographic data. Here, data often exhibit spatial effects, which can be difficult to learn for neural networks. We propose SXL, a method for embedding information on the autoregressive nature of spatial data directly into the learning process using auxiliary tasks. We utilize the local Moran's I, a measure of local spatial autocorrelation, to "nudge" the model to learn the direction and magnitude of local spatial effects, complementing learning of the primary task. We further introduce a novel expansion of Moran's I to multiple resolutions, capturing spatial interactions over longer and shorter distances simultaneously. The novel multi-resolution Moran's I can be constructed easily and offers seamless integration into existing machine learning frameworks. Over a range of experiments using real-world data, we highlight how our method consistently improves the training of neural networks in unsupervised and supervised learning tasks. In generative spatial modeling experiments, we propose a novel loss for auxiliary task GANs utilizing task uncertainty weights. SXL outperforms domain-specific spatial interpolation benchmarks, highlighting its potential for downstream applications.

Standard spatial cluster detection methods used in public health surveillance assign each disease case to a single location (typically, the patient's home address), aggregate locations to small areas, and monitor the number of cases in each area over time. However, such methods cannot detect clusters of disease resulting from visits to non-residential locations, such as a park or a university campus. Thus we develop two new spatial scan methods, the unconditional and conditional spatial logistic models, to search for spatial clusters of increased infection risk. We use mobility data from two sets of individuals, disease cases and healthy individuals, where each individual is represented by a sparse sample of geographical locations (e.g., from geo-tagged social media data). The methods account for the multiple, varying number of spatial locations observed per individual, either by non-parametric estimation of the odds of being a case, or by matching case and control individuals with similar numbers of observed locations. Applying our methods to synthetic and real-world scenarios, we demonstrate robust performance on detecting spatial clusters of infection risk from mobility data, outperforming competing baselines.

Identifying changes in model parameters is fundamental in machine learning and statistics. However, standard changepoint models are limited in expressiveness, often addressing unidimensional problems and assuming instantaneous changes. We introduce change surfaces as a multidimensional and highly expressive generalization of changepoints. We provide a model-agnostic formalization of change surfaces, illustrating how they can provide variable, heterogeneous, and non-monotonic rates of change across multiple dimensions. Additionally, we show how change surfaces can be used for counterfactual prediction. As a concrete instantiation of the change surface framework, we develop Gaussian Process Change Surfaces (GPCS). We demonstrate counterfactual prediction with Bayesian posterior mean and credible sets, as well as massive scalability by introducing novel methods for additive non-separable kernels. Using two large spatio-temporal datasets we employ GPCS to discover and characterize complex changes that can provide scientific and policy relevant insights. Specifically, we analyze twentieth century measles incidence across the United States and discover previously unknown heterogeneous changes after the introduction of the measles vaccine. Additionally, we apply the model to requests for lead testing kits in New York City, discovering distinct spatial and demographic patterns.

Identifying anomalous patterns in real-world data is essential for understanding where, when, and how systems deviate from their expected dynamics. Yet methods that separately consider the anomalousness of each individual data point have low detection power for subtle, emerging irregularities. Additionally, recent detection techniques based on subset scanning make strong independence assumptions and suffer degraded performance in correlated data. We introduce methods for identifying anomalous patterns in non-iid data by combining Gaussian processes with novel log-likelihood ratio statistic and subset scanning techniques. Our approaches are powerful, interpretable, and can integrate information across multiple streams. We illustrate their performance on numeric simulations and three open source spatiotemporal datasets of opioid overdose deaths, 311 calls, and storm reports.

Inferring causal relationships in observational data is crucial for understanding scientific and social processes. We develop the first statistical machine learning approach for automatically discovering regression discontinuity designs (RDDs), a quasi-experimental setup often used in econometrics. Our method identifies interpretable, localized RDDs in arbitrary dimensional data and can seamlessly compute treatment effects without expert supervision. By applying the technique to a variety of synthetic and real datasets, we demonstrate robust performance under adverse conditions including unobserved variables, substantial noise, and model

We present a scalable Gaussian process model for identifying and characterizing smooth multidimensional changepoints, and automatically learning changes in expressive covariance structure. We use Random Kitchen Sink features to flexibly define a change surface in combination with expressive spectral mixture kernels to capture the complex statistical structure. Finally, through the use of novel methods for additive non-separable kernels, we can scale the model to large datasets. We demonstrate the model on numerical and real world data, including a large spatio-temporal disease dataset where we identify previously unknown heterogeneous changes in space and time.

In applied fields, practitioners hoping to apply causal structure learning or causal orientation algorithms face an important question: which independence test is appropriate for my data? In the case of real-valued iid data, linear dependencies, and Gaussian error terms, partial correlation is sufficient. But once any of these assumptions is modified, the situation becomes more complex. Kernel-based tests of independence have gained popularity to deal with nonlinear dependencies in recent years, but testing for conditional independence remains a challenging problem. We highlight the important issue of non-iid observations: when data are observed in space, time, or on a network, “nearby” observations are likely to be similar. This fact biases estimates of dependence between variables. Inspired by the success of Gaussian process regression for handling non-iid observations in a wide variety of areas and by the usefulness of the Hilbert-Schmidt Independence Criterion (HSIC), a kernel-based independence test, we propose a simple framework to address all of these issues: first, use Gaussian process regression to control for certain variables and to obtain residuals. Second, use HSIC to test for independence. We illustrate this on two classic datasets, one spatial, the other temporal, that are usually treated as iid. We show how properly accounting for spatial and temporal variation can lead to more reasonable causal graphs. We also show how highly structured data, like images and text, can be used in a causal inference framework using a novel structured input/output Gaussian process formulation. We demonstrate this idea on a dataset of translated sentences, trying to predict the source language.

Gaussian processes (GPs) are a flexible class of methods with state of the art performance on spatial statistics applications. However, GPs require O(n^3) computations and O(n^2) storage, and popular GP kernels are typically limited to smoothing and interpolation. To address these difficulties, Kronecker methods have been used to exploit structure in the GP covariance matrix for scalability, while allowing for expressive kernel learning (Wilson et al., 2014). However, fast Kronecker methods have been confined to Gaussian likelihoods. We propose new scalable Kronecker methods for Gaussian processes with non-Gaussian likelihoods, using a Laplace approximation which involves linear conjugate gradients for inference, and a lower bound on the GP marginal likelihood for kernel learning. Our approach has near linear scaling, requiring O(D n^(D+1)/D) operations and O(D n^2/D) storage, for n training data-points on a dense D > 1 dimensional grid. Moreover, we introduce a log Gaussian Cox process, with highly expressive kernels, for modelling spatiotemporal count processes, and apply it to a point pattern (n = 233,088) of a decade of crime events in Chicago. Using our model, we discover spatially varying multiscale seasonal trends and produce highly accurate long-range local area forecasts.

Recent research has highlighted the need for upstream planning in healthcare service delivery systems, patient scheduling, and resource allocation in the hospital inpatient setting. This study examines the value of upstream planning within hospital-wide resource allocation decisions based on machine learning (ML) and mixed-integer programming (MIP), focusing on prediction of diagnosis-related groups (DRGs) and the use of these predictions for allocating scarce hospital resources. DRGs are a payment scheme employed at patients’ discharge, where the DRG and length of stay determine the revenue that the hospital obtains. We show that early and accurate DRG classification using ML methods, incorporated into an MIP-based resource allocation model, can increase the hospital’s contribution margin, the number of admitted patients, and the utilization of resources such as operating rooms and beds. We test these methods on hospital data containing more than 16,000 inpatient records and demonstrate improved DRG classification accuracy as compared to the hospital’s current approach. The largest improvements were observed at and before admission, when information such as procedures and diagnoses is typically incomplete, but performance was improved even after a substantial portion of the patient’s length of stay, and under multiple scenarios making different assumptions about the available information. Using the improved DRG predictions within our resource allocation model improves contribution margin by 2.9% and the utilization of scarce resources such as operating rooms and beds from 66.3% to 67.3% and from 70.7% to 71.7%, respectively. This enables 9.0% more nonurgent elective patients to be admitted as compared to the baseline.

Event detection in social media is an important but challenging problem. Most existing approaches are based on burst detection, topic modeling, or clustering techniques, which cannot naturally model the implicit heterogeneous network structure in social media. As a result, only limited information, such as terms and geographic locations, can be used. This paper presents Non-Parametric Heterogeneous Graph Scan (NPHGS), a new approach that considers the entire heterogeneous network for event detection: we first model the network as a "sensor" network, in which each node senses its "neighborhood environment" and reports an empirical p-value measuring its current level of anomalousness for each time interval (e.g., hour or day). Then, we efficiently maximize a nonparametric scan statistic over connected subgraphs to identify the most anomalous network clusters. Finally, the event represented by each cluster is summarized with information such as type of event, geographical locations, time, and participants. As a case study, we consider two applications using Twitter data, civil unrest event detection and rare disease outbreak detection, and present empirical evaluations illustrating the effectiveness and efficiency of our proposed approach.

We explore scalable and accurate dynamic pattern detection methods in graph-based data sets. We apply our proposed Dynamic Subset Scan method to the task of detecting, tracking, and source-tracing contaminant plumes spreading through a water distribution system equipped with noisy, binary sensors. While static patterns affect the same subset of data over a period of time, dynamic patterns may affect different subsets of the data at each time step. These dynamic patterns require a new approach to define and optimize penalized likelihood ratio statistics in the subset scan framework, as well as new computational techniques that scale to large, real-world networks. To address the first concern, we develop new subset scan methods that allow the detected subset of nodes to change over time, while incorporating temporal consistency constraints to reward patterns that do not dramatically change between adjacent time steps. Second, our Additive Graph Scan algorithm allows our novel scan statistic to process small graphs (500 nodes) in 4.1 seconds on average while maintaining an approximation ratio over 99% compared to an exact optimization method, and to scale to large graphs with over 12,000 nodes in 30 minutes on average. Evaluation results across multiple detection, tracking, and source-tracing tasks demonstrate substantial performance gains achieved by the Dynamic Subset Scan approach.

We present Generalized Fast Subset Sums (GFSS), a new Bayesian framework for scalable and accurate detection of irregularly shaped spatial clusters using multiple data streams. GFSS extends the previously proposed Multivariate Bayesian Scan Statistic (MBSS) and Fast Subset Sums (FSS) approaches for detection of emerging events. The detection power of MBSS is primarily limited by computational considerations, which limit it to searching over circular spatial regions. GFSS enables more accurate and timely detection by defining a hierarchical prior over all subsets of the N locations, first selecting a local neighborhood consisting of a center location and its neighbors, and introducing a sparsity parameter p to describe how likely each location in the neighborhood is to be affected. This approach allows us to consider all possible subsets of locations (including irregularly-shaped regions) but also puts higher weight on more compact regions. We demonstrate that MBSS and FSS are both special cases of this general framework (assuming p = 1 and p = 0.5 respectively), but substantially higher detection power can be achieved by choosing an appropriate value of p. Thus we show that the distribution of the sparsity parameter p can be accurately learned from a small number of labeled events. Our evaluation results (on synthetic disease outbreaks injected into real-world hospital data) show that the GFSS method with learned sparsity parameter has higher detection power and spatial accuracy than MBSS and FSS, particularly when the affected region is irregular or elongated. We also show that the learned models can be used for event characterization, accurately distinguishing between two otherwise identical event types based on the sparsity of the affected spatial region.

OBJECTIVE: Evidence suggests that the medication lists of patients are often incomplete and could negatively affect patient outcomes. In this article, the authors propose the application of collaborative filtering methods to the medication reconciliation task. Given a current medication list for a patient, the authors employ collaborative filtering approaches to predict drugs the patient could be taking but are missing from their observed list. DESIGN: The collaborative filtering approach presented in this paper emerges from the insight that an omission in a medication list is analogous to an item a consumer might purchase from a product list. Online retailers use collaborative filtering to recommend relevant products using retrospective purchase data. In this article, the authors argue that patient information in electronic medical records, combined with artificial intelligence methods, can enhance medication reconciliation. The authors formulate the detection of omissions in medication lists as a collaborative filtering problem. Detection of omissions is accomplished using several machine-learning approaches. The effectiveness of these approaches is evaluated using medication data from three long-term care centers. The authors also propose several decision-theoretic extensions to the methodology for incorporating medical knowledge into recommendations. RESULTS: Results show that collaborative filtering identifies the missing drug in the top-10 list about 40-50% of the time and the therapeutic class of the missing drug 50%-65% of the time at the three clinics in this study. CONCLUSION: Results suggest that collaborative filtering can be a valuable tool for reconciling medication lists, complementing currently recommended process-driven approaches. However, a one-size-fits-all approach is not optimal, and consideration should be given to context (eg, types of patients and drug regimens) and consequence (eg, the impact of omission on outcomes).

T-Cube Web Interface is a generic tool to visualize and manipulate large scale multivariate time series datasets. The interface allows the user to execute complex queries quickly and to run various types of statistical tests on the loaded data. We show its utility in an important application scenario: real-time bio-surveillance system designed to support rapid detection and mitigation of bio-medical threats in developing countries.

We propose a new method for detecting patterns of anomalies in categorical datasets. We assume that anomalies are generated by some underlying process which affects only a particular subset of the data. Our method consists of two steps: we first use a "local anomaly detector" to identify individual records with anomalous attribute values, and then detect patterns where the number of anomalous records is higher than expected. Given the set of anomalies flagged by the local anomaly detector, we search over all subsets of the data defined by any set of fixed values of a subset of the attributes, in order to detect self-similar patterns of anomalies. We wish to detect any such subset of the test data which displays a significant increase in anomalous activity as compared to the normal behavior of the system (as indicated by the training data). We perform significance testing to determine if the number of anomalies in any subset of the test data is significantly higher than expected, and propose an efficient algorithm to perform this test over all such subsets of the data. We show that this algorithm is able to accurately detect anomalous patterns in real-world hospital, container shipping and network intrusion data.

We propose a new class of spatio-temporal cluster detection methods designed for the rapid detection of emerging space-time clusters. We focus on the motivating application of prospective disease surveillance: detecting space-time clusters of disease cases resulting from an emerging disease outbreak. Automatic, real-time detection of outbreaks can enable rapid epidemiological response, potentially reducing rates of morbidity and mortality. Building on the prior work on spatial and space-time scan statistics, our methods combine time series analysis (to determine how many cases we expect to observe for a given spatial region in a given time interval) with new "emerging cluster" space-time scan statistics (to decide whether an observed increase in cases in a region is significant), enabling fast and accurate detection of emerging outbreaks. We evaluate these methods on two types of simulated outbreaks: aerosol release of inhalational anthrax (e.g. from a bioterrorist attack) and FLOO ("Fictional Linear Onset Outbreak"), injected into actual baseline data (Emergency Department records and over-the-counter drug sales data from Allegheny County). We demonstrate that our methods are successful in rapidly detecting both outbreak types while keeping the number of false positives low, and show that our new "emerging cluster" scan statistics consistently outperform the standard "persistent cluster" scan statistics approach.

We propose a new Bayesian method for spatial cluster detection, the “Bayesian spatial scan statistic,” and compare this method to the standard (frequentist) scan statistic approach. We demonstrate that the Bayesian statistic has several advantages over the frequentist approach, including increased power to detect clusters and (since randomization testing is unnecessary) much faster runtime. We evaluate the Bayesian and fre- quentist methods on the task of prospective disease surveillance: detect- ing spatial clusters of disease cases resulting from emerging disease out- breaks. We demonstrate that our Bayesian methods are successful in rapidly detecting outbreaks while keeping number of false positives low.

Given an N x N grid of squares, where each square has a count cij and an underlying population pij, our goal is to find the rectangular region with the highest density, and to calculate its significance by randomization. An arbitrary density function D, dependent on a region's total count C and total population P, can be used. For example, if each count represents the number of disease cases occurring in that square, we can use Kulldorff's spatial scan statistic DK to find the most significant spatial disease cluster. A naive approach to finding the maximum density region requires O(N4) time, and is generally computationally infeasible. We present a multiresolution algorithm which partitions the grid into overlapping regions using a novel overlap-kd tree data structure, bounds the maximum score of subregions contained in each region, and prunes regions which cannot contain the maximum density region. For sufficiently dense regions, this method finds the maximum density region in O((N log N)2) time, in practice resulting in significant (20-2000x) speedups on both real and simulated datasets.

Assume a uniform, multidimensional grid of bivariate data, where each cell of the grid has a count ci and a baseline bi. Our goal is to find spatial regions (d-dimensional rectangles) where the ci are significantly higher than expected given bi. We focus on two applications: detection of clusters of disease cases from epidemiological data (emergency depart- ment visits, over-the-counter drug sales), and discovery of regions of in- creased brain activity corresponding to given cognitive tasks (from fMRI data). Each of these problems can be solved using a spatial scan statistic (Kulldorff, 1997), where we compute the maximum of a likelihood ratio statistic over all spatial regions, and find the significance of this region by randomization. However, computing the scan statistic for all spatial regions is generally computationally infeasible, so we introduce a novel fast spatial scan algorithm, generalizing the 2D scan algorithm of (Neill and Moore, 2004) to arbitrary dimensions. Our new multidimensional multiresolution algorithm allows us to find spatial clusters up to 1400x faster than the naive spatial scan, without any loss of accuracy. 1 Introduction One of the core goals of modern statistical inference and data mining is to discover patterns and relationships in data. In many applications, however, it is important not only to discover patterns, but to distinguish those patterns that are significant from those that are likely to have occurred by chance. This is particularly important in epidemiological applications, where a rise in the number of disease cases in a region may or may not be indicative of an emerging epidemic. In order to decide whether further investigation is necessary, epidemiologists must know not only the location of a possible outbreak, but also some measure of the likelihood that an outbreak is occurring in that region. Similarly, when investigating brain imaging data, we want to not only find regions of increased activity, but determine whether these increases are significant or due to chance fluctuations. More generally, we are interested in spatial data mining problems where the goal is detec- tion of overdensities: spatial regions with high counts relative to some underlying baseline. In the epidemiological datasets, the count is some quantity (e.g. number of disease cases, or units of cough medication sold) in a given area, where the baseline is the expected value of that quantity based on historical data. In the brain imaging datasets, our count is the total fMRI activation in a given set of voxels under the experimental condition, while our baseline is the total activation in that set of voxels under the null or control condition. We consider the case in which data has been aggregated to a uniform, d-dimensional grid. For the fMRI data, we have three spatial dimensions; for the epidemiological data, we have two spatial dimensions but also use several other quantities (time, patients' age and gender) as "pseudo-spatial" dimensions; this is discussed in more detail below. In the general case, let G be a d-dimensional grid of cells, with size N1 N2 ... Nd. Each cell si G (where i is a d-dimensional vector) is associated with a count ci and a baseline bi. Our goal is to search over all d-dimensional rectangular regions S G, and find regions where the total count C(S) = S ci is higher than expected, given the baseline B(S) = S bi. In addition to discovering these high-density regions, we must also perform statistical testing to determine whether these regions are significant. As is necessary in the scan statistics framework, we focus on finding the single, most significant region; the method can be iterated (removing each significant cluster once it is found) to find multiple significant regions. 1.1 Likelihood ratio statistics Our basic model assumes that counts ci are generated by an inhomogeneous Poisson pro- cess with mean qbi, where q (the underlying ratio of count to baseline) may vary spatially. We wish to detect hyper-rectangular regions S such that q is significantly higher inside S than outside S. To do so, for a given region S, we assume that q = qin uniformly for cells si S, and q = qout uniformly for cells si G-S. We then test the null hypothesis H0(S): qin (1+)qout against the alternative hypothesis H1(S): qin > (1+)qout. If = 0, this is equivalent to the classical spatial scan statistic [1-2]: we are testing for regions where qin is greater than qout . However, in many real-world applications (including the epidemiological and fMRI datasets discussed later) we expect some fluctuation in the underlying baseline; thus, we do not want to detect all deviations from baseline, but only those where the amount of deviation is greater than some threshold. For example, a 10% increase in disease cases in some region may not be interesting to epidemiologists, even if the underlying population is large enough to conclude that this is a "real" (statistically significant) increase in q. By increasing , we can focus the scan statistic on regions with larger ratios of count to base- line. For example, we can use the scan statistic with = 0.25 to test for regions where qin is more than 25% higher than qout . Following Kulldorff [1], our spatial scan statistic is the maximum, over all regions S, of the ratio of the likelihoods under the alternative and null hypotheses. Taking logs for convenience, we have: sup q s D i (S) = log in>(1+)qout S P(ci Po(qinbi))siG-S P(ci Po(qoutbi)) sup qin(1+)qout siS P(ci Po(qinbi))siG-S P(ci Po(qoutbi)) C(S) C C = ( tot tot sgn) C(S) log + (C -C(S) ( tot 1 + )B(S) -C(S))log Btot -B(S) -CtotlogBtot+B(S) where C(S) and B(S) are the count and baseline of the region S under consideration, Ctot and Btot are the total count and baseline of the entire grid G, and sgn = +1 if C(S) > (1 + B(S) )Ctot-C(S) and -1 otherwise. Then the scan statistic D B ,max is equal to the maximum D(S) tot -B(S) over all spatial regions (d-dimensional rectangles) under consideration. We note that our statistical and computational methods are not limited to the Poisson model given here; any model of null and alternative hypotheses such that the resulting statistic D(S) satisfies the conditions given in [4] can be used for the fast spatial scan. 1.2 Randomization testing Once we have found the highest scoring region S = arg maxS D(S) of grid G, we must still determine the statistical significance of this region. Since the exact distribution of the test statistic Dmax is only known in special cases, in general we must find the region's p-value by randomization. To do so, we run a large number R of random replications, where a replica has the same underlying baselines bi as G, but counts are randomly drawn from the null hypothesis H0(S). More precisely, we pick ci Po(qbi), where q = qin = (1+) Ctot Btot +B(S) for si S, and q = qout = Ctot for s B i tot +B(S) G - S. The number of replicas G with Dmax(G ) Dmax(G), divided by the total number of replications R, gives us the p-value for our most significant region S. If this p-value is less than (where is the false positive rate, typically chosen to be 0.05 or 0.1), we can conclude that the discovered region is statistically significant at level . 1.3 The naive spatial scan The simplest method of finding Dmax is to compute D(S) for all rectangular regions of sizes k1 k2 ...kd, where 1 kj Nj. Since there are a total of d (N j=1 j - kj + 1) regions of each size, there are a total of O(d N2) j=1 regions to examine. We can compute D(S) j for any region S in constant time, by first finding the count C(S) and baseline B(S), then computing D.1 This allows us to compute Dmax of a grid G in O(d N2) j=1 time. However, j significance testing by randomization also requires us to find Dmax for each replica G , and compare this to Dmax(G); thus the total complexity is multiplied by the number of replications R. When the size of the grid is large, as is the case for the epidemiological and fMRI datasets we are considering, this naive approach is computationally infeasible. Instead, we apply our "overlap-multiresolution partitioning" algorithm [3-4], generalizing this method from two-dimensional to d-dimensional datasets. This reduces the complexity to O(d N j=1 j log N j ) in cases where the most significant region S has a sufficiently high ra- tio of count to baseline, and (as we show in Section 3) typically results in tens to thousands of times speedup over the naive approach. We note that this fast spatial scan algorithm is exact (always finds the correct value of Dmax and the corresponding region S); the speedup results from the observation that we do not need to search a given set of regions if we can prove that none of them have score > Dmax. Thus we use a top-down, branch-and-bound approach: we maintain the current maximum score of the regions we have searched so far, calculate upper bounds on the scores of subregions contained in a given region, and prune regions whose upper bounds are less than the current value of Dmax. When searching a replica grid, we care only whether Dmax of the replica grid is greater than Dmax(G). Thus we can use Dmax of the original grid for pruning on the replicas, and can stop searching a replica if we find a region with score > Dmax(G). 2 Overlap-multiresolution partitioning As in [4], we use a multiresolution search method which relies on an overlap-kd tree data structure. The overlap-kd tree, like kd-trees [5] and quadtrees [6], is a hierarchical, space- partitioning data structure. The root node of the tree represents the entire space under consideration (i.e. the entire grid G), and each other node represents a subregion of the grid. Each non-leaf node of a d-dimensional overlap-kd tree has 2d children, an "upper" and a "lower" child in each dimension. For example, in three dimensions, a node has six children: upper and lower children in the x, y, and z dimensions. The overlap-kd tree is different from the standard kd-tree and quadtree in that adjacent regions overlap: rather than splitting the region in half along each dimension, instead each child contains more than half the area of the parent region. For example, a 64 64 64 grid will have six children: two of size 48 6464, two of size 644864, and two of size 646448. 1An old trick makes it possible to compute the count and baseline of any rectangular region in time constant in N: we first form a d-dimensional array of the cumulative counts, then compute each region's count by adding/subtracting at most 2d cumulative counts. Note that because of the exponential dependence on d, these techniques suffer from the "curse of dimensionality": neither the naive spatial scan, nor the fast spatial scan discussed below, are appropriate for very high dimensional datasets. In general, let region S have size k1 k2...kd. Then the two children of S in dimension j (for j = 1 . . . d) have size k1 ...kj-1 fjkj kj+1 ...kd, where 1 < f 2 j < 1. This partitioning (for the two-dimensional case, where f1 = f2 = 3 ) is illustrated in Figure 1. 4 Note that there is a region SC common to all of these children; we call this region the center of S. When we partition region S in this manner, it can be proved that any subregion of S either a) is contained entirely in (at least) one of S1 . . . S2d, or b) contains the center region SC. Figure 1 illustrates each of these possibilities, for the simple case of d = 2. S Figure 1: Overlap-multires partitioning of region S (for d = 2). Any subregion of S either a) is contained in some S S_1 S_2 S_3 i, S_4 S_C i = 1 . . . 4, or b) contains SC. Now we can search all subregions of S by recursively searching S1 . . . S2d, then searching all of the regions contained in S which contain the center SC. There may be a large number of such "outer regions," but since we know that each such region contains the center, we can place very tight bounds on the score of these regions, often allowing us to prune most or all of them. Thus the basic outline of our search procedure (ignoring pruning, for the moment) is: overlap-search(S) { call base-case-search(S) define child regions S1..S2d, center SC as above call overlap-search(Si) for i=1..2d for all S' such that S' is contained in S and contains S_C, call base-case-search(S') } The fractions fi are selected based on the current sizes ki of the region being searched: if ki = 2m, then fi = 3 , and if k . For simplicity, we assume that 4 i = 3 2m, then fi = 23 all Ni are powers of two, and thus all region sizes ki will fall into one of these two cases. Repeating this partitioning recursively, we obtain the overlap-kd tree structure. For d = 2, the first two levels of the overlap-kd tree are shown in Figure 2. Figure 2: The first two levels of the two- dimensional overlap-kd tree. Each node represents a gridded region (denoted by a thick rectangle) of the entire dataset (thin square and dots). The overlap-kd tree has several useful properties, which we present here without proof. First, for every rectangular region S G, either S is a gridded region (contained in the overlap-kd tree), or there exists a unique gridded region S such that S is an outer region of S (i.e. S is contained in S , and contains the center region of S ). This means that, if overlap-search is called exactly once for each gridded region2, and no pruning is done, then base-case-search will be called exactly once for every rectangular region S G. In practice, we will prune many regions, so base-case-search will be called at most once for every rect- angular region, and every region will be either searched or pruned. The second nice prop- erty of our overlap-kd tree is that the total number of gridded regions is O(d N j=1 j log N j ). This implies that, if we are able to prune (almost) all outer regions, we can find Dmax of the grid in O(d N N2) j=1 j log N j ) time rather than O(dj=1 . In fact, we may not even need to j search all gridded regions, so in many cases the search will be even faster. 2As in [4], we use "lazy expansion" to ensure that gridded regions are not multiply searched. 2.1 Score bounds and pruning We now consider which regions can be pruned (discarded without searching) during our multiresolution search procedure. First, given some region S, we must calculate an upper bound on the scores D(S ) for regions S S. More precisely, we are interested in two upper bounds: a bound on the score of all subregions S S, and a bound on the score of the outer subregions of S (those regions contained in S and containing its center SC). If the first bound is less than or equal to Dmax, we can prune region S completely; we do not need to search any (gridded or outer) subregion of S. If only the second bound is less than or equal to Dmax, we do not need to search the outer subregions of S, but we must recursively call overlap-search on the gridded children of S. If both bounds are greater than Dmax, we must both recursively call overlap-search and search the outer regions. Score bounds are calculated based on various pieces of information about the subregions of S, including: upper and lower bounds bmax, bmin on the baseline of subregions S ; an upper bound dmax on the ratio C of S ; an upper bound d of S B inc on the ratio C B -SC; and a lower bound dmin on the ratio C of S B -S . We also know the count C and baseline B of region S, and the count ccenter and baseline bcenter of region SC. Let cin and bin be the count and baseline of S . To find an upper bound on D(S ), we must calculate the values of cin and bin which maximize D subject to the given constraints: cin-ccenter bin-bcenter dinc, cin bin dmax, C-cin B-bin dmin, and bmin bin bmax. The solution to this maximization problem is derived in [4], and (since scores are based only on count and baseline rather than the size and shape of the region) it applies directly to the multidimensional case. The bounds on baselines and ratios C are first calculated using global values (as a fast, "first-pass" pruning technique). B For the remaining, unpruned regions, we calculate tighter bounds using the quartering method of [4], and use these to prune more regions. 2.2 Related work Our work builds most directly on the results of Kulldorff [1], who presents the two- dimensional spatial scan framework and the classical ( = 0) likelihood ratio statistic. It also extends [4], in which we present the two-dimensional fast spatial scan. Our major extensions in the present work are twofold: the d-dimensional fast spatial scan, and the generalized likelihood ratio statistics D. A variety of other cluster detection techniques exist in the literature on epidemiology [1-3, 7-8], brain imaging [9-11], and machine learn- ing [12-15]. The machine learning literature focuses on heuristic or approximate cluster- finding techniques, which typically cannot deal with spatially varying baselines, and more importantly, give no information about the statistical significance of the clusters found. Our technique is exact (in that it calculates the maximum of the likelihood ratio statistic over all hyper-rectangular spatial regions), and uses a powerful statistical test to determine significance. Nevertheless, other methods in the literature have some advantages over the present approach, such as applicability to high-dimensional data and fewer assumptions on the underlying model. The fMRI literature generally tests significance on a per-voxel basis (after applying some method of spatial smoothing); clusters must then be inferred by grouping individually significant voxels, and (with the exception of [10]) no per-cluster false positive rate is guaranteed. The epidemiological literature focuses on detecting signif- icant circular, two-dimensional clusters, and thus cannot deal with multidimensional data or elongated regions. Detection of elongated regions is extremely important in both epi- demiology (because of the need to detect windborne or waterborne pathogens) and brain imaging (because of the "folded sheet" structure of the brain); the present work, as well as [4], allow detection of such clusters.

Given an N(cid:2)N grid of squares, where each square has a count and an un- derlying population, our goal is to find the square region with the highest density, and to calculate its significance by randomization. Any density measure D, dependent on the total count and total population of a re- gion, can be used. For example, if each count represents the number of disease cases occurring in that square, we can use Kulldorff’s spatial scan statistic DK to find the most significant spatial disease cluster. A naive approach to finding the maximum density region requires O(N 3) time, and is generally computationally infeasible. We present a novel algorithm which partitions the grid into overlapping regions, bounds the maximum score of subregions contained in each region, and prunes re- gions which cannot contain the maximum density region. For sufficiently dense regions, this method finds the maximum density region in optimal O(N2) time, in practice resulting in significant (10-200x) speedups.

In many real-world situations, "cooperation" in the simple sense of the Prisoner's Dilemma is not sufficient for success: instead, cooperators must precisely coordinate more complex behaviors in a noisy environment. We investigate one such model, the Turn-Taking Dilemma, a variant of the repeated Prisoner's Dilemma in which the highest total payoff is achieved not by simultaneous mutual cooperation, but by taking turns defecting (alternating temptation and sucker payoffs). The Turn-Taking Dilemma more accurately models interactions where players must take short-term losses for long-term gains: situations marked by the intricate give-and-take of bargaining and compromise. Using "evolutionary dominance" as a general measure of performance, we investigated which strategies are most successful in Turn-Taking Dilemma interactions. Our experiments demonstrate that turn-taking can be achieved in a noisy environment, even when agents have strict resource constraints (limited memory strategies). Top strategies such as EXALT2 can effectively coordinate turn-taking under noise, while exploiting cooperators and resisting exploitation by defectors; these strategies are likely to achieve success in the variety of real-world interactions modeled by the Turn-Taking Dilemma.