2. Examples

To use OpenBT in R, install Ropenbt as described in Getting Started with R. This example assumes that the command line tools were built with MPI support.

Let’s create a test function. A popular one is the Branin function:

# Test Branin function, rescaled
braninsc <- function(xx)
{
  x1 <- xx[1]
  x2 <- xx[2]

  x1bar <- 15*x1 - 5
  x2bar <- 15 * x2

  term1 <- x2bar - 5.1*x1bar^2/(4*pi^2) + 5*x1bar/pi - 6
  term2 <- (10 - 10/(8*pi)) * cos(x1bar)

  y <- (term1^2 + term2 - 44.81) / 51.95
  return(y)
}


# Simulate Branin data for testing
set.seed(99)
n=500
p=2
x = matrix(runif(n*p),ncol=p)
y=rep(0,n)
for(i in 1:n) y[i] = braninsc(x[i,])

And then we can load the Ropenbt package and fit a BART model. Here we set the model type as model="bart", which ensures that we fit a homoscedastic BART model. The number of MPI processes to use is specified as tc=4. For a list of all optional parameters, see args(openbt).

library(Ropenbt)
fit=openbt(x,y,tc=4,model="bart",modelname="branin")

Next we can construct predictions and make a simple plot. Here, we are calculating the in-sample predictions since we passed the same x matrix to the predict.openbt() function.

# Calculate in-sample predictions
fitp=predict.openbt(fit,x,tc=4)

# Make a simple plot
plot(y,fitp$mmean,xlab="observed",ylab="fitted")
abline(0,1)

To save the model, use the openbt.save() function. Similarly, load the model using openbt.load(). Because the posterior can be large in sample-based models such as these, the fitted model is saved in a compressed file format with the extension .obt.

# Save fitted model as test.obt in the working directory
openbt.save(fit,"test")

# Load fitted model to a new object.
fit2=openbt.load("test")

The standard variable activity information, calculated as the proportion of splitting rules involving each variable, can be computed using the vartivity.openbt() function.

# Calculate variable activity information
fitv=vartivity.openbt(fit2)

# Plot variable activity
plot(fitv)

A more accurate alternative is to calculate the Sobol’ indices.

# Calculate Sobol' indices
fits=sobol.openbt(fit2)
fits$msi
fits$mtsi
fits$msij