3
$\begingroup$

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: testdata

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: testdata fit

$\endgroup$
  • $\begingroup$ Your model description should make $y_i$ conditional on $z_i$. $\endgroup$ – Xi'an Apr 23 '15 at 6:28
  • $\begingroup$ 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] ? $\endgroup$ – Xi'an Apr 23 '15 at 6:31
  • $\begingroup$ and your conditional posteriors on $\beta$ and $\tau$ are incorrect. $\endgroup$ – Xi'an Apr 23 '15 at 6:54
  • $\begingroup$ @Xi'an, thanks for the comments. Can you elaborate on what is incorrect about the $\beta$ and $\tau$ conditional posteriors? $\endgroup$ – theochli Apr 23 '15 at 12:26
  • $\begingroup$ 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. $\endgroup$ – Xi'an Apr 23 '15 at 13:41
0
$\begingroup$

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)
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.