How to simulate survival data with the sim.survdata function

Jonathan Kropko and Jeffrey J. Harden

2019-04-12

1 Introduction

Simulating data with known properties is an essential step in the development of new statistical methods. Simulating survival data, however, is more challenging than most simulation tasks. To describe this difficulty and our approach to overcoming it, we quote from Harden and Kropko (2018):

Typical methods for generating . . . durations often employ known distributions - such as the exponential, Weibull, or Gompertz — that imply specific shapes for the baseline hazard function. This approach is potentially problematic because it contradicts a key advantage of the Cox model — the ability to leave the distribution of the baseline hazard unspecified. By restricting the simulated data to a known (usually parametric) form, these studies impose an assumption that is not required in applications of the Cox model. In applied research, the baseline hazard in data can exhibit considerable heterogeneity, both across the many fields which employ the Cox model and within a given field . . . . Thus, simulating durations from one specific parametric distribution may not adequately approximate the data that many applied researchers analyze, reducing the simulation’s generalizability.

Here we address these problems by introducing a novel method for simulating durations without specifying a particular distribution for the baseline hazard function. Instead, the method generates — at each iteration of the simulation — a unique baseline hazard by fitting a cubic spline to randomly-drawn points. This produces a wide variety of shapes for the baseline hazard, including those that are unimodal, multimodal, monotonically increasing or decreasing, and other shapes. The method then generates a density function based on each baseline hazard and draws durations accordingly. Because the shape of the baseline hazard can vary considerably, this approach matches the Cox model’s inherent flexibility. By remaining agnostic about the distribution of the baseline hazard, our method better constructs the assumed data generating process (DGP) of the Cox model. Moreover, repeating this process over many iterations in a simulation yields more heterogeneity in the simulated samples of data, which matches the fact that applied researchers analyze a wide variety of data with the Cox model. This increases the generalizability of the simulation results.

In this vignette we demonstrate the uses of the sim.survdata() function for generating survival data. To begin, we load the coxed library:

library(coxed)

2 Simulating survival data using the flexible-hazard method

The flexible-hazard method described by Harden and Kropko (2018) first generates a baseline failure CDF: it plots points at (0, 0) and (T+1, 1), and it plots knots additional points with x-coordinates drawn uniformly from integers in [2, T] and y-coordinates drawn from U[0, 1]. It sorts these coordinates in ascending order (because a CDF must be non-decreasing) and if spline=TRUE it fits a spline using Hyman’s (1983) cubic smoothing function to preserve the CDF’s monotonicity. Next it constructs the failure-time PDF by computing the first differences of the CDF at each time point. It generates the survivor function by subtracting the failure CDF from 1. Finally, it computes the baseline hazard by dividing the failure PDF by the survivor function. Full technical details are listed in Harden and Kropko (2018).

2.1 Simulating a single dataset

To generate a survival dataset, use the sim.survdata() function, with the num.data.frames argument set to 1. Here we generate a single survival dataset with 1000 observations, in which durations can fall on any integer between 1 and 100:

2.2 Attributes of a simulation

The simulated data has several important attributes:

## $names
## [1] "data"             "xdata"            "baseline"        
## [4] "xb"               "exp.xb"           "betas"           
## [7] "ind.survive"      "marg.effect"      "marg.effect.data"
## 
## $class
## [1] "simSurvdata"

The full survival dataset, including the durations (y), the censoring variable (failed), and the covariates, is contained in the data frame object stored in the data attribute:

##             X1          X2          X3  y failed
## 1  -0.07392768 -0.16681023  0.95278612 85   TRUE
## 2  -0.44604373  0.76796254  0.22352026 19   TRUE
## 3   0.12489742 -0.57420780 -0.99545858 61   TRUE
## 4  -0.71203795  0.04409812  1.78181886 68   TRUE
## 5   1.79128542 -0.92812647 -1.24928618 20  FALSE
## 6   0.06870889  1.62915198 -0.76242481 67   TRUE
## 7  -0.54569924  1.66808802 -1.72180296 91   TRUE
## 8  -0.03863737 -0.44085091  0.45146624 85   TRUE
## 9   0.90402590  0.46646237 -1.23922206 69  FALSE
## 10  0.33471996 -0.77802688 -0.06973362 12   TRUE

The xdata attribute only contains the covariates:

##             X1          X2          X3
## 1  -0.07392768 -0.16681023  0.95278612
## 2  -0.44604373  0.76796254  0.22352026
## 3   0.12489742 -0.57420780 -0.99545858
## 4  -0.71203795  0.04409812  1.78181886
## 5   1.79128542 -0.92812647 -1.24928618
## 6   0.06870889  1.62915198 -0.76242481
## 7  -0.54569924  1.66808802 -1.72180296
## 8  -0.03863737 -0.44085091  0.45146624
## 9   0.90402590  0.46646237 -1.23922206
## 10  0.33471996 -0.77802688 -0.06973362

The baseline attribute contains the baseline functions – failure PDF, failure CDF, survivor, and hazard – generated by the flexible-hazard method:

##    time failure.PDF failure.CDF  survivor     hazard
## 1     1  0.01540280  0.01540280 0.9845972 0.01540280
## 2     2  0.01655095  0.03195375 0.9680463 0.01680987
## 3     3  0.01768918  0.04964293 0.9503571 0.01827307
## 4     4  0.01891866  0.06856159 0.9314384 0.01990690
## 5     5  0.02026468  0.08882627 0.9111737 0.02175633
## 6     6  0.02172725  0.11055352 0.8894465 0.02384534
## 7     7  0.02330636  0.13385988 0.8661401 0.02620322
## 8     8  0.02500201  0.15886189 0.8411381 0.02886601
## 9     9  0.02681420  0.18567609 0.8143239 0.03187848
## 10   10  0.02874294  0.21441903 0.7855810 0.03529669

As we demonstrate below, we can use the survsim.plot() function to plot these baseline functions.

The betas attribute contains the coefficients used in the simulation to generate the durations:

##            [,1]
## [1,] 0.06027158
## [2,] 0.12146851
## [3,] 0.05836502

These coefficients are the “true” coefficients in a data generating process, and a survival model such as the Cox proportional hazards model produces estimates of these coefficients. In this case, the Cox model estimates the coefficients as follows:

##          X1          X2          X3 
##  0.06740749  0.14893910 -0.02236168

The coefficients and the covariates together form values of the linear predictor, contained in the xb attribute, and exponentiated values of the linear predictor, contained in the exp.xb attribute:

##             [,1]      [,2]
## [1,]  0.03089145 1.0313735
## [2,]  0.07944527 1.0826863
## [3,] -0.12032036 0.8866363
## [4,]  0.06643678 1.0686934
## [5,] -0.07768915 0.9252520
## [6,]  0.15753291 1.1706193

In the course of generating simulations, the flexible-hazard method requires expressing each individual observation’s survivor function. These functions are contained in the rows of the ind.survive attribute, a matrix with N rows and T columns. For example, here’s the survivor function for observation 1:

##   [1] 9.841178e-01 9.670604e-01 9.488401e-01 9.293652e-01 9.085184e-01
##   [6] 8.861832e-01 8.622438e-01 8.365851e-01 8.090932e-01 7.796555e-01
##  [11] 7.481604e-01 7.121531e-01 6.714498e-01 6.293351e-01 5.890709e-01
##  [16] 5.538938e-01 5.270165e-01 5.116381e-01 5.034025e-01 4.958037e-01
##  [21] 4.888156e-01 4.824121e-01 4.765674e-01 4.712555e-01 4.664507e-01
##  [26] 4.621272e-01 4.582594e-01 4.548218e-01 4.517887e-01 4.491348e-01
##  [31] 4.468345e-01 4.448625e-01 4.431935e-01 4.418021e-01 4.406630e-01
##  [36] 4.397509e-01 4.390407e-01 4.385069e-01 4.381245e-01 4.378682e-01
##  [41] 4.377127e-01 4.376329e-01 4.376034e-01 4.375992e-01 4.375941e-01
##  [46] 4.375578e-01 4.374594e-01 4.372677e-01 4.369516e-01 4.364802e-01
##  [51] 4.358223e-01 4.349469e-01 4.338229e-01 4.324193e-01 4.307052e-01
##  [56] 4.286495e-01 4.262215e-01 4.233901e-01 4.201246e-01 4.163944e-01
##  [61] 4.095196e-01 3.974201e-01 3.809676e-01 3.610353e-01 3.384958e-01
##  [66] 3.142194e-01 2.890718e-01 2.639125e-01 2.395928e-01 2.169549e-01
##  [71] 1.968303e-01 1.800402e-01 1.673976e-01 1.597118e-01 1.557904e-01
##  [76] 1.535670e-01 1.525427e-01 1.522204e-01 1.521037e-01 1.516961e-01
##  [81] 1.505016e-01 1.480246e-01 1.437715e-01 1.300289e-01 1.045007e-01
##  [86] 7.480869e-02 4.842829e-02 3.251940e-02 2.484392e-02 1.849891e-02
##  [91] 1.335264e-02 9.274594e-03 6.135654e-03 3.808197e-03 2.166322e-03
##  [96] 1.086105e-03 4.459604e-04 1.271887e-04 1.489449e-05 0.000000e+00

The last two attributes, marg.effect and marg.effect.data, are discussed further below.

2.3 Simulating multiple datasets

In a simulation, researchers generally want to create many iterations of simulated data to draw inferences from the sampling variation that produces each dataset. To generate multiple datasets, set the num.data.frames argument to the desired number of iterations. For example, to create 2 datasets, type:

##      Length Class  Mode
## [1,] 9      -none- list
## [2,] 9      -none- list

The output is now a list of two objects, where each object contains another version of the attributes described above. For example, to the see the data for the second iteration, type

##            X1          X2         X3  y failed
## 1   0.7629271  0.03402121  1.1743157 46   TRUE
## 2  -0.4677598  0.80173462 -1.2249148 19   TRUE
## 3  -1.8897361  1.19052235 -0.4565519 46   TRUE
## 4  -1.4026512 -0.55020073 -0.1670356 46   TRUE
## 5  -1.1870153 -1.78614674 -2.1705478 47   TRUE
## 6   1.0341610 -0.76363408  1.7103576 17   TRUE
## 7  -0.2395232 -0.30959740 -0.3095540 45   TRUE
## 8   1.6391908 -0.03897204  0.3109840 46   TRUE
## 9  -0.7036793 -1.64006238  1.5113857 41  FALSE
## 10  0.7618128 -0.19042504 -1.5539722 46   TRUE

In the above example we created two datasets to reduce the time needed to compile this vignette, but for research applications this number can be set as large as a researcher wants, even to num.data.frames=1000 for example.

2.4 Plotting the baseline functions and histograms

The flexible-hazard method is designed to avoid assuming a particular parametric form for the hazard function, but it still employs a hazard function. To see the baseline functions, including hazard, use the survsim.plot() function. The first argument is the output of the sim.survdata() function. Every simulation iteration produces different durations (and if fixed.hazard=FALSE, different hazard functions), so it is necessary to use the df argument specify which simulation iteration to use when plotting histograms and baseline functions. To see the baseline functions for the first simulation iteration, specify type="baseline" as follows:

To see the histograms of the simulated durations, linear predictors, and exponentiated linear predictors, specify type="hist":

And to see both the baseline functions and the histograms, specify type="both":

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]

The survsim.plot() function also has a bins argument, which allows the user to change the number of bins in these histograms (the default is 30).

3 Changing simulation parameters

The sim.survdata() function has a number of arguments that allow the user to change parameters of the simulation.

3.1 Changing the number of observations, time points, and covariates, and the censoring rate

The N argument changes the number of observations in each simulated dataset, the T argument changes the maximum possible survival time, xvars sets the number of covariates, and the censor argument specifies the proportion of observations that should be designated as right-censored. By default, N=1000, T=100, xvars=3, and censor=.1. To generate data with 700 observations, with durations ranging from 1 to 250, with 5 covariates and 20% of observations right-censored, type:

##        X1                 X2                 X3           
##  Min.   :-3.32251   Min.   :-2.86431   Min.   :-3.285649  
##  1st Qu.:-0.67918   1st Qu.:-0.64491   1st Qu.:-0.671033  
##  Median : 0.01627   Median : 0.02412   Median :-0.004759  
##  Mean   : 0.00826   Mean   : 0.04375   Mean   : 0.005663  
##  3rd Qu.: 0.66174   3rd Qu.: 0.69897   3rd Qu.: 0.674779  
##  Max.   : 3.84799   Max.   : 2.74023   Max.   : 3.164353  
##        X4                 X5                  y         failed       
##  Min.   :-2.90819   Min.   :-2.834889   Min.   :  1   Mode :logical  
##  1st Qu.:-0.74946   1st Qu.:-0.696081   1st Qu.: 11   FALSE:153      
##  Median :-0.08525   Median :-0.012920   Median :103   TRUE :547      
##  Mean   :-0.08710   Mean   :-0.003722   Mean   : 86                  
##  3rd Qu.: 0.61666   3rd Qu.: 0.678233   3rd Qu.:113                  
##  Max.   : 3.21455   Max.   : 4.502816   Max.   :250

3.2 User-supplied coefficients

By default the coefficients are drawn from normal distributions, however the user may specify different coefficients with the beta argument if he or she wishes:

##      [,1]
## [1,]  1.0
## [2,] -1.5
## [3,]  2.0

Another way that a user can specify coefficients is to let these coefficients be dependent on time. An example of that is provided below.

3.3 User-supplied data

By default, the covariates are drawn from standard normal distributions, but the sim.survdata() function allows a user to specify different data for the covariates using the X argument. Suppose, for example, that we want to use the data on the length of time needed for governing coalitions in European democracies to conclude negotiations from Martin and Vanberg (2003). The data are stored in the coxed package as martinvanberg:

##     formdur           postel          prevdef            cont       
##  Min.   :  1.00   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:  6.00   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:0.0000  
##  Median : 19.00   Median :1.0000   Median :1.0000   Median :0.0000  
##  Mean   : 28.03   Mean   :0.5764   Mean   :0.8522   Mean   :0.3251  
##  3rd Qu.: 37.50   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :205.00   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##      ident           rgovm            pgovno         tpgovno      
##  Min.   :1.000   Min.   :0.0000   Min.   :1.000   Min.   : 0.000  
##  1st Qu.:1.000   1st Qu.:0.0000   1st Qu.:1.000   1st Qu.: 2.080  
##  Median :2.000   Median :0.4200   Median :2.000   Median : 5.890  
##  Mean   :2.049   Mean   :0.8009   Mean   :2.261   Mean   : 6.552  
##  3rd Qu.:3.000   3rd Qu.:1.4050   3rd Qu.:3.000   3rd Qu.: 8.600  
##  Max.   :3.000   Max.   :3.2700   Max.   :6.000   Max.   :25.470  
##     minority     
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.3596  
##  3rd Qu.:1.0000  
##  Max.   :1.0000

We will use the postel, rgovm, and pgovno variables to create a data frame that we pass to the X argument:

The data now contain the exogenous covariates and the simulated duration and censoring variables:

##   postel rgovm pgovno  y failed
## 1      1  3.01      2 14   TRUE
## 2      1  1.95      2 47   TRUE
## 3      1  2.36      2 85   TRUE
## 4      0  2.36      2 98   TRUE
## 5      0  2.36      2 95  FALSE
## 6      1  1.46      2 41   TRUE

3.4 Simulated marginal effects

The sim.survdata() function also generates a marginal change in duration, conditional on a particular covariate, so that researchers can compare the performance of estimators of this statistic using simulated data. The function calculates simulated durations for each observation conditional on a baseline hazard function and exogenous covariates and coefficients. The argument specifies the variable in the X matrix to vary so as to measure the marginal effect. First the covariate is set to the value specified in for all observations, then to the value specified in for all observations. Given each value, new durations are drawn. The durations when the covariate equals the low value are subtracted from the durations when the covariate equals the high value. The marginal effect is calculated by employing the statistic given by , which is by default. The marginal effect itself is stored in the marg.effect attribute, and the durations and covariates for the high and low profiles are stored as a list in the marg.effect.data attribute.

For example, suppose we are interested in comparing the duration when X1 is 1 to the duration when X1 is 0. We specify covariate=1 to indicate that X1 is the variable whose values should be fixed to the ones specified with low=0 and high=1. In this example we fix the coefficients so that X1 has a large effect, to ensure that we see a larger marginal effect. To simulate with these parameters, we type:

The simulated marginal effect is

## [1] -26

The covariates and durations for each of the low and high conditions are stored in marg.effect.data:

##   X1           X2          X3
## 1  0  1.345293789  2.08304792
## 2  0 -0.110570447 -1.35776100
## 3  0  0.723268017  0.45365921
## 4  0  0.862820298 -1.09684624
## 5  0 -1.009181418  0.10102341
## 6  0  0.004834229 -0.03386995
## [1] 16 43  8 88 92 48
##   X1           X2          X3
## 1  1  1.345293789  2.08304792
## 2  1 -0.110570447 -1.35776100
## 3  1  0.723268017  0.45365921
## 4  1  0.862820298 -1.09684624
## 5  1 -1.009181418  0.10102341
## 6  1  0.004834229 -0.03386995
## [1]  7 11  8 10 11  8

3.5 Fixed or varying hazard

A researcher might want to hold the hazard function fixed from one iteration to the next in order to isolate the effect of different factors. Alternatively, a researcher might want the hazard to vary from one iteration to the next in order to examine a range of possibilities within one analysis and to avoid making any restrictive assumptions about hazard. If fixed.hazard=TRUE, one baseline hazard is generated and the same function is used to generate all of the simulated datasets. If fixed.hazard=FALSE (the default), a new hazard function is generated with each simulation iteration. To illustrate, we create two simulated datasets, twice - once with fixed.hazard=TRUE and once with fixed.hazard=FALSE. If fixed.hazard=TRUE, both data frames yield the same baseline functions:

If fixed.hazard=FALSE, the two data frames yield different baseline functions:

3.6 Covariate-dependent censoring

An important assumption of many survival models is that the mechanism that causes some observations to be right-censored is independent from the covariates in the model. Some researchers, studying the effects of violations of the assumptions of a model, may want to dispel this assumption in the data generating process to see the effect on the model estimates. In the sim.survdata() function, if censor.cond=FALSE then a proportion of the observations specified by censor is randomly and uniformly selected to be right-censored. If censor.cond=TRUE then censoring depends on the covariates as follows: new coefficients are drawn from normal distributions with mean 0 and standard deviation of 0.1, and these new coefficients are used to create a new linear predictor using the X matrix. The observations with the largest (100 x censor) percent of the linear predictors are designated as right-censored. In other words, if censor.cond=FALSE, then a logit model in which the censoring indicator is regressed on the covariates should yield null results:

## 
## Call:
## glm(formula = failed ~ X1 + X2 + X3, family = binomial(link = "logit"), 
##     data = simdata$data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3010   0.4362   0.4617   0.4854   0.5730  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.15540    0.10401  20.722   <2e-16 ***
## X1          -0.07288    0.10159  -0.717    0.473    
## X2           0.10982    0.10354   1.061    0.289    
## X3          -0.06481    0.10196  -0.636    0.525    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 667.57  on 999  degrees of freedom
## Residual deviance: 665.62  on 996  degrees of freedom
## AIC: 673.62
## 
## Number of Fisher Scoring iterations: 5

If, however, censor.cond=TRUE, then a logit model in which the censoring indicator is regressed on the covariates should yield significant results:

## 
## Call:
## glm(formula = failed ~ X1 + X2 + X3, family = binomial(link = "logit"), 
##     data = simdata$data)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.80228   0.00820   0.05056   0.16435   2.63233  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   5.8652     0.4875  12.030  < 2e-16 ***
## X1            1.1229     0.1882   5.967 2.42e-09 ***
## X2           -0.4619     0.1847  -2.501   0.0124 *  
## X3            4.0505     0.3739  10.833  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 650.17  on 999  degrees of freedom
## Residual deviance: 245.80  on 996  degrees of freedom
## AIC: 253.8
## 
## Number of Fisher Scoring iterations: 8

3.7 User-specified hazard functions

For some applications, a researcher may want to specify his or her own hazard function rather than relying on the flexible-hazard method of generating hazard functions and drawing durations. The sim.survdata() allows a user to specify a customized hazard function with the hazard.fun argument, which takes a function with one argument that represents time. For example, suppose we want a lognormal hazard function with a mean of 50 and a standard deviation of 10. We can specify that function as follows:

We can then pass that function to sim.survdata(). Note that the lognormal hazard function does not equal 0 when T=100. That means that many observations will be naturally right-censored above and beyond the amount of right-censoring we’ve specified with the censor argument (a proportion of .1 by default), simply because we haven’t allowed the simulated durations to extend to later time points. To reduce the number of additional right-censored observations, we increase T to 1000 to allow for the longer durations that the log-normal hazard function implies. Note, the sim.survdata() produces a warning with the number of additional right-censored observations:

## Warning in FUN(X[[i]], ...): 151 additional observations right-censored because the user-supplied hazard function
##                                   is nonzero at the latest timepoint. To avoid these extra censored observations, increase T

To see the hazard, we can use the simsurv.plot() function with type="baseline":

3.8 Time-varying covariates

Time-varying covariates require a different data structure. The Surv() function in the survival package sets up the dependent variable for a Cox model. Generally it has two arguments: duration time and a censoring variable. But for time-varying covariates it replaces the duration argument with two time arguments, representing the start and end of discrete intervals, which allows a covariate to take on different values across different intervals for the same observation.

The sim.survdata() function generates data with a time-varying covariate if type="tvc". Durations are drawn again using proportional hazards, and are passed to the permalgorithm() function in the PermAlgo package to generate the time-varying data structure (Sylvestre and Abrahamowicz 2008). If type="tvc", the sim.survdata() function cannot accept user-supplied data for the covariates, as a time-varying covariate is expressed over time frames which themselves convey part of the variation of the durations, and we are generating these durations. If user-supplied X data is provided, the function passes a warning and generates random data instead as if . To generate survival data with 1000 observations, a maximum duration of 100, and time-varying covariates, type:

##    id failed start end          X1         X2         X3
## 1   1      0     0   1  0.80659639 -0.4126163 -0.2763447
## 2   1      0     1   2 -0.05678167 -0.4126163 -0.2763447
## 3   1      0     2   3 -0.64024574 -0.4126163 -0.2763447
## 4   1      0     3   4 -0.29441903 -0.4126163 -0.2763447
## 5   1      0     4   5  0.79861165 -0.4126163 -0.2763447
## 6   1      0     5   6 -1.71966224 -0.4126163 -0.2763447
## 7   1      0     6   7  0.77871246 -0.4126163 -0.2763447
## 8   1      0     7   8 -1.81961063 -0.4126163 -0.2763447
## 9   1      0     8   9 -1.60665020 -0.4126163 -0.2763447
## 10  1      0     9  10 -0.79584193 -0.4126163 -0.2763447
## 11  1      0    10  11 -0.61168950 -0.4126163 -0.2763447
## 12  1      0    11  12  1.04714193 -0.4126163 -0.2763447
## 13  1      0    12  13  1.59199971 -0.4126163 -0.2763447
## 14  1      0    13  14 -0.47626976 -0.4126163 -0.2763447
## 15  1      0    14  15  0.14671021 -0.4126163 -0.2763447
## 16  1      0    15  16  0.85775883 -0.4126163 -0.2763447
## 17  1      0    16  17  0.25678265 -0.4126163 -0.2763447
## 18  1      0    17  18  1.11156161 -0.4126163 -0.2763447
## 19  1      0    18  19  0.53144000 -0.4126163 -0.2763447
## 20  1      0    19  20 -0.87317660 -0.4126163 -0.2763447

3.9 Time-varying coefficients

The Cox model assumes proportional hazards, and one way to dispel this assumption is to allow the coefficients to vary over time. The sim.survdata() function generates data using time-varying coefficients if the type="tvbeta" argument is specified. If this option is specified, the first coefficient, whether coefficients are user-supplied or randomly generated, is interacted with the natural log of the time counter from 1 to T (the maximum possible duration). Durations are generated via proportional hazards, and coefficients are saved as a matrix to illustrate their dependence on time. To generate data with time-varying coefficients, type:

The coefficients are saved as a matrix with a column to represent time:

##    time        beta1       beta2      beta3
## 1     1  0.000000000 -0.03250039 0.09471359
## 2     2 -0.001096576 -0.03250039 0.09471359
## 3     3 -0.001738031 -0.03250039 0.09471359
## 4     4 -0.002193151 -0.03250039 0.09471359
## 5     5 -0.002546170 -0.03250039 0.09471359
## 6     6 -0.002834607 -0.03250039 0.09471359
## 7     7 -0.003078477 -0.03250039 0.09471359
## 8     8 -0.003289727 -0.03250039 0.09471359
## 9     9 -0.003476063 -0.03250039 0.09471359
## 10   10 -0.003642745 -0.03250039 0.09471359

An alternative approach to specifying time-varying coefficients is for the user to supply the matrix of time-varying coefficients. With this approach the coefficient matrix must have the same number of rows as the maximum duration specified with T, which in this case is 100. Suppose that we specify three coefficients for three covariates, over the time frame from 1 to 100. We want the first coefficient to be given by \[\beta_{1_t} = \frac{(t - 25)^2}{2500},\] and the second and third coefficients to be given by .5 and -.25 respectively:

##    beta1 beta2 beta3
## 1 0.2304   0.5 -0.25
## 2 0.2116   0.5 -0.25
## 3 0.1936   0.5 -0.25
## 4 0.1764   0.5 -0.25
## 5 0.1600   0.5 -0.25
## 6 0.1444   0.5 -0.25

We pass this matrix to the sim.survdata() function through the beta argument:

The data from this simulation are as follows:

##            X1          X2          X3  y failed
## 1  -0.1019964  0.07308832 -0.18423534 82   TRUE
## 2   0.4375385 -0.64651322  1.33692129 92  FALSE
## 3   1.1065059 -0.19794039 -0.63576681 22   TRUE
## 4   1.0991005  0.52697166  0.08448124 35   TRUE
## 5   1.4395873  0.66778150 -0.44694843 23   TRUE
## 6  -0.4807752 -0.46896444 -0.90295821 95   TRUE
## 7   0.7856846  0.97196353 -1.15458192 37   TRUE
## 8  -0.3084329  2.09888610 -1.51031320 26   TRUE
## 9  -1.5353525 -0.41373213 -0.65690715 99   TRUE
## 10 -0.2935970  0.43872654  0.67964543 27  FALSE

And the coefficients are as we specified earlier:

##    time  beta1 beta2 beta3
## 1     1 0.2304   0.5 -0.25
## 2     2 0.2116   0.5 -0.25
## 3     3 0.1936   0.5 -0.25
## 4     4 0.1764   0.5 -0.25
## 5     5 0.1600   0.5 -0.25
## 6     6 0.1444   0.5 -0.25
## 7     7 0.1296   0.5 -0.25
## 8     8 0.1156   0.5 -0.25
## 9     9 0.1024   0.5 -0.25
## 10   10 0.0900   0.5 -0.25

4 References