hrf: Random forest for longitudinal data

Description

Fits a random forest of historical regression trees to longitudinal data.

Usage

hrf(x,
    time=NULL,
    id=NULL,
    yindx,
    ntrees = 100,
    method="freq",
    mtry=NULL,
    se=FALSE,
    B=100,
    R=10,
    nsamp=5,
    historical=TRUE,
    vh=NULL,
    vc=NULL,
    delta=NULL,
    classify=FALSE,
    control=list())

Arguments

A data frame containing response and predictors

time

A vector of observation times associated with rows of x.

Subject identifier, if NULL observations are assumed independent

yindx

Name of response variable, alt. its column number in x

ntrees

Number of trees in ensemble

method

Historical summary method, can be freq, frac, mean0, freqw, fracw and mean0w (see below).

mtry

Number of predictors sampled at each split

If TRUE then bootstrap standard errors are computed. Total number of trees fit for bootstrapping is B*R.

Only used if se=TRUE, number of bootstrap samples, defaults to 100.

Only used if se=TRUE, forest size for each bootstrap sample, defaults to R=10.

nsamp

Number of sampled delta values, see below

historical

If TRUE then historical splitting is done, else only standard (ie concurrent predictor) splitting.

Optional vector of variable names to be used as historical predictors, can be a numeric vector giving column numbers of x.

Optional vector of variable names to be used as concurrent predictors, can be a numeric vector giving column numbers of x.

delta

An optional vector of time-lags to be used (see below).

classify

If TRUE then a classification tree is built (using gini-impurity node splitting).

control

A list of control parameters (see below). All arguments, except those describing the data, can be set in control. Arguments in control are used if both are given.

Value

Returns a list with elements, the most important being: h is a list returned by parLapply for fitting trees in parallel, error gives the OOB error (mse for regression, misclassification rate for classification), control a list containing control arguments, boot gives bootstrapped model estimates (if se=TRUE), x,id,time give the training data.

Details

The hrf 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 vector of predictors at time $t_{ij}$ are $x_{ij}$. The number of observations can vary across subjects, and the sampling in time can be irregular.

hrf 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 (i.e. prior to 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 number (and observation times) 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 using the best split among all splits of concurrent and historical predictors.

Different summary functions are available in hrf, specified by the argument method. If method="freq" the summary function summarizing a covariate history is: $$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="freq". 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 within-subject 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. See below on control for futher arguments governing the historical splitting. See below on control for futher arguments governing the historical splitting.

Setting se=TRUE performs standard error estimation. The number of bootstrap samples (sampling subjects with replacement) is determined by B. For each bootstrap sample a random forest with R trees is built, which defaults to R=10. The bias induced by using smaller bootstrap ensemble sizes is corrected for in the estimate. Using se=TRUE will influence summaries from the fitted model, such as providing approximate confidence intervals for partial dependence plots (when running partdep_hrf), and give standard errors for predictions when predict_hrf is used.

All arguments (except those decribing the data x, yindx,time and id) can be set in the control list. The arguments supplied in control are used if both are supplied. So if ntrees=300 and control=list(ntrees=500) then 500 trees are fit. Besides the arguments described above, a number of other parameters can be set in control. These are: nodesize giving the minimum number of training observations in a terminal node; sample_fraction giving the fraction of data sample to train each tree; dtry the number of sampled delta values used when splitting based on a historical variable (note this is alternatively controlled by the argument nsamp above ; ndelta the number of delta values to use if delta is not supplied, these are taken as the quantiles from the distribution of observed delta values in the data; qtry the number of sampled values of $\eta_2$ for method=freq/frac, $\eta_3$ for method=freqw/fracw; quantiles is a vector of probabilities, and is used when method=frac or method=fracw, ie when covariate histories are split on their windowed empirical distributions. Splitting is restricted to these quantiles.

References

L. Breiman (2001). “Random Forests,” Machine Learning 45(1):5-32.

Zhang and Singer (2010) “Recursive Partioning and Applications” Springer.

Sexton and Laake (2009) “Standard errors for bagged and random forest estimators,” Computational Statistics and Data Analysis.

Examples

Run this code

# NOT RUN {
# }
# NOT RUN {
# ----------------------------------------------------------------------------------------- ##
# Mother's stress on child illness:
# 	Investigate whether mother's stress is (Granger) causal for 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))


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

control=list(vh=historical_predictors,vc=concurrent_predictors)


# -- fit historical random forest
#	(NOTE: response is 0/1 so a regression tree is
# 	 the same as  a classification tree with Gini-index splitting) 
ff=hrf(x=mscm,id=mscm$id,time=mscm$day,yindx="illness",control=control)

# out-of-bag error 
plot(1:length(ff$error),ff$error,type="l",main="OOB error",xlab="forest size",ylab="mse")

# .. larger nodesize works slightly better
control$nodesize=20
ff=hrf(x=mscm,id=mscm$id,time=mscm$day,yindx="illness",control=control)
points(1:length(ff$error),ff$error,type="l",col="blue")


# -- variable importance table 
vi=varimp_hrf(ff)
vi



# -- fit historical random forest with 'se=TRUE'
control$se=TRUE
ff=hrf(x=mscm,id=mscm$id,time=mscm$day,yindx="illness",control=control)

# -- partial dependence for top 4 predictors (with +/-2 SE estimates) 
par(mfrow=c(2,2))
for(k in 1:4)
	pd=partdep_hrf(ff,xindx=as.character(vi$Predictor[k]))

par(mfrow=c(1,1))




## -- Classification trees 

## setting classify=TRUE builds classification tree (gini-impurity node splitting)
control$classify=TRUE
## ... standard error estimation not implemented .. turn off bootstrapping 
control$se=FALSE

ff=hrf(x=mscm,id=mscm$id,time=mscm$day,yindx="illness",control=control)

# -- plot oob classification error 
plot(1:length(ff$error),ff$error,type="l",xlab="forest size",ylab="oob classification error")
abline(mean(mscm$illness),0,lty=2)  ## error of constant model

p=predict_hrf(ff)

## variable importance table (model error measured by gini-impurity)
vi=varimp_hrf(ff)
vi


# -------------------------------- #
# Data w/irregular observation times 
# ------------------------------- #
data(cd4)

control=list(se=TRUE)
ff=hrf(x=cd4,id=cd4$id,time=cd4$time,yindx="count",control=control)

vi=varimp_hrf(ff)

# -- partial dependence for top 4 predictors (with +/-2 SE estimates) 
par(mfrow=c(2,2))
for(k in 1:4)
	pd=partdep_hrf(ff,xindx=as.character(vi$Predictor[k]))
par(mfrow=c(1,1))

plot(1:length(ff$error),ff$error,xlab="forest size",ylab="oob mse",type="l")

## by default, the number of delta values (parameter 'eta_1' above) is 20
## can set this using 'ndelta'  
control$ndelta=50

control$se=FALSE # -- turning off bootstrapping .. 
ff=hrf(x=cd4,id=cd4$id,time=cd4$time,yindx="count",control=control)
points(1:length(ff$error),ff$error,type="l",lty=2)

# the grid of delta values 
ff$control$delta


# --------------------------------------- ##
# Boston Housing data (not longitudinal)
# --------------------------------------- ##
# library(htree)
library(mlbench)
library(randomForest)

data(BostonHousing)
dat=as.data.frame(na.omit(BostonHousing))


## omitting arguments time/id assumes rows are iid 
control=list(ntrees=500,sample_fraction=.5,nodesize=1)
h=hrf(x=dat,yindx="medv",control=control)

## randomForest comparison 
##	(by default, randomForest samples with replacement, while hrf samples without) 
r=randomForest(medv~.,data=dat,replace=F,sampsize=ceiling(.5*nrow(dat)),nodesize=1)

## plot oob-error for both
plot(1:length(r$mse),r$mse,type="l",ylim=c(min(r$mse,h$error),max(r$mse,h$error)),
	main="BostonHousing",xlab="forest size",ylab="out-of-bag mse")
points(1:length(h$error),h$error,type="l",col="blue")

## -- variable importance table
vi=varimp_hrf(h)
vi

## -- partial dependence plots with approximate 95<!-- % C.I  -->
control$se=TRUE
h=hrf(x=dat,yindx="medv",control=control)

par(mfrow=c(2,2))
for(k in 1:4)
	pd=partdep_hrf(h,xindx=as.character(vi$Predictor[k]))

par(mfrow=c(1,1))



# }

Run the code above in your browser using DataLab