Maximum Likelihood Estimation
As discussed in the references, MPN is estimated by maximizing likelihood. Combining the notaton of Blodgett (2) and Jarvis et al. (7), we write the likelihood function as:
\[L = L(\lambda; x_i, n_i, z_i, i = 1,...,k) = \prod_{i=1}^k \binom{n_i}{x_i} {(1-exp(-{\lambda}{z_i}))} ^ {x_i} {(exp(-{\lambda}{z_i}))} ^ {n_i-x_i}\]
where
- \(\lambda\) is the microbial density (concentration) to be estimated
- \(k\) is the number of dilution levels
- \(x_i\) is the number of positive tubes at the \(i^{th}\) dilution level
- \(n_i\) is the total number of tubes at the \(i^{th}\) dilution level
- \(z_i\) is the amount of inoculum per tube at the \(i^{th}\) dilution level
As an R function:
#likelihood
L <- function(lambda, positive, tubes, amount) {
binom_coef <- choose(tubes, positive)
exp_term <- exp(-lambda * amount)
prod(binom_coef * ((1 - exp_term) ^ positive) * exp_term ^ (tubes - positive))
}
L_vec <- Vectorize(L, "lambda")As is typical of maximum likelihood approaches, mpn()
uses the score function (derivative of the log-likelihood) to solve for
\(\hat{\lambda}\), the maximum
likelihood estimate (MLE) of \(\lambda\) (i.e., the point estimate of
MPN). However, let’s demonstrate what is happening in terms of
the likelihood function itself. Assume we have 10g of undiluted inoculum
in each of 3 tubes. Now we use a 10-fold dilution twice (i.e., the
relative dilution levels are 1, .1, .01). Also assume that exactly 1 of
the 3 tubes is positive at each dilution level:
#MPN calculation
library(MPN)
my_positive <- c(1, 1, 1) #xi
my_tubes <- c(3, 3, 3) #ni
my_amount <- 10 * c(1, .1, .01) #zi
(my_mpn <- mpn(my_positive, my_tubes, my_amount))
#> $MPN
#> [1] 0.1118255
#>
#> $MPN_adj
#> [1] 0.08725827
#>
#> $variance
#> [1] 0.00431064
#>
#> $var_log
#> [1] 0.3447151
#>
#> $conf_level
#> [1] 0.95
#>
#> $CI_method
#> [1] "Jarvis"
#>
#> $LB
#> [1] 0.03538176
#>
#> $UB
#> [1] 0.3534288
#>
#> $RI
#> [1] 0.005814327If we plot the likelihood function, we see that \(\hat{\lambda}\) maximizes the likelihood:
my_mpn$MPN
#> [1] 0.1118255
lambda <- seq(0, 0.5, by = .001)
my_L <- L_vec(lambda, my_positive, my_tubes, my_amount)
plot(lambda, my_L, type = "l", ylab = "Likelihood", main = "Maximum Likelihood")
abline(v = my_mpn$MPN, lty = 2, col = "red")
If none of the tubes are positive, the MLE is zero:
no_positive <- c(0, 0, 0) #xi
(mpn_no_pos <- mpn(no_positive, my_tubes, my_amount)$MPN)
#> [1] 0
L_no_pos <- L_vec(lambda, no_positive, my_tubes, my_amount)
plot(lambda, L_no_pos, type = "l", xlim = c(-0.02, 0.2), ylab = "Likelihood",
main = "No Positives")
abline(v = mpn_no_pos, lty = 2, col = "red")
If all of the tubes are positive, then no finite MLE exists:
all_positive <- c(3, 3, 3) #xi
mpn(my_tubes, all_positive, my_amount)$MPN
#> [1] Inf
lambda <- seq(0, 200, by = .1)
L_all_pos <- L_vec(lambda, all_positive, my_tubes, my_amount)
plot(lambda, L_all_pos, type = "l", xlim = c(0, 100), ylim = c(0, 1.1),
ylab = "Likelihood", main = "All Positives")
abline(h = 1, lty = 2)
From a practical perspective, if all the tubes are positive, then the scientist should probably further dilute the sample until some tubes are negative.