Baseline Drift Estimation for Air Quality Data Using Quantile Trend Filtering

TL;DR

This paper introduces a novel quantile trend filtering method with a parallelizable ADMM algorithm for estimating baseline trends in long time series data, especially from low-cost air quality sensors, improving accuracy and computational efficiency.

Contribution

The study extends quantile trend filtering to multiple quantiles with non-crossing constraints and develops a scalable ADMM algorithm for large datasets.

Findings

Better quantile trend estimates than existing methods

Improved signal classification in air quality data

Enhanced computational efficiency for long time series

Abstract

We address the problem of estimating smoothly varying baseline trends in time series data. This problem arises in a wide range of fields, including chemistry, macroeconomics, and medicine; however, our study is motivated by the analysis of data from low cost air quality sensors. Our methods extend the quantile trend filtering framework to enable the estimation of multiple quantile trends simultaneously while ensuring that the quantiles do not cross. To handle the computational challenge posed by very long time series, we propose a parallelizable alternating direction method of moments (ADMM) algorithm. The ADMM algorthim enables the estimation of trends in a piecewise manner, both reducing the computation time and extending the limits of the method to larger data sizes. We also address smoothing parameter selection and propose a modified criterion based on the extended Bayesian Information Criterion. Through simulation studies and our motivating application to low cost air quality sensor data, we demonstrate that our model provides better quantile trend estimates than existing methods and improves signal classification of low-cost air quality sensor output.

Figures

Equations

$$\mathbf{y}$$

$$\underset{\theta }{\mathrm{minimize}}\frac{1}{2}\Vert \mathbf{y}-\theta {\Vert }_2^2+\lambda \Vert {\mathbf{D}}^{(k+1)}\theta {\Vert }_1,$$

$${\mathbf{D}}^{(1)}$$

$${\mathbf{D}}^{(k+1)}$$

$$\underset{\theta }{\mathrm{minimize}}{\rho }_{\tau }(\mathbf{y}-\theta )+\lambda \Vert {\mathbf{D}}^{(k+1)}\theta {\Vert }_1,$$

$${\rho }_{\tau }(\mathbf{r})$$

$$\sum _{i=1}^n{\rho }_{\tau }(y(t_i)-\theta (t_i))+\lambda \int {\theta }^{\prime \prime }(t{)}^{2}dt,$$

$$\underset{\mathbf{\Theta }\in \mathcal{C}}{\mathrm{minimize}}\sum _{j=1}^J\left[{\rho }_{{\tau }_j}(\mathbf{y}-{\theta }_j)+{\lambda }_j\Vert {\mathbf{D}}^{(k+1)}{\theta }_j{\Vert }_1\right]$$

$$1=l_1<l_2<u_1<l_3<u_2<l_4<u_3<\cdots <u_W=n.$$

$${\mathbf{U}}^{(w)}$$

$$\begin{array}{rl}& \text{minimize}\sum _{w=1}^W\left\{\sum _{j=1}^J\left[{\rho }_{{\tau }_j}({\mathbf{y}}^{(w)}-{\theta }_j^{(w)})+{\lambda }_j\Vert {\mathbf{D}}^{(k+1)}{\theta }_j^{(w)}{\Vert }_1\right]+{\iota }_{\mathcal{C}}({\mathbf{\Theta }}^{(w)})\right\}\\ & \text{subject to}{\mathbf{\Theta }}^{(w)}={\mathbf{U}}^{(w)}\mathbf{\Theta }\text{ for }w=1,\dots ,W.\end{array}$$

$$\begin{array}{rl}& \text{minimize}f(\varphi )+g(\tilde{\varphi })\\ & \text{subject to}\mathbf{A}\varphi +\mathbf{B}\tilde{\varphi }=\mathbf{c}.\end{array}$$

$${\mathcal{L}}_{\gamma }(\varphi ,\tilde{\varphi },\omega )$$

$$\begin{array}{rl}{\varphi }_{m+1}& =\underset{\varphi }{\mathrm{argmin}}{\mathcal{L}}_{\gamma }(\varphi ,{\tilde{\varphi }}_m,{\omega }_m)\\ {\tilde{\varphi }}_{m+1}& =\underset{\tilde{\varphi }}{\mathrm{argmin}}{\mathcal{L}}_{\gamma }({\varphi }_{m+1},\tilde{\varphi },{\omega }_m)\\ {\omega }_{m+1}& ={\omega }_m+\gamma (\mathbf{c}-\mathbf{A}{\varphi }_{m+1}-\mathbf{B}{\tilde{\varphi }}_{m+1}).\end{array}$$

$$\mathcal{L}(\mathbf{\Theta },\{{\mathbf{\Theta }}^{(w)}\},\{{\mathbf{\Omega }}^{(w)}\})$$

$${\mathcal{L}}_w(\mathbf{\Theta },{\mathbf{\Theta }}^{(w)},{\mathbf{\Omega }}^{(w)})$$

$${\mathbf{\Theta }}_{m+1}$$

$${\mathbf{\Theta }}_{m+1}^{(w)}$$

$${\theta }_{ij}$$

$$r_{\text{primal}}^m$$

$$r_{\text{dual}}^m$$

$$r_{\text{primal}}^m$$

$$r_{\text{dual}}^m$$

$${\tilde{\rho }}_{\tau }(\mathbf{r})$$

$$\underset{\mathbf{\Theta }\in \mathcal{C}}{\mathrm{minimize}}\sum _{j=1}^J\left[{\tilde{\rho }}_{{\tau }_j}(\mathbf{y}-{\theta }_j)+{\lambda }_j\Vert {\mathbf{D}}^{(k+1)}{\theta }_j{\Vert }_1\right],$$

$${\breve{\rho }}_{\tau }(\mathbf{r})$$

$$\text{mbox}SIC(p_{\lambda },{\tau }_j)$$

$$\text{mbox}BIC(\nu )=-2\log (L\{\hat{\theta }\})+\nu \log n$$

$$L(\mathbf{y}\mid \theta )$$

$$\text{mbox}BIC(\nu ,{\tau }_j)=\frac{2}{\sigma }{\rho }_{{\tau }_j}(\mathbf{y}-{\hat{\theta }}_j)+\nu \log n$$

$$\text{mbox}{eBIC}_{\gamma }(\nu )$$

Full text

air quality,

non-parametric quantile regression,

trend estimation,

,

and

1 Introduction

In the last decade, low cost and portable air quality sensors have enjoyed dramatically increased usage. These sensors can provide an un-calibrated measure of a variety of pollutants in near real time, but deriving meaningful information from sensor data remains a challenge (Snyder et al., 2013). The “SPod” is a low-cost sensor currently being investigated by researchers at the U.S. Environmental Protection Agency to detect volatile organic compound (VOC) emissions from industrial facilities (Thoma et al., 2016). Due to changes in temperature and relative humidity the output signal exhibits a slowly varying baseline drift on the order of minutes to hours. Figure 1 provides an example of measurements from three SPod sensors co-located at the border of an industrial facility. All of the sensors respond to the pollutant signal, which is illustrated by the three sharp transient spikes at 11:32, 14:10, and 16:03. However, the baseline drift varies from one sensor to another, obscuring the detection of the peaks that alert the intrusion of pollutants. We show later that by estimating the baseline drift in each sensor and removing it from the observed signals, peaks can be reliably detected from concordant residual signals from a collection of SPods using a simple data-driven thresholding strategy. Thus, accurately demixing a noisy observed time series into a slowly varying component and a transient component can lead to greatly improved and simplified downstream analysis.

While this work is motivated by the analysis of data from low cost air quality sensors, the problem of demixing nosy time series into trends and transients is ubiquitous across many fields of study. In a wide range of applications that spans chemistry (Ning, Selesnick and Duval, 2014), macroeconomics (Yamada, 2017), environmental science (Brantley et al., 2014), and medical sciences (Pettersson et al., 2013; Marandi and Sabzpoushan, 2015), scalar functions of time $y(t)$ are observed and assumed to be a superposition of an underlying slowly varying baseline trend $\theta(t)$, other more rapidly varying components $s(t)$, and noise. In practice, $y(t)$ is observed at discrete time points $t_{1},\ldots,t_{n}$, and we model the vector of samples ${y}_{i}=y(t_{i})$ as

[TABLE]

where ${\theta}_{i}=\theta(t_{i})$, ${s}_{i}=s(t_{i})$, and ${\bm{\mathbf{{\varepsilon}}}}\in\mathbb{R}^{n}$ is a vector of uncorrelated noise. For notational simplicity, for the rest of the paper, we assume that the time points take on the values $t_{i}=i$, but it is straightforward to generalize to an arbitrary grid of time points.

In some applications, the slowly varying component ${\bm{\mathbf{{\theta}}}}$ is the signal of interest, and the transient component ${\bm{\mathbf{{s}}}}$ is a vector of nuisance parameters. In our air quality application, the roles of ${\bm{\mathbf{{\theta}}}}$ and ${\bm{\mathbf{{s}}}}$ are reversed; ${\bm{\mathbf{{s}}}}$ represents the signal of interest and ${\bm{\mathbf{{\theta}}}}$ represents a baseline drift that obscures the identification of the important transient events encoded in ${\bm{\mathbf{{s}}}}$.

To tackle demixing problems, we introduce a scalable baseline estimation framework by building on $\ell_{1}$-trend filtering, a relatively new nonparametric estimation framework. Our contributions are three-fold.

• •

Kim et al. (2009) proposed using the check function as a possible extension of $\ell_{1}$-trend filtering but did not investigate it further. Here, we develop the basic $\ell_{1}$-quantile-trend-filtering framework and extend it to model multiple quantiles simultaneously with non-crossing constraints to ensure validity and improve trend estimates.

• •

To reduce computation time and extend the method to long time series, we develop a parallelizable ADMM algorithm. The algorithm proceeds by splitting the time domain into overlapping windows, fitting the model separately for each of the windows and reconciling estimates from the overlapping intervals.

• •

Finally, we propose a modified criterion for performing model selection.

In the rest of the paper, we detail our quantile trend filtering algorithms (Section 2) as well as how to choose the smoothing parameter (Section 3). We demonstrate through simulation studies that our proposed model provides better or comparable estimates of non-parametric quantile trends than existing methods (Section 4). We further show that quantile trend filtering is a more effective method of drift removal for low-cost air quality sensors and results in improved signal classification compared to quantile smoothing splines (Section 5). Finally, we discuss potential extensions of quantile trend filtering (Section 6).

2 Baseline Trend Estimation

2.1 Background

Kim et al. (2009) originally proposed $\ell_{1}$-trend filtering to estimate trends with piecewise polynomial functions, assuming that the observed time series ${\bm{\mathbf{{y}}}}$ consists of a trend ${\bm{\mathbf{{\theta}}}}$ plus uncorrelated noise ${\bm{\mathbf{{\varepsilon}}}}$, namely ${\bm{\mathbf{{y}}}}={\bm{\mathbf{{\theta}}}}+{\bm{\mathbf{{\varepsilon}}}}$. The estimated trend is the solution to the following convex optimization problem

[TABLE]

where $\lambda$ is a nonnegative regularization parameter, and the matrix ${\bm{\mathbf{{D}}}}^{(k+1)}\in\mathbb{R}^{(n-k-1)\times n}$ is the discrete difference operator of order $k+1$. To understand the purpose of penalizing the 1-norm of ${\bm{\mathbf{{D}}}}^{(k+1)}{\bm{\mathbf{{\theta}}}}$ consider the difference operator when $k=0$.

[TABLE]

Thus, $\lVert{\bm{\mathbf{{D}}}}^{(1)}{\bm{\mathbf{{\theta}}}}\rVert_{1}=\sum_{i=1}^{n-1}\lvert\theta_{i}-\theta_{i+1}\rvert$, which is known as the total variation denoising penalty in one dimension in the signal processing literature (Rudin, Osher and Fatemi, 1992) or the fused lasso penalty in the statistics literature (Tibshirani et al., 2005). The penalty term incentivizes solutions which are piecewise constant. For $k\geq 1$, the difference operator ${\bm{\mathbf{{D}}}}^{(k+1)}\in\mathbb{R}^{(n-k-1)\times n}$ is defined recursively as follows

[TABLE]

Penalizing the 1-norm of the vector ${\bm{\mathbf{{D}}}}^{(k+1)}{\bm{\mathbf{{\theta}}}}$ produces estimates of ${\bm{\mathbf{{\theta}}}}$ that are piecewise polynomials of order $k$.

Tibshirani (2014) proved that with a judicious choice of $\lambda$ the trend filtering estimate converges to the true underlying function at the minimax rate for functions whose $k$th derivative is of bounded variation and showed that trend filtering is locally adaptive when the time series consists of only the trend and random noise, which is illustrated in Figure 2a. As noted earlier, in some applications, such as the air quality monitoring problem considered in this paper, the data contain a rapidly varying signal in addition to the slowly varying trend and noise. Figure 2b shows that standard trend filtering is not designed to distinguish between the slowly varying trend and the rapidly-varying signal, as the smooth component estimate ${\bm{\mathbf{{\theta}}}}$ is biased towards the peaks of the transient components.

To account for the presence of transient components in the observed time series ${\bm{\mathbf{{y}}}}$, we propose quantile trend filtering Figure 2b. To estimate the trend in the $\tau$th quantile, we solve the convex optimization problem

[TABLE]

where $\rho_{\tau}({\bm{\mathbf{{r}}}})$ is the check function

[TABLE]

and $\mathbbm{1}(A)$ is 1 if its input $A$ is true and 0 otherwise. Note that we do not explicitly model ${\bm{\mathbf{{s}}}}$. Rather, we focus on estimating ${\bm{\mathbf{{\theta}}}}$. We then estimate ${\bm{\mathbf{{s}}}}+{\bm{\mathbf{{\varepsilon}}}}$ as the difference ${\bm{\mathbf{{y}}}}-{\bm{\mathbf{{\theta}}}}$.

Before elaborating on how we compute our proposed $\ell_{1}$-quantile trend filtering estimator, we discuss similarities and differences between our proposed estimator and existing quantile trend estimators.

2.1.1 Relationship to Prior Work

In this application, as well as those described in Ning, Selesnick and Duval (2014), Marandi and Sabzpoushan (2015), and Pettersson et al. (2013), the goal is to estimate the trend in the baseline not the mean. We can define the trend in the baseline as the trend in a low quantile of the data. A variety of methods for estimating quantile trends have already been proposed. Koenker and Bassett (1978) were the first to propose substituting the sum-of-squares term with the check function (2) to estimate a conditional quantile instead of the conditional mean. Later, Koenker, Ng and Portnoy (1994) proposed quantile trend filtering with $k=1$ producing quantile trends that are piecewise linear, but they did not consider extensions to higher order differences. Rather than using the $\ell_{1}$-norm to penalize the discrete differences, Nychka et al. (1995) used the smoothing spline penalty based on the square of the $\ell_{2}$-norm:

[TABLE]

where $\theta(t)$ is a smooth function of time and $\lambda$ is a tuning parameter that controls the degree of smoothing. Oh, Lee and Nychka (2011) proposed an algorithm for solving the quantile smoothing spline problem by approximating the check function with a differentiable function. Racine and Li (2017) propose a method for estimating quantile trends that does employ the check function. In their method, the response is constrained to follow a location scale model and the conditional quantiles are estimated by combining Gaussian quantile functions with a kernel smoother and solving a local-linear least squares problem.

2.2 Quantile Trend Filtering

We combine the ideas of quantile regression and trend filtering. For a single quantile level $\tau$, the quantile trend filtering problem is given in (1). As with classic quantile regression, the quantile trend filtering problem is a linear program which can be solved by a number of methods. We want to estimate multiple quantiles simultaneously and to ensure that our quantile estimates are valid by enforcing the constraint that if $\tau_{2}>\tau_{1}$ then $Q(\tau_{2})\geq Q(\tau_{1})$ where $Q$ is the quantile function of ${\bm{\mathbf{{y}}}}$. Even if a single quantile is ultimately desired, ensuring non-crossing allows information from nearby quantiles to be used to improve the estimates as we will see in the peak detection experiments in Section 4.2. Given quantiles $\tau_{1}<\tau_{2}<\ldots<\tau_{J}$, the optimization problem becomes

[TABLE]

where ${\bm{\mathbf{{\Theta}}}}=\begin{pmatrix}{\bm{\mathbf{{\theta}}}}_{1}&{\bm{\mathbf{{\theta}}}}_{2}&\cdots&{\bm{\mathbf{{\theta}}}}_{J}\end{pmatrix}\in\mathbb{R}^{n\times J}$ is a matrix whose $j$th column corresponds to the $j$th quantile signal ${\bm{\mathbf{{\theta}}}}_{j}$ and the set \mathcal{C}=\{{\bm{\mathbf{{\Theta}}}}\in\mathbb{R}^{n\times J}:{\theta}_{ij}\leq{\theta}_{ij^{\prime}}\;\text{for j\leq j^{\prime}}\} encodes the non-crossing quantile constraints. The additional non-crossing constraints are linear inequalities involving the parameters, so the non-crossing quantile trends can still be estimated by a number of available linear programming solvers. We allow for the possibility that the degree of smoothness in the trends varies by quantile by allowing the smoothing parameter to vary with quantile as well. In the rest of this paper, we use $k=2$ to produce piecewise quadratic polynomials and report numerical results using the commercial solver Gurobi (Gurobi Optimization, 2018) and its R package implementation. However, we could easily substitute a free solver such as the Rglpk package by Theussl and Hornik (2017).

2.3 ADMM for Big Data

As the size of the data increases, computation time becomes prohibitive. In our application to air quality sensor data, measurements are recorded every second resulting in 86,400 observations per day. This number of observations is already too large to use with currently available R packages used for estimating quantile trends (Douglas Nychka et al., 2017; Koenker, 2018). To our knowledge, no one has addressed the problem of finding smooth quantile trends of series that are too large to be processed simultaneously. We propose a divide-and-conquer approach via an ADMM algorithm for solving large problems in a piecewise fashion.

2.3.1 Formulation

To decrease computation time and extend our method to larger problems, we divide our observed series ${y}_{i}$ with $i=1,\ldots,n$ into $W$ overlapping windows of observations, defining the vector of sequential elements indexed from $l_{w}$ to $u_{w}$ as ${\bm{\mathbf{{y}}}}^{(w)}=\{y_{l_{w}},y_{l_{w}+1},\ldots,y_{u_{w}-1},y_{u_{w}}\}$, with

[TABLE]

We define $n_{w}=u_{w}-l_{w}+1$ so that ${\bm{\mathbf{{y}}}}^{(w)}\in\mathbb{R}^{n_{w}}$. Figure 3 shows an example of 1200 observations being mapped into three equally sized overlapping windows of observations. While the overlapping trend estimates between $l_{2}$ and $u_{2}$ do not vary dramatically, the difference is more pronounced in the trend in the 5th quantile between $l_{3}$ and $u_{2}$. Thus, we need a way of enforcing estimates to be identical in the overlapping regions.

Given quantiles $\tau_{1}<\cdots<\tau_{J}$, we introduce dummy variables ${\bm{\mathbf{{\theta}}}}^{(w)}_{j}\in\mathbb{R}^{n_{w}}$ as the value of the $\tau_{j}$th quantile trend in window $w$. We then “stitch” together the $W$ quantile trend estimates into consensus over the overlapping regions by introducing the constraint ${\theta}^{(w)}_{ij}={\theta}_{i+l_{w}-1,j}$ for $i=1,\ldots,n_{w}$ and for all $j$. Let ${\bm{\mathbf{{\Theta}}}}^{(w)}$ be the matrix whose $j$th column is ${\bm{\mathbf{{\theta}}}}^{(w)}_{j}$. Then we can write these constraints more concisely as ${\bm{\mathbf{{\Theta}}}}^{(w)}={\bm{\mathbf{{U}}}}^{(w)}{\bm{\mathbf{{\Theta}}}}$, where ${\bm{\mathbf{{U}}}}^{(w)}\in\{0,1\}^{n_{w}\times n}$ is a matrix that selects rows of ${\bm{\mathbf{{\Theta}}}}$ corresponding to the $w$th window, namely

[TABLE]

where ${\bm{\mathbf{{e}}}}_{i}\in\mathbb{R}^{n}$ denotes the $i$th standard basis vector. Furthermore, let $\iota_{\mathcal{C}}$ denote the indicator function of the non-crossing quantile constraint, namely $\iota_{\mathcal{C}}({\bm{\mathbf{{\Theta}}}})$ is zero if ${\bm{\mathbf{{\Theta}}}}\in\mathcal{C}$ and infinity otherwise. Our windowed quantile trend optimization problem can then be written as

[TABLE]

The solution to (4) is not identical to the solution to (3) because of double counting of the overlapping sections. The solutions are very close, however, and the differences are essentially immaterial concerning downstream analysis. Figure 4 provides an illustration of the trends estimated using multiple windows compared with the trends estimated using a single window; estimates using multiple and single windows are nearly indistinguishable.

2.3.2 Algorithm

The ADMM algorithm (Gabay and Mercier, 1975; Glowinski and Marroco, 1975) is described in greater detail by Boyd et al. (2011), but we briefly review how it can be used to iteratively solve the following equality constrained optimization problem which is a more general form of (4).

[TABLE]

Recall that finding the minimizer to an equality constrained optimization problem is equivalent to the identifying the saddle point of the Lagrangian function associated with the problem (5). ADMM seeks the saddle point of a related function called the augmented Lagrangian,

[TABLE]

where the dual variable ${\bm{\mathbf{{\omega}}}}$ is a vector of Lagrange multipliers and $\gamma$ is a nonnegative tuning parameter. When $\gamma=0$, the augmented Lagrangian coincides with the ordinary Lagrangian.

ADMM minimizes the augmented Lagrangian one block of variables at a time before updating the dual variable ${\bm{\mathbf{{\omega}}}}$. This yields the following sequence of updates at the $(m+1)$th ADMM iteration

[TABLE]

Returning to our constrained windows problem giving in (4), let ${\bm{\mathbf{{\Omega}}}}^{(w)}$ denote the Lagrange multiplier matrix for the $w$th consensus constraint, namely ${\bm{\mathbf{{\Theta}}}}^{(w)}={\bm{\mathbf{{U}}}}^{(w)}{\bm{\mathbf{{\Theta}}}}$, and let ${\bm{\mathbf{{\omega}}}}^{(w)}_{j}$ denote its $j$th column.

The augmented Lagrangian is given by

[TABLE]

where

[TABLE]

where $\gamma$ is a positive tuning parameter.

The ADMM algorithm alternates between updating the consensus variable ${\bm{\mathbf{{\Theta}}}}$, the window variables $\{{\bm{\mathbf{{\Theta}}}}^{(w)}\}$, and the Lagrange multipliers $\{{\bm{\mathbf{{\Omega}}}}^{(w)}\}$. At the $(m+1)$th iteration, we perform the following sequence of updates

[TABLE]

**Updating ${\bm{\mathbf{{\Theta}}}}$: ** Some algebra shows that, defining $i_{w}=i-l_{w}+1$, updating the consensus variable step is computed as follows.

[TABLE]

The consensus update (7) is rather intuitive. We essentially average the trend estimates in overlapping sections of the windows, subject to some adjustment by the Lagrange multipliers, and leave the trend estimates in non-overlapping sections of the windows untouched. For notational ease, we write the consensus update (7) compactly as ${\bm{\mathbf{{\Theta}}}}=\psi(\{{\bm{\mathbf{{\Theta}}}}^{(w)}\},\{{\bm{\mathbf{{\Omega}}}}^{(w)}\})$.

**Updating $\{{\bm{\mathbf{{\Theta}}}}^{(w)}\}$: ** We then estimate the trend separately in each window, which can be done in parallel, while penalizing the differences in the overlapping pieces of the trends as outlined in Algorithm 1. The use of the Augmented Lagrangian converts the problem of solving a potentially large linear program into a solving a collection of smaller quadratic programs. The gurobi R package (Gurobi Optimization, 2018) can solve quadratic programs in addition to linear programs, but we can also use the free R package quadprog (Weingessel and Turlach, 2013).

Algorithm 1 has the following convergence guarantees.

Proposition 2.1**.**

Let $\{\{{\bm{\mathbf{{\Theta}}}}^{(w)}_{m}\},{\bm{\mathbf{{\Theta}}}}_{m}\}$ denote the $m$th collection of iterates generated by Algorithm 1. Then (i) $\lVert{\bm{\mathbf{{\Theta}}}}^{(w)}_{m}-{\bm{\mathbf{{U}}}}^{(w)}{\bm{\mathbf{{\Theta}}}}_{m}\rVert_{\text{F}}\rightarrow 0$ and (ii) $p_{m}\rightarrow p^{\star}$, where $p^{\star}$ is the optimal objective function value of (4) and $p_{m}$ is the objective function value of (4) evaluated at $\{\{{\bm{\mathbf{{\Theta}}}}^{(w)}_{m}\},{\bm{\mathbf{{\Theta}}}}_{m}\}$.

The proof of Property 2.1 is a straightforward application of the convergence result presented in Section 3.2 of Boyd et al. (2011).

To terminate our algorithm, we use the stopping criteria described by Boyd et al. (2011). The criteria are based on the primal and dual residuals, which represent the residuals for primal and dual feasibility, respectively. The primal residual at the $m$th iteration,

[TABLE]

represents the difference between the trend values in the windows and the consensus trend value. The dual residual at the $m$th iteration,

[TABLE]

represents the change in the consensus variable from one iterate to the next. The algorithm is stopped when

[TABLE]

The quantile trend filtering problem for a single window is a linear program with $N\times J$ parameters (number of observations by number of quantiles), which can be solved in computational time proportional to $(NJ)^{3}$. Consequently, solving a large problem using Algorithm 1 should require less computational time than solving (3), even if the sub-problems are solved sequentially. We demonstrate the advantages of Algorithm 1 through timing experiments (Figure 5). For each data size, 25 datasets were simulated using the peaks simulation design described below. Trends for the fifth, tenth, and fifteenth quantiles were fit simultaneously using $\lambda_{j}=n/5$ for all $j$. We use from one to four windows for each data size with an overlap of 500. Algorithm 1 was the stopped when (8) was satisfied, defining $\epsilon_{abs}=0.01$ and $\epsilon_{rel}=0.001$. Figure 5 shows that using 4 windows instead of one on data sizes of 55,000 provides a factor of 3 decrease in computation time. The timing experiments were conducted on an Intel Xeon based Linux cluster using two processor cores.

3 Model Selection

An important practical issue in baseline estimation is the choice of the regularization parameter $\lambda$, which controls the degree of smoothness in ${\bm{\mathbf{{\theta}}}}$. In this section, we introduce four methods for choosing $\lambda$. The first is a validation based approach; the latter three are based on information criteria. Each of the criteria we compare is calculated for a single quantile ($\tau_{j}$). Rather than combine results over quantiles, we allow the value of $\lambda$ to vary by quantile resulting in ${\bm{\mathbf{{\lambda}}}}=\{\lambda_{1},...,\lambda_{J}\}$. To choose the best value for each $\lambda_{j}$, we first estimate all of the quantile trends using $\lambda_{j}=\lambda$ for all $j$ over a grid of values for $\lambda$. We then determine the $\lambda_{j}$ that maximizes the criteria chosen evaluated using $\tau_{j}$. Finally, we re-estimate the non-crossing trends with the optimal values for $\lambda_{j}$. A more thorough approach would involve fitting the model on a J dimensional grid of values for ${\bm{\mathbf{{\lambda}}}}$ but this is computationally infeasible.

3.1 Validation

Our method can easily handle missing data by defining the check loss function to output 0 for missing values. Specifically, we use the following modified function $\tilde{\rho}_{\tau}$ in place of the $\rho_{\tau}$ function given in (2)

[TABLE]

where $\mathcal{V}$ is a held-out validation subset of $\{1,\ldots,n\}$ and solve the problem

[TABLE]

which can be solved via Algorithm 1 with trivial modification to the quadratic program sub-problems. For each quantile level $j$, we select the $\lambda_{j}$ that minimizes the hold-out prediction error quantified by $\breve{\rho}_{\tau_{j}}({\bm{\mathbf{{y}}}}-{\bm{\hat{\mathbf{{\theta}}}}}_{j}(\lambda_{j}))$ where ${\bm{\hat{\mathbf{{\theta}}}}}_{j}(\lambda_{j})$ is the solution to (10) and

[TABLE]

3.2 Information Criteria

Koenker, Ng and Portnoy (1994) addressed the choice of regularization parameter by proposing the Schwarz criterion for the selection of $\lambda$

[TABLE]

where $p_{\lambda}=\sum_{i}\mathbbm{1}({y}_{i}=\widehat{\theta}_{i})$ is the number of non-interpolated points, which can be thought of as active knots. Equivalently, $p_{\lambda}$ can be substituted with the number of non-zero components in ${\bm{\mathbf{{D}}}}^{(k+1)}{\bm{\hat{\mathbf{{\theta}}}}}_{j}$ which we denote $\nu$ and have found to be more numerically stable. The SIC is based on the traditional Bayesian Information Criterion (BIC) which is given by

[TABLE]

where $L$ is the likelihood function. If we take the approach used in Bayesian quantile regression (Yu and Moyeed, 2001), and view minimizing the check function as maximizing the asymmetric Laplace likelihood,

[TABLE]

we can compute the BIC as

[TABLE]

where ${\bm{\hat{\mathbf{{\theta}}}}}$ is the estimated trend, and $\nu$ is the number of non-zero elements of ${\bm{\mathbf{{D}}}}^{(k+1)}{\bm{\hat{\mathbf{{\theta}}}}}$. We can choose any $\sigma>0$ and have found empirically that $\sigma=\frac{1-|1-2\tau|}{2}$ produces stable estimates.

Chen and Chen (2008) proposed the extended Bayesian Information Criteria (eBIC), specifically designed for large parameter spaces.

[TABLE]

where $P$ is the total number of possible parameters and $\nu$ is the number of non-zero parameters included in given mod. We used this criteria with $\gamma=1$, and $P=n-k-1$. In the simulation study, we compare the performance of the SIC, scaled eBIC (with $\sigma$ defined above), and validation methods.

4 Simulation Studies

We conduct two simulation studies to compare the performance of our quantile trend filtering method and regularization parameter selection criteria with previously published methods. The first study compares the method’s ability to estimate quantiles when the observed series consists of a smooth trend plus independent error, but does not contain transient components. The second study is based on our application and compares the method’s ability to estimate baseline trends and enable peak detection when the time series contains a non-negative, transient signal in addition to the trend and random component.

We compare three criteria for choosing the smoothing parameter for quantile trend filtering: ${\bm{\mathbf{{\lambda}}}}$ chosen using SIC (3.2) (detrendr_SIC); ${\bm{\mathbf{{\lambda}}}}$ chosen using the validation method with the validation set consisting of every 5th observation (detrendr_valid); and ${\bm{\mathbf{{\lambda}}}}$ chosen using the proposed eBIC criterion (3.2) (detrendr_eBIC). For the second study we also examine the effect of the non-crossing quantile constraint by estimating the quantile trends separately and choosing ${\bm{\mathbf{{\lambda}}}}$ using eBIC (detrendr_Xing). We do not include detrendr_Xing in the first study because the difference in quantiles is large enough that we would not expect the non-crossing constraint to make a difference.

We also compare the performance of our quantile trend filtering method with three previously published methods, none of which guarantee non-crossing quantiles:

• •

npqw: The local linear quantile method (quantile-ll) described in Racine and Li (2017). Code was obtained from the author.

• •

qsreg: Quantile smoothing splines described in Oh et al. (2004) and Nychka et al. (1995) and implemented in the fields R package (Douglas Nychka et al., 2017). The regularization parameter was chosen using generalized cross-validation.

• •

rqss: Quantile trend filtering with $k=1$ available in the quantreg package and described in Koenker, Ng and Portnoy (1994). The regularization parameter is chosen using a grid search and minimizing the SIC (3.2) as described in Koenker, Ng and Portnoy (1994).

4.1 Estimating Quantiles

To compare performance in estimating quantile trends in the absence of a signal component, three simulation designs from Racine and Li (2017) were considered. For all designs $t=1,\ldots,n$, $x(t)=t/n$, and the response ${\bm{\mathbf{{y}}}}$ was generated as

[TABLE]

The errors were simulated as independent draws from the following distributions:

• •

Gaussian: $\epsilon(x(t))\sim N\left(0,\left(\frac{1+x(t)^{2}}{4}\right)^{2}\right)$

• •

Beta: $\epsilon(x(t))\sim Beta(1,11-10x(t))$

• •

Mixed normal: $\epsilon(x(t))$ is simulated from a mixture of $N(-1,1)$ and $N(1,1)$ with mixing probability $x(t)$.

The true quantile trends and an example simulated data set is show in Figure 6. One hundred datasets were generated of sizes 300, 500 and 1000.

Quantile trends were estimated for $\tau=\{0.05,0.25,0.5,0.75,0.95\}$ and the root mean squared error was calculated as $\mbox{RMSE}(\tau_{j})=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\hat{\theta}_{ij}-\theta_{ij})^{2}}$, where $\theta_{ij}$ is the true value of the $\tau_{j}$th quantile of $y$ at $t=i$. Figure 7 shows the mean RMSE plus or minus twice the standard error for each method, quantile level, and sample size. In all three designs the proposed detrend methods are either better than or comparable to existing methods. Overall the detrend_eBIC performs best. In the mixed normal design, specifically, our methods have lower RMSEs for the 5th and 95th quantiles. The npqw method performs particularly poorly on the mixed normal design due to the violation of the assumption that the data come from a scale-location model.

4.2 Peak Detection

The second simulation design is closely motivated by our air quality analysis problem. We assume that the measured data can be represented by

[TABLE]

where $t_{i}=i$ for $i=1,...,n$, $\theta(t)$ is the drift component that varies smoothly over time, $s(t)$ is the true signal at time $t$, and ${\varepsilon}_{i}$ are i.i.d. errors distributed as $N(0,0.25^{2})$. We generate $\theta(t)$ using cubic natural spline basis functions with degrees of freedom sampled from a Poisson distribution with mean parameter equal to $n/100$, and coefficients drawn from an exponential distribution with rate 1. The true signal function $s(t)$ is assumed to be zero with peaks generated using the Gaussian density function. The number of peaks is sampled from a binomial distribution with size equal to $n$ and probability equal to $0.005$ with location parameters uniformly distributed between $1$ and $n-1$ and bandwidths uniformly distributed between $2$ and $12$. The simulated peaks were multiplied by a factor that was randomly drawn from a normal distribution with mean 20 and standard deviation of 4. An example dataset with 4 signal peaks is shown in Figure 9. One hundred datasets were generated for each $n\in\{500,1000,2000,4000\}$.

We compare the ability of the methods to estimate the true quantiles of $y(t_{i})-s(t_{i})$ for $\tau\in\{0.01,0.05,0.1\}$ and calculate the RMSE (Figure 8). In this simulation study, our proposed method detrend_eBIC method substantially outperforms the others. The detrend_Xing method, which is the detrend_eBIC method fit without the non-crossing constraints, performs similarly for larger quantiles and larger datasets. However, detrend_Xing produces significantly worse estimates for more extreme quantiles ($\tau=0.01$) and smaller data sets ($\tau=0.05$ and $n=500$). These results indicate that even when a single quantile is of interest, simultaneously fitting nearby quantiles and utilizing the non-crossing constraint can improve estimates when data is sparse by using information from nearby quantiles. The qsreg method is comparable to the detrend_eBIC method on the larger datasets, but its performance deteriorates as the data size shrinks. The npqw and rqss methods both perform poorly on this design.

While minimizing RMSE is desirable in general, in our application, the primary metric of success is accurately classifying observations ${y}_{i}$ into signal present or absent. To evaluate the accuracy of our method compared to other methods we define true signal as any time point when the simulated peak value is greater than 0.5. We compare three different quantiles for the baseline estimation and four different thresholds for classifying the signal after subtracting the estimated baseline from the observations. Figure 9 illustrates the observations classified as signal after subtracting the baseline trend compared to the “true signal.”

To compare the resulting signal classifications, we calculate the class averaged accuracy (CAA), which is defined as

[TABLE]

where $\delta_{i}\in\{0,1\}$ is the vector of true signal classifications and $\hat{\delta}_{i}\in\{0,1\}$ is the vector of estimated signal classifications, namely $\hat{\delta}_{i}=\mathbbm{1}({y}_{i}-\hat{\theta}_{i}>0.5)$. We use this metric because our classes tend to be very imbalanced with many more zeros than ones. The CAA metric should give a score close to 0.5 both for random guessing and also for trivial classifiers such as $\hat{\delta}_{i}=0$ for all $i$.

Our detrend_eBIC method results in the largest CAA values (Figure 10) in addition to the smallest RMSE values (Figure 8). While qsreg was competitive with our method in some cases, in the majority of cases the largest CAA values for each threshold were produced using the detrend_eBIC method with the 1st or 5th quantiles.

5 Analysis of Air Quality Data

The low-cost “SPod” air quality sensors output a time series that includes a slowly varying baseline, the sensor response to pollutants, and high frequency random noise. These sensors record measurements every second and are used to monitor pollutant concentrations at the perimeter of industrial facilities. Time points with high concentrations are identified and compared with concurrent wind direction and speed. Ideally, three co-located and time aligned sensors (as shown in Figure 1) responding to a pollutant plume would result in the same signal classification after baseline trend removal and proper threshold choice. We first illustrate the difference between our detrend_eBIC method, hereafter referred to as detrendr, and qsreg using data from 13:10 to 15:10 from Figure 1 (Section 5.1). We then compare the methods on the complete day shown in Figure 1, estimating trends by applying the qsreg method to 2 hour non-overlapping windows of the data and Algorithm 1 to the entire day. We focus on this day because data from three sensors was available. Finally, we examine an entire week of measurements from two co-located sensors (Section 5.2).

5.1 Short series of air quality measurements

We compare our detrendr method with the qsreg method on a two-hour subset of one-second SPod data (n=7200) both to facilitate visualization and because the qsreg method cannot handle all 24 hours simultaneously. We estimate the baseline trend using $\tau=\{0.01,0.05,0.1\}$ and compare three thresholds for classifying the signal. The thresholds are calculated as the 90th, 95th, and 99th quantiles of the de-trended series for each SPod. If there is signal present in the dataset, values above these thresholds should occur simultaneously on all three SPods. We do not use class-averaged accuracy to compare the signal classifications because we do not have a reference value to define as the “true” signal. Instead, we compute the variation of information (VI) which compares the similarity between two classifications. Given the signal classifications for SPods a and b, $\delta^{(a)}_{i}\in\{0,1\}$ and $\delta^{(b)}_{i}\in\{0,1\}$, for $i\in\{1,...,n\}$ the VI is defined as:

[TABLE]

where $(j,k)\in\{(0,0),(0,1),(1,0),(1,1)\}$. The VI is a distance metric for measuring similarity of classifications and will be 0 if the classifications are identical and increase as the classifications become more different.

Figure 11 shows the estimated 5th quantile trends from each method for each SPod. The detrendr method results in a smoother baseline estimate while the qsreg method absorbs more of the peaks obscuring some of the signal. Figure 12 shows the series after subtracting the detrendr estimate of the 5th quantile and classifying the signal using the 95th quantile of the detrended data. The 90th and 99th thresholds are also shown for comparison in blue and orange, respectively. The largest peaks at 14:10 are easily identified as signal, but good baseline estimates also enable proper classification of the smaller peaks like the one at 15:12. The under-smoothing of the qsreg method results in less similar signal classifications and higher VI values for the 90th and 95th quantile thresholds (Figure 13). However, when the 99th threshold is used only the highest observations are classified as signal and the baseline estimation method isn’t as important (Figure 13).

5.2 Long series of air quality measurements

Algorithm 1 was used to remove the baseline drift from the full day of data (Figure 1) consisting of 86,400 observations per SPod and compared to the series detrended using the qsreg trends estimated using non-overlapping 2 hour windows. As in the shorter illustration, the detrendr method results in generally lower VI scores than the qsreg method (Figure 14). The detrendr method also results in better correlation in the detrended series as is illustrated in (Figure 15). The Spearman correlation coefficients for SPods a and b, SPods a and c, and SPods b and c after removing the 5th quantile trend using detrendr were 0.37, 0.75, and 0.43, compared with 0.07, 0.24, and 0.16 using qsreg. The noise variance was higher for SPod b than for SPods a and c resulting in lower correlation and higher VI values for the ab and bc metrics compared with ac.

Finally, we estimated the quantile trends for 7 days of measurements of two co-located SPods. Figure 16 demonstrates the improvement in classification similarity when using detrendr. Each point represents a day of measurements and all points that fall below the dashed line have more similar classifications using detrendr compared to qsreg. The improvement of detrendr over qsreg is more severe at lower thresholds. This indicates that detrendr gives greater agreement on signal classification when the methods are tuned to deliver positive classifications more frequently.

6 Conclusion and Discussion

We have expanded the quantile trend filtering method by implementing a non-crossing constraint, a new algorithm for processing large series, and proposing a modified criteria for smoothing parameter selection. Furthermore, we have demonstrated the utility of quantile trend filtering in both simulations and applied settings. Our ADMM algorithm for large series both reduces the computing time and allows trends to be estimated on series that cannot be estimated simultaneously while our scaled extended BIC criterion was shown to provide better estimates of quantile trends in series with and without a signal component.

In our application to low cost air quality sensor data, we have shown that the baseline drift in low cost air quality sensors can be removed through estimating quantile trends, but the data size was too large for existing methods to be computationally feasible. While qsreg cannot feasibly handle more than a few hour windows of data, our new methods were able to process 24 hours simultaneously and deliver signal classifications that were more consistent between the two sensors for a week of data (168 hours).

In the future, quantile trend filtering could be extended to observations measured at non-uniform spacing by incorporating the distance in covariate spacing into the differencing matrix. It could also be extended to estimate smooth spatial trends by a similar adjustment to the differencing matrix based on spatial distances between observations.

7 Acknowledgments

The authors greatly appreciate Eben Thoma at the US Environmental Protection Agency for providing the SPod datasets. The authors are also thankful for the Oak Ridge Institute of Science Fellowship that partially supported this work. Plots were created using ggplot2 (Wickham, 2016).

SUPPLEMENTARY MATERIAL

R-package for detrend routine:

R-package detrendr containing code to perform the methods described in the article. (GNU zipped tar file also available at https://github.com/halleybrantley/detrendr)