Hazard of death by time since randomization and time since recurrence
The two time variables are \(t\) = time since randomization and \(s\) = time since recurrence. Additionally, \(u\) = time at recurrence, is the fixed time at entry in the risk set recorded as number of days from randomization to recurrence. The analysis is performed over the \((u,s)\)-plane, so we need to bin the data over \(u\) and \(s\). We consider bins of width 30 days on both axes, and we obtain 77 bins on the \(u\) axis and \(91\) bins on the \(s\) axis.
dt2ts <- prepare_data(data = reccolon2ts,
u = "timer",
s_out = "timesr",
events = "status",
ds = 30, du = 30)> print(dt2ts)
An object of class 'data2ts'
Data:
List of 2
$ bins :List of 6
$ bindata:List of 2
- attr(*, "class")= chr "data2ts"
NULL
Range covered by the bins:
$bins_u
[1] 8 2288
$bins_s
[1] 0 2730
Number of bins:
$nu
[1] 76
$ns
[1] 91
Overview of the binned data:
Total exposure time: 246018
Total number of events: 409Here, we do not use the artificially created left truncated entry times on the \(s\)-axis, therefore the function returns a message informing the user that an entry time of 0 is imputed to all observations.
After having prepared the data we can estimate the model. In the analysis presented in Carollo et al. (2024), we build 23 \(B\)-splines over each of the two dimensions, for a total of 529 \(\alpha\) parameters. The optimal smoothing parameters are chosen by numerical optimization of the AIC of the model as function of the smoothing parameters. We use cubic \(B\)-splines bases and a second order penalty. Here we explicitly specify all these parameters, even though some of them are the default options.
mod1 <- fit2ts(data2ts = dt2ts,
Bbases_spec = list(bdeg = 3,
nseg_s = 20,
min_s = 0,
max_s = 2730,
nseg_u = 20,
min_u = 0,
max_u = 2300),
lrho = c(2, 0),
pord = 2,
optim_method = "ucminf",
optim_criterion = "aic") The object returned by fit2ts() is of class
'haz2ts'. The optimal smoothing parameters are \(\varrho_u = 10^{2.4}\) and \(\varrho_s=10^{0.3}\), and the effective
dimension of the model is 11.1.
> summary(mod1)
Number of events = 409
Model specifications:
nu = 76
ns = 91
cu = 23
cs = 23
Optimal smoothing:
log10(rho_u) = 2.402191
log10(rho_s) = 0.3104794
rho_u = 252.4593
rho_s = 2.043993
Model with no covariates
Model diagnostics:
AIC = 1249.297
BIC = 1314.651
ED = 11.13843We can change the optimization criterion to BIC, and compare the results in terms of smoothing parameters and effective dimensions:
mod2 <- fit2ts(data2ts = dt2ts,
Bbases_spec = list(bdeg = 3,
nseg_s = 20,
min_s = 0,
max_s = 2730,
nseg_u = 20,
min_u = 0,
max_u = 2300),
lrho = c(2, 0),
pord = 2,
optim_method = "ucminf",
optim_criterion = "bic")> mod2$optimal_logrho
[1] 6.137452 1.550218
> mod2$optimal_model$ed
[1] 5.451275As expected, choosing BIC as optimization criterion results in larger smoothing parameters and a smaller effective dimension, as BIC penalizes model complexity more strongly than AIC. The following code-chunk shows how to use the grid-search method to select the optimal pair of smoothing parameters and, at the same time, to produce plots of the AIC and BIC values of the grid of \(\log_{10}\) values of both smoothing parameters.
mod3 <- fit2ts(data2ts = dt2ts,
Bbases_spec = list(bdeg = 3,
nseg_s = 20,
min_s = 0,
max_s = 2730,
nseg_u = 20,
min_u = 0,
max_u = 2300),
optim_method = "grid_search",
optim_criterion = "aic",
lrho = list(seq(-1, 3, by = .2),
seq(-1, 3, by = .2)),
par_gridsearch = list(
plot_aic = TRUE,
plot_bic = TRUE,
mark_optimal = TRUE,
plot_contour = TRUE
))
Alternatively, we can ask the function to return the matrices of AIC and/or BIC values as part of the fitted object and then plot them separately (here not shown).
Lastly, we show how to fit the same model by using the package
LMMsolver Martin P. Boer (2023), that uses
the connection between linear mixed models and P-splines, and uses
sparse representation to speed-up calculations. It is possible to fit
this model by using the same function fit2ts() and
specifying the option optim_method = "LMMsolver". This
returns an object of class haz2tsLMM, which differ in
structure from objects of class haz2ts, but has the same
methods implemented.
mod_LMM <- fit2ts(data2ts = dt2ts,
Bbases_spec = list(bdeg = 3,
nseg_s = 20,
min_s = 0,
max_s = 2730,
nseg_u = 20,
min_u = 0,
max_u = 2300),
pord = 2,
optim_method = "LMMsolver",
optim_criterion = "aic")> summary(mod_LMM)
Number of events = 409
Model specifications:
nu = 76
ns = 91
cu = 23
cs = 23
Optimal smoothing:
log10(rho_u) = 2.14465
log10(rho_s) = 0.5217892
rho_u = 139.5243
rho_s = 3.324981
Model with no covariates
Model diagnostics:
AIC = 1249.449
BIC = 1314.151
ED = 11.02726Note: the functions from LMMsolver needed to fit the two time scales
hazard model are incorporated in the TwoTimeScales package.
Nevertheless, we do recommend interested readers to check out the
excellent R-package LMMsolver and its accompanying webpage:
https://biometris.github.io/LMMsolver/
