5
$\begingroup$

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.

$\endgroup$
1
  • 4
    $\begingroup$ Consistent is easy; unbiased not so much. One ad hoc method I have used, with some trepidation -- it was necessary to be consistent with the OLS estimator in cases where no censoring actually occurs -- is to use the MLE for the mean but adjust the MLE estimate of the error variance by a factor of $(n-p-1)/n$ for the purpose of computing confidence and prediction limits. There's also an interesting issue here: would it suffice just to be unbiased, or should it be unbiased conditional on the number of censored values? $\endgroup$
    – whuber
    Commented Jul 21, 2023 at 19:19

3 Answers 3

2
$\begingroup$

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: histogram of means

$\endgroup$
19
  • $\begingroup$ If not all observations are subject to censoring and you have a decent number of unccensored observations you can use multiple imputation to derive a more distribution-free estimate of the mean. $\endgroup$ Commented Jul 22, 2023 at 12:06
  • 1
    $\begingroup$ @FrankHarrell I think that approach will be flawed. If you use a parametric distribution to impute the censored values, you're circling back to the potential bias of using a parametric distribution. If you impute from the KM curves, you will have a big issue of properly sampling from the tails (since it's a step function) which will likely lead to much more bias than just using an incorrect parametric distribution. $\endgroup$
    – Cliff AB
    Commented Jul 22, 2023 at 15:22
  • $\begingroup$ I should have clarified. Use predictive mean matching or redistribution to the left/right, not parametric imputation. This imputation doesn't use KM. $\endgroup$ Commented Jul 22, 2023 at 16:11
  • 2
    $\begingroup$ @FrankHarrell: after reading up, I learned this reproduces the KM curves...but those can have a very biased estimate of the mean. See my edits. $\endgroup$
    – Cliff AB
    Commented Jul 23, 2023 at 16:06
  • 1
    $\begingroup$ @MattF: in my simulation, both the censoring time and event time are exp(1) distributed. It's right censored, not left. It's gotten off topic from the question, but it does show that a very large bias (ie 20% with a moderate sample size + moderately strong censoring) can be induced from using the KM curves to estimate the mean. Note that in the case with a constant cutoff, this approach fails much worse, since nothing is greater than the cutoff. . $\endgroup$
    – Cliff AB
    Commented Jul 23, 2023 at 19:41
3
$\begingroup$

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.

enter image description here

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))
$\endgroup$
2
  • $\begingroup$ Nice analysis. But I believe this estimator will be biased, though. $\endgroup$
    – whuber
    Commented Jul 22, 2023 at 17:34
  • $\begingroup$ +1, this answer makes sense to me in the context of my problem. Before I do a direct comparison with the Bayesian method though, I'll have to try some simulations to see how much bias I get under different parameters. Thank you! I'll wait a bit longer before accepting to see if I can get some more explanation on the other answer. $\endgroup$
    – wzbillings
    Commented Jul 25, 2023 at 2:13
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.