Simple unbiased mean estimator for censored data

Question

Suppose I have an i.i.d. random sample of size $n$ s.t. $X_i \sim \mathcal{N}(\mu, \sigma^2)$ for all $i$. But suppose my observed sample is left censored, such that any $x_i$ observation is replaced with some constant value $U$ for all $x_i \leq U$. So $X_i$ is left-censored at $U$. I am interesting in estimating $\bar{x}$, the sample mean of the $x_i$.

I followed the Stan user manual to estimate the mean of a sample like this by treating the observed values as data, and the censored values as parameters to be estimated. The Stan code that I wrote was able to correctly estimate the mean with a reasonable amount of precision. But unfortunately this kind of Bayesian estimator with Stan is kind of cumbersome to use and implement. Is there a frequentist/MLE way to estimate the mean? I would want such an estimator to be unbiased and consistent. I'm specifically talking about the left-censored case, but I guess an MLE estimator for the right- or interval-censored case would be similar.

Answer 1

In the case of censored data, I would strongly suggest not caring about the mean so much but rather consider quantiles, i.e. the survival function if this ok in terms of the research question being asked. Outliers can have a strong influence on the mean and are more likely to be censored. So if your normal assumption is incorrect, you can add a lot of bias to the estimate of the mean. For example, if you assume your data is normal, but it's really a t-distribution, the left censoring will introduce a right-side bias to your estimate because on average, you'll assume your left censored values are higher than the really are.

However, estimators like the Kaplan Meier estimator will be a consistent estimator for the median and other quantiles. As long as your outcome is independent of your censoring mechanism, this will be consistent and as long as the censoring isn't too extreme, the estimate of the median should be minimal bias.

In your case of a constant left censoring at $U$, this will be dependent on whether $U < \mu$, in which case you case you can just use the the median directly, ignoring censoring. You can estimate whether $U < \mu$ by simply checking if less than half the values are censored. On the other hand, if you have constant censoring and more than half of your data is censored, you don't really have any hope of getting a reliable estimate of the median (or mean!) without extremely strong assumptions of the distribution which you can't really check!

To circle back to your question, if you still want to estimate the mean and if you are willing to make the strong assumption that your data is normal, the median will be equal to the mean, so you can just use your standard estimator of the median given above. But this is leaning for on assumptions than I would be comfortable with.

EDIT: In the comments, there was discussion of using the "redistribute to the right" (or left in this case) method. Reading up about this, redistributing to the right will result in the Kaplan Meier curves (see background here). Unfortunately, the Kaplan Meier curves can lead to heavily biased estimates of the mean! As a simulation, see this R code:

library(survival)

# Simulating right censored data
sim_surv_data = function(n = 100){
  t_time = rexp(n)
  c_time = rexp(n)
  is_observed = t_time < c_time
  y = c_time
  y[is_observed] = t_time[is_observed]
  ans = data.frame(y = y, obs = is_observed)
  return(ans) 
}

# Compute mean from KM curve
km_to_mean = function(km_fit){
  sum_obj = summary(km_fit)
  times = sum_obj$time
  # Convert survival function to PMF
  surv = sum_obj$surv
  shifted_survs = c(1, head(surv, length(surv) - 1))
  pmf = shifted_survs - surv 
  ans = sum(pmf * times)
  return(ans)
}

sim_once = function(n = 100){
  df = sim_surv_data(n)
  km_fit = survfit(Surv(df$y, event=df$obs, type='right')~1)
  ans = km_to_mean(km_fit)
  return(ans)
}

res = NULL 
for(i in 1:400){
  res[i] = sim_once()
}

estimated_mean = round(mean(res), 3)
se = round(sd(res) / sqrt(length(res) - 1), 3)
msg = "Estimated Mean From KM Curve\n True mean = 1"
msg = paste0(msg, "\nEstimate = ", estimated_mean, " SE = ", se)
hist(res, 
     main = msg, 
     xlab = "Estimated Mean")

This results in the following:

Answer 2

Here is a consistent estimator based on method of moments. Suppose

\begin{align} X &\sim N(\mu,\sigma^2)\\ Y &=\max(U,X)\\ p &=P[Y>U]\\ q &=E[Y-U]\\ r &=E[(Y-U)Y] \end{align} where we estimate $p,q,r$ from averages of the data in the problem.

Those variables are defined by integrals, and the integrals lead to an equation which we can solve numerically for $\mu$: $$\sqrt{\frac{r-q\mu}{2\pi p}}\,\exp\left(\frac{-p(\mu-U)^2}{2(r-q\mu)}\right)=q-p(\mu-U).$$

This equation has the following interpretation: Let $s=\sqrt{(r-q\mu)/p}$, which is equal to $\sigma$ if we have the exact values for $\mu,p,q,r$. Let $Z\sim N(\mu,s)$. Then $$f_Z(U) = \frac{q - p(\mu-U)}{s^2}.$$

Example: Suppose we use this procedure with 100 standard normal variables with a left cutoff at $u=1$. Simulating this 1000 times, I found estimated means that are only slightly biased, averaging around -0.006 with a standard deviation around 0.3. The accuracy is even better for lower cutoffs.

roots = c()
for(i in 1:100000){
  u = 1
  x = rnorm(100)
  y = pmax(x,u)
  p = mean(y>u)
  q = mean(y-u)
  r = mean((y-u)*y)
  
  f = function(mu){
    s2 = max(0,(r-q*mu)/p)
    ex = exp(-(u-mu)^2/(2*s2))
    return(sqrt(s2/(2*pi)) * ex - (q-p*(mu-u)))
  }
  
  if(f(-r/q)>0)
    root = r/q
  else
    root = uniroot(f, c(-r/q, r/q))$root
  roots = c(roots, root)
}
hist(roots)
c(mean(roots), sd(roots))

Answer 3

Here is maybe the simplest answer: If

\begin{align} p &= P[X>U] \\ m &= \text{median}\{X:X>U\}\\ q_u &= \Phi^{-1}(1-p)\\ q_m &= \Phi^{-1}(1-p/2)\\ \end{align} then $u=\mu+\sigma q_u$ and $m = \mu + \sigma q_m$ lead to the estimate of: $$\mu \simeq \frac{uq_m-mq_u}{q_m-q_u}.$$ This does as well as my moments-based answer for $u<0$, is more volatile and more biased for $0<u<2$, and is more volatile but less biased for $u>2$.

Since this estimator is simple, maybe I or maybe someone else can figure out how to correct its bias.