joint_models_report_part4
Abderrahim Oussama Batouche
2024-02-15
Intro:
This document demonstrate how to use the predict()
function from INLAJoint library. The aim is to create a
prediction model that:
Giving patient i
data(type of treatment, time to treatment, psa results) and prediction time (5 years) (horizon), it returns therecurrence-free probability(Risk).Or, giving patient i
dataandrecurrence risk, returns thetime to recurrence.
The process is as follow :
fit (using
joint()function ) a joint model that relates the risk of recurrence (accounting for different covariates) and the longitudinal PSA data.Once the model is fitted, we can use it to
predict()both the survival data (at N data points) and the longitudinal data (also at N data points). We can then use splines to interpolate and get the function that links (approximate) those predicted data points in order to obtain the desired prediction data point(at timetor with riskr). It is important to keep in mind that the predict function horizon cannot (will lead to false results) predict more then the seen data (max(data$drfs_time))
Libraries
library(readr)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(INLA)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: sp## This is INLA_24.02.09 built 2024-02-09 03:43:24 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - List available models/likelihoods/etc with inla.list.models()
## - Use inla.doc(<NAME>) to access documentation#inla.update(testing=TRUE)
library(INLAjoint)## Package 'INLAjoint' version 24.2.5
## Type 'citation("INLAjoint")' for citing this R package in publications.library(splines)Importing data
surv_long_ <- read.csv("processed/surv_long_.csv")
# We can remove the following after getting the new treatment data
# The new treatment data is supposed to contain only RP or RT when status=1 and UNK when status=0
# Now we filter the available data
survlong_rprt <- surv_long_[surv_long_$tx1_type %in% c('rp','radio','tfs-unk'),]
# Rename attributes
names(survlong_rprt)[names(survlong_rprt) == "tx1_type"] <- "tfs_tx"
names(survlong_rprt)[names(survlong_rprt) == "after_tx"] <- "drfs_tx"
# Change values
survlong_rprt$tfs_tx[survlong_rprt$tfs_tx %in% c('radio')] <- 'rt'
survlong_rprt$tfs_tx[survlong_rprt$tfs_tx %in% c('tfs-unk')] <- 'unk'
survlong_rprt$drfs_tx[survlong_rprt$drfs_tx %in% c('bcr-unk')] <- 'unk'
# There shouldn't be any time to event = 0 (checked status 0, no time 0, good!) ,
# drfs_time min = 1 , good!
survlong_rprt$tfs_time[(survlong_rprt$tfs_status == 1) & (survlong_rprt$tfs_time == 0)] <- 0.5
# Convert time to event from Months to Years
survlong_rprt$tfs_time_m <- survlong_rprt$tfs_time
survlong_rprt$drfs_time_m <- survlong_rprt$drfs_time
survlong_rprt$time_to_psa_m <- survlong_rprt$time_to_psa
survlong_rprt$tfs_time <- round(survlong_rprt$tfs_time_m/12, 2)
survlong_rprt$drfs_time <- round(survlong_rprt$drfs_time_m/12, 2)
survlong_rprt$time_to_psa <- round(survlong_rprt$time_to_psa_m/12, 2)Functions
reassigning_ids <- function(data){
id_rle <- rle(data$id)
new_id <- rep(seq_along(id_rle$lengths), id_rle$lengths)
data$id <- new_id
return(data)
}
bf_sampling <- function(data, sample_size=0.3){
#data <- survlong_rprt
#sample_size <- 1
# Uncomment when needed, Sampling
unique_ids <- unique(data$id)
sample_size <- round(sample_size * length(unique_ids))
# Sample X.% of the unique IDs
sampled_ids <- sample(unique_ids, sample_size, replace = FALSE)
# Subset the data based on sampled IDs
sampled_data <- data[data$id %in% sampled_ids, ]
# Factors
sampled_data$mri <- as.factor(sampled_data$mri)
sampled_data$grade_group <- as.factor(sampled_data$grade_group)
sampled_data$bx_age <- as.factor(sampled_data$bx_age)
sampled_data$tfs_status <- as.factor(sampled_data$tfs_status)
sampled_data$drfs_status <- as.factor(sampled_data$drfs_status)
sampled_data$drfs_tx <- as.factor(sampled_data$drfs_tx)
sampled_data$tfs_tx <- as.factor(sampled_data$tfs_tx)
# Order factors
sampled_data$tfs_tx <- factor(sampled_data$tfs_tx, levels = c('rt', 'rp', 'unk'))
sampled_data$drfs_tx <- factor(sampled_data$drfs_tx, levels = c('rt', 'rp', 'unk'))
# Set reference
sampled_data$tfs_tx <- relevel(sampled_data$tfs_tx, ref = "rt")
sampled_data$drfs_tx <- relevel(sampled_data$drfs_tx, ref = "rt")
# Remove Nas
sampled_data <- sampled_data[!is.na(sampled_data$time_to_psa),]
# Reassigning ID
# We have to do it here, as all will be reordered after the previous cmd
# -> Calculate run-length encoding of the original id column
id_rle <- rle(sampled_data$id)
# -> Generate a sequence of IDs based on the lengths of runs in the original id column
new_id <- rep(seq_along(id_rle$lengths), id_rle$lengths)
# -> Reassign the new IDs to the id column in sampled_data
sampled_data$id <- new_id
sampled_data <- sampled_data[with(sampled_data, order(id)), ]
# Order it by :
# 1- Type of treatment at recurrence drfs_tx, so the unk ones will be at the end.
# 2- Time to psa , so the INLAjoint Long model works properly
# We already have PSA at t=0 for all these patients from part 1
#sampled_data <- sampled_data[with(sampled_data, order(drfs_tx, time_to_psa)), ]
sampled_data <- sampled_data[order(sampled_data$drfs_tx,
sampled_data$id,
sampled_data$time_to_psa), ]
# Quick reordering
id_rle <- rle(sampled_data$id)
new_id <- rep(seq_along(id_rle$lengths), id_rle$lengths)
sampled_data$id <- new_id
return(sampled_data)
}
interp <- function(jm){
n_sample <- 1e4 # number of samples for uncertainty quantification
smp_H <- inla.hyperpar.sample(n_sample, jm) # sample values for frailty and association parameters
sigma <- sqrt(1/smp_H[, which(colnames(smp_H)=="Precision for IDIntercept_S1")]) # sampled sd of frailty
assoc <- smp_H[, which(colnames(smp_H)=="Beta for IDIntercept_S1_S2")] # sampled association parameter
SMP <- sapply(sigma, function(x) rnorm(1e3, mean = 0, sd = x)) # sample realizations of each frailty
SMP2 <- rbind(sigma, SMP) # add sigmas there to vectorize computations and avoid loop
mean_low15 <- apply(SMP, 2, function(x) mean(x[-1][x[-1]<(-x[1])])) # mean frailty deviation for lowest 15% , for 2.5% -2x[1]
mean_up15 <- apply(SMP, 2, function(x) mean(x[-1][x[-1]>x[1]])) # mean frailty deviation for top 15%, same here
HR_low15 <- exp(assoc*mean_low15) # hazard ratios for lowest 15% vs. average
HR_up15 <- exp(assoc*mean_up15) # hazard ratios for top 15% vs. average
Q_low15 <- quantile(HR_low15, c(0.025, 0.5, 0.975), na.rm=T) # should not have NAs, here it's a bad example on unstable model
Q_up15 <- quantile(HR_up15, c(0.025, 0.5, 0.975), na.rm=T)
# Interpretation:
res_l15 <- paste0("Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the lowest risk of receiving treatment (i.e., longest time-to-treatment) is associated to a ", round(Q_low15[2], 2), " [", round(Q_low15[1], 2), ",", round(Q_low15[3], 2), "] increased risk of recurrence (i.e., shorter time-to-recurrence) compared to the average individual.")
res_u15 <- paste0("Conditional on covariates included in the model (i.e., grade group, mri, ...) and on the time-dependent PSA level, the top 15% individuals with the highest risk of receiving treatment (i.e., shortest time-to-treatment) is associated to a ", round(Q_up15[2], 2), " [", round(Q_up15[1], 2), ",", round(Q_up15[3], 2), "] decreased risk of recurrence (i.e., longer time-to-recurrence) compared to the average individual.")
return(c(res_l15, res_u15))
}Modeling
Sampling
# Sampling
s_survlong_rprt <- bf_sampling(survlong_rprt, 0.9)
# As we are modeling only BCR (no TFS), we need only treated patients (cannot truncate when no tx)
s_survlong_rprt_bcr <- s_survlong_rprt[s_survlong_rprt$tfs_status==1,]
s_survlong_rprt_bcr$drfs_tx <- droplevels(s_survlong_rprt_bcr$drfs_tx)
s_survlong_rprt_bcr$tfs_tx <- droplevels(s_survlong_rprt_bcr$tfs_tx)
# Surv (only) data to be used in the latter join model
s_surv_rprt_bcr <- s_survlong_rprt_bcr[!duplicated(s_survlong_rprt_bcr$id), ]
summary(s_survlong_rprt_bcr)## id bx_age mri grade_group tfs_tx
## Min. : 1 <60 :20358 0:96632 1:35675 rt:65615
## 1st Qu.:1494 >70 :34209 1:11106 2:31501 rp:42123
## Median :2949 60-70:53171 3:22500
## Mean :2985 4:10054
## 3rd Qu.:4479 5: 8008
## Max. :6039
## tx_date tfs_time tfs_status drfs_tx drfs_time
## Length:107738 Min. : 0.040 0: 0 rt:65615 Min. : 0.080
## Class :character 1st Qu.: 0.250 1:107738 rp:42123 1st Qu.: 3.580
## Mode :character Median : 0.330 Median : 6.330
## Mean : 1.204 Mean : 6.816
## 3rd Qu.: 0.920 3rd Qu.: 9.670
## Max. :20.000 Max. :23.670
## drfs_status psa_results psa_date time_to_psa
## 0:66600 Min. : 0.00 Length:107738 Min. :-20.170
## 1:41138 1st Qu.: 0.04 Class :character 1st Qu.: 0.330
## Median : 0.78 Mode :character Median : 2.250
## Mean : 11.66 Mean : 2.802
## 3rd Qu.: 5.38 3rd Qu.: 5.250
## Max. :76800.00 Max. : 23.670
## log_psa tfs_time_m drfs_time_m time_to_psa_m
## Min. : 0.00000 Min. : 0.50 Min. : 1.00 Min. :-242.00
## 1st Qu.: 0.03922 1st Qu.: 3.00 1st Qu.: 43.00 1st Qu.: 4.00
## Median : 0.57661 Median : 4.00 Median : 76.00 Median : 27.00
## Mean : 1.03175 Mean : 14.44 Mean : 81.79 Mean : 33.62
## 3rd Qu.: 1.85395 3rd Qu.: 11.00 3rd Qu.:116.00 3rd Qu.: 63.00
## Max. :11.24897 Max. :240.00 Max. :284.00 Max. : 284.00Fit
# Setup the number of used threads
inla.setOption(num.threads='8:1')
NSplines <- ns(s_survlong_rprt_bcr$time_to_psa, knots=c(3)) # natural cubic splines (always check)
f1 <- function(x) predict(NSplines, x)[,1] # first basis
f2 <- function(x) predict(NSplines, x)[,2] # second basis
JM3 <- joint(formSurv=list(inla.surv(drfs_time, drfs_status) ~ drfs_tx+mri),
id="id",
formLong=log_psa ~ (f1(time_to_psa) + f2(time_to_psa)) * drfs_tx+mri + (1 + f1(time_to_psa) + f2(time_to_psa)| id),
timeVar="time_to_psa",
basRisk=c("rw2"),
family="gaussian",
assoc=c("CV"),
dataSurv = s_surv_rprt_bcr,
dataLong = s_survlong_rprt_bcr,
control=list(int.strategy="eb"))## Warning in joint(formSurv = list(inla.surv(drfs_time, drfs_status) ~ drfs_tx +
## : Internal correlation between hyperparameters is abnormally high, this is a
## sign of identifiability issues / ill-defined model.summary(JM3)## Longitudinal outcome (gaussian)
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept_L1 12.3704 0.3201 11.7431 12.3704 12.9978
## f1time_to_psa_L1 -20.1218 0.5254 -21.1516 -20.1218 -19.0920
## f2time_to_psa_L1 -1.0193 0.2301 -1.4703 -1.0193 -0.5683
## drfs_txrp_L1 3.2850 0.5104 2.2846 3.2850 4.2853
## mri1_L1 0.0157 0.0271 -0.0375 0.0157 0.0688
## f1time_to_psa:drfs_txrp_L1 -7.6629 0.8350 -9.2994 -7.6629 -6.0263
## f2time_to_psa:drfs_txrp_L1 -1.2883 0.3687 -2.0109 -1.2883 -0.5658
## Res. err. (variance) 0.3922 0.0017 0.3885 0.3924 0.3949
##
## Random effects variance-covariance (L1)
## mean sd 0.025quant 0.5quant
## Intercept_L1 278.1286 7.3039 263.0206 278.5591
## f1time_to_psa_L1 764.5609 18.7760 726.2997 765.6753
## f2time_to_psa_L1 136.5785 4.4717 128.7658 136.3160
## Intercept_L1:f1time_to_psa_L1 -455.7412 11.6738 -477.3622 -456.5812
## Intercept_L1:f2time_to_psa_L1 118.1179 5.1603 108.1461 118.1931
## f1time_to_psa_L1:f2time_to_psa_L1 -155.9696 8.2385 -171.6086 -156.2902
## 0.975quant
## Intercept_L1 291.4873
## f1time_to_psa_L1 799.7444
## f2time_to_psa_L1 145.7816
## Intercept_L1:f1time_to_psa_L1 -431.7903
## Intercept_L1:f2time_to_psa_L1 128.1304
## f1time_to_psa_L1:f2time_to_psa_L1 -139.4775
##
## Survival outcome
## mean sd 0.025quant 0.5quant 0.975quant
## Baseline risk (variance)_S1 0.0528 0.0008 0.0513 0.0528 0.0544
## drfs_txrp_S1 0.6798 0.0242 0.6323 0.6798 0.7272
## mri1_S1 0.9709 0.0373 0.8977 0.9709 1.0441
##
## Association longitudinal - survival
## mean sd 0.025quant 0.5quant 0.975quant
## CV_L1_S1 0.4928 0.0099 0.4704 0.4937 0.5085
##
## log marginal-likelihood (integration) log marginal-likelihood (Gaussian)
## -252982.6 -252974.0
##
## Deviance Information Criterion: 55313.84
## Widely applicable Bayesian information criterion: 49190.38
## Computation time: 135.34 secondsPredict
ND1 <- s_survlong_rprt_bcr[s_survlong_rprt_bcr$id==1,] # make prediction for this individual
ND2 <- s_survlong_rprt_bcr[c(1,1,1,1),] # first line 4 time
ND2$id <- 2:5
ND2$drfs_tx <- c("rp", "rt","rp", "rt")
ND2$mri <- c(0, 1, 1,0)
ND2$log_psa <- NA # make predictions for average trajectories conditional on treatment
ND <- rbind(ND1, ND2)
ND## id bx_age mri grade_group tfs_tx tx_date tfs_time tfs_status drfs_tx
## 36 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 37 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 35 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 30 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 44 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 41 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 32 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 39 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 38 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 40 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 33 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 43 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 31 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 42 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 34 1 60-70 0 3 rt 2016-04-11 0.25 1 rt
## 361 2 60-70 0 3 rt 2016-04-11 0.25 1 rp
## 36.1 3 60-70 1 3 rt 2016-04-11 0.25 1 rt
## 36.2 4 60-70 1 3 rt 2016-04-11 0.25 1 rp
## 36.3 5 60-70 0 3 rt 2016-04-11 0.25 1 rt
## drfs_time drfs_status psa_results psa_date time_to_psa log_psa
## 36 4.83 0 63.00 2015-08-21 -0.42 4.1588831
## 37 4.83 0 36.00 2015-10-01 -0.25 3.6109179
## 35 4.83 0 36.00 2016-01-05 0.00 3.6109179
## 30 4.83 0 8.30 2016-04-11 0.25 2.2300144
## 44 4.83 0 8.20 2016-12-12 0.92 2.2192035
## 41 4.83 0 7.70 2017-02-14 1.08 2.1633230
## 32 4.83 0 7.50 2017-04-20 1.25 2.1400662
## 39 4.83 0 2.40 2017-07-26 1.50 1.2237754
## 38 4.83 0 1.40 2018-01-16 2.00 0.8754687
## 40 4.83 0 1.10 2018-06-11 2.42 0.7419373
## 33 4.83 0 0.92 2018-09-06 2.67 0.6523252
## 43 4.83 0 0.52 2019-04-02 3.25 0.4187103
## 31 4.83 0 0.41 2019-10-15 3.75 0.3435897
## 42 4.83 0 0.35 2020-04-17 4.25 0.3001046
## 34 4.83 0 0.35 2020-11-03 4.83 0.3001046
## 361 4.83 0 63.00 2015-08-21 -0.42 NA
## 36.1 4.83 0 63.00 2015-08-21 -0.42 NA
## 36.2 4.83 0 63.00 2015-08-21 -0.42 NA
## 36.3 4.83 0 63.00 2015-08-21 -0.42 NA
## tfs_time_m drfs_time_m time_to_psa_m
## 36 3 58 -5
## 37 3 58 -3
## 35 3 58 0
## 30 3 58 3
## 44 3 58 11
## 41 3 58 13
## 32 3 58 15
## 39 3 58 18
## 38 3 58 24
## 40 3 58 29
## 33 3 58 32
## 43 3 58 39
## 31 3 58 45
## 42 3 58 51
## 34 3 58 58
## 361 3 58 -5
## 36.1 3 58 -5
## 36.2 3 58 -5
## 36.3 3 58 -5# survival is : to calculate the survival in addition to risk (integrate area uner the curve),
# Nsample is: number of samples from the covariates dist.
# NsampleRE is: for each sample, this is the number of samples for the random effect.
PRED <- predict(JM3, ND, horizon=max(s_survlong_rprt_bcr$time_to_psa), survival=TRUE, Nsample=500, NsampleRE=100) ## Start sampling## Sampling done.## Computing longitudinal predictions for individual 1 on PID: 71810## Computing survival predictions for individual 1## Computing longitudinal predictions for individual 2 on PID: 71810## Computing survival predictions for individual 2## Computing longitudinal predictions for individual 3 on PID: 71810## Computing survival predictions for individual 3## Computing longitudinal predictions for individual 4 on PID: 71810## Computing survival predictions for individual 4## Computing longitudinal predictions for individual 5 on PID: 71810## Computing survival predictions for individual 5# This function will faster when the code is written in C++ (later this year!)Results
PRED_1 <- PRED
PRED_1$PredL <- PRED_1$PredL[PRED_1$PredL$id==1, ]# longitudinal prediction for individual 1
PRED_1$PredS <- PRED_1$PredS[PRED_1$PredS$id==1, ]# survival prediction for individual 1
# plot observed vs. predicted longitudinal for individual 1
plot(PRED_1$PredL$time_to_psa, PRED_1$PredL$quant0.5, ylim=range(c(s_survlong_rprt_bcr$log_psa,
PRED_1$PredL$quant0.025,
PRED_1$PredL$quant0.975)), type="o", pch=19)
lines(PRED_1$PredL$time_to_psa, PRED_1$PredL$quant0.025, type="o", pch=19, cex=0.5)
lines(PRED_1$PredL$time_to_psa, PRED_1$PredL$quant0.975, type="o", pch=19, cex=0.5)
points(ND1$time_to_psa, ND1$log_psa, pch=19, col="red")
# plot survival probability for individual 1
plot(PRED_1$PredS$time_to_psa, PRED_1$PredS$Surv_quant0.5, ylim=c(0,1), type="o", pch=19)
lines(PRED_1$PredS$time_to_psa, PRED_1$PredS$Surv_quant0.025, type="o", pch=19, cex=0.5)
lines(PRED_1$PredS$time_to_psa, PRED_1$PredS$Surv_quant0.975, type="o", pch=19, cex=0.5)
‘Predictions + Interpretation’ for avg individual conditional on treatment
# now we make predictions for average individual conditional on treatment
PRED_T <- PRED$PredS # only interested in survival part of predictions?
PRED_T <- PRED_T[PRED_T$id!=1,]
# If we have more scenarios, i.e more cov, it's better to predict inside the function and use the results, instead of predicting all cases and play with large if else.
FCT <- function(X, TRT, MRI){ # function with arguments X (time) and TRT (treatment)
if(TRT=="rp" & MRI==1) ID <- 4
else if(TRT=="rt" & MRI==1) ID <- 3
else if(TRT=="rp" & MRI==0) ID <- 2
else if(TRT=="rt" & MRI==0) ID <- 1 # choose the individual associated to the given treatment line for the prediction (this could be grade_group for you)
PRED_T_i <- PRED_T[PRED_T$id==ID,] # select the prediction for the given category of the covariate(s), here only treatment
fpred_0025 <- splinefun(PRED_T_i$time_to_psa, PRED_T_i$Surv_quant0.025, method="monoH.FC")
fpred_05 <- splinefun(PRED_T_i$time_to_psa, PRED_T_i$Surv_quant0.5, method="monoH.FC")
fpred_0975 <- splinefun(PRED_T_i$time_to_psa, PRED_T_i$Surv_quant0.975, method="monoH.FC")
return(c(fpred_0025(X), fpred_05(X), fpred_0975(X)))
}
X <- 1 # time of prediction
P1_C <- round((1-FCT(X, "rp", 1))*100, 2)
P1_S <- round((1-FCT(X, "rt", 1))*100, 2)
paste0("An average individual has ", P1_C[2], "% [", P1_C[3], ", ", P1_C[1], "] probability of event at time ", X, " if he receives treatment RP (MRI=1).")## [1] "An average individual has 41.93% [24.85, 64.74] probability of event at time 1 if he receives treatment RP (MRI=1)."paste0("An average individual has ", P1_S[2], "% [", P1_S[3], ", ", P1_S[1], "] probability of event at time ", X, " if he receives treatment RT (MRI=1).")## [1] "An average individual has 31.54% [17.87, 51.46] probability of event at time 1 if he receives treatment RT (MRI=1)."# important point, predictions are always conditional on covariates included in the model, here the treatment only
# if you want more general predictions (unconditional) then do not include covariates in the model (and you may have higher uncertainty)
FCT_inv <- function(X, TRT, MRI){ # function with arguments X (survival probability) and TRT (treatment)
if(TRT=="rp" & MRI==1) ID <- 4
else if(TRT=="rt" & MRI==1) ID <- 3
else if(TRT=="rp" & MRI==0) ID <- 2
else if(TRT=="rt" & MRI==0) ID <- 1 # choose the individual associated to the given treatment line for the prediction (this could be grade_group for you)
PRED_T_i <- PRED_T[PRED_T$id==ID,] # select the prediction for the given category of the covariate(s), here only treatment
i_fpred_025 <- splinefun(PRED_T_i$Surv_quant0.025, PRED_T_i$time_to_psa, method="monoH.FC")
i_fpred_05 <- splinefun(PRED_T_i$Surv_quant0.5, PRED_T_i$time_to_psa, method="monoH.FC")
i_fpred_0975 <- splinefun(PRED_T_i$Surv_quant0.975, PRED_T_i$time_to_psa, method="monoH.FC")
return(c(i_fpred_025(X), i_fpred_05(X), i_fpred_0975(X)))
}
X <- 0.6 # survival probability (i.e., 1 - probability of event)
P1_C <- round(FCT_inv(X, "rp", 1), 2)## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' valuesP1_S <- round(FCT_inv(X, "rt", 1), 2)## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' valuespaste0("An average individual has ", round((1-X)*100,2), "% probability of event at time ", P1_C[2], " [", P1_C[1], ", ", P1_C[3], "] if he receives treatment RP (MRI=1).")## [1] "An average individual has 40% probability of event at time 0.93 [0.36, 1.86] if he receives treatment RP (MRI=1)."paste0("An average individual has ", round((1-X)*100,2), "% probability of event at time ", P1_S[2], " [", P1_S[1], ", ", P1_S[3], "] if he receives treatment RT (MRI=1).")## [1] "An average individual has 40% probability of event at time 1.36 [0.66, 2.62] if he receives treatment RT (MRI=1)."