BayesTree package plot extraction in R

Question

I am using the BayesTree package in R. Using the author's example:

##simulate data (example from Friedman MARS paper)
f = function(x){
10*sin(pi*x[,1]*x[,2]) + 20*(x[,3]-.5)^2+10*x[,4]+5*x[,5]
}
sigma = 1.0 #y = f(x) + sigma*z , z~N(0,1)
n = 100 #number of observations
set.seed(99)
x=matrix(runif(n*10),n,10) #10 variables, only first 5 matter
Ey = f(x)
y=Ey+sigma*rnorm(n)

##run BART
set.seed(99)
bartFit = bart(x,y)
plot(bartFit) # plot bart fit

When the plot(bartFit) is executed it produces two plots "The plot method sets mfrow to c(1,2) and makes two plots." I am looking to only use the second plot produced from this command. Is there an easy way to extract only the second plot?

This is provided as well:

## S3 method for class 'bart':
plot((
   x,
   plquants=c(.05,.95), cols =c('blue','black'),
   ...))

I am relatively new to R, in general if a command gives more then 1 plot how can you choose a specific plot you would like? For example in the linear model (lm) plot command you also get multiple plots.

Answer 1

Ok, I don't have a solution, but maybe I'm onto something. First, I tried to figure which list elements from bartFit are plotted on the second plot. You can look at the structure of the bart using str(bartFit). It seems that second plot is created from $y and $yhat.train.mean, which seems right when you look at it:

plot(bartFit$yhat.train.mean ~ bartFit$y)

but it's obviously missing the bars. So I decided to change the strategy and traceback() to see how those plots are created. No luck, until I changed the name of x-scale to force an error.

plot(bartFit, xlab="something")

Which gave me only the first plot, error message, but also ability to traceback(). There, I could see that the second plot is created using:

4: plot.default(x$y, qm, ylim = range(ql, qu), xlab = "y", ylab = "posterior interval for E(Y|x)", ...)
3: plot(x$y, qm, ylim = range(ql, qu), xlab = "y", ylab = "posterior interval for E(Y|x)", ...)

but that's not good. I have no idea what is qm, or where ylim ranges come from. So I started playing with par. There I discovered very, very crude solution to use par argument mfg to "force" plotting only the second plot. This will still leave empty space for one plot, but won't plot it. But it behaves strangely as well. Tricky one.

plot(bartFit)
par(mfg=c(1,2))
plot(bartFit)