Multivariate Survival Trees Package

Description

This package constructs trees for multivariate survival data using marginal and frailty models

Details

Decision trees require few statistical assumptions, handle a variety of data structures, and provide meaningful interpretations. There are several R packages that provide functions to construct survival trees (see rpart, partykit, and DStree); this package extends the implementation to multivariate survival data. There are two main approaches to analyzing correlated failure times. One is the marginal approach studied by authors Wei et al. (1989) and Liang et al. (1993). In the marginal model, the correlation is modeled implicitly using generalized estimating equations on the marginal distribution formulated by the Cox (1972) proportional hazards model. The other approach is the frailty model studied by Clayton (1978) and Clayton and Cuzick (1985). In the frailty model, the correlation is modeled explicitly by a multiplicative random effect called frailty, which corresponds to some common unobserved characteristics shared by all correlated times.

The construction of the tree adopts a modified CART procedure controlling for the correlated failure times. The procedure consists of three stages: growing the initial tree, pruning the tree, and selecting the best-sized subtree; details of these steps are described elsewhere (Fan et al. [2006], Su and Fan [2004], and Fan et al. [2009]). There are two methods for selecting the best-sized subtree. When the dataset is large, one may divide the dataset into a training sample to grow and prune the initial tree and a test sample to select the best-sized tree. When the dataset is small, one can resample the dataset to choose the best-sized subtree.

Author(s)

Xiaogang Su, Peter Calhoun, & Juanjuan Fan

Maintainer: Peter Calhoun <calhoun.peter@gmail.com>

References

Calhoun P., Su X., Nunn M., Fan J. (2018) Constructing Multivariate Survival Trees: The MST Package for R. Journal of Statistical Software, 83(12), 1–21.

Clayton D.G. (1978) A model for association in bivariate life tables and its application in epidemiologic studies of familial tendency in chronic disease incidence. Biometrika, 65(1), 141–151

Clayton D.G. and Cuzick J. (1985) Multivariate generalization of the proportional hazards model. Journal of the Royal Statistical Society Series A, 148(2), 82–108

Cox D.R. (1972) Regression models and life-tables (with discussion). Journal of the Royal Statistical Society Series B, 34(2), 187–220.

Fan J., Su X., Levine R., Nunn M., LeBlanc M. (2006) Trees for Correlated Survival Data by Goodness of Split, With Applications to Tooth Prognosis. Journal of American Statistical Association, 101(475), 959–967.

Fan J., Nunn M., Su X. (2009) Multivariate exponential survival trees and their application to tooth prognosis. Computational Statistics and Data Analysis, 53(4), 1110–1121.

Liang K.Y., Self S.G., Chang Y. (1993) Modeling marginal hazards in multivariate failure time data. Journal of the Royal Statistical Society Series B, 55(2), 441–453

Su X., Fan J. (2004) Multivariate Survival Trees: A Maximum Likelihood Approach Based on Frailty Models. Biometrics, 60(1), 93–99.

Su X., Fan J., Wang A., Johnson M. (2006) On Simulating Multivariate Failure Times. International Journal of Applied Mathematics & Statistics, 5, 8–18

Wei L.J., Lin D.Y., Weissfeld L. (1989) Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. Journal of the American Statistical Association, 84(408), 1065–1073

Multivariate Survival Trees

Description

Constructs trees for multivariate survival data using marginal and frailty models. A wrapper function that grows a large initial tree, prunes the tree, and selects the best sized tree.

Usage

MST(formula, data, test = NULL, weights_data, weights_test, subset,
  method = c("marginal", "gamma.frailty", "exp.frailty", "stratified", "independence"),
  minsplit = 20, minevents = 3, minbucket = round(minsplit/3), maxdepth = 10,
  mtry = NULL, distinct = TRUE, delta = 0.05, nCutPoints = 50,
  selection.method = c("test.sample", "bootstrap"),
  B = 30, LeBlanc = TRUE, min.boot.tree.size = 1,
  plot.Ga = TRUE, filename = NULL, horizontal = TRUE, details = FALSE, sortTrees = TRUE)

Arguments

Details

Marginal and frailty models are the two main ways to analyze correlated failure times. Let X_{ij} represent the covariate vector for the jth member in the ith cluster.

The marginal model uses the Cox (1972) proportional hazards model:

\lambda_{ij}(t|X_{ij})=\lambda_{0}(t) \exp(\beta \cdot I(X_{ij} \leq c))

where \lambda_{0}(t) is an unspecified baseline hazard function and I(\cdot) is the indicator function.

The gamma frailty model uses the proportional hazards model:

\lambda_{ij}(t|X_{ij}, w_{i})=\lambda_{0}(t) \exp(\beta \cdot I(X_{ij} \leq c)) w_{i}

where \lambda_{0}(t) is an unspecified baseline hazard function, I(\cdot) is the indicator function, and w_{i} is the frailty term for the ith cluster.

The exponential frailty model uses the proportional hazards model:

\lambda_{ij}(t|X_{ij}, w_{i})=\exp(\beta_{0} + \beta_{1} \cdot I(X_{ij} \leq c)) w_{i}

where I(\cdot) is the indicator function and w_{i} is the frailty term for the ith cluster.

For the marginal model, a robust logrank statistic is calculated for each covariate X and possible cutpoint c. The estimate of the score function and likelihood of \beta can be obtained assuming independence. However, the variance-covariance structure adjusts for the dependence using a sandwich-type estimator. The best split is the one with the largest robust logrank statistic.

For the frailty models, a score test statistic is calculated from the maximum integrated log likelihood for each covariate X and possible cutpoint c. The frailty term must follow some known positive distribution; one common choice is w_{i} \sim \Gamma(1/\nu, 1/\nu) where \nu represents an unknown variance. Note, the exponential frailty model replaces the baseline hazard function with a constant, yielding different score test statistics and typically computationally faster splits. The best split is the one with the largest score test statistic.

Stratified model grows a tree by minimizing the within-strata variation. This method should be used with care because the tree will not split on variables with a fixed value within each stratum. The independence model ignores the dependence and uses the logrank statistic as the splitting rule.

For continuous variables with many distinct cutpoints, the number of cutpoints considered can be reduced to percentiles. Using percentiles increases efficiency at the expense of less accuracy.

Growing the initial tree is done by splitting nodes (as described above) reiteratively until the maximum depth of the tree is reached or a small number of observations remain at terminal node. However, as the final tree model can be any subtree of the initial tree, the number of subtrees can become massive. A goodness-of-fit with an added penalty for the number of internal nodes is used to prune the trees (i.e. reduce the number of subtrees considered). The best-sized tree is selected by the largest goodness-of-fit with the added penalty using either the test sample or bootstrap samples.

Value

Note, the constparty object requires a constant fit from each terminal node. Thus, the predict and plot functions ignore the dependence, so users are recommended to fit their own model when making predictions (see example)

Warning

Error messages in the gamma frailty models sometimes occur when using the bootstrap method. Increasing minsplit may help fix these errors. The exponential frailty model can have problems for large, extremely unbalanced designs. Currently weights can only be applied to marginal and gamma frailty models.

Note

Code may take awhile to implement large datasets. To decrease computation time, user should use test sample (selection.method = "test.sample"). User can also split continuous variables based on percentiles (distinct = FALSE) at the expense of slightly less accuracy. Gamma frailty models are more computationally intensive

Author(s)

Xiaogang Su, Peter Calhoun, and Juanjuan Fan

References

Calhoun P., Su X., Nunn M., Fan J. (2018) Constructing Multivariate Survival Trees: The MST Package for R. Journal of Statistical Software, 83(12), 1–21.

Cox D.R. (1972) Regression models and life-tables (with discussion). Journal of the Royal Statistical Society Series B, 34(2), 187–220.

Fan J., Su X., Levine R., Nunn M., LeBlanc M. (2006) Trees for Correlated Survival Data by Goodness of Split, With Applications to Tooth Prognosis. Journal of American Statistical Association, 101(475), 959–967.

Fan J., Nunn M., Su X. (2009) Multivariate exponential survival trees and their application to tooth prognosis. Computational Statistics and Data Analysis, 53(4), 1110–1121.

Su X., Fan J. (2004) Multivariate Survival Trees: A Maximum Likelihood Approach Based on Frailty Models. Biometrics, 60(1), 93–99.

See Also

rpart

Examples

set.seed(186117)
data <- rmultime(N = 200, K = 4, beta = c(-1, 0.8, 0.8, 0, 0), cutoff = c(0.5, 0.3, 0, 0),
    model = "marginal.multivariate.exponential", rho = 0.65)$dat
test <- rmultime(N = 100, K = 4, beta = c(-1, 0.8, 0.8, 0, 0), cutoff = c(0.5, 0.3, 0, 0),
    model = "marginal.multivariate.exponential", rho = 0.65)$dat

#Construct Multivariate Survival Tree:
fit <- MST(formula = Surv(time, status) ~ x1 + x2 + x3 + x4 | id, data, test,
    method = "marginal", minsplit = 100, minevents = 20, selection.method = "test.sample")

(tree_final <- getTree(fit, 4))
plot(tree_final)

#Fit a model from the final tree
data$term_nodes <- as.factor(predict(tree_final, newdata = data, type = 'node'))
coxph(Surv(time, status) ~ term_nodes + cluster(id), data = data)

Tooth Loss Data

Description

Survival of teeth with various predictors.

Usage

data("Teeth")

Format

A data frame with 65,890 teeth on the following 56 variables.

x1

numeric. mobil Mobility score (on a scale 0–5).

x2

numeric. bleed Bleeding on Probing (percentage).

x3

numeric. plaque Plaque Score (percentage).

x4

numeric. pocket_mean Periodontal Probing Depth (tooth-level mean).

x5

numeric. pocket_max Periodontal Probing Depth (tooth-level mean).

x6

numeric. cal_mean Clinical Attachment Level (tooth-level mean).

x7

numeric. cal_max Clinical Attachment Level (tooth-level max).

x8

numeric. fgm_mean Free Gingival Margin (tooth-level mean).

x9

numeric. fgm_max Free Gingival Margin (tooth-level max).

x10

numeric. mg Mucogingival Defect.

x11

numeric. filled Filled Surfaces.

x12

numeric. decay_new Decayed Surfaces – new.

x13

numeric. decay_recur Decayed Surfaces – recurrent.

x14

numeric. dfs Decayed and Filled Surfaces.

x15

factor. crown Crown.

x16

factor. endo Endodontic Therapy.

x17

factor. implant Tooth Implant.

x18

factor. pontic Bridge Pontic.

x19

factor. missing_tooth Missing Tooth.

x20

factor. filled_tooth Filled Tooth.

x21

factor. decayed_tooth Decayed Tooth.

x22

factor. furc_max Furcation Involvement for Molars.

x23

numeric. bleed_ave Bleeding on Probing (mean percentage).

x24

numeric. plaque_ave Plaque Index (mean percentage).

x25

numeric. pocket_mean_ave Periodontal Probing Depth (mean of tooth mean).

x26

numeric. pocket_max_ave Periodontal Probing Depth (mean of tooth max).

x27

numeric. cal_mean_ave Clinical Attachment Level (mean of tooth mean).

x28

numeric. cal_max_ave Clinical Attachment Level (mean of tooth max).

x29

numeric. fgm_mean_ave Free Gingival Margin (mean of tooth max).

x30

numeric. fgm_max_ave Free Gingival Margin (mean of tooth max).

x31

numeric. mg_ave Mucogingival Defect (mean).

x32

numeric. filled_sum Filled Surfaces (total).

x33

numeric. filled_ave Filled Surfaces (mean).

x34

numeric. decay_new_sum New Decayed Surfaces (total).

x35

numeric. decay_new_ave New Decayed Surfaces (mean).

x36

numeric. decay_recur_sum Recurrent Decayed Surfaces (total).

x37

numeric. decay_recur_ave Recurrent Decayed Surfaces (mean).

x38

numeric. dfs_sum Decayed and Filled Surfaces (total).

x39

numeric. dfs_ave Decayed and Filled Surfaces (mean).

x40

numeric. filled_tooth_sum Number of Filled Teeth.

x41

numeric. filled_tooth_ave Percentage of Filled Teeth.

x42

numeric. decayed_tooth_sum Number of Decayed Teeth.

x43

numeric. decayed_tooth_ave Percentage of Decayed Teeth.

x44

numeric. missing_tooth_sum Number of Missing Teeth.

x45

numeric. missing_tooth_ave Percentage of Missing Teeth.

x46

numeric. total_tooth Number of Teeth.

x47

numeric. dft Number of Decayed and Filled Teeth.

x48

numeric. baseline_age Patient Age at Baseline (years).

x49

factor. gender Gender.

x50

factor. diabetes Diabetes Mellitus.

x51

factor. tobacco_ever Tobacco Use.

molar

logical. Molar.

id

numeric. Patient ID.

tooth

numeric. Tooth ID.

event

numeric. Tooth Loss Status.

time

numeric. Follow Up Time.

Details

Patients were treated at the Creighton University School of Dentistry from August 2007 to March 2013. This is a subset of the original data.

The goal is to estimate the survival time of teeth (molars or non-molars) using 51 predictors (22 tooth-level factors (x1–x22) and 29 patient-level factors (x23–x51)).

Examples

data(Teeth)

Extract initial or best-sized tree

Description

This function extracts the tree based on the split penalty.

Usage

getTree(mstObj, Ga = c("0", "2", "3", "4", "log_n"))

Arguments

Value

The tree of object class "constparty"

Author(s)

Peter Calhoun <calhoun.peter@gmail.com>

See Also

MST

Random Multivariate Survival Data

Description

Generates multivariate survival data

Usage

rmultime(N = 100, K = 4, beta = c(-1, 2, 1, 0, 0), cutoff = c(0.5, 0.5, 0, 0),
  digits = 1, icensor = 1, model = c("gamma.frailty", "log.normal.frailty",
  "marginal.multivariate.exponential", "marginal.nonabsolutely.continuous",
  "nonPH.weibull"), v = 1, rho = 0.65, a = 1.5, lambda = 0.1)

Arguments

Details

This function generates multivariate survival data. Letting i=1,...,N number of clusters, j=1,...,K number of units per cluster, and X_{ij} be a candidate covariate, the following multivariate survival models can be used:

gamma.frailty: \hspace{2mm} \lambda_{ij}(t)=\exp(\beta_{0}+\beta_{1} \cdot I(X_{ij} \leq c)) w_{i} with w_{i} \sim \Gamma(1/v, 1/v)

log.normal.frailty: \hspace{2mm} \lambda_{ij}(t)=\exp(\beta_{0}+\beta_{1} \cdot I(X_{ij} \leq c) + w_{i}) with w_{i} \sim N(0, v)

marginal.multivariate.exponential: \hspace{2mm} \lambda_{ij}(t)=\exp(\beta_{0}+\beta_{1} \cdot I(X_{ij} \leq c)) absolutely continuous

marginal.nonabsolutely.continuous: \hspace{2mm} \lambda_{ij}(t)=\exp(\beta_{0}+\beta_{1} \cdot I(X_{ij} \leq c)) not absolutely continuous

nonPH.weibull: \hspace{2mm} \lambda_{ij}(t)=\lambda_{0}(t) \exp(\beta_{0}+\beta_{1} \cdot I(X_{ij} \leq c)) w_{i} with w_{i} \sim \Gamma(1/v ,1/v) and

\hspace{96mm} \lambda_{0}(t)=\alpha \lambda t^{\alpha-1}

The user specifies the coefficients (\beta_{0} and \beta_{1}), the cutoff values, the censoring rate, and the model with the respective parameters.

Value

Author(s)

Xiaogang Su, Peter Calhoun, Juanjuan Fan

References

Fan J., Nunn M., Su X. (2009) Multivariate exponential survival trees and their application to tooth prognosis. Computational Statistics and Data Analysis, 53(4), 1110–1121.

Su X., Fan J., Wang A., Johnson M. (2006) On Simulating Multivariate Failure Times. International Journal of Applied Mathematics & Statistics, 5, 8–18

See Also

genSurv, complex.surv.dat.sim, survsim

Examples

randMarginalExp <- rmultime(N = 200, K = 4, beta = c(-1, 2, 2, 0, 0), cutoff = c(0.5, 0.5, 0, 0),
    digits = 1, icensor = 1, model = "marginal.multivariate.exponential", rho = .65)$dat

randFrailtyGamma <- rmultime(N = 200, K = 4, beta = c(-1, 1, 3, 0), cutoff = c(0.4, 0.6, 0),
    digits = 1, icensor = 1, model = "gamma.frailty", v = 1)$dat