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:
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).
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:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=1)
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.5240913 -1.6020903 0.2527589 96 TRUE
## 2 0.1239700 -0.9921704 1.7907474 88 TRUE
## 3 0.1639584 -0.5359740 -0.4428468 17 TRUE
## 4 -0.3686127 1.1937109 -1.4668718 17 TRUE
## 5 0.3049723 -0.5755462 -0.8669821 20 TRUE
## 6 -1.4262383 0.5881105 -0.2705132 14 TRUE
## 7 -1.2382806 -1.4294450 -1.2519687 84 TRUE
## 8 -0.7567318 0.2353116 0.4250513 14 TRUE
## 9 0.6222265 -0.1076480 -0.1095822 1 TRUE
## 10 0.4952521 -0.1975973 -0.4783299 17 TRUE
The xdata attribute only contains the covariates:
## X1 X2 X3
## 1 0.5240913 -1.6020903 0.2527589
## 2 0.1239700 -0.9921704 1.7907474
## 3 0.1639584 -0.5359740 -0.4428468
## 4 -0.3686127 1.1937109 -1.4668718
## 5 0.3049723 -0.5755462 -0.8669821
## 6 -1.4262383 0.5881105 -0.2705132
## 7 -1.2382806 -1.4294450 -1.2519687
## 8 -0.7567318 0.2353116 0.4250513
## 9 0.6222265 -0.1076480 -0.1095822
## 10 0.4952521 -0.1975973 -0.4783299
The baseline attribute contains the baseline functions – failure PDF, failure CDF, survivor, and hazard – generated by the flexible-hazard method:
head(simdata$baseline, 10)
## time failure.PDF failure.CDF survivor hazard
## 1 1 7.254861e-04 0.0007254861 0.9992745 7.254861e-04
## 2 2 6.141108e-04 0.0013395969 0.9986604 6.145567e-04
## 3 3 5.120168e-04 0.0018516138 0.9981484 5.127037e-04
## 4 4 4.192041e-04 0.0022708179 0.9977292 4.199818e-04
## 5 5 3.356727e-04 0.0026064906 0.9973935 3.364367e-04
## 6 6 2.614225e-04 0.0028679131 0.9971321 2.621057e-04
## 7 7 1.964536e-04 0.0030643667 0.9969356 1.970186e-04
## 8 8 1.407660e-04 0.0032051326 0.9967949 1.411986e-04
## 9 9 9.435960e-05 0.0032994922 0.9967005 9.466301e-05
## 10 10 5.723451e-05 0.0033567267 0.9966433 5.742398e-05
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.06176803
## [2,] -0.02686945
## [3,] -0.12833991
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:
model <- coxph(Surv(y, failed) ~ X1 + X2 + X3, data=simdata$data)
model$coefficients
## X1 X2 X3
## -0.08519733 -0.06808770 -0.08059583
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:
head(cbind(simdata$xb, simdata$exp.xb))
## [,1] [,2]
## [1,] -0.02176386 0.9784713
## [2,] -0.21082267 0.8099177
## [3,] 0.06110886 1.0630146
## [4,] 0.17895233 1.1959637
## [5,] 0.10789548 1.1139313
## [6,] 0.10701136 1.1129469
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] 0.99929013 0.99868922 0.99818821 0.99777802 0.99744955 0.99719374
## [7] 0.99700151 0.99686376 0.99677143 0.99671542 0.99668666 0.99667606
## [13] 0.99667455 0.98877923 0.93347248 0.79457398 0.59731217 0.38424921
## [19] 0.37310532 0.36312713 0.35425135 0.34641445 0.33955260 0.33360174
## [25] 0.32849759 0.32417566 0.32057127 0.31761958 0.31525561 0.31341428
## [31] 0.31203040 0.31103872 0.31037393 0.30997068 0.30976360 0.30968731
## [37] 0.30967641 0.30967372 0.30965489 0.30960377 0.30950424 0.30934014
## [43] 0.30909533 0.30875366 0.30836576 0.30799289 0.30763489 0.30729156
## [49] 0.30696273 0.30664823 0.30634787 0.30606148 0.30578888 0.30552989
## [55] 0.30528433 0.30505202 0.30483279 0.30462646 0.30443284 0.30425176
## [61] 0.30408304 0.30392651 0.30378198 0.30364927 0.30352820 0.30341861
## [67] 0.30332030 0.30323310 0.30315682 0.30309130 0.30303635 0.30299180
## [73] 0.30295745 0.30293315 0.30291869 0.30291391 0.30259871 0.30161663
## [79] 0.29991282 0.29743233 0.29412002 0.28992053 0.28477827 0.27863731
## [85] 0.27144133 0.26313357 0.25368628 0.24307126 0.23123002 0.21810302
## [91] 0.20362935 0.18774637 0.17038902 0.15148891 0.13097268 0.10875920
## [97] 0.08475388 0.05883515 0.03081070 0.00000000
The last two attributes, marg.effect and marg.effect.data, are discussed further below.
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:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=2)
summary(simdata)
## 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
head(simdata[[2]]$data, 10)
## X1 X2 X3 y failed
## 1 -0.6637338 1.4575280 0.11949721 60 TRUE
## 2 -0.4591399 -1.0422882 0.32166959 74 TRUE
## 3 -0.6664545 0.3410277 -0.94298409 6 TRUE
## 4 -0.2997064 0.2396492 -1.15308714 82 TRUE
## 5 -0.2285684 -1.3844324 0.42070486 71 FALSE
## 6 0.5835292 -0.2554817 -0.50814882 32 TRUE
## 7 -1.2534555 0.1194795 0.41693526 85 TRUE
## 8 0.1740976 -1.5212277 1.04338891 65 TRUE
## 9 0.8090130 -0.5800915 -0.35023788 69 TRUE
## 10 0.6006085 -1.6530343 0.07130339 47 FALSE
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.
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:
survsim.plot(simdata, df=1, type="baseline")
To see the histograms of the simulated durations, linear predictors, and exponentiated linear predictors, specify type="hist":
survsim.plot(simdata, df=1, type="hist")
And to see both the baseline functions and the histograms, specify type="both":
survsim.plot(simdata, df=1, 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).
Changing simulation parameters
The sim.survdata() function has a number of arguments that allow the user to change parameters of the simulation.
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:
simdata <- sim.survdata(N=700, T=250, xvars=5, censor=.2, num.data.frames = 1)
summary(simdata$data)
## X1 X2 X3
## Min. :-3.44724 Min. :-3.17257 Min. :-2.886606
## 1st Qu.:-0.63587 1st Qu.:-0.64283 1st Qu.:-0.707874
## Median : 0.07265 Median :-0.02332 Median : 0.006778
## Mean : 0.05814 Mean :-0.02137 Mean : 0.006236
## 3rd Qu.: 0.73802 3rd Qu.: 0.59765 3rd Qu.: 0.702254
## Max. : 2.68327 Max. : 3.05147 Max. : 3.590737
## X4 X5 y failed
## Min. :-2.777387 Min. :-3.327179 Min. : 1.0 Mode :logical
## 1st Qu.:-0.681695 1st Qu.:-0.638326 1st Qu.: 41.0 FALSE:152
## Median :-0.013661 Median :-0.002426 Median :191.5 TRUE :548
## Mean : 0.002193 Mean : 0.027107 Mean :153.5
## 3rd Qu.: 0.673026 3rd Qu.: 0.689157 3rd Qu.:247.0
## Max. : 3.246752 Max. : 3.169282 Max. :249.0
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:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=1, beta=c(1, -1.5, 2))
simdata$betas
## [,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.
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:
mv.data <- dplyr::select(martinvanberg, postel, rgovm, pgovno)
simdata <- sim.survdata(T=100, X=mv.data, num.data.frames = 1)
The data now contain the exogenous covariates and the simulated duration and censoring variables:
## postel rgovm pgovno y failed
## 1 1 3.01 2 16 TRUE
## 2 1 1.95 2 15 TRUE
## 3 1 2.36 2 60 TRUE
## 4 0 2.36 2 15 TRUE
## 5 0 2.36 2 53 FALSE
## 6 1 1.46 2 15 TRUE
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:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=1, covariate=1, low=0, high=1,
beta = c(2, .1, .1))
The simulated marginal effect is
## [1] -19
The covariates and durations for each of the low and high conditions are stored in marg.effect.data:
head(simdata$marg.effect.data$low$x)
## X1 X2 X3
## 1 0 0.3534795 0.04459034
## 2 0 -0.2432656 0.98062479
## 3 0 1.8719443 -0.64737419
## 4 0 -0.4859762 1.05314175
## 5 0 -0.1651501 1.81490247
## 6 0 0.1051242 -1.19036680
head(simdata$marg.effect.data$low$y)
## [1] 35 34 34 83 34 34
head(simdata$marg.effect.data$high$x)
## X1 X2 X3
## 1 1 0.3534795 0.04459034
## 2 1 -0.2432656 0.98062479
## 3 1 1.8719443 -0.64737419
## 4 1 -0.4859762 1.05314175
## 5 1 -0.1651501 1.81490247
## 6 1 0.1051242 -1.19036680
head(simdata$marg.effect.data$high$y)
## [1] 15 14 4 20 33 34
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:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=2, fixed.hazard=TRUE)
survsim.plot(simdata, df=1, type="baseline")
survsim.plot(simdata, df=2, type="baseline")
If fixed.hazard=FALSE, the two data frames yield different baseline functions:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=2, fixed.hazard=FALSE)
survsim.plot(simdata, df=1, type="baseline")
survsim.plot(simdata, df=2, type="baseline")
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:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=1, censor.cond=FALSE)
logit <- glm(failed ~ X1 + X2 + X3, data=simdata$data, family=binomial(link="logit"))
summary(logit)
##
## Call:
## glm(formula = failed ~ X1 + X2 + X3, family = binomial(link = "logit"),
## data = simdata$data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3333 0.4100 0.4424 0.4694 0.5746
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.25358 0.10847 20.776 <2e-16 ***
## X1 -0.03975 0.10706 -0.371 0.710
## X2 -0.11717 0.11337 -1.033 0.301
## X3 -0.13811 0.10667 -1.295 0.195
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 632.41 on 999 degrees of freedom
## Residual deviance: 629.60 on 996 degrees of freedom
## AIC: 637.6
##
## 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:
simdata <- sim.survdata(N=1000, T=100, num.data.frames=1, censor.cond=TRUE)
logit <- glm(failed ~ X1 + X2 + X3, data=simdata$data, family=binomial(link="logit"))
summary(logit)
##
## Call:
## glm(formula = failed ~ X1 + X2 + X3, family = binomial(link = "logit"),
## data = simdata$data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.62651 0.01108 0.05474 0.18931 2.89615
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.6137 0.4723 11.885 < 2e-16 ***
## X1 3.2675 0.3236 10.098 < 2e-16 ***
## X2 1.6181 0.2091 7.739 1.00e-14 ***
## X3 1.5614 0.2162 7.221 5.15e-13 ***
## ---
## 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: 277.27 on 996 degrees of freedom
## AIC: 285.27
##
## Number of Fisher Scoring iterations: 8
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:
my.hazard <- function(t){
dnorm((log(t) - log(50))/log(10)) /
(log(10)*t*(1 - pnorm((log(t) - log(50))/log(10))))
}
We can then pass that function to sim.survdata():
simdata <- sim.survdata(N=1000, T=100, num.data.frames=1, hazard.fun = my.hazard)
To see the hazard, we can use the simsurv.plot() function with type="baseline":
survsim.plot(simdata, df=1, type="baseline")
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:
simdata <- sim.survdata(N=1000, T=100, type="tvc", num.data.frames=1)
head(simdata$data, 20)
## id failed start end X1 X2 X3
## 1 1 0 0 1 -1.0309368 1.312018 -0.2360551
## 2 1 0 1 2 0.1688022 1.312018 -0.2360551
## 3 1 0 2 3 -1.3901208 1.312018 -0.2360551
## 4 1 0 3 4 -0.5943993 1.312018 -0.2360551
## 5 1 0 4 5 -1.3798346 1.312018 -0.2360551
## 6 1 0 5 6 0.1614912 1.312018 -0.2360551
## 7 1 0 6 7 1.1600534 1.312018 -0.2360551
## 8 1 0 7 8 0.0714863 1.312018 -0.2360551
## 9 1 0 8 9 0.2376392 1.312018 -0.2360551
## 10 1 0 9 10 -1.1128862 1.312018 -0.2360551
## 11 1 0 10 11 0.2292287 1.312018 -0.2360551
## 12 1 0 11 12 1.1123070 1.312018 -0.2360551
## 13 1 0 12 13 0.5223005 1.312018 -0.2360551
## 14 1 0 13 14 0.5592167 1.312018 -0.2360551
## 15 1 0 14 15 0.6726185 1.312018 -0.2360551
## 16 1 0 15 16 -0.3910532 1.312018 -0.2360551
## 17 1 0 16 17 -1.5024762 1.312018 -0.2360551
## 18 1 0 17 18 0.2488101 1.312018 -0.2360551
## 19 1 0 18 19 0.2354240 1.312018 -0.2360551
## 20 1 0 19 20 -0.6097164 1.312018 -0.2360551
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:
simdata <- sim.survdata(N=1000, T=100, type="tvbeta", num.data.frames = 1)
The coefficients are saved as a matrix with a column to represent time:
## time beta1 beta2 beta3
## 1 1 0.000000000 0.08060892 0.1203639
## 2 2 0.008545802 0.08060892 0.1203639
## 3 3 0.013544776 0.08060892 0.1203639
## 4 4 0.017091605 0.08060892 0.1203639
## 5 5 0.019842739 0.08060892 0.1203639
## 6 6 0.022090579 0.08060892 0.1203639
## 7 7 0.023991101 0.08060892 0.1203639
## 8 8 0.025637407 0.08060892 0.1203639
## 9 9 0.027089553 0.08060892 0.1203639
## 10 10 0.028388541 0.08060892 0.1203639
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:
beta.mat <- data.frame(beta1 = (1:100 - 25)^2/2500,
beta2 = .5,
beta3 = -.25)
head(beta.mat)
## 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:
simdata <- sim.survdata(N=1000, T=100, type="tvbeta", beta=beta.mat, num.data.frames = 1)
The data from this simulation are as follows:
## X1 X2 X3 y failed
## 1 -0.009923238 0.3301969 -0.2191167 4 TRUE
## 2 1.992899129 -0.1943596 0.6749138 62 TRUE
## 3 -1.001732940 -1.0630761 -1.5919218 99 TRUE
## 4 -1.595622186 -1.1243103 -0.2534580 5 TRUE
## 5 -0.245100625 -1.1127584 3.1503928 32 TRUE
## 6 -1.722751180 1.3392698 0.8861647 10 TRUE
## 7 -0.358534629 -0.8264453 0.1461627 33 TRUE
## 8 1.748877594 -0.8875988 1.1383078 55 FALSE
## 9 1.452307649 0.6467879 -0.2992108 1 FALSE
## 10 -0.003008743 0.7989141 -0.1984989 3 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