0
$\begingroup$

I'm trying to fit the mean $\mu$ of right-censored gaussian data ($n$ samples) in a toy example (let's assume $\sigma^2=1$ is known), and the censoring happens always at the same value $s$. As far as I can understand the log likelihood is $$L(\hat{\mu})=\sum_{i=1}^n y_i~\log \mathcal{N}(x_i;\hat{\mu})+(1-y_i)~\log(1-F(s;\hat{\mu})),$$ where $y_i=1$ for observed and $y_i=0$ for censored, $\mathcal{N}$ is the gaussian density function, and $F$ is the gaussian cumulative density function.

Translating this into R code in which the true $\mu=0$, censoring at $s=0.5$, and plotting the log likelihood over a grid of possible $\mu$s:

n <- 3e3
mu <- 0
x1 <- rnorm(n, mean=mu)

x2 <- x1
s <- 0.5

# Censoring status, 1=event (not censored), 0=censored
y <- as.integer(x2 > s)
table(y)
x2[y == 0] <- s


# Grid of means
mv <- seq(-5, 5, length=50)

# log likelihoods for each mean, original data
ll1 <- sapply(mv, function(g) sum(dnorm(x1, mean=g, log=TRUE)))

# log likelihoods for each mean, censored data
ll2 <- sapply(mv, function(g) sum(dnorm(x2, mean=g, log=TRUE)))

# log likelihood: sum of log likelihoods at non-censored observations plus
# sum of log(1 - F(s, mu)) for censored observations
ll3 <- sapply(mv, function(g) {
   sum(y * dnorm(x2, mean=g, log=TRUE) +
      (1 - y) * pnorm(x2, mean=g, lower=FALSE, log=TRUE))
})

png("censored.png", width=900, height=900, res=200)
matplot(mv, cbind(ll1, ll2, ll3), type="l", lwd=1.5,
   xlab=expression(hat(mu)), ylab="log likelihood")
abline(v=mv[which.max(ll1)], col=1, lty=1, lwd=1.5)
abline(v=mv[which.max(ll2)], col=2, lty=2, lwd=1.5)
abline(v=mv[which.max(ll3)], col=3, lty=3, lwd=1.5)
dev.off()

enter image description here

If this is correct, I'm not clear about why the maximum likelihood of the censored data but ignoring the censoring structure (ll2 in red) produces a better estimate of $\mu$ than ll3 in green which accounts for censoring (the true $\mu$ is 0). Unless I've made a mistake somewhere...

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.