# Create nested directories if they don't exist
if (!dir.exists("data/confounding")) {
dir.create("data/confounding", recursive = TRUE)
}
# Importing and transforming RHC data
rhc_data <-
read.csv("https://hbiostat.org/data/repo/rhc.csv", header = TRUE) |>
# Generating `age_category`
transform(age_category = cut(age,
breaks = c(0, 50, 60, 70, 80, Inf),
include.lowest = TRUE, right = FALSE)) |>
# Factoring `race`
transform(race = factor(race,
levels = c('white', 'black', 'other'),
labels = c('White', 'Black', 'Other'))) |>
# Factoring `sex`
transform(sex = factor(sex, levels = c('Male', 'Female'))) |>
# Factoring `primary_diagnosis`
transform(primary_diagnosis = factor(cat1)) |>
within(levels(primary_diagnosis) <- list(ARF = 'ARF',
CHF = 'CHF',
MOSF = c('MOSF w/Malignancy', 'MOSF w/Sepsis'),
Other = c('Cirrhosis', 'Colon Cancer', 'Coma',
'COPD', 'Lung Cancer'))) |>
# Factoring `cancer_status`
transform(cancer_status = factor(ca,
levels = c('No', 'Yes', 'Metastatic'),
labels = c('None', 'Localized (Yes)', 'Metastatic'))) |>
# Generating `num_comorbidities`
{\(.x) transform(.x, num_comorbidities = rowSums(subset(.x, select =
cardiohx:amihx)))}() |>
# Factoring `rhc_status`
transform(rhc_status = factor(swang1, levels = c('No RHC', 'RHC'))) |>
# Converting `death` to logical
transform(death_status = ifelse(death == 'Yes', TRUE, FALSE)) |>
# Renaming other variables to be more meaningful
rename(dnr_status = dnr1,
aps_score = aps1,
two_month_survival_probability = surv2md1,
activities_of_daily_living = adld3p,
duke_activity_status_index = das2d3pc,
temperature = temp1,
heart_rate = hrt1,
mean_blood_pressure = meanbp1,
respiratory_rate = resp1,
white_blood_cell_count = wblc1,
pafi_ratio = pafi1,
partial_pressure_of_co2 = paco21,
blood_ph = ph1,
creatinine = crea1,
albumin = alb1,
glasgow_coma_scale = scoma1) |>
# Subsetting relevant columns
subset(select = c(age_category, sex, race, primary_diagnosis,
cancer_status, dnr_status, aps_score,
two_month_survival_probability, activities_of_daily_living,
duke_activity_status_index, temperature, heart_rate,
mean_blood_pressure, respiratory_rate, white_blood_cell_count,
pafi_ratio, partial_pressure_of_co2, blood_ph, creatinine,
albumin, glasgow_coma_scale,
rhc_status, death_status, num_comorbidities)) |>
# Subsetting complete case observations
na.omit()
# Save the transformed dataset
saveRDS(rhc_data, "data/confounding/rhc_data.rds")Exercise 2 Solution (R)
About the data
The following lab builds upon the RHC data used previously.
We begin by importing the RHC data, which has already been prepared for this assignment. What is different in this exercise is that we are changing the reference levels of RHC and DNR variables to No RHC and No, respectively.
# Loading the RHC data
rhc_data <- readRDS("data/confounding/rhc_data.rds")
# Convert 'rhc_status' to binary (1 for 'No RHC', 0 for 'RHC')
rhc_data$no_rhc_status <- ifelse(rhc_data$rhc_status == "No RHC", 1, 0)
# Convert 'dnr_status' to binary (1 for 'No', 0 for 'Yes')
rhc_data$no_dnr_status <- ifelse(rhc_data$dnr_status == "No", 1, 0)The variables included within this data extract are:
str(rhc_data)
#> 'data.frame': 1439 obs. of 26 variables:
#> $ age_category : Factor w/ 5 levels "[0,50)","[50,60)",..: 4 5 1 1 1 3 1 4 3 2 ...
#> $ sex : Factor w/ 2 levels "Male","Female": 1 2 2 1 1 1 2 1 2 1 ...
#> $ race : Factor w/ 3 levels "White","Black",..: 1 1 1 2 3 1 2 1 1 1 ...
#> $ primary_diagnosis : Factor w/ 4 levels "ARF","CHF","MOSF",..: 4 4 3 3 2 2 1 2 3 1 ...
#> $ cancer_status : Factor w/ 3 levels "None","Localized (Yes)",..: 2 1 1 1 1 1 1 1 1 2 ...
#> $ dnr_status : chr "No" "No" "No" "No" ...
#> $ aps_score : int 46 38 55 60 30 26 46 35 48 103 ...
#> $ two_month_survival_probability: num 0.641 0.665 0.672 0.777 0.889 ...
#> $ activities_of_daily_living : int 0 0 0 2 0 2 0 5 0 3 ...
#> $ duke_activity_status_index : num 23.5 17.5 13.5 21 22 13 29 11 17 12 ...
#> $ temperature : num 38.7 39.2 35.7 39 36.9 ...
#> $ heart_rate : int 124 134 110 130 122 115 104 54 65 170 ...
#> $ mean_blood_pressure : num 41 115 77 53 47 63 144 77 41 35 ...
#> $ respiratory_rate : num 10 36 40 12 20 22 8 42 20 42 ...
#> $ white_blood_cell_count : num 22.1 18 7.3 5.4 10.1 ...
#> $ pafi_ratio : num 68 184 171 390 333 ...
#> $ partial_pressure_of_co2 : num 40 68 25 31 40 40 39 42 40 35 ...
#> $ blood_ph : num 7.36 7.3 7.47 7.32 7.4 ...
#> $ creatinine : num 1.2 1.4 1.9 15 1.3 ...
#> $ albumin : num 3.5 3.1 3.5 3.4 4.5 ...
#> $ glasgow_coma_scale : int 0 0 37 9 0 0 0 0 0 9 ...
#> $ rhc_status : Factor w/ 2 levels "No RHC","RHC": 1 1 1 1 2 1 2 1 2 1 ...
#> $ death_status : logi FALSE FALSE TRUE FALSE TRUE TRUE ...
#> $ num_comorbidities : num 2 2 1 1 1 2 0 2 2 2 ...
#> $ no_rhc_status : num 1 1 1 1 0 1 0 1 0 1 ...
#> $ no_dnr_status : num 1 1 1 1 1 1 1 1 1 1 ...
#> - attr(*, "na.action")= 'omit' Named int [1:4296] 2 3 4 5 7 8 9 10 11 12 ...
#> ..- attr(*, "names")= chr [1:4296] "2" "3" "4" "5" ...For information regarding the variables included in the data, review the variable descriptions here and here.
The Duke Activity Status Index (DASI) is a self-administered questionnaire that measures a patient’s functional capacity or physical activity level. It is commonly used in cardiology to assess a patient’s ability to perform daily activities and is predictive of cardiovascular outcomes. The DASI score ranges from 0 to 58.2, with higher scores indicating better physical functioning and higher activity levels. The DASI score in this dataset ranges from 11 to 33. The median value of 19 and the mean value of 20.36 suggest that the typical person in this dataset has a moderate activity level. In this exercise, we will only work with the subset of the data where the DASI score is greater than 20.
About the assignment
Effect Modification refers to how the effect of an exposure on an outcome differs across levels of a third variable (the effect modifier). It focuses on how the association between the exposure and the outcome changes based on this modifying factor. This lab adapts exercises from the following tutorial on “Interaction in Epidemiology”.
For Problem Set #1, we will explore in the high_dasi_data dataset whether the association between right-heart-catheterization (no_rhc_status; treatment or exposure) on mortality (death_status; outcome) is modified by whether a participant signed a do-not-resuscitate order (no_dnr_status; potential effect modifier).
Note that we will assume that confounders for RHC and Mortality are the following:
- Age Category (
age_category) - APS Score (
aps_score) - Cancer Status (
cancer_status) - Number of Comorbidities (
num_comorbidities) - Glasgow Coma Scale (
glasgow_coma_scale) - Heart Rate (
heart_rate) - Mean Blood Pressure (
mean_blood_pressure) - Primary Diagnosis (
primary_diagnosis) - Sex (
sex) - Respiratory Rate (
respiratory_rate) - White Blood Cell Count (
white_blood_cell_count)
Interaction, particularly in the statistical context, often refers to a departure from multiplicativity when using models like logistic regression. It examines how the combined effect of two exposures differs from what would be expected if their effects were simply multiplied. Interaction is typically tested on a multiplicative scale (e.g., odds ratios) and is interpreted through synergy (super-multiplicative) or antagonism (sub-multiplicative).
For Problem Set #2, we will explore in the high_dasi_data dataset whether binary DNR status (“Do Not Resuscitate”: no_dnr_status) interacts with right-heart-catheterization (no_rhc_status) in its association with mortality (death_status).
Note that we will assume that confounders for DNR Status and Mortality are the following:
- Age Category (
age_category) - APS Score (
aps_score) - Cancer Status (
cancer_status) - Number of Comorbidities (
num_comorbidities) - Glasgow Coma Scale (
glasgow_coma_scale) - Heart Rate (
heart_rate) - Mean Blood Pressure (
mean_blood_pressure) - Primary Diagnosis (
primary_diagnosis) - White Blood Cell Count (
white_blood_cell_count) - Activities of Daily Living (
activities_of_daily_living)
Problem #1: Exploring Effect Modification
1(a) Fit the effect modification model
Generate a multivariable logistic regression model (effect_modification_model), adhering to the following specifications:
- Set the outcome as death.
- Set the exposure as no RHC.
- Include a multiplicative interaction term of the exposure by no DNR status.
- Adjust for the stated confounders in the relationship (no RHC vs death) in the model.
# Logistic regression with interaction between no_rhc_status and no_dnr_status in the subset data
effect_modification_model <- glm(death_status ~ no_rhc_status * no_dnr_status +
age_category + aps_score + cancer_status +
num_comorbidities + glasgow_coma_scale +
heart_rate + mean_blood_pressure +
primary_diagnosis + sex +
respiratory_rate + white_blood_cell_count,
data = high_dasi_data,
family = binomial)
summary(effect_modification_model)
#>
#> Call:
#> glm(formula = death_status ~ no_rhc_status * no_dnr_status +
#> age_category + aps_score + cancer_status + num_comorbidities +
#> glasgow_coma_scale + heart_rate + mean_blood_pressure + primary_diagnosis +
#> sex + respiratory_rate + white_blood_cell_count, family = binomial,
#> data = high_dasi_data)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -2.3712376 1.4019207 -1.691 0.090756 .
#> no_rhc_status 3.4969691 1.6575046 2.110 0.034877 *
#> no_dnr_status 0.6228925 1.3008762 0.479 0.632063
#> age_category[50,60) 0.5615264 0.3007867 1.867 0.061921 .
#> age_category[60,70) 1.0150295 0.2615213 3.881 0.000104 ***
#> age_category[70,80) 0.6505600 0.2780222 2.340 0.019286 *
#> age_category[80,Inf] 0.8609830 0.3498095 2.461 0.013844 *
#> aps_score -0.0032209 0.0070546 -0.457 0.647977
#> cancer_statusLocalized (Yes) 0.6516821 0.2753557 2.367 0.017948 *
#> cancer_statusMetastatic 1.8068081 0.4031984 4.481 7.42e-06 ***
#> num_comorbidities 0.2526729 0.0888800 2.843 0.004471 **
#> glasgow_coma_scale -0.0006699 0.0053385 -0.125 0.900136
#> heart_rate -0.0030546 0.0026816 -1.139 0.254659
#> mean_blood_pressure -0.0049564 0.0027455 -1.805 0.071034 .
#> primary_diagnosisCHF 0.7685691 0.3003703 2.559 0.010505 *
#> primary_diagnosisMOSF 0.4425051 0.2653825 1.667 0.095430 .
#> primary_diagnosisOther 0.6229054 0.2796315 2.228 0.025908 *
#> sexFemale -0.1207658 0.1982125 -0.609 0.542343
#> respiratory_rate 0.0119828 0.0085488 1.402 0.161005
#> white_blood_cell_count 0.0153771 0.0081258 1.892 0.058442 .
#> no_rhc_status:no_dnr_status -3.2992008 1.6709649 -1.974 0.048333 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 839.22 on 622 degrees of freedom
#> Residual deviance: 712.15 on 602 degrees of freedom
#> AIC: 754.15
#>
#> Number of Fisher Scoring iterations: 51(b) Calculate Relative Excess Risk due to Interaction (RERI)
Effect modification looks at how the effect of one exposure on an outcome changes based on the level of another exposure or factor (the effect modifier). RERI is useful in this context because it quantifies the additive interaction (whether the combined effect is larger or smaller than expected based on individual effects).
# Extract the coefficients
beta_rhc <- as.numeric(effect_modification_model$coefficients['no_rhc_status']) # Coefficient for no_rhc_status
beta_dnr <- as.numeric(effect_modification_model$coefficients['no_dnr_status']) # Coefficient for no_dnr_status
beta_interaction <- as.numeric(effect_modification_model$coefficients['no_rhc_status:no_dnr_status']) # Interaction term
# Calculate the Odds Ratios
OR_rhc_only <- exp(beta_rhc) # OR when no_rhc_status = 1, no_dnr_status = 0
OR_dnr_only <- exp(beta_dnr) # OR when no_rhc_status = 0, no_dnr_status = 1
OR_both <- exp(beta_rhc + beta_dnr + beta_interaction)
# OR when both no_rhc_status and no_dnr_status are 1
# Calculate RERI (Relative Excess Risk due to Interaction)
RERI <- OR_both - OR_rhc_only - OR_dnr_only + 1
# Output the results
cat("OR for no RHC only (no_rhc_status = 1, no_dnr_status = 0):", OR_rhc_only, "\n")
#> OR for no RHC only (no_rhc_status = 1, no_dnr_status = 0): 33.01524
cat("OR for no DNR only (no_rhc_status = 0, no_dnr_status = 1):", OR_dnr_only, "\n")
#> OR for no DNR only (no_rhc_status = 0, no_dnr_status = 1): 1.864313
cat("OR for both no RHC and no DNR (no_rhc_status = 1, no_dnr_status = 1):", OR_both, "\n")
#> OR for both no RHC and no DNR (no_rhc_status = 1, no_dnr_status = 1): 2.272001
cat("Relative Excess Risk due to Interaction (RERI):", RERI, "\n")
#> Relative Excess Risk due to Interaction (RERI): -31.60755Interpretation of RERI in Effect Modification:
- RERI > 0: Indicates a positive additive interaction, meaning the combined effect of two exposures (e.g., no RHC and no DNR) is greater than the sum of their individual effects. In other words, there is an excess risk attributable to the interaction between the two exposures.
- RERI = 0: Indicates no additive interaction. The combined effect of the two exposures is equal to the sum of their individual effects. There is no additional risk due to their interaction.
- RERI < 0: Indicates a negative additive interaction (also known as antagonism), meaning the combined effect of the two exposures is less than the sum of their individual effects. This suggests that the exposures may be protective when combined, compared to when considered separately.
1(c) Exploring Effect Modification using interactionR [15%] and Publish (optional)
Based on the output of the above regression, implement the following:
- Use
interactionR::interactionR()to obtain a complete reporting of effect modification estimates.
# Load the necessary package
library(interactionR)
# Obtain the interaction effect modification estimates
interaction_result <- interactionR(effect_modification_model,
exposure_names = c("no_rhc_status", "no_dnr_status"),
ci.type = "delta", # Confidence interval method
em = TRUE, # Effect modification estimates
recode = FALSE) # Do not recode the exposure variables
# Display the complete effect modification report
interaction_result$dframe- Use
Publish::publish()to obtain stratum-specific estimates of the main effect. NOTE: Only return those quantities requested (i.e., not the full regression table).
# Load the necessary package
library(Publish)
# Recode variables for better interpretation
high_dasi_data$no_dnr_factor <- factor(high_dasi_data$no_dnr_status, levels = c(1, 0), labels = c("No", "Yes"))
high_dasi_data$no_dnr_factor <- relevel(high_dasi_data$no_dnr_factor, ref = "Yes")
high_dasi_data$no_rhc_factor <- factor(high_dasi_data$no_rhc_status, levels = c(1, 0), labels = c("No RHC", "RHC"))
high_dasi_data$no_rhc_factor <- relevel(high_dasi_data$no_rhc_factor, ref = "RHC")
# Fit the model with recoded factors
effect_modification_model2 <- glm(death_status ~ no_rhc_factor * no_dnr_factor +
age_category + aps_score + cancer_status +
num_comorbidities + glasgow_coma_scale +
heart_rate + mean_blood_pressure +
primary_diagnosis + sex +
respiratory_rate + white_blood_cell_count,
data = high_dasi_data,
family = binomial)
# Obtain stratum-specific estimates from the model
stratum_estimates <- publish(effect_modification_model2)
#> Variable Units OddsRatio CI.95 p-value
#> age_category [0,50) Ref
#> [50,60) 1.75 [0.97;3.16] 0.0619213
#> [60,70) 2.76 [1.65;4.61] 0.0001039
#> [70,80) 1.92 [1.11;3.31] 0.0192860
#> [80,Inf] 2.37 [1.19;4.70] 0.0138438
#> aps_score 1.00 [0.98;1.01] 0.6479771
#> cancer_status None Ref
#> Localized (Yes) 1.92 [1.12;3.29] 0.0179479
#> Metastatic 6.09 [2.76;13.42] < 1e-04
#> num_comorbidities 1.29 [1.08;1.53] 0.0044712
#> glasgow_coma_scale 1.00 [0.99;1.01] 0.9001364
#> heart_rate 1.00 [0.99;1.00] 0.2546595
#> mean_blood_pressure 1.00 [0.99;1.00] 0.0710337
#> primary_diagnosis ARF Ref
#> CHF 2.16 [1.20;3.89] 0.0105053
#> MOSF 1.56 [0.93;2.62] 0.0954302
#> Other 1.86 [1.08;3.23] 0.0259076
#> sex Male Ref
#> Female 0.89 [0.60;1.31] 0.5423427
#> respiratory_rate 1.01 [1.00;1.03] 0.1610054
#> white_blood_cell_count 1.02 [1.00;1.03] 0.0584416
#> no_rhc_factor(RHC): no_dnr_factor(No vs Yes) 1.86 [0.15;23.87] 0.6320628
#> no_rhc_factor(No RHC): no_dnr_factor(No vs Yes) 0.07 [0.01;0.55] 0.0113453
#> no_dnr_factor(Yes): no_rhc_factor(No RHC vs RHC) 33.02 [1.28;850.32] 0.0348774
#> no_dnr_factor(No): no_rhc_factor(No RHC vs RHC) 1.22 [0.79;1.87] 0.3682362
# Print only the relevant stratum-specific estimates (for no_rhc_status and no_dnr_status)
tail(stratum_estimates$rawTable,2)Note: Effect modification occurs when the effect of an exposure on an outcome differs across levels of a third variable (the effect modifier). It does not require the interaction term (no_rhc_status:no_dnr_status) to be statistically significant; rather, we are looking at whether the association between the primary exposure (RHC) and the outcome (death) changes depending on the level of the effect modifier (DNR).
Problem #2: Exploring Interaction
2(a) Interaction Model
Generate a multivariable logistic regression model, adhering to the following specifications:
- Set the outcome as death.
- Set two exposures as RHC and DNR status.
- Include a multiplicative interaction term of the exposures.
- Adjust for the stated confounders in the relationship (no RHC vs death, as well as no DNR vs death) in the model. That means the union of the two confounder sets.
# Logistic regression with interaction between no_rhc_status and no_dnr_status in the subset data
interaction_model <- glm(death_status ~ no_rhc_status * no_dnr_status +
age_category + aps_score + cancer_status +
num_comorbidities + glasgow_coma_scale +
heart_rate + mean_blood_pressure +
primary_diagnosis + sex +
respiratory_rate + white_blood_cell_count +
activities_of_daily_living,
data = high_dasi_data,
family = binomial)
summary(interaction_model)
#>
#> Call:
#> glm(formula = death_status ~ no_rhc_status * no_dnr_status +
#> age_category + aps_score + cancer_status + num_comorbidities +
#> glasgow_coma_scale + heart_rate + mean_blood_pressure + primary_diagnosis +
#> sex + respiratory_rate + white_blood_cell_count + activities_of_daily_living,
#> family = binomial, data = high_dasi_data)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -2.2800128 1.4080440 -1.619 0.105388
#> no_rhc_status 3.3958648 1.6658528 2.039 0.041499 *
#> no_dnr_status 0.5238144 1.3076850 0.401 0.688740
#> age_category[50,60) 0.5760202 0.3013259 1.912 0.055925 .
#> age_category[60,70) 1.0173738 0.2632735 3.864 0.000111 ***
#> age_category[70,80) 0.6498421 0.2793077 2.327 0.019986 *
#> age_category[80,Inf] 0.8040125 0.3523068 2.282 0.022481 *
#> aps_score -0.0037222 0.0070865 -0.525 0.599411
#> cancer_statusLocalized (Yes) 0.6443974 0.2770866 2.326 0.020039 *
#> cancer_statusMetastatic 1.8888507 0.4064679 4.647 3.37e-06 ***
#> num_comorbidities 0.2408454 0.0895108 2.691 0.007131 **
#> glasgow_coma_scale -0.0008445 0.0053673 -0.157 0.874971
#> heart_rate -0.0032735 0.0026971 -1.214 0.224856
#> mean_blood_pressure -0.0050871 0.0027627 -1.841 0.065575 .
#> primary_diagnosisCHF 0.7932514 0.3007341 2.638 0.008347 **
#> primary_diagnosisMOSF 0.3831509 0.2685658 1.427 0.153679
#> primary_diagnosisOther 0.5884193 0.2812664 2.092 0.036435 *
#> sexFemale -0.1391314 0.1993885 -0.698 0.485308
#> respiratory_rate 0.0123729 0.0086003 1.439 0.150246
#> white_blood_cell_count 0.0159359 0.0081423 1.957 0.050326 .
#> activities_of_daily_living 0.2249761 0.1006831 2.234 0.025450 *
#> no_rhc_status:no_dnr_status -3.2074335 1.6793158 -1.910 0.056138 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 839.22 on 622 degrees of freedom
#> Residual deviance: 707.11 on 601 degrees of freedom
#> AIC: 751.11
#>
#> Number of Fisher Scoring iterations: 5Note: Check if the interaction term (no_rhc_status:no_dnr_status) has a low p-value. If this p-value is below the acceptable cutpoint (commonly used significance threshold of 0.05 or something reasonable/justifiable), we could conclude that the interaction between RHC and DNR status on the outcome (death) is statistically significant. This could suggest that the association between RHC status and death depends on DNR status.
2(b) Exploring Interaction using interactionR
Based on the output of the above regression, implement the following:
- Use
interactionR::interactionR()to obtain a complete reporting of interaction estimates.
# Load the necessary package
library(interactionR)
# Obtain the interaction estimates
interaction_result <- interactionR(interaction_model,
exposure_names = c("no_rhc_status", "no_dnr_status"),
ci.type = "delta", # Confidence interval method
em = FALSE, # Interaction estimates
recode = FALSE) # Do not recode the exposure variables
# Display the complete interaction report
interaction_result$dframe2(c) Calculate RERI, AP, and SI
Based on the output of the above regression, calculate RERI, the attributable proportion due to interaction (AP), and the synergy index (SI).
# Extract the coefficients
beta_rhc <- coef(interaction_model)["no_rhc_status"]
beta_dnr <- coef(interaction_model)["no_dnr_status"]
beta_interaction <- coef(interaction_model)["no_rhc_status:no_dnr_status"]
# Calculate the odds ratios
OR_rhc_only <- exp(beta_rhc) # OR when no_rhc_status = 1, no_dnr_status = 0
OR_dnr_only <- exp(beta_dnr) # OR when no_rhc_status = 0, no_dnr_status = 1
OR_both <- exp(beta_rhc + beta_dnr + beta_interaction)
# OR when both no_rhc_status and no_dnr_status are 1
# Calculate RERI
RERI <- OR_both - OR_rhc_only - OR_dnr_only + 1
# Calculate AP
AP <- RERI / OR_both
# Calculate SI
SI <- (OR_both - 1) / ((OR_rhc_only - 1) + (OR_dnr_only - 1))
# Output the results
cat("Relative Excess Risk due to Interaction (RERI):", RERI, "\n")
#> Relative Excess Risk due to Interaction (RERI): -28.49034
cat("Attributable Proportion (AP):", AP, "\n")
#> Attributable Proportion (AP): -13.97569
cat("Synergy Index (SI):", SI, "\n")
#> Synergy Index (SI): 0.03517111Interpretation of Results:
- RERI: A positive RERI indicates an excess risk due to the interaction between RHC and DNR. A negative RERI suggests antagonism.
-
AP: Represents the proportion of the risk due to the interaction.
- AP>0: A positive proportion of the risk is due to the interaction.
- AP=0: No additional risk is due to the interaction.
- AP<0: A negative interaction, indicating antagonism between the two exposures.
-
SI: A synergy index greater than 1 indicates a synergistic effect, while an SI less than 1 indicates antagonism.
- SI>1: Synergistic effect. The combined effect of the two exposures is greater than expected based on multiplicative effects.
- SI=1: No multiplicative interaction.
- SI<1: Antagonism. The combined effect of the two exposures is less than expected based on multiplicative effects.
Note: If the p-values for the interaction term, RERI, AP, and SI are all greater than a reasonable cut-point (e.g., 0.05 is one example), that indicates no statistically significant interaction. In that case, this would mean we do not have strong evidence to claim that the combined effects of RHC and DNR on mortality are significantly different from their individual effects.