# =============================================================================
# LEGACY STUDY HARMONIZATION - ANALYSIS
# Meta-analysis, patient-level pooling, endpoint re-derivation and prediction
#
# Developed by: Planimeter Kft.
# Version:      1.0.0
# Date:         2026-09-24
# Copyright:    (c) 2026 Planimeter Kft.
# SPDX-License-Identifier: GPL-3.0-or-later
#
# Purpose:
#   Reproduce the final analytical story from the harmonized patient-level
#   legacy dataset created by legacy_data_preparation_planimeter_v1.0.0.R.
# =============================================================================

required_packages <- c("dplyr","readr","ggplot2","survival","broom",
                       "ranger","yardstick","tidyr","tibble","scales","grid")
missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages)>0) stop(paste("Missing required packages:",paste(missing_packages,collapse=", ")))

library(dplyr); library(readr); library(ggplot2); library(survival); library(broom)
library(ranger); library(yardstick); library(tidyr); library(tibble); library(scales); library(grid)

dir.create("output",recursive=TRUE,showWarnings=FALSE)
dir.create("figures",recursive=TRUE,showWarnings=FALSE)

integrated <- read_csv(file.path("data","generated","integrated_legacy_breast_cancer.csv"), show_col_types=FALSE) %>%
  mutate(treatment=factor(treatment,levels=c("Control","Treatment")), grade=factor(grade))

# ---- 1. Study-level treatment effects --------------------------------------
study_meta <- integrated %>% group_by(study_id) %>% group_modify(~{
  fit <- coxph(Surv(analysis_time_years,analysis_event)~treatment,data=.x)
  tidy(fit) %>% filter(term=="treatmentTreatment") %>% transmute(logHR=estimate,SE=std.error)
}) %>% ungroup() %>% mutate(HR=exp(logHR),conf.low=exp(logHR-1.96*SE),conf.high=exp(logHR+1.96*SE))

fixed_meta <- study_meta %>% mutate(weight=1/SE^2) %>% summarise(
  logHR=sum(weight*logHR)/sum(weight), SE=sqrt(1/sum(weight))) %>% mutate(
  study_id="Fixed-effect meta-analysis", HR=exp(logHR),
  conf.low=exp(logHR-1.96*SE), conf.high=exp(logHR+1.96*SE), method="Meta-analysis")

k <- nrow(study_meta)
q_df <- study_meta %>% mutate(weight_fixed=1/SE^2)
fixed_logHR <- sum(q_df$weight_fixed*q_df$logHR)/sum(q_df$weight_fixed)
Q <- sum(q_df$weight_fixed*(q_df$logHR-fixed_logHR)^2)
C <- sum(q_df$weight_fixed)-sum(q_df$weight_fixed^2)/sum(q_df$weight_fixed)
tau2 <- max(0,(Q-(k-1))/C)
random_meta <- study_meta %>% mutate(weight_random=1/(SE^2+tau2)) %>% summarise(
  logHR=sum(weight_random*logHR)/sum(weight_random), SE=sqrt(1/sum(weight_random))) %>% mutate(
  study_id="Random-effects meta-analysis", HR=exp(logHR),
  conf.low=exp(logHR-1.96*SE), conf.high=exp(logHR+1.96*SE), method="Meta-analysis")

pooled_fit <- coxph(Surv(analysis_time_years,analysis_event)~treatment+strata(study_id),data=integrated)
patient_pooled <- tidy(pooled_fit) %>% filter(term=="treatmentTreatment") %>% transmute(
  study_id="Patient-level pooled analysis", logHR=estimate, SE=std.error,
  HR=exp(logHR), conf.low=exp(logHR-1.96*SE), conf.high=exp(logHR+1.96*SE), method="Patient-level pooled")

forest_meta_df <- bind_rows(
  study_meta %>% mutate(method="Individual legacy study") %>% select(study_id,logHR,SE,HR,conf.low,conf.high,method),
  fixed_meta %>% select(study_id,logHR,SE,HR,conf.low,conf.high,method),
  random_meta %>% select(study_id,logHR,SE,HR,conf.low,conf.high,method),
  patient_pooled %>% select(study_id,logHR,SE,HR,conf.low,conf.high,method)
)
write_csv(forest_meta_df,file.path("output","meta_analysis_vs_patient_level.csv"))

# ---- 2. Common 3-year and 5-year endpoints ---------------------------------
integrated_endpoints <- integrated %>% mutate(
  event_3y=ifelse(analysis_time_years<=3 & analysis_event==1,1,0), time_3y=pmin(analysis_time_years,3),
  event_5y=ifelse(analysis_time_years<=5 & analysis_event==1,1,0), time_5y=pmin(analysis_time_years,5))

fit_3y <- coxph(Surv(time_3y,event_3y)~treatment+strata(study_id),data=integrated_endpoints)
fit_5y <- coxph(Surv(time_5y,event_5y)~treatment+strata(study_id),data=integrated_endpoints)
endpoint_results <- bind_rows(
  tidy(fit_3y,exponentiate=TRUE,conf.int=TRUE) %>% filter(term=="treatmentTreatment") %>% mutate(endpoint="Re-derived 3-year endpoint"),
  tidy(fit_5y,exponentiate=TRUE,conf.int=TRUE) %>% filter(term=="treatmentTreatment") %>% mutate(endpoint="Re-derived 5-year endpoint")
) %>% transmute(endpoint,HR=estimate,conf.low,conf.high,p.value)
write_csv(endpoint_results,file.path("output","rederived_endpoint_results.csv"))

# ---- 3. Patient-level prediction -------------------------------------------
ml_df <- integrated_endpoints %>% mutate(
  endpoint_3y=factor(event_3y,levels=c(0,1)), endpoint_5y=factor(event_5y,levels=c(0,1)),
  treatment=factor(treatment), study_id=factor(study_id), grade=factor(grade)) %>%
  select(study_id,age,tumor_size_mm,nodes,grade,treatment,endpoint_3y,endpoint_5y) %>% na.omit()

set.seed(20260610)
idx <- sample(seq_len(nrow(ml_df)),size=floor(0.70*nrow(ml_df)))
train_ml <- ml_df[idx,]; test_ml <- ml_df[-idx,]
rf_3y <- ranger(endpoint_3y~study_id+age+tumor_size_mm+nodes+grade+treatment,
                data=train_ml,probability=TRUE,importance="impurity",num.trees=500,seed=20260610)
rf_5y <- ranger(endpoint_5y~study_id+age+tumor_size_mm+nodes+grade+treatment,
                data=train_ml,probability=TRUE,importance="impurity",num.trees=500,seed=20260610)
pred_3y <- predict(rf_3y,data=test_ml)$predictions[,"1"]
pred_5y <- predict(rf_5y,data=test_ml)$predictions[,"1"]
auc_3y <- roc_auc_vec(test_ml$endpoint_3y,pred_3y,event_level="second")
auc_5y <- roc_auc_vec(test_ml$endpoint_5y,pred_5y,event_level="second")
write_csv(tibble(endpoint=c("3-year endpoint","5-year endpoint"),auc=c(auc_3y,auc_5y)),
          file.path("output","prediction_auc.csv"))

vip_df <- bind_rows(
  tibble(variable=names(rf_3y$variable.importance),importance=as.numeric(rf_3y$variable.importance),endpoint="3-year endpoint"),
  tibble(variable=names(rf_5y$variable.importance),importance=as.numeric(rf_5y$variable.importance),endpoint="5-year endpoint")
) %>% group_by(endpoint) %>% mutate(importance_scaled=importance/max(importance)) %>% ungroup()
write_csv(vip_df,file.path("output","random_forest_variable_importance.csv"))

calibration_df <- bind_rows(
  tibble(endpoint="3-year endpoint",pred=pred_3y,observed=as.integer(as.character(test_ml$endpoint_3y))),
  tibble(endpoint="5-year endpoint",pred=pred_5y,observed=as.integer(as.character(test_ml$endpoint_5y)))
) %>% group_by(endpoint) %>% mutate(risk_decile=ntile(pred,10)) %>%
  group_by(endpoint,risk_decile) %>% summarise(mean_predicted=mean(pred),observed_rate=mean(observed),n=n(),.groups="drop")
write_csv(calibration_df,file.path("output","prediction_calibration.csv"))

# ---- 4. Final figures -------------------------------------------------------
bg="#07111f"; gridc="#24384F"; txt="white"; subtxt="#A9BED6"; cyan="#55DDE0"; gold="#FFD166"; purple="#7B61FF"; muted="#132238"
theme_dark <- function(base_size=16) theme_minimal(base_size=base_size)+theme(
  plot.background=element_rect(fill=bg,color=NA), panel.background=element_rect(fill=bg,color=NA),
  panel.grid.major=element_line(color=gridc,linewidth=.35), panel.grid.minor=element_blank(),
  plot.title=element_text(color=txt,face="bold",size=base_size+8), plot.subtitle=element_text(color=subtxt),
  axis.title=element_text(color=subtxt), axis.text=element_text(color=subtxt),
  legend.background=element_rect(fill=bg,color=NA), legend.key=element_rect(fill=bg,color=NA),
  legend.text=element_text(color=subtxt), legend.title=element_text(color=txt),
  strip.text=element_text(color=txt,face="bold"), strip.background=element_rect(fill=muted,color=NA))
save_fig <- function(p,n,w=15,h=8) ggsave(file.path("figures",n),p,width=w,height=h,dpi=300,bg=bg)

legacy_map <- read_csv(file.path("data","legacy_study_specification.csv"),show_col_types=FALSE) %>%
  select(study_id,source_endpoint,source_time_unit,source_treatment_coding,source_tumor_size_representation,source_node_representation) %>%
  pivot_longer(-study_id,names_to="feature",values_to="value")
p1 <- ggplot(legacy_map,aes(study_id,feature))+geom_tile(fill=muted,color=bg)+geom_text(aes(label=value),color=txt,size=3.6)+
  labs(title="Six legacy studies are similar, but not identical",subtitle="Endpoint definitions, time units and coding conventions must be harmonized",x=NULL,y=NULL)+theme_dark(14)+theme(panel.grid=element_blank())
save_fig(p1,"figure_01_legacy_heterogeneity_matrix.png",16,8)

plot_meta <- forest_meta_df %>% mutate(study_id=factor(study_id,levels=rev(study_id)))
p2 <- ggplot(plot_meta,aes(study_id,HR,color=method))+geom_hline(yintercept=1,linetype="dashed",color=subtxt)+
  geom_errorbar(aes(ymin=conf.low,ymax=conf.high),width=.15)+geom_point(size=3)+scale_y_log10()+coord_flip()+
  labs(title="Meta-analysis summarizes fragmented evidence",subtitle="Study-level synthesis and patient-level pooling can be compared on the same effect scale",x=NULL,y="Hazard ratio, log scale",color=NULL)+theme_dark()
save_fig(p2,"figure_02_meta_analysis_and_patient_level.png",16,9)

capability_df <- tribble(~approach,~capability,~available,
 "Study-level meta-analysis","Common treatment effect",1,"Study-level meta-analysis","Forest plot",1,
 "Study-level meta-analysis","Re-derived 3-year endpoint",0,"Study-level meta-analysis","Re-derived 5-year endpoint",0,
 "Study-level meta-analysis","Covariate-adjusted patient-level model",0,"Study-level meta-analysis","Predictive modelling",0,
 "Patient-level harmonization","Common treatment effect",1,"Patient-level harmonization","Forest plot",1,
 "Patient-level harmonization","Re-derived 3-year endpoint",1,"Patient-level harmonization","Re-derived 5-year endpoint",1,
 "Patient-level harmonization","Covariate-adjusted patient-level model",1,"Patient-level harmonization","Predictive modelling",1)
p3 <- ggplot(capability_df,aes(approach,capability,fill=factor(available)))+geom_tile(color=bg)+geom_text(aes(label=ifelse(available==1,"YES","-")),color=txt,fontface="bold")+
  scale_fill_manual(values=c("0"="#301820","1"="#1E3A34"),guide="none")+labs(title="Patient-level harmonization enables analyses beyond meta-analysis",subtitle="The value is a reusable integrated evidence dataset",x=NULL,y=NULL)+theme_dark()+theme(panel.grid=element_blank())
save_fig(p3,"figure_03_beyond_meta_analysis_capability_map.png",15,8)

p4 <- ggplot(endpoint_results,aes(endpoint,HR))+geom_hline(yintercept=1,linetype="dashed",color=subtxt)+geom_errorbar(aes(ymin=conf.low,ymax=conf.high),width=.15,color=cyan)+geom_point(size=4,color=gold)+coord_flip()+scale_y_log10()+labs(title="Common endpoints can be re-derived consistently",subtitle="Integrated patient-level data support harmonized 3-year and 5-year analyses",x=NULL,y="Hazard ratio, log scale")+theme_dark()
save_fig(p4,"figure_04_rederived_endpoints.png",13,7)

auc_labels <- tibble(endpoint=c("3-year endpoint","5-year endpoint"),auc=c(auc_3y,auc_5y),x=.08,y=c(.88,.78),label=paste0("AUC = ",round(c(auc_3y,auc_5y),3)))
p5 <- ggplot(calibration_df,aes(mean_predicted,observed_rate))+geom_abline(slope=1,intercept=0,linetype="dashed",color=subtxt)+geom_line(color=cyan,linewidth=1.1)+geom_point(aes(size=n),color=gold)+geom_text(data=auc_labels,aes(x=x,y=y,label=label),inherit.aes=FALSE,color=txt,hjust=0,size=5)+facet_wrap(~endpoint)+scale_x_continuous(labels=percent_format(),limits=c(0,1))+scale_y_continuous(labels=percent_format(),limits=c(0,1))+labs(title="Integrated patient-level data enable prediction",subtitle="Calibration for re-derived 3-year and 5-year endpoints",x="Mean predicted risk",y="Observed event rate",size="Patients")+theme_dark()
save_fig(p5,"figure_05_prediction_calibration.png",15,8)

cat("Legacy analysis completed.\n")
print(forest_meta_df)
print(endpoint_results)
cat("AUC 3Y:",round(auc_3y,3),"\nAUC 5Y:",round(auc_5y,3),"\n")