# 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.4.0 (2024-04-24)
## Platform: x86_64-apple-darwin23.4.0
## Running under: macOS Sonoma 14.5
##
## Matrix products: default
## BLAS:   /usr/local/Cellar/openblas/0.3.27/lib/libopenblasp-r0.3.27.dylib
## LAPACK: /usr/local/Cellar/r/4.4.0_1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base
##
## other attached packages:
## [1] ggplot2_3.5.1 CVXR_1.0-13
##
## loaded via a namespace (and not attached):
##  [1] gmp_0.7-4         sass_0.4.9        utf8_1.2.4        generics_0.1.3
##  [5] slam_0.1-50       blogdown_1.19     lattice_0.22-6    digest_0.6.35
##  [9] magrittr_2.0.3    evaluate_0.23     grid_4.4.0        bookdown_0.39
## [13] fastmap_1.2.0     jsonlite_1.8.8    Matrix_1.7-0      Rmosek_10.2.0
## [17] fansi_1.0.6       scales_1.3.0      codetools_0.2-20  jquerylib_0.1.4
## [21] cli_3.6.2         Rmpfr_0.9-5       rlang_1.1.3       Rglpk_0.6-5.1
## [25] bit64_4.0.5       munsell_0.5.1     withr_3.0.0       cachem_1.1.0
## [29] yaml_2.3.8        tools_4.4.0       osqp_0.6.3.2      Rcplex_0.3-6
## [33] rcbc_0.1.0.9001   dplyr_1.1.4       colorspace_2.1-0  gurobi_11.0-0
## [37] assertthat_0.2.1  vctrs_0.6.5       R6_2.5.1          lifecycle_1.0.4
## [41] bit_4.0.5         pkgconfig_2.0.3   cccp_0.3-1        pillar_1.9.0
## [45] bslib_0.7.0       gtable_0.3.5      glue_1.7.0        Rcpp_1.0.12
## [49] highr_0.11        xfun_0.44         tibble_3.2.1      tidyselect_1.2.1
## [53] knitr_1.47        farver_2.1.2      htmltools_0.5.8.1 labeling_0.4.3
## [57] rmarkdown_2.27    compiler_4.4.0

R Markdown