htree (version 1.0.0)

htb: Tree boosting for longitudinal data

Description

Fits a boosted ensemble of historical regression trees to longitudinal data.

Usage

htb(x,
    time=NULL,
    id=NULL,
    yindx,
    ntrees = 100,
    method="freqw",
    nsplit=1,
    lambda=.1,
    family="gaussian",
    cv.fold=0,
    cv.rep=NULL,
    nsamp=5,
    historical=TRUE,
    keep_data=TRUE,
    vh=NULL,
    vc=NULL,
    delta=NULL)

Arguments

x

a data frame containing response and predictors

time

vector of observation times

id

subject identifier, if NULL observations are assumed independent

yindx

Column number in x corresponding response variable

ntrees

Number of trees in ensemble

method

Historical summary method, can be freq, frac, mean or sum

nsplit

Number of splits in each regression tree.

lambda

Shrinkage parameter.

family

Either "gaussian" (default) or "bernoulli".

cv.fold

Number of cross-validation folds, if cv<=1 no cross-validation is run.

cv.rep

Number of times to repeat the cross-validation. If not given set to cv.fold.

historical

If TRUE then historical splitting is done, else standard splitting.

nsamp

Number of sampled delta values, see below

keep_data

If TRUE training data is returned in fit object.

vh

Optional vector of indexes giving the historical predictors.

vc

Optional vector of indexes giving strictly concurrent effects.

delta

A vector of history lags to be used (see below), defaults to NULL in which case all possible observed lags are available for splitting.

Value

Returns a list with elements: trees giving the tree ensemble, cv (if cv.fold>1) giving cross-validation model estimates, cv_error cross-validated error, as well as arguments of the function call.

Details

The htb function fits a random forest model to longitudinal data. Data is assumed to be of form: \(z_{ij}=(y_{ij},t_{ij},x_{ij})\) for \(i=1,..,n\) and \(j=1,..,n_i\), with \(y_{ij}\) being the response for the \(i\)-th subject at the \(j\)-th observation time \(t_{ij}\). The predictors at time \(t_{ij}\) are \(x_{ij}\). The number of observations can vary across subjects, and the sampling in time can be irregular.

htb estimates a model for the response \(y_{ij}\) using both \((t_{ij},x_{ij})\) (the observations concurrent with \(y_{ij}\)) and all preceeding observations of the \(i\)-th subject up to (but not including) time \(t_{ij}\). The model is fit using historical regression (alt. classification) trees. Here a predictor is one of two types, either concurrent or historical. The concurrent predictors for \(y_{ij}\) are the elements of the vector (\((t_{ij},x_{ij})\)), while a historic predictor is the set of all preceeding values (preceeding time \(t_{ij}\)) of a given element of \((y_{ij},t_{ij},x_{ij})\), for subject \(i\). In a historic regression tree, node splitting on a concurrent predictor follows the approach in standard regression (classification) trees. For historical predictors the splitting is modified since, associated with each observed response \(y_{ij}\), the set of observations of a historical predictor will vary according to \(i\) and \(j\). For these, the splitting is done by first transforming the preceeding values of a predictor using a summary function. This summary is invertible, in the sense that knowledge of it is equivalent to knowledge of the covariate history. Letting \(\bar{z}_{ijk}\) denote the set of historical values of the \(k\)-th element of \(z_{ij}\), the summary function is denoted \(s(\eta;\bar{z}_{ijk})\) where \(\eta\) is the argument vector of the summary function. Node splitting based on a historical predictor is done by solving $$\mbox{argmin}\sum_{(ij)\in Node} (y_{ij}-\mu_L I(s(\eta;\bar{z}_{ijk})<c)-\mu_R I(s(\eta;\bar{z}_{ijk})\geq c))^2,$$ where the minimization is over the vector \((k,\mu,c,\eta)\). Each node of historical regression tree is split according to the best split among all possible splits based on both concurrent and historical predictors.

Different summary functions are available in hrf, specified by the argument method, where method="freq" and corresponds the summary $$s(\eta;\bar{z}_{ijk})=\sum_{h: t_{ij}-\eta_1\leq t_{ih}<t_{ij}} I(z_{ihk}<\eta_2)$$; method="frac": $$s(\eta;\bar{z}_{ijk})=\sum_{h: t_{ij}-\eta_1\leq t_{ih}<t_{ij}} I(z_{ihk}<\eta_2)/n_{ij}(\eta)$$; method="mean0": $$s(\eta;\bar{z}_{ijk})=\sum_{h: t_{ij}-\eta_1\leq t_{ih}<t_{ij}} z_{ihk}/n_{ij}(\eta)$$; method="freqw": $$s(\eta;\bar{z}_{ijk})=\sum_{h: t_{ij}-\eta_1\leq t_{ih}<t_{ij}-\eta_2} I(z_{ihk}<\eta_3)$$; method="fracw": $$s(\eta;\bar{z}_{ijk})=\sum_{h: t_{ij}-\eta_1\leq t_{ih}<t_{ij}-\eta_2} I(z_{ihk}<\eta_3)/n_{ij}(\eta)$$; method="meanw0": $$s(\eta;\bar{z}_{ijk})=\sum_{h: t_{ij}-\eta_1\leq t_{ih}<t_{ij}-\eta_2} z_{ihk}/n_{ij}(\eta)$$. Here \(n_{ij}(\eta)\) denotes the number of observations of subject \(i\) in the time window \([t_{ij}-\eta_1,t_{ij}-\eta_2)\). In the case \(n_{ij}(\eta)=0\), the summary function is set to zero, i.e \(s(\eta;\bar{z}_{ijk})=0\). The default is method="freqw". The possible values of \(\eta_1\) in the summary function can be set by the argument delta. If not supplied, the set of possible values of \(\eta_1\) is determined by the difference in time between successive observations in the data. When a split is attempted on a historical predictor, a sample of this set is taken from which the best split is selected. The size of this set equals that of the nsamp argument.

When cv.fold>1 then cross-validation is performed. In subsequent summaries of the model, say the partial dependence plots from partdep_htb, these cross-validation model fits are used to estimate the standard error. This is done using the delete-d jackknife estimator (where deletion refers to subjects, instead of rows of the training data). Each cross-validation model fit is performed by removing a random sample of 1/cv.fold of the subjects. The number of cross-validation runs is determined by cv.rep which defaults to the value of cv.fold.

References

J.H. Friedman (2001). “Greedy Function Approximation: A Gradient Boosting Machine,” Annals of Statistics 29(5):1189-1232.

J.H. Friedman (2002). “Stochastic Gradient Boosting,” Computational Statistics and Data Analysis 38(4):367-378.

See Also

predict_htb, partdep_htb, varimp_htb.

Examples

# NOT RUN {
# }
# NOT RUN {

# ----------------------------------------------------------------------------------------- ##
# Mother's stress on child illness:
# 	Investigate whether mother's stress is (Granger) causal on child illness 
#	'hrf' model is fit using previous observations of mother's stress to predict 
#	child's illness at given time point, but not mother's stress at that time point
#
#	Predictor variables are classified into "historical" and "concurrent"
#
#	A predictor is "historical" if its prior realizations can be used to predict 
#	the outcome.   
#
#	A predictor is "concurrent" if its realization at the same timepoint as the outcome
#	can be used to predict the outcome at that timepoint
#
#	A predictor can be both "concurrent" and "historical", the status of the predictors 
#	can be set by the 'vh' and 'vc' arguments of 'hrf'. 
#	(if not set these are automatically determined) 
#	 
# ------------------------------------------------------------------------------------------- ##


data(mscm) 
mscm=as.data.frame(na.omit(mscm))

## first 2 weeks
mscm=mscm[mscm$day<=14,] 

# -- set concurrent and historical predictors 
historical_predictors=match(c("illness","stress"),names(mscm))
concurrent_predictors=which(names(mscm)!="stress")


# logistic regression
ff=htb(x=mscm,id=mscm$id,time=mscm$day,yindx="illness",ntrees=100,lambda=.5,method="freqw",
	vh=historical_predictors,vc=concurrent_predictors,cv.fold=10,family="bernoulli")

# cross-validated negative log-likelihood
plot(1:ff$ntrees,ff$cv_error,type="l",col="blue",ylab="",
	xlab="iterations",main="Cross-validated deviance")

vi=varimp_htb(ff)
vi

par(mfrow=c(2,2))
for(k in 1:4)
	partdep_htb(ff,xindx=rownames(vi)[k])

par(mfrow=c(1,1))

# -- random intercept example --- #

p=5;sigma_e=.5;sigma_a=.5;v=rep(1,p);n=500;pnoise=2
random_intercept=rnorm(n,sd=sigma_a^.5)
random_intercept=as.numeric(matrix(random_intercept,nrow=p,ncol=n,byrow=TRUE))
x=random_intercept+rnorm(n*p,sd=sigma_e^.5)
id=sort(rep(1:n,p))
time<-rep(1:p,n)
znoise=matrix(rnorm(n*p*pnoise),ncol=pnoise)
xx=cbind(time,x,znoise)

# fit historical random forest
ff=htb(x=xx,time=time,id=id,yindx=2,ntrees=100,nsamp=5,lambda=.1,cv.fold=10)

# plot cv-error 
plot(1:ff$ntrees,ff$cv_error,type="l",col="blue",xlab="iterations",
	ylab="",main="Cross-validated mse")


# plot partial dependence of response on its past 
pd=partdep_htb(ff,xindx=2,ngrid=25,subsample=.1)

# variable importance table
varimp_htb(ff,nperm=20)




# ------------------------------------------ #
# Non-longitudinal data examples
# ------------------------------------------ #


# ------------------------------------------ #
# Boston Housing data 
# ------------------------------------------ #

library(htree)
library(mlbench)
data(BostonHousing)
dat=as.data.frame(na.omit(BostonHousing))


# omitting arguments 'time' and 'id' assumes rows are iid 
h=htb(x=dat,yindx="medv",ntrees=500,cv.fold=10)

# -- plot cross-validated Mean-squared error --- #
plot(1:length(h$cv_error),h$cv_error,type="l",xlab="Boosting iterations",
	ylab="MSE",main="Cross-validated Mean-squared error")

# -- iteration with lowest error
iter_hat=order(h$cv_error)[1]

# -- variable importance table
vi=varimp_htb(h,nperm=20)
vi 

# -- partial dependence of top 4 predictors with +/-2 S.E.
#	(standard errors are based on the 
#	delete-d jack-knife estimator, using the cross-validation runs)
par(mfrow=c(2,2))
for(k in 1:4)
	partdep_htb(h,xindx=rownames(vi)[k])


par(mfrow=c(1,1))

# }
# NOT RUN {
# }