MLE for Normally Distributed Data Censored by Limit of Detection

Description

Values below the limit of detection (LOD) are a problem in several fields of science, and there are numerous approaches for replacing the missing data. We present a new mathematical solution for maximum likelihood estimation that allows us to estimate the true values of the mean and standard deviation for normal distributions and is significantly faster than previous implementations. The article with the details was submitted to JSS and can be currently seen on <https://www2.arnes.si/~tverbo/LOD/Verbovsek_Sega_2_Manuscript.pdf>.

Details

The DESCRIPTION file:

Index of help topics:

mlelod                  Estimates the parameters of the normal
                        distribution.
mlelod-package          MLE for Normally Distributed Data Censored by
                        Limit of Detection

Values below the limit of detection (LOD) are a problem in several fields of science, and there are numerous approaches for replacing the missing data. Thic package uses a new mathematical solution for maximum likelihood estimation that allows us to estimate the true values of the mean and standard deviation for normal distributions and is significantly faster than previous implementations. The core function is the function mlelod with three parameters: the size of the sample, the values of the sample and the value of limit of detection. The function returns two estimates: for mu (expected value) and sigma (standard deviation).

Author(s)

Gregor Sega [aut, cre]

Maintainer: Gregor Sega <gregor.sega@fmf.uni-lj.si>

References

The article with the derived method is submited to Journal of Statistical Software

Estimates the parameters of the normal distribution.

Description

The function returns the two estimates for the parameters of the normal distribution.

Usage

mlelod(n, censoreddata, lod)

Arguments

Value

Author(s)

Gregor Sega

References

The article with the derived method is submited to Journal of Statistical Software

Examples

## The function is currently defined as
mlelod <- function (n, censoreddata, lod) 
{
    k <- length(censoreddata)
    s <- sum(censoreddata)
    s2 <- sum(censoreddata^2)
    g <- function(x, n, k, s, s2, lod) {
        a <- (k * (lod^2 - x^2) + s2 - 2 * lod * s)/(x * (k * 
            lod - s))
        b <- s - k * (k * x^2 + lod * s - s2)/(k * lod - s) - 
            x * (n - k) * dnorm(a, 0, 1)/pnorm(a, 0, 1)
        return(b)
    }
    h <- function(x) {
        return(g(x, n, k, s, s2, lod))
    }
    lsi <- sqrt(s2/k - (s/k)^2)
    sigmaEst <- uniroot(h, lower = lsi/4, upper = 4 * lsi, tol = 1e-06)$root
    muEst <- (k * sigmaEst^2 + lod * s - s2)/(k * lod - s)
    return(list(muEst = muEst, sigmaEst = sigmaEst))
  }
##define the parameters of the normal distribution
mu <- 5
sigma <- 4
##define the size of the sample and the value of lod
n <- 100
lod <- 2
##generate normally distributed data and extract the observable values (the ones exceeding lod)
data <- rnorm(n, mu, sigma)
data2 <-  Filter(function(x) x>lod,data)
##run the function to obtain the estimates
mlelod(n, data2, lod)