# 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.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin19.5.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.10_1/lib/libopenblasp-r0.3.10.dylib
##
## locale:
##  en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
##  stats     graphics  grDevices datasets  utils     methods   base
##
## other attached packages:
##  ggplot2_3.3.2 CVXR_1.0-9
##
## loaded via a namespace (and not attached):
##   gmp_0.6-0        Rcpp_1.0.5       compiler_4.0.2   pillar_1.4.6
##   tools_4.0.2      digest_0.6.25    bit_1.1-15.2     evaluate_0.14
##   lifecycle_0.2.0  tibble_3.0.3     gtable_0.3.0     lattice_0.20-41
##  pkgconfig_2.0.3  rlang_0.4.7      Matrix_1.2-18    gurobi_9.0.3.1
##  Rglpk_0.6-4      yaml_2.2.1       blogdown_0.19    xfun_0.15
##  osqp_0.6.0.3     cccp_0.2-4       withr_2.2.0      dplyr_1.0.0
##  Rmpfr_0.8-1      stringr_1.4.0    knitr_1.28       generics_0.0.2
##  vctrs_0.3.2      tidyselect_1.1.0 bit64_0.9-7      grid_4.0.2
##  glue_1.4.1       R6_2.4.1         rmarkdown_2.3    bookdown_0.19
##  farver_2.0.3     purrr_0.3.4      magrittr_1.5     codetools_0.2-16
##  rcbc_0.1.0.9001  scales_1.1.1     htmltools_0.5.0  ellipsis_0.3.1
##  assertthat_0.2.1 colorspace_1.4-1 labeling_0.3     Rcplex_0.3-3
##  stringi_1.4.6    Rmosek_9.2.3     munsell_0.5.0    slam_0.1-47
##  crayon_1.3.4

R Markdown