#check the tutorial
library(EMJMCMC)
library(clusterGeneration)
library(bindata)
library(ggplot2)
#***********************IMPORTANT******************************************************
# if a multithreaded backend openBLAS for matrix multiplications
# is installed on your machine, please force it to use 1 thread explicitly
# unless ncores in LogrRegr is reasonably small, in the latter case
# you might want to experiment with the combinations of blas_set_num_threads and ncores
library(RhpcBLASctl)
blas_set_num_threads(1)
omp_set_num_threads(1)
#***********************IMPORTANT******************************************************
# set the seed
set.seed(040590)
# construct a correlation matrix for M = 50 variables
M = 50
m = clusterGeneration::rcorrmatrix(M,alphad=2.5)
# simulate 1000 binary variables with this correlation matrix
X = bindata::rmvbin(1000, margprob = rep(0.5,M), bincorr = m)
# prepare the correlation matrix in the melted format
melted_cormat = reshape2::melt(cor(X))
# plot the heat-map of the correlations
ggplot2::ggplot(data = melted_cormat,
ggplot2::aes(x=Var1, y=Var2, fill=value)) +
ggplot2::geom_tile() +
ggplot2::theme(axis.title.x = ggplot2::element_blank(),
axis.title.y = ggplot2::element_blank(),
axis.text.x = ggplot2::element_blank())
# simulate Gaussian responses from a model with two-way interactions
# which is considered in S.4 of the paper
df = data.frame(X)
df$Y=rnorm(n = 1000,mean = 1+1.43*(df$X5*df$X9)+
0.89*(df$X8*df$X11)+0.7*(df$X1*df$X4),sd = 1)
help("LogicRegr")
# specify the initial formula
formula1 = as.formula(paste(colnames(df)[M+1],"~ 1 + ",
paste0(colnames(df)[-c(M+1)],collapse = "+")))
# Bayesian logic regression with the robust-g-prior
res4G = LogicRegr(formula = formula1, data = df,
family = "Gaussian", prior = "G", report.level = 0.5,
d = 15,cmax = 2,kmax = 15, p.and = 0.9, p.not = 0.1, p.surv = 0.2,
ncores = 32)
# Bayesian logic regression with the Jeffreys prior
res4J = LogicRegr(formula = formula1, data = df,
family = "Gaussian", prior = "J", report.level = 0.5,
d = 15, cmax = 2,kmax = 15, p.and = 0.9, p.not = 0.1, p.surv = 0.2,
ncores = 32)
# print the results for the robust g-prior
print(base::rbind(c("expressions","probabilities"),res4G$feat.stat))
#print the results for the Jeffreys prior
print(base::rbind(c("expressions","probabilities"),res4J$feat.stat))
# simulate Gaussian responses from a model with two-way interactions
# and an age effect which is an extension of S.4 of the paper
Xp = data.frame(X)
Xp$age = rpois(1000,lambda = 34)
Xp$Y=rnorm(n = 1000,mean = 1+0.7*(Xp$X1*Xp$X4) +
0.89*(Xp$X8*Xp$X11)+1.43*(Xp$X5*Xp$X9) + 2*Xp$age, sd = 1)
teid = sample.int(size =100,n = 1000,replace = F)
test = Xp[teid,]
train = Xp[-teid,]
# specify the initial formula
formula1 = as.formula(paste("Y ~ 1 +",
paste0(colnames(test)[-c(51,52)],collapse = "+")))
# specify the link function
g = function(x) x
# specify the parameters of the custom estimator function
estimator.args = list(data = train, n = dim(train)[1],
m =stri_count_fixed(as.character(formula1)[3],"+"),k.max = 15)
# specify the parameters of gmjmcmc algorithm
gmjmcmc.params = list(allow_offsprings=1,mutation_rate = 250,
last.mutation=10000, max.tree.size = 5, Nvars.max =15,
p.allow.replace=0.9,p.allow.tree=0.01,p.nor=0.01,p.and = 0.9)
# specify some advenced parameters of mjmcmc
mjmcmc.params = list(max.N.glob=10, min.N.glob=5, max.N=3, min.N=1,
printable = F)
# run the inference of BLR with a non-binary covariate and predicions
res.alt = pinferunemjmcmc(n.cores = 30, report.level = 0.2,
num.mod.best = 100,simplify = T,predict = T,test.data = test,
link.function = g,
runemjmcmc.params = list(formula = formula1,latnames = c("I(age)"),
data = train,estimator = estimate.logic.lm.jef,
estimator.args =estimator.args,
recalc_margin = 249, save.beta = T,interact = T,outgraphs=F,
interact.param = gmjmcmc.params,
n.models = 10000,unique = F,max.cpu = 4,max.cpu.glob = 4,
create.table = F,create.hash = T,pseudo.paral = T,burn.in = 100,
print.freq = 1000,
advanced.param = mjmcmc.params))
print(base::rbind(c("expressions","probabilities"),res.alt$feat.stat))
print(sqrt(mean((res.alt$predictions-test$Y)^2)))
print(mean(abs((res.alt$predictions-test$Y))))
library(HDeconometrics)
ridge = ic.glmnet(x = train[,-51],y=train$Y,family = "gaussian",
alpha = 0)
predict.ridge = predict(ridge$glmnet,newx = as.matrix(test[,-51]),
type = "response")[,which(ridge$glmnet$lambda == ridge$lambda)]
print(sqrt(mean((predict.ridge-test$Y)^2)))
print(mean(abs((predict.ridge-test$Y))))
tmean = 1+2*test$age+0.7*(test$X1*test$X4) +
0.89*(test$X8*test$X11)+1.43*(test$X5*test$X9)
print(sqrt(mean((tmean -test$Y)^2)))
print(mean(abs((tmean -test$Y))))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.