Measurement Error in Regression
Classical Measurement Error Problems
Consider a linear regression model, \(y = X\beta + \epsilon\), where \(y\) is an \(n\)-dimensional response vector, \(X\) is an \(n \times p\) design matrix, and \(\epsilon\) is normally distributed error. With \(n>p\) (and \(X\) positive definite), unbiased regression coefficients are obtained as
\[ \hat{\beta} = (X^{T}X )^{-1} X^{T} y. \]
In many cases, however, the true covariates \(X\) are not directly observed. Instead, we have noisy measurements
\[ W = X + U, \]
where \(W\) is the \(n\times p\) measurement matrix, and \(U\) is the \(n\times p\) matrix of measurement errors. Assume for the moment that the \(n\) rows of \(U\), \(u_{i}\) are identically and independently distributed with covariance matrix \(\Sigma_{uu}\), i.e.,
\[ u_{i} \sim N(0, \Sigma_{uu}), ~ i = 1,\dots,n. \]
Since we do not have measurements of \(X\), using a classical linear regression model now yields coefficient estimates
\[ \hat{\beta}_{naive} = (W^{T}W )^{-1} W^{T} y, \] which is referred to as the naive estimate in the measurement error literature. Naive estimates are typically biased. For example, when all measurement errors are uncorrelated and each have variance \(\sigma_{u}\), i.e., \(\Sigma_{uu} = \sigma_{u} I_{p\times p}\), the expected value of \(\hat{\beta}\) is
\[ E(\hat{\beta}_{naive}) = \frac{\sigma_{x}^{2}}{\sigma_{u}^{2} + \sigma_{x}^2} \beta = \gamma \beta, \] where \(\gamma\) is the attenuation factor. Unbiased estimates for the case of linear regression can be obtained by minimizing the corrected loss function \[ L_{corr}(\beta) = \|y - W\beta\|^{2} - \beta^{T}\Sigma_{uu} \beta. \] This \(L_{corr}(\beta)\) is not always convex. If the Hessian \(W^{T}W - \Sigma_{uu}\) is positive definite, \(L_{corr}(\beta)\) is convex, and the estimates can be found using
\[ \hat{\beta}_{corr} = (W^{T}W - \Sigma_{uu} )^{-1} W^{T} y. \] Otherwise, iteration methods must be used to find a minimum of \(L(\beta)\).
An estimate \(\hat{\Sigma}_{uu}\) of \(\Sigma_{uu}\) can be typically obtained with replicate measurements.
Similar results hold for generalized linear models, e.g., logistic and Poisson regression: Measurement error typically leads to bias in the naive estimates, and correction methods can be used to obtained unbiased estimates. We refer to the book Carroll et al. (2006) for a thorough introduction to measurement error modeling.
Measurement Error in High-Dimensional Regression Problems
Now consider the same setup as in the previous section, but with high-dimensional data. That is, the number of covariates \(p\) is either on the same order as \(n\), or even larger than \(n\). In this case, even in the absence of measurement error, regularization is typically needed (see, e.g., Chapter 3 of Hastie, Tibshirani, and Friedman (2009)).
The lasso
The lasso (Tibshirani (1996)) performs
\(L1\)-regularization, which induces
sparsity. A popular R implementation can be found in the
glmnet package (Simon et al.
(2011)). For linear regression models, the lasso finds \(\hat{\beta}\) minimizing the loss \[
L(\beta) = \|y - X \beta\|^{2} + \lambda \|\beta\|_{1},
\] where \(\lambda\) is a
regularization parameter. The lasso can also be formulated as a
constrained optimization problem, \[
\text{minimize } ~\|y - X\beta\|^{2}~ \text{ subject to }~
\|\beta\|_{1} \leq R,
\] where \(R\) is in a
one-to-one relationship with \(\lambda\).
A lot of effort has been put into understanding the statistical properties of the lasso, both in the classical \(p<n\) case and in the high-dimensional \(p>n\) case, summarized in Bühlmann and Geer (2011). The results are typically either probabilistic upper bounds on or asymptotic limits for
- Prediction error \(E(\|y_{new} - X_{new}\hat{\beta}\|)\) where \((X_{new}, y_{new})\) are a new sample of data.
- Sum of squared estimation error \(\|\hat{\beta} - \beta\|_{2}^{2}\) or sum of absolute estimation error \(\|\hat{\beta} - \beta\|_{1}\)
- Covariate selection: If true coefficients in the index set \(J \subseteq \{1, \dots, p\}\) are nonzero, while the rest are zero, i.e., \(\beta_{J} \neq 0\) and \(\beta_{J^{C}} = 0\), under which conditions will the lasso recover the true set of relevant covariates while correctly setting the irrelevant covariates to zero?
The impact of measurement error in the lasso for linear models has been studied by Sorensen, Frigessi, and Thoresen (2015). The authors show that estimation of \(\hat{\beta}\) in the asymptotic limit suffers the same bias as for a multivariate linear model described above. Consistent covariate selection, on the other hand, requires stricter conditions than in the case without measurement error.
Using the glmnet package, we can illustrate the impact
of measurement error in a small experiment.
First we set up the parameters and generate random data. In order to
repeat the procedure, we create the function
create_example_data for doing this. For illustration
purposes, we set the true coefficient vector to \(\beta = (-2, -1, 0.5, 1, 2, 0, \dots,
0)^{T}\).
create_example_data <- function(n, p, s = 5, sdX = 1, sdU = 0.5,
sdEpsilon = 0.1, family = "gaussian") {
# Independent true covariates with mean zero and standard deviation sdX
X <- matrix(rnorm(n * p, sd = sdX), nrow = n, ncol = p)
# if Gaussian, response has standard deviation sdEpsilon, and zero intercept
# if binomial, response is binomial with mean (1 + exp(-X %*% beta))^(-1)
beta <- c(-2, -1, 0.5, 1, 2, rep(0, p - s))
if(family == "gaussian") {
# True coefficient vector has s non-zero elements and p-s zero elements
y <- X %*% beta + rnorm(n, sd = sdEpsilon)
} else if (family == "binomial") {
# Need an amplification in the binomial case
beta <- beta * 3
y <- rbinom(n, size = 1, prob = (1 + exp(-X %*% beta))**(-1))
}
# The measurements W have mean X and standard deviation sdU.
# We assume uncorrelated measurement errors
W <- X + matrix(rnorm(n * p, sd = sdU), nrow = n, ncol = p)
return(list(X = X, W = W, y = y, beta = beta, sigmaUU = diag(p) * sdU))
}We then call the function to get the data and put the result in the
list ll.
n <- 100
p <- 500
set.seed(1000)
ll <- create_example_data(n, p)Next we run the lasso with cross-validation on the true covariates
and on the noisy measurements. We pick the coefficient estimate at \(\lambda_{1se}\), i.e., on standard error on
the sparse side of the cross-validation minimum. This is the default of
the coef function for objects of class
cv.glmnet.
library(glmnet)
library(dplyr)
# Lasso with cross-validation on data without measurement error
fit1 <- cv.glmnet(ll$X, ll$y)
# Lasso with cross-validation on data with measurement error
fit2 <- cv.glmnet(ll$W, ll$y)
# Create a data frame with results ([-1] because we drop the intercept)
lassoEstimates <- tibble(
index = rep(1:p, times = 3),
beta = c(ll$beta, as.numeric(coef(fit1)[-1]), coef(fit2)[-1]),
label = c(rep("True values", p), rep("No measurement error", p), rep("Measurement error", p))
)By plotting all the estimated regression coefficients, we see that when the data are subject to measurement error, the number of false positives may increase. Note that this is not necessarily the case for all choices of parameters.
library(ggplot2)
ggplot(lassoEstimates, aes(x = index, y = beta, color = label)) +
geom_point() +
xlab("p") +
theme(legend.title=element_blank()) +
ggtitle("Measurement error leading to false positives")
We can also focus on the 5 parameters which are truly nonzero. In this case, we see that in the absence of measurement error, the lasso estimates values quite close to truth. With measurement error, on the other hand, the attenuation is quite clear: the estimates are biased toward zero.
library(tidyr)
estimatesOfNonzero <- lassoEstimates %>%
spread(key = label, value = beta) %>%
filter(`True values` != 0) %>%
gather(key = label, value = beta, -index)
ggplot(estimatesOfNonzero, aes(x = index, y = beta, color = label)) +
geom_point() +
xlab("p") +
theme(legend.title=element_blank()) +
ggtitle("Measurement error leading to attenuation")
The Dantzig Selector
The Dantzig selector Candes and Tao (2007) is closely related to the lasso. It is defined as the solution to the optimization problem \[ \text{minimize } ~ \|\beta\|_{1}, ~ \text{ subject to } ~ (1/n)\| X^{T} (y - X\beta)\|_{\infty} \leq \lambda, \] where \(\|\cdot\|_{\infty}\) denotes the maximum component norm. A generalized Dantzig selector for generalized linear models was introduced by James and Radchenko (2009). It is defined as the solution to the optimization problem \[ \text{minimize } ~ \|\beta\|_{1}, ~ \text{ subject to } ~ (1/n)\| X^{T} (y - \mu(X\beta))\|_{\infty} \leq \lambda, \] where \(\mu(\cdot) \in \mathbb{R}^{n}\) is the vector valued mean function of the generalized linear model. Examples include logistic regression with \(\mu(x) = (1+\exp(-x))^{-1}\) and Poisson regression with \(\mu(x) = \exp(x)\). Using an iterative reweighing algorithm, the generalized Dantzig selector can be solved as a sequence of linear optimization problems of the same form as the Dantzig selector for linear models.
Generalized Dantzig Selector
To our knowledge, no R package to date contains an implementation of
the Generalized Dantzig Selector (GDS). Since the Generalized Matrix
Uncertainty Selector (GMUS) presented later in this vignette is a
generalization of the GDS, we have included a function for computing the
GDS, called gds.
# Number of samples
n <- 1000
# Number of covariates
p <- 50
# Create example data
ll <- create_example_data(n, p, family = "binomial")Arguments to gds. Only \(X\) and \(y\) are required.
args(gds)
#> function (X, y, lambda = NULL, family = "gaussian", weights = NULL)
#> NULLFit the generalized Dantzig selector on the data.
# Fit the Generalized Dantzig Selector
gds_estimate <- gds(ll$X, ll$y, family = "binomial")The result is a list of class gds.
class(gds_estimate)
#> [1] "gds"The list contains the intercept (which has not been penalized), the
regression coefficients beta, the family, the value of
lambda and the number of nonzero components of
beta. By default, gds uses the value \(\lambda_{min}\) corresponding to the
minimum cross-validated loss for the lasso, using
cv.glmnet.
str(gds_estimate)
#> List of 5
#> $ intercept : num 0.189
#> $ beta : num [1:50, 1] -3.971 -2.058 0.942 1.839 4.04 ...
#> $ family : chr "binomial"
#> $ lambda : num 0.00378
#> $ num_non_zero: num 26
#> - attr(*, "class")= chr "gds"At the moment, only a single value of lambda is accepted
by GDS. We will fix this in the future, as well as adding a cross
validation function.
