joint_models_report_part2
Abderrahim Oussama Batouche
2024-02-13
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)Importing data
surv <- read.csv("processed/surv.csv")
long <- read.csv("processed/long.csv")
surv_long_ <- read.csv("processed/surv_long_.csv")Some updates:
surv$mri <- as.factor(surv$mri)
surv$bx_age <- as.factor(surv$bx_age)
surv$grade_group <- as.factor(surv$grade_group)
surv$tx1_type[surv$tx1_type %in% c('chemo','horm')] <- 'CTHT'
surv$tx1_type <- factor(surv$tx1_type, levels = c('radio', 'rp', 'CTHT', 'tfs-unk'))
#surv$tx1_type <- as.factor(surv$tx1_type)
surv$after_tx <- as.factor(surv$after_tx)
surv$tx1_type <- relevel(surv$tx1_type, ref = "radio")
surv$after_tx <- relevel(surv$after_tx, ref = "rt")Modeling
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.M-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(surv))
summary(M1)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 1.4285 0.436 0.7681 1.3623 2.4719
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -25557.82 -25557.82
##
## Deviance Information Criterion: 51026.39
## Widely applicable Bayesian information criterion: 51026.79
## Computation time: 2.28 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+tx1_type),
basRisk = c("rw2"),
dataSurv = list(surv))
summary(M1_cov)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.9142 0.3158 0.4510 0.8618 1.6828
## bx_age70 0.2229 0.0312 0.1617 0.2229 0.2840
## bx_age6070 0.1336 0.0297 0.0753 0.1336 0.1919
## mri1 0.5044 0.0296 0.4464 0.5044 0.5624
## grade_group2 0.6476 0.0287 0.5914 0.6476 0.7038
## grade_group3 0.8895 0.0303 0.8301 0.8895 0.9489
## grade_group4 0.7855 0.0382 0.7105 0.7855 0.8605
## grade_group5 1.1679 0.0353 1.0985 1.1679 1.2372
## tx1_typerp 0.5423 0.0261 0.4912 0.5423 0.5935
## tx1_typeCTHT 0.1041 0.0257 0.0538 0.1041 0.1544
## tx1_typetfsunk -14.9541 2.8728 -20.5874 -14.9541 -9.3208
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -18546.38 -18546.38
##
## Deviance Information Criterion: 36923.4
## Widely applicable Bayesian information criterion: 36926.76
## Computation time: 2.96 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(surv))
summary(M2)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.1443 0.0793 0.0461 0.1259 0.3513
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -13683.63 -13683.63
##
## Deviance Information Criterion: 27311.55
## Widely applicable Bayesian information criterion: 27311.41
## Computation time: 2.28 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+tx1_type+tfs_time+after_tx),
basRisk = c("rw2"),
dataSurv = list(surv))
summary(M2_cov)##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance) 0.1437 0.0794 0.0458 0.1253 0.3513
## bx_age70 -0.0148 0.0575 -0.1275 -0.0148 0.0978
## bx_age6070 -0.0175 0.0527 -0.1207 -0.0175 0.0858
## mri1 0.2475 0.0641 0.1218 0.2475 0.3733
## grade_group2 0.0281 0.0592 -0.0879 0.0281 0.1442
## grade_group3 0.5772 0.0580 0.4634 0.5772 0.6909
## grade_group4 0.9350 0.0667 0.8043 0.9350 1.0658
## grade_group5 1.4646 0.0627 1.3416 1.4646 1.5875
## tx1_typerp -0.7584 0.3389 -1.4230 -0.7584 -0.0938
## tx1_typeCTHT 0.0028 0.0523 -0.0999 0.0028 0.1054
## tx1_typetfsunk -10.5665 6.0405 -22.4113 -10.5665 1.2783
## tfs_time -0.0057 0.0008 -0.0073 -0.0057 -0.0042
## after_txbcrunk -15.0724 3.9290 -22.7768 -15.0724 -7.3679
## after_txrp 1.0025 0.3360 0.3436 1.0025 1.6614
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -12447 -12447
##
## Deviance Information Criterion: 24738.26
## Widely applicable Bayesian information criterion: 24744.82
## Computation time: 3.78 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(surv))
summary(JM1)##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 3.0008 0.5424 2.0931 2.9462 4.2205
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.3743 0.0301 0.3204 0.3723 0.4388
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.041 0.0073 0.028 0.0406 0.0565
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 0.5208 0.1182 0.3058 0.5154 0.7695
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -39927.88 -39927.88
##
## Deviance Information Criterion: 79373.3
## Widely applicable Bayesian information criterion: 79721.81
## Computation time: 13.38 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+tx1_type+(1|id),
inla.surv(time = drfs_time, event = drfs_status) ~ bx_age+mri+grade_group+tx1_type+tfs_time+after_tx), id="id",
basRisk=c("rw2", "rw2"), assocSurv=TRUE, NbasRisk = 30,
dataSurv = list(surv))## 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_cov)##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 1.9428 0.4866 1.1443 1.8917 3.0485
## bx_age70_S1 0.1859 0.0327 0.1218 0.1859 0.2499
## bx_age6070_S1 0.1123 0.0311 0.0514 0.1123 0.1732
## mri1_S1 0.4302 0.0308 0.3699 0.4302 0.4905
## grade_group2_S1 0.5650 0.0302 0.5058 0.5650 0.6242
## grade_group3_S1 0.7690 0.0320 0.7064 0.7690 0.8318
## grade_group4_S1 0.6896 0.0401 0.6110 0.6896 0.7683
## grade_group5_S1 1.0828 0.0377 1.0089 1.0828 1.1568
## tx1_typerp_S1 0.4188 0.0271 0.3657 0.4188 0.4719
## tx1_typeCTHT_S1 0.1892 0.0274 0.1355 0.1891 0.2429
## tx1_typetfsunk_S1 -11.1001 2.8748 -16.7372 -11.1001 -5.4630
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.0533 0.0129 0.0314 0.0523 0.0817
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.2057 0.0830 0.0808 0.1939 0.4025
## bx_age70_S2 -0.0002 0.0628 -0.1233 -0.0002 0.1230
## bx_age6070_S2 -0.0087 0.0578 -0.1221 -0.0087 0.1046
## mri1_S2 0.2846 0.0686 0.1501 0.2846 0.4192
## grade_group2_S2 0.0933 0.0637 -0.0316 0.0933 0.2183
## grade_group3_S2 0.6984 0.0636 0.5738 0.6984 0.8231
## grade_group4_S2 1.1112 0.0746 0.9651 1.1112 1.2576
## grade_group5_S2 1.7332 0.0724 1.5914 1.7332 1.8754
## tx1_typerp_S2 -0.8311 0.3872 -1.5903 -0.8311 -0.0719
## tx1_typeCTHT_S2 -0.0895 0.0575 -0.2024 -0.0895 0.0232
## tx1_typetfsunk_S2 -2.9149 5.8964 -14.4770 -2.9149 8.6474
## tfs_time_S2 -0.0012 0.0009 -0.0030 -0.0012 0.0006
## after_txbcrunk_S2 -7.1790 4.0801 -15.1795 -7.1790 0.8217
## after_txrp_S2 1.1507 0.3842 0.3975 1.1507 1.9042
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 2.8234 0.3543 2.1486 2.8159 3.543
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -34035.75 -34035.75
##
## Deviance Information Criterion: 68841.93
## Widely applicable Bayesian information criterion: 7.096363e+12
## Computation time: 25.42 secondsplot(JM1_cov)$Baseline+scale_y_log10()
# Survival curves
onePatient <- surv[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 4 on PID: 9394## 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 4# 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 4 on PID: 9394## 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 4lines(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)
#################################### SAMPLING + CHECK THE DATA
sampling <- function(data, sample_size=0.3){
# Uncomment when needed, Sampling
unique_ids <- unique(surv_long_$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
surv_long_data <- surv_long_[surv_long_$id %in% sampled_ids, ]
# Reassigning ID
# -> Calculate run-length encoding of the original id column
id_rle <- rle(surv_long_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 surv_long_data
surv_long_data$id <- new_id
surv_long_data$tx1_type[surv_long_data$tx1_type %in% c('chemo','horm')] <- 'CTHT'
# factors
surv_long_data$mri <- as.factor(surv_long_data$mri)
surv_long_data$grade_group <- as.factor(surv_long_data$grade_group)
surv_long_data$bx_age <- as.factor(surv_long_data$bx_age)
surv_long_data$tfs_status <- as.factor(surv_long_data$tfs_status)
surv_long_data$drfs_status <- as.factor(surv_long_data$drfs_status)
surv_long_data$after_tx <- as.factor(surv_long_data$after_tx)
surv_long_data$tx1_type <- as.factor(surv_long_data$tx1_type)
# order factors
surv_long_data$tx1_type <- factor(surv_long_data$tx1_type, levels = c('radio', 'rp', 'CTHT', 'tfs-unk'))
surv_long_data$after_tx <- factor(surv_long_data$after_tx, levels = c('rt', 'rp', 'bcr-unk'))
# set reference
surv_long_data$tx1_type <- relevel(surv_long_data$tx1_type, ref = "radio")
surv_long_data$after_tx <- relevel(surv_long_data$after_tx, ref = "rt")
# Remove Nas
surv_long_data <- surv_long_data[!is.na(surv_long_data$time_to_psa),]
return(surv_long_data)
}
#Sampling
surv_long_data <- sampling(surv_long_, 0.4)
# Surv data to be used in the latter join model
surv_nolong_data <- surv_long_data[!duplicated(surv_long_data$id), ]
# Rename my dfs
fsurvlong <- surv_long_data
fsurv <- surv_nolong_data
summary(fsurvlong)## id bx_age mri grade_group tx1_type
## Min. : 1 <60 :13076 0:68070 1:26196 radio :29501
## 1st Qu.:1160 >70 :29032 1: 9699 2:20483 rp :18204
## Median :2352 60-70:35661 3:15869 CTHT :17356
## Mean :2348 4: 7267 tfs-unk:12708
## 3rd Qu.:3519 5: 7954
## Max. :4701
## tx_date tfs_time tfs_status after_tx
## Length:77769 Min. : 0.00 0:12708 rt :43197
## Class :character 1st Qu.: 2.00 1:65061 rp :18399
## Mode :character Median : 5.00 bcr-unk:16173
## Mean : 29.23
## 3rd Qu.: 36.00
## Max. :256.00
## drfs_time drfs_status psa_results psa_date
## Min. : 1.00 0:55046 Min. : 0.00 Length:77769
## 1st Qu.: 40.00 1:22723 1st Qu.: 0.08 Class :character
## Median : 70.00 Median : 1.21 Mode :character
## Mean : 79.51 Mean : 13.95
## 3rd Qu.:113.00 3rd Qu.: 6.00
## Max. :256.00 Max. :76800.00
## time_to_psa log_psa
## Min. :-242.00 Min. : 0.00000
## 1st Qu.: 1.00 1st Qu.: 0.07696
## Median : 24.00 Median : 0.79299
## Mean : 30.22 Mean : 1.12905
## 3rd Qu.: 59.00 3rd Qu.: 1.94591
## Max. : 265.00 Max. :11.24897library(splines)
# We applied log transformation to PSA with shift=1
# Longitudinal Data should be sorted for each patient.
fsurvlong <- fsurvlong[with(fsurvlong, order(id, time_to_psa)), ]
# 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=fsurvlong
)
summary(M3)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 1.3695 0.0123 1.3454 1.3695 1.3936
## time_to_psa -0.0094 0.0004 -0.0101 -0.0094 -0.0086
## Res. err. (variance) 0.5212 0.0028 0.5157 0.5212 0.5267
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 0.6314 0.0152 0.6020 0.6313 0.6618
## time_to_psa 0.0007 0.0000 0.0006 0.0007 0.0007
## Intercept:time_to_psa 0.0001 0.0004 -0.0007 0.0001 0.0008
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -98899.51 -98895.54
##
## Deviance Information Criterion: 178094
## Widely applicable Bayesian information criterion: 179633.6
## Computation time: 10.46 seconds# patientID = 25 # pick one
patient_counts <- table(fsurvlong$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 <- fsurvlong[fsurvlong$id==patientID,]
P1 <- predict(M3, ND, id="id", horizon=max(fsurvlong$time_to_psa))$PredL # Make Linear prediction## Start sampling## Sampling done.## Computing longitudinal predictions for individual 1612 on PID: 9394# 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) # to HERE , pch=size of the dots
# use splines?
NSplines <- ns(fsurvlong$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(fsurvlong$time_to_psa), ylim=c(-1,1))
curve(f2, xlim=range(fsurvlong$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+after_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=fsurvlong
)## 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.1724 0.5575 8.0797 9.1724 10.2650
## f1time_to_psa -15.1937 0.9491 -17.0538 -15.1937 -13.3335
## f2time_to_psa -3.4542 0.2843 -4.0115 -3.4542 -2.8970
## grade_group2 3.7586 0.8222 2.1470 3.7586 5.3701
## grade_group3 8.3358 0.8989 6.5740 8.3358 10.0977
## grade_group4 6.1956 1.1775 3.8877 6.1956 8.5036
## grade_group5 11.3536 1.1483 9.1030 11.3536 13.6042
## mri1 0.0629 0.0293 0.0054 0.0629 0.1204
## bx_age70 0.2535 0.0308 0.1932 0.2535 0.3139
## bx_age6070 0.1582 0.0294 0.1007 0.1582 0.2158
## after_txrp -0.4781 0.0269 -0.5308 -0.4781 -0.4254
## after_txbcrunk -0.1055 0.0260 -0.1564 -0.1055 -0.0546
## f1time_to_psa:grade_group2 -7.7064 1.3983 -10.4471 -7.7064 -4.9658
## f1time_to_psa:grade_group3 -15.2830 1.5280 -18.2779 -15.2830 -12.2882
## f1time_to_psa:grade_group4 -10.0502 1.9985 -13.9672 -10.0502 -6.1333
## f1time_to_psa:grade_group5 -16.1900 1.9322 -19.9770 -16.1900 -12.4030
## f2time_to_psa:grade_group2 -1.1106 0.4288 -1.9509 -1.1106 -0.2702
## f2time_to_psa:grade_group3 0.8112 0.4685 -0.1071 0.8112 1.7294
## f2time_to_psa:grade_group4 3.1437 0.6368 1.8955 3.1437 4.3918
## f2time_to_psa:grade_group5 8.9940 0.6532 7.7137 8.9940 10.2743
## Res. err. (variance) 0.3957 0.0020 0.3923 0.3955 0.4001
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 371.3493 8.6521 355.1144 371.1286 389.0664
## f1time_to_psa 1091.2794 20.4212 1053.8269 1090.6857 1134.8672
## f2time_to_psa 88.5020 1.4730 86.0277 88.3538 91.5593
## Intercept:f1time_to_psa -631.9603 13.3805 -659.4835 -631.5585 -607.1706
## Intercept:f2time_to_psa 77.4405 5.5729 66.9838 77.1757 89.1395
## f1time_to_psa:f2time_to_psa -99.4742 9.5514 -119.8387 -99.0524 -81.6505
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -93005.05 -92998.33
##
## Deviance Information Criterion: 159477.4
## Widely applicable Bayesian information criterion: 161035.1
## Computation time: 38.4 seconds# Again, pick one patient (here the same)
ND2 <- fsurvlong[fsurvlong$id==patientID,] # observed vs. fitted for a couple individuals.
P2 <- predict(M3_cov, ND2, id="id", horizon=max(fsurvlong$time_to_psa))$PredL## Start sampling## Sampling done.## Computing longitudinal predictions for individual 1612 on PID: 9394plot(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
#Sampling
sample_size <- 0.5
surv_long_data <- sampling(surv_long_, sample_size)
# Surv data to be used in the latter join model
surv_nolong_data <- surv_long_data[!duplicated(surv_long_data$id), ]
# Rename my dfs
fsurvlong <- surv_long_data
fsurv <- surv_nolong_data
# Setup the number of used threads
inla.setOption(num.threads='8:1')
# 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(fsurv[fsurv$after_tx=='rp',]),
dataLong=fsurvlong
)
summary(JM2_RP)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 1.9945 0.0923 1.8136 1.9945 2.1755
## f1time_to_psa_L1 -2.1169 0.1496 -2.4100 -2.1169 -1.8237
## f2time_to_psa_L1 -3.3618 0.1319 -3.6202 -3.3618 -3.1034
## grade_group2_L1 0.5886 0.1167 0.3600 0.5886 0.8173
## grade_group3_L1 0.5955 0.1215 0.3573 0.5955 0.8337
## grade_group4_L1 0.3550 0.1439 0.0730 0.3550 0.6369
## grade_group5_L1 -0.2506 0.1464 -0.5376 -0.2506 0.0364
## mri1_L1 -0.9503 0.0494 -1.0470 -0.9503 -0.8535
## bx_age70_L1 -0.2221 0.0512 -0.3224 -0.2221 -0.1218
## bx_age6070_L1 -0.0697 0.0495 -0.1667 -0.0697 0.0273
## f1time_to_psa:grade_group2_L1 -1.9387 0.2083 -2.3469 -1.9387 -1.5306
## f1time_to_psa:grade_group3_L1 -1.7789 0.2158 -2.2019 -1.7789 -1.3559
## f1time_to_psa:grade_group4_L1 -0.4603 0.2575 -0.9651 -0.4603 0.0445
## f1time_to_psa:grade_group5_L1 1.6230 0.2543 1.1245 1.6230 2.1215
## f2time_to_psa:grade_group2_L1 -1.3820 0.2004 -1.7748 -1.3820 -0.9892
## f2time_to_psa:grade_group3_L1 -1.1631 0.2206 -1.5955 -1.1631 -0.7307
## f2time_to_psa:grade_group4_L1 0.1368 0.2900 -0.4316 0.1368 0.7053
## f2time_to_psa:grade_group5_L1 0.6179 0.2999 0.0301 0.6179 1.2057
## Res. err. (variance) 0.5670 0.0030 0.5613 0.5670 0.5730
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 0.0922 0.0213 0.0581 0.0892 0.1430
## f1time_to_psa_L1 4.4444 0.1381 4.1897 4.4398 4.7254
## f2time_to_psa_L1 21.1880 0.7285 19.8907 21.1359 22.7089
## Intercept_L1:f1time_to_psa_L1 -0.1081 0.0570 -0.2157 -0.1121 0.0060
## Intercept_L1:f2time_to_psa_L1 -0.0062 0.0996 -0.2098 -0.0032 0.1776
## f1time_to_psa_L1:f2time_to_psa_L1 0.4345 0.2228 0.0205 0.4288 0.8769
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0555 0.0337 0.0144 0.0477 0.1436
## grade_group2_S1 0.6624 0.0558 0.5530 0.6624 0.7718
## grade_group3_S1 0.7969 0.0630 0.6734 0.7969 0.9203
## grade_group4_S1 0.9265 0.0836 0.7626 0.9265 1.0905
## grade_group5_S1 1.0039 0.0986 0.8106 1.0039 1.1972
## mri1_S1 0.2852 0.0529 0.1815 0.2852 0.3888
## bx_age70_S1 0.0437 0.0645 -0.0827 0.0437 0.1701
## bx_age6070_S1 0.0032 0.0499 -0.0947 0.0032 0.1011
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.0258 0.0111 0.0094 0.0241 0.0525
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0566 0.0257 0.0214 0.0516 0.1210
## grade_group2_S2 0.4158 0.0666 0.2852 0.4158 0.5464
## grade_group3_S2 0.6255 0.0735 0.4813 0.6255 0.7696
## grade_group4_S2 0.7841 0.0950 0.5980 0.7841 0.9703
## grade_group5_S2 1.0928 0.1104 0.8764 1.0928 1.3092
## mri1_S2 0.6213 0.0674 0.4893 0.6213 0.7533
## bx_age70_S2 0.2976 0.0776 0.1456 0.2976 0.4497
## bx_age6070_S2 0.0115 0.0607 -0.1074 0.0115 0.1305
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -0.8859 0.0586 -1.0104 -0.8828 -0.7809
## CV_L1_S2 0.3354 0.0905 0.1595 0.3346 0.5158
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 1.159 0.6735 -0.1754 1.162 2.4764
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -169485.6 -169473.4
##
## Deviance Information Criterion: 173728.2
## Widely applicable Bayesian information criterion: 172561.1
## Computation time: 143.79 secondsEffect post RT
Assessing the effect of time to treatment on risk of relapse for patients having RadioTherapy as there first Currative Treatment
fsurvlong_rt <- fsurvlong[fsurvlong$after_tx=='rt',]
# Reassigning ID
# -> Calculate run-length encoding of the original id column
id_rle <- rle(fsurvlong_rt$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 fsurvlong_rt
fsurvlong_rt$id <- new_id
fsurv_rt <- fsurvlong_rt[!duplicated(fsurvlong_rt$id), ]
# Final joint model, combining 2 survival modals and the longitudinal one
JM2_RT <- 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=fsurv_rt,
dataLong=fsurvlong_rt
)## 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_RT)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 11.0554 0.8704 9.3495 11.0554 12.7613
## f1time_to_psa_L1 -18.5539 1.4853 -21.4650 -18.5539 -15.6428
## f2time_to_psa_L1 -3.8064 0.4797 -4.7465 -3.8064 -2.8662
## grade_group2_L1 2.2260 1.2626 -0.2486 2.2260 4.7006
## grade_group3_L1 8.1609 1.3084 5.5964 8.1609 10.7253
## grade_group4_L1 8.2078 1.6502 4.9735 8.2078 11.4421
## grade_group5_L1 9.2173 1.5080 6.2617 9.2173 12.1728
## mri1_L1 -0.0253 0.0427 -0.1090 -0.0253 0.0584
## bx_age70_L1 0.0663 0.0468 -0.0254 0.0663 0.1580
## bx_age6070_L1 0.0731 0.0472 -0.0193 0.0731 0.1656
## f1time_to_psa:grade_group2_L1 -5.0578 2.1538 -9.2791 -5.0578 -0.8365
## f1time_to_psa:grade_group3_L1 -14.5008 2.2273 -18.8662 -14.5008 -10.1355
## f1time_to_psa:grade_group4_L1 -12.1650 2.8000 -17.6529 -12.1650 -6.6771
## f1time_to_psa:grade_group5_L1 -12.3382 2.5456 -17.3275 -12.3382 -7.3488
## f2time_to_psa:grade_group2_L1 -1.3800 0.7011 -2.7542 -1.3800 -0.0059
## f2time_to_psa:grade_group3_L1 1.3862 0.7341 -0.0527 1.3862 2.8250
## f2time_to_psa:grade_group4_L1 6.5915 0.9482 4.7330 6.5915 8.4499
## f2time_to_psa:grade_group5_L1 8.9553 0.8947 7.2018 8.9553 10.7089
## Res. err. (variance) 0.4566 0.0028 0.4509 0.4567 0.4621
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant
## Intercept_L1 569.7622 26.4613 520.6850 568.4971
## f1time_to_psa_L1 1689.8895 77.6847 1547.0745 1684.7344
## f2time_to_psa_L1 150.6147 5.9828 139.9116 150.0624
## Intercept_L1:f1time_to_psa_L1 -973.8521 45.2016 -1069.3380 -971.5855
## Intercept_L1:f2time_to_psa_L1 116.6054 9.1175 99.0093 116.4421
## f1time_to_psa_L1:f2time_to_psa_L1 -144.8517 15.2844 -175.2515 -144.6202
## 0.975quant
## Intercept_L1 625.4419
## f1time_to_psa_L1 1852.0858
## f2time_to_psa_L1 163.0072
## Intercept_L1:f1time_to_psa_L1 -889.8816
## Intercept_L1:f2time_to_psa_L1 135.0325
## f1time_to_psa_L1:f2time_to_psa_L1 -114.8852
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0828 0.0045 0.0743 0.0827 0.0919
## grade_group2_S1 0.7520 0.0595 0.6353 0.7520 0.8686
## grade_group3_S1 1.1036 0.0609 0.9842 1.1036 1.2230
## grade_group4_S1 1.1709 0.0757 1.0225 1.1709 1.3193
## grade_group5_S1 1.7286 0.0679 1.5955 1.7286 1.8616
## mri1_S1 0.5757 0.0597 0.4587 0.5757 0.6928
## bx_age70_S1 0.5436 0.0659 0.4145 0.5436 0.6727
## bx_age6070_S1 0.3656 0.0664 0.2355 0.3656 0.4957
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.7148 0.0137 0.6873 0.7151 0.7412
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0656 0.0067 0.0542 0.0650 0.0804
## grade_group2_S2 0.2659 0.0682 0.1321 0.2659 0.3997
## grade_group3_S2 0.6806 0.0697 0.5441 0.6806 0.8172
## grade_group4_S2 0.7807 0.0853 0.6134 0.7807 0.9479
## grade_group5_S2 1.2444 0.0766 1.0942 1.2444 1.3947
## mri1_S2 1.4247 0.0719 1.2838 1.4247 1.5656
## bx_age70_S2 0.4965 0.0745 0.3505 0.4965 0.6426
## bx_age6070_S2 0.2618 0.0744 0.1159 0.2618 0.4076
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -0.1651 0.0232 -0.2078 -0.1661 -0.1167
## CV_L1_S2 0.5386 0.0109 0.5161 0.5388 0.5591
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 0.9944 0.0447 0.8923 0.9995 1.0639
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -151462.5 -151450.2
##
## Deviance Information Criterion: 18572.42
## Widely applicable Bayesian information criterion: 15867.08
## Computation time: 200.64 secondsInterpretation
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))
}
cat('\nAFTER RP\n')##
## AFTER RPint_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.89 [0.51,1.01] 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.13 [0.99,1.99] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."cat('\n\nAFTER RT\n')##
##
## AFTER RTint_res <- interp(JM2_RT)
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.52 [0.14,0.95] 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.91 [1.05,7.04] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."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 differente risk associations
# CV, CV_CS, SRE
JM2_CVCS <- joint(formSurv=list(inla.surv(tfs_time, tfs_status) ~ grade_group+mri+bx_age + (1|id),
inla.surv(drfs_time, drfs_status) ~ grade_group+mri+bx_age+after_tx),
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_CS", "CV_CS"), # what is cv : Current value = risk event at t depends on longit.marker(psa) at t
control=list(int.strategy="eb"),
dataSurv=list(fsurv),
dataLong=fsurvlong
)## 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_CVCS)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 9.9303 0.5002 8.9499 9.9303 10.9107
## f1time_to_psa_L1 -16.5933 0.8455 -18.2504 -16.5933 -14.9363
## f2time_to_psa_L1 -2.7804 0.2983 -3.3651 -2.7804 -2.1958
## grade_group2_L1 3.5162 0.7407 2.0645 3.5162 4.9679
## grade_group3_L1 7.0230 0.8007 5.4536 7.0230 8.5924
## grade_group4_L1 6.9018 1.0342 4.8749 6.9018 8.9287
## grade_group5_L1 7.1755 1.0107 5.1946 7.1755 9.1565
## mri1_L1 0.0500 0.0278 -0.0045 0.0500 0.1044
## bx_age70_L1 0.3144 0.0294 0.2567 0.3144 0.3721
## bx_age6070_L1 0.1471 0.0286 0.0911 0.1471 0.2031
## f1time_to_psa:grade_group2_L1 -7.3263 1.2501 -9.7764 -7.3263 -4.8762
## f1time_to_psa:grade_group3_L1 -13.0958 1.3494 -15.7405 -13.0958 -10.4510
## f1time_to_psa:grade_group4_L1 -11.3416 1.7399 -14.7517 -11.3416 -7.9315
## f1time_to_psa:grade_group5_L1 -9.6368 1.6886 -12.9464 -9.6368 -6.3273
## f2time_to_psa:grade_group2_L1 -1.1762 0.4492 -2.0566 -1.1762 -0.2957
## f2time_to_psa:grade_group3_L1 0.0142 0.4891 -0.9445 0.0142 0.9728
## f2time_to_psa:grade_group4_L1 3.1124 0.6424 1.8532 3.1124 4.3715
## f2time_to_psa:grade_group5_L1 6.4065 0.6555 5.1217 6.4065 7.6914
## Res. err. (variance) 0.4015 0.0021 0.3972 0.4016 0.4055
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant
## Intercept_L1 359.2457 13.2238 335.4928 358.3589
## f1time_to_psa_L1 1044.8932 36.8220 980.7592 1041.6253
## f2time_to_psa_L1 125.0094 4.0971 117.1070 124.9376
## Intercept_L1:f1time_to_psa_L1 -606.4905 22.0310 -654.5778 -604.7612
## Intercept_L1:f2time_to_psa_L1 89.7997 5.4363 79.2790 89.6607
## f1time_to_psa_L1:f2time_to_psa_L1 -107.0782 8.9311 -125.0449 -106.9978
## 0.975quant
## Intercept_L1 388.2346
## f1time_to_psa_L1 1125.8172
## f2time_to_psa_L1 133.1679
## Intercept_L1:f1time_to_psa_L1 -567.4804
## Intercept_L1:f2time_to_psa_L1 100.8109
## f1time_to_psa_L1:f2time_to_psa_L1 -90.0251
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0586 0.0052 0.0501 0.0580 0.0704
## grade_group2_S1 1.0523 0.0539 0.9466 1.0523 1.1579
## grade_group3_S1 1.4455 0.0579 1.3319 1.4455 1.5590
## grade_group4_S1 1.4812 0.0741 1.3360 1.4812 1.6265
## grade_group5_S1 2.1544 0.0694 2.0184 2.1544 2.2904
## mri1_S1 0.6676 0.0560 0.5578 0.6676 0.7773
## bx_age70_S1 0.5589 0.0591 0.4430 0.5589 0.6748
## bx_age6070_S1 0.3822 0.0573 0.2699 0.3822 0.4945
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 1.729 0.0281 1.6822 1.726 1.7911
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0843 0.0031 0.0785 0.0842 0.0906
## grade_group2_S2 1.3773 0.0958 1.1895 1.3773 1.5650
## grade_group3_S2 2.1363 0.1032 1.9341 2.1363 2.3385
## grade_group4_S2 2.1768 0.1316 1.9189 2.1768 2.4347
## grade_group5_S2 3.0010 0.1238 2.7583 3.0010 3.2437
## mri1_S2 2.0109 0.1008 1.8133 2.0109 2.2086
## bx_age70_S2 0.9002 0.1056 0.6932 0.9002 1.1072
## bx_age6070_S2 0.4965 0.1012 0.2981 0.4965 0.6949
## after_txrp_S2 1.1393 0.0600 1.0218 1.1393 1.2568
## after_txbcrunk_S2 3.7512 0.0600 3.6336 3.7512 3.8687
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 -0.1180 0.0159 -0.1486 -0.1182 -0.0858
## CS_L1_S1 -0.9334 0.0945 -1.1196 -0.9334 -0.7476
## CV_L1_S2 0.7596 0.0217 0.7122 0.7610 0.7968
## CS_L1_S2 0.2949 0.0936 0.1083 0.2956 0.4770
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 1.9096 0.0299 1.8436 1.9119 1.96
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -449889.1 -449875.1
##
## Deviance Information Criterion: -336425.9
## Widely applicable Bayesian information criterion: -354139.8
## Computation time: 823.29 secondsint_res <- interp(JM2_CVCS)
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.13 [0,0.86] 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 7.57 [1.15,334.67] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."JM2_SRE <- joint(formSurv=list(inla.surv(tfs_time, tfs_status) ~ grade_group+mri+bx_age + (1|id),
inla.surv(drfs_time, drfs_status) ~ 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("SRE", "SRE"), # what is cv : Current value = risk event at t depends on longit.marker(psa) at t
control=list(int.strategy="eb"),
dataSurv=list(fsurv),
dataLong=fsurvlong
)## 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_SRE)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 10.7521 0.4451 9.8798 10.7521 11.6245
## f1time_to_psa_L1 -17.8839 0.7797 -19.4120 -17.8839 -16.3558
## f2time_to_psa_L1 -2.2508 0.2345 -2.7103 -2.2508 -1.7912
## grade_group2_L1 3.2705 0.6030 2.0886 3.2705 4.4525
## grade_group3_L1 5.3621 0.6899 4.0099 5.3621 6.7144
## grade_group4_L1 4.8984 0.8906 3.1528 4.8984 6.6440
## grade_group5_L1 4.4416 0.9010 2.6756 4.4416 6.2075
## mri1_L1 0.0652 0.0277 0.0110 0.0652 0.1194
## bx_age70_L1 0.3315 0.0293 0.2741 0.3315 0.3888
## bx_age6070_L1 0.1630 0.0284 0.1074 0.1630 0.2187
## f1time_to_psa:grade_group2_L1 -6.8883 1.0805 -9.0060 -6.8883 -4.7706
## f1time_to_psa:grade_group3_L1 -10.5905 1.2166 -12.9751 -10.5905 -8.2059
## f1time_to_psa:grade_group4_L1 -8.3643 1.5690 -11.4393 -8.3643 -5.2892
## f1time_to_psa:grade_group5_L1 -5.7811 1.5614 -8.8414 -5.7811 -2.7207
## f2time_to_psa:grade_group2_L1 -1.2166 0.3004 -1.8054 -1.2166 -0.6278
## f2time_to_psa:grade_group3_L1 -1.1436 0.3554 -1.8402 -1.1436 -0.4470
## f2time_to_psa:grade_group4_L1 1.6228 0.4650 0.7114 1.6228 2.5342
## f2time_to_psa:grade_group5_L1 3.8734 0.5075 2.8788 3.8734 4.8680
## Res. err. (variance) 0.4050 0.0021 0.4008 0.4050 0.4090
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant
## Intercept_L1 372.4317 15.8792 343.2991 371.3711
## f1time_to_psa_L1 1086.1295 42.6344 1008.9172 1083.2561
## f2time_to_psa_L1 126.9316 5.0105 117.0003 127.0075
## Intercept_L1:f1time_to_psa_L1 -629.5014 25.9088 -685.9196 -627.8332
## Intercept_L1:f2time_to_psa_L1 90.7197 6.9261 77.0161 90.7235
## f1time_to_psa_L1:f2time_to_psa_L1 -106.9700 10.9185 -128.8973 -106.9050
## 0.975quant
## Intercept_L1 406.3460
## f1time_to_psa_L1 1180.2903
## f2time_to_psa_L1 136.6337
## Intercept_L1:f1time_to_psa_L1 -581.9810
## Intercept_L1:f2time_to_psa_L1 104.7369
## f1time_to_psa_L1:f2time_to_psa_L1 -85.2963
##
## Survival outcome (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0817 0.0119 0.0657 0.0792 0.1111
## grade_group2_S1 0.8996 0.0371 0.8269 0.8996 0.9723
## grade_group3_S1 1.2459 0.0398 1.1680 1.2459 1.3239
## grade_group4_S1 1.2484 0.0510 1.1483 1.2484 1.3484
## grade_group5_S1 1.7823 0.0478 1.6887 1.7823 1.8759
## mri1_S1 0.5094 0.0386 0.4338 0.5094 0.5851
## bx_age70_S1 0.3422 0.0408 0.2623 0.3422 0.4221
## bx_age6070_S1 0.2626 0.0395 0.1852 0.2626 0.3400
##
## Frailty term variance (S1)
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1 0.4885 0.0083 0.4761 0.4873 0.5075
##
## Survival outcome (S2)
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S2 0.0524 0.0045 0.0453 0.0518 0.0627
## grade_group2_S2 0.3712 0.0635 0.2467 0.3712 0.4957
## grade_group3_S2 0.8058 0.0668 0.6749 0.8058 0.9366
## grade_group4_S2 1.1572 0.0837 0.9932 1.1572 1.3213
## grade_group5_S2 1.9796 0.0773 1.8282 1.9796 2.1310
## mri1_S2 1.5237 0.0603 1.4055 1.5237 1.6420
## bx_age70_S2 0.6411 0.0619 0.5198 0.6411 0.7624
## bx_age6070_S2 0.3141 0.0596 0.1974 0.3141 0.4309
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## SRE_L1_S1 -0.0111 0.0117 -0.0319 -0.0117 0.0139
## SRE_L1_S2 0.5943 0.0126 0.5722 0.5936 0.6215
##
## Association survival - survival
## mean sd 0.025quant 0.5quant 0.975quant
## IDIntercept_S1_S2 1.6605 0.0167 1.6294 1.6599 1.6953
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -305512.4 -305500.1
##
## Deviance Information Criterion: -24368.16
## Widely applicable Bayesian information criterion: -31228.5
## Computation time: 253.72 secondsint_res <- interp(JM2_SRE)
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.54 [1.08,14.96] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual."