Censored Regression
Introduction
Data collected from an experimental study is sometimes censored, so that only partial information is known about a subset of observations. For instance, when measuring the lifespan of mice, we may find a number of subjects live beyond the duration of the project. Thus, all we know is the lower bound on their lifespan. This right censoring can be incorporated into a regression model via convex optimization.
Suppose that only \(K\) of our observations \((x_i,y_i)\) are fully observed, and the remaining are censored such that we observe \(x_i\), but only know \(y_i \geq D\) for \(i = K+1,\ldots,m\) and some constant \(D \in {\mathbf R}\). We can build an OLS model using the uncensored data, restricting the fitted values \(\hat y_i = x_i^T\beta\) to lie above \(D\) for the censored observations:
\[ \begin{array}{ll} \underset{\beta}{\mbox{minimize}} & \sum_{i=1}^K (y_i - x_i^T\beta)^2 \\ \mbox{subject to} & x_i^T\beta \geq D, \quad i = K+1,\ldots,m. \end{array} \]
This avoids the bias introduced by standard OLS, while still utilizing
all of the data points in the regression. The constraint requires only
one more line in CVXR.
Example
We will generate synthetic data for this example, censoring observations beyond a value \(D\).
## Problem data
n <- 30
M <- 50
K <- 200
set.seed(n * M * K)
X <- matrix(stats::rnorm(K * n), nrow = K, ncol = n)
beta_true <- matrix(stats::rnorm(n), nrow = n, ncol = 1)
y <- X %*% beta_true + 0.3 * sqrt(n) * stats::rnorm(K)
## Order variables based on y
idx <- order(y, decreasing = FALSE)
y_ordered <- y[idx]
X_ordered <- X[idx, ]
## Find cutoff and censor
D <- (y_ordered[M] + y_ordered[M + 1]) / 2
censored <- (y_ordered > D)
y_censored <- pmin(y_ordered, D)We now fit regular OLS, OLS using just the censored data and finally the censored regression.
## Regular OLS
beta <- Variable(n)
obj <- sum((y_censored - X_ordered %*% beta)^2)
prob <- Problem(Minimize(obj))
result <- solve(prob)
beta_ols <- result$getValue(beta)
## OLS using uncensored data
obj <- sum((y_censored[1:M] - X_ordered[1:M,] %*% beta)^2)
prob <- Problem(Minimize(obj))
result <- solve(prob)
beta_unc <- result$getValue(beta)
## Censored regression
obj <- sum((y_censored[1:M] - X_ordered[1:M,] %*% beta)^2)
constr <- list(X_ordered[(M+1):K,] %*% beta >= D)
prob <- Problem(Minimize(obj), constr)
result <- solve(prob)
beta_cens <- result$getValue(beta)Here’s are some plots comparing the results. The blue diamond points are estimates.
plot_results <- function(beta_res, title) {
d <- data.frame(index = seq_len(K),
y_ordered = y_ordered,
y_fit = as.numeric(X_ordered %*% beta_res),
censored = as.factor(censored))
ggplot(data = d) +
geom_point(mapping = aes(x = index, y = y_ordered, color = censored)) +
geom_point(mapping = aes(x = index, y = y_fit), color = "blue", shape = 23) +
geom_abline(intercept = D, slope = 0, lty = "dashed") +
labs(x = "Observations", y = "y") +
ggtitle(title)
}plot_results(beta_ols, "Regular OLS.")
plot_results(beta_unc, "OLS using uncensored data.")
plot_results(beta_cens, "Censored Regression.")
Session Info
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin21.6.0 (64-bit)
## Running under: macOS Ventura 13.0
##
## Matrix products: default
## BLAS: /usr/local/Cellar/openblas/0.3.21/lib/libopenblasp-r0.3.21.dylib
## LAPACK: /usr/local/Cellar/r/4.2.1_4/lib/R/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] ggplot2_3.3.6 CVXR_1.0-11
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 xfun_0.34 bslib_0.4.0 slam_0.1-50
## [5] lattice_0.20-45 Rmosek_10.0.25 colorspace_2.0-3 vctrs_0.5.0
## [9] generics_0.1.3 htmltools_0.5.3 yaml_2.3.6 gmp_0.6-6
## [13] utf8_1.2.2 rlang_1.0.6 jquerylib_0.1.4 pillar_1.8.1
## [17] glue_1.6.2 Rmpfr_0.8-9 withr_2.5.0 DBI_1.1.3
## [21] Rcplex_0.3-5 bit64_4.0.5 lifecycle_1.0.3 stringr_1.4.1
## [25] munsell_0.5.0 blogdown_1.13 gtable_0.3.1 gurobi_9.5-2
## [29] codetools_0.2-18 evaluate_0.17 labeling_0.4.2 knitr_1.40
## [33] fastmap_1.1.0 fansi_1.0.3 cccp_0.2-9 highr_0.9
## [37] Rcpp_1.0.9 scales_1.2.1 cachem_1.0.6 osqp_0.6.0.5
## [41] jsonlite_1.8.3 farver_2.1.1 bit_4.0.4 digest_0.6.30
## [45] stringi_1.7.8 bookdown_0.29 dplyr_1.0.10 Rglpk_0.6-4
## [49] grid_4.2.1 cli_3.4.1 tools_4.2.1 magrittr_2.0.3
## [53] sass_0.4.2 tibble_3.1.8 pkgconfig_2.0.3 rcbc_0.1.0.9001
## [57] Matrix_1.5-1 assertthat_0.2.1 rmarkdown_2.17 R6_2.5.1
## [61] compiler_4.2.1