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.
suppressWarnings(suppressMessages(library(CVXR)))
## 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 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/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] methods stats graphics grDevices datasets utils base
##
## other attached packages:
## [1] CVXR_0.95 ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] gmp_0.5-13.1 Rcpp_0.12.15 pillar_1.1.0
## [4] compiler_3.4.3 plyr_1.8.4 R.methodsS3_1.7.1
## [7] R.utils_2.6.0 tools_3.4.3 digest_0.6.15
## [10] bit_1.1-12 evaluate_0.10.1 tibble_1.4.2
## [13] gtable_0.2.0 lattice_0.20-35 rlang_0.2.0
## [16] Matrix_1.2-12 yaml_2.1.16 blogdown_0.5.4
## [19] xfun_0.1 Rmpfr_0.7-0 ECOSolveR_0.4
## [22] stringr_1.3.0 knitr_1.20 rprojroot_1.3-2
## [25] bit64_0.9-7 grid_3.4.3 R6_2.2.2
## [28] rmarkdown_1.8.10 bookdown_0.7 magrittr_1.5
## [31] backports_1.1.2 scales_0.5.0 htmltools_0.3.6
## [34] scs_1.1-1 colorspace_1.3-2 labeling_0.3
## [37] stringi_1.1.6 lazyeval_0.2.1 munsell_0.4.3
## [40] R.oo_1.21.0