joint_models_report_part3
Abderrahim Oussama Batouche
2024-02-15
Background
This document is a continuation of the previous parts. The treatment free survival in this document is defined as the time to the first curative treatment (RP/RT). We will model the :
Effect of Time to curative treatment
(rp+rt)on Time to relapse from curative treatment.In this case 1 in the
tfs_statusmeans either rp or rt occurred.and
tfs_status= 0 means non of them occurred
Effect of Time to
specificcurative treatment(rp or rt)on Time to relapse fromanycurative treatment.These are 2 different models, in the first one status is = 1 if treatment is rp and 0 otherwise, and the second model status is =1 if treatment is rt and 0 otherwise.
So this data should contain only RP and RT.
Effect of Time to
specificcurative treatment(rp or rt)on Time to relapse fromspecificcurative treatment.- These are also 2 models, as the previous ones, but in these models we’re assessing the effect of time to specific type of curative treatment(let’s say RT) on the risk of recurrence from that type of treatment (RT).
It is important to keep in mind that in the joint survival model,
patients are on risk of relapse only if they get a curative treatment ,
this is why we use left truncation to take into account the fact the the
recurrence risk starts at treatment time
(truncation = tfs_time). Also, as the
INLAjoint survival function tries to match the TFS data
with RFS data, it’s crucial to keep the RFS patients with no treatment
(tfs_status = 0 and rfs_status = 0) at the end
of the data. These patients are important for the first model (represent
the censored ones) but they aren’t important for the second model (no
treatment = no risk for bcr). Moreover, in order for the
INLAjoint survival function to work properly, time to event
must be greater then 0
(tfs_time >0andrfs_time>0).
We will try simpler models for PSA data, the first one is a simple
linear model (y=ax+b, y:psa ; x=time_to_psa) . this will
give us a simple linear model, that only captures the trend.
Next, we will try splines, 2 and 3 splines to capture to PSA
variation and give the model more freedom. Later we can also combine
quadratic and cubic functions (f1:y=x*x and f2:y=x*x*x). It
would be interesting to try the combination between linear + cubic.
Keep in mind that reference in factors is important.
The interpretation is hold only for Surv models and the
surv.surv joint model. The results from JM-1 are the only
needed ones to assess the effect of time to curative treatment on risk
of relapse (time to relapse).
join model including surv.surv.long, will be used to
explain the JM-1 results. While the JM-1 is showing
negative effect of time to CTx on recurrence
(lowest risk of treatment = highest risk of recurrence -accounting for cov-),
the JM-2 that includes the shared effect of the
longitudinal PSA data, showed a positive effect. This could be due to
the role of PSA, as we are using it as confounder (not on the causal
pathway), while it could be seen as a mediator or a moderator. In other
words, while accounting for PSA, time to treatment has positive effect
on risk of recurrence (increasing PSA
(slop>0) => higher risk of treatment = higher risk of recurrence
).
In JM-1 we’re assuming CTx as RP or RT for both treatment survival and recurrence survival. There are other variations that we can test, i.e, keep only one of the treatments, and modeling the effect of time to CTxi on risk of recurrence.
In JM-2 CTx we split the cohort into RP or RT, keep it in mind while interpreting. We can use both later if needed for interpretation (and for explaining the effect of PSA).
Libraries
library(readr)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(INLA)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: sp## This is INLA_24.02.09 built 2024-02-09 03:43:24 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - List available models/likelihoods/etc with inla.list.models()
## - Use inla.doc(<NAME>) to access documentation#inla.update(testing=TRUE)
library(INLAjoint)## Package 'INLAjoint' version 24.2.5
## Type 'citation("INLAjoint")' for citing this R package in publications.library(splines)Importing data
surv_long_ <- read.csv("processed/surv_long_.csv")
# We can remove the following after getting the new treatment data
# The new treatment data is supposed to contain only RP or RT when status=1 and UNK when status=0
# Now we filter the available data
survlong_rprt <- surv_long_[surv_long_$tx1_type %in% c('rp','radio','tfs-unk'),]
# Rename attributes
names(survlong_rprt)[names(survlong_rprt) == "tx1_type"] <- "tfs_tx"
names(survlong_rprt)[names(survlong_rprt) == "after_tx"] <- "drfs_tx"
# Change values
survlong_rprt$tfs_tx[survlong_rprt$tfs_tx %in% c('radio')] <- 'rt'
survlong_rprt$tfs_tx[survlong_rprt$tfs_tx %in% c('tfs-unk')] <- 'unk'
survlong_rprt$drfs_tx[survlong_rprt$drfs_tx %in% c('bcr-unk')] <- 'unk'
# There shouldn't be any time to event = 0 (checked status 0, no time 0, good!) ,
# drfs_time min = 1 , good!
survlong_rprt$tfs_time[(survlong_rprt$tfs_status == 1) & (survlong_rprt$tfs_time == 0)] <- 0.5
# Convert time to event from Months to Years
survlong_rprt$tfs_time_m <- survlong_rprt$tfs_time
survlong_rprt$drfs_time_m <- survlong_rprt$drfs_time
survlong_rprt$time_to_psa_m <- survlong_rprt$time_to_psa
survlong_rprt$tfs_time <- round(survlong_rprt$tfs_time_m/12, 2)
survlong_rprt$drfs_time <- round(survlong_rprt$drfs_time_m/12, 2)
survlong_rprt$time_to_psa <- round(survlong_rprt$time_to_psa_m/12, 2)Functions
reassigning_ids <- function(data){
id_rle <- rle(data$id)
new_id <- rep(seq_along(id_rle$lengths), id_rle$lengths)
data$id <- new_id
return(data)
}
bf_sampling <- function(data, sample_size=0.3){
#data <- survlong_rprt
#sample_size <- 1
# Uncomment when needed, Sampling
unique_ids <- unique(data$id)
sample_size <- round(sample_size * length(unique_ids))
# Sample X.% of the unique IDs
sampled_ids <- sample(unique_ids, sample_size, replace = FALSE)
# Subset the data based on sampled IDs
sampled_data <- data[data$id %in% sampled_ids, ]
# Factors
sampled_data$mri <- as.factor(sampled_data$mri)
sampled_data$grade_group <- as.factor(sampled_data$grade_group)
sampled_data$bx_age <- as.factor(sampled_data$bx_age)
sampled_data$tfs_status <- as.factor(sampled_data$tfs_status)
sampled_data$drfs_status <- as.factor(sampled_data$drfs_status)
sampled_data$drfs_tx <- as.factor(sampled_data$drfs_tx)
sampled_data$tfs_tx <- as.factor(sampled_data$tfs_tx)
# Order factors
sampled_data$tfs_tx <- factor(sampled_data$tfs_tx, levels = c('rt', 'rp', 'unk'))
sampled_data$drfs_tx <- factor(sampled_data$drfs_tx, levels = c('rt', 'rp', 'unk'))
# Set reference
sampled_data$tfs_tx <- relevel(sampled_data$tfs_tx, ref = "rt")
sampled_data$drfs_tx <- relevel(sampled_data$drfs_tx, ref = "rt")
# Remove Nas
sampled_data <- sampled_data[!is.na(sampled_data$time_to_psa),]
# Reassigning ID
# We have to do it here, as all will be reordered after the previous cmd
# -> Calculate run-length encoding of the original id column
id_rle <- rle(sampled_data$id)
# -> Generate a sequence of IDs based on the lengths of runs in the original id column
new_id <- rep(seq_along(id_rle$lengths), id_rle$lengths)
# -> Reassign the new IDs to the id column in sampled_data
sampled_data$id <- new_id
sampled_data <- sampled_data[with(sampled_data, order(id)), ]
# Order it by :
# 1- Type of treatment at recurrence drfs_tx, so the unk ones will be at the end.
# 2- Time to psa , so the INLAjoint Long model works properly
# We already have PSA at t=0 for all these patients from part 1
#sampled_data <- sampled_data[with(sampled_data, order(drfs_tx, time_to_psa)), ]
sampled_data <- sampled_data[order(sampled_data$drfs_tx,
sampled_data$id,
sampled_data$time_to_psa), ]
# Quick reordering
id_rle <- rle(sampled_data$id)
new_id <- rep(seq_along(id_rle$lengths), id_rle$lengths)
sampled_data$id <- new_id
return(sampled_data)
}
interp <- function(jm){
n_sample <- 1e4 # number of samples for uncertainty quantification
smp_H <- inla.hyperpar.sample(n_sample, jm) # sample values for frailty and association parameters
sigma <- sqrt(1/smp_H[, which(colnames(smp_H)=="Precision for IDIntercept_S1")]) # sampled sd of frailty
assoc <- smp_H[, which(colnames(smp_H)=="Beta for IDIntercept_S1_S2")] # sampled association parameter
SMP <- sapply(sigma, function(x) rnorm(1e3, mean = 0, sd = x)) # sample realizations of each frailty
SMP2 <- rbind(sigma, SMP) # add sigmas there to vectorize computations and avoid loop
mean_low15 <- apply(SMP, 2, function(x) mean(x[-1][x[-1]<(-x[1])])) # mean frailty deviation for lowest 15% , for 2.5% -2x[1]
mean_up15 <- apply(SMP, 2, function(x) mean(x[-1][x[-1]>x[1]])) # mean frailty deviation for top 15%, same here
HR_low15 <- exp(assoc*mean_low15) # hazard ratios for lowest 15% vs. average
HR_up15 <- exp(assoc*mean_up15) # hazard ratios for top 15% vs. average
Q_low15 <- quantile(HR_low15, c(0.025, 0.5, 0.975), na.rm=T) # should not have NAs, here it's a bad example on unstable model
Q_up15 <- quantile(HR_up15, c(0.025, 0.5, 0.975), na.rm=T)
# Interpretation:
res_l15 <- paste0("Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the lowest risk of receiving treatment (i.e., longest time-to-treatment) is associated to a ", round(Q_low15[2], 2), " [", round(Q_low15[1], 2), ",", round(Q_low15[3], 2), "] increased risk of recurrence (i.e., shorter time-to-recurrence) compared to the average individual.")
res_u15 <- paste0("Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the highest risk of receiving treatment (i.e., shortest time-to-treatment) is associated to a ", round(Q_up15[2], 2), " [", round(Q_up15[1], 2), ",", round(Q_up15[3], 2), "] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual.")
return(c(res_l15, res_u15))
}Modeling
Sampling
# Sampling
s_survlong_rprt <- bf_sampling(survlong_rprt, 1)
# Surv (only) data to be used in the latter join model
s_surv_rprt <- s_survlong_rprt[!duplicated(s_survlong_rprt$id), ]
summary(s_survlong_rprt)## id bx_age mri grade_group tfs_tx
## Min. : 1 <60 :29259 0:135874 1:58040 rt :73137
## 1st Qu.:2106 >70 :49074 1: 16264 2:42686 rp :46638
## Median :4200 60-70:73805 3:28561 unk:32363
## Mean :4289 4:12730
## 3rd Qu.:6399 5:10121
## Max. :9030
## tx_date tfs_time tfs_status drfs_tx drfs_time
## Length:152138 Min. : 0.040 0: 32363 rt :73137 Min. : 0.080
## Class :character 1st Qu.: 0.250 1:119775 rp :46638 1st Qu.: 3.580
## Mode :character Median : 0.500 unk:32363 Median : 6.500
## Mean : 2.647 Mean : 7.081
## 3rd Qu.: 3.330 3rd Qu.:10.000
## Max. :21.580 Max. :23.670
## drfs_status psa_results psa_date time_to_psa
## 0:106475 Min. : 0.00 Length:152138 Min. :-20.330
## 1: 45663 1st Qu.: 0.05 Class :character 1st Qu.: 0.250
## Median : 1.06 Mode :character Median : 2.170
## Mean : 11.57 Mean : 2.741
## 3rd Qu.: 5.69 3rd Qu.: 5.170
## Max. :76800.00 Max. : 23.670
## log_psa tfs_time_m drfs_time_m time_to_psa_m
## Min. : 0.00000 Min. : 0.50 Min. : 1.00 Min. :-244.00
## 1st Qu.: 0.04879 1st Qu.: 3.00 1st Qu.: 43.00 1st Qu.: 3.00
## Median : 0.72271 Median : 6.00 Median : 78.00 Median : 26.00
## Mean : 1.07824 Mean : 31.77 Mean : 84.97 Mean : 32.89
## 3rd Qu.: 1.90061 3rd Qu.: 40.00 3rd Qu.:120.00 3rd Qu.: 62.00
## Max. :11.24897 Max. :259.00 Max. :284.00 Max. : 284.00M-1: Time to Treatment (SURV)
# No Covariates
M1 <- joint(formSurv =list(inla.surv(time = tfs_time, event = tfs_status) ~ -1),
basRisk = c("rw2"),
dataSurv = list(s_surv_rprt))
summary(M1)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.9567 0.2999 0.5055 0.9103 1.6771
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -36491.23 -36491.23
##
## Deviance Information Criterion: 72893.69
## Widely applicable Bayesian information criterion: 72898.28
## Computation time: 2.06 secondsplot(M1)## $Baseline
##
## attr(,"class")
## [1] "plot.INLAjoint" "list"# With Covariates
M1_cov <- joint(formSurv =list(inla.surv(time = tfs_time, event = tfs_status) ~ bx_age+mri+grade_group),
basRisk = c("rw2"),
dataSurv = list(s_surv_rprt))
summary(M1_cov)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.8920 0.2821 0.4694 0.8477 1.5712
## bx_age70 0.1048 0.0230 0.0598 0.1048 0.1499
## bx_age6070 0.1326 0.0213 0.0908 0.1326 0.1744
## mri1 0.2587 0.0242 0.2111 0.2587 0.3062
## grade_group2 0.5619 0.0202 0.5223 0.5619 0.6014
## grade_group3 0.8328 0.0226 0.7884 0.8328 0.8771
## grade_group4 0.7681 0.0305 0.7082 0.7681 0.8279
## grade_group5 0.8577 0.0320 0.7949 0.8577 0.9204
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -35447.77 -35447.77
##
## Deviance Information Criterion: 70736.79
## Widely applicable Bayesian information criterion: 70752.54
## Computation time: 2.58 secondsplot(M1_cov)## $Outcomes
## $Outcomes$S1
##
##
## $Baseline
##
## attr(,"class")
## [1] "plot.INLAjoint" "list"M-2: Time to Recurrence (SURV)
# We need separate df for BCR (only treated patients)
s_surv_rprt_bcr <- s_surv_rprt[s_surv_rprt$tfs_status==1,]
s_surv_rprt_bcr$drfs_tx <- droplevels(s_surv_rprt_bcr$drfs_tx)
s_surv_rprt_bcr$tfs_tx <- droplevels(s_surv_rprt_bcr$tfs_tx)
# No Covariates
M2 <- joint(formSurv =list(inla.surv(time = drfs_time, event = drfs_status, truncation = tfs_time) ~ -1),
basRisk = c("rw2"),
dataSurv = s_surv_rprt_bcr) # Make sure only treated patients are included in this model
summary(M2)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.0741 0.0564 0.0142 0.0587 0.2269
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -27668.26 -27668.26
##
## Deviance Information Criterion: 55283.86
## Widely applicable Bayesian information criterion: 55287
## Computation time: 1.97 secondsplot(M2)## $Baseline
##
## attr(,"class")
## [1] "plot.INLAjoint" "list"# With Covariates
M2_cov <- joint(formSurv =list(inla.surv(time = drfs_time, event = drfs_status, truncation = tfs_time) ~ bx_age+mri+grade_group+drfs_tx),
basRisk = c("rw2"),
dataSurv = s_surv_rprt_bcr) # Make sure only treated patients are included in this model
summary(M2_cov)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.0653 0.0512 0.0121 0.0510 0.2048
## bx_age70 0.2406 0.0319 0.1780 0.2406 0.3031
## bx_age6070 0.0389 0.0286 -0.0171 0.0389 0.0950
## mri1 0.8323 0.0359 0.7619 0.8323 0.9027
## grade_group2 0.0880 0.0282 0.0327 0.0880 0.1432
## grade_group3 0.3035 0.0305 0.2437 0.3035 0.3634
## grade_group4 0.4739 0.0388 0.3978 0.4739 0.5501
## grade_group5 0.9123 0.0402 0.8334 0.9123 0.9912
## drfs_txrp 0.3229 0.0242 0.2756 0.3229 0.3703
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -27044.39 -27044.39
##
## Deviance Information Criterion: 53959.53
## Widely applicable Bayesian information criterion: 53966.98
## Computation time: 2.56 secondsplot(M2_cov)## $Outcomes
## $Outcomes$S1
##
##
## $Baseline
##
## attr(,"class")
## [1] "plot.INLAjoint" "list"JM-1: Treatment.X.Recurrence (SURV.SURV)
# Because of:
# 1 - the truncation time must be only related to RP/RT , and
# 2 - the treatments in the first surv model are not all RP/RT
# We have to make sure that those patients with first treatment as NOT a curative treatment (RP/RT)
# should be at the end of the dataframe. (to avoid mismatch)
#surv <- surv[with(surv, order(tx1_type)), ]
# No covariates
JM1 <- joint(formSurv=list(inla.surv(time = tfs_time, event = tfs_status) ~ -1+(1|id),
inla.surv(time = drfs_time, event = drfs_status, truncation=tfs_time) ~ -1),
id="id",
basRisk=c("rw2", "rw2"), assocSurv=TRUE, NbasRisk = 15,
dataSurv = list(s_surv_rprt,
s_surv_rprt_bcr# Make sure only treated patients are included in this model
)
)
summary(JM1)##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.3696 0.028 0.3235 0.3663 0.4327
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 2.8158 0.0972 2.6311 2.8137 3.013
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0512 0.0036 0.045 0.0509 0.059
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 0.0215 0.0142 -0.0037 0.0206 0.0519
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -59399.57 -59399.57
##
## Deviance Information Criterion: 114891.1
## Widely applicable Bayesian information criterion: 123807.8
## Computation time: 11.03 secondsplot(JM1)$Baseline+scale_y_log10()
# With Covariates
JM1_cov <- joint(formSurv=list(inla.surv(time = tfs_time , event = tfs_status) ~ bx_age+mri+grade_group+(1|id),
inla.surv(time = drfs_time, event = drfs_status, truncation=tfs_time) ~ bx_age+mri+grade_group+drfs_tx), id="id",
basRisk=c("rw2", "rw2"), assocSurv=TRUE, NbasRisk = 15,
dataSurv = list(s_surv_rprt,
s_surv_rprt_bcr # Make sure only treated patients are included in this model
)
)
summary(JM1_cov)##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.2238 0.0333 0.1691 0.2199 0.2995
## bx_age70_S1 0.1817 0.0517 0.0804 0.1817 0.2832
## bx_age6070_S1 0.2021 0.0483 0.1074 0.2021 0.2970
## mri1_S1 0.3768 0.0543 0.2704 0.3768 0.4834
## grade_group2_S1 1.0156 0.0460 0.9255 1.0156 1.1060
## grade_group3_S1 1.3612 0.0524 1.2585 1.3611 1.4641
## grade_group4_S1 1.3197 0.0706 1.1814 1.3196 1.4582
## grade_group5_S1 1.4324 0.0731 1.2893 1.4324 1.5760
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 2.3502 0.1441 2.09 2.3418 2.6559
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0479 0.0074 0.0337 0.0479 0.0624
## bx_age70_S2 0.2390 0.0326 0.1751 0.2390 0.3029
## bx_age6070_S2 0.0288 0.0293 -0.0286 0.0288 0.0862
## mri1_S2 0.8063 0.0366 0.7346 0.8063 0.8780
## grade_group2_S2 0.0379 0.0289 -0.0188 0.0379 0.0946
## grade_group3_S2 0.2337 0.0315 0.1720 0.2337 0.2953
## grade_group4_S2 0.4230 0.0400 0.3446 0.4230 0.5013
## grade_group5_S2 0.8498 0.0415 0.7683 0.8498 0.9312
## drfs_txrp_S2 0.3444 0.0244 0.2966 0.3444 0.3921
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 -0.1293 0.021 -0.1705 -0.1293 -0.0877
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -58246.52 -58246.52
##
## Deviance Information Criterion: 110012.2
## Widely applicable Bayesian information criterion: 112891.6
## Computation time: 11.9 secondsplot(JM1_cov)$Baseline+scale_y_log10()
int_res <- interp(JM1_cov)
int_res[1]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the lowest risk of receiving treatment (i.e., longest time-to-treatment) is associated to a 1.16 [1.01,1.6] increased risk of recurrence (i.e., shorter time-to-recurrence) compared to the average individual."int_res[2]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the highest risk of receiving treatment (i.e., shortest time-to-treatment) is associated to a 0.86 [0.63,0.99] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."# Survival curves
onePatient <- s_surv_rprt[1, ]
P <- predict(M1, onePatient, id="id", horizon=30, surv=TRUE)$PredS## Warning in predict.INLAjoint(M1, onePatient, id = "id", horizon = 30, surv =
## TRUE): The fitted model has baseline risk information up until value 21.58 for
## survival outcome 1. Since you ask for prediction at horizon 30 I will assume
## constant baseline hazard beyond the maximum available value. Alternatively, you
## can use baselineHaz='smooth' to use splines to predict the baseline hazard (for
## each sample). Alternatively, adding 'horizon' in the control options of the
## inla() call allows to extend the baseline beyond the last observed event time
## (linear extension, less flexible than the smooth method).## Start sampling## Sampling done.## Computing longitudinal predictions for individual 1 on PID: 27140## Warning in assign(paste0(object$dataSurv), SdataPred): only the first element
## is used as variable name## Warning in assign(paste0(object$dataSurv), SdataPred): only the first element
## is used as variable name## Computing survival predictions for individual 1# plot survival curve for the two time to event outcomes
plot(P$time, P$Surv_quant0.5, type="l", lwd=2, ylim=c(0, 1), xlab="time", ylab="survival probability")
lines(P$time, P$Surv_quant0.025, lty=2)
lines(P$time, P$Surv_quant0.975, lty=2)
# add observed event times
# sapply(surv_data[surv_data$tfs_status==1, "tfs_time"], function(x) abline(v=x, lty=3, lwd=0.5))
P2 <- predict(M2, onePatient, id="id", horizon=30, surv=TRUE)$PredS## Warning in predict.INLAjoint(M2, onePatient, id = "id", horizon = 30, surv =
## TRUE): The fitted model has baseline risk information up until value 23.67 for
## survival outcome 1. Since you ask for prediction at horizon 30 I will assume
## constant baseline hazard beyond the maximum available value. Alternatively, you
## can use baselineHaz='smooth' to use splines to predict the baseline hazard (for
## each sample). Alternatively, adding 'horizon' in the control options of the
## inla() call allows to extend the baseline beyond the last observed event time
## (linear extension, less flexible than the smooth method).## Start sampling## Sampling done.## Computing longitudinal predictions for individual 1 on PID: 27140## Computing survival predictions for individual 1lines(P2$time, P2$Surv_quant0.5, lwd=2, col=2)
lines(P2$time, P2$Surv_quant0.025, lty=2, col=2)
lines(P2$time, P2$Surv_quant0.975, lty=2, col=2)
# sapply(surv_data[surv_data$drfs_time==1, "drfs_time"], function(x) abline(v=x, lty=3, lwd=0.5, col=2))
legend("topright", c("Treatment", "Relapse"), lwd=2, col=1:2)
M-3: Time to PSA (LONG)
Linear
# First model for longitudinal marker PSA
M3 <- joint(formLong=log_psa ~ time_to_psa + (1 + time_to_psa | id),
timeVar="time_to_psa",
id="id",
family="gaussian",
control=list(int.strategy="eb"), # To make calculation faster, use only the mean for the hyper-parameter distributions
dataLong=s_survlong_rprt
)
summary(M3)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 1.3149 0.0085 1.2983 1.3149 1.3314
## time_to_psa -0.0933 0.0021 -0.0975 -0.0933 -0.0891
## Res. err. (variance) 0.4583 0.0018 0.4548 0.4583 0.4618
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.5816 0.0101 0.5617 0.5816 0.6018
## time_to_psa 0.0326 0.0008 0.0311 0.0326 0.0340
## Intercept:time_to_psa 0.0161 0.0021 0.0124 0.0159 0.0208
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -180184.8 -180180.9
##
## Deviance Information Criterion: 327981.8
## Widely applicable Bayesian information criterion: 330535.3
## Computation time: 17.79 seconds# patientID = 25 # pick one
patient_counts <- table(s_survlong_rprt$id)
max_count <- max(patient_counts) # Pick the patient with Max PSA results
patients_with_max_count <- names(patient_counts[patient_counts == max_count])
patientID = patients_with_max_count
ND <- s_survlong_rprt[s_survlong_rprt$id==patientID,]
P1 <- predict(M3, ND, id="id", horizon=max(s_survlong_rprt$time_to_psa))$PredL # Make Linear prediction## Start sampling## Sampling done.## Computing longitudinal predictions for individual 2106 on PID: 27140# To plot, should run from HRE
plot(P1$time_to_psa, P1$quant0.5, type="l", lwd=2, ylim=range(c(P1$quant0.025, P1$quant0.975)), xlab="time", ylab="log (PSA+1)")
lines(P1$time_to_psa, P1$quant0.025, lty=2)
lines(P1$time_to_psa, P1$quant0.975, lty=2)
points(ND$time_to_psa, ND$log_psa, pch=19) # pch=size of the dots
# to HERE Splines
# use splines?
NSplines <- ns(s_survlong_rprt$time_to_psa, knots=c(3)) # natural cubic splines , knots, the starting rise time, look at time to psa summary
f1 <- function(x) predict(NSplines, x)[,1] # first basis
f2 <- function(x) predict(NSplines, x)[,2] # second basis
# check splines
curve(f1, xlim=range(s_survlong_rprt$time_to_psa), ylim=c(-1,1))
curve(f2, xlim=range(s_survlong_rprt$time_to_psa), add=T)
# Second joint model for the longitudinal data, this time accounting for all
M3_cov <- joint(formLong=log_psa ~ (f1(time_to_psa) + f2(time_to_psa)) * grade_group+mri+bx_age+drfs_tx + (1 + f1(time_to_psa) + f2(time_to_psa)| id),
timeVar="time_to_psa",
id="id",
family="gaussian",
control=list(int.strategy="eb"),
dataLong=s_survlong_rprt
)## Warning in joint(formLong = log_psa ~ (f1(time_to_psa) + f2(time_to_psa)) * :
## Internal correlation between hyperparameters is abnormally high, this is a sign
## of identifiability issues / ill-defined model.summary(M3_cov)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 9.0064 0.3502 8.3200 9.0064 9.6929
## f1time_to_psa -14.7400 0.5809 -15.8786 -14.7400 -13.6015
## f2time_to_psa -2.5818 0.2469 -3.0657 -2.5818 -2.0978
## grade_group2 3.6562 0.5271 2.6231 3.6562 4.6893
## grade_group3 7.0605 0.6039 5.8769 7.0605 8.2441
## grade_group4 6.7732 0.8227 5.1607 6.7732 8.3857
## grade_group5 9.3474 0.9187 7.5468 9.3474 11.1481
## mri1 0.1474 0.0211 0.1061 0.1474 0.1887
## bx_age70 0.2369 0.0204 0.1969 0.2369 0.2769
## bx_age6070 0.1086 0.0188 0.0717 0.1086 0.1455
## drfs_txrp -0.5275 0.0172 -0.5612 -0.5275 -0.4939
## drfs_txunk -0.1844 0.0179 -0.2194 -0.1844 -0.1493
## f1time_to_psa:grade_group2 -7.3400 0.8732 -9.0515 -7.3400 -5.6286
## f1time_to_psa:grade_group3 -12.7019 0.9993 -14.6604 -12.7019 -10.7434
## f1time_to_psa:grade_group4 -10.0302 1.3595 -12.6948 -10.0302 -7.3657
## f1time_to_psa:grade_group5 -11.3804 1.5019 -14.3241 -11.3804 -8.4368
## f2time_to_psa:grade_group2 -0.4800 0.3780 -1.2208 -0.4800 0.2607
## f2time_to_psa:grade_group3 1.6475 0.4335 0.7980 1.6475 2.4971
## f2time_to_psa:grade_group4 5.8928 0.6008 4.7153 5.8928 7.0704
## f2time_to_psa:grade_group5 11.6067 0.7072 10.2206 11.6067 12.9929
## Res. err. (variance) 0.3318 0.0006 0.3307 0.3318 0.3330
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 307.4347 8.9810 288.5149 308.0869 323.6489
## f1time_to_psa 867.0852 24.9862 813.9520 869.0045 910.5119
## f2time_to_psa 145.1463 4.1883 135.5651 145.6648 151.7853
## Intercept:f1time_to_psa -509.6712 14.9502 -536.3252 -510.6209 -478.2912
## Intercept:f2time_to_psa 113.3721 4.5694 103.8422 113.6088 121.9315
## f1time_to_psa:f2time_to_psa -141.5667 7.3315 -154.9002 -141.9427 -125.9167
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -170272.5 -170265.8
##
## Deviance Information Criterion: 288060
## Widely applicable Bayesian information criterion: 291625.1
## Computation time: 52.98 seconds# Again, pick one patient (here the same)
ND2 <- s_survlong_rprt[s_survlong_rprt$id==patientID,] # observed vs. fitted for a couple individuals.
P2 <- predict(M3_cov, ND2, id="id", horizon=max(s_survlong_rprt$time_to_psa))$PredL## Start sampling## Sampling done.## Computing longitudinal predictions for individual 2106 on PID: 27140plot(P2$time_to_psa, P2$quant0.5, type="l", lwd=2, ylim=range(c(P2$quant0.025, P2$quant0.975)), xlab="time", ylab="log (PSA+1)")
lines(P2$time_to_psa, P2$quant0.025, lty=2)
lines(P2$time_to_psa, P2$quant0.975, lty=2)
points(ND$time_to_psa, ND$log_psa, pch=19)
Linear + account for all
# use Linear
NSplines <- ns(s_survlong_rprt$time_to_psa, knots=c(3)) # natural cubic splines , knots, the starting rise time, look at time to psa summary
# Second joint model for the longitudinal data, this time accounting for all
M3_cov <- joint(formLong=log_psa ~ (time_to_psa) * grade_group+mri+bx_age+drfs_tx + (1 + time_to_psa| id),
timeVar="time_to_psa",
id="id",
family="gaussian",
control=list(int.strategy="eb"),
dataLong=s_survlong_rprt
)
summary(M3_cov)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 1.4384 0.0227 1.3938 1.4384 1.4829
## time_to_psa -0.0864 0.0034 -0.0931 -0.0864 -0.0797
## grade_group2 -0.2388 0.0199 -0.2777 -0.2388 -0.1999
## grade_group3 -0.2880 0.0227 -0.3325 -0.2880 -0.2434
## grade_group4 -0.1373 0.0306 -0.1972 -0.1373 -0.0775
## grade_group5 0.3171 0.0317 0.2549 0.3171 0.3792
## mri1 0.1448 0.0234 0.0989 0.1448 0.1907
## bx_age70 0.2543 0.0227 0.2098 0.2543 0.2988
## bx_age6070 0.1188 0.0210 0.0777 0.1188 0.1599
## drfs_txrp -0.4839 0.0191 -0.5214 -0.4839 -0.4464
## drfs_txunk -0.1035 0.0199 -0.1424 -0.1035 -0.0645
## time_to_psa:grade_group2 -0.0374 0.0052 -0.0475 -0.0374 -0.0273
## time_to_psa:grade_group3 -0.0332 0.0059 -0.0448 -0.0332 -0.0217
## time_to_psa:grade_group4 0.0389 0.0081 0.0230 0.0389 0.0548
## time_to_psa:grade_group5 0.0992 0.0091 0.0814 0.0992 0.1169
## Res. err. (variance) 0.4585 0.0018 0.4550 0.4585 0.4620
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.4909 0.0088 0.4742 0.4910 0.5090
## time_to_psa 0.0312 0.0007 0.0298 0.0312 0.0326
## Intercept:time_to_psa -0.0009 0.0020 -0.0049 -0.0009 0.0028
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -179414.1 -179410.1
##
## Deviance Information Criterion: 327946.3
## Widely applicable Bayesian information criterion: 330513
## Computation time: 25.8 seconds# Again, pick one patient (here the same)
ND2 <- s_survlong_rprt[s_survlong_rprt$id==patientID,] # observed vs. fitted for a couple individuals.
P2 <- predict(M3_cov, ND2, id="id", horizon=max(s_survlong_rprt$time_to_psa))$PredL## Start sampling## Sampling done.## Computing longitudinal predictions for individual 2106 on PID: 27140plot(P2$time_to_psa, P2$quant0.5, type="l", lwd=2, ylim=range(c(P2$quant0.025, P2$quant0.975)), xlab="time", ylab="log (PSA+1)")
lines(P2$time_to_psa, P2$quant0.025, lty=2)
lines(P2$time_to_psa, P2$quant0.975, lty=2)
points(ND$time_to_psa, ND$log_psa, pch=19)
Quadratic
f1 <- function(x) x*x
#f2 <- function(x) predict(NSplines, x)[,2] # second basis
# check splines
curve(f1, xlim=range(s_survlong_rprt$time_to_psa))
#curve(f2, xlim=range(s_survlong_rprt$time_to_psa), add=T)
# Second joint model for the longitudinal data, this time accounting for all
M3_cov <- joint(formLong=log_psa ~ (time_to_psa + f1(time_to_psa)) * grade_group+mri+bx_age+drfs_tx + (1 + time_to_psa + f1(time_to_psa)| id),
timeVar="time_to_psa",
id="id",
family="gaussian",
control=list(int.strategy="eb"),
dataLong=s_survlong_rprt
)
summary(M3_cov)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 1.6173 0.0226 1.5729 1.6173 1.6616
## time_to_psa -0.1980 0.0060 -0.2097 -0.1980 -0.1863
## f1time_to_psa 0.0097 0.0009 0.0081 0.0097 0.0114
## grade_group2 -0.2123 0.0200 -0.2514 -0.2123 -0.1731
## grade_group3 -0.2397 0.0228 -0.2844 -0.2397 -0.1949
## grade_group4 -0.1057 0.0307 -0.1659 -0.1057 -0.0455
## grade_group5 0.3486 0.0319 0.2862 0.3486 0.4110
## mri1 0.1387 0.0232 0.0932 0.1387 0.1841
## bx_age70 0.2319 0.0225 0.1879 0.2319 0.2760
## bx_age6070 0.1072 0.0207 0.0666 0.1072 0.1479
## drfs_txrp -0.5222 0.0189 -0.5593 -0.5222 -0.4851
## drfs_txunk -0.1848 0.0197 -0.2234 -0.1848 -0.1463
## time_to_psa:grade_group2 -0.0952 0.0089 -0.1127 -0.0952 -0.0778
## time_to_psa:grade_group3 -0.1384 0.0102 -0.1584 -0.1384 -0.1185
## time_to_psa:grade_group4 -0.0635 0.0139 -0.0908 -0.0635 -0.0363
## time_to_psa:grade_group5 -0.0180 0.0153 -0.0479 -0.0180 0.0119
## f1time_to_psa:grade_group2 0.0062 0.0013 0.0037 0.0062 0.0088
## f1time_to_psa:grade_group3 0.0153 0.0015 0.0124 0.0153 0.0182
## f1time_to_psa:grade_group4 0.0208 0.0020 0.0168 0.0208 0.0248
## f1time_to_psa:grade_group5 0.0313 0.0024 0.0266 0.0313 0.0359
## Res. err. (variance) 0.3442 0.0014 0.3414 0.3442 0.3470
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.4971 0.0093 0.4790 0.4971 0.5156
## time_to_psa 0.0936 0.0039 0.0873 0.0932 0.1030
## f1time_to_psa 0.0017 0.0001 0.0016 0.0017 0.0018
## Intercept:time_to_psa -0.0292 0.0048 -0.0401 -0.0288 -0.0209
## Intercept:f1time_to_psa 0.0048 0.0010 0.0032 0.0047 0.0069
## time_to_psa:f1time_to_psa -0.0076 0.0006 -0.0090 -0.0075 -0.0066
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -171290.5 -171283.8
##
## Deviance Information Criterion: 290278.9
## Widely applicable Bayesian information criterion: 293071.7
## Computation time: 73.95 seconds# Again, pick one patient (here the same)
ND2 <- s_survlong_rprt[s_survlong_rprt$id==patientID,] # observed vs. fitted for a couple individuals.
P2 <- predict(M3_cov, ND2, id="id", horizon=max(s_survlong_rprt$time_to_psa))$PredL## Start sampling## Sampling done.## Computing longitudinal predictions for individual 2106 on PID: 27140plot(P2$time_to_psa, P2$quant0.5, type="l", lwd=2, ylim=range(c(P2$quant0.025, P2$quant0.975)), xlab="time", ylab="log (PSA+1)")
lines(P2$time_to_psa, P2$quant0.025, lty=2)
lines(P2$time_to_psa, P2$quant0.975, lty=2)
points(ND$time_to_psa, ND$log_psa, pch=19)
JM-2: Treatment.X.Recurrence.X.PSA (SURV.SURV.LONG)
Effect post RP
############################################## FINAL JOINT MODEL
# Get only RP-Caused Recurrence cases
fm_surlong_rp <- s_survlong_rprt[s_survlong_rprt$drfs_tx == 'rp',]
fm_surlong_rp <- reassigning_ids(fm_surlong_rp)
fm_sur_rp <- fm_surlong_rp[!duplicated(fm_surlong_rp$id), ]
# We need separate df for BCR (only treated patients)
fm_sur_rp_bcr <- fm_sur_rp[fm_sur_rp$tfs_status==1,]
fm_sur_rp_bcr$drfs_tx <- droplevels(fm_sur_rp_bcr$drfs_tx)
fm_sur_rp_bcr$tfs_tx <- droplevels(fm_sur_rp_bcr$tfs_tx)
# Setup the number of used threads
inla.setOption(num.threads='8:1')Splines
NSplines <- ns(fm_surlong_rp$time_to_psa, knots=c(3)) # natural cubic splines , knots, the starting rise time, look at time to psa summary
f1 <- function(x) predict(NSplines, x)[,1] # first basis
f2 <- function(x) predict(NSplines, x)[,2] # second basis
# Final joint model, combining 2 survival modals and the longitudinal one
JM2_RP_splines <- joint(formSurv=list(inla.surv(tfs_time, tfs_status) ~ grade_group+mri+bx_age + (1|id),
inla.surv(drfs_time, drfs_status, truncation = tfs_time) ~ grade_group+mri+bx_age),
formLong=log_psa ~ (f1(time_to_psa) + f2(time_to_psa)) * grade_group+mri+bx_age + (1 + f1(time_to_psa) + f2(time_to_psa)| id),
id="id",
timeVar="time_to_psa",
basRisk=c("rw2", "rw2"), # random walk
assocSurv=TRUE,
family="gaussian",
assoc=c("CV", "CV"), # what is cv : Current value = risk event at t depends on longit.marker(psa) at t
control=list(int.strategy="eb"),
dataSurv=list(fm_sur_rp, fm_sur_rp_bcr),
dataLong=fm_surlong_rp
)## Warning in joint(formSurv = list(inla.surv(tfs_time, tfs_status) ~ grade_group
## + : Stupid local search strategy used: This can be a sign of a ill-defined
## model and/or non-informative data.summary(JM2_RP_splines)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 2.8404 0.1680 2.5110 2.8404 3.1697
## f1time_to_psa_L1 -5.8638 0.2479 -6.3496 -5.8638 -5.3780
## f2time_to_psa_L1 -6.5888 0.2383 -7.0558 -6.5888 -6.1218
## grade_group2_L1 -0.1482 0.2242 -0.5875 -0.1482 0.2912
## grade_group3_L1 -0.4321 0.2387 -0.8999 -0.4321 0.0357
## grade_group4_L1 -0.4179 0.3126 -1.0306 -0.4179 0.1948
## grade_group5_L1 -0.2429 0.3696 -0.9673 -0.2429 0.4815
## mri1_L1 0.0070 0.0188 -0.0297 0.0070 0.0438
## bx_age70_L1 0.0891 0.0225 0.0450 0.0891 0.1332
## bx_age6070_L1 0.0788 0.0174 0.0446 0.0788 0.1130
## f1time_to_psa:grade_group2_L1 -0.5336 0.3265 -1.1737 -0.5336 0.1064
## f1time_to_psa:grade_group3_L1 0.2047 0.3431 -0.4678 0.2047 0.8772
## f1time_to_psa:grade_group4_L1 0.5658 0.4372 -0.2912 0.5658 1.4228
## f1time_to_psa:grade_group5_L1 0.9263 0.5102 -0.0738 0.9263 1.9263
## f2time_to_psa:grade_group2_L1 -0.7008 0.3312 -1.3500 -0.7008 -0.0516
## f2time_to_psa:grade_group3_L1 -0.0641 0.3649 -0.7794 -0.0641 0.6512
## f2time_to_psa:grade_group4_L1 0.7569 0.5063 -0.2355 0.7569 1.7492
## f2time_to_psa:grade_group5_L1 2.1118 0.6133 0.9097 2.1118 3.3140
## Res. err. (variance) 0.4780 0.0040 0.4690 0.4785 0.4842
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant
## Intercept_L1 7.7197 0.5584 6.7550 7.6813
## f1time_to_psa_L1 8.3224 0.8141 6.8731 8.2681
## f2time_to_psa_L1 33.2906 1.8899 30.1685 33.0984
## Intercept_L1:f1time_to_psa_L1 -7.9264 0.6729 -9.3501 -7.8904
## Intercept_L1:f2time_to_psa_L1 15.9719 0.9523 14.3591 15.8867
## f1time_to_psa_L1:f2time_to_psa_L1 -16.5483 1.1098 -18.8516 -16.4667
## 0.975quant
## Intercept_L1 8.8874
## f1time_to_psa_L1 9.9952
## f2time_to_psa_L1 37.5049
## Intercept_L1:f1time_to_psa_L1 -6.7319
## Intercept_L1:f2time_to_psa_L1 18.0517
## f1time_to_psa_L1:f2time_to_psa_L1 -14.5888
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0261 0.0033 0.0197 0.0262 0.0324
## grade_group2_S1 0.4642 0.0461 0.3738 0.4642 0.5546
## grade_group3_S1 0.5874 0.0514 0.4866 0.5874 0.6881
## grade_group4_S1 0.7491 0.0691 0.6137 0.7491 0.8845
## grade_group5_S1 0.9394 0.0799 0.7828 0.9394 1.0960
## mri1_S1 0.3502 0.0441 0.2637 0.3502 0.4366
## bx_age70_S1 0.1975 0.0530 0.0935 0.1975 0.3015
## bx_age6070_S1 0.1086 0.0407 0.0288 0.1086 0.1884
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.1352 0.0073 0.121 0.1352 0.1499
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0001 0.0000 0.0000 0.0001 0.0001
## grade_group2_S2 1.0486 0.1412 0.7719 1.0486 1.3252
## grade_group3_S2 1.6545 0.1565 1.3477 1.6545 1.9612
## grade_group4_S2 2.3972 0.2135 1.9788 2.3972 2.8157
## grade_group5_S2 3.6706 0.2484 3.1837 3.6706 4.1574
## mri1_S2 1.4868 0.1413 1.2098 1.4868 1.7638
## bx_age70_S2 0.5215 0.1680 0.1923 0.5215 0.8506
## bx_age6070_S2 -0.0788 0.1285 -0.3308 -0.0788 0.1731
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -1.5816 0.0495 -1.6703 -1.5840 -1.4762
## CV_L1_S2 -1.2876 0.0745 -1.4431 -1.2848 -1.1501
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 7.1552 0.1177 6.956 7.1464 7.4137
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -138621.9 -138609.7
##
## Deviance Information Criterion: -12158.52
## Widely applicable Bayesian information criterion: -16407.13
## Computation time: 387.12 secondsLinear
# Final joint model, combining 2 survival modals and the longitudinal one
JM2_RP_linear <- joint(formSurv=list(inla.surv(tfs_time, tfs_status) ~ grade_group+mri+bx_age + (1|id),
inla.surv(drfs_time, drfs_status, truncation = tfs_time) ~ grade_group+mri+bx_age),
formLong=log_psa ~ (time_to_psa) * grade_group+mri+bx_age + (1 + time_to_psa| id),
id="id",
timeVar="time_to_psa",
basRisk=c("rw2", "rw2"), # random walk
assocSurv=TRUE,
family="gaussian",
assoc=c("CV", "CV"), # what is cv : Current value = risk event at t depends on longit.marker(psa) at t
control=list(int.strategy="eb"),
dataSurv=list(fm_sur_rp, fm_sur_rp_bcr),
dataLong=fm_surlong_rp
)## Warning in joint(formSurv = list(inla.surv(tfs_time, tfs_status) ~ grade_group
## + : Stupid local search strategy used: This can be a sign of a ill-defined
## model and/or non-informative data.summary(JM2_RP_linear)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 1.1292 0.0193 1.0914 1.1292 1.1670
## time_to_psa_L1 -0.1816 0.0048 -0.1910 -0.1816 -0.1721
## grade_group2_L1 -0.3130 0.0206 -0.3533 -0.3130 -0.2727
## grade_group3_L1 -0.3372 0.0227 -0.3817 -0.3372 -0.2926
## grade_group4_L1 -0.3127 0.0309 -0.3732 -0.3127 -0.2522
## grade_group5_L1 -0.2448 0.0361 -0.3156 -0.2448 -0.1739
## mri1_L1 0.0240 0.0190 -0.0132 0.0240 0.0612
## bx_age70_L1 0.1013 0.0227 0.0569 0.1013 0.1457
## bx_age6070_L1 0.0801 0.0176 0.0456 0.0801 0.1146
## time_to_psa:grade_group2_L1 0.0070 0.0066 -0.0059 0.0070 0.0200
## time_to_psa:grade_group3_L1 0.0260 0.0072 0.0119 0.0260 0.0402
## time_to_psa:grade_group4_L1 0.0481 0.0099 0.0287 0.0481 0.0675
## time_to_psa:grade_group5_L1 0.0804 0.0119 0.0571 0.0804 0.1037
## Res. err. (variance) 0.4957 0.0037 0.4895 0.4953 0.5038
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 0.1157 0.0050 0.1066 0.1154 0.1265
## time_to_psa_L1 0.0121 0.0006 0.0109 0.0122 0.0132
## Intercept_L1:time_to_psa_L1 -0.0128 0.0011 -0.0152 -0.0128 -0.0108
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0164 0.0019 0.0133 0.0162 0.0207
## grade_group2_S1 0.4583 0.0446 0.3708 0.4583 0.5457
## grade_group3_S1 0.5852 0.0497 0.4877 0.5852 0.6827
## grade_group4_S1 0.7586 0.0666 0.6281 0.7586 0.8892
## grade_group5_S1 0.9518 0.0772 0.8005 0.9518 1.1030
## mri1_S1 0.2780 0.0426 0.1944 0.2780 0.3615
## bx_age70_S1 0.1468 0.0512 0.0465 0.1468 0.2471
## bx_age6070_S1 0.0912 0.0395 0.0139 0.0912 0.1686
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.0964 0.0045 0.0885 0.096 0.1062
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0915 0.0170 0.0657 0.0888 0.1321
## grade_group2_S2 1.4151 0.1553 1.1108 1.4151 1.7195
## grade_group3_S2 2.0347 0.1722 1.6971 2.0347 2.3722
## grade_group4_S2 2.8047 0.2355 2.3430 2.8047 3.2663
## grade_group5_S2 4.1424 0.2747 3.6041 4.1424 4.6807
## mri1_S2 1.7381 0.1562 1.4321 1.7381 2.0442
## bx_age70_S2 0.5866 0.1855 0.2231 0.5866 0.9501
## bx_age6070_S2 -0.1030 0.1417 -0.3808 -0.1030 0.1748
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -1.5290 0.0710 -1.6457 -1.5354 -1.3710
## CV_L1_S2 -0.1765 0.1371 -0.4629 -0.1712 0.0762
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 9.5878 0.1598 9.3084 9.5779 9.9324
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -139230.8 -139221.4
##
## Deviance Information Criterion: -10427.88
## Widely applicable Bayesian information criterion: -14421.25
## Computation time: 254.66 secondsMore Combinations
f1 <- function(x) x*x
# Final joint model, combining 2 survival modals and the longitudinal one
JM2_RP_quad <- joint(formSurv=list(inla.surv(tfs_time, tfs_status) ~ grade_group+mri+bx_age + (1|id),
inla.surv(drfs_time, drfs_status, truncation = tfs_time) ~ grade_group+mri+bx_age),
formLong=log_psa ~ (time_to_psa + f1(time_to_psa)) * grade_group+mri+bx_age + (1 + time_to_psa + f1(time_to_psa)| id),
id="id",
timeVar="time_to_psa",
basRisk=c("rw2", "rw2"), # random walk
assocSurv=TRUE,
family="gaussian",
assoc=c("CV", "CV"), # what is cv : Current value = risk event at t depends on longit.marker(psa) at t
control=list(int.strategy="eb"),
dataSurv=list(fm_sur_rp, fm_sur_rp_bcr),
dataLong=fm_surlong_rp
)
summary(JM2_RP_quad)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 1.3375 0.0209 1.2966 1.3375 1.3784
## time_to_psa_L1 -0.3726 0.0108 -0.3937 -0.3726 -0.3514
## f1time_to_psa_L1 0.0208 0.0021 0.0167 0.0208 0.0249
## grade_group2_L1 -0.3325 0.0229 -0.3774 -0.3325 -0.2877
## grade_group3_L1 -0.3560 0.0253 -0.4056 -0.3560 -0.3063
## grade_group4_L1 -0.3479 0.0345 -0.4155 -0.3479 -0.2803
## grade_group5_L1 -0.3017 0.0405 -0.3811 -0.3017 -0.2222
## mri1_L1 0.0131 0.0195 -0.0251 0.0131 0.0514
## bx_age70_L1 0.0889 0.0235 0.0429 0.0889 0.1349
## bx_age6070_L1 0.0689 0.0181 0.0333 0.0689 0.1044
## time_to_psa:grade_group2_L1 -0.0231 0.0146 -0.0516 -0.0231 0.0054
## time_to_psa:grade_group3_L1 -0.0063 0.0160 -0.0377 -0.0063 0.0250
## time_to_psa:grade_group4_L1 0.0459 0.0217 0.0034 0.0459 0.0883
## time_to_psa:grade_group5_L1 0.0675 0.0258 0.0169 0.0675 0.1181
## f1time_to_psa:grade_group2_L1 0.0071 0.0029 0.0015 0.0071 0.0128
## f1time_to_psa:grade_group3_L1 0.0079 0.0032 0.0017 0.0079 0.0141
## f1time_to_psa:grade_group4_L1 0.0031 0.0043 -0.0054 0.0031 0.0115
## f1time_to_psa:grade_group5_L1 0.0086 0.0052 -0.0017 0.0086 0.0188
## Res. err. (variance) 0.3882 0.0029 0.3826 0.3882 0.3939
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 0.1456 0.0064 0.1338 0.1454 0.1590
## time_to_psa_L1 0.0606 0.0028 0.0553 0.0605 0.0665
## f1time_to_psa_L1 0.0025 0.0001 0.0022 0.0025 0.0027
## Intercept_L1:time_to_psa_L1 -0.0406 0.0037 -0.0480 -0.0406 -0.0339
## Intercept_L1:f1time_to_psa_L1 0.0028 0.0005 0.0018 0.0028 0.0039
## time_to_psa_L1:f1time_to_psa_L1 -0.0083 0.0005 -0.0092 -0.0083 -0.0075
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.2088 0.0449 0.1348 0.2041 0.3106
## grade_group2_S1 0.4295 0.0444 0.3424 0.4295 0.5165
## grade_group3_S1 0.5595 0.0496 0.4623 0.5595 0.6567
## grade_group4_S1 0.7118 0.0662 0.5820 0.7118 0.8415
## grade_group5_S1 0.8723 0.0767 0.7221 0.8723 1.0226
## mri1_S1 0.2830 0.0420 0.2007 0.2830 0.3653
## bx_age70_S1 0.1208 0.0506 0.0217 0.1208 0.2199
## bx_age6070_S1 0.0704 0.0388 -0.0057 0.0704 0.1465
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.0844 0.0055 0.0741 0.0843 0.0957
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.2316 0.0515 0.1485 0.2254 0.3503
## grade_group2_S2 1.5509 0.1628 1.2318 1.5509 1.8700
## grade_group3_S2 2.1782 0.1806 1.8242 2.1782 2.5322
## grade_group4_S2 2.9578 0.2472 2.4734 2.9578 3.4423
## grade_group5_S2 4.3312 0.2883 3.7661 4.3312 4.8962
## mri1_S2 1.8561 0.1639 1.5349 1.8561 2.1773
## bx_age70_S2 0.6074 0.1946 0.2260 0.6074 0.9888
## bx_age6070_S2 -0.1139 0.1487 -0.4053 -0.1139 0.1775
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -1.4882 0.0567 -1.5996 -1.4883 -1.3764
## CV_L1_S2 0.2540 0.1125 0.0328 0.2540 0.4757
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 10.6123 0.2132 10.1753 10.6181 11.0141
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -137558.5 -137546.3
##
## Deviance Information Criterion: -20813.53
## Widely applicable Bayesian information criterion: -24849.32
## Computation time: 439.54 secondsEffect post RT
Assessing the effect of time to treatment on risk of relapse for patients having RadioTherapy as there first Curative Treatment
# Get only RP-Caused Recurrence cases
fm_surlong_rt <- s_survlong_rprt[s_survlong_rprt$drfs_tx == 'rt',]
fm_surlong_rt <- reassigning_ids(fm_surlong_rt)
fm_sur_rt <- fm_surlong_rt[!duplicated(fm_surlong_rt$id), ]
# We need separate df for BCR (only treated patients)
fm_sur_rt_bcr <- fm_sur_rt[fm_sur_rt$tfs_status==1,]
fm_sur_rt_bcr$drfs_tx <- droplevels(fm_sur_rt_bcr$drfs_tx)
fm_sur_rt_bcr$tfs_tx <- droplevels(fm_sur_rt_bcr$tfs_tx)
# Setup the number of used threads
inla.setOption(num.threads='8:1')Splines
NSplines <- ns(fm_surlong_rt$time_to_psa, knots=c(3)) # natural cubic splines , knots, the starting rise time, look at time to psa summary
f1 <- function(x) predict(NSplines, x)[,1] # first basis
f2 <- function(x) predict(NSplines, x)[,2] # second basis
# Final joint model, combining 2 survival modals and the longitudinal one
JM2_RP_splines <- joint(formSurv=list(inla.surv(tfs_time, tfs_status) ~ grade_group+mri+bx_age + (1|id),
inla.surv(drfs_time, drfs_status, truncation = tfs_time) ~ grade_group+mri+bx_age),
formLong=log_psa ~ (f1(time_to_psa) + f2(time_to_psa)) * grade_group+mri+bx_age + (1 + f1(time_to_psa) + f2(time_to_psa)| id),
id="id",
timeVar="time_to_psa",
basRisk=c("rw2", "rw2"), # random walk
assocSurv=TRUE,
family="gaussian",
assoc=c("CV", "CV"), # what is cv : Current value = risk event at t depends on longit.marker(psa) at t
control=list(int.strategy="eb"),
dataSurv=list(fm_sur_rt, fm_sur_rt_bcr),
dataLong=fm_surlong_rt
)## Warning in joint(formSurv = list(inla.surv(tfs_time, tfs_status) ~ grade_group
## + : Stupid local search strategy used: This can be a sign of a ill-defined
## model and/or non-informative data.## Warning in joint(formSurv = list(inla.surv(tfs_time, tfs_status) ~ grade_group
## + : Internal correlation between hyperparameters is abnormally high, this is a
## sign of identifiability issues / ill-defined model.summary(JM2_RP_splines)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 9.7818 0.5154 8.7716 9.7818 10.7919
## f1time_to_psa_L1 -16.2564 0.8656 -17.9529 -16.2564 -14.5598
## f2time_to_psa_L1 -3.0294 0.3889 -3.7916 -3.0294 -2.2671
## grade_group2_L1 1.8496 0.7663 0.3476 1.8496 3.3516
## grade_group3_L1 6.5645 0.8517 4.8953 6.5645 8.2338
## grade_group4_L1 7.1007 1.1279 4.8901 7.1007 9.3113
## grade_group5_L1 10.4734 1.2225 8.0773 10.4734 12.8695
## mri1_L1 0.1564 0.0450 0.0682 0.1564 0.2446
## bx_age70_L1 0.0917 0.0371 0.0189 0.0917 0.1645
## bx_age6070_L1 0.0457 0.0369 -0.0267 0.0457 0.1181
## f1time_to_psa:grade_group2_L1 -4.1136 1.2881 -6.6382 -4.1136 -1.5891
## f1time_to_psa:grade_group3_L1 -11.0015 1.4289 -13.8020 -11.0015 -8.2010
## f1time_to_psa:grade_group4_L1 -8.0894 1.8845 -11.7830 -8.0894 -4.3958
## f1time_to_psa:grade_group5_L1 -11.1488 2.0192 -15.1064 -11.1488 -7.1912
## f2time_to_psa:grade_group2_L1 -0.6221 0.5856 -1.7699 -0.6221 0.5257
## f2time_to_psa:grade_group3_L1 3.3077 0.6518 2.0302 3.3077 4.5852
## f2time_to_psa:grade_group4_L1 11.1206 0.8783 9.3993 11.1206 12.8420
## f2time_to_psa:grade_group5_L1 16.4294 0.9868 14.4953 16.4294 18.3634
## Res. err. (variance) 0.4034 0.0023 0.3993 0.4033 0.4082
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant
## Intercept_L1 287.1197 11.9725 265.0268 286.6572
## f1time_to_psa_L1 838.3426 32.2381 778.5791 837.3761
## f2time_to_psa_L1 155.0039 8.6327 137.7685 154.9908
## Intercept_L1:f1time_to_psa_L1 -482.4813 19.4490 -523.6496 -481.6402
## Intercept_L1:f2time_to_psa_L1 91.1177 8.4901 75.0479 90.8070
## f1time_to_psa_L1:f2time_to_psa_L1 -96.1587 12.7745 -121.0214 -95.7986
## 0.975quant
## Intercept_L1 312.3467
## f1time_to_psa_L1 905.8890
## f2time_to_psa_L1 171.1777
## Intercept_L1:f1time_to_psa_L1 -446.7611
## Intercept_L1:f2time_to_psa_L1 107.7523
## f1time_to_psa_L1:f2time_to_psa_L1 -71.6411
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0659 0.0065 0.0522 0.0665 0.0769
## grade_group2_S1 0.6398 0.0498 0.5421 0.6398 0.7375
## grade_group3_S1 0.7672 0.0547 0.6600 0.7672 0.8743
## grade_group4_S1 0.7212 0.0705 0.5830 0.7212 0.8593
## grade_group5_S1 0.8704 0.0706 0.7322 0.8704 1.0087
## mri1_S1 0.3909 0.0694 0.2549 0.3909 0.5269
## bx_age70_S1 0.4966 0.0574 0.3841 0.4966 0.6092
## bx_age6070_S1 0.2976 0.0571 0.1857 0.2976 0.4094
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.8512 0.0478 0.7483 0.8559 0.9319
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0613 0.0055 0.0530 0.0604 0.0742
## grade_group2_S2 0.3420 0.0632 0.2181 0.3420 0.4659
## grade_group3_S2 0.6039 0.0690 0.4687 0.6039 0.7390
## grade_group4_S2 0.4757 0.0886 0.3020 0.4757 0.6494
## grade_group5_S2 0.9594 0.0887 0.7855 0.9594 1.1333
## mri1_S2 1.8715 0.0918 1.6915 1.8715 2.0514
## bx_age70_S2 0.5820 0.0729 0.4391 0.5820 0.7250
## bx_age6070_S2 0.2427 0.0720 0.1016 0.2427 0.3838
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -0.1620 0.0181 -0.1999 -0.1613 -0.1286
## CV_L1_S2 0.9776 0.0159 0.9465 0.9775 1.0092
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 1.2546 0.0279 1.206 1.2529 1.315
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -206920.1 -206907.8
##
## Deviance Information Criterion: 11587.01
## Widely applicable Bayesian information criterion: 7455.239
## Computation time: 302.71 secondsInterpretation
cat('\nAFTER RP - splines: \n')##
## AFTER RP - splines:int_res <- interp(JM2_RP_splines)
int_res[1]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the lowest risk of receiving treatment (i.e., longest time-to-treatment) is associated to a 0.39 [0.07,0.94] increased risk of recurrence (i.e., shorter time-to-recurrence) compared to the average individual."int_res[2]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the highest risk of receiving treatment (i.e., shortest time-to-treatment) is associated to a 2.54 [1.07,14.94] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."cat('\nAFTER RP - linear: \n')##
## AFTER RP - linear:int_res <- interp(JM2_RP_linear)
int_res[1]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the lowest risk of receiving treatment (i.e., longest time-to-treatment) is associated to a 0.09 [0,0.83] increased risk of recurrence (i.e., shorter time-to-recurrence) compared to the average individual."int_res[2]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the highest risk of receiving treatment (i.e., shortest time-to-treatment) is associated to a 11.24 [1.2,1078.93] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."cat('\nAFTER RP - Quadratic: \n')##
## AFTER RP - Quadratic:int_res <- interp(JM2_RP_quad)
int_res[1]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the lowest risk of receiving treatment (i.e., longest time-to-treatment) is associated to a 0.09 [0,0.83] increased risk of recurrence (i.e., shorter time-to-recurrence) compared to the average individual."int_res[2]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the highest risk of receiving treatment (i.e., shortest time-to-treatment) is associated to a 11.73 [1.21,1398.03] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."cat('\n-----------------------------\n')##
## -----------------------------cat('\nAFTER RT - splines: \n')##
## AFTER RT - splines:int_res <- interp(JM2_RP_splines)
int_res[1]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the lowest risk of receiving treatment (i.e., longest time-to-treatment) is associated to a 0.39 [0.07,0.93] increased risk of recurrence (i.e., shorter time-to-recurrence) compared to the average individual."int_res[2]## [1] "Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the highest risk of receiving treatment (i.e., shortest time-to-treatment) is associated to a 2.53 [1.07,14.92] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."# cat('\nAFTER RT - linear: \n')
# int_res <- interp(JM2_RP_linear)
# int_res[1]
# int_res[2]
# cat('\nAFTER RT - Quadratic: \n')
# int_res <- interp(JM2_RP_quad)
# int_res[1]
# int_res[2]NOTE:
The comparison assume the same value of covariates and the same PSA level. This interpretation is valid under some assumptions:
Proportional hazards (i.e., interpretation is valid regardless of the category of covariates and PSA level, as long as we compare individuals with the same)
Model choice is appropriate for the data (e.g., you want enough flexibility to capture complex trajectories of PSA)
There are no missing confounders (i.e., factors that affect both the time-to-treatment and the time-to-recurrence and that are not included in the model)
Trying different risk associations
# Soon!