A Bayesian model using random-walk Metropolis method for data augmentation

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)

}

#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?