0
$\begingroup$

I have following model,

z=beta1+X1*beta2+e; e~N(0,sigma2)

prob=exp(z)/(1+exp(z)); and y= 1 with prob, 0 with (1-prob).

I have the following prior:

beta~N(betahat,A^(-1)); betahat=c(0,0) and A=0.01*I

sigma2~inv-gamma(nu,ssq); nu=0.01*N and ssq=1

I want to use a data augmentation method to solve this model. The steps go like,

  1. Based on the likelihood function, draw the augmented value of z for each observation.

  2. With the draws of z, obtain beta by standard univariate regression of z on X=[1, X1].

So, I write the following R code,

#simulate dataset
N=200
X=cbind(runif(N),runif(N))

beta_true=c(1,2)
k=length(beta_true)

z_true=X%*%beta_true+rnorm(N,sd=2)
prob=exp(z_true)/(1+exp(z_true)) 

y=rbinom(N,1,prob)

#prior 
betabar=c(1,2)
A=0.01*diag(k)
nu=0.01*N
ssq=1

#univeriate regression of y on X
uni_regression_one_gibbs<-function(y,X,nvar,last_beta,last_sigmasq,betabar,A,nu,ssq){

  n=length(y)
  tXX=crossprod(X)
  tXy=crossprod(X,y)
  IR=backsolve(chol(tXX/sqrt(last_sigmasq)+A),diag(nvar))
  beta_tilde=crossprod(t(IR))%*%(tXy/sqrt(last_sigmasq) + A%*%betabar)
  beta=beta_tilde+IR%*%rnorm(nvar)

  #draw sigmasq | beta
  res=y-X%*%beta
  s=crossprod(res)

  sigmasq=(s+nu*ssq)/rchisq(1,n+nu)
  sigmasq=as.vector(sigmasq)

  return(list(betadraw=beta,sigmasqdraw=sigmasq))
}

#metropolis functions----
#likelihood
logit_z_likelihood<-function(y,z){
  prob = exp(z)/(1+exp(z))
  prob = prob*y + (1-prob)*(1-y)
  sum(log(prob))
}

#random walk draw
RwDAMetroplis<-function(y,X,beta,zr,rej){
  z_old = zr
  #print(length(z_old))
  z_new = z_old + rnorm(1,sd=sqrt(0.015))
  #print(length(z_old))
  # data        
  lognew = logit_z_likelihood(y,z_new)
  logold = logit_z_likelihood(y,z_old)

  # prior
  logknew = -.5*(1/sigma2)*t(z_new-X%*%beta) %*% (z_new-X%*%beta)
  logkold = -.5*(1/sigma2)*t(z_old-X%*%beta) %*% (z_old-X%*%beta)
  # MH step
  alpha = exp(lognew + logknew - logold - logkold)
  if(alpha=="NaN") alpha=-1
  u = runif(n=1,min=0, max=1)
  if(u < alpha) { 
    zr=z_new
  } else {
    rej = rej+1  }

  return(list(z=zr,rej=rej))
}

##MCMC
iter=10000
rej=0

result.beta=matrix(0,nrow=iter,ncol=k)
result.sigma2=rep(0,iter)

#starting value
beta=rep(0,k)
sigma2=4
z=rep(0,N)

for(i in 1:iter){

  for(r in 1:N){
    #augmenting z for each observation
    oneMetropDraw=RwDAMetroplis(y[r],X[r,],beta,z[r],rej)
    z[r]=oneMetropDraw$z
    rej=oneMetropDraw$rej

  }

  # regresssion of z on x
  onedraw<-uni_regression_one_gibbs(z,X,k,beta,sigma2,betabar,A,nu,ssq)
  # store results
  beta<-onedraw$betadraw
  sigma2<-onedraw$sigmasqdraw
  result.beta[i,]<-beta
  result.sigma2[i]<-sigma2
}

However, it can not recover the true value of beta=[1,2]. The problem may be with the likelihood function, I'm not sure.

This model is similar to logit, but not the same: we have an error term for z here, instead of z=X*beta in logit.

I need latent variable z for later study so I have to use the data augmentation method. Is Metropolis here able to achieve this goal?

Any comments and suggestions welcomed!

$\endgroup$
  • $\begingroup$ What do you mean by "it cannot recover the true value of beta"? $\endgroup$ – Xi'an Aug 20 '17 at 15:01
  • $\begingroup$ @Xi'an Hi, I mean the estimated results of beta (store in result.beta, accessed by command apply(result.beta,2,mean)) is very very different from the true value of beta, which should be [1,2]. $\endgroup$ – Ding Li Aug 21 '17 at 1:40

Your Answer

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

Browse other questions tagged or ask your own question.