dynaTree (version 1.2-17)

dynaTrees: Fitting Dynamic Tree Models

Description

A function to initialize and fit dynamic tree models to regression and classification data by the sequential Monte Carlo (SMC) method of particle learning (PL)

Usage

dynaTree(X, y, N = 1000, model = c("constant", "linear", "class", "prior"),
         nu0s20 = c(0,0), ab = c(0.95, 2), minp = NULL, sb = NULL, 
	       nstart = minp, icept = c("implicit", "augmented", "none"), 
         rprop = c("luvar", "luall", "reject"), verb = round(length(y)/10))
dynaTrees(X, y, N = 1000,  R = 10, sub = length(y),
          model = c("constant", "linear", "class", "prior"), nu0s20 = c(0,0),
          ab=c(0.95, 2), minp = NULL, sb = NULL, nstart = minp,
          icept =  c("implicit", "augmented", "none"),
          rprop = c("luvar", "luall", "reject"), XX = NULL, yy = NULL,
	  varstats = FALSE, lhs = NULL, plotit = FALSE, proj = 1,
          rorder = TRUE, verb = round(sub/10), pverb=round(N/10), ...)

Value

Both functions return an object of class "dynaTree", which is a list containing the following fields

m

ncol(X)

T

nrow(X)

N

the number of particles used

X

a copy of the design matrix X

y

a copy of the responses y

model

a copy of the specified leaf model

params

a vector containing c(nu0s20, alpha, beta, minp, sb, icept, rprop), where the latter two are in integer form

verb

a copy of the verbosity argument

lpred

a vector of log posterior probabilities for each observation, conditional on the ones previous, for all time (2*minp):T; see getBF for calculating Bayes factors from these

icept

a copy of the intercept argument

time

the total computing time used to build the particle cloud

num

a “pointer” to the C-side particle cloud; see the note below

-

The dynaTrees function can obtain predictive samples (via predict.dynaTree) at each of the R

repeats. Therefore, the "dynaTree" object returned contains extra fields collecting these predictive samples, primarily comprising of R columns of information for each of the fields returned by predict.dynaTree; see that function for more details. Likewise, when varstats = TRUE the returned object also contains vpu, vpt and parde[ fields whose columns contain the varpropuse and

varproptotal outputs.

Likewise, dynaTrees, can provide variable usage summaries if varstats = TRUE, in which case the output includes

vpu and vpt fields; See varpropuse

and varproptotal for more details

The dynaTrees function does not return num since it does not leave any allocated particle clouds on the C-side

Arguments

X

A design matrix of real-valued predictors

y

A vector of length nrow(X) containing real-valued responses (for regression) or positive integer-valued class labels (for classification)

N

a positive scalar integer indicating the number of particles to be used

R

a scalar integer >= 2 indicating the number of “repeats” or passes through the data, as facilitated by dynaTrees; see details below

sub

Optional argument allowing only a subset of the length(y) X-y pairs to be used in each repeat of dynaTrees; each repeat will use a different random subset of size sub

model

indicates the type of model to be used at the leaves of the tree; "constant" and "linear" apply to regression, and "class" to multinomial classification; finally "prior" was recently added to explore sampled without data

nu0s20

a two-vector indicating Inverse Gamma prior parameters c(nu0, sigma20) for the variance in each leaf node, \(\sigma^2\). A c(0,0) setting indicates a default, scale-invariant, prior; does not apply to the "class" model

ab

tree prior parameter c(alpha, beta); see details below

minp

a positive scalar integer describing the smallest allowable region in the treed partition; if NULL (default) a suitable minimum is calculated based on dim(X) and the type of model being fit

sb

an optional two-vector of positive integers indicating c(splitmin, basemax) for the "linear" model. It gives the first column of X on which treed partitioning is allowed, and the last column of X to use as covariates in the linear model at the leaves, respectively

nstart

a positive scalar integer >= minp indicating the time index at which treed partitioning is allowed to start

icept

indicates the type of intertcept term used (only applies to model="linear"). The default, "implicit" causes the inputs X to be centered so the intercept can be implied as an afterthought; "augmented" causes the inputs X to automatically gain a leading column of ones in a way that is transparent to the user; and "none" assumes that no intercept is being used, or that the user has pre-treated X to have a column of ones. The main advantage of "implicit" over "augmented" is that the former can default to a constant model fit if leaf design matrices become rank deficient. The latter defaults to the zero-model in such cases

XX

a design matrix of predictive locations (where ncol(XX) == ncol(X)) for dynaTrees; also see predict.dynaTree and some explanation in the details below

yy

an optional vector of “true” responses at the XX predictive locations at which the log posterior probability are to be reported

varstats

if TRUE causes the varpropuse, varproptotal, and relevance.dynaTree functions to be called on after each repetition to collect the usage proportions of each input variable (column of X); see those documentation files for more details

lhs

an optional lhs argument to sens.dynaTree if a sensitivity analysis step is desired after each restart (XX="sens")

plotit

a scalar logical indicating if the fit should be plotted after each of the R repeats; only applies to 1-d data and dynaTrees

proj

when ncol(x$X) > 1 and plotit = TRUE this argument is passed to plot.dynaTree to make a 1-d projection using x$X[,proj]

rorder

a scalar logical indicating if the rows of X (and corresponding components of y) should be randomly re-ordered for repeats 2:R in order to assess the how the time-ordering of the SMC effects the Monte Carlo error; only applies to dynaTrees. Alternatively, one can specify an nrow(X)-by-(R-1) matrix of orderings (permutations of 1:nrow(X))

rprop

indicates the scheme used to construct a grow proposal. The best setting, "luall" uses the lower (L) and upper (U) rectangle method based on minp (above) as described in the seminal reference in which the growing location and dimension is sampled uniformly. It can be computationally intensive for large input spaces. A thriftier option (the default) in this case is "luvar" which uniformly chooses the splitting variable first and then uses the LU method marginally. Thriftier still is "reject" which just proposes uniformly in the bounding leaf rectangle and rejects subsequent grows that lead to partitions with too few data points; (see the minp argument)

verb

a positive scalar integer indicating how many time steps (iterations) should pass before a progress statement is printed to the console; a value of verb = 0 is quiet

pverb

a positive scalar integer indicating after many particles should be processed for prediction before a progress statement is printed to the console; a value of pverb = 0 is quiet

...

extra arguments to predict.dynaTree passed from dynaTrees

Author

Robert B. Gramacy ,
Matt Taddy and Christoforos Anagnostopoulos

Details

The dynaTree function processes the X and y pairs serially via PL. It builds up a particle cloud which is stored as an object in C. A “pointer” to that object is the primary return value. The dynaTrees function fits several (R) different dynamic tree models on different time-orderings of the data indices and also obtains samples from the posterior predictive distribution at new XX locations. These predictions can be averaged over each repeat, or used to assess the Monte Carlo predictive error.

Three different leaf models are supported: two for regression and one for classification. If model == "class" then the y values must contain representatives from every class (1:max(y)). For details of these models and the complete description of their use at the leaves of the dynamic trees, see the Taddy, et al., (2009) reference, below.

The tree prior is specified by ab=c(alpha, beta) via the and minp. It was originally described by Chipman et al., (1998, 2002) $$p_{\mbox{\tiny split}}(\eta, \mathcal{T}) = \alpha*(1+\eta)^\beta$$ and subsequently augmented to enforce a minimum number of points (minp) in each region.

Once a "dynaTree"-class object has been built (by dynaTree), predictions and estimates of sequential design and optimization criteria can be obtained via predict.dynaTree, a generic prediction method. These values can be used to augment the design, and the update.dynaTree function can be used to quickly update the fit with the augmenting data

References

Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586

Gramacy, R.B., Taddy, M.A., and S. Wild (2011). “Variable Selection and Sensitivity Analysis via Dynamic Trees with an Application to Computer Code Performance Tuning” arXiv:1108.4739

Carvalho, C., Johannes, M., Lopes, H., and Polson, N. (2008). “Particle Learning and Smoothing”. Discussion Paper 2008-32, Duke University Dept. of Statistical Science.

Chipman, H., George, E., & McCulloch, R. (1998). Bayesian CART model search (with discussion). Journal of the American Statistical Association, 93, 935--960.

Chipman, H., George, E., & McCulloch, R. (2002). Bayesian treed models. Machine Learning, 48, 303--324.

https://bobby.gramacy.com/r_packages/dynaTree/

See Also

predict.dynaTree, update.dynaTree, plot.dynaTree, deletecloud, copy.dynaTree, getBF, varpropuse, varproptotal, sens.dynaTree, relevance.dynaTree

Examples

## simple parabolic data
n <- 100
Xp <- sort(runif(n,-3,3))
Yp <- Xp + Xp^2 + rnorm(n, 0, .2)

## fit a piece-wise linear model
parab.fit <- dynaTree(Xp, Yp, model="linear")

## obtain predictions at a new set of locations
## and plot
parab.fit <- predict(parab.fit, XX=seq(-3, 3, length=100))
plot(parab.fit)

## try duplicating the object
parab.fit.copy <- copy(parab.fit)

## must delete the cloud or memory may leak
deletecloud(parab.fit); parab.fit$num <- NULL
## to delete all clouds, do:
deleteclouds()

## for more examples of dynaTree see update.dynaTree

## Motorcycle accident data
if(require("MASS")) {
  data(mcycle)
  Xm <- mcycle[,1]
  Ym <- mcycle[,2]
  XXm <- seq(min(mcycle[,1]), max(mcycle[,1]), length=100)
  
  R <- 2 ## use R >= 10 for better results
  ## small R is for faster CRAN checks
  ## fit constant model with R=2 repeats and predictions
  moto.fit <- dynaTrees(Xm, Ym, XX=XXm, R=R, plotit=TRUE)
  
  ## plot the averages
  plot(moto.fit, ptype="mean")
  
  ## clouds automatically deleted by dynaTrees
}

if (FALSE) {
## 2-d/3-class classification data
library(plgp)
library(tgp)
xx <- seq(-2, 2, length=20)
XX <- expand.grid(xx, xx)
X <- dopt.gp(125, Xcand=XX)$XX
C <- exp2d.C(X)

## fit a classification model with R=10 repeats, 
class.fit <- dynaTrees(X, C, XX=XX, model="class")

## for plot the output (no generic plotting available)
cols <- c(gray(0.85), gray(0.625), gray(0.4))
par(mfrow=c(1,2))
library(interp)

## plot R-averaged predicted class
mclass <- apply(class.fit$p, 1, which.max)
image(interp(XX[,1], XX[,2], mclass), col=cols,
      xlab="x1", ylab="x2", main="repeated class mean")
points(X)
## plot R-averaged entropy
ment <-  apply(class.fit$entropy, 1, mean)
image(interp(XX[,1], XX[,2], ment),
      xlab="x1", ylab="x2", main="repeated entropy mean")
}