joint_models_report_part3
Abderrahim Oussama Batouche
2024-02-13
Intro:
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 occured
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.
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 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.5Functions
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%
mean_up15 <- apply(SMP, 2, function(x) mean(x[-1][x[-1]>x[1]])) # mean frailty deviation for top 15%
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, 0.5)
# Surv 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 :14856 0:67618 1:29181 rt :36149
## 1st Qu.:1076 >70 :24079 1: 8503 2:21484 rp :23482
## Median :2103 60-70:37186 3:14683 unk:16490
## Mean :2149 4: 5909
## 3rd Qu.:3196 5: 4864
## Max. :4512
## tx_date tfs_time tfs_status drfs_tx drfs_time
## Length:76121 Min. : 0.50 0:16490 rt :36149 Min. : 1.00
## Class :character 1st Qu.: 3.00 1:59631 rp :23482 1st Qu.: 44.00
## Mode :character Median : 6.00 unk:16490 Median : 78.00
## Mean : 32.16 Mean : 85.87
## 3rd Qu.: 41.00 3rd Qu.:121.00
## Max. :259.00 Max. :284.00
## drfs_status psa_results psa_date time_to_psa
## 0:53416 Min. : 0.00 Length:76121 Min. :-242.00
## 1:22705 1st Qu.: 0.05 Class :character 1st Qu.: 3.00
## Median : 1.05 Mode :character Median : 26.00
## Mean : 12.19 Mean : 33.63
## 3rd Qu.: 5.56 3rd Qu.: 63.00
## Max. :76800.00 Max. : 284.00
## log_psa
## Min. : 0.00000
## 1st Qu.: 0.04879
## Median : 0.71784
## Mean : 1.06614
## 3rd Qu.: 1.88099
## Max. :11.24897M-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.9693 0.3181 0.4951 0.9189 1.7368
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -18310.98 -18310.98
##
## Deviance Information Criterion: 36538.26
## Widely applicable Bayesian information criterion: 36542.92
## Computation time: 1.81 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+tfs_tx),
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.5157 0.2054 0.2266 0.4780 1.0257
## bx_age70 0.2391 0.0331 0.1741 0.2391 0.3041
## bx_age6070 0.0376 0.0300 -0.0213 0.0376 0.0964
## mri1 0.4453 0.0345 0.3776 0.4453 0.5130
## grade_group2 0.5104 0.0287 0.4542 0.5104 0.5666
## grade_group3 0.7391 0.0320 0.6763 0.7391 0.8019
## grade_group4 0.6102 0.0444 0.5232 0.6102 0.6973
## grade_group5 0.8813 0.0464 0.7904 0.8813 0.9722
## tfs_txrp 0.6602 0.0268 0.6076 0.6602 0.7127
## tfs_txunk -2.1442 0.0365 -2.2158 -2.1442 -2.0727
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -14345.67 -14345.67
##
## Deviance Information Criterion: 28534.59
## Widely applicable Bayesian information criterion: 28546.39
## Computation time: 3.07 secondsplot(M1_cov)## $Outcomes
## $Outcomes$S1
##
##
## $Baseline
##
## attr(,"class")
## [1] "plot.INLAjoint" "list"M-2: Time to Recurrence (SURV)
# No covariates
M2 <- joint(formSurv =list(inla.surv(time = drfs_time, event = drfs_status) ~ -1),
basRisk = c("rw2"),
dataSurv = list(s_surv_rprt))
summary(M2)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.2648 0.1243 0.0991 0.2394 0.5805
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -18427.11 -18427.11
##
## Deviance Information Criterion: 36788.79
## Widely applicable Bayesian information criterion: 36791.31
## Computation time: 2.2 secondsplot(M2)## $Baseline
##
## attr(,"class")
## [1] "plot.INLAjoint" "list"# With Covariates
M2_cov <- joint(formSurv =list(inla.surv(time = drfs_time, event = drfs_status) ~ bx_age+mri+grade_group+tfs_tx+tfs_time+drfs_tx),
basRisk = c("rw2"),
dataSurv = list(s_surv_rprt))
summary(M2_cov)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.2996 0.1364 0.1156 0.2724 0.6444
## bx_age70 0.2669 0.0401 0.1882 0.2669 0.3455
## bx_age6070 0.0501 0.0359 -0.0203 0.0501 0.1205
## mri1 0.7919 0.0447 0.7043 0.7919 0.8795
## grade_group2 0.0996 0.0345 0.0319 0.0996 0.1673
## grade_group3 0.3159 0.0384 0.2406 0.3159 0.3912
## grade_group4 0.4556 0.0512 0.3552 0.4556 0.5559
## grade_group5 0.9199 0.0521 0.8177 0.9199 1.0221
## tfs_txrp 0.2147 7.0710 -13.6509 0.2147 14.0803
## tfs_txunk 0.6850 7.0711 -13.1807 0.6850 14.5506
## tfs_time -0.0162 0.0005 -0.0171 -0.0162 -0.0152
## drfs_txrp 0.2147 7.0710 -13.6509 0.2147 14.0803
## drfs_txunk 0.6850 7.0711 -13.1807 0.6850 14.5506
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -17130.59 -17130.59
##
## Deviance Information Criterion: 34092.34
## Widely applicable Bayesian information criterion: 34101.34
## Computation time: 3.88 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 = 30,
dataSurv = list(s_surv_rprt))## Warning in joint(formSurv = list(inla.surv(time = tfs_time, event = tfs_status)
## ~ : Stupid local search strategy used: This can be a sign of a ill-defined
## model and/or non-informative data.summary(JM1)##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.9067 0.1222 0.7148 0.8893 1.1913
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 2.048 0.089 1.8945 2.0399 2.2424
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0608 0.0077 0.0471 0.0603 0.0772
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 -0.6063 0.0176 -0.6424 -0.6058 -0.573
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -36921.4 -36921.4
##
## Deviance Information Criterion: 70589.07
## Widely applicable Bayesian information criterion: 72489.89
## Computation time: 10.3 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+tfs_tx+(1|id),
inla.surv(time = drfs_time, event = drfs_status) ~ bx_age+mri+grade_group+tfs_tx+tfs_time+drfs_tx), id="id",
basRisk=c("rw2", "rw2"), assocSurv=TRUE, NbasRisk = 30,
dataSurv = list(s_surv_rprt))
summary(JM1_cov)##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 1.3340 0.4119 0.6640 1.2920 2.2671
## bx_age70_S1 0.2966 0.0382 0.2216 0.2966 0.3715
## bx_age6070_S1 0.0607 0.0346 -0.0071 0.0607 0.1285
## mri1_S1 0.5342 0.0398 0.4562 0.5342 0.6123
## grade_group2_S1 0.4165 0.0331 0.3515 0.4165 0.4815
## grade_group3_S1 0.6053 0.0374 0.5320 0.6053 0.6785
## grade_group4_S1 0.5302 0.0513 0.4295 0.5302 0.6309
## grade_group5_S1 0.7454 0.0535 0.6405 0.7454 0.8503
## tfs_txrp_S1 0.5082 0.0312 0.4469 0.5082 0.5694
## tfs_txunk_S1 -2.3560 0.0423 -2.4390 -2.3560 -2.2733
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.1831 0.0171 0.1522 0.1822 0.2193
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 1.6888 0.3444 1.1104 1.6560 2.4606
## bx_age70_S2 0.1300 0.2078 -0.2779 0.1302 0.5372
## bx_age6070_S2 -0.0403 0.1880 -0.4091 -0.0403 0.3283
## mri1_S2 1.7378 0.2176 1.3113 1.7378 2.1646
## grade_group2_S2 -0.0738 0.1796 -0.4262 -0.0737 0.2782
## grade_group3_S2 0.4976 0.2050 0.0955 0.4976 0.8994
## grade_group4_S2 1.4930 0.2826 0.9390 1.4929 2.0476
## grade_group5_S2 2.8737 0.2930 2.2996 2.8735 3.4487
## tfs_txrp_S2 0.4240 7.0716 -13.4427 0.4240 14.2906
## tfs_txunk_S2 4.6636 7.0723 -9.2044 4.6636 18.5316
## tfs_time_S2 -0.1401 0.0031 -0.1462 -0.1400 -0.1341
## drfs_txrp_S2 0.4240 7.0716 -13.4427 0.4240 14.2906
## drfs_txunk_S2 4.6636 7.0723 -9.2044 4.6636 18.5316
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 -10.8992 0.2987 -11.4681 -10.9054 -10.2926
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -35180.23 -35180.23
##
## Deviance Information Criterion: 61120.48
## Widely applicable Bayesian information criterion: 64888.68
## Computation time: 13.56 secondsplot(JM1_cov)$Baseline+scale_y_log10()
# Survival curves
onePatient <- s_surv_rprt[1, ]
P <- predict(M1, onePatient, id="id", horizon=300, surv=TRUE)$PredS## Warning in predict.INLAjoint(M1, onePatient, id = "id", horizon = 300, surv =
## TRUE): The fitted model has baseline risk information up until value 259 for
## survival outcome 1. Since you ask for prediction at horizon 300 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: 71801## 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=300, surv=TRUE)$PredS## Warning in predict.INLAjoint(M2, onePatient, id = "id", horizon = 300, surv =
## TRUE): The fitted model has baseline risk information up until value 284 for
## survival outcome 1. Since you ask for prediction at horizon 300 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: 71801## 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 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
)## Warning in joint(formLong = log_psa ~ time_to_psa + (1 + time_to_psa | id), :
## Stupid local search strategy used: This can be a sign of a ill-defined model
## and/or non-informative data.summary(M3)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 1.3014 0.0118 1.2784 1.3014 1.3244
## time_to_psa -0.0086 0.0004 -0.0094 -0.0086 -0.0079
## Res. err. (variance) 0.4447 0.0024 0.4398 0.4448 0.4493
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.5543 0.0140 0.5269 0.5544 0.5811
## time_to_psa 0.0006 0.0000 0.0006 0.0006 0.0006
## Intercept:time_to_psa 0.0002 0.0004 -0.0005 0.0001 0.0009
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -90910.24 -90906.27
##
## Deviance Information Criterion: 162109.6
## Widely applicable Bayesian information criterion: 163719.8
## Computation time: 11.34 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 1677 on PID: 71801# 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(25)) # 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
)
summary(M3_cov)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 8.6140 0.4145 7.8016 8.6140 9.4264
## f1time_to_psa -14.1207 0.6941 -15.4812 -14.1207 -12.7602
## f2time_to_psa -2.9951 0.2964 -3.5759 -2.9951 -2.4142
## grade_group2 3.4230 0.6270 2.1941 3.4230 4.6518
## grade_group3 5.8318 0.7129 4.4345 5.8318 7.2290
## grade_group4 5.2901 1.0067 3.3170 5.2901 7.2633
## grade_group5 7.9869 1.1145 5.8025 7.9869 10.1712
## mri1 0.1197 0.0320 0.0570 0.1197 0.1825
## bx_age70 0.2431 0.0316 0.1811 0.2431 0.3051
## bx_age6070 0.1165 0.0290 0.0597 0.1165 0.1733
## drfs_txrp -0.5109 0.0267 -0.5632 -0.5109 -0.4586
## drfs_txunk -0.1673 0.0275 -0.2213 -0.1673 -0.1134
## f1time_to_psa:grade_group2 -6.8841 1.0499 -8.9418 -6.8841 -4.8264
## f1time_to_psa:grade_group3 -10.8915 1.1925 -13.2287 -10.8915 -8.5542
## f1time_to_psa:grade_group4 -7.9107 1.6771 -11.1978 -7.9107 -4.6237
## f1time_to_psa:grade_group5 -9.3998 1.8454 -13.0168 -9.3998 -5.7829
## f2time_to_psa:grade_group2 -0.5049 0.4543 -1.3954 -0.5049 0.3855
## f2time_to_psa:grade_group3 0.4763 0.5158 -0.5346 0.4763 1.4871
## f2time_to_psa:grade_group4 4.9015 0.7523 3.4269 4.9015 6.3760
## f2time_to_psa:grade_group5 10.4643 0.8533 8.7918 10.4643 12.1367
## Res. err. (variance) 0.3403 0.0014 0.3377 0.3403 0.3433
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 211.6698 5.8471 200.4594 211.5442 223.3859
## f1time_to_psa 612.0800 16.5756 582.2212 611.1390 647.6566
## f2time_to_psa 104.4675 2.9179 98.3298 104.7137 109.4553
## Intercept:f1time_to_psa -354.8559 9.8550 -375.1419 -354.4858 -336.6491
## Intercept:f2time_to_psa 65.1665 3.7962 57.6439 65.1025 72.8412
## f1time_to_psa:f2time_to_psa -72.9334 6.3988 -85.3398 -72.9183 -60.5386
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -84645.17 -84638.45
##
## Deviance Information Criterion: 144162.6
## Widely applicable Bayesian information criterion: 145954.1
## Computation time: 29.01 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 1677 on PID: 71801plot(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(25)) # 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.4095 0.0311 1.3485 1.4095 1.4704
## time_to_psa -0.0077 0.0006 -0.0089 -0.0077 -0.0064
## grade_group2 -0.2555 0.0273 -0.3089 -0.2555 -0.2020
## grade_group3 -0.2881 0.0311 -0.3490 -0.2881 -0.2272
## grade_group4 -0.2412 0.0429 -0.3252 -0.2412 -0.1571
## grade_group5 0.3069 0.0445 0.2198 0.3069 0.3941
## mri1 0.1627 0.0317 0.1005 0.1627 0.2249
## bx_age70 0.2825 0.0313 0.2211 0.2825 0.3439
## bx_age6070 0.1372 0.0287 0.0809 0.1372 0.1936
## drfs_txrp -0.4516 0.0264 -0.5034 -0.4516 -0.3999
## drfs_txunk -0.0992 0.0273 -0.1527 -0.0992 -0.0458
## time_to_psa:grade_group2 -0.0034 0.0009 -0.0052 -0.0034 -0.0015
## time_to_psa:grade_group3 -0.0041 0.0011 -0.0062 -0.0041 -0.0020
## time_to_psa:grade_group4 0.0031 0.0015 0.0001 0.0031 0.0061
## time_to_psa:grade_group5 0.0102 0.0017 0.0069 0.0102 0.0134
## Res. err. (variance) 0.4452 0.0024 0.4403 0.4452 0.4498
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.4616 0.0124 0.4380 0.4615 0.4865
## time_to_psa 0.0006 0.0000 0.0006 0.0006 0.0006
## Intercept:time_to_psa -0.0009 0.0005 -0.0020 -0.0008 -0.0002
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -90577.47 -90573.51
##
## Deviance Information Criterion: 162124.8
## Widely applicable Bayesian information criterion: 163671.6
## Computation time: 13.86 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 1677 on PID: 71801plot(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
# use
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.6126 0.0307 1.5524 1.6126 1.6727
## time_to_psa -0.0190 0.0010 -0.0210 -0.0190 -0.0171
## f1time_to_psa 0.0000 0.0004 -0.0008 0.0000 0.0007
## grade_group2 -0.2041 0.0281 -0.2592 -0.2041 -0.1490
## grade_group3 -0.2073 0.0320 -0.2701 -0.2073 -0.1446
## grade_group4 -0.1665 0.0443 -0.2532 -0.1665 -0.0797
## grade_group5 0.3497 0.0459 0.2597 0.3497 0.4398
## mri1 0.1254 0.0307 0.0652 0.1254 0.1856
## bx_age70 0.2684 0.0302 0.2092 0.2684 0.3275
## bx_age6070 0.1216 0.0276 0.0675 0.1216 0.1758
## drfs_txrp -0.5045 0.0254 -0.5543 -0.5045 -0.4547
## drfs_txunk -0.1710 0.0264 -0.2227 -0.1710 -0.1193
## time_to_psa:grade_group2 -0.0125 0.0015 -0.0154 -0.0125 -0.0095
## time_to_psa:grade_group3 -0.0175 0.0017 -0.0209 -0.0175 -0.0142
## time_to_psa:grade_group4 -0.0170 0.0024 -0.0218 -0.0170 -0.0123
## time_to_psa:grade_group5 -0.0248 0.0027 -0.0301 -0.0248 -0.0195
## f1time_to_psa:grade_group2 0.0001 0.0006 -0.0010 0.0001 0.0012
## f1time_to_psa:grade_group3 0.0002 0.0006 -0.0010 0.0002 0.0015
## f1time_to_psa:grade_group4 0.0005 0.0009 -0.0013 0.0005 0.0022
## f1time_to_psa:grade_group5 0.0023 0.0009 0.0005 0.0023 0.0041
## Res. err. (variance) 0.3161 0.0018 0.3126 0.3161 0.3196
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.4921 0.0132 0.4669 0.4922 0.5192
## time_to_psa 0.0014 0.0000 0.0013 0.0014 0.0014
## f1time_to_psa 0.0002 0.0000 0.0002 0.0002 0.0002
## Intercept:time_to_psa -0.0093 0.0006 -0.0107 -0.0093 -0.0082
## Intercept:f1time_to_psa 0.0006 0.0003 0.0001 0.0005 0.0011
## time_to_psa:f1time_to_psa 0.0000 0.0000 -0.0001 0.0000 0.0000
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -101178.9 -101172.2
##
## Deviance Information Criterion: 139952.1
## Widely applicable Bayesian information criterion: 141615.5
## Computation time: 19.32 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 1677 on PID: 71801plot(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)
#### Cubic
# use
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.6122 0.0307 1.5520 1.6122 1.6723
## time_to_psa -0.0190 0.0010 -0.0210 -0.0190 -0.0171
## f1time_to_psa 0.0000 0.0004 -0.0008 0.0000 0.0007
## grade_group2 -0.2039 0.0281 -0.2590 -0.2039 -0.1487
## grade_group3 -0.2070 0.0321 -0.2699 -0.2070 -0.1442
## grade_group4 -0.1662 0.0443 -0.2530 -0.1662 -0.0793
## grade_group5 0.3502 0.0460 0.2600 0.3502 0.4403
## mri1 0.1250 0.0307 0.0648 0.1250 0.1852
## bx_age70 0.2688 0.0302 0.2096 0.2688 0.3279
## bx_age6070 0.1219 0.0276 0.0677 0.1219 0.1761
## drfs_txrp -0.5051 0.0254 -0.5549 -0.5051 -0.4553
## drfs_txunk -0.1701 0.0264 -0.2218 -0.1701 -0.1184
## time_to_psa:grade_group2 -0.0125 0.0015 -0.0154 -0.0125 -0.0095
## time_to_psa:grade_group3 -0.0175 0.0017 -0.0209 -0.0175 -0.0142
## time_to_psa:grade_group4 -0.0170 0.0024 -0.0218 -0.0170 -0.0123
## time_to_psa:grade_group5 -0.0248 0.0027 -0.0302 -0.0248 -0.0195
## f1time_to_psa:grade_group2 0.0001 0.0006 -0.0010 0.0001 0.0012
## f1time_to_psa:grade_group3 0.0002 0.0006 -0.0010 0.0002 0.0015
## f1time_to_psa:grade_group4 0.0005 0.0009 -0.0013 0.0005 0.0022
## f1time_to_psa:grade_group5 0.0023 0.0009 0.0005 0.0023 0.0041
## Res. err. (variance) 0.3161 0.0018 0.3125 0.3161 0.3196
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.4930 0.0128 0.4684 0.4927 0.5189
## time_to_psa 0.0014 0.0000 0.0013 0.0014 0.0014
## f1time_to_psa 0.0002 0.0000 0.0002 0.0002 0.0002
## Intercept:time_to_psa -0.0094 0.0006 -0.0107 -0.0093 -0.0082
## Intercept:f1time_to_psa 0.0002 0.0002 -0.0001 0.0002 0.0006
## time_to_psa:f1time_to_psa 0.0000 0.0000 0.0000 0.0000 0.0000
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -101179.4 -101172.7
##
## Deviance Information Criterion: 139951.8
## Widely applicable Bayesian information criterion: 141617.6
## Computation time: 19.43 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 1677 on PID: 71801plot(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
Splines
############################################## 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), ]
# Setup the number of used threads
inla.setOption(num.threads='8:1')
NSplines <- ns(fm_surlong_rp$time_to_psa, knots=c(25)) # 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 <- 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),
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)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 3.1723 0.2572 2.6683 3.1723 3.6763
## f1time_to_psa_L1 -5.6375 0.4170 -6.4548 -5.6375 -4.8201
## f2time_to_psa_L1 -5.2345 0.2268 -5.6791 -5.2345 -4.7899
## grade_group2_L1 -0.1976 0.3458 -0.8754 -0.1976 0.4802
## grade_group3_L1 -0.7161 0.3651 -1.4317 -0.7161 -0.0005
## grade_group4_L1 -0.6153 0.5048 -1.6047 -0.6153 0.3740
## grade_group5_L1 -0.3035 0.5785 -1.4372 -0.3035 0.8303
## mri1_L1 0.0147 0.0259 -0.0360 0.0147 0.0655
## bx_age70_L1 0.1052 0.0310 0.0444 0.1052 0.1660
## bx_age6070_L1 0.0903 0.0235 0.0442 0.0903 0.1363
## f1time_to_psa:grade_group2_L1 -0.2643 0.5588 -1.3595 -0.2643 0.8309
## f1time_to_psa:grade_group3_L1 0.6955 0.5859 -0.4528 0.6955 1.8439
## f1time_to_psa:grade_group4_L1 0.8197 0.7973 -0.7430 0.8197 2.3825
## f1time_to_psa:grade_group5_L1 0.8980 0.9024 -0.8708 0.8980 2.6668
## f2time_to_psa:grade_group2_L1 -0.4232 0.3145 -1.0395 -0.4232 0.1932
## f2time_to_psa:grade_group3_L1 -0.1907 0.3419 -0.8607 -0.1907 0.4794
## f2time_to_psa:grade_group4_L1 0.5876 0.5049 -0.4019 0.5876 1.5771
## f2time_to_psa:grade_group5_L1 2.0700 0.5954 0.9031 2.0700 3.2370
## Res. err. (variance) 0.4761 0.0047 0.4667 0.4762 0.4850
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant
## Intercept_L1 9.6634 0.9925 7.8592 9.6426
## f1time_to_psa_L1 19.2094 2.1657 15.2558 19.0958
## f2time_to_psa_L1 14.6321 1.2237 12.4279 14.5778
## Intercept_L1:f1time_to_psa_L1 -13.5350 1.4596 -16.4976 -13.4714
## Intercept_L1:f2time_to_psa_L1 11.8154 0.9036 10.2010 11.8059
## f1time_to_psa_L1:f2time_to_psa_L1 -16.6036 1.3015 -19.3044 -16.5808
## 0.975quant
## Intercept_L1 11.7842
## f1time_to_psa_L1 23.7223
## f2time_to_psa_L1 17.1299
## Intercept_L1:f1time_to_psa_L1 -10.8899
## Intercept_L1:f2time_to_psa_L1 13.6813
## f1time_to_psa_L1:f2time_to_psa_L1 -14.2454
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0804 0.0786 0.0103 0.0564 0.2988
## grade_group2_S1 0.4829 0.0604 0.3645 0.4829 0.6013
## grade_group3_S1 0.6297 0.0672 0.4980 0.6297 0.7614
## grade_group4_S1 0.7777 0.0924 0.5966 0.7777 0.9588
## grade_group5_S1 0.9599 0.1088 0.7467 0.9599 1.1731
## mri1_S1 0.3603 0.0584 0.2458 0.3603 0.4748
## bx_age70_S1 0.2061 0.0701 0.0686 0.2061 0.3435
## bx_age6070_S1 0.0952 0.0523 -0.0073 0.0952 0.1978
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.0636 0.0543 0.0145 0.0469 0.2146
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0827 0.0795 0.0144 0.0578 0.3043
## grade_group2_S2 0.3383 0.0717 0.1978 0.3383 0.4787
## grade_group3_S2 0.5594 0.0777 0.4071 0.5594 0.7117
## grade_group4_S2 0.8885 0.1046 0.6834 0.8885 1.0935
## grade_group5_S2 1.4280 0.1222 1.1886 1.4280 1.6675
## mri1_S2 0.5626 0.0713 0.4229 0.5626 0.7023
## bx_age70_S2 0.4058 0.0828 0.2436 0.4058 0.5681
## bx_age6070_S2 0.1432 0.0627 0.0203 0.1432 0.2660
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -1.5137 0.0834 -1.6785 -1.5136 -1.3499
## CV_L1_S2 -0.8018 0.1601 -1.1753 -0.7802 -0.5698
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 1.7344 0.5997 0.6553 1.7043 3.0059
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -70292.20 -70279.96
##
## Deviance Information Criterion: -4860.322
## Widely applicable Bayesian information criterion: -7444.014
## Computation time: 59.65 secondsLinear
# 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), ]
# Setup the number of used threads
inla.setOption(num.threads='8:1')
# Final joint model, combining 2 survival modals and the longitudinal one
JM2_RP2 <- 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),
dataLong=fm_surlong_rp
)
summary(JM2_RP2)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 1.0923 0.0296 1.0343 1.0923 1.1504
## time_to_psa_L1 -0.0172 0.0017 -0.0206 -0.0172 -0.0139
## grade_group2_L1 -0.2780 0.0311 -0.3389 -0.2780 -0.2170
## grade_group3_L1 -0.3135 0.0340 -0.3802 -0.3135 -0.2468
## grade_group4_L1 -0.3184 0.0477 -0.4120 -0.3184 -0.2249
## grade_group5_L1 -0.2493 0.0567 -0.3604 -0.2493 -0.1382
## mri1_L1 0.1140 0.0304 0.0543 0.1140 0.1737
## bx_age70_L1 0.1074 0.0366 0.0358 0.1074 0.1790
## bx_age6070_L1 0.1028 0.0280 0.0479 0.1028 0.1577
## time_to_psa:grade_group2_L1 -0.0008 0.0023 -0.0053 -0.0008 0.0037
## time_to_psa:grade_group3_L1 0.0001 0.0025 -0.0048 0.0001 0.0050
## time_to_psa:grade_group4_L1 0.0040 0.0036 -0.0030 0.0040 0.0110
## time_to_psa:grade_group5_L1 0.0070 0.0043 -0.0013 0.0070 0.0154
## Res. err. (variance) 0.4881 0.0048 0.4790 0.4880 0.4979
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 0.1371 0.0101 0.1184 0.1367 0.1575
## time_to_psa_L1 0.0010 0.0000 0.0009 0.0010 0.0011
## Intercept_L1:time_to_psa_L1 -0.0025 0.0004 -0.0033 -0.0025 -0.0019
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0633 0.0523 0.0106 0.0483 0.2067
## grade_group2_S1 0.5056 0.0621 0.3839 0.5056 0.6273
## grade_group3_S1 0.6343 0.0690 0.4991 0.6343 0.7695
## grade_group4_S1 0.8118 0.0945 0.6265 0.8118 0.9971
## grade_group5_S1 1.0067 0.1114 0.7884 1.0067 1.2250
## mri1_S1 0.3442 0.0602 0.2261 0.3442 0.4622
## bx_age70_S1 0.1568 0.0722 0.0153 0.1568 0.2984
## bx_age6070_S1 0.0807 0.0543 -0.0257 0.0807 0.1872
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.0909 0.0079 0.0771 0.0904 0.108
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.3235 0.3735 0.0406 0.2034 1.3617
## grade_group2_S2 1.2870 0.2251 0.8458 1.2870 1.7282
## grade_group3_S2 1.8993 0.2473 1.4145 1.8993 2.3840
## grade_group4_S2 2.8986 0.3479 2.2167 2.8986 3.5805
## grade_group5_S2 4.4124 0.4136 3.6018 4.4124 5.2230
## mri1_S2 1.5610 0.2277 1.1147 1.5610 2.0074
## bx_age70_S2 1.0287 0.2716 0.4965 1.0287 1.5610
## bx_age6070_S2 0.3093 0.2021 -0.0869 0.3093 0.7055
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -1.4389 0.0956 -1.6192 -1.4414 -1.2432
## CV_L1_S2 -1.0454 0.2090 -1.4416 -1.0504 -0.6189
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 9.8593 0.4227 8.9446 9.8871 10.5943
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -71552.14 -71542.67
##
## Deviance Information Criterion: -6404.239
## Widely applicable Bayesian information criterion: -8402.502
## Computation time: 68.75 secondsMore Combinations
# 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), ]
# Setup the number of used threads
inla.setOption(num.threads='8:1')
f1 <- function(x) x*x
# Final joint model, combining 2 survival modals and the longitudinal one
JM2_RP3 <- 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),
dataLong=fm_surlong_rp
)
summary(JM2_RP3)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 1.3610 0.0331 1.2962 1.3610 1.4259
## time_to_psa_L1 -0.0370 0.0025 -0.0419 -0.0370 -0.0321
## f1time_to_psa_L1 0.0000 0.0014 -0.0028 0.0000 0.0028
## grade_group2_L1 -0.2940 0.0358 -0.3642 -0.2940 -0.2239
## grade_group3_L1 -0.3176 0.0392 -0.3944 -0.3176 -0.2407
## grade_group4_L1 -0.3772 0.0552 -0.4853 -0.3772 -0.2690
## grade_group5_L1 -0.3015 0.0660 -0.4309 -0.3015 -0.1721
## mri1_L1 0.0634 0.0327 -0.0007 0.0634 0.1275
## bx_age70_L1 0.0750 0.0394 -0.0021 0.0750 0.1522
## bx_age6070_L1 0.0706 0.0298 0.0122 0.0706 0.1291
## time_to_psa:grade_group2_L1 -0.0084 0.0034 -0.0150 -0.0084 -0.0018
## time_to_psa:grade_group3_L1 -0.0057 0.0037 -0.0129 -0.0057 0.0015
## time_to_psa:grade_group4_L1 -0.0021 0.0052 -0.0123 -0.0021 0.0082
## time_to_psa:grade_group5_L1 0.0008 0.0063 -0.0116 0.0008 0.0132
## f1time_to_psa:grade_group2_L1 0.0003 0.0019 -0.0034 0.0003 0.0041
## f1time_to_psa:grade_group3_L1 0.0003 0.0021 -0.0037 0.0003 0.0044
## f1time_to_psa:grade_group4_L1 0.0005 0.0030 -0.0053 0.0005 0.0063
## f1time_to_psa:grade_group5_L1 0.0005 0.0035 -0.0064 0.0005 0.0075
## Res. err. (variance) 0.3667 0.0037 0.3595 0.3666 0.3741
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 0.1893 0.0103 0.1710 0.1888 0.2105
## time_to_psa_L1 0.0020 0.0001 0.0018 0.0020 0.0022
## f1time_to_psa_L1 0.0007 0.0000 0.0007 0.0007 0.0008
## Intercept_L1:time_to_psa_L1 -0.0080 0.0006 -0.0092 -0.0080 -0.0069
## Intercept_L1:f1time_to_psa_L1 0.0001 0.0001 -0.0002 0.0001 0.0003
## time_to_psa_L1:f1time_to_psa_L1 0.0000 0.0000 -0.0001 0.0000 0.0000
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0346 0.0137 0.0147 0.0323 0.0681
## grade_group2_S1 0.4728 0.0629 0.3495 0.4728 0.5960
## grade_group3_S1 0.6244 0.0699 0.4873 0.6244 0.7615
## grade_group4_S1 0.7183 0.0960 0.5302 0.7183 0.9064
## grade_group5_S1 0.9380 0.1131 0.7163 0.9380 1.1596
## mri1_S1 0.3080 0.0613 0.1879 0.3080 0.4281
## bx_age70_S1 0.1224 0.0734 -0.0214 0.1224 0.2661
## bx_age6070_S1 0.0552 0.0549 -0.0523 0.0552 0.1628
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.104 0.011 0.0842 0.1034 0.1273
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.2751 0.1175 0.1156 0.2516 0.5718
## grade_group2_S2 1.5128 0.2191 1.0835 1.5128 1.9422
## grade_group3_S2 2.1041 0.2407 1.6323 2.1041 2.5758
## grade_group4_S2 2.9376 0.3386 2.2739 2.9376 3.6013
## grade_group5_S2 4.3968 0.4026 3.6077 4.3968 5.1860
## mri1_S2 1.7228 0.2219 1.2879 1.7228 2.1577
## bx_age70_S2 0.9171 0.2646 0.3984 0.9171 1.4357
## bx_age6070_S2 0.2365 0.1969 -0.1495 0.2365 0.6225
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -1.4579 0.0793 -1.6137 -1.4579 -1.3016
## CV_L1_S2 0.4948 0.1357 0.2205 0.4971 0.7547
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 9.0791 0.3855 8.3096 9.0827 9.8277
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -75860.42 -75848.19
##
## Deviance Information Criterion: -12025.31
## Widely applicable Bayesian information criterion: -14092.12
## Computation time: 128.19 secondsEffect post RT
Assessing the effect of time to treatment on risk of relapse for patients having RadioTherapy as there first Currative Treatment
Interpretation
cat('\nAFTER RP - splines: \n')##
## AFTER RP - splines:int_res <- interp(JM2_RP)
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.76 [0.18,0.99] 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 1.31 [1.01,5.43] 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_RP2)
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.1 [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 10.49 [1.2,1010.85] 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_RP3)
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.1 [0,0.85] 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 10.22 [1.18,898.21] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."#cat('\n\nAFTER RT\n')
#int_res <- interp(JM2_RT)
#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)