CoxRFX {ebmstate}

Empirical Bayes, multi-state Cox model

Description

This function estimates a multi-state Cox model with one or more Gaussian priors imposed on the regression coefficients (see Therneau et al., 2003). Multiple groups of coefficients can be defined: coefficients within a group share the same (possibly unknown) mean and variance. The parameters and hyperparameters are efficiently estimated by an EM-type algorithm built around the function survival::coxph.

Usage

CoxRFX(
  Z,
  surv,
  groups = rep(1, ncol(Z)),
  which.mu = unique(groups),
  tol = 0.001,
  max.iter = 50,
  sigma0 = 0.1,
  sigma.hat = c("df", "MLE", "REML", "BLUP"),
  verbose = FALSE,
  ...
)

Arguments

Z

A data frame consisting of the covariate columns of a data set in 'long format', and two extra columns: one named 'trans', with the transition that each row refers to, and another named 'strata', with the stratum of each transition (transitions belonging to the same stratum are assumed to have the same baseline hazard function).

surv

A ‘survival’ object created with survival::Surv.

groups

A character or numeric vector whose ith element gives the group of the regression coefficient associated with the ith covariate column of Z (coefficients belonging to the same group share the same Gaussian prior).

which.mu

A vector with names or numbers of coefficient groups (see argument groups). If the name or number of a group of coefficients is given in this argument, CoxRFX will estimate the mean of its Gaussian distribution; otherwise the mean will be fixed at zero.

tol

Convergence criterium of the EM algorithm. The algorithm stops unless there is at least one parameter (or hyperparameter) for which it holds that the current estimate differs in absolute terms by more than tol from the previous estimate.

max.iter

The maximum number of iterations in the EM algorithm.

sigma0

A vector with the initial value of the variance hyperparameter for each group of coefficients. Or a single value, in case the initial value of the variance hyperparameter is meant to be the same for all groups.

sigma.hat

Which estimator to use for the variance hyperparameters (see details).

verbose

Gives more output.

...

Further arguments passed to the function survival::coxph.

Details

Different estimators exist for the variance hyperparameters: the default is "df", as used by Perperoglou (2014) and introduced by Schall (1991). Alternatives are MLE, REML, and BLUP, as defined by Therneau et al. (2003). Simulations suggest that the 'df' method is the most accurate.

The model can also be fitted using package coxme; the coxme routine numerically optimises the integrated partial likelihood, which may be more accurate, but is computationally expensive.

Value

An object of class c(coxrfx,coxph.penal,coxph), which is essentially a coxph object with a few extra fields [the inputs $groups, $Z and $surv, and the hyperparameters $sigma2 (variances) and $mu (means)]. See survival::coxph.object.

Author(s)

Moritz Gerstung & Rui Costa, extending the work of Terry Therneau et al. in the package survival.

References

Terry M Therneau, Patricia M Grambsch & V. Shane Pankratz (2003) Penalized Survival Models and Frailty, Journal of Computational and Graphical Statistics, 12:1, 156-175, http://dx.doi.org/10.1198/1061860031365

A. Perperoglou (2014). Cox models with dynamic ridge penalties on time-varying effects of the covariates. Stat Med, 33:170-80. http://dx.doi.org/10.1002/sim.5921

R. Schall (1991). Estimation in generalized linear models with random effects. Biometrika, 78:719-727. http://dx.doi.org/10.1093/biomet/78.4.719

See Also

Package survival survival::coxph.object; survival::Surv; package coxme.

Examples

# Fit an empirical Bayes Cox model using
# simulated, illness-death data from 250
# patients ('mstate_data_sample').

#load simulated data
data("mstate_data_sample")

# Set class of ‘mstate_data_sample’
class(mstate_data_sample)<-c("data.frame","msdata")

# add transition matrix as attribute
tmat<-mstate::transMat(x=list(c(2,3),c(4),c(),c()),
      names=c("health","illness","death",
     "death_after_illness"))
attr(mstate_data_sample,"trans")<-tmat 

# expand covariates by transition:
covariates.expanded<-mstate::expand.covs(
      mstate_data_sample,
      covs=names(mstate_data_sample)
      [!names(mstate_data_sample)%in%c("id","from",
      "to","trans","Tstart","Tstop","time","status",
      "strata")],append=FALSE)


# argument ‘Z’ of coxrfx
Z<-data.frame(covariates.expanded,
   trans=mstate_data_sample$trans,
   strata=mstate_data_sample$trans)

# argument ‘surv’ for a non-homogeneous 
# Markov model
surv<-survival::Surv(mstate_data_sample$Tstart,
           mstate_data_sample$Tstop,
           mstate_data_sample$status)

# argument ‘groups’ of coxrfx
groups<-paste0(rep("group", ncol(Z)-2),c("_1","_2","_3"))

#fit random effects model
coxrfx_object<-CoxRFX(Z,surv,groups)

#show point estimates
summary(coxrfx_object)