LowCon: A design-based subsampling approach in a misspecified linear model

Suppose the underlying true model has the form

where $y_{i}$’s are the responses, $\bm{x}_{i}$’s are the predictors, $\bm{\beta}_{0}\in\mathbb{R}^{p}$ ($p\ll n$) is the vector of unknown coefficients, the random errors $\{u_{i}\}_{i=1}^{n}$ are independently distributed, and $u_{i}$ follows a non-centered normal distribution $N(h(\bm{x}_{i}),\sigma^{2})$, $i=1,\ldots,n$. Let $\mathcal{X}$ be the design space. In this paper, we assume that the unknown multivariate function $h$ satisfies

where $\alpha$ is a finite positive constant. When $\bm{x}_{i}=(x_{i1},\ldots,x_{ip})^{\intercal}$ has finite values, some examples of $h$ include $h(\bm{x}_{i})=\sin(x_{i1})$ and $h(\bm{x}_{i})=x_{i1}x_{i2}$. Let $\bm{y}=(y_{1},\ldots,y_{n})^{\intercal}$ be the response vector, $\mathbf{X}=(\bm{x}_{1},\ldots,\bm{x}_{n})^{\intercal}$ be the predictor matrix, and $\bm{h}_{X}=(h(\bm{x}_{1}),\ldots,h(\bm{x}_{n}))^{\intercal}$ be the misspecification term. For model identifiability, we assume the matrix $[\mathbf{X};\bm{h}_{X}]$ has a full column rank. Under this assumption, we exclude the case that $h(\bm{x})$ is a linear function of $\bm{x}$, i.e., $h(\bm{x}_{i})$ cannot be a linear combination of $\bm{x}_{i1},\ldots,\bm{x}_{ip}$.

We consider the scenario that practitioners have no prior information on the true model (1) and postulate a classical linear model,

where the random errors $\{\epsilon_{i}\}_{i=1}^{n}$ are i.i.d. and follow a normal distribution with mean zero and constant variance $\sigma^{2}$, i.e., $N(0,\sigma^{2})$. Model (3) is thus a misspecified linear model of the true model (1). Fitting model (3) without taking into account the model misspecification may result in the degenerated performance of the coefficient estimation and model prediction. For example, the full-sample ordinary least squares (OLS) estimate, known as the best linear unbiased estimate, is a biased estimate of the true coefficient when the model is misspecified [5]. More discussion on misspecified linear models can be found in kiefer1975optimal [20] and sacks1978linear [34].

In our measurement-constrained setting, practitioners are given the full sample of predictors $\{\bm{x}_{i}\}_{i=1}^{n}$. The responses $\{y_{i}\}_{i=1}^{n}$ in model (1), however, are hidden unless explicitly requested. Practitioners are then allowed to reveal a subset of $\{y_{i}\}_{i=1}^{n}$, denoted by $\bm{y}^{*}=(y^{*}_{1},\ldots,y^{*}_{r})^{\intercal}$, where $p<r\ll n$. The goal is to estimate the true coefficient $\bm{\beta}_{0}$ using $(\bm{x}_{i}^{*},y_{i}^{*})$, where $i=1,\ldots,r$, and $\bm{x}_{i}^{*}$ is the corresponding predictor for $y_{i}^{*}$. A natural estimator for the coefficient $\bm{\beta}_{0}$ is the subsample least squares estimator [50],

where $\mathbf{X}^{*}=(\bm{x}_{1}^{*},\ldots,\bm{x}_{r}^{*})^{\intercal}$. We derive the mean squared error (MSE) and the worst-case MSE of this estimator, in the next subsection.

Let $\mathbf{Q}=(\mathbf{X}^{*\intercal}\mathbf{X}^{*})^{-1}\mathbf{X}^{*\intercal}$ and $\bm{h}=(h(\bm{x}^{*}_{1}),\ldots,h(\bm{x}^{*}_{r}))^{\intercal}\in\mathbb{R}^{r}$. The MSE of the estimator $\widetilde{\bm{\beta}}_{\mathbf{X}^{*}}$ (conditional on $\mathbf{X}$) thus can be decomposed as

where the bias term $\bm{h}^{\intercal}\mathbf{Q}^{\intercal}\mathbf{Q}\bm{h}$ is associated with the model misspecification. Note that when the bias term vanishes, $\bm{h}_{X}=\mathbf{0}$, i.e., when the model is correctly specified, minimizing MSE is equivalent to minimizing the variance term $\sigma^{2}\text{tr}[(\mathbf{X}^{*\intercal}\mathbf{X}^{*})^{-1}]$. Further discussion following this line of thinking can be found in wang2017computationally [50] and wang2018information [48], in which the authors focused on selecting the subsample that minimizes the variance term. In our setting, where the model is misspecified, however, minimizing the variance term does not necessarily lead to a small MSE.

Recall that our goal is to select a subsample such that the corresponding SLS estimator is robust to various model misspecification. Since the misspecification term $\bm{h}_{X}$ is unknown to practitioners, a natural and intuitive approach is to find the “minimax” subsample that minimizes the so-called “worst-case” MSE, i.e., the maximum value of $\text{MSE}(\widetilde{\bm{\beta}}_{\mathbf{X}^{*}})$ with respect to all the possible choices of the misspecification term $\bm{h}_{X}$. The following lemma gives an explicit form of the worst-case MSE; the proof can be found in the Supplementary Material.

Two conclusions can be made from Lemma 2.1. First, the worst-case MSE of an SLS estimator can be inflated to arbitrarily large values by a very small value of $\lambda_{min}(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$. It is thus very easy to construct a “worst-case” sample and a misspecification term for which an SLS estimator will have unacceptable performance. Second, $\widetilde{\bm{\beta}}_{\mathbf{X}^{*}}$ is the most robust SLS estimator if the selected subsample minimizes the worst-case MSE. Such a subsample, however, is impossible to obtain in real practice, since both values of $\sigma^{2}$ and $\alpha^{2}$ are unknown to practitioners.

In this paper, we are more interested in the setting where the misspecified term $h(\bm{x})$ is large enough. In particular, the value of $\alpha^{2}$ is large enough such that, on the right-hand side of the Inequality (5), the second term dominates the first term. Under this setting, the desired subsample $\mathbf{X}^{*}$ should yield a relatively small value of $\text{tr}(\mathbf{X}^{*\intercal}\mathbf{X}^{*})/\lambda_{min}(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$. Notice that

where the equality holds when the condition number of the subsample information matrix, i.e., $\kappa(\mathbf{X}^{*\intercal}\mathbf{X}^{*})\stackrel{{\scriptstyle\text{def}}}{{=}}\lambda_{max}(\mathbf{X}^{*\intercal}\mathbf{X}^{*})/\lambda_{min}(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$, takes the minimum value $1$. Inequality (6) thus suggests the desired subsample $\mathbf{X}^{*}$ is the one with a relatively small value of $\kappa(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$.

We now give another intuition about how $\kappa(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$ is related to the robustness of the SLS estimator. casella1985condition [6] showed that

where $\delta\widehat{\bm{\beta}}_{ols}$ and $\delta{\mathbf{X}^{\intercal}\bm{y}}$ are perturbations of $\widehat{\bm{\beta}}_{ols}$ and ${\mathbf{X}^{\intercal}\bm{y}}$ respectively. Analogously, one can also show that

Inequality (7) thus suggests that a smaller value of $\kappa(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$ associates with a more robust estimator $\widetilde{\bm{\beta}}_{\mathbf{X}^{*}}$.

It is worth noting that, if the subsample matrix $\mathbf{X}^{*}$ minimizes the worst-case MSE, it does not necessarily minimize $\kappa(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$ simultaneously since both the value of $\sigma^{2}$ and $\alpha^{2}$ are not available in practice. A robust subsample $\mathbf{X}^{*}$ should at least yield a relatively small value of $\kappa(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$ and balance the trade-off between the bias and the variance in the Equation (4). Following this line of thinking, we propose a novel subsampling algorithm, the details of which are presented in the next section.

Taking a subsample with some specific characteristics has many similarities to the design of experiments, which aims to place design points in a continuous design space, so that resulting design points have certain properties [51]. The theory and methods in the design of experiments are potentially useful for solving subsampling problems. The fundamental difference between the design of experiments and the subsampling is that, in subsampling, the selected points cannot be freely designed in a continuous space as the design of experiments but must be taken from the given finite sample $\{\bm{x}_{i}\}_{i=1}^{n}$. To borrow the strength of the design of experiments, we focus on space-filling designs, which aims to place the design points that cover a continuous design space as uniformly as possible [11, 21, 18, 29, 47]. In other words, for any point in the experimental region, space-filling designs have a design point close to it. We thus propose to round the design point to its nearest neighbor in the sample. Details are provided in Section 3.2.

We now introduce a specific space-filling design that is of interest, the Latin hypercube design (LHD) [37, 26, 46].

Intuitively, if one divides the design space $[-1,1]^{p}$ into $r$ equally-sized slices in the $j$th ($j=1,\ldots,p$) dimension, a Latin hypercube design ensures that there is exactly one design point in each slice. The left panel of Figure 2 shows an example of a set of Latin hypercube design points (black dots). Although uniformly distributed on the marginal, the Latin hypercube design points do not necessarily spread out in the whole design space. That is to say, a set of LHD points may not be “space-filling” enough. To improve the “space-filling” property of LHD, various methods have been developed [40, 31, 12, 19]. Of particular interest in this paper is the orthogonal Latin hypercube design (OLHD), which achieves the goal by reducing the pairwise correlations of LHD [53]; see the right panel of Figure 2 for an example.

Consider the information matrix $\mathbf{L}^{\intercal}\mathbf{L}$, where $\mathbf{L}$ is an OLHD matrix. Intuitively, the matrix $\mathbf{L}^{\intercal}\mathbf{L}$ has a relatively small condition number, since all of the diagonal elements of $\mathbf{L}^{\intercal}\mathbf{L}$ are the same and all of the off-diagonal elements of $\mathbf{L}^{\intercal}\mathbf{L}$ have relatively small absolute value. Although there is a lack of theoretical guarantee, empirically, it is known that $\kappa(\mathbf{L}^{\intercal}\mathbf{L})$ is in general no greater than 1.13 [7]. Such a fact motivates us to select the subsample that approximates a set of orthogonal Latin hypercube design points.

Without loss of generality, we assume the data points $\{\bm{x}_{i}\}_{i=1}^{n}$ are first scaled to $[-1,1]^{p}$. The proposed algorithm works as follows. We first generate a set of orthogonal Latin hypercube design points from a design space $\mathcal{X}\subseteq[-1,1]^{p}$. We then search and select the nearest neighbor from the sample for every design point.

The key to success is that the selected subsample can well-represent the set of design points, i.e., each selected subsample point is close-enough to its nearest design point, respectively. We provide more discussion in Section 3.3 about when such a requirement is met in practice. Empirically, we find $[-1,1]^{p}$ may not be a good choice for the design space $\mathcal{X}$. This is because, in such a scenario, the design points, which are close to the boundary of $[-1,1]^{p}$, may be too far away from their nearest neighbors, especially when the population density function has a heavy tail. As a result, a design space that is slightly smaller than $[-1,1]^{p}$ would be a safer choice. We opt to set the design space as $\mathcal{X}_{\theta}=[\theta_{j1},\theta_{j2}]^{p}$, where $\theta_{j1}$ and $\theta_{j2}$ are the $\theta$-percentile and $(100-\theta)$-percentile of the $j$th column of the scaled data points, respectively. The algorithm is summarized below.

Figure 3 illustrates LowCon algorithm. The synthetic data points in the left panel were generated from a bivariate normal distribution and are scaled to $[-1,1]^{2}$. A set of orthogonal Latin hypercube design points are then generated, labeled as black triangles in the middle panel. For each design point, the nearest data point is selected, marked as black dots in the right panel. The selected points can well-approximate the design points.

Note that the set of design points generated by the orthogonal Latin hypercube design technique is not unique; different sets of design points may result in different subsamples. Algorithm 1 thus is a random subsampling method instead of a deterministic subsampling method. In practice, the set of design points $\{\bm{l}_{i}\}_{i=1}^{r}$ in Algorithm 1 can be randomly generated.

We now present the theoretical property of the proposed subsample least squares estimator, obtained by the LowCon algorithm. Some notations are needed before we show our main theorem. Recall that $\mathbf{L}$ represents an orthogonal Latin hypercube design matrix. Let $\mathbf{X}^{*}_{L}$ be the subsample matrix obtained by the proposed algorithm. One thus can decompose $\mathbf{X}^{*}_{L}$ into a sum of the design matrix $\mathbf{L}$ and a matrix $\mathbf{D}=(\bm{d}_{1},\cdots,\bm{d}_{r})^{\intercal}$, i.e., $\mathbf{X}^{*}_{L}=\mathbf{L}+\mathbf{D}$.

Following the notations in Algorithm 1, one can write $\mathbf{L}=(\bm{l}_{1},\ldots,\bm{l}_{r})^{\intercal}$ and $\mathbf{X}^{*}_{L}=(\bm{l}^{*}_{1},\ldots,\bm{l}^{*}_{r})^{\intercal}$, where $\bm{l}_{i}$ and $\bm{l}^{*}_{i}$ represent the $i$th design point and its corresponding nearest neighbor from the sample, respectively. One thus has $\bm{d}_{i}=\bm{l}^{*}_{i}-\bm{l}_{i}$, for $i=1,\ldots,r$. Intuitively, $\mathbf{D}$ is a random perturbation matrix, and the selected data points can well-approximate the design points if $\mathbf{D}$ is “negligible”. In such a case, $MSE(\widetilde{\bm{\beta}}_{\mathbf{X}^{*}_{L}})$, which is a function of $\mathbf{X}^{*}_{L}$, can be expanded around $MSE(\widetilde{\bm{\beta}}_{\mathbf{L}})$ through Taylor expansion. From this, we can establish our main theorem below; the proof is relegated to the appendix.

When the Taylor expansion in Theorem 3.1 is valid, three significant conclusions can be made. First, the theorem indicates that the MSE of the proposed estimator is finite. Specifically, following the Definition 3.1, we have

Moreover, the value of $\kappa(\mathbf{L}^{\intercal}\mathbf{L})$ is in general no greater than 1.13, as discussed in Section 3.1. Combining these two facts yields an informal but finite upper bound for $MSE(\widetilde{\bm{\beta}}_{\mathbf{X}^{*}_{L}})$, i.e.,

Recall that Lemma 2.1 shows that the worst-case MSE of an SLS estimator can be inflated to an arbitrarily large value by a very small value of $\lambda_{min}(\mathbf{X}^{*\intercal}\mathbf{X}^{*})$. The fact that the proposed estimator has a finite MSE thus indicates the proposed estimator is robust.

Second, the upper bound of the squared bias of the proposed estimator, which equals $\alpha^{2}p\kappa(\mathbf{L}^{\intercal}\mathbf{L})$, is very close to the minimum value of the worst-case squared bias. This is because the worst-case squared bias has the minimum value of $\alpha^{2}p$, and the value of $\kappa(\mathbf{L}^{\intercal}\mathbf{L})$ is close to 1. Consider the common situation when the value of $\alpha^{2}$ is large enough such that, in Inequality (5), the bias term dominates the variance term. Under such a situation, the second conclusion thus indicates the proposed estimator is very close to the “most robust” estimator, which minimizes the worst-case squared bias.

Third, the proposed estimator has a finite variance. Recall that in Algorithm 1, sometimes we may choose a design space $\mathcal{X}_{\theta}\subset[-1,1]^{p}$. The value of $tr(\mathbf{L}^{\intercal}\mathbf{L})$ will decrease in such cases, compared to the case when the design space equals $[-1,1]^{p}$. The variance of the proposed estimator thus will increase in such cases. Nevertheless, the variance term will not be inflated to be arbitrarily large, as long as the design space is not too small. More discussion on the impact of the design space to the Inequality (8) is relegated to the Supplementary Material.

There are two essential assumptions in Theorem 3.1. One is that $s_{p}(\mathbf{L})>s_{1}(\mathbf{D})$, and the other is that the Taylor expansion is valid, i.e., when $s_{1}(\mathbf{D})$ is “small”. Although we will evaluate the quality of the proposed estimator empirically in the next section, a precise theoretical characterization of when these two assumptions are valid is currently not available. Here, we simply give an example such that $s_{1}(\mathbf{D})$ converges to zero as $n$ goes to infinity, in which case the desired Taylor expansion is apparently valid. The assumption $s_{p}(\mathbf{L})>s_{1}(\mathbf{D})$ is also satisfied in such a case, as $n$ goes to infinity, since the value of $s_{p}(\mathbf{L})$ is not relevant to $n$. Consider the case when the non-zero support of the population distribution is $[-1,1]^{p}$, i.e., the sample and the design points have the same domain. In such a case, the distance between each design point and its nearest neighbor converges to zero, as $n$ goes to infinity. As a result, each entry of the matrix $\mathbf{D}$ converges to zero, and thus $s_{1}(\mathbf{D})$ converges to zero as well, as $n$ goes to infinity. Consequently, the desired Taylor expansion is valid in such a case.

Soil functional properties refer to the properties related to a soil’s capacity to support essential ecosystem services, which include primary productivity, nutrient and water retention, and resistance to soil erosion [16]. The soil functional properties are thus important for planning sustainable agricultural intensification and natural resource management. To measure the soil functional properties in a target area, a natural paradigm is to first collect a sample of soil in this area, then analyze the sample using the technique of diffuse reflectance infrared spectroscopy [36]. Such a paradigm might be time-consuming or even impractical if the desired sample of soil from the target area is difficult to obtain. Predicting the soil functional properties is thus a measurement-constrained problem.

With the help of greater availability of Earth remote sensing data, the practitioners are provided new opportunities to predict soil functional properties at unsampled locations. One of the Earth remote sensing databases is provided by the Shuttle Radar Topography Mission (SRTM), which aims to generate the most complete high-resolution digital topographic database of Earth [13]. In this section, we consider the Africa Soil Property Prediction dataset, which contains the soil samples from 1157 different areas ($n=1157$). We aim to analyze the relationship between the sand content, one of the soil functional properties, and the five features ($p=5$) derived from the SRTM data. The features include compound topographic index calculated from SRTM elevation data (CTI), SRTM elevation data (ELEV), topographic Relief calculated from SRTM elevation data (RELI), mean annual precipitation of average long-term Tropical Rainfall Monitoring Mission data (TMAP), and modified Fournier index of average long-term Tropical Rainfall Monitoring Mission data (TMFI). We assume the data follow the model,

where the random errors $u_{i}$ are i.i.d. and follow a non-centered normal distribution $N(h(\bm{x}_{i}),\sigma^{2})$. Here, $\bm{x}_{i}=(1,CTI_{i},ELEV_{i},RELI_{i},TMAP_{i},TMFI_{i})^{\intercal}$ and $h(\cdot)$ represents a multivariate function that is unknown to the practitioner. The postulated model is thus a misspecified linear model. In our measurement-constrained setting, we assume the response vector is hidden unless explicitly requested. We then estimate the true coefficient of Model (11), i.e., $(\beta_{0},\beta_{1},\beta_{2},\beta_{3},\beta_{4},\beta_{5})^{\intercal}$, using subsampling methods.

The subsampling methods considered here are uniform subsampling (UNIF), basic leverage subsampling (BLEV), shrinkage leverage subsampling (SLEV) with parameter $\alpha=0.9$, unweighted-leverage subsampling (LEVUNW) [22, 23], information-based optimal subset selection (IBOSS) [48] and the proposed LowCon method. Table 1 summarizes the EMSEs for all six SLS estimators, and the best result in each row is in bold letter. We observe that the proposed LowCon method yields the best result in every row.

The second real-data example we consider is the Diamond Price Prediction dataset ¹¹1The dataset can be downloaded from https://www.kaggle.com/shivam2503/diamonds., which contains the prices and the features of around 54,000 diamonds. Of interest is to analyze the relationship between the price of the diamond, and three continuous features ($p$=3): weight of the diamond (caret), total depth percentage (depth), and width of top of diamond relative to widest point (table).

As the same setting used in Section 5.1, we assume the data follow a misspecified linear model,

Here, the random errors $u_{i}$ are i.i.d. and follow a non-centered normal distribution $N(h(\bm{x}_{i}),\sigma^{2})$, where $\bm{x}_{i}=(1,caret_{i},depth_{i},table_{i})^{\intercal}$, and $h(\cdot)$ is a multivariate function that is unknown to the practitioner. Note that the price of a diamond might be time-consuming or even impossible to obtain if the diamond has not been on the market yet. We thus assume the value of the response vector is hidden unless explicitly requested, and we estimate the true coefficient using subsampling methods.

Table 2 summarizes the EMSEs for all the subsample estimators, and the best result in each row is in bold letter. From Table 2, we observe that the proposed LowCon algorithm yields decent performance in all the cases and the best result in most of the cases.