## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  error=TRUE,
  warning=FALSE
)

## ----echo=T,results='hide'----------------------------------------------------
library(bmabart)

## ----echo=T,results='hide'----------------------------------------------------
data(weight_behavior)
try0= bma.bart(pred=weight_behavior[,3], m=weight_behavior[,c(2,4:14)], 
               y=weight_behavior[,15], refy = 0, predref = "F",nskip=10,ndpost=100)
summary(try0)

## ----echo=T-------------------------------------------------------------------
summary(try0)

## -----------------------------------------------------------------------------
 summary(try0,method=2,RE=F)

## ----echo=T,results='hide'----------------------------------------------------
try1= bma.bart(pred=weight_behavior[,3], m=weight_behavior[,c(5:12)], 
               cova=weight_behavior[,2], y=weight_behavior[,1], 
               refy = 0, predref = "F",nskip=10,ndpost=20)
summary(try1)

## ----echo=T, results='hide'---------------------------------------------------
try2= bma.bart(pred=weight_behavior[,3], m=weight_behavior[,c(5:12)], 
               cova=weight_behavior[,2], mcov=weight_behavior[,13:14], 
               mclist=c(list(var=1:7),rep(NA,6),list(1)),
               y=weight_behavior[,1], refy = 0, predref = "F",nskip=0,ndpost=2)

## ----fig.show='hold'----------------------------------------------------------
summary(try2)

## ----echo=T, results='hide'---------------------------------------------------
try3= bma.bart(pred=weight_behavior[,c(1,4)], m=weight_behavior[,c(2:3,5:14)], 
               y=weight_behavior[,15], refy = 0, predref = "OTHER",nskip=0,ndpost=2)

## -----------------------------------------------------------------------------
summary(try3)

## -----------------------------------------------------------------------------
 try0$DIC$D_bar
 try0$DIC$var_D
 try0$DIC$p_D #calculated using two methods
 try0$DIC$DIC

## ----results='hide'-----------------------------------------------------------
#generating the data
set.seed(123)
n <- 1000
X <- rnorm(n)
M <- 0.5 + 0.5 * X^2 + rnorm(n)     
Y <- 0.6 * M^3 + 0.5* X + rnorm(n)  

data <- data.frame(X = X, M = M, Y = Y)
try4= bma.bart(pred=X, m=M, y=Y, ndpost=10L, nskip=100L)  #use the sd/10 for X and M
summary(try4)


## -----------------------------------------------------------------------------
summary(try4)
#Associations between X and M
plot(X, M, main = "BART Model: X vs M", ylab = "M", xlab = "", pch = 16, col = "blue") 
points(X, try4$m.models[[1]]$yhat.train.mean, col = "red", pch = 16)

ie2.apart=try4$apart.ie[,,1]  
plot(X,apply(ie2.apart,2,mean), main = "X vs apart", ylab = "apart",xlab = "", pch = 16, col = "red") 
#Association between M and Y
plot(M, Y, main = "BART Model: M vs Y", ylab = "Y", xlab = "", pch = 16, col = "blue") #
points(M, try4$y.model$yhat.train.mean, col = "red", pch = 16)
legend("topleft", legend = c("Observed", "BART Predictions"), col = c("blue", "red"), 
       pch = c(16,16), bty="n")

ie2.bpart=try4$bpart.ie[,,1]
plot(M, apply(ie2.bpart,2,mean), main = "M vs bpart", ylab = "bpart", xlab = "", col = "red", pch = 16)