# NOT RUN {
burnin <- 100
nmc <- 500
thin <- 1
p <- 100 # number of predictors
ntrain <- 250 # training size
ntest <- 100 # test size
n <- ntest + ntrain # sample size
q <- 10 # number of true predictos
beta.t <- c(sample(x = c(1, -1), size = q, replace = TRUE), rep(0, p - q))
x <- mvtnorm::rmvnorm(n, mean = rep(0, p), sigma = diag(p))
zmean <- x %*% beta.t
z <- rbinom(n, size = 1, prob = boot::inv.logit(zmean))
X <- scale(as.matrix(x)) # standarization
# Training set
ztrain <- z[1:ntrain]
Xtrain <- X[1:ntrain, ]
# Test set
ztest <- z[(ntrain + 1):n]
Xtest <- X[(ntrain + 1):n, ]
posterior.fit <- logiths(z = ztrain, X = Xtrain, method.tau = "halfCauchy",
burn = burnin, nmc = nmc, thin = 1,
Xtest = Xtest)
posterior.fit$BetaHat
# Posterior processing to recover the true predictors
cluster <- kmeans(abs(posterior.fit$BetaHat), centers = 2)$cluster
cluster1 <- which(cluster == 1)
cluster2 <- which(cluster == 2)
min.cluster <- ifelse(length(cluster1) < length(cluster2), 1, 2)
which(cluster == min.cluster) # this matches with the true variables
# }
Run the code above in your browser using DataLab