<p>Concerning your last question about estimating CI ranges, there are three common methods for ML estimators:</p>
<ol>
<li>Variance estimation from the inverted Hessian matrix.</li>
<li>Jackknife estimator for the variance (simpler and more stable, if the Hessian is estimated numerically, but computationally more expensive)</li>
<li>Bootstrap CIs (the computatianally most expensive approach).</li>
</ol>
<p>For bootstrap CIs, you do not need to implement them yourself (bias correction, e.g. can be tricky), but can rely on the R library <em>boot</em>.</p>
<p>Incidentally, I have written a summary with R code for all three approaches two years ago: <a href="https://arxiv.org/abs/1807.03582" rel="nofollow noreferrer">Construction of Confidence Intervals</a> (see section 5). For the method utilizing the Hessian Matrix, e.g., the outline is as follows:</p>
<pre><code>lnL &lt;- function(theta1, theta2, ...) {
  # definition of the negative (!)
  # log-likelihood function...
}

# starting values for the optimization
theta0 &lt;- c(start1, start2, ...)

# optimization
p &lt;- optim(theta0, lnL, hessian=TRUE)
if (p$convergence == 0) {
  theta &lt;- p$par
  covmat &lt;- solve(p$hessian)
  sigma &lt;- sqrt(diag(covmat))
}
</code></pre>
<p>The function <em>mle</em> from <em>stats4</em> already wraps the covrainace matrix estimation and retruns it in <em>vcov</em>. In the practical use cases in which I have tried this (paired comparison models), though, this estimation was rather unstable, and I have resorted to the jackknife method instead.</p>
R code for maximum likelihood estimate from a specific likelihood function
<p>I have been trying to generate R code for maximum likelihood estimation from a log likelihood function in a <a href="https://onlinelibrary.wiley.com/doi/full/10.1111/risa.12880" rel="nofollow noreferrer">paper</a> (equation 9 in page 609). Authors in the paper estimated it using MATLAB, which I am not familiar with. So I tried to generate codes in R.</p>
<p>Here is the snapshot of the log likelihood function in the paper:</p>
<p><a href="https://i.stack.imgur.com/rrZz1.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/rrZz1.png" alt="enter image description here" /></a></p>
<p>, where</p>
<p><em>r</em>: Binary decision (0 or 1) indicating infested plant(s) detection (1) or not (0).</p>
<p><em>e</em>: Inspection efficiency. This is known.</p>
<p><em>n</em>: Sample size</p>
<p>The overall objective is to estimate plant infestation rate (gamma: γ) and epsilon (<em>e</em>) based on binary decision of presence and absence of infested plants instead of using infested plant(s) detected. So, the function has only binary information (<em>r</em>) of infested plant detection and sample size. Since epsilon (<em>e</em>) is known or fixed, the actual goal is to estimate gamma (γ) in a population.</p>
<p>Another objective is to compare estimated infestation rates from above with ones in hypergeometric sampling formula in <a href="https://pubmed.ncbi.nlm.nih.gov/19139495/" rel="nofollow noreferrer">another paper</a> (in page 6). The formula is:</p>
<p><a href="https://i.stack.imgur.com/ZRk30.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ZRk30.png" alt="enter image description here" /></a></p>
<p>This formula generates required sample size to detect infested plants with selected probability (e.g., 95) given an infested rate. For example:</p>
<pre><code># Sample size calculation function
fosgate.sample1 &lt;- function(box, p, ci){ # Note: box represent total plant number
  ninf &lt;- p*box
  sample.size &lt;- round(((1-(1-ci)^(1/ninf))*(box-(ninf-1)/2)))
  #sample.size &lt;- ceiling(((1-(1-ci)^(1/ninf))*(box-(ninf-1)/2)))
  sample.size
}

fosgate.sample1(box=100, p = .05, ci = .95) # where box: population or total plants, p: infestation rate, and ci: probability of detection
## 44
</code></pre>
<p>The idea is if sample size (e.g., 44) and binary decision data are provided the log-likelihood function can be used to estimate infestation rate and the rate may be close to anticipated rate (e.g., .05). <strong>Ultimately, I would like to compare plant infestation rates (gamma: γ) estimated from the log likelihood function above and D/N in the sample size calculation formula (second) or p in the sample size code below.</strong></p>
<p>I generated R code for the log-likelihood described above.</p>
<pre><code>### MLE with stat4
library(stats4)
# Log-likelihood function
plant.inf.lik &lt;- function(inf.rate){
logl &lt;- suppressWarnings(
        sum((1-insp.result)*n*log(1-inf.rate) + 
            insp.result*log(1-(1-inf.rate)^n))
        )
return(-logl)
}
</code></pre>
<p>Using the sample size function (i.e., fosgate.sample1) I generated sample sizes for various cases of total plant (or box) and anticipated detection rate (p) in the function. Since I am also interested in error/confidence ranges of estimated plant infestation rates, I used bootstrapping to calculate range of estimates (I am not sure if this is appropriate/acceptable). Here is the final code I generated:</p>
<pre><code>### MLE and CI with bootstrapping with multiple scenarios
plant &lt;- c(100, 500, 1000, 5000, 10000, 100000) # Total plant number
ir &lt;- seq(.01, .2, by = .01) # Plant infestation rate
df.result &lt;- data.frame(expand.grid(plant=plant, inf.rate = ir))
df.result$sample.size &lt;- fosgate.sample1(box=df.result$plant, p=df.result$inf.rate, ci=.95) # Sample size
df.result$insp.result &lt;- 1000 # Shipment number (can be replaced with random integers)
df.result &lt;- df.result[order(df.result$plant, df.result$inf.rate, df.result$sample.size), ]
rownames(df.result) &lt;- 1:nrow(df.result)
df.result$est.mean &lt;- 0
#df.result$est.median &lt;- 0
df.result$est.lower.ci &lt;- 0
df.result$est.upper.ci &lt;- 0
df.result$nsim &lt;- 0
str(df.result)
head(df.result)

# Looping
    est &lt;- rep(NA, 1000)
for(j in 1:nrow(df.result)){
    for(i in 1:1000){
        insp.result &lt;- sample(c(rep(1, df.result$insp.result[j]-df.result$insp.result[j]*df.result$inf.rate[j]), 
                    rep(0, df.result$insp.result[j]*df.result$inf.rate[j])))
        ir &lt;- df.result$inf.rate[j]
        n &lt;- df.result$sample.size[j]
        insp.result &lt;- sample(insp.result, replace = TRUE)
        est[i] &lt;- mle(plant.inf.lik, start = list(inf.rate = ir*.9), method = &quot;BFGS&quot;, nobs = length(insp.result))@coef
    df.result$est.mean[j] &lt;- mean(est, na.rm = TRUE)
#   df.result$est.median[j] &lt;- median(est, na.rm = TRUE)
    df.result$est.lower.ci[j] &lt;- quantile(est, prob = .025, na.rm = TRUE)
    df.result$est.upper.ci[j] &lt;- quantile(est, prob = .975, na.rm = TRUE)
    df.result$nsim[j] &lt;- length(est)
    }
}

# Significance test result
sig &lt;- ifelse(df.result$inf.rate &gt;= df.result$est.lower.ci &amp; df.result$inf.rate &lt;= df.result$est.upper.ci, &quot;no sig&quot;, &quot;sig&quot;)
table(sig)

# Plot
library(ggplot2)
library(reshape2)
df.result$num &lt;- ave(df.result$inf.rate, df.result$plant, FUN=seq_along)
df.result.m &lt;- melt(df.result, id.vars=c(&quot;plant&quot;, &quot;sample.size&quot;, &quot;insp.result&quot;, &quot;est.lower.ci&quot;, &quot;est.upper.ci&quot;, &quot;nsim&quot;, &quot;num&quot;))
df.result.m$est.lower.ci &lt;- ifelse(df.result.m$variable == &quot;inf.rate&quot;, NA, df.result.m$est.lower.ci)
df.result.m$est.upper.ci &lt;- ifelse(df.result.m$variable == &quot;inf.rate&quot;, NA, df.result.m$est.upper.ci)
str(df.result.m)

ggplot(data = df.result.m, aes(x = num, y = value, group=variable, color=variable, shape=variable))+
    geom_point()+
    geom_errorbar(aes(ymin = est.lower.ci, ymax = est.upper.ci), width=.5)+
    scale_y_continuous(breaks = seq(0, .2, .02))+
    xlab(&quot;Index&quot;)+
    ylab(&quot;Plant infestation rate&quot;)+
    facet_wrap(~plant, ncol = 3)
</code></pre>
<p>When I ran the code, I was able to obtain results and to compare estimated (est.mean) and anticipated (inf.rate) infestation rates as shown in the plot below.</p>
<p><a href="https://i.stack.imgur.com/8yBvU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8yBvU.png" alt="enter image description here" /></a></p>
<p>If results are correct, plot indicates that estimation looks fine but off for greater infestation rates.</p>
<p>Also, I always got warning messages without &quot;suppressWarnings&quot; function and occasionally error messages below. I have no clue how to fix them.</p>
<pre><code>## Warning messages
## 29: In log(1 - (1 - inf.rate)^n) : NaNs produced
## 30: In log(1 - inf.rate) : NaNs produced

## Error message (occasionally)
## Error in solve.default(oout$hessian) : 
## Lapack routine dgesv: system is exactly singular: U[1,1] = 0
</code></pre>
<p>My questions are:</p>
<ul>
<li>Is R function (plant.inf.lik) for maximum likelihood estimation of the log-likelihood function appropriate?</li>
<li>Should I take care of warning and error messages? If yes, how? Again, I have no clue how to fix...</li>
<li>Is bootstrapping (resampling?) method appropriate to estimate CI ranges and/or standard error?</li>
</ul>
<p>I found <a href="https://math.stackexchange.com/questions/40319/maximum-likelihood-estimate-of-hypergeometric-distribution-parameter">this link</a> useful for alternative approach. Although I am still working both approaches together, results seem different (maybe following question).</p>
<p>Any suggestion would be greatly appreciated.</p>