Last active
August 12, 2025 02:43
-
-
Save sashagusev/8755afd308e7a039a4d6c0e04b386bdc to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library(survival) | |
| # --- | |
| # Parameters | |
| # --- | |
| # number of individuals | |
| n = 50e3 | |
| # desired hazard ratio | |
| HR = 1.64 | |
| beta = log(HR) | |
| risk_reduction = 0.5 | |
| # --- | |
| # Simulate a normally distributed risk factor | |
| X = rnorm(n, mean = 0, sd = 1) | |
| U = runif(n) | |
| # Generate survival times from Weibull proportional hazards model | |
| shape = 4 | |
| scale = 140 | |
| T <- scale * (-log(U) / exp(beta * X))^(1 / shape) | |
| # Simulate an intervention that lowers risk by 0.5SDs | |
| X_new = X - risk_reduction | |
| T_new <- scale * (-log(U) / exp(beta * X_new))^(1 / shape) | |
| par(mfrow=c(1,4)) | |
| # --- Plot the cumulative incidence | |
| # Fit a survival curve | |
| fit <- survfit(Surv(T, rep(1,length(T))) ~ 1) | |
| plot(fit$time, 1 - fit$surv, type = "s", lwd = 2, | |
| xlab = "Age", ylab = "Cumulative incidence", | |
| main = "Cumulative incidence curve", | |
| col = "#08519c",las=1,bty="n",xlim=c(0,150)) | |
| fit <- survfit(Surv(T_new, rep(1,length(T_new))) ~ 1) | |
| lines( fit$time, 1 - fit$surv , col="#9ecae1",lwd=2) | |
| fit <- survfit(Surv(death_times, rep(1,length(death_times))) ~ 1) | |
| legend("bottomright",pch=19,legend=c("Baseline","Treated"),col=c("#08519c","#9ecae1"),bty="n") | |
| # --- Plot years gained | |
| max_age = 100 | |
| min_age = 50 | |
| ages = seq(min_age,max_age,by=1) | |
| plot(T[T>min_age&T<max_age],(T_new-T)[T>min_age&T<max_age],xlim=c(min_age,max_age),type="p",xlab="Age at diagnosis",ylab="Disease-free years gained",cex=0.5,pch=19,bty="n",las=1,main="Years gained post-diagnosis") | |
| grid() | |
| # Add a death curve | |
| # Sample mortality from a Weibull | |
| death_times <- rweibull(n, shape = 5, scale = 88) | |
| # Truncate mortality at age 100 | |
| death_times[ death_times > 100 ] = 100 | |
| # Calculate disease free time prevented on average | |
| dtime_a = death_times - T | |
| dtime_a[ dtime_a < 0 ] = 0 | |
| dtime_b = death_times - T_new | |
| dtime_b[ dtime_b < 0 ] = 0 | |
| disease_time = rep(NA,length(ages)) | |
| for ( i in 1:length(ages) ) { | |
| dtime_a = ages[i] - T | |
| dtime_a[ dtime_a < 0 ] = 0 | |
| dtime_b = ages[i] - T_new | |
| dtime_b[ dtime_b < 0 ] = 0 | |
| disease_time[i] = mean(dtime_a - dtime_b) | |
| cat( ages[i] , disease_time[i] , '\n' ) | |
| } | |
| plot(ages,disease_time*12,ylab="Disease-free months gained",xlab="Age at censoring",type="o",pch=19,cex=0.5,lwd=2,bty="n",las=1,main="Population expected months gained") | |
| # Calculate the classic relative risk reduction | |
| # --- Resimulate with more data | |
| n <- 10e6 | |
| X <- rnorm(n, mean = 0, sd = 1) | |
| U <- runif(n) | |
| T <- scale * (-log(U) / exp(beta * X))^(1 / shape) | |
| X_new = X - risk_reduction | |
| T_new <- scale * (-log(U) / exp(beta * X_new))^(1 / shape) | |
| # Calculate the relative risk reduction by age cutoff | |
| RR = rep(NA,length(ages)) | |
| for ( i in 1:length(ages) ) { | |
| RR[i] = (mean(T_new<ages[i]) / mean(T<ages[i])) | |
| cat( ages[i] , mean(T_new<ages[i]) , mean(T<ages[i]) , 1-RR[i] , '\n' ) | |
| } | |
| # Plot risk reduction | |
| plot(ages,(1-RR)*100,type="o",bty="n",xlim=c(min_age,max_age),cex=0.5,pch=19,xlab="Age at censoring",ylab="Relative risk reduction (%)",las=1,main="Relative risk reduction") | |
| grid() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment