How to know dispersion if $\mu$ is close to or below 0 (chance-corrected beta-binomial model)

Background

In sensory science, "replication" means having a panelist in a taste panel do multiple rounds of the same test. You cannot just count those additional rounds as additional panelists, because panelists might show consistent performance. In that case, replications would just confirm the abilities of individual panelists rather than provide more information about the difference between the two stimulants. So, theory suggests using a corrected beta-binomial model in which the abilities of the panelists are assumed to be beta-distributed, not between 0 and 1, but between $C$ (the guessing probability of the taste test used) and 1:

$$P(x|n,a,b,C)=\frac{(1-C)^n}{B(a,b)}\binom{n}{x}\sum_{i=0}^{x}\binom{x}{i}\left(\frac{C}{1-C}\right)^{x-i}B(a+i,n+b-x)$$

$C$ is the guessing probability of the test method used (in my case $C=\frac{1}{3}$),
$n$ is the number of replications (in my case $n=2$ for all panelists),
$x$ is the number of successes.

Log likelihood:

$$L=\sum_{i=1}^{k}log[Pr(x_i|n_i,\mu,\gamma)]$$

$k$ is the number of panelists (in my case $k=54$),
$n_i$ and $x_i$ are the number of replications and successes of the $i$th panelist,
$\mu$ is the mean,
$\gamma$ is the over-dispersion parameter (sometimes also referred to as $\rho$).

$a$ and $b$ can be calculated from $\mu$ and $\gamma$:

$$a=\mu\left(\frac{1}{\gamma}-1\right)\text{ and }b=(1-\mu)\left(\frac{1}{\gamma}-1\right)$$

(Bi, Sensory discrimination tests and measurements, 2006, as shown in Gacula et al., Statistical Methods in Food and Consumer Research, 2009 § 10.3. See also Brockhoff, The statistical power of replications in difference tests, 2003 and Næs et al., Statistics for Sensory and Consumer Science, 2010 § 7.7.)

My test results

I had 54 panelists perform 2 replications of a sensory test with a guessing probability of 33% (n=2, k=54, C=1/3), all comparing the same two pancake recipes. Depending on the taste difference between the pancake recipes, you would expect the observed proportion of successes $p_c$ to run up from 33% in case of no taste difference to 100% in case of obvious taste difference.

However, my test results show panelists correctly grouped the pancakes in only 29 of the 108 rounds (27%). Here is the R-code for my data:

round1 <- c(0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1,
1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0)
round2 <- c(1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0)
x <- round1 + round2
n <- c(rep(2,length(x)))
dat <- data.frame(id = factor(1:length(x)), x, n)


I have no reason to assume the panelists gave wrong answers on purpose, nor that I mixed up samples, so I guess it was just chance that led to this 'low' performance.

Maximum Likelihood Estimation

As described in Gacula (see link above), I created an R-function of the above-mentioned log-likelihood of the corrected beta-binomial (CBB) model, called it llcbb, then ran nlminb to find MLE-estimates of $\mu$ and $\gamma$. The estimate for $\mu$ is close to 0 and the estimate for $\gamma$ will usually stay near its start value:

fit<-nlminb(start=c(.1, .9), function(x) -1*llcbb(dat,mu=x[1],gamma=x[2],1/3)
,control=list(trace=1))
fit$par # [1] 5.986091e-13 8.974503e-01 fit<-nlminb(start=c(.1, .1), function(x) -1*llcbb(dat,mu=x[1],gamma=x[2],1/3) ,control=list(trace=1)) fit$par
# [1] 3.040535e-12 9.976926e-02


The following picture shows my findings in black and some possible configurations of the CBB-model. The cyan line shows that when $\mu$ approaches 0, it does not really matter anymore what $\gamma$ is. I guess this explains why nlminb leaves $\gamma$ untouched.

Moment estimation

Using other formulas provided by Gacula (see link above) I calculated the following moment estimates: $$\hat{\mu}=-0.097222\text{, }\hat{\gamma}= 0.349006$$

Wikipedia states the following about negative moment estimates:

Note that these estimates can be non-sensically negative which is evidence that the data is either undispersed or underdispersed relative to the binomial distribution. In this case, the binomial distribution and the hypergeometric distribution are alternative candidates respectively.

(Beta-binomial distribution. (2014, February 7). In Wikipedia, The Free Encyclopedia. Retrieved 14:42, February 14, 2014, from http://en.wikipedia.org/w/index.php?title=Beta-binomial_distribution&oldid=594290522)

I understand the suggestion to revert to the binomial distribution, but even then I still feel I need to be sure about the dispersion, because the upper limit of the confidence interval on $p_c$ varies quite a bit depending on whether I can use the replication data (in case of no dispersion) or not (in case of over-dispersion):

binom.test(29,108,p=1/3,alternative="less")
# 95 percent confidence interval:
#  0.0000000 0.3475966

binom.test(15,54,p=1/3,alternative="less")
# 95 percent confidence interval:
#  0.0000000 0.3950352


Question

How to deal with possible dispersion ($\gamma$ or $\rho$) in beta-binomial distribution in case $\mu$ is estimated to approach 0, or to be below 0?

Rune Haubo Bojesen Christensen, author of the sensR-package, provided the following solution:

# Solution by Rune Haubo Bojesen Christensen
# dat contains the data (see question)
library(sensR)
X <- dat[, -1]

prof.mu <- function(mu) {
## Profile likelihood of mu: optimize the likelihood over gamma for
## given value of mu. Return value is a nlminb fit.
fun <- function(gamma)
sensR:::setParBB(rho=rho, par=c(mu, gamma))
fit <- nlminb(.5, objective=fun, lower=0+1e-3, upper=1-1e-3)
fit
}

## Sequence of values at which to evaluate the profile likelihood of
## mu:
(museq <- c(1e-6, seq(.01, .99, by=.01)))

## Fit profile likelihood
fit.list <- lapply(museq, prof.mu)
## Check convergence for all mu:
sapply(fit.list, "[[", "convergence")
all(sapply(fit.list, "[[", "convergence") == 0) ## Ok
sapply(fit.list, "[[", "par") ## gamma estimates - seems ok.
## Negative log-likelihood:
(nll <- sapply(fit.list, "[[", "objective"))


## Profile negative log-likelihood
## I think zero corresponds to the binomial model.
par(mfrow=c(2,2))
plot(museq, -2 * nll, type = "l")
abline(h=-3.84)
plot(museq, -2 * nll, type = "l", ylim=c(-10, 0), xlim=c(0, .2))
abline(h=-3.84)

## So the upper 95% CI limit is around pd = 0.07 corresponding to

The MLE for $\mu$ is as I found earlier ($\mu\to0$, unable to determine $\gamma$), but by calculating the log-likelihood values of a range of other possible $\mu$'s (between 0 to 1) combined with their respective MLE's of $\gamma$, it is possible to determine a likelihood ratio confidence interval by comparing the likelihood ratio-statistic with the value of a chi-square distribution with 1 degree of freedom. The value 3.84 that is compared against, is the 95% centile of the $\chi^2$-distribution on one degree of freedom. We’re looking at the profile likelihood of $\mu$, where $\gamma$ is regarded as a nuisance parameter, and since we’re only ‘testing’ one parameter we only use df=1.
Then, regarding the MLE of $\gamma$ for $\mu\to0$ I was originally looking for, I now understand that this does not matter, as the outcome of the log-likelihood is what matters for creating the profile, and that (in this case) for the log-likelihood for $\mu\to0$, it does not matter what $\gamma$ is. Actually, even at the found upper 95% CI limit of around $\mu$ = $p_d$ = 0.07, $\gamma$ still hit the lower boundary of 0+1e-3 set in Christensen's prof.mu function.
If I still want to say something about the MLE of $\gamma$ for $\mu\to0$, looking at the MLE's for $\gamma$ where it can be determined (for $\mu\geq0.1$), we can say that the MLE’s of $\mu$ and $\gamma$ are both essentially zero:
plot(museq[11:100], sapply(fit.list, "[[", "par")[11:100], type = "l")