Change Surfaces for Expressive Multidimensional
Changepoints and Counterfactual Prediction

In this paper, we introduce change surfaces as expressive, multidimensional generalizations of changepoints. We present a model-agnostic formulation of change surfaces and instantiate this framework with scalable Gaussian process models. The resulting model is capable of automatically learning expressive covariance functions and a sophisticated continuous change surface. Additionally, we derive massively scalable inference procedures, as well as counterfactual prediction techniques. Finally, we apply the proposed methods to a wide variety of numerical data and complex human systems. In particular, we:

1. 1.

   Introduce change surfaces as multidimensional and highly flexible generalizations of changepoint modeling.

2. 2.

   Introduce a procedure which allows one to specify background functions and change functions, for more powerful inductive biases and added interpretability.

3. 3.

   Provide a new framework for counterfactual prediction using change surfaces.

4. 4.

   Present the Gaussian Process Change Surface model (GPCS) which models change surfaces with highly flexible Random Kitchen Sink (Rahimi and Recht, 2007) features.

5. 5.

   Develop massively scalable additive, non-stationary, non-separable kernels by using the Weyl inequality (Weyl, 1912) and novel Kronecker methods. In addition we integrate our approach into the recent KISS-GP framework (Wilson and Nickisch, 2015). The resulting approach is the first scalable Gaussian process multidimensional changepoint model.

6. 6.

   Describe a novel initialization method for spectral mixture kernels (Wilson and Adams, 2013) by fitting a Gaussian mixture model to the Fourier transform of the data. This method provides good starting values for hyperparameters of expressive stationary kernels, allowing for successful optimization over a multimodal parameter space.

7. 7.

   Demonstrate that the GPCS approach is robust to misspecification, and automatically discourages extraneous model complexity, leading to the discovery of interpretable generative hypotheses for the data.

8. 8.

   Perform counterfactual prediction in complex real world data with posterior mean and covariance estimates for each point in the input domain.

9. 9.

   Use GPCS for discovering and characterizing continuous changes in large observational data. We demonstrate our approach on a recently released public health dataset providing new insight that suggests how the effect of the 1963 measles vaccine may have varied over space and time in the United States. Additionally, we apply the model to requests for lead testing kits in New York City from 2014-2016. The results illustrate distinct spatial patterns in increased concern about lead-tainted water.

The paper is divided into three main units.

Section 2 formally introduces the notion of change surfaces as a multidimensional, expressive generalization of changepoints. We discuss a variant of change surfaces in section 2.1 and detail how to use change surfaces for counterfactual prediction in section 2.2. The discussion of change surfaces in this unit is method-agnostic, and should be relevant to experts from a wide variety of statistical and machine learning disciplines. We emphasize the novel contribution of this framework to the general field of change detection.

Section 3 presents the Gaussian Process Change Surface (GPCS) as a scalable method for change surface modeling. We review Gaussian process basics in section 3.1. We specify the GPCS model in section 3.2. Counterfactual predictions with GPCS are derived in section 3.3. Scalable inference using novel Kronecker methods are presented in section 3.4, and we describe a novel initialization technique for expressive Gaussian process kernels in section 3.5.

Section 4 demonstrates GPCS on out-of-class numerical data and complex spatio-temporal data. We describe our numerical setup in section 4.1 presenting results for posterior prediction, change surface identification, and counterfactual prediction. We present a one-dimensional application of GPCS on coal mining data in section 4.2 including a comparison to state-of-the-art changepoint methods. Moving to spatio-temporal data, we apply GPCS to model requests for lead testing kits in New York City in section 4.3 and discuss the policy relevant conclusions. Additionally, we use GPCS to model measles incidence in the United States in section 4.4 and discuss scientifically relevant insights.

Finally, we conclude with summary remarks in section 5.

In certain applications we are interested in modeling how a change occurs concurrent with a background function which is common to all regimes. For example, consider urban crime. If a police department staged a prolonged intervention in one sector of the city, we expect that some of the crime dynamics in that sector might change. However, seasonal and other weather-related patterns may remain the same throughout the entire city. In this case we want a model to identify and isolate those general background patterns as well as one or more clearly interpretable functions representing regions of change from the background distribution.

We can accommodate such a model as a special case of the generalized change surface from Eq. (2). Each latent function is modeled as $f_{0}(x)+f_{i}(x)$ where $f_{0}(x)$ models “background” dynamics, and $f_{i}(x)$ models each change function. Since changes do not necessarily persist over the entire domain, we fix $f_{r}(x)=0$, and allow $\sum_{i=1}^{r-1}s_{i}(x)\leq 1$. This approach results in the following change surface background model:

Figure 3 presents a two-dimensional representation of the change surface and change surface background models. The data depicted comes from the numerical experiments in section 4.1.

The explicit decomposition into background and change functions is valuable, for instance, if we wish to model counterfactuals: we want to know what the data in a region might look like had there been no change. The decomposition also enables us to interpret the precise effect of each change. Moreover, from a statistical perspective, the decomposition allows us to naturally encode inductive biases into the change surface model, allowing meaningful a priori statistical dependencies between each region. In the particular case of $r=2$, the change surface background model has the form $y(x)=f_{0}(x)+s_{1}(x)f_{1}(x)$, where $f_{1}(x)$ is the only change function modulated by a change surface, $s_{1}(x)\in[0,1]$. This corresponds to observation studies or natural experiments where a single change is observed in the data. We explore this special case further in our discussion of counterfactual prediction, in section 2.2.

Finally, for any change surface or change surface background model, it is critical that the model not overfit the data due to a proliferation of parameters, which could lead to erroneously detected changes even when no dynamic change is present. We discuss one strategy for preventing overfitting through the use of Gaussian processes in section 3.

By simultaneously characterizing the change surface, $s(x)$, and the underlying generative functions, $f(x)$, change surface models allow us to ask questions about how the data would have looked had there been only one latent function. In other words, change surface models allow us to consider counterfactual questions.

For example, in section 4.4 we consider measles disease incidence in the United States in the twentieth century. The measles vaccine was introduced in 1963, radically changing the dynamics of disease incidence. Counterfactual studies such as van Panhuis et al. (2013) attempt to estimate how many cases of measles there would have been in the absence of the vaccine. To be clear, since change surface models do not consider explicit indicators of an intervention, they do not directly estimate the counterfactual with respect to a particular treatment variable such as vaccination. Instead, they identify and characterize changes in the data generating process that may or may not correspond to a known intervention. The change surface counterfactuals estimate the $y$ values for each functional regime in the absence of the change identified by the change surface model. In cases where the discovered change surface does correspond to a known intervention of interest, domain experts may interpret the change surface predictions as a counterfactual “what if” that intervention and any contemporaneous changes in the data generating process (note that we cannot disentangle these causal factors without explicit intervention labels) did not occur.

Counterfactuals are typically studied in econometrics. In observational studies econometricians try to measure the effect of a “treatment” over some domain. Econometric models often measure simple features of the intervention effect, such as the expected value of the treatment over the entire domain, also known as the average treatment effect. A nascent body of work considers machine learning approaches to provide counterfactual prediction in complex data (Athey and Imbens, 2006; Brodersen et al., 2015; Johansson et al., 2016; Hartford et al., 2016), as well as richer measures of the intervention effect (Athey and Imbens, 2006; McFowland et al., 2016). Recent work by Schulam and Saria (2017) uses Gaussian processes for trajectory counterfactual prediction over time. However, these methods generally follow a common framework using the potential outcomes model, which assumes that each observation is observed with a discrete treatment (Rubin, 2005; Holland, 1986). With discrete treatments a unit, $x$, is either intervened upon or not intervened upon — there are no partial interventions. For example, in a medical study a patient may be given a vaccination, or given a placebo. Such discretization is similar to a traditional changepoint model where $s(x)\in\{0,1\}$ can only be in one of two states. Yet discrete states prove challenging in practical applications where units may be partially treated or affected through spillover. For example, there may be herd effects in vaccinations whereby a person’s neighbor being vaccinated reduces the risk of infection to the person. Certain econometric models attempt to account for partial treatment such as treatment eligibility (Abadie et al., 2002), where partial treatments are induced by defining proportions of the population that could potentially be treated. Yet a model that directly enables and estimates continuous levels of treatment may be more natural in such cases.

We provide a brief review of Gaussian processes. More detail can be found in Rasmussen and Williams (2006), Schölkopf and Smola (2002), and MacKay (1998).

Consider data, $(x,y)$, as in section 2, where $x=\{x_{1},\dots,x_{n}\},x_{i}\in\mathbb{R}^{D}$, are inputs or covariates, and $y=\{y_{1},\dots,y_{n}\}$, $y_{i}\in\mathbb{R}$ are outputs or response variables indexed by $x$. We assume that $y$ is generated from $x$ by a latent function with a Gaussian process prior (GP) and Gaussian noise. In particular,

A Gaussian process is a nonparametric prior over functions completely specified by mean and covariance functions. The mean function, $\mu(x)$, is the prior expectation of $f(x)$, while the covariance function, $k(x,x^{\prime})$, is a positive semidefinite kernel that defines the covariance between function values $f(x)$ and $f(x^{\prime})$.

Any finite collection of function values is normally distributed $[f(x_{1})...f(x_{p})]\sim\mathcal{N}(\mu(x),K)$ where $p\times p$ matrix $K_{i,j}=k(x_{i},x_{j})$. Thus we can draw samples from a Gaussian process at a finite set of points by sampling from a multivariate Gaussian distribution. In this paper we generally consider $\mu(x)=0$ and concentrate on the covariance function. The choice of kernel is particularly important in Gaussian process applications since the kernel defines the types of correlations encoded in the Gaussian process. For example, a common kernel choice is a Radial Basis Function (RBF), also known as a Gaussian kernel,

where $s^{2}$ is the signal variance and $V$ is a diagonal matrix of bandwidths. The RBF kernel implies that nearby values are more highly correlated. While this may be true in many applications, it would be inappropriate for data with significant periodicity. In such cases a periodic kernel would be more fitting. We consider more expressive kernel representations in section 3.2.2. This formulation of Gaussian processes naturally accommodates inputs $x$ of arbitrary dimensionality.

Change surface data consists of latent functions $f_{1},\dots,f_{r}$ defining $r$ regimes in the data. The change surface defines the transitions between these functions. We could initially consider an input-dependent mixture model such as in Wilson et al. (2012),

where the weighting functions, $w_{i}(x):\mathbb{R}^{D}\rightarrow\mathbb{R}^{1}$, describe the mixing proportions over the input domain. However, for data with changing regimes we are particularly interested in latent functions that exhibit some amount of mutual exclusivity.

We induce this partial discretization with $\sigma(z):\mathbb{R}^{r}\rightarrow[0,1]^{r}$. These functions have support over the entire real line, but a range in $[0,1]$ and concentrated towards $0$ and $1$. Thus, each $w_{i}(x)$ in Eq. (13) becomes $\sigma_{i}(w(x))$, where $w(x)=[w_{1}(x),..,w_{r}(x)]$. Additionally, we choose $\sigma(z)$ such that it produces a convex combination over the weighting functions, $\sum_{i=1}^{r}\sigma_{i}(w(x))=1$. In this way, each $w_{i}(x)$ defines the strength of latent $f_{i}$ over the domain, while $\sigma(z)$ normalizes these weights to induce weak mutual exclusivity. Thus considering the general model of change surfaces in Eq. (1) we define each warping function as $s_{i}(x)=\sigma_{i}(w(x))$.

A natural choice for flexible change surfaces is to let $\sigma(z)$ be the softmax function. In this way the change surface can approximate a Heaviside step function, corresponding to the sharp transitions of standard changepoints, or more gradual changes. For $r$ latent functions, the resulting warping function is:

The Gaussian process change surface (GPCS) model is thus

where each $f_{i}$ is drawn from a Gaussian process. Importantly, we expect that each Gaussian process, $f_{i}(x)$, will have different hyperparameter values corresponding to different dynamics in the various regimes.

Since a sum of Gaussian processes is a Gaussian process, we can re-write Eq. (15) as $y(x)=f(x)+\epsilon$, where $f(x)$ has a single Gaussian process prior with covariance function,

In this form we can see that $\sigma_{1}(w(x))\dots\sigma_{r}(w(x))$ induce non-stationarity since they are dependent on the input $x$. Thus, even if we use stationary kernels for all $k_{i}$, GPCS observations follow a Gaussian process with a flexible, non-stationary kernel.

We consider counterfactuals when using two latent functions in a GPCS, $f_{1}(x)$ and $f_{2}(x)$. This two-function setup addresses a typical setting for counterfactual prediction when considering two alternatives. The derivations below can be extended to multiple functional regimes. As discussed above, we note that change surface counterfactuals are only valid with respect to the regimes of the data as identified by GPCS. Subsequent analysis and domain expertise are necessary to make any further claims about the relationship between an identified change surface and some latent intervention.

In counterfactual prediction we wish to infer the value of $f_{1}(x)$ and $f_{2}(x)$ in the absence of the other function. Therefore we condition on the observations, $(x,y)$, and GPCS model parameters in order to compute the conditional distribution $p\big{(}[f_{1}(x),f_{2}(x)]|y\big{)}$ from the multivariate Gaussian joint distribution $p\big{(}[f_{1}(x),f_{2}(x)],y\big{)}$. For notational convenience we omit explicit reference to the model parameters in the subsequent derivations but note that all distributions are conditional on these parameters.

To recall, for two latent functions, $f_{1}(x)$ and $f_{2}(x)$, GPCS specifies

where for notational simplicity we let $K_{1}=k_{1}(x,x^{\prime})$, $K_{2}=k_{2}(x,x^{\prime})$, $\sigma_{1}=\sigma_{1}(w(x))$, and $\sigma_{2}=\sigma_{2}(w(x))$.

We consider the most general case when we want to predict counterfactuals for both $f_{1}(x)$ and $f_{2}(x)$ over the domain $X$. No restrictions are placed over $X$. It can include the entire original domain, parts of the original domain, or different inputs entirely. We concatenate $f(X)$ and $g(X)$ together,

Since in section 3.2 we assumed that $f_{1}(x)$ and $f_{2}(x)$ have independent Gaussian process priors, we know that,

Considering the observed data, $y$, we know that $u$ and $y$ are jointly Gaussian,

and using multivariate Gaussian identities, we find that $u$ has the conditional Gaussian distribution

Thus in order to derive counterfactuals for both $f(X)$ and $g(X)$ we only need to compute $\Sigma_{u,y},\Sigma_{y,y}$, and $\Sigma_{u,u}$. Note that with respect to $\Sigma_{u,u}$ we have already derived the covariance structure for $u$ in Equation (29).

Analytic optimization and inference for Gaussian processes requires computation of the log marginal likelihood from Eq. (12). Yet solving linear systems and computing log determinants over $n\times n$ covariance matrices, using standard approaches such as the Cholesky decomposition, requires $O(n^{3})$ computations and $O(n^{2})$ memory, which is impractical for large datasets. Recent advances in scalable Gaussian processes (Wilson, 2014) have reduced this computational burden by exploiting Kronecker structure under two assumptions: (1) the inputs lie on a grid formed by a Cartesian product, $x\in X=X^{(1)}\times...\times X^{(D)}$; and, (2) the kernel is multiplicative across each dimension. Multiplicative kernels are commonly employed in spatio-temporal Gaussian process modeling (Martin, 1990; Majumdar et al., 2005; Flaxman et al., 2015), corresponding to a soft a priori assumption of independence across input dimensions, without ruling out posterior correlations. The popular RBF and ARD kernels, for instance, already have this multiplicative structure. Under these assumptions, the $n\times n$ covariance matrix $K=K_{1}\otimes\dots\otimes K_{D}$, where each $K_{d}$ is $n_{d}\times n_{d}$ such that $\prod_{1}^{D}n_{d}=n$.

Using efficient Kronecker algebra, Saatçi (2011) shows how one can solve linear systems and compute log determinants in $O(Dn^{\frac{D+1}{D}})$ operations using $O(Dn^{\frac{2}{D}})$ memory. Furthermore, Wilson et al. (2014) extends the Kronecker methods for incomplete grids. Yet for additive compositions of kernels, such as those needed for change surface modeling in Eq. (16), the resulting sum of matrix Kronecker products does not decompose as a Kronecker product. Thus, the standard Kronecker approaches for scalable inference and learning are inapplicable. Instead, solving linear systems for the kernel inverse can be efficiently carried out through linear conjugate gradients as in Flaxman et al. (2015) that only rely on matrix vector multiplications, which can be performed efficiently with sums of Kronecker matrices.

However, there is no exact method for efficient computation of the log determinant of the sum of Kronecker products. Instead, Flaxman et al. (2015) upper bound the log determinant using the Fiedler bound (Fiedler, 1971) which says that for $n\times n$ Hermitian matrices $A$ and $B$ with sorted eigenvalues $\alpha_{1},\dots,\alpha_{n}$ and $\beta_{1},\dots,\beta_{n}$ respectively,

While efficient, the Fiedler bound does not generalize to more than two matrices.

Since GPCS uses flexible spectral mixture kernels, as well as RKS features for the change surface, the parameter space is highly multimodal. Therefore, it is essential to initialize the model hyperparameters appropriately. Below we present an approach where we first initialize the $w(x)$ RKS features and then use those values in a novel initialization method for the spectral mixture kernels. Like most GP optimization problems, GPCS hyperparameter optimization is non-convex and there are no provable guarantees that the proposed initialization will result in optimal solutions. However, it is our experience that this initialization procedure works well in practice for the GPCS as well as spectral mixture kernels in general.

To initialize $w(x)$ defined by RKS features we first simplify the change surface model by assuming that each latent function, $f_{1},...,f_{r}$, from Eq. (15) is drawn from a Gaussian process with an RBF kernel. Since RBF kernels have far fewer hyperparameters than spectral mixture kernels, starting with RBF kernels helps our approach find good starting vaglues for $w(x)$. Algorithm 1 provides the procedure for initializing this simplified change surface model. Note that depending on the application domain, a model with latent functions defined by RBF kernels may be sufficient as a terminal model.

In the algorithm, we test multiple possible sets of values for $w(x)$ by drawing the hyperparameters $a$, $\omega$, and $b$ from their respective prior distributions (see section 3.2.1) $m_{1}$ number of times. We set reasonable values for hyperparameters in those prior distributions. Specifically, we let $\Lambda=(\frac{\text{range}(x)}{2})^{2}$, $\sigma_{0}=\mbox{std}(y)$, and $\sigma_{n}=\frac{\mbox{mean}(|y|)}{10}$. These choices are similar to those employed in Lázaro-Gredilla et al. (2010).

For each sampled set of $w(x)$ hyperparameters, we sample $m_{2}$ sets of hyperparameters for the RBF kernels and select the set with the highest marginal likelihood. Then we run an abbreviated optimization procedure over the combined $w(x)$ and RBF hyperparameters and select the joint set that achieves the highest marginal likelihood. Finally, we optimize the resulting hyperparameters until convergence.

In order to initialize the spectral mixture kernels, we use the initialized $w(x)$ from above to define the subset $\{x:\sigma_{i}(w(x))>0.5\}$ where each latent function, $f_{i}$ from Eq. (15), is dominant. We then take a Fourier transform of $y(x)$ over each dimension, $x^{(d)}$, of $\{x:\sigma_{i}(w(x))>0.5\}$ to obtain the empirical spectrum in that dimension. Note that we consider each dimension of $x$ individually since we have a multiplicative Q-component spectral mixture kernel over each dimension (Wilson, 2014). Since spectral mixture kernels model the spectral density with $Q$ Gaussians on $\mathbb{R}^{1}$, we fit a 1-dimensional Gaussian mixture model,