2
$\begingroup$

I'm trying to estimate using MLE the parameters of a (Gamma) probability distribution given some data, but the data contain a good number of samples below the minimum detection limit (MDL) of the measurement method. These show up as zeroes in the dataset. I came up with a piecewise likelihood function where all 0s are assigned the mean probability density over [0,MDL], i.e. CDF(MDL)/MDL. And then all values above the MDL just have their likelihood calculated as normal.

Questions:

  1. Is there anything WRONG with this approach? Are there better approaches?
  2. Are there any non-obvious situations where it will fail to reliably estimate the parameters?
  3. Has this been done in the literature before? If so, does somebody have a citation I can use?

Reprex below...

x_actual <- rgamma(100, 1,1)
MDL <- .25
x_with_0s <- x_actual
x_with_0s[x_actual<MDL] <- 0

likfunc <- function(k, theta){
  mean_lik_below_MDL <- pgamma(MDL, k, theta)/MDL
  liks <- log(dgamma(x, k, theta))
  liks[x==0] <- log(mean_lik_below_MDL)
  loglik <- -sum(liks)
  return(loglik)
}

This seemingly returns pretty good parameter estimates:

> x <- x_actual

> mle(likfunc, start=list(k=.5, theta=.5), method = "L-BFGS-B", lower = c(0,0))

Call:
mle(minuslogl = likfunc, start = list(k = 0.5, theta = 0.5), 
    method = "L-BFGS-B", lower = c(0, 0))

Coefficients:
       k    theta 
1.082315 1.077583 

> x <- x_with_0s

> mle(likfunc, start=list(k=.5, theta=.5), method = "L-BFGS-B", lower = c(0,0))

Call:
mle(minuslogl = likfunc, start = list(k = 0.5, theta = 0.5), 
    method = "L-BFGS-B", lower = c(0, 0))

Coefficients:
       k    theta 
1.031814 1.029893 
>

One thing I don't like is that there's a discontinuity in the PDF for the likelihood function, but I don't know if this is really an issue.

x <- seq(0,5, .01)
pdf <- function(k, theta){
  mean_lik_below_MDL <- pgamma(MDL, k, theta)/MDL
  liks <- dgamma(x, k, theta)
  liks[x<MDL] <- mean_lik_below_MDL
  return(liks)
}
plot(pdf(1,1)~x)

enter image description here

$\endgroup$
1

2 Answers 2

4
$\begingroup$

Your data is left-censored because you only know these measurements have a value below MDL. So, using Survival Analysis to incorporate left-censored data into the likelihood function you would have to replace all zeros for the MDL, and because you only know that $X < MDL$, their contribution to the likelihood is $F(MDL)$.

The code could be something like

set.seed(1996)

x_actual <- rgamma(100, 1,1)
MDL <- .25
x_censored <- x_actual
x_censored[x_actual<MDL] <- MDL

dt <- data.frame(x_obs = x_censored)
dt$censored <- ifelse(dt$x_obs == MDL, 1, 0)

likfunc <- function(k, theta){
  liks <- log(dgamma(x$x_obs, k, theta))
  liks[x$censored == 1] <- log(pgamma(x[x$censored == 1, 1], k, theta))
  loglik <- -sum(liks)
  return(loglik)
}

This yields the estimates

> x <- data.frame(x_obs = x_actual, censored = rep(0, length(x_actual)))
> 
> mle(likfunc, start=list(k=.5, theta=.5), method = "L-BFGS-B", lower = c(0,0))

Call:
mle(minuslogl = likfunc, start = list(k = 0.5, theta = 0.5), 
    method = "L-BFGS-B", lower = c(0, 0))

Coefficients:
       k    theta 
1.069465 1.275040 
> 
> x <- dt
> 
> mle(likfunc, start=list(k=.5, theta=.5), method = "L-BFGS-B", lower = c(0,0))

Call:
mle(minuslogl = likfunc, start = list(k = 0.5, theta = 0.5), 
    method = "L-BFGS-B", lower = c(0, 0))

Coefficients:
       k    theta 
1.045510 1.247177 

Klein's Survival Analysis explains this contribution to the likelihood function.

$\endgroup$
2
  • 1
    $\begingroup$ Survival analysis typically concerns right-censored data, although in various ways it can be pressed into service for analyzing left-censored data. But it does not specifically deal with truncated data--nor is the present question about truncated data. It's also not about replacing the censored data with zeros. A closer consideration of the question is advisable. $\endgroup$
    – whuber
    Mar 4, 2020 at 23:21
  • 2
    $\begingroup$ You are right, I apologize. I've edited my answer to show how to incorporate left-censored data in the likelihood function. I think that this approach with Survival Analysis has the theoretical foundations you need. $\endgroup$ Mar 5, 2020 at 0:11
2
$\begingroup$

What you want to use for the likelihood is the following. Suppose $f(x)$ is the density function and $F(x)$ is the cumulative distribution function at $x$. If from a sample of size $n$ with a threshold of $x_0$, there are $n_0$ below that threshold and $n-n_0$ at or above $x_0$, the likelihood is

$$L=F(x_0)^{n_0} \times \prod_{i=1}^{n-n_0} f(x_i)$$

where $x_1, x_2,\ldots,x_{n-n_0}$ are the values at or above the threshold. As stated by @Ben-ReinstateMonica in a comment, the log of the likelihood is what is typically maximized because of better numeric stability:

$$\ln{L}=n_0 \ln{F(x_0)}+\sum_{i=1}^{n-n_0}\ln f(x_i)$$

Note that the R function mle in the stats4 package minimizes $-\ln{L}$ rather than maximizing $\ln{L}$ which can cause a bit of confusion. (I'm guessing the reason for that is that mle uses the optim function which by default minimizes an objective function.)

That can be implemented for your example as

x_actual <- rgamma(100, 1, 1)
MDL <- 0.25
nCensored <- sum(x_actual < MDL)
xNotCensored <- x_actual[x_actual >= MDL]

library(stats4)

likfunc <- function(k, theta){
  -nCensored*log(pgamma(MDL, k, theta)) - sum(log(dgamma(xNotCensored, k, theta)))
}

(results <- mle(likfunc, start=list(k=.5, theta=.5), method = "L-BFGS-B", lower = c(0,0)))

Call:
mle(minuslogl = likfunc, start = list(k = 0.5, theta = 0.5), 
    method = "L-BFGS-B", lower = c(0, 0))

Coefficients:
        k     theta 
1.0068535 0.9319662 

Saving the results in an object gets you the other essential piece for an estimate: a measure of precision (here in the form of a standard error).

summary(results)
Maximum likelihood estimation

Call:
mle(minuslogl = likfunc, start = list(k = 0.5, theta = 0.5), 
    method = "L-BFGS-B", lower = c(0, 0))

Coefficients:
       Estimate Std. Error
k     0.8912489  0.1402007
theta 0.9177082  0.1721334

-2 log L: 269.0074 

You can also get the estimated covariance matrix and correlation between the maximum likelihood estimators of k and theta with the following:

# Covariance matrix
  results@vcov
               k      theta
k     0.01965624 0.01991023
theta 0.01991023 0.02962990

# Correlation between estimators
results@vcov[1,2]/sqrt(results@vcov[1,1]*results@vcov[2,2])
[1] 0.8250135

In general, suppose you had "complete" censoring in the sense that all of the data is binned into $b$ bins with boundaries $x_0,x_1,\ldots,x_b$ and associated freqencies $n_1,n_2,\ldots,n_b$. The likelihood function is

$$L=\prod_{i=1}^{b} (F(x_i)-F(x_{i-1}))^{n_i}$$

with the log of the likelihood being

$$\ln L=\sum_{i=1}^{b} n_i \ln(F(x_i)-F(x_{i-1}))$$

The maximum likelihood estimates would be the values of the parameters that maximize $\ln L$.

$\endgroup$
4
  • 1
    $\begingroup$ Good answer (+1). May I suggest you add an additional equation to show the corresponding form of the log-likelihood $\ell = \ln L$, since this function is often the one that is maximised in the MLE procedure (for numerical stability reasons). $\endgroup$
    – Ben
    Mar 5, 2020 at 0:35
  • 1
    $\begingroup$ @Ben-ReinstateMonica. Thanks. Good idea. I've added that and a blurb about including a measure of precision. $\endgroup$
    – JimB
    Mar 6, 2020 at 0:36
  • $\begingroup$ OP here: So it would appear to me, then, that the only difference between my approach and the two answers here is that, for the truncated data, I assign a probability of F(x0)/x0 (a probability density, specifically the mean probability density of x < x0) , whereas you give the probability as the probability mass F(x0). In your example, how can you multiply both probability masses (for x < x0) and probability densities (for x > x0) into a coherent likelihood for for all x? $\endgroup$
    – user278411
    Mar 6, 2020 at 0:54
  • 1
    $\begingroup$ This doesn't invalidate your argument but $F(x_0)/x_0$ is not the mean probability density given that $x<x_0$. For a gamma distribution with parameters $k=1$ and $\theta=1$ and $x_0=1/4$, the mean probability density given that $x<1/4$ is $\frac{1}{2} \left(1+\frac{1}{\sqrt[4]{e}}\right)\approx 0.8894$ and $F(1/4)/(1/4)=\frac{1}{2} \left(1+\frac{1}{\sqrt[4]{e}}\right)\approx 0.8848$. For those parameters, the two values are similar but not necessarily for other sets of parameters. How would your approach fare if the gamma distribution was right censored? Would the mean density be zero? $\endgroup$
    – JimB
    Mar 6, 2020 at 1:38

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.