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

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$.

• Have you tried finding the posterior distribution for this? – Greenparker Jul 12 '16 at 21:11
• @Greenparker are you referring to the conditional posteriors for each parameter? I'm afraid I'm fairly new to bayesianism, and did not get very far at all trying to work the math out, to the point where I felt it was better omitted from the post. – John Madden Jul 12 '16 at 21:21
• I am referring to the posterior, not the conditional posterior. This will help you understand (1) Is the posterior intractable. Because if the posterior turns out to be a known distribution, you don't need MCMC(!). It is likely, the posterior is not a known distribution. (2) It will force you to write down the parameters of the posterior, and understand the quantities you are dealing with. (3) Most often we use Metropolis-Hastings to sample from the posterior, which needs the kernel os the posterior. So you will have to write it down. That would be the first step. – Greenparker Jul 12 '16 at 21:26
• Alright, I'll give her another go and edit the OP with my progress. – John Madden Jul 12 '16 at 21:30
• @Greenparker OK I gave it a go – John Madden Jul 12 '16 at 22:39

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

• Very cool work. Playing around with your code, it seems that this "adjustment" gets the parameter estimates closer to the truth. Thanks a bunch! This serves me as a very cool intro to EM type algorithms. – John Madden Jul 13 '16 at 2:12