# 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.1364151 -0.5776570  0.69888490 58   TRUE
## 2   0.6696157 -1.2650216  0.59843683 24  FALSE
## 3  -0.3537196 -0.2893007  0.32638838 25   TRUE
## 4   1.1079940 -0.3852589 -0.89485350 63   TRUE
## 5  -0.5263717  0.5411580  0.01939824 23   TRUE
## 6   0.2871423 -0.3336064  0.30500820 26   TRUE
## 7  -0.4093448  0.8434817  0.36242347 24   TRUE
## 8  -0.3699496 -0.2268843 -0.12054188 33   TRUE
## 9   0.2419329  0.8151564 -0.22749104 43   TRUE
## 10 -0.6697677 -0.6267988 -0.05804856 24  FALSE

The xdata attribute only contains the covariates:

##            X1         X2          X3
## 1   0.1364151 -0.5776570  0.69888490
## 2   0.6696157 -1.2650216  0.59843683
## 3  -0.3537196 -0.2893007  0.32638838
## 4   1.1079940 -0.3852589 -0.89485350
## 5  -0.5263717  0.5411580  0.01939824
## 6   0.2871423 -0.3336064  0.30500820
## 7  -0.4093448  0.8434817  0.36242347
## 8  -0.3699496 -0.2268843 -0.12054188
## 9   0.2419329  0.8151564 -0.22749104
## 10 -0.6697677 -0.6267988 -0.05804856

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 5.479626e-06 5.479626e-06 0.9999945 5.479626e-06
## 2     2 3.835739e-05 4.383701e-05 0.9999562 3.835760e-05
## 3     3 1.041129e-04 1.479499e-04 0.9998521 1.041175e-04
## 4     4 2.027462e-04 3.506961e-04 0.9996493 2.027762e-04
## 5     5 3.342572e-04 6.849533e-04 0.9993150 3.343745e-04
## 6     6 4.986460e-04 1.183599e-03 0.9988164 4.989878e-04
## 7     7 6.959126e-04 1.879512e-03 0.9981205 6.967372e-04
## 8     8 9.260569e-04 2.805569e-03 0.9971944 9.278007e-04
## 9     9 1.189079e-03 3.994648e-03 0.9960054 1.192424e-03
## 10   10 1.484979e-03 5.479626e-03 0.9945204 1.490935e-03

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.232870415
## [2,] -0.124517232
## [3,]  0.002520564

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.34600951 -0.12134438  0.09152857

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.10545687 1.1112182
## [2,]  0.31495908 1.3702032
## [3,] -0.04552523 0.9554955
## [4,]  0.30373487 1.3549098
## [5,] -0.18991099 0.8270327
## [6,]  0.10917548 1.1153581

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] 1.000000e+00 9.999939e-01 9.999513e-01 9.998356e-01 9.996103e-01
##   [6] 9.992389e-01 9.986848e-01 9.979117e-01 9.968829e-01 9.955621e-01
##  [11] 9.939128e-01 9.918988e-01 9.894837e-01 9.866313e-01 9.833057e-01
##  [16] 9.794707e-01 9.750905e-01 9.701296e-01 9.645523e-01 9.583234e-01
##  [21] 9.514078e-01 9.437708e-01 9.353778e-01 8.613998e-01 7.275470e-01
##  [26] 6.364083e-01 6.141035e-01 5.927056e-01 5.648511e-01 5.341769e-01
##  [31] 5.042568e-01 4.785865e-01 4.574273e-01 4.386304e-01 4.220228e-01
##  [36] 4.074359e-01 3.947046e-01 3.836671e-01 3.741651e-01 3.660430e-01
##  [41] 3.591476e-01 3.533279e-01 3.484349e-01 3.443208e-01 3.408393e-01
##  [46] 3.378448e-01 3.351926e-01 3.327384e-01 3.303386e-01 3.278499e-01
##  [51] 3.251295e-01 3.220349e-01 3.184246e-01 3.141578e-01 3.090949e-01
##  [56] 3.030978e-01 2.960305e-01 2.877019e-01 2.785055e-01 2.691784e-01
##  [61] 2.604481e-01 2.530330e-01 2.476438e-01 2.449886e-01 2.440328e-01
##  [66] 2.432516e-01 2.426185e-01 2.421074e-01 2.416922e-01 2.413467e-01
##  [71] 2.410448e-01 2.407605e-01 2.404676e-01 2.401401e-01 2.397520e-01
##  [76] 2.392771e-01 2.386897e-01 2.379636e-01 2.370730e-01 2.359920e-01
##  [81] 2.346949e-01 2.331562e-01 2.313502e-01 2.292517e-01 2.268356e-01
##  [86] 1.952117e-01 1.476629e-01 1.153397e-01 8.832711e-02 6.608761e-02
##  [91] 4.809845e-02 3.385256e-02 2.285949e-02 1.464676e-02 8.761209e-03
##  [96] 4.770909e-03 2.267444e-03 8.690285e-04 2.249089e-04 2.230882e-05
## [101] 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  -1.09314491 -1.4504598028  0.1737110   5   TRUE
## 2   1.79619320 -0.3117859865  3.0575759  95   TRUE
## 3   0.17812434 -1.2098908917  1.4683695   5   TRUE
## 4   0.05424818  0.5908027294 -0.2074044  98   TRUE
## 5   0.91923027  0.0001077636 -0.7186776 100   TRUE
## 6   0.02050577  0.5912798242 -1.7256100   9   TRUE
## 7  -1.24611510 -0.7122088797  0.5182174   1   TRUE
## 8   0.55985820 -0.3225662877  0.7284522   4   TRUE
## 9  -1.26028511  0.3279342471  0.8787280  64   TRUE
## 10 -0.22411793  1.3395373036 -0.8451259  49   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.   :-1.54675   Min.   :-1.667013   Min.   :-1.80326
##  1st Qu.:-0.37239   1st Qu.:-0.335963   1st Qu.:-0.30612
##  Median :-0.01206   Median :-0.005102   Median : 0.02276
##  Mean   :-0.02194   Mean   :-0.009655   Mean   : 0.01370
##  3rd Qu.: 0.34929   3rd Qu.: 0.351629   3rd Qu.: 0.33478
##  Max.   : 1.44054   Max.   : 1.456765   Max.   : 1.78528
##        X4                  X5                  y           failed
##  Min.   :-1.444284   Min.   :-1.571971   Min.   :  1.0   Mode :logical
##  1st Qu.:-0.334935   1st Qu.:-0.308472   1st Qu.: 72.0   FALSE:151
##  Median : 0.009969   Median : 0.000606   Median : 83.0   TRUE :549
##  Mean   : 0.011683   Mean   : 0.021234   Mean   :100.6
##  3rd Qu.: 0.347150   3rd Qu.: 0.357980   3rd Qu.:111.0
##  Max.   : 1.694506   Max.   : 1.596067   Max.   :250.0

## 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:

## Warning in FUN(X[[i]], ...): 194 observations have drawn durations
##                                     at the minimum or maximum possible value. Generating coefficients
##                                     and other quantities of interest are unlikely to be returned
##                                     by models due to truncation. Consider making user-supplied coefficients
##                                     smaller, making T bigger, or decreasing the variance of the X variables.
##      [,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 60   TRUE
## 2      1  1.95      2 37   TRUE
## 3      1  2.36      2 59   TRUE
## 4      0  2.36      2 94   TRUE
## 5      0  2.36      2 53   TRUE
## 6      1  1.46      2 57   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:

## Warning in FUN(X[[i]], ...): 68 observations have drawn durations
##                                     at the minimum or maximum possible value. Generating coefficients
##                                     and other quantities of interest are unlikely to be returned
##                                     by models due to truncation. Consider making user-supplied coefficients
##                                     smaller, making T bigger, or decreasing the variance of the X variables.

The simulated marginal effect is

## [1] -41

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

##   X1         X2          X3
## 1  0 -0.0140624 -0.02211878
## 2  0  0.3877006  0.22235036
## 3  0  0.4844482 -0.27843449
## 4  0 -0.2344638  0.25951588
## 5  0  0.7505747  0.64858703
## 6  0  0.2811640  0.62776229
## [1] 36 45  9  2  3 53
##   X1         X2          X3
## 1  1 -0.0140624 -0.02211878
## 2  1  0.3877006  0.22235036
## 3  1  0.4844482 -0.27843449
## 4  1 -0.2344638  0.25951588
## 5  1  0.7505747  0.64858703
## 6  1  0.2811640  0.62776229
## [1]  1  2  1 11  1 10

## 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.3340 0.4001 0.4329 0.4632 0.5874 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.287101 0.110091 20.775 <2e-16 *** ## X1 0.003771 0.226858 0.017 0.9867 ## X2 0.197537 0.217842 0.907 0.3645 ## X3 -0.367358 0.219328 -1.675 0.0939 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 618.85 on 999 degrees of freedom ## Residual deviance: 615.24 on 996 degrees of freedom ## AIC: 623.24 ## ## 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.7233   0.2188   0.3335   0.4831   1.3369
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   2.6013     0.1423  18.274  < 2e-16 ***
## X1           -0.5812     0.2355  -2.468   0.0136 *
## X2            0.4152     0.2206   1.883   0.0597 .
## X3           -1.9690     0.2460  -8.004 1.21e-15 ***
## ---
## 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: 564.95  on 996  degrees of freedom
## AIC: 572.95
##
## Number of Fisher Scoring iterations: 6

## 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]], ...): 116 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 -1.3015315 0.4078455 0.3337465
## 2   1      0     1   2 -1.5982472 0.4078455 0.3337465
## 3   1      0     2   3 -1.5859503 0.4078455 0.3337465
## 4   1      0     3   4  1.4230914 0.4078455 0.3337465
## 5   1      0     4   5  2.0093171 0.4078455 0.3337465
## 6   1      0     5   6 -0.6421900 0.4078455 0.3337465
## 7   1      0     6   7  0.1502766 0.4078455 0.3337465
## 8   1      0     7   8  1.5255397 0.4078455 0.3337465
## 9   1      0     8   9  0.5023510 0.4078455 0.3337465
## 10  1      0     9  10 -1.2281639 0.4078455 0.3337465
## 11  1      0    10  11 -1.2572803 0.4078455 0.3337465
## 12  1      0    11  12 -0.1439728 0.4078455 0.3337465
## 13  1      0    12  13  0.2278704 0.4078455 0.3337465
## 14  1      0    13  14 -0.3037370 0.4078455 0.3337465
## 15  1      0    14  15  0.8369139 0.4078455 0.3337465
## 16  1      0    15  16 -0.1408818 0.4078455 0.3337465
## 17  1      0    16  17  1.0990125 0.4078455 0.3337465
## 18  1      0    17  18  2.2711663 0.4078455 0.3337465
## 19  1      0    18  19  0.3489595 0.4078455 0.3337465
## 20  1      1    19  20 -1.5322524 0.4078455 0.3337465

## 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.0000000 -0.1469655 0.1669951
## 2     2 0.1099974 -0.1469655 0.1669951
## 3     3 0.1743418 -0.1469655 0.1669951
## 4     4 0.2199948 -0.1469655 0.1669951
## 5     5 0.2554061 -0.1469655 0.1669951
## 6     6 0.2843392 -0.1469655 0.1669951
## 7     7 0.3088018 -0.1469655 0.1669951
## 8     8 0.3299922 -0.1469655 0.1669951
## 9     9 0.3486835 -0.1469655 0.1669951
## 10   10 0.3654035 -0.1469655 0.1669951

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:

## Warning in FUN(X[[i]], ...): 29 observations have drawn durations
##                                     at the minimum or maximum possible value. Generating coefficients
##                                     and other quantities of interest are unlikely to be returned
##                                     by models due to truncation. Consider making user-supplied coefficients
##                                     smaller, making T bigger, or decreasing the variance of the X variables.

The data from this simulation are as follows:

##              X1         X2           X3  y failed
## 1  -0.006828448  0.4729367 -0.195018482 60   TRUE
## 2  -0.335344072  0.1123093  0.234364749 42   TRUE
## 3   0.331864056 -0.9283091  0.344384499 39   TRUE
## 4  -0.437932466 -0.3667950 -0.440302465 16   TRUE
## 5  -0.004573757 -0.2354934  0.454745396 38   TRUE
## 6   0.934768177  1.0345235  0.002044735 45   TRUE
## 7  -0.382439326  0.8293747 -0.509321980 17   TRUE
## 8   0.711626280  0.4647250 -0.601334698 21   TRUE
## 9   0.280833582 -0.4734335  0.286422890 80   TRUE
## 10 -0.287489203  0.2305915  1.641656222 39   TRUE

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

• Harden, J. J. and Kropko, J. (2018). “Simulating Duration Data for the Cox Model.”" Political Science Research and Methods https://doi.org/10.1017/psrm.2018.19

• Hyman, J. M. (1983) “Accurate Monotonicity Preserving Cubic Interpolation.” SIAM J. Sci. Stat. Comput. 4, 645–654. https://doi.org/10.1137/0904045

• Martin, L. W and Vanberg, G. (2003) “Wasting Time? The Impact of Ideology and Size on Delay in Coalition Formation.” British Journal of Political Science 33 323-344 https://doi.org/10.1017/S0007123403000140

• Sylvestre M.-P., Abrahamowicz M. (2008) :“Comparison of Algorithms to Generate Event Times Conditional on Time-Dependent Covariates.” Statistics in Medicine 27(14):2618–34. https://doi.org/10.1002/sim.3092