Coronary Anomalies Survival Analysis

Author

Sarah Urbut

Published

December 11, 2024

Overview

This analysis compares survival outcomes in patients with right coronary anomalies. We use age as the time scale and survival as the primary outcome, measuring time from birth to avoid lead time bias. The analysis compares patients with and without intervention, incorporating matching by age and sex for sensitivity analysis.

Initial Survival Analysis

We define survival time as the minimum of today’s date and last known alive date, and calculate age at the end of follow-up. We then categorize patients as cases (with intervention) or controls (without intervention) based on the intervention status.

We perform matching as described below, matching cases and controls first 1:1 based on age, sex, htn, dm, smoking and hyperlipidemia, then in 2:1 matching, and with skeptical, flat, and optimistic priors.

Code

# Load required libraries
# Load required libraries
library(data.table)
library(lubridate)
library(MatchIt)
library(survival)
library(survminer)

# 1. Load and prepare case data
#df <- fread("~/Dropbox (Personal)/jack_data.csv")
df <- fread("~/Dropbox (Personal)/arca_jack_2.csv")
# First, let's check the actual column names in df
print("Column names in case data:")

[1] "Column names in case data:"

Code

print(names(df))

 [1] "Record ID"                                            
 [2] "Last name"                                            
 [3] "First name"                                           
 [4] "MRN"                                                  
 [5] "DOB"                                                  
 [6] "Sex"                                                  
 [7] "Comorbidities (choice=DM)"                            
 [8] "Comorbidities (choice=Smoker (ever))"                 
 [9] "Comorbidities (choice=Hyperlipidemia)"                
[10] "Comorbidities (choice=Hypertension)"                  
[11] "Comorbidities (choice=Stroke)"                        
[12] "Comorbidities (choice=Atrial arrhythmias (ever))"     
[13] "Comorbidities (choice=Ventricular arrhythmias (ever))"
[14] "Comorbidities (choice=Others)"                        
[15] "Comorbidities (choice=None)"                          
[16] "Comorbidities (choice=Unknown)"                       
[17] "Date of diagnosis"                                    
[18] "Age at diagnosis"                                     
[19] "Alive?"                                               
[20] "Date of death"                                        
[21] "Alive as of ---"                                      
[22] "Intervention"                                         
[23] "Date of intervention"                                 
[24] "Age at last followup"                                 
[25] "Age at death"                                         
[26] "Survival after diagnosis"                             
[27] "Duration of followup"

Code

# Adjust the column names based on your actual data
cases <- df[, .(
  hyperlip=ifelse(`Comorbidities (choice=Hyperlipidemia)`=="Checked",1,0),
    htn = ifelse(`Comorbidities (choice=Hypertension)`=="Checked",1,0),
    dm = ifelse(`Comorbidities (choice=DM)`=="Checked",1,0),
    smoking = ifelse(`Comorbidities (choice=Smoker (ever))`=="Checked",1,0),
    MRN = MRN,
    sexnum = ifelse(Sex == "Male", 1, 0),
    DOB = as.Date(DOB, format = "%Y-%m-%d"),
    Date_of_death = as.Date(`Date of death`),  # Adjust if column name is different
    Alive_as_of = as.Date(`Alive as of ---`),      # Adjust if column name is different
    alive = ifelse(`Alive?`=="Yes",1,0),
    Int = ifelse(Intervention == "Yes", 1, 0),
    case_control = 1  # Mark as cases
)]

# Set today's date explicitly
today_date <- as.Date("2024-11-20")

# Calculate end of follow-up and survival time for cases
cases[, `:=`(
    end_followup = fifelse(!is.na(Date_of_death), 
                          Date_of_death, 
                          pmin(Alive_as_of, today_date, na.rm = TRUE)),
    survival_time = as.numeric(difftime(
        fifelse(!is.na(Date_of_death),
               Date_of_death,
               pmin(Alive_as_of, today_date, na.rm = TRUE)),
        DOB, 
        units = "days"
    ))/365.25
)]

# 2. Load and prepare control data (MGBB data)
controls <- fread("~/Downloads/MGBB_Phenos_2CODE 2022-06-21.csv")
mgbb_pce=fread("~/Dropbox (Personal)/PRS_ethnicity/mgbb_data/mgbb_pce_collection_date.csv")
controls=merge(controls,mgbb_pce,by="EMPI")
# Check control data column names
print("Column names in control data:")

[1] "Column names in control data:"

Code

print(names(controls))

  [1] "EMPI"                           "Biobank_Subject_ID"            
  [3] "Sex"                            "Age_Genotyping"                
  [5] "Date_of_Birth"                  "Vital_Status"                  
  [7] "Date_of_Death"                  "Collect_Date"                  
  [9] "Last_Enc_Date"                  "Censor_Date"                   
 [11] "Time_Censor_Date"               "htn_Date"                      
 [13] "Incd_htn"                       "Prev_htn"                      
 [15] "Has_htn"                        "Time_Censor_htn"               
 [17] "dm_Date"                        "Incd_dm"                       
 [19] "Prev_dm"                        "Has_dm"                        
 [21] "Time_Censor_dm"                 "cvaTia_stroke_Date"            
 [23] "Incd_cvaTia_stroke"             "Prev_cvaTia_stroke"            
 [25] "Has_cvaTia_stroke"              "Time_Censor_cvaTia_stroke"     
 [27] "cad_Date"                       "Incd_cad"                      
 [29] "Prev_cad"                       "Has_cad"                       
 [31] "Time_Censor_cad"                "mi_Date"                       
 [33] "Incd_mi"                        "Prev_mi"                       
 [35] "Has_mi"                         "Time_Censor_mi"                
 [37] "hyperlipidemia_Date"            "Incd_hyperlipidemia"           
 [39] "Prev_hyperlipidemia"            "Has_hyperlipidemia"            
 [41] "Time_Censor_hyperlipidemia"     "cerebralAthero_Date"           
 [43] "Incd_cerebralAthero"            "Prev_cerebralAthero"           
 [45] "Has_cerebralAthero"             "Time_Censor_cerebralAthero"    
 [47] "pad_Date"                       "Incd_pad"                      
 [49] "Prev_pad"                       "Has_pad"                       
 [51] "Time_Censor_pad"                "ischemic_stroke_Date"          
 [53] "Incd_ischemic_stroke"           "Prev_ischemic_stroke"          
 [55] "Has_ischemic_stroke"            "Time_Censor_ischemic_stroke"   
 [57] "hemorrhagic_stroke_Date"        "Incd_hemorrhagic_stroke"       
 [59] "Prev_hemorrhagic_stroke"        "Has_hemorrhagic_stroke"        
 [61] "Time_Censor_hemorrhagic_stroke" "nontraumatic_SAH_Date"         
 [63] "Incd_nontraumatic_SAH"          "Prev_nontraumatic_SAH"         
 [65] "Has_nontraumatic_SAH"           "Time_Censor_nontraumatic_SAH"  
 [67] "TIA_Date"                       "Incd_TIA"                      
 [69] "Prev_TIA"                       "Has_TIA"                       
 [71] "Time_Censor_TIA"                "cerebral_aneurysm_Date"        
 [73] "Incd_cerebral_aneurysm"         "Prev_cerebral_aneurysm"        
 [75] "Has_cerebral_aneurysm"          "Time_Censor_cerebral_aneurysm" 
 [77] "all_stroke_noTIA_Date"          "Incd_all_stroke_noTIA"         
 [79] "Prev_all_stroke_noTIA"          "Has_all_stroke_noTIA"          
 [81] "Time_Censor_all_stroke_noTIA"   "mgbb_collection_date"          
 [83] "prior_ascvd"                    "enrollment_age"                
 [85] "sex"                            "sex_imputed"                   
 [87] "race"                           "race_imputed"                  
 [89] "tc"                             "tc_imputed"                    
 [91] "hdlc"                           "hdlc_imputed"                  
 [93] "sbp"                            "sbp_imputed"                   
 [95] "current_smoking"                "current_smoking_imputed"       
 [97] "dm"                             "antihtn"                       
 [99] "lipidmed"                       "pce_10yr_risk"                 
[101] "pce_10yr_risk_imputed"

Code

controls <- controls[!which(controls$Vital_Status=="Deceased" & is.na(controls$Date_of_Death)), .(
    MRN = Biobank_Subject_ID,
    sexnum = ifelse(Sex == "Male", 1, 0),
    htn=Has_htn,
    hyperlip=Has_hyperlipidemia,
    dm=Has_dm,
    smoking=ifelse(current_smoking_imputed=="yes",1,0),
    DOB = as.Date(Date_of_Birth),
    Date_of_death = as.Date(Date_of_Death),
    Censor_Date = as.Date(Censor_Date),
    alive = ifelse(Vital_Status == "Living", 1, 0),
    Int = 0,  # Controls have no intervention
    case_control = 0  # Mark as controls
)]

# Calculate survival time for controls
controls[, `:=`(
    end_followup = fifelse(!is.na(Date_of_death),
                          Date_of_death,
                          Censor_Date),
    survival_time = as.numeric(difftime(
        fifelse(!is.na(Date_of_death),
               Date_of_death,
               Censor_Date),
        DOB,
        units = "days"
    ))/365.25
)]

# 3. Combine cases and controls
combined_df <- rbind(cases, controls, fill = TRUE)

# 4. Add age calculation as of a specific date
combined_df[, age.today := as.numeric(
    difftime(today_date, DOB, units = "days")/365.25
)]

# 5. Add death indicator
combined_df[, death := ifelse(alive == 0, 1, 0)]

# Verify the data
print("Summary of survival times:")

[1] "Summary of survival times:"

Code

print(summary(combined_df$survival_time))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.856  43.822  60.033  57.553  71.288 106.913

Code

print("Distribution of cases/controls:")

[1] "Distribution of cases/controls:"

Code

print(table(combined_df$case_control))


    0     1 
63521   291

Code

print("Distribution of deaths:")

[1] "Distribution of deaths:"

Code

print(table(combined_df$death))


    0     1 
59443  4369

Code

print("Distribution of intervention by case/control status:")

[1] "Distribution of intervention by case/control status:"

Code

print(table(combined_df$case_control, combined_df$Int))

   
        0     1
  0 63521     0
  1   207    84

Code

# Approach 1: Direct Case-Control Matching
perform_direct_matching <- function(combined_df) {
    # Match cases to controls based on age and sex
    m.out <- matchit(case_control ~ age.today + sexnum+dm+htn+smoking+hyperlip, 
                    data = combined_df,
                    distance = "glm",
                    replace = TRUE)
    
    matched_data <- match.data(m.out)
   
    return(matched_data)
}

## keep more untreated among the treated 
# Plot survival curves
plot_survival <- function(matched_data) {
    fit <- survfit(Surv(time = survival_time, event = death) ~ case_control, 
                  data = matched_data)
    
    ggsurvplot(
        fit,
        data = matched_data,
        size = 1,
        conf.int = TRUE,
        pval = TRUE,
        risk.table = TRUE,
        risk.table.col = "strata",
        risk.table.height = 0.25,
        xlab="Age (years)",
        ggtheme = theme_bw()
    
    )
}

Direct Case-Control Matching (Approach 1): - Matches cases and controls directly Preserves the original intervention distribution in cases. Uses Sex, Smoking, Hypertension, and Hyperlipidemia.

Code

# For Approach 1
matched_direct <- perform_direct_matching(combined_df)
plot_survival(matched_direct)

using more controls:

Code

perform_direct_matching <- function(combined_df, ratio = 2) {
    # Match cases to controls based on age and sex
    m.out <- matchit(case_control ~ age.today + sexnum + dm + htn + smoking+hyperlip, 
                    data = combined_df,
                    distance = "glm",
                    ratio = ratio,     # Specify number of controls per case
                    replace = FALSE,   # Sample without replacement
                    method = "nearest" # Use nearest neighbor matching
    )
    
    matched_data <- match.data(m.out)
    
    # Print matching summary
    print(summary(m.out))
   
    return(matched_data)
}

# Try with 4 controls per case
matched_direct <- perform_direct_matching(combined_df, ratio = 2)


Call:
matchit(formula = case_control ~ age.today + sexnum + dm + htn + 
    smoking + hyperlip, data = combined_df, method = "nearest", 
    distance = "glm", replace = FALSE, ratio = ratio)

Summary of Balance for All Data:
          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance         0.0267        0.0045          0.5668    19.4105    0.2616
age.today       67.6184       61.7963          0.3303     1.0148    0.0771
sexnum           0.6289        0.4448          0.3811          .    0.1841
dm               0.1924        0.1007          0.2328          .    0.0918
htn              0.5430        0.3015          0.4846          .    0.2414
smoking          0.3058        0.0286          0.6017          .    0.2772
hyperlip         0.4605        0.2039          0.5148          .    0.2566
          eCDF Max
distance    0.3757
age.today   0.1196
sexnum      0.1841
dm          0.0918
htn         0.2414
smoking     0.2772
hyperlip    0.2566

Summary of Balance for Matched Data:
          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance         0.0267        0.0266          0.0006     1.0044    0.0000
age.today       67.6184       67.8270         -0.0118     1.3566    0.0320
sexnum           0.6289        0.6306         -0.0036          .    0.0017
dm               0.1924        0.2148         -0.0567          .    0.0223
htn              0.5430        0.5481         -0.0103          .    0.0052
smoking          0.3058        0.3076         -0.0037          .    0.0017
hyperlip         0.4605        0.4502          0.0207          .    0.0103
          eCDF Max Std. Pair Dist.
distance    0.0086          0.0019
age.today   0.0739          0.6425
sexnum      0.0017          0.6651
dm          0.0223          0.7627
htn         0.0052          0.6243
smoking     0.0017          0.0037
hyperlip    0.0103          0.2413

Sample Sizes:
          Control Treated
All         63521     291
Matched       582     291
Unmatched   62939       0
Discarded       0       0

Code

# Plot survival curves with more controls
plot_survival <- function(matched_data) {
    fit <- survfit(Surv(time = survival_time, event = death) ~ case_control, 
                  data = matched_data)
    
    ggsurvplot(
        fit,
        data = matched_data,
        size = 1,
       
        conf.int = TRUE,
        pval = TRUE,
        risk.table = TRUE,
        risk.table.col = "strata",
        risk.table.height = 0.25,
        ggtheme = theme_bw(),
        xlab="Age (years)",
        title = "Survival Analysis with 2:1 Control Matching"
    )
}

plot_survival(matched_direct)

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.3)’
using SDK: ‘MacOSX15.0.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000183 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.83 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.504 seconds (Warm-up)
Chain 1:                1.459 seconds (Sampling)
Chain 1:                2.963 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 9.6e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.96 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.484 seconds (Warm-up)
Chain 2:                1.438 seconds (Sampling)
Chain 2:                2.922 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 9.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.532 seconds (Warm-up)
Chain 3:                1.655 seconds (Sampling)
Chain 3:                3.187 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 9.4e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.94 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.499 seconds (Warm-up)
Chain 4:                1.383 seconds (Sampling)
Chain 4:                2.882 seconds (Total)
Chain 4:

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.3)’
using SDK: ‘MacOSX15.0.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000188 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.88 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.595 seconds (Warm-up)
Chain 1:                1.638 seconds (Sampling)
Chain 1:                3.233 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 9.7e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.97 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.592 seconds (Warm-up)
Chain 2:                1.495 seconds (Sampling)
Chain 2:                3.087 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000103 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.635 seconds (Warm-up)
Chain 3:                1.442 seconds (Sampling)
Chain 3:                3.077 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 9.5e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.95 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.634 seconds (Warm-up)
Chain 4:                1.692 seconds (Sampling)
Chain 4:                3.326 seconds (Total)
Chain 4:

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.3)’
using SDK: ‘MacOSX15.0.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000175 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.75 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.522 seconds (Warm-up)
Chain 1:                1.484 seconds (Sampling)
Chain 1:                3.006 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 9.6e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.96 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.494 seconds (Warm-up)
Chain 2:                1.546 seconds (Sampling)
Chain 2:                3.04 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000108 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.564 seconds (Warm-up)
Chain 3:                1.517 seconds (Sampling)
Chain 3:                3.081 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000108 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.49 seconds (Warm-up)
Chain 4:                1.758 seconds (Sampling)
Chain 4:                3.248 seconds (Total)
Chain 4:

[1] "Treatment Effect Estimates:"

           prior        mean       lower       upper prob_positive
2.5%   Skeptical -0.05661100 -0.09310435 -0.02175616       0.00125
2.5%1       Flat -0.05671999 -0.09431859 -0.02153903       0.00000
2.5%2 Optimistic -0.05581678 -0.09447028 -0.02202158       0.00125

[1] "Model Diagnostics:"

         mean       se_mean            sd          2.5%           25% 
-5.661100e-02  3.441805e-04  1.856050e-02 -9.310435e-02 -6.873806e-02 
          50%           75%         97.5%         n_eff          Rhat 
-5.622896e-02 -4.410580e-02 -2.175616e-02  2.908085e+03  1.000661e+00

For each prior (Skeptical, Flat, and Optimistic), we have several key metrics:

Mean: The average treatment effect from the posterior distribution
Skeptical: -0.057
Flat: -0.057
Optimistic: -0.056

These negative values suggest that being a case (having coronary anomaly) is associated with worse survival outcomes compared to controls.

Lower and Upper: These represent the 95% credible interval (similar to confidence intervals in frequentist statistics)
For all three priors, the intervals are entirely negative (from about -0.09 to -0.02)
The fact that these intervals don’t cross zero suggests strong evidence for a negative effect
prob_positive: The probability that the treatment effect is positive
All three priors show very low probability (~0.1%) of a positive effect
This reinforces the strong evidence for a negative effect

Key Interpretations:

The consistency across different priors (skeptical, flat, and optimistic) suggests that the data is strong enough to overcome prior beliefs
The negative effect is approximately a 5.7% increase in the hazard rate for cases vs controls
There is very strong evidence (>99.8% probability) that having a coronary anomaly is associated with worse survival