Superset model problem

Regression analysis is a statistical method widely used in many research areas. It is often specified as the normal linear model, where coefficients are linear and the error term follows the normal distribution to simplify the analysis. This specification aims to approximate the state of nature and is often useful in prediction as well as discussing the causality.

Among its specifications, variable selection is a central issue in the regression analysis. It is important to select an appropriate set of explanatory variables partly because of the cost of collecting variables. Many methods are proposed for the variable selection problem.

In relation to the variable selection problem, this paper focuses on the superset model problem where the linear regression model selects a larger set of variables than the state of nature does (see the example provided in Section 2). The linear regression model tends to choose smaller set of variables when the state of nature is linear in variables due to its least squares loss. However, the nonlinear relationship is more likely and this case may lead to select a larger set of variables, which will be a deficiency of the linear regression model in terms of the data collection cost.

To evaluate the superset model problem, this paper utilizes the Bayesian approach, which provides a measure of uncertainty in the form of the posterior probability, and proposes an alternative model that will select the true set of variables when the sample size is large.

This paper is organized as follows. Section 2 describes the superset model problem by providing an example. Bayes’ Theorem is adopted to evaluate the problem in Section 3 and two regression models for specification is explained in Section 4. Section 6 illustrates the proposed method and discusses its robustness.

Suppose the continuous response $Y$ is associated with the set of explanatory variables $\bm{x}_{T}$. We are interested in its mean response conditional on $\bm{x}_{T}$. We often assume it to be linear in practice, although it is more likely to be nonlinear in reality. To this end, a regression model is specified as

where $\epsilon$ is an additive error term with mean zero. The functional form $\phi(\cdot)$, the distribution of the error term, and the true set of explanatory variables $\bm{x}_{T}$ are all unknown. With this model, we make statistical inferences about the conditional mean response by estimating $\phi(\cdot)$ and the set of explanatory variables. Among problems about how this regression model should be specified, the variable selection problem focuses on the set of explanatory variables, based on the dataset $\{y_{i},\bm{x}_{i}\}_{i=1}^{n}$.

When the set of explanatory variables is known, the best fit is $E(Y\mid\bm{x}_{T})$ as an estimator of $\phi(\bm{x}_{T})$ when we use the squared loss. The linear regression model assumes $E(Y\mid\bm{x}_{T})=\bm{x}_{T}^{\prime}\bm{\beta}$, where $\bm{\beta}$ is referred to as the regression coefficient vector. However, in general, $E(Y\mid\bm{x}_{T})\neq\bm{x}_{T}^{\prime}\bm{\beta}$, contrary to the linearity assumption.

For example, suppose

The linear regression with $(x_{T},x_{U})$, where $x_{U}=x_{T}^{2}$, is better than the one with $x_{T}$, even though the latter selects the true explanatory variable. This is an example of the superset model problem. On the other hand, when $x_{U}^{\prime}$ is independent of $x_{T}$, $x_{U}^{\prime}$ should not be included in the regression to improve the fit.

Above example suggests that the knowledge about association among variables is helpful to examine the superset model problem, and hence the variable selection problem. One approach is to estimate the conditional expectation without linearity and compare it with the one implied by the normal linear model. If they are different and the latter contains more explanatory variables, there exists the superset model problem. Because the dataset at hand is limited, it is difficult to determine whether the superset model problem exists or not. Rather, it is evaluated in a probability form, which is explained in the next section.

Suppose $\mathcal{M}^{\ast}$ is the set of explanatory variables in the state of nature, which is $\bm{x}_{T}$ when vectorized, and is known for a moment. The current dataset $\{y_{i},\bm{x}_{i}\}_{i=1}^{n}$ is generated from this state of nature independently for each observation $i$ and is observed. Depending on a context, the normal linear model with $\mathcal{M}^{\ast}$ explanatory variables would be a choice if it approximates the state of nature well. In this case, we do not have the superset model problem. On the other hand, a normal linear model with a set of explanatory variables indexed by $\mathcal{M}(\neq\mathcal{M}^{\ast})$ is chosen independent of the state of nature in terms of, say, prediction, where the superset model problem arises when $\mathcal{M}\supset\mathcal{M}^{\ast}$.

However, the state of nature is usually unknown and is inferred from the dataset. The uncertainty from inference is evaluated by the posterior probability over possible subsets of explanatory variables. This paper approximates it by assuming a flexible model (see Model (5) in Section 4), which is denoted by $H_{0}$. Then, this posterior probability is calculated via Bayes’ theorem, which is given by

The numerator is the cross product of the marginal likelihood and the prior belief about the set of explanatory variables. When the latter is uniform (which is assumed in the following empirical illustration), the posterior probability is proportional to the marginal likelihood under $H_{0}$.

When uncertainty from inference about the normal linear model is evaluated from its posterior probability as well, the overall superset model probability is calculated as

where $H_{1}$ denotes the normal linear model (see Model (4) in Section 4).

We note that the above expression is general enough to include variables (the response and explanatory variables) that are continuous or discrete. The next section specifies two regression models $H_{0}$ and $H_{1}$, where the response is assumed to be continuous for simplicity.

First, the linear regression model $H_{1}$ is specified as

where each of explanatory variables is standardized without loss of generality. To estimate model parameters $(\bm{\beta},\lambda^{2})$, the hyper-$g$ prior is assumed. Then, the marginal likelihood is analytically tractable (see Miyawaki and MacEachern (na) for example).

Second, the model $H_{1}$ that is alternative to the normal linear regression model is specified as

given $\bm{x}_{i}=\bm{x}$. The normal error assumption is made because we have no other knowledge on it. Further, it makes the conditional mean estimation simpler, in terms of the number of parameters as well as the computational burden. The main purpose of this semiparametric model is to estimate conditional means of $Y$ in a flexible manner, and to captures the association between $Y$ and $\bm{x}$ in the state of nature as the sample size increases, which can be viewed as as an extreme of the local constant estimation (see, e.g., Fan and Gijbels (2003) for the local estimation).

When the dataset is fixed, the covariate space becomes sparse as its dimension gets larger. Then, the marginal likelihood (hence, the superset model probability) under this alternative model is strongly dependent on the prior specification. To mitigate this influence, this paper takes the $m$-fold cross-validation approach, which is described in details below.

The dataset is randomly divided into $m$ groups. One of them is used as the test set, while the remainings belong to the training set. Explanatory variables in the training set are standardized, and those in the test set are standardized by the mean and standard deviation of those in the training set. Let $\mathcal{D}_{0}$ and $\mathcal{D}_{1}$ be the sets of identification numbers of observations which belongs to the training and test sets. More precisely, $\mathcal{D}_{0}=\{i\mid\text{the $i$-th observation is in the training set},i=1,\dots,n\}$ and $\mathcal{D}_{1}=\{i\mid\text{the $i$-th observation is in the test set},i=1,\dots,n\}$. Given a choice of $\mathcal{D}_{0}$ and $\mathcal{D}_{1}$, we construct the prior and conditional marginal likelihood in the following manner.

The prior for $\theta_{x}$ in the model (5) given $\bm{x}_{i}=\bm{x}$ and $i\in\mathcal{D}_{1}$ is assumed as

When $\mathcal{A}$ is a set, $|\mathcal{A}|$ is the number of elements in the set. This prior is constructed from classical OLS estimates of mean and standard deviation of $Y_{i}$ at $\bm{x}_{i}=\bm{x}$. By using the prior that is obtained from the linear model and using the model that focuses on the local observation, we are able to combine local and global information. It is possible to use other estimates such as corresponding normal linear regression estimates under the hyper-g prior. However, to keep the methodology as simple as possible, we take the above prior specification.

Then, we are able to derive the marginal likelihood conditional on the nuisance parameter $\sigma_{x}^{2}$ for each $\bm{x}_{i}=\bm{x}$ and $i\in\mathcal{D}_{1}$. Let $\mathcal{M}_{x}=\{i\mid\bm{x}_{i}=\bm{x},i\in\mathcal{D}_{1}\}$ and $n_{x}=|\mathcal{M}_{x}|$. Then, this conditional marginal likelihood is given by

The full Bayes analysis specifies a prior on the nuisance parameter $\sigma_{x}^{2}$ as well. However, because the data are sparse at $\bm{x}$, how we specify it affects the (unconditional) marginal likelihood much. To mitigate this problem, this paper takes the empirical Bayes approach. The marginal likelihood for each $\bm{x}_{i}=\bm{x}$ is the conditional marginal likelihood $m^{\ast}(\{Y_{i}\}_{i\in\mathcal{M}_{x}}\mid\{\bm{x}_{i}\}_{i\in\mathcal{M}_{x}},\sigma_{x}^{2},\{y_{i},\bm{x}_{i}\}_{i\in\mathcal{D}_{0}})$ maximized over $\sigma_{x}^{2}$. More precisely,

See the next section for this maximization in details. By multiplying it over all distinct $\bm{x}$ in $\mathcal{D}_{1}$ and taking the geometric mean, we have the marginal likelihood for $\mathcal{D}_{1}$ per one observation, which is given by

The geometric mean is to take care of different sample sizes in different test sets.

Finally, we repeat above process until all $m$ groups are used as the test set and calculate above marginal likelihood for each test group selection. After averaging $m$ marginal likelihoods, we raise it to the power of $n$ to obtain the final marginal likelihood estimate for the model (5).

The robustness of this approach is of interest. The approach will be more useful if we know the upper and lower bounds of the superset model probability under $H_{0}$ when its specification changes. Two points are discussed regarding robustness.

First, we consider the robustness to the number of folds in the cross-validation. In the following empirical analysis in Section 6, we use the 10-fold cross-validation. When the number of folds changes from 2 to 15, we see the resulting probability does not change much with the diabetes dataset (see Figure 1).

Second, it is important to check the robustness to the prior. One approach would be to use the $\epsilon$-contamination class prior (see Section 4.7.4 of Berger (1985)), and to show the sensitivity of the marginal likelihood, which will be our future analysis.

Letting

the local marginal likelihood function is characterized by Theorem 1.

See Appendix A for its proof. Table 1 summarizes the result.

By this theorem, we are able to set a value of $\sigma_{x}^{2}$ to maximize the marginal likelihood, instead of placing a prior on it.

The diabetes data (see Efron et al. (2004)) are used to illustrate our method. This dataset contains 442 observations. For the analysis below, we use the logarithm of the diabetes progression measure as the response and use remaining 10 variables are used as exlanatory variables. The proposed method is applied, and the superset model probability for this dataset with 10-fold cross-validation is estimated to be 22.37%.

As discussed by MacEachern and Miyawaki (2022), the diabetes dataset seems to be collected from least two different sources. In particular, the precision of two explanatory variables (the blood pressure and fourth blood serum measurement) consists of a mix of finer and coarser observations. When the data are divided into two groups by this precision, MacEachern and Miyawaki (2022) suggests these two datasets have different characteristics.

This conclusion is also confirmed in terms of the superset model probability. When the dataset for observations with finer variables is used, it is 22.19%. When, on the other hand, that for observations with coarser variables is used, it is 10.72%. The superset model problem is more likely to occur with the former dataset than the latter one, which would be due to the difference in characteristics of these two datasets.