Classical Measurement Error Problems

Consider a linear regression model, y = Xβ + ϵ, where y is an n-dimensional response vector, X is an n × p design matrix, and ϵ is normally distributed error. With n > p (and X positive definite), unbiased regression coefficients are obtained as

β̂ = (X^TX)⁻¹X^Ty.

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 × p measurement matrix, and U is the n × 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 Σ_uu, i.e.,

u_i ∼ N(0, Σ_uu), i = 1, …, n.

Since we do not have measurements of X, using a classical linear regression model now yields coefficient estimates

β̂_naive = (W^TW)⁻¹W^Ty, 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 σ_u, i.e., Σ_uu = σ_uI_p × p, the expected value of β̂ is

$$ E(\hat{\beta}_{naive}) = \frac{\sigma_{x}^{2}}{\sigma_{u}^{2} + \sigma_{x}^2} \beta = \gamma \beta, $$ where γ is the attenuation factor. Unbiased estimates for the case of linear regression can be obtained by minimizing the corrected loss function L_corr(β) = ∥y − Wβ∥² − β^TΣ_uuβ. This L_corr(β) is not always convex. If the Hessian W^TW − Σ_uu is positive definite, L_corr(β) is convex, and the estimates can be found using

β̂_corr = (W^TW − Σ_uu)⁻¹W^Ty. Otherwise, iteration methods must be used to find a minimum of L(β).

An estimate Σ̂_uu of Σ_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)).

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))  
}

n <- 100
p <- 500
set.seed(1000)
ll <- create_example_data(n, p)

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))
  )

library(ggplot2)
theme_set(theme_bw())
theme_update(
  panel.grid = element_blank()
)
ggplot(lassoEstimates, aes(x = index, y = beta, color = label)) +
  geom_point() +
  xlab("p") +
  scale_color_brewer(type = "qualitative", palette = "Paired") +
  theme(legend.title=element_blank()) + 
  ggtitle("Measurement error leading to false positives")

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") +
  scale_color_brewer(type = "qualitative", palette = "Paired") +
  theme(legend.title=element_blank()) + 
  ggtitle("Measurement error leading to attenuation")

# Number of samples
n <- 1000
# Number of covariates
p <- 50
# Create example data
ll <- create_example_data(n, p, family = "binomial")

args(gds)
#> function (X, y, lambda = NULL, family = "gaussian", weights = NULL) 
#> NULL

# Fit the Generalized Dantzig Selector
gds_estimate <- gds(ll$X, ll$y, family = "binomial")

class(gds_estimate)
#> [1] "gds"

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"

Corrected Lasso for Linear Regression

When an estimate of the measurement error covariance matrix Σ_uu is available, a corrected lasso can be defined as minimizing the loss minimize  L(β) = ∥y − Wβ∥² − β^TΣ_uuβ + λ∥β∥₁  subject to ∥β∥₁ ≤ R. Because we subtract the positive semidefinite matrix Σ_uu from the convex function ∥y − Wβ∥², this corrected lasso may be non-convex. In fact, when p > n is is always non-convex. Hence, in order to avoid non-trivial solutions, we must constrain the solution to lie in some L1-ball with radius R. Since this problem involves two regularization parameters, λ and R, it is more convenient to use the constrained version of the lasso

minimize  L(β) = ∥y − Wβ∥² − β^TΣ_uuβ  subject to ∥β∥₁ ≤ R. Loh and Wainwright (2012) analyze the properties of these two closely related versions of the lasso. They show that the bounds for the estimates of ∥β̂ − β∥₂² and ∥β̂ − β∥₁ are of the same order as for the standard lasso without measurement error, where β̂ is the global minimum of the optimization problem. More remarkably, they show that despite non-convexity, under mild assumptions a projected gradient descent algorithm will converge with high probability to a local optimum which is very close to the global optimum.

Sorensen, Frigessi, and Thoresen (2015) analyzed the covariate selection properties of the same model, and similarly showed that results very similar to those for covariate selection with the standard lasso in the absence of measurement, also hold for this corrected lasso.

set.seed(1000)
# Generate example data
ll <- create_example_data(n, p)
# Fit the corrected lasso
corrected_fit <- corrected_lasso(W = ll$W, y = ll$y, sigmaUU = ll$sigmaUU)

# Class of the object
class(corrected_fit)
#> [1] "corrected_lasso"
# The coef() method prints the number of nonzero estimates as a function of the radius
coef(corrected_fit)
#> Number of nonzero coefficient estimates
#> as a function of regularization parameter (radius):
#>       radius nonzeros
#>   0.01006214        1
#>   0.53911872        2
#>   1.06817530        2
#>   1.59723188        2
#>   2.12628846        2
#>   2.65534505        2
#>   3.18440163        4
#>   3.71345821        4
#>   4.24251479        4
#>   4.77157137        4
#>   5.30062795        4
#>   5.82968454        5
#>   6.35874112        5
#>   6.88779770        5
#>   7.41685428        6
#>   7.94591086        8
#>   8.47496744       14
#>   9.00402403       17
#>   9.53308061       22
#>  10.06213719       24

args(corrected_lasso)
#> function (W, y, sigmaUU, family = c("gaussian", "binomial", "poisson"), 
#>     radii = NULL, no_radii = NULL, alpha = 0.1, maxits = 5000, 
#>     tol = 1e-12) 
#> NULL

plot(corrected_fit)

plot(corrected_fit, type = "path")

Corrected Lasso for Generalized Linear Models

Sorensen, Frigessi, and Thoresen (2015) show how we can also correct for measurement error in the lasso for generalized linear models, using the conditional score method introduced by Stefanski and Carroll (1987), combined with the projected gradient descent algorithm used by Loh and Wainwright (2012) for the corrected lasso for linear models.

set.seed(323)
n <- 100
p <- 50
ll <- create_example_data(n, p, sdU = 0.2, family = "binomial")

corrected_fit <- corrected_lasso(ll$W, ll$y, ll$sigmaUU, family = "binomial")

plot(corrected_fit)

plot(corrected_fit, type = "path")

Model Tuning

The corrected lasso for linear regression has a clearly defined loss function, L(β) = ∥y − Wβ∥₂² − β^TΣ_uuβ. In order to find the optimal value of the regularization parameter R, we can hence use cross-validation to minimize L(β). This is implemented in the function cv_corrected_lasso. We illustrate its use below.

For generalized linear models, we are using the conditional score approach, which does not yield a well defined loss function. Hence, cross-validation is not straightforward in these cases.

set.seed(1000)
# Generate example data
ll <- create_example_data(n, p)
# Run lasso with cross-validation
cv_corrected_fit <- cv_corrected_lasso(W = ll$W, y = ll$y, sigmaUU = ll$sigmaUU)

The result is a list of class cv_corrected_lasso.

class(cv_corrected_fit)
#> [1] "cv_corrected_lasso"

The print below shows all the elements of the resulting list.

str(cv_corrected_fit)
#> List of 6
#>  $ cv        :'data.frame':  100 obs. of  5 variables:
#>   ..$ radii    : num [1:100] 0.0119 0.1324 0.2528 0.3732 0.4937 ...
#>   ..$ mean_loss: num [1:100] 8.74 8.48 8.22 7.97 7.7 ...
#>   ..$ sd_loss  : num [1:100] 9.07 8.74 8.42 8.1 7.8 ...
#>   ..$ upper_1se: num [1:100] 11.6 11.2 10.9 10.5 10.2 ...
#>   ..$ lower_1se: num [1:100] 5.88 5.71 5.56 5.41 5.23 ...
#>  $ radius_min: num 9.28
#>  $ loss_min  : num -4.11
#>  $ radius_1se: num 7.48
#>  $ loss_1se  : num -3.79
#>  $ family    : chr "gaussian"
#>  - attr(*, "class")= chr "cv_corrected_lasso"

The element cv contains all the details of the cross-validation runs.

• radii contains all the constraint radii R that were used.
• mean_loss contains the mean of the loss function at the given radius over all folds.
• sd_loss gives the standard deviation of the loss at the given radius.
• upper_1se and lower_1se contain the upper and lower standard error of mean_loss at the given radius.

Next, loss_min gives the minimum of the mean loss, and radius_min is the corresponding radius. loss_1se gives the smallest loss within one standard error of loss_min, and radius_1se is a smallest radius given a loss less than or equal to loss_1se.

The result of performing cross-validation can be illustrated using the plot function. It shows the cross-validated loss over the grid of radii. R_min and R_1se are shown with labels, and the corresponding loss as horizontal red lines.

plot(cv_corrected_fit)

Having used cv_corrected_lasso to find the right value of the constraint parameter R, we can use corrected_lasso on this single value. The snippet below shows how to compute the solution using R_1se.

corrected_fit <- corrected_lasso(ll$W, ll$y, ll$sigmaUU, radii = cv_corrected_fit$radius_1se)

The final parameter estimates can be found in corrected_fit$betaCorr.

str(corrected_fit)
#> List of 3
#>  $ betaCorr: num [1:50, 1] -2.8 -1.2 0 1.28 1.76 ...
#>  $ radii   : num 7.48
#>  $ family  : chr "gaussian"
#>  - attr(*, "class")= chr "corrected_lasso"

Generalized Matrix Uncertainty Selector

The Generalized Matrix Uncertainty Selector (GMUS) is an extension of the MUS to generalized linear models, introduced by Sorensen et al. (2018). It is defined as the solution to the optimization problem

$$ \text{minimize } ~ \|\beta\|_{1}, ~ \text{ subject to } ~ (1/n)\|W^{T} (y - \mu(W \beta))\|_{\infty} \leq \lambda + \sum_{r=1}^{R}\delta^{r}(r! \sqrt{n})^{-1} \|\beta\|_{1}^{r} \|\mu^ {(r)}(W\beta)\|_{2}, $$ where μ(⋅) is the vector valued mean function of the generalized linear model and μ^(r)(⋅) is its rth derivative. R is a parameter controlling the number of Taylor expansion terms which are included. When R → ∞, the true solution is a member of the feasible set given that the bounds (1/n)∥W^Tϵ∥_∞ ≤ λ and ∥U∥_∞ ≤ δ hold. In practice, we set R = 1 for computational reasons. When R = 1, the GMUS can be solved using a sequence of linear programming problems on the same form as the MUS, with an iterative reweighing algorithm.

The gmus function in the hdme package implements the GMUS for R = 1, which is defined as minimize  ∥β∥₁,  subject to  (1/n)∥W^T(y − μ(Wβ))∥_∞ ≤ λ + δn^−(1/2)∥β∥₁∥μ⁽¹⁾(Wβ)∥₂. As the MUS, the GMUS does not require an estimate of Σ_uu.

We illustrate gmus using logistic regression. First, we generate some sample data.

set.seed(323)
n <- 100
p <- 50
ll <- create_example_data(n, p, sdU = 0.2, family = "binomial")
gmus_fit <- gmus(ll$W, ll$y, family = "binomial")

The returned object is again a list with class gmus.

class(gmus_fit)
#> [1] "gmus"
str(gmus_fit)
#> List of 6
#>  $ intercept   : num [1:26] -0.587 -0.322 -0.257 -0.225 -0.204 ...
#>  $ beta        : num [1:50, 1:26] -3.223 -1.406 0.885 1.696 2.729 ...
#>  $ family      : chr "binomial"
#>  $ delta       : num [1:26] 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 ...
#>  $ lambda      : num 0.00926
#>  $ num_non_zero: num [1:26] 29 14 10 8 7 5 5 5 5 5 ...
#>  - attr(*, "class")= chr "gmus"

The model fitting works the same way as for the MUS. A λ value is found by running the lasso with the appropriate link function, and taking the value yielding minimum cross-validation loss. A range of δ values are tried. We can call the plot function to study the behavior.

plot(gmus_fit)

Again we see that the number of nonzero coefficients decreases in δ. Following the elbow rule, a reasonable choice for δ might here be 0.1. To obtain a final estimate, we might therefore run gmus again with this value.

gmus_fit <- gmus(ll$W, ll$y, delta = 0.1, family = "binomial")

The plot will now show the estimated coefficients.

plot(gmus_fit)