If missing data process is known and it is MNAR, is it possible to get an unbiased estimate of parameter?

Question

I was reading on missing data handling techniques from "Statistical Analysis with Missing Data Second Edition" by Little and Rubin. In example 6.22 they consider that $y_1, y_2, ... y_n$ are i.i.d. exponential random variables. i.e. $f(y_i) = \theta^{-1}e^{-\dfrac y \theta}$. Further $y_i$ is only recorded if it is less than a cutoff $c$. Suppose $r$ out of $n$ observations are recorded and rest are missing. A corresponding missingness indicator $M_i$ is 0 if it is missing and 1 if observed. They show that the estimate based on the joint likelihood of missing data process and observed data gives

$\hat{\theta} = \dfrac {\Sigma_{i=1}^{r} y_i + (n-r)c} r $

My question is if this estimate is unbiased?

Here is what I've tried so far:

Since the observed $y_i$ belong to a truncated exponential distribution with support $[0, c)$, $E(y_i) = \dfrac \theta {1-e^{-\dfrac c \theta}} (1-e^{-\dfrac c \theta}(1+\dfrac c \theta))$.

$\Rightarrow E(\hat{\theta}) = \dfrac \theta {1-e^{-\dfrac c \theta}} (1-e^{-\dfrac c \theta}(1+\dfrac c \theta)) + nc * E(\dfrac {1} {r}) - c$

Since $r$, the number of observed data which are less than $c$ follows a binomial distribution with probability of success $p(y < c)$, $\dfrac {1} {r}$ is the inverse of binomial distribution. I tried to calculate this expected value but I ended up concluding that it is infinity. So it seems to me that the bias also infinity?

Link: Inverse Binomial Expected value

Simulations that I performed to check unbiasedness of estimator

I also performed some simualations to get a rough idea about how biased the estimator is, and from the simulations it seems to me it may be unbiased and its variance in inversely proportional to the cutoff $c$. However I can't conclude from simulations. Nevertheless here is my R code:

library(ggplot2)

numObs = 100
len = 100
cutoff = 30
origEstimator = vector("numeric", len)
newEstimator = vector("numeric", len)
for(i in 1:len){
  sample = rexp(n = numObs, rate = 0.05)
  origEstimator[i] = mean(sample)

  truncSample = sample[sample<cutoff]
  r = length(truncSample)
  newEstimator[i] = (sum(truncSample) + (numObs-r)*cutoff)/r
}

ggplot(data = data.frame(mean=c(origEstimator, newEstimator), type=c(rep("Complete data", len), rep("Adjusted", len)))) + 
  geom_density(aes(x=mean, color=type))

Screenshot from the book:

Answer 1

This is an interesting question.

First I will show that $r>0$. If $r=0$, then there is no observed data, and the likelihood function is no longer concave, thus this statistical problem is not well defined. Given $r>0$, $E(1/r)$ should be finite.

Let $p=\mathrm{Pr}(y_i<c) = 1-\mathrm{exp}(-c/\theta)$. We have $$E(1/r) = \sum_{r=1}^n \binom{n}{r} p^r (1-p)^{n-r}/r = n p (1-p)^{n-1} F(1, 1, 1-n;2, 2;p/(p-1)),$$ where $F(\cdot)$ is the generalized hypergeometric function. $E(\hat{\theta})$ can be computed by plugging $E(1/r)$ into the equation in the question.

The following numerical computation shows that $E(\hat{\theta}) = \theta$, i.e., $\hat\theta$ that considers the MNAR (missing not at random) mechanism is unbiased. Note that genhypergeo() used to calculate the generalized hypergeometric function has numerical error, but the above summation can be computed exactly.

library(hypergeo)

theta_hat = function(theta, n = 100) {
  c = .5*theta # arbitrary c can be used
  p = 1 - exp(-c/theta)
  theta*(1-exp(-c/theta)*(c/theta+1))/(1-exp(-c/theta)) + (n *c) *
    (genhypergeo(U=c(1,1,1-n), L=c(2,2), z=p/(p-1)) * n * p * (1-p)^(n-1)) - c
}

theta = 1:10
plot(theta, sapply(theta, theta_hat), ylab='theta_hat')
abline(a=0, b=1)