2
$\begingroup$

I am trying to understand the effect of censored data on parameter estimation, in particular with respect to Bayesian modeling. Here is an illustrating example with some code:

Assume you have $N = n+m$ i.i.d observations from an exponential distribution with parameter $\lambda=3$. We will also define a value $L$ such that any observation less than $L$ is censored. With this in mind, the following is the likelihood for the data \begin{align} \ell(x|\lambda) &= \prod_{i=1}^mF(x)\prod_{i=1}^nf(x)\\ &=\prod_{i=1}^m(1-e^{-\lambda L})\prod_{i=1}^n\lambda e^{-\lambda x_i} \end{align} where $m$ and $n$ are the number of censored and noncensored observations.

We also place a prior such that $\lambda\sim\text{Gamma}(0.01,0.01)$. Now, no matter what the level of censoring is, as long as I have a single uncensored data point my models seems to do well recovering the true $\lambda$. However, as soon as I make the entire data set censored, I cannot estimate $\lambda$ well. This doesn't seem too shocking to me given that you don't have any information other than the censoring point, however, I would have thought that the estimate of $\lambda$ would be whatever the prior mean is, which is not the case. Should you be recovering the prior if you have all censored data?

Here is come code that estimates $\lambda$ in an MCMC. To change the percent of censoring you just need to change the value for prob in the quantile() function.

likelihood <- function(x, censor, lambda){
  
  llik1 <- log(1 - exp(-lambda * x[censor == 1])) 
  llik2 <- log(lambda) - lambda * x[censor == 0]
  lik  <- sum(llik1) + sum(llik2)
  
  return(lik)
}

prior <- function(lambda, sigma){
  
  prior <- dgamma(lambda, .01, .01, log = TRUE)
  
  return(prior)
}

posterior <- function(x, censor, lambda){
  
  post <- likelihood(x, censor, lambda) + prior(lambda)
  
  return(post)
}


set.seed(4)
N = 100
x = rexp(N, 3)
L = quantile(x, prob = 1) # Censoring point
censor = ifelse(x <= L, 1, 0) # Censoring indicator

x[censor == 1] <- L

S <- 10000
lambda <- rep(NA, S)
lambda[1] <- 1

for(i in 2:S){
  
  # MCMC for lambda
  lambda.star = lambda[i-1] + rnorm(1,0)
  
  if(lambda.star < 0){
    alpha = 0
  }else{
    ratio <- exp(posterior(x, censor, lambda.star) - posterior(x, censor, lambda[i-1]))
    alpha <- min(1, ratio)
  }
  
  if(runif(1) < alpha){
    lambda[i] <- lambda.star
  }else{
    lambda[i] <- lambda[i - 1]
  }
  
}


lambda <- lambda[5000:S]
mean(lambda)
$\endgroup$

2 Answers 2

2
$\begingroup$

It is not the case that you should recover the prior if all your data is censored; there is still information in the data.

Consider the following situation. The data is i.i.d Gaussian with unknown mean $\mu$ (whether or not we know the variance is irrelevant to this example.) My prior on the mean is that $\mu \sim N(0,10)$. My censoring point is $\gamma = 10$. I draw 100 observations; every one is above the censoring point. Without getting into the details of what my posterior should look like, why, given that I have observed 100 values $> 10$ and none $< 10$, and the underlying distribution is Normal, would I hold to the belief that there's a 50-50 chance that the true mean $< 0$ (which is a straightforward implication of the prior)? For any proposed value of $\mu \leq 10$, there is only a $2^{-100}$ probability that all 100 observations would lie above $\mu$ (thanks to the symmetry of the Gaussian); clearly our posterior should be substantially different than our prior if we want to capture this information.

$\endgroup$
6
  • $\begingroup$ Thanks for that very useful answer. I think that makes sense to me. An additional question then would be is the parameter estimate under 100% censoring reliable? In the sense that will the MCMC converge to a stationary distribution if left to run long enough? $\endgroup$
    – John Smith
    Oct 15, 2021 at 21:58
  • $\begingroup$ If you have proper priors on all your parameters, writing informally, almost certainly, otherwise, possibly not. For the latter, consider the standard improper prior on $\mu \sim U(-\infty, \infty)$; if all we know is that we have 100 observations $>10$, the posterior will still be improper because there's no information to cause the right tail to start dropping off at a sufficiently great rate to be integrable. Also note that the Gamma(0.01,0.01) distribution is no longer considered a good choice of uninformative prior for scale parameters. $\endgroup$
    – jbowman
    Oct 15, 2021 at 22:11
  • $\begingroup$ I strongly suspect that you can construct counterexamples even with proper priors, but I also strongly suspect you have to work to do it, so your chances of doing so by accident are pretty low. $\endgroup$
    – jbowman
    Oct 15, 2021 at 22:13
  • $\begingroup$ Thanks for the useful information. Can you recommend a useful noninformative prior for scale parameters? $\endgroup$
    – John Smith
    Oct 15, 2021 at 22:13
  • $\begingroup$ stat.columbia.edu/~gelman/research/published/taumain.pdf , although not directly applicable to your problem, contains suggestions that I have found helpful in the past (early days when I still used the Gamma dist'n.) I'd use the half-t if I were in your shoes. $\endgroup$
    – jbowman
    Oct 15, 2021 at 22:16
0
$\begingroup$

If all the data points are censored then this is just like the case where you have an observable set of binary data (i.e., all data points zero or one) and you observe all the values are zeroes (or equivalently, all the values are ones). This is the famous "sunrise problem" examined by Laplace. Under the sunrise problem, you can still get an inference for the probability of the outcome that has occurred in the data.

Transferring over to your model, this is the probability of a value being censored, which is a function of $\lambda$. Consequently, you still get an inference for $\lambda$ in this case, and the posterior is not equivalent to the prior --- i.e., the data is still informative about the parameter. Of course, the inferences is not great (compared to the case where you have non-censored data), but that is hardly surprising, since your data essentially reduces to one outcome from a binary.

$\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.