# Bayesian Mixture Model Gibbs Sampler for two linear relationships

I am attempting to use a Gibbs Sampler to model a mixture of two groups, where the group membership is defined by a linear relationship conditional on x. Both groups have the same slope and intercept, but one group has an additional, positive constant added to the common intercept. My question is that my set up does not seem to be correctly separating out the two groups. (see plot with result below)

$$z_i = \begin{cases} 1, & \text{if observation i's response is in group 1 (with positive shift)}, \\ 0, & \text{if observation i's response is in group 2}. \end{cases}$$

$$y_i = \begin{cases} \text{N} (\gamma + \beta x_i + \tau z_i\, , \, \sigma^2 ) & \text{ with prob. } \lambda, \\ \text{N} (\gamma + \beta x_i \,,\, \sigma^2) & \text{ with prob. } 1- \lambda. \end{cases}$$

The (simulated) data look like this:

And here is my code:

#synthetic data a synthetic dataset of 100 points
truebeta <-0.003
truegamma <-0
truetau <-15
truesigma2 <-1
truez <- rbinom(100, 1, 0.50)
truex <- seq(from=1, to= 1000, by = 10)
truey <- rnorm(100, truegamma+truebeta*truex+truetau*truez, sqrt(truesigma2+truesigma2*truez))
plot(truex, truey)

xs <-truex
ys <-truey
n <-length(xs)

niters <-1000

logcondfxn <- function(y, z, gamma, beta, x, tau, sigma2){
ndist1 <-dnorm(y[z==1], gamma + beta*x[z==1] + tau, sqrt(sigma2))
ndist2 <-dnorm(y[z==0],gamma + beta*x[z==0], sqrt(sigma2))
sum(log(ndist1)[is.finite(log(ndist1))]) + sum(log(ndist2)[is.finite(log(ndist2))])
}

library(MCMCpack)

#initialize values
lambda <- rep(0.5, niters)
beta <- rep(0.001, niters)
tau <- rep(0, niters)
gamma <- rep(0, niters)
sigma2 <- rep(5, niters)

### Gibbs Sampler with M-H ####
for (i in 2: niters) {
#sample zs from bernoulli with prob
znum <-lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
zdenom1 <- lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
zdenom2 <- (1-lambda[i-1])*dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs, sd=sqrt(sigma2[i-1]))
prob <-znum/(zdenom1+zdenom2)
prob[is.na(prob)] <-0
z <- rbinom(n,1, prob)

#sample gamma
ystar <- ys - beta[i-1]*xs- tau[i-1]*z
#what is left from this is the intercept gamma regress ystar on x
gamma[i] <- rnorm(1, summary(lm(ystar~1))$coef[1], sqrt(sigma2[i-1]/n )) #sample beta ystar <- ys - gamma[i] - tau[i-1]*z beta[i] <- rnorm(1, summary(lm(ystar~xs))$coef[2], sqrt(sigma2[i-1]/n ))

#sample tau
ystar <- ys[z==1] - beta[i]*xs[z==1] - gamma[i]
tau[i] <- rnorm(1, min(mean(ystar), 0, na.rm = TRUE), sqrt(sigma2[i-1]/n))

#sample lambda
lambda[i] <-rbeta(1, sum(z)+1, n - sum(z)+1)

#sample sigma2 using metropolis-hastings
sigma2.prop <- rgamma(1, sigma2[i-1], 1)
Mnum <- logcondfxn(ys, z, gamma[i], beta[i], xs, tau[i], sigma2.prop)
Mnum[is.na(Mnum)] <-0
Mden <- logcondfxn(ys, z, gamma[i], beta[i], xs, tau[i], sigma2[i-1])
Mden[is.na(Mden)] <-0

#Hastings step to correct of assymmetry of proposed sigma2 dist
logHr <-log(dgamma(sigma2[i-1], sigma2.prop, 1)) - log(dgamma(sigma2.prop, sigma2[i-1], 1))

#Calc logr
logr <- Mnum - Mden + logHr

logu <- log(runif(1, 0, 1))
sigma2[i] <- ifelse(logr>=logu, sigma2.prop, sigma2[i-1])
}


I tried plotting the iterations to gauge convergence, and the chain is not converging.

plot(sigma2, type = "l")
plot(lambda, type = "l", xlim = c(0,100))
plot(beta, type = "l")
plot(gamma, type = "l", xlim= c(0,50))
plot(tau, type = "l", xlim = c(0,50))


And to diagnose autocorrelation:

acf(beta)
acf(lambda)
acf(gamma)
acf(tau)


I haven't thinned or thrown out burn-in yet, but from these preliminary plots, I am not getting what I hope to (regression lines through each group of points). Instead, group membership goes virtually all to one group (0) in the z vector, and both lines are between the two groups.

plot(xs,ys)
abline(mean(gamma), mean(beta))
abline(mean(gamma) + mean(tau), mean(beta))


Here is a plot of what I am getting right now:

• Your model description should make $y_i$ conditional on $z_i$. – Xi'an Apr 23 '15 at 6:28
• In logcondfxn, why do you write gamma - beta*x[z==1] - tau and gamma - beta*x[z==0] instead of gamma + beta*x[z==1] + tau and gamma + beta*x[z==0] ? – Xi'an Apr 23 '15 at 6:31
• and your conditional posteriors on $\beta$ and $\tau$ are incorrect. – Xi'an Apr 23 '15 at 6:54
• @Xi'an, thanks for the comments. Can you elaborate on what is incorrect about the $\beta$ and $\tau$ conditional posteriors? – theochli Apr 23 '15 at 12:26
• How did you compute those full conditionals? it seems you replicated the line for the $\gamma$'s but the posterior variance on $\beta$ depends on $x$ and the posterior on $\tau$ depends on the number of $z$'s equal to 1. I do not have time to re-derive those for myself but advise you to check them afresh. – Xi'an Apr 23 '15 at 13:41

Here is the corrected code to run this problem:

lambda<- rep(0.5, niters)
beta<- rep(0.001, niters)
tau<- rep(0, niters)
gamma<- rep(0, niters)
sigma2<- rep(5, niters)

for (i in 2: niters){
## 1. Sample z. sample zs from bernoulli with prob
znum<-lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
zdenom1<- lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
zdenom2<- (1-lambda[i-1])*dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs, sd=sqrt(sigma2[i-1]))

prob<-znum/(zdenom1+zdenom2)

#prob[is.na(prob)]<-0

z<- rbinom(n,1, prob)

## 2. Sample lamda from Beta dist
lambda[i]<-rbeta(1, sum(z)+1, n - sum(z)+1)

## 3. Sample mus (beta and gamma are common, tau is only for group1)
#### 3a. GAMMA
ystar<- ys - beta[i-1]*xs- tau[i-1]*z
#what is left from this is the intercept gamma regress ystar on x
gamma[i]<- rnorm(1, summary(lm(ystar~1))$coef[1], sqrt(sigma2[i-1]/n )) #### 3b. BETA ystar<- ys - gamma[i] - tau[i-1]*z beta[i]<- rnorm(1, summary(lm(ystar~xs-1))$coef[1],vcov(lm(ystar~xs-1)))
#### 3c. TAU
n1<- max(sum(z), 0.00001)
ystar<- ys[z==1] - beta[i]*xs[z==1] - gamma[i]
tau[i]<- rnorm(1, max(mean(ystar), 0, na.rm = TRUE), sqrt(sigma2[i-1]/n1))

## 4. Sample sigma2 from invgamma
sigma2[i]<- 1/rgamma(1, n/2 +1 , sum((ys - gamma[i]- beta[i]*xs - tau[i]*z)^2, na.rm = TRUE)/ 2)
}