MCMC for Probit/Logit model with some 1's flipped to 0's

Question

I would like help constructing a sampler for the following model, which is the latent variable interpretation of either logistic or probit glm (doesn't matter which one to me), with a small twist: there is a probability $p$ that "successes" are flipped to failures after the fact.

Model:

Let $X$ represent the design matrix. Let $\vec{y}$ represent a vector of binary responses

$z = X \vec{\beta} + \vec{\epsilon}$

$\epsilon_i \sim Logis(0,s)$ or $\epsilon_i \sim Normal(0,\phi)$

$y_i = 1$ if $z_i > 0$ AND $rbernoulli(p) == 0$, $y_i = 0$ otherwise

$p \sim beta(\alpha,\beta)$

$\beta_j \sim Normal(0,0.001)$

Whether or not an observation is flipped from 1 to 0 is independent of an observation's "ability" (X values), conditional on our other parameters. Observations may not be flipped from 0 to 1.

Toy Example:

At the University of Foo, statistics students are selected to receive the Bar Award based on things such as academic achievement, extra-curricular participation, and community service. Professors convene and decide to give a certain number of students the award, writing the awardees' names on separate slips of paper, and handing them to the department head to give out the awards. However, as they often are, this department head is quite clumsy and, at random, loses some of the deserving students' names.

We need to build a model to predict student success, taking into account that some of the "losers" are deserving of the award. Last year, we found that 3 of the 20 original awardee's names were lost. Using a Jeffreys' prior with binomial likelihood to estimate the probability of droppage, we develop a $Beta(3.5,17.5)$ posterior on $p$.

I would love to see references to papers where something like this is done, or some help in constructing a sampler for this problem. My own efforts have not lead anywhere.

My Math:

As requested, I am putting here my attempt at a solution.

Above, I have the prior for $\beta$ as a normal. It doesn't really matter to me what it is, so long as it is not very informative. My attempts will use marginal Jeffreys' priors for $\beta$ and $\phi$:

$P(\beta) \propto 1$

$P(\phi) \propto 1/\phi$

These are the Jeffreys priors for regular old linear models, I assume I can use them in the latent varaible model as well. Please inform me if I am wrong.

$$(\vec{z} - X \vec{beta}) \sim N(0,\phi I)$$

$$P(\vec{z},\vec{\beta},\phi | X,y) \propto | \phi I| ^{1/2} e^{-1/2 (z - X \beta)' (\phi I)(z - X\beta)}$$

$$\propto | \phi I| ^{1/2} e^{-\phi/2 (z' - \beta' X')(z - X\beta)}$$

$$\propto | \phi I| ^{1/2} e^{-\phi/2 (z'z - \beta'X'z - z'X\beta + \beta'X'X\beta)}$$

$$\propto | \phi I| ^{1/2} e^{-\phi/2(\beta'X'X\beta - 2z'X\beta )}$$

It would seem $\beta$ is normally distributed conditional on $z$ and $\phi$.

I am not sure how to tie all of this to $y$ and $p$.

Answer 1

There are two steps to this problem: 1) Obtaining an unbiased estimate of $\beta$ from the fully observed model and 2) Estimating the probability model for the error generating mechanism.

In general the approach to handling these types of problems is using the EM algorithm. A few assumptions are necessary to gain any traction. Assumptions might be along the lines of 1) what is the actual probability model for so called random flipping? Is it an unobserved additive model? Is it a mixture model? Does it depend on other observed factors? Is it even an estimable model? 2) what is the expected number of true positives in the $y$ variable? Is this known? Should it have been a fixed amount, or consistent with some target? 3) What is the known conditional distribution of the $X$ variables given $y$?

So the problem as stated is not defined sufficiently, but below an approach is outlined.

If we assume $X$ is normally distributed conditional upon $Y$, and that the label flipping is done completely at random, a clustering based approach is outlined as follows. Using the distribution of $X$ where $y$ is 0, one would look for probabilistic outliers estimated from a normal mixture model, and presume a small cluster of high risk X values here should belong to the $Y=1$ group. After flipping labels, one would re-estimate $\beta$ for predictions of risk, then re-run the algorithm iteratively until convergence.

An example in practice is simulated here:

## example
set.seed(1234)
n <- 1000
p <- 0.3

e.prob <- 0.05 ## error rate, completely at random model
y.orig <- rbinom(1000, 1, p)
x <- numeric(1000)
x[y.orig == 1] <- rnorm(sum(y.orig), 10, 2)
x[y.orig == 0] <- rnorm(sum(1-y.orig), 3, 3)

## true beta
f <- glm(y ~ x, family=binomial)
# int: -9.2, slope: 1.2

err <- rbinom(1000, 1, e.prob)
y.obs <- y.orig
y.obs[err==1] <- 0

hist(x[y.obs==1], col=rgb(1,0,0,.5), breaks = -10:20, ylim=c(0, 120))
hist(x[y.obs==0], col=rgb(0,1,0,.5), add=T, breaks=-10:20, ylim=c(0, 120))

## initialize em:
y.working <- y.obs

mu <- tapply(x, y.working, mean)
sig <- tapply(x, y.working, var)

## loglikelihood
loglik.prev <- sum(dnorm(x[y.working==0], mean=mu[1], sd = sig[1], log=T))+sum(dnorm(x[y.working==1], mean=mu[2], sd = sig[2], log=T))
itr <- 0
repeat({
  itr <- itr+1
  ## classify least likely negative X as a false negative value
  minLL.y0 <- which.min(dnorm(x[y.working==0], mean=mu[1], sd=sig[1], log=T))
  index<-which(y.working==0)[minLL.y0]
  y.working[index] <- 1

  ## evaluate updated likelihood, stop if it gets worse
  mu <- tapply(x, y.working, mean)
  sig <- tapply(x, y.working, sd)

  loglik.curr <- sum(dnorm(x[y.working==0], mean=mu[1], sd = sig[1], log=T))+sum(dnorm(x[y.working==1], mean=mu[2], sd = sig[2], log=T))

  if(loglik.curr < loglik.prev) {
    y.working[index] <- 0 ## set the old value back since that was the ML value
    break()
  } else {
    loglik.prev <- loglik.curr
  }
})

table(y.working, y.obs)
table(y.working, y.orig)

par(mfrow=c(2,1))
hist(x[y.obs==1], col=rgb(1,0,0,.5), breaks = -10:20, ylim=c(0, 120), main='Observed')
hist(x[y.obs==0], col=rgb(0,1,0,.5), add=T, breaks=-10:20, ylim=c(0, 120))
hist(x[y.working==1], col=rgb(1,0,0,.5), breaks = -10:20, ylim=c(0, 120), main='Corrected')
hist(x[y.working==0], col=rgb(0,1,0,.5), add=T, breaks=-10:20, ylim=c(0, 120))
legend('topright', bty='n', pch=22, pt.bg = c(rgb(1,0,0,.5), rgb(0,1,0,.5)), c('Positives', 'Negatives'))

3 / 12, not very good really