Simple case of MNAR missing data

Question

I believe I have a very simple problem with missing data, but I'm a bit lost because all the materials I read seem to be focused on much more complicated cases.

I have a random variable $X$ which has a Binomial distribution with parameters $n$ and $p$ (i.e., $X\sim B(n,p)$), where $p$ is unknown ($n$ is known). I have $K$ independent samples of $X$, let them be $X_1,\dotsc,X_K$ but some of them are missing, in the sense that I observe the values $R_iX_i$, where $R_i\in\{0,1\}$. I know that $R_i = 0$ when $X_i< C$, for some known constant $C$, otherwise $R_i=1$. This means that my missing data are MNAR (Missing Not At Random), but I know the "mechanism" leading to their being missing, which I believe is what makes my case easier.

Edit: it's a case of left-censoring with Binomial data.

How can I estimate $p$? I'm also interested in finding a confidence interval for $p$.

References and links are much appreciated.

Thanks in advance.

Answer 1

You can still use maximum likelihood fairly easily without having to resort to the E-M algorithm, which is commonly used for missing data problems. Just put a wrapper around the binomial distribution so that it calculates $p(x=\text{NA}),\space p(x=C),\space p(x=C+1), \dots$ instead of $p(x=0),\space p(x=1),\space p(x=2), \dots$ and maximize the likelihood using a univariate maximization routine.

In R, I choose to minimize -2*LL (log likelihood) instead. The aforementioned "wrapper" is the calculation of term1 inside the m2.log.like function below.

# Random generation of data
n <- 25
p <- 0.7
C <- 16
K <- 100
x <- rbinom(K,n,p)
x[x<C] <- NA 

sum(is.na(x))
[1] 19 # 19 of the 100 samples were censored

m2.log.lik <- function(p, n, C, x)
{
  term1 <- sum(dbinom(x[!is.na(x)],n,p,log=TRUE))
  term1 <- term1 + sum(is.na(x)) * log(pbinom(C-1,n,p))
  -2*term1
}

opt.ll <- optimize(m2.log.lik, lower=0.01, upper=0.99, n=n, C=C, x=x)
opt.ll
$minimum
[1] 0.7004906  # Our MLE of p.

$objective
[1] 408.2332

To calculate a confidence interval, you can use the asymptotic chi-square distribution of -2*LL and do some root finding to find values for $p$ such that the associated values of -2*LL equal the value associated with the MLE plus some appropriate number based on the chi-square distribution:

m2.log.lik.rf <- function(p, n, C, x, tgt)
{
  log.lik(p, n, C, x) - tgt
}

# Calculate a 95% confidence interval
tgt <- m2.log.lik(phat, n, C, x) + qchisq(0.975, 1)  
ub95 <- uniroot(m2.log.lik.rf, lower=phat, upper=phat+0.1, n=n, C=C, x=x, tgt=tgt)
lb95 <- uniroot(m2.log.lik.rf, lower=phat-0.1, upper=phat, n=n, C=C, x=x, tgt=tgt)

lb95$root
[1] 0.6792
ub95$root
[1] 0.721157

The sample interval was pretty tight, but then, we did have 100 observations from a binomial(25,0.7), so even with the censoring we should do pretty well.