diff --git a/analysis/apply_across_MI/apply_across_MI.R b/analysis/apply_across_MI/apply_across_MI.R new file mode 100644 index 0000000..a1a296c --- /dev/null +++ b/analysis/apply_across_MI/apply_across_MI.R @@ -0,0 +1,692 @@ +# ------------------------------------------------------------------------------ +# +# apply_across_MI.R +# +# This file applies multiple imputation to the BMI and Smoking covariates +# MI is conducted in "across" fashion, meaning that the result is a dataframe +# 10x larger than the origina containing all datasets +# +# Arguments: +# - cohort - string, defines which of three opensafely cohorts to describe +# (prevax, vax, unvax) +# - preex - boolean/string, defines preexisting conditions +# for the replication preex = FALSE always +# ("All", TRUE, or FALSE) +# +# Returns: +# - output/apply_across_MI/input__clean_across_MI_ami.rds +# - output/apply_across_MI/input__clean_across_MI_sahhs.rds +# +# Authors: Emma Tarmey +# +# ------------------------------------------------------------------------------ + + +# Load libraries --------------------------------------------------------------- +print("Load libraries") + +library(magrittr) +library(mice) +library(here) +library(dplyr) +library(fs) +library(survival) + + +# Source common functions ------------------------------------------------------ +print("Source common functions") + +source("analysis/utility.R") + +# only needed here, generates unique patient IDs for stacked data +# unique IDs needed for cox_ipw action +new_patient_ids <- function(ids = NULL, num_needed = NULL, shift = NULL) { + if (num_needed == 1) { + return (ids + shift) + } else { + new_ids <- c( + ids, + new_patient_ids(ids = (ids + shift), num_needed = (num_needed - 1), shift = shift) + ) + } + return (new_ids) +} + + +# Specify arguments ------------------------------------------------------------ +print("Specify arguments") + +args <- commandArgs(trailingOnly = TRUE) +print(length(args)) + +if (length(args) == 0) { + # default argument values + cohort <- "prevax" + preex <- "All" + +} else { + # YAML arguments + cohort <- args[[1]] + + # optional argument + if (length(args) < 2) { + preex <- "All" + } else { + preex <- args[[2]] + } # allow an empty input for the preex variable +} + + +# Load data -------------------------------------------------------------------- +print("Load data") + +df <- readr::read_rds(paste0( + "output/dataset_clean/input_", + cohort, + "_clean_prehoc.rds" +)) + +df <- as.data.frame(df) + + +# Define variables ------------------------------------------------------------- +print("Define variables") + +df$cov_bin_covid <- !is.na(df$exp_date_covid) +df$cov_bin_sahhs <- !is.na(df$out_date_stroke_sahhs) + + +# Check all covariates types --------------------------------------------------- +print("Check all covariate types") + +df$cov_bin_ami <- as.factor(df$cov_bin_ami) # outcome +df$cov_bin_sahhs <- as.factor(df$cov_bin_sahhs) # outcome +df$cov_bin_covid <- as.factor(df$cov_bin_covid) # exposure + +df$cov_num_age <- as.numeric(df$cov_num_age) +df$cov_cat_sex <- as.factor(df$cov_cat_sex) +df$cov_num_bmi <- as.numeric(df$cov_num_bmi) +df$cov_cat_ethnicity <- as.factor(df$cov_cat_ethnicity) +df$cov_cat_imd <- as.factor(df$cov_cat_imd) +df$cov_cat_smoking <- as.factor(df$cov_cat_smoking) + +df$cov_bin_carehome <- as.factor(df$cov_bin_carehome) +df$cov_bin_hcworker <- as.factor(df$cov_bin_hcworker) +df$cov_bin_dementia <- as.factor(df$cov_bin_dementia) +df$cov_bin_liver_disease <- as.factor(df$cov_bin_liver_disease) +df$cov_bin_ckd <- as.factor(df$cov_bin_ckd) + +df$cov_bin_cancer <- as.factor(df$cov_bin_cancer) +df$cov_bin_hypertension <- as.factor(df$cov_bin_hypertension) +df$cov_bin_diabetes <- as.factor(df$cov_bin_diabetes) +df$cov_bin_obesity <- as.factor(df$cov_bin_obesity) +df$cov_bin_copd <- as.factor(df$cov_bin_copd) + +df$cov_bin_depression <- as.factor(df$cov_bin_depression) +df$cov_bin_stroke_all <- as.factor(df$cov_bin_stroke_all) +df$cov_bin_other_ae <- as.factor(df$cov_bin_other_ae) +df$cov_bin_vte <- as.factor(df$cov_bin_vte) +df$cov_bin_hf <- as.factor(df$cov_bin_hf) + +df$cov_bin_angina <- as.factor(df$cov_bin_angina) +df$cov_bin_lipidmed <- as.factor(df$cov_bin_lipidmed) +df$cov_bin_antiplatelet <- as.factor(df$cov_bin_antiplatelet) +df$cov_bin_anticoagulant <- as.factor(df$cov_bin_anticoagulant) +df$cov_bin_cocp <- as.factor(df$cov_bin_cocp) + +df$cov_bin_hrt <- as.factor(df$cov_bin_hrt) +df$strat_cat_region <- as.factor(df$strat_cat_region) + + +# Applying multiple imputation to BMI and smoking covariates ------------------- +print("Applying multiple imputation to BMI and smoking covariates") + +# set random seed +set.seed(2026) + +# re-cast missing smoking to NA +smoking_missing <- df$cov_cat_smoking == "Missing" +df$cov_cat_smoking[smoking_missing] <- NA + +# check missingness of smoking and bmi variables +percent_smoking_missing <- signif(100 * (sum(is.na(df$cov_cat_smoking)) / length(df$cov_cat_smoking)), digits = 4) +percent_bmi_missing <- signif(100 * (sum(is.na(df$cov_num_bmi)) / length(df$cov_num_bmi)), digits = 4) + +print(paste0("The variable smoking is ", percent_smoking_missing, "% missing")) +print(paste0("The variable bmi is ", percent_bmi_missing, "% missing")) + + +# Specify imputation methods for each outcome (ami and sahhs) ------------------ +print("Specify imputation methods for each outcome (ami and sahhs)") + +# The below ensures that bmi and smoking are handled with specific +# imputation methods and that all other covariates are left alone +# See: https://www.rdocumentation.org/packages/mice/versions/3.17.0/topics/mice +imp_method <- rep("", length.out = length(colnames(df))) +names(imp_method) <- colnames(df) +imp_method["cov_cat_smoking"] <- "polyreg" # smoking is categorical with 3 levels, Polytomous logistic regression +imp_method["cov_num_bmi"] <- "norm" # bmi is numerical, Bayesian linear regression + + +# Specify imputation formulas for each outcome (ami and sahhs) ----------------- +print("Specify imputation formulas for outcome") + +# Specify imputation formulas, exclude variable such as index date +all_var_names <- c( + "cov_bin_ami", + "cov_bin_sahhs", + "cov_bin_covid", + + "cov_num_age", + "cov_cat_sex", + "cov_num_bmi", + "cov_cat_ethnicity", + "cov_cat_imd", + "cov_cat_smoking", + + "cov_bin_carehome", + "cov_bin_hcworker", + "cov_bin_dementia", + "cov_bin_liver_disease", + "cov_bin_ckd", + + "cov_bin_cancer", + "cov_bin_hypertension", + "cov_bin_diabetes", + "cov_bin_obesity", + "cov_bin_copd", + + "cov_bin_depression", + "cov_bin_stroke_all", + "cov_bin_other_ae", + "cov_bin_vte", + "cov_bin_hf", + + "cov_bin_angina", + "cov_bin_lipidmed", + "cov_bin_antiplatelet", + "cov_bin_anticoagulant", + "cov_bin_cocp", + + "cov_bin_hrt", + "strat_cat_region", + "vax_cat_jcvi_group" + # "cens_status" # excluded +) + +all_var_names_except_bmi <- c( + "cov_bin_ami", + "cov_bin_sahhs", + "cov_bin_covid", + + "cov_num_age", + "cov_cat_sex", + # "cov_num_bmi", + "cov_cat_ethnicity", + "cov_cat_imd", + "cov_cat_smoking", + + "cov_bin_carehome", + "cov_bin_hcworker", + "cov_bin_dementia", + "cov_bin_liver_disease", + "cov_bin_ckd", + + "cov_bin_cancer", + "cov_bin_hypertension", + "cov_bin_diabetes", + "cov_bin_obesity", + "cov_bin_copd", + + "cov_bin_depression", + "cov_bin_stroke_all", + "cov_bin_other_ae", + "cov_bin_vte", + "cov_bin_hf", + + "cov_bin_angina", + "cov_bin_lipidmed", + "cov_bin_antiplatelet", + "cov_bin_anticoagulant", + "cov_bin_cocp", + + "cov_bin_hrt", + "strat_cat_region", + "vax_cat_jcvi_group" + # "cens_status" # excluded +) + +all_var_names_except_smoking <- c( + "cov_bin_ami", + "cov_bin_sahhs", + "cov_bin_covid", + + "cov_num_age", + "cov_cat_sex", + "cov_num_bmi", + "cov_cat_ethnicity", + "cov_cat_imd", + # "cov_cat_smoking", + + "cov_bin_carehome", + "cov_bin_hcworker", + "cov_bin_dementia", + "cov_bin_liver_disease", + "cov_bin_ckd", + + "cov_bin_cancer", + "cov_bin_hypertension", + "cov_bin_diabetes", + "cov_bin_obesity", + "cov_bin_copd", + + "cov_bin_depression", + "cov_bin_stroke_all", + "cov_bin_other_ae", + "cov_bin_vte", + "cov_bin_hf", + + "cov_bin_angina", + "cov_bin_lipidmed", + "cov_bin_antiplatelet", + "cov_bin_anticoagulant", + "cov_bin_cocp", + + "cov_bin_hrt", + "strat_cat_region", + "vax_cat_jcvi_group" + # "cens_status" # excluded +) + +my_formulas <- list( + cov_cat_smoking = as.formula(paste0("cov_cat_smoking ~ ", paste(all_var_names_except_smoking, collapse = " + "), " + H0")), + cov_num_bmi = as.formula(paste0("cov_num_bmi ~ ", paste(all_var_names_except_bmi, collapse = " + "), " + H0")) +) + +# Calculate Nelson-Aalen Estimator for outcome ----------- +print("Calculate Nelson-Aalen Estimator for outcome") + +basehaz_with_SE <- function(fit, newdata, centered = TRUE) +{ + if (inherits(fit, "coxphms")) + stop("the basehaz function is not implemented for multi-state models") + if (!inherits(fit, "coxph")) + stop("must be a coxph object") + if (!missing(newdata)) { + sfit <- survfit(fit, newdata=newdata, se.fit=TRUE) + chaz <- sfit$cumhaz + } + else { + sfit <- survfit(fit, se.fit=TRUE) + if (!centered) { + # The right thing to do here is to call survfit with a vector of + # all zeros for the "subject to predict". But if there is a factor + # in the model, there may be no subject at all who will give all + # zeros, so we post process instead + zcoef <- ifelse(is.na(coef(fit)), 0, coef(fit)) + offset <- sum(fit$means * zcoef) + chaz <- sfit$cumhaz * exp(-offset) + } + else { + chaz <- sfit$cumhaz + } + } + + new <- data.frame( + hazard = chaz, + time = sfit$time, + std_err_cumhaz = sfit$std.err + ) + + strata <- sfit$strata + if (!is.null(strata)) { + new$strata <- factor(rep(names(strata), strata), levels = names(strata)) + } + + return (new) +} + +nelsonaalen_with_SE <- function(data, timevar, statusvar, ...) { + if (!is.data.frame(data)) { + stop("Data must be a data frame") + } + timevar <- as.character(substitute(timevar)) + statusvar <- as.character(substitute(statusvar)) + time <- data[, timevar, drop = TRUE] + status <- data[, statusvar, drop = TRUE] + + coxph_obj <- survival::coxph(survival::Surv(time, status) ~ 1, ...) + hazard <- basehaz_with_SE(coxph_obj) + + # Adjust depending on near-tie correction + idx <- if (coxph_obj$timefix) { + match(coxph_obj$y[, "time"], hazard[, "time"]) + } else match(time, hazard[, "time"]) + + nelsonaalen_estimates_with_SE <- data.frame( + nelsonaalen_estimates = hazard[idx, "hazard"], + nelsonaalen_se = hazard[idx, "std_err_cumhaz"] + ) + + return (nelsonaalen_estimates_with_SE) +} + +df_ami <- df +df_sahhs <- df + +outcome_cox_dates_ami <- rep(as.Date(NA), times = nrow(df_ami)) +cens_status_ami <- rep(NA, times = nrow(df_ami)) + +# 0 = censoring time = date of end of study +# 1 = failure time = time of outcome event +# See: https://www.rdocumentation.org/packages/survival/versions/3.8-3/topics/Surv +# and: https://glmnet.stanford.edu/articles/Coxnet.html#basic-usage-for-right-censored-data +for (i in c(1:nrow(df_ami))) { + if (is.na(df_ami$out_date_ami[i])) { + # right-hand censorship takes place + cens_status_ami[i] <- 0 + outcome_cox_dates_ami[i] <- df_ami$end_date_outcome[i] + } else { + # event takes place (failure) + cens_status_ami[i] <- 1 + outcome_cox_dates_ami[i] <- df_ami$out_date_ami[i] + } +} + +outcome_cox_dates_sahhs <- rep(as.Date(NA), times = nrow(df_sahhs)) +cens_status_sahhs <- rep(NA, times = nrow(df_sahhs)) + +# 0 = censoring time = date of end of study +# 1 = failure time = time of outcome event +# See: https://www.rdocumentation.org/packages/survival/versions/3.8-3/topics/Surv +# and: https://glmnet.stanford.edu/articles/Coxnet.html#basic-usage-for-right-censored-data +for (i in c(1:nrow(df_sahhs))) { + if (is.na(df_sahhs$out_date_stroke_sahhs[i])) { + # right-hand censorship takes place + cens_status_sahhs[i] <- 0 + outcome_cox_dates_sahhs[i] <- df_sahhs$end_date_outcome[i] + } else { + # event takes place (failure) + cens_status_sahhs[i] <- 1 + outcome_cox_dates_sahhs[i] <- df_sahhs$out_date_stroke_sahhs[i] + } +} + +# add data to dataframes +df_ami$outcome_cox_dates_ami <- as.numeric(outcome_cox_dates_ami) +df_ami$cens_status_ami <- cens_status_ami + +df_sahhs$outcome_cox_dates_sahhs <- as.numeric(outcome_cox_dates_sahhs) +df_sahhs$cens_status_sahhs <- cens_status_sahhs + +# ami +df_ami_nelsonaalen <- data.frame(time = df_ami$outcome_cox_dates_ami, status = df_ami$cens_status_ami) +H0_ami <- nelsonaalen_with_SE(df_ami_nelsonaalen, time, status) +df_ami_nelsonaalen$H0 <- H0_ami$nelsonaalen_estimates +df_ami_nelsonaalen$se <- H0_ami$nelsonaalen_se +df_ami$H0 <- H0_ami$nelsonaalen_estimates + +# sahhs +df_sahhs_nelsonaalen <- data.frame(time = df_sahhs$outcome_cox_dates_sahhs, status = df_sahhs$cens_status_sahhs) +H0_sahhs <- nelsonaalen_with_SE(df_sahhs_nelsonaalen, time, status) +df_sahhs_nelsonaalen$H0 <- H0_sahhs$nelsonaalen_estimates +df_sahhs_nelsonaalen$se <- H0_sahhs$nelsonaalen_se +df_sahhs$H0 <- H0_sahhs$nelsonaalen_estimates + + +# Applying multiple imputation to BMI and smoking covariates for outcome --- +print("Applying multiple imputation to BMI and smoking covariates for outcome") + +# Apply multiple imputation for ami outcome +imp_ami <- mice::mice( + data = df_ami, + m = get_number_of_imputed_datasets(), + maxit = 20, + formulas = my_formulas, + imp_method = unname(imp_method) +) + +df_post_imputation_ami <- mice::complete( + imp_ami, + action = "long", + include = FALSE +) + +df_post_imputation_ami <- subset( + df_post_imputation_ami, + select = -c(.imp, .id) +) + +# Apply multiple imputation for sahhs outcome +imp_sahhs <- mice::mice( + data = df_sahhs, + m = get_number_of_imputed_datasets(), + maxit = 20, + formulas = my_formulas, + imp_method = unname(imp_method) +) + +df_post_imputation_sahhs <- mice::complete( + imp_sahhs, + action = "long", + include = FALSE +) + +df_post_imputation_sahhs <- subset( + df_post_imputation_sahhs, + select = -c(.imp, .id) +) + + +# Remove now unused level 'missing' from smoking covariate ----- +print("Remove now unused level 'missing' from smoking covariate") + +df_post_imputation_ami$cov_cat_smoking <- factor( + df_post_imputation_ami$cov_cat_smoking, + levels = levels(droplevels(df_post_imputation_ami$cov_cat_smoking)) +) + +df_post_imputation_sahhs$cov_cat_smoking <- factor( + df_post_imputation_sahhs$cov_cat_smoking, + levels = levels(droplevels(df_post_imputation_sahhs$cov_cat_smoking)) +) + + +# Re-assign unique identifiers to imputed dataset for outcome ----- +print("Re-assign unique identifiers to imputed dataset for outcome ") + +patient_id_1 <- as.numeric(unique(df_post_imputation_ami$patient_id)) +shift <- max(patient_id_1) + 1 + +new_patient_id <- new_patient_ids( + ids = patient_id_1, + num_needed = get_number_of_imputed_datasets(), + shift = shift +) + +df_post_imputation_ami$patient_id <- new_patient_id +df_post_imputation_sahhs$patient_id <- new_patient_id + +# Re-convert date variables to date type +all_date_vars <- grep("date", colnames(df_post_imputation_ami)) +for (var in all_date_vars) { + df_post_imputation_ami[, var] <- as.Date(format(as.Date(df_post_imputation_ami[, var], origin = lubridate::origin), "%Y-%m-%d")) + df_post_imputation_sahhs[, var] <- as.Date(format(as.Date(df_post_imputation_sahhs[, var], origin = lubridate::origin), "%Y-%m-%d")) +} + +# check missingness of smoking and bmi variables +percent_smoking_missing_ami <- signif(100 * (sum(is.na(df_post_imputation_ami$cov_cat_smoking)) / length(df_post_imputation_ami$cov_cat_smoking)), digits = 4) +percent_bmi_missing_ami <- signif(100 * (sum(is.na(df_post_imputation_ami$cov_num_bmi)) / length(df_post_imputation_ami$cov_num_bmi)), digits = 4) + +print(paste0("The variable smoking is ", percent_smoking_missing_ami, "% missing (for ami)")) +print(paste0("The variable bmi is ", percent_bmi_missing_ami, "% missing (for ami)")) + +percent_smoking_missing_sahhs <- signif(100 * (sum(is.na(df_post_imputation_sahhs$cov_cat_smoking)) / length(df_post_imputation_sahhs$cov_cat_smoking)), digits = 4) +percent_bmi_missing_sahhs <- signif(100 * (sum(is.na(df_post_imputation_sahhs$cov_num_bmi)) / length(df_post_imputation_sahhs$cov_num_bmi)), digits = 4) + +print(paste0("The variable smoking is ", percent_smoking_missing_sahhs, "% missing (for sahhs)")) +print(paste0("The variable bmi is ", percent_bmi_missing_sahhs, "% missing (for sahhs)")) + + +# Remove unnecessary variables --- +print("Remove unnecessary_variables ") + +print(colnames(df_post_imputation_ami)) + +df_post_imputation_ami_clean <- ( + df_post_imputation_ami %>% + select( + # outcome_cox_dates_ami, # excluded + patient_id, + index_date, + end_date_exposure, + end_date_outcome, + sub_bin_covidhistory, + sub_cat_covidhospital, + exp_date_covid, + out_date_ami, + out_date_stroke_sahhs, + cov_num_age, + cov_cat_sex, + cov_num_bmi, + cov_cat_ethnicity, + cov_cat_imd, + cov_cat_smoking, + cov_bin_carehome, + cov_bin_hcworker, + cov_bin_dementia, + cov_bin_liver_disease, + cov_bin_ckd, + cov_bin_cancer, + cov_bin_hypertension, + cov_bin_diabetes, + cov_bin_obesity, + cov_bin_copd, + cov_bin_ami, + cov_bin_depression, + cov_bin_stroke_all, + cov_bin_other_ae, + cov_bin_vte, + cov_bin_hf, + cov_bin_angina, + cov_bin_lipidmed, + cov_bin_antiplatelet, + cov_bin_anticoagulant, + cov_bin_cocp, + cov_bin_hrt, + cens_date_dereg, + cens_date_death, + strat_cat_region, + vax_cat_jcvi_group, + cov_bin_covid, + cov_bin_sahhs, + # cens_status_ami, # excluded + # H0, # excluded + vax_date_eligible, + vax_date_covid_1, + vax_date_covid_2, + vax_date_covid_3, + vax_date_Pfizer_1, + vax_date_Pfizer_2, + vax_date_Pfizer_3, + vax_date_AstraZeneca_1, + vax_date_AstraZeneca_2, + vax_date_AstraZeneca_3, + vax_date_Moderna_1, + vax_date_Moderna_2, + vax_date_Moderna_3 + ) +) + +df_post_imputation_sahhs_clean <- ( + df_post_imputation_sahhs %>% + select( + # outcome_cox_dates_sahhs, # excluded + patient_id, + index_date, + end_date_exposure, + end_date_outcome, + sub_bin_covidhistory, + sub_cat_covidhospital, + exp_date_covid, + out_date_ami, + out_date_stroke_sahhs, + cov_num_age, + cov_cat_sex, + cov_num_bmi, + cov_cat_ethnicity, + cov_cat_imd, + cov_cat_smoking, + cov_bin_carehome, + cov_bin_hcworker, + cov_bin_dementia, + cov_bin_liver_disease, + cov_bin_ckd, + cov_bin_cancer, + cov_bin_hypertension, + cov_bin_diabetes, + cov_bin_obesity, + cov_bin_copd, + cov_bin_ami, + cov_bin_depression, + cov_bin_stroke_all, + cov_bin_other_ae, + cov_bin_vte, + cov_bin_hf, + cov_bin_angina, + cov_bin_lipidmed, + cov_bin_antiplatelet, + cov_bin_anticoagulant, + cov_bin_cocp, + cov_bin_hrt, + cens_date_dereg, + cens_date_death, + strat_cat_region, + vax_cat_jcvi_group, + cov_bin_covid, + cov_bin_sahhs, + # cens_status_ami, # excluded + # H0, # excluded + vax_date_eligible, + vax_date_covid_1, + vax_date_covid_2, + vax_date_covid_3, + vax_date_Pfizer_1, + vax_date_Pfizer_2, + vax_date_Pfizer_3, + vax_date_AstraZeneca_1, + vax_date_AstraZeneca_2, + vax_date_AstraZeneca_3, + vax_date_Moderna_1, + vax_date_Moderna_2, + vax_date_Moderna_3 + ) +) + + +# Save data after 'across' multiple imputation --- +print("Save data after 'across' multiple imputation ") + +saveRDS( + df_post_imputation_ami_clean, + file = paste0("output/dataset_clean/input_", cohort, "_clean_ami.rds"), + compress = TRUE +) + +saveRDS( + df_ami_nelsonaalen, + file = paste0("output/dataset_clean/nelson_aalen_", cohort, "_ami.rds"), + compress = TRUE +) + +saveRDS( + df_post_imputation_sahhs_clean, + file = paste0("output/dataset_clean/input_", cohort, "_clean_sahhs.rds"), + compress = TRUE +) + +saveRDS( + df_sahhs_nelsonaalen, + file = paste0("output/dataset_clean/nelson_aalen_", cohort, "_sahhs.rds"), + compress = TRUE +) diff --git a/analysis/cox_ipw/cox-ipw.R b/analysis/cox_ipw/cox-ipw.R new file mode 100644 index 0000000..f676c76 --- /dev/null +++ b/analysis/cox_ipw/cox-ipw.R @@ -0,0 +1,720 @@ +# # # # # # # # # # # # # # # # # # # # # +# This script: +# - defines its arguments +# - samples data and applies inverse probability weights +# - performs survival data setup +# - checks covariate variation +# - fits Cox model +# # # # # # # # # # # # # # # # # # # # # + +# Define flag style arguments using the optparse package ---- +library(optparse) +option_list <- list( + make_option( + "--df_input", + type = "character", + default = "input.csv", + help = "Input dataset. csv, csv.gz, rds, arrow, or a feather filename (this is assumed to be within the output directory but can be within a subdirectory) [default %default]", + metavar = "filename.csv" + ), + make_option( + "--ipw", + type = "logical", + default = TRUE, + help = "Logical, indicating whether sampling and IPW are to be applied [default %default]", + metavar = "TRUE/FALSE" + ), + make_option( + "--sample_exposed", + type = "logical", + default = FALSE, + help = "Logical, indicating whether exposed individuals should be sampled [default %default]", + metavar = "TRUE/FALSE" + ), + make_option( + "--exposure", + type = "character", + default = "exp_date_covid19_confirmed", + help = "Exposure variable name [default %default]", + metavar = "exposure_varname" + ), + make_option( + "--outcome", + type = "character", + default = "out_date_vte", + help = "Outcome variable name [default %default]", + metavar = "outcome_varname" + ), + make_option( + "--strata", + type = "character", + default = "cov_cat_region", + help = "Semi-colon separated list of variable names to be included as strata in the regression model [default %default]", + metavar = "varname_1;varname_2;..." + ), + make_option( + "--covariate_sex", + type = "character", + default = "cov_cat_sex", + help = "Variable name for the sex covariate; specify argument as NULL to model without sex covariate [default %default]", + metavar = "sex_varname" + ), + make_option( + "--covariate_age", + type = "character", + default = "cov_num_age", + help = "Variable name for the age covariate; specify argument as NULL to model without age covariate [default %default]", + metavar = "age_varname" + ), + make_option( + "--covariate_other", + type = "character", + default = "cov_cat_ethnicity;cov_num_consultation_rate;cov_bin_healthcare_worker;cov_bin_carehome_status", + help = "Semi-colon separated list of other covariates or filename of text file containing semi-colon separated list of other covariates to be included in the regression model; specify argument as NULL to run age, age squared, sex adjusted model only. Note that if you are including only a single covariate please include the semi-colon after it [default %default]", + metavar = "varname_1;varname_2;... or covariate-other.txt" + ), + make_option( + "--cox_start", + type = "character", + default = "pat_index_date", + help = "Semi-colon separated list of variable names used to define start of patient follow-up or single variable if already defined [default %default]", + metavar = "varname_1;varname_2;..." + ), + make_option( + "--cox_stop", + type = "character", + default = "death_date;out_date_vte;vax_date_covid_1", + help = "Semi-colon separated list of variable names used to define end of patient follow-up or single variable if already defined [default %default]", + metavar = "varname_1;varname_2;..." + ), + make_option( + "--study_start", + type = "character", + default = "2021-06-01", + help = "Study start date; this is used to remove events outside study dates [default %default]", + metavar = "YYYY-MM-DD" + ), + make_option( + "--study_stop", + type = "character", + default = "2021-12-14", + help = "Study end date; this is used to remove events outside study dates [default %default]", + metavar = "YYYY-MM-DD" + ), + make_option( + "--cut_points", + type = "character", + default = "28;197", + help = "Semi-colon separated list of cut points to be used to define time post exposure [default %default]", + metavar = "cutpoint_1;cutpoint_2" + ), + make_option( + "--controls_per_case", + type = "integer", + default = 20L, + help = "Number of controls to retain per case in the analysis [default %default]", + metavar = "integer" + ), + make_option( + "--total_event_threshold", + type = "integer", + default = 50L, + help = "Number of events that must be present for any model to run [default %default]", + metavar = "integer" + ), + make_option( + "--episode_event_threshold", + type = "integer", + default = 5L, + help = "Number of events that must be present in a time period; if threshold is not met, time periods are collapsed [default %default]", + metavar = "integer" + ), + make_option( + "--covariate_threshold", + type = "integer", + default = 5L, + help = "Minimum number of individuals per covariate level for covariate to be retained [default %default]", + metavar = "integer" + ), + make_option( + "--age_spline", + type = "logical", + default = TRUE, + help = "Logical, if age should be included in the model as a spline with knots at 0.1, 0.5, 0.9 [default %default]", + metavar = "TRUE/FALSE" + ), + make_option( + "--df_output", + type = "character", + default = "results.csv", + help = "Output data csv filename (this is assumed to be within the output directory but can be within a subdirectory) [default %default]", + metavar = "filename.csv" + ), + make_option( + "--seed", + type = "integer", + default = 137L, + help = "Random number generator seed passed to IPW sampling [default %default]", + metavar = "integer" + ), + make_option( + "--save_analysis_ready", + type = "character", + default = "", + help = "If provided, analysis ready Stata dataset filename (this is assumed to be within the output directory but can be within a subdirectory) [default %default]", + metavar = "filename.dta" + ), + make_option( + "--run_analysis", + type = "logical", + default = TRUE, + help = "Logical, if analysis should be run [default %default]", + metavar = "TRUE/FALSE" + ), + make_option( + "--config", + type = "character", + default = "", + help = "Config parsed from the YAML", + metavar = "" + ) +) +opt_parser <- OptionParser( + usage = "cox-ipw:[version] [options]", + option_list = option_list +) +opt <- parse_args(opt_parser) + +# Parse the config from YAML +if (opt$config != "") { + config <- jsonlite::fromJSON(opt$config) + opt <- modifyList(opt, config) +} + +# Record input arguments -------------------------------------------------------- +print("Record input arguments") + +record_args <- data.frame( + argument = names(opt), + value = unlist(opt), + stringsAsFactors = FALSE +) + +row.names(record_args) <- NULL + +print(record_args) + +write.csv( + record_args, + file = paste0("output/", gsub(".csv", "-args.csv", opt$df_output)), + row.names = FALSE +) + +# Import libraries ------------------------------------------------------------- +print("Import libraries") + +library(survival) +library(magrittr) + +# Import functions ------------------------------------------------------------- +print("Import functions") + +lapply( + list.files("analysis/cox_ipw/", full.names = TRUE, pattern = "fn-"), + source +) +source("analysis/utility.R") + +# Read in covariate_other if a filename ---- +# Assuming it's a filename if the character string contains no semi-colons +# Note to future self/developer - could change the second condition here to look for say .txt in string +if (!is.null(opt$covariate_other) && !grepl(";", opt$covariate_other)) { + # Check if file exists + if (!file.exists(opt$covariate_other)) { + stop(paste("The file", opt$covariate_other, "does not exist, please check how you have specified the covariate_other option.")) + } + + # Read in text file contents + # Get the size of the file in bytes + covariate_other_file_size <- file.info(opt$covariate_other)$size + + # Read the entire file content as a single character string + covariate_other_single_string <- readChar(opt$covariate_other, covariate_other_file_size) + # Strip newlines from (likely the end of) string (but could be multi-line text file) + covariate_other_single_string <- gsub("[\r\n]", "", covariate_other_single_string) + + # Overwrite opt$covariate_other so it can be processed as the direct specification of this option + opt$covariate_other <- covariate_other_single_string +} + +# If last character of opt$covariate_other is a ; replace with nothing +opt$covariate_other <- trimws(opt$covariate_other, which = "right", whitespace = ";") + +# Separate list arguments ------------------------------------------------------ +print("Separate list arguments") + +optlistargs <- c( + "strata", + "covariate_other", + "cut_points", + "cox_start", + "cox_stop" +) +for (i in seq_len(length(optlistargs))) { + tmp <- opt[optlistargs[i]] + if (tmp[1] == "NULL") { + assign(optlistargs[i], NULL) + } else { + tmp <- stringr::str_split(as.vector(tmp), ";")[[1]] + assign(optlistargs[i], tmp) + } +} +rm(tmp) + +# Make numeric arguments numeric ----------------------------------------------- +print("Make numeric arguments numeric") + +cut_points <- as.numeric(cut_points) +controls_per_case <- opt$controls_per_case +total_event_threshold <- opt$total_event_threshold +episode_event_threshold <- opt$episode_event_threshold +covariate_threshold <- opt$covariate_threshold + +# Load data -------------------------------------------------------------------- +print("Load data") + +if (grepl(".csv.gz", opt$df_input)) { + R.utils::gunzip(paste0("output/", opt$df_input), remove = FALSE) + opt$df_input <- substr(opt$df_input, 1, nchar(opt$df_input) - 3) +} +if (grepl(".csv", opt$df_input)) { + data <- readr::read_csv(paste0("output/", opt$df_input)) +} else if (grepl(".rds", opt$df_input)) { + data <- readr::read_rds(paste0("output/", opt$df_input)) +} else if (grepl(".feather", opt$df_input) || grepl(".arrow", opt$df_input)) { + data <- arrow::read_feather(paste0("output/", opt$df_input)) +} + +print(summary(data)) + +# Make binary variables logical ------------------------------------------------ +print("Make binary variables logical") + +var_bin <- colnames(data)[grepl("_bin_", colnames(data))] +data[, var_bin] <- lapply(data[, var_bin], as.logical) + +# Make date variables dates ---------------------------------------------------- +print("Make date variables dates") + +var_date <- colnames(data)[grepl("_date", colnames(data))] +data[, var_date] <- lapply( + data[, var_date], + function(x) as.Date(x, origin = "1970-01-01") +) + +# Make categorical variables factors ------------------------------------------- +print("Make categorical variables factors") + +var_cat <- colnames(data)[grepl("_cat_", colnames(data))] +data[, var_cat] <- lapply(data[, var_cat], as.factor) + +# Make numerical variables numerical ------------------------------------------- +print(" Make numerical variables numerical") + +var_num <- colnames(data)[grepl("_num_", colnames(data))] +data[, var_num] <- lapply(data[, var_num], as.numeric) + +# Restrict to core variables --------------------------------------------------- +print("Restrict to core variables") + +core <- c("patient_id", opt$exposure, opt$outcome, cox_start, cox_stop) +input <- data[, core] +print(paste0("Core variables: ", paste0(core, collapse = ", "))) + +# Give generic names to variables ---------------------------------------------- +print("Give generic names to variables") + +input <- dplyr::rename( + input, + "outcome" = tidyselect::all_of(opt$outcome), + "exposure" = tidyselect::all_of(opt$exposure) +) + +cox_start <- gsub(opt$outcome, "outcome", cox_start) +cox_start <- gsub(opt$exposure, "exposure", cox_start) + +cox_stop <- gsub(opt$outcome, "outcome", cox_stop) +cox_stop <- gsub(opt$exposure, "exposure", cox_stop) + +print(summary(input)) + +# Specify study dates ---------------------------------------------------------- +print("Specify study dates") + +input$study_start <- as.Date(opt$study_start) +input$study_stop <- as.Date(opt$study_stop) + +print(summary(input)) + +# Specify follow-up dates ------------------------------------------------------ +print("Specify follow-up dates") + +input$fup_start <- do.call( + pmax, + c(input[, c("study_start", cox_start)], list(na.rm = TRUE)) +) + +input$fup_stop <- do.call( + pmin, + c(input[, c("study_stop", cox_stop, "outcome")], list(na.rm = TRUE)) +) + +input <- input[input$fup_stop >= input$fup_start, ] + +print(summary(input)) + +# Remove exposures and outcomes outside follow-up ------------------------------ +print("Remove exposures and outcomes outside follow-up") + +print(paste0( + "Exposure data range: ", + min(input$exposure, na.rm = TRUE), + " to ", + max(input$exposure, na.rm = TRUE) +)) +print(paste0( + "Outcome data range: ", + min(input$outcome, na.rm = TRUE), + " to ", + max(input$outcome, na.rm = TRUE) +)) + +input <- input %>% + dplyr::mutate( + exposure = replace( + exposure, + which(exposure > fup_stop | exposure < fup_start), + NA + ), + outcome = replace( + outcome, + which(outcome > fup_stop | outcome < fup_start), + NA + ) + ) + +print(paste0( + "Exposure data range: ", + min(input$exposure, na.rm = TRUE), + " to ", + max(input$exposure, na.rm = TRUE) +)) +print(paste0( + "Outcome data range: ", + min(input$outcome, na.rm = TRUE), + " to ", + max(input$outcome, na.rm = TRUE) +)) + +# Make indicator variable for outcome status ----------------------------------- +print("Make indicator variable for outcome status") + +input$outcome_status <- input$outcome == input$fup_stop & + !is.na(input$outcome) & + !is.na(input$fup_stop) + +print(table(input$outcome_status)) + +# Sample control population ---------------------------------------------------- + +N_total <- nrow(input) +print(paste0("N_total = ", N_total)) + +N_exposed <- nrow(input[!is.na(input$exposure), ]) +print(paste0("N_exposed = ", N_exposed)) + +if (opt$ipw == TRUE) { + print("Sample control population") + input <- ipw_sample( + df = input, + controls_per_case = controls_per_case, + seed = opt$seed, + sample_exposed = opt$sample_exposed + ) + print(paste0("After sampling, N_total = ", nrow(input))) +} + +print(summary(input)) + +# Define episode labels -------------------------------------------------------- +print("Define episode labels") + +if (length(cut_points) == 1L) { + time_period_labels <- paste0("days0_", cut_points) +} else { + time_period_labels <- paste0( + "days", + c("0", cut_points[1:(length(cut_points) - 1)]), + "_", + cut_points + ) +} + +episode_labels <- data.frame( + episode = 0:length(cut_points), + time_period = c("days_pre", time_period_labels), + stringsAsFactors = FALSE +) + +print(episode_labels) + +# Survival data setup ---------------------------------------------------------- +print("Survival data setup") + +data_surv <- survival_data_setup( + df = input, + cut_points = cut_points, + episode_labels = episode_labels +) + +print(summary(data_surv)) + +# Calculate events in each time period ----------------------------------------- +print("Calculate events in each time period") + +if (nrow(data_surv[data_surv$outcome_status == 1, ]) > 0) { + episode_info <- get_episode_info( + df = data_surv, + cut_points = cut_points, + episode_labels = episode_labels, + ipw = opt$ipw + ) + + print(episode_info) +} + +# STOP if the total number of events is insufficient --------------------------- + +if ( + sum(episode_info[episode_info$time_period != "days_pre", ]$N_events) < + total_event_threshold +) { + results <- data.frame( + error = paste0( + "The total number of post-exposure events is less than the prespecified limit (limit = ", + total_event_threshold, + ")." + ) + ) + print(results$error) +} else { + # Add strata information to data --------------------------------------------- + print("Add strata information to data") + + data_strata <- data[, c("patient_id", strata)] + data_surv <- merge(data_surv, data_strata, by = "patient_id", all.x = TRUE) + + print(summary(data_surv)) + + # Add age covariates --------------------------------------------------------- + print("Add age covariates") + + if (opt$covariate_age != "NULL") { + data_covar <- data[, c("patient_id", opt$covariate_age)] + + data_covar <- dplyr::rename( + data_covar, + "cov_num_age" = tidyselect::all_of(opt$covariate_age) + ) + + data_covar$cov_num_age_sq <- data_covar$cov_num_age^2 + + data_surv <- merge(data_surv, data_covar, by = "patient_id", all.x = TRUE) + + print(summary(data_surv)) + } + + # Add sex covariate ---------------------------------------------------------- + print("Add sex covariate") + + if (opt$covariate_sex != "NULL") { + data_covar <- data[, c("patient_id", opt$covariate_sex)] + + data_covar <- dplyr::rename( + data_covar, + "cov_cat_sex" = tidyselect::all_of(opt$covariate_sex) + ) + + data_surv <- merge(data_surv, data_covar, by = "patient_id", all.x = TRUE) + + print(summary(data_surv)) + } + + # If additional covariates are specified, add covariate data ----------------- + + covariate_removed <- NULL + covariate_collapsed <- NULL + + if (!is.null(covariate_other)) { + # Add covariate information to data ---------------------------------------- + print("Additional covariates specified: Add covariate information to data") + + data_covar <- data[, c("patient_id", covariate_other)] + + data_surv <- merge(data_surv, data_covar, by = "patient_id", all.x = TRUE) + + # Remove covariates with insufficient variation ---------------------------- + print( + "Additional covariates specified: Remove covariates with insufficient variation" + ) + + tmp <- check_covariates( + df = data_surv, + covariate_threshold = covariate_threshold, + strata = strata + ) + + data_surv <- tmp$df + covariate_removed <- tmp$covariate_removed + covariate_collapsed <- tmp$covariate_collapsed + strata_warning <- tmp$strata_warning + rm(tmp) + + print(summary(data_surv)) + print(paste0("Removed covariates: ", paste0(covariate_removed))) + print(paste0("Collapsed covariates: ", paste0(covariate_collapsed))) + } + + # Save analysis ready dataset ------------------------------------------------ + print("Save analysis ready dataset") + + data_surv[, c("study_start", "study_stop")] <- NULL + + if (opt$save_analysis_ready != "") { + foreign::write.dta(data_surv, paste0("output/", opt$save_analysis_ready)) + } else { + analysis_ready_empty <- data.frame() + foreign::write.dta(analysis_ready_empty, "output/analysis_ready_empty.dta") + } + + if (opt$run_analysis == TRUE) { + # Perform Cox modelling ---------------------------------------------------- + print("Perform Cox modelling") + + results <- fit_model( + df = data_surv, + time_periods = episode_info[ + episode_info$time_period != "days_pre", + ]$time_period, + covariates = covariate_other, + strata = strata, + age_spline = opt$age_spline, + covariate_removed = covariate_removed, + covariate_collapsed = covariate_collapsed, + ipw = opt$ipw + ) + + # Merge results with number of events and person time ---------------------- + print("Merge results with number of events and person time") + + results <- merge( + results, + episode_info[, c( + "time_period", + "N_events", + "person_time_total", + "outcome_time_median" + )], + by.x = "term", + by.y = "time_period", + all.x = TRUE + ) + + print(summary(results)) + + # Add dummy row for days_pre term -------------------------------------------- + print("Add dummy row for days_pre term") + + tmp <- data.frame( + term = "days_pre", + lnhr = NA, + se_lnhr = NA, + model = c("mdl_age_sex", "mdl_max_adj"), + surv_formula = c( + results[results$model == "mdl_age_sex", ]$surv_formula[1], + results[results$model == "mdl_max_adj", ]$surv_formula[1] + ), + covariate_removed = "", + covariate_collapsed = "", + obs_warning = "", + N_events = episode_info[ + episode_info$time_period == "days_pre", + ]$N_events, + person_time_total = episode_info[ + episode_info$time_period == "days_pre", + ]$person_time_total, + outcome_time_median = episode_info[ + episode_info$time_period == "days_pre", + ]$outcome_time_median, + stringsAsFactors = FALSE + ) + + results <- rbind(results, tmp) + + print(summary(results)) + + # Tidy variables for outputting -------------------------------------------- + print("Tidy variables for outputting") + + results$N_total <- N_total + results$N_exposed <- N_exposed + + results$exposure <- opt$exposure + results$outcome <- opt$outcome + + results$input <- opt$df_input + + results$hr <- exp(results$lnhr) + results$conf_low <- exp(results$lnhr - qnorm(0.975) * results$se_lnhr) + results$conf_high <- exp(results$lnhr + qnorm(0.975) * results$se_lnhr) + + results$strata_warning <- strata_warning + + results$cox_ipw <- "v0.0.41" + + results <- results[ + order(results$model), + c( + "model", + "exposure", + "outcome", + "term", + "lnhr", + "se_lnhr", + "hr", + "conf_low", + "conf_high", + "N_total", + "N_exposed", + "N_events", + "person_time_total", + "outcome_time_median", + "covariate_collapsed", + "strata_warning", + "obs_warning", + "surv_formula", + "input", + "cox_ipw" + ) + ] + } else { + results <- NULL + } +} + +# Save output ------------------------------------------------------------------ +print("Save output") + +write.csv(results, file = paste0("output/", opt$df_output), row.names = FALSE) + +print(summary(results)) diff --git a/analysis/cox_ipw/fn-check_covariates.R b/analysis/cox_ipw/fn-check_covariates.R new file mode 100644 index 0000000..013dbc8 --- /dev/null +++ b/analysis/cox_ipw/fn-check_covariates.R @@ -0,0 +1,210 @@ +check_covariates <- function(df, covariate_threshold, strata) { + library(magrittr) + + # Identify non-numeric covariates to remove ---------------------------------- + print("Identify non-numeric covariates to remove") + + covariate_removed <- NULL + + for (i in unique(c(colnames(df)[grepl("cov_", colnames(df))], strata))) { + # Consider non-numeric covariates ------------------------------------------ + + if (grepl("cov_bin_", i) || grepl("cov_cat_", i)) { + print(paste0("Covariate: ", i)) + + # Calculate frequency for each level ------------------------------------- + print("Calculate frequency for each level") + + tmp <- unique(df[ + df$exposure_status == 1 & df$outcome_status == 1, + c("patient_id", i) + ]) + freq <- data.frame(table(tmp[, i])) + + print(freq) + + # Add covariates to removal list if they fall below covariate threshold -- + + if ( + nrow(freq[ + freq$Freq <= covariate_threshold & freq$Var1 != "Missing", + ]) > + 0 + ) { + print("Add covariate to removal list") + covariate_removed <- c(covariate_removed, i) + } + + # Add covariates to removal list if the covariate has a single level ----- + + if (nrow(freq) == 1) { + print("Add covariate to removal list") + covariate_removed <- c(covariate_removed, i) + } + } + } + + # Collapse special case covariates ------------------------------------------- + + covariate_collapsed <- NULL + + ## Region + + if ("cov_cat_region" %in% covariate_removed) { + print("Collapsing region as special case") + + df <- df %>% + dplyr::mutate( + cov_cat_region = dplyr::case_when( + cov_cat_region == "North East" ~ "Northern England", + cov_cat_region == "North West" ~ "Northern England", + cov_cat_region == "Yorkshire and The Humber" ~ "Northern England", + cov_cat_region == "London" ~ "Southern England", + cov_cat_region == "South East" ~ "Southern England", + cov_cat_region == "East Midlands" ~ "Midlands", + cov_cat_region == "West Midlands" ~ "Midlands", + cov_cat_region == "South West" ~ "Southern England", + cov_cat_region == "East" ~ "Southern England" + ) + ) + + df$cov_cat_region <- factor(df$cov_cat_region) + df$cov_cat_region <- relevel(df$cov_cat_region, ref = "Southern England") + + covariate_removed <- setdiff(covariate_removed, "cov_cat_region") + covariate_collapsed <- c(covariate_collapsed, "cov_cat_region") + } + + ## Deprivation + + if ("cov_cat_deprivation" %in% covariate_removed) { + print("Collapsing deprivation as special case") + + df <- df %>% + dplyr::mutate( + cov_cat_deprivation = dplyr::case_when( + cov_cat_deprivation == "1-2 (most deprived)" ~ "1-4", + cov_cat_deprivation == "3-4" ~ "1-4", + cov_cat_deprivation == "5-6" ~ "5-6", + cov_cat_deprivation == "7-8" ~ "7-10", + cov_cat_deprivation == "9-10 (least deprived)" ~ "7-10" + ) + ) + + df$cov_cat_deprivation <- ordered( + df$cov_cat_deprivation, + levels = c("1-4", "5-6", "7-10") + ) + + covariate_removed <- setdiff(covariate_removed, "cov_cat_deprivation") + covariate_collapsed <- c(covariate_collapsed, "cov_cat_deprivation") + } + + # Smoking status ------------------------------------------------------------- + + if ("cov_cat_smoking_status" %in% covariate_removed) { + print("Collapsing smoking status as special case") + + df <- df %>% + dplyr::mutate( + cov_cat_smoking_status = dplyr::case_when( + cov_cat_smoking_status == "Never smoker" ~ "Never smoker", + cov_cat_smoking_status == "Ever smoker" ~ "Ever smoker", + cov_cat_smoking_status == "Current smoker" ~ "Ever smoker", + cov_cat_smoking_status == "Missing" ~ "Missing" + ) + ) + + df$cov_cat_smoking_status <- ordered( + df$cov_cat_smoking_status, + levels = c("Never smoker", "Ever smoker", "Missing") + ) + + covariate_removed <- setdiff(covariate_removed, "cov_cat_smoking_status") + } + + # Check special case collapsed covariates ------------------------------------ + + for (i in c( + "cov_cat_deprivation", + "cov_cat_smoking_status", + "cov_cat_region" + )) { + if (i %in% colnames(df)) { + print(paste0("Rechecking covariate: ", i)) + + # Calculate frequency for each level among exposed with outcome ---------- + print("Calculate frequency for each level among exposed with outcome") + + tmp <- unique(df[ + df$exposure_status == 1 & df$outcome_status == 1, + c("patient_id", i) + ]) + freq <- data.frame(table(tmp[, i])) + + print(freq) + + # Add covariates to removal list if they fall below covariate threshold -- + + if ( + nrow(freq[ + freq$Freq <= covariate_threshold & freq$Var1 != "Missing", + ]) > + 0 + ) { + print("Add covariate to removal list") + covariate_removed <- c(covariate_removed, i) + covariate_collapsed <- setdiff(covariate_collapsed, i) + } + + # Add covariates to removal list if the covariate has a single level ----- + + if (nrow(freq) == 1) { + print("Add covariate to removal list") + covariate_removed <- c(covariate_removed, i) + covariate_collapsed <- setdiff(covariate_collapsed, i) + } + } + } + + # Check strata variables meet covariate threshold ---------------------------- + print("Check strata variables meet covariate threshold") + + strata_warning <- "" + + if (length(intersect(covariate_removed, strata)) > 0) { + strata_warning <- paste0( + intersect(covariate_removed, strata), + collapse = ";" + ) + for (i in intersect(covariate_removed, strata)) { + tmp <- unique(df[ + df$exposure_status == 1 & df$outcome_status == 1, + c("patient_id", i) + ]) + freq <- data.frame(table(tmp[, i])) + print(paste0( + "Warning: strata variable '", + i, + "' does not meet covariate threshold" + )) + print(freq) + } + } + + # Remove covariates ---------------------------------------------------------- + print("Remove covariates") + + df <- df[, !(colnames(df) %in% setdiff(covariate_removed, strata))] + + # Return data and list of removed covariates --------------------------------- + + output <- list( + df = df, + covariate_removed = setdiff(covariate_removed, strata), + covariate_collapsed = covariate_collapsed, + strata_warning = strata_warning + ) + + return(output) +} diff --git a/analysis/cox_ipw/fn-fit_model.R b/analysis/cox_ipw/fn-fit_model.R new file mode 100644 index 0000000..622a4f6 --- /dev/null +++ b/analysis/cox_ipw/fn-fit_model.R @@ -0,0 +1,231 @@ +fit_model <- function( + df, + time_periods, + covariates, + strata, + age_spline, + covariate_removed, + covariate_collapsed, + ipw +) { + + # Calculate stacked weights ----- + print("Calculate stacked weights") + df$cox_stacked_weight <- generate_weights(initial_weights = df$cox_weight) + + # Check weights in logs ----- + print("Initial weights:") + print(summary(df$cox_weight)) + print("Stacked weights for across MI method:") + print(summary(df$cox_stacked_weight)) + + # Define model formula ------------------------------------------------------- + print("Define model formula") + + surv_formula <- paste0( + "Surv(tstart, tstop, outcome_status) ~ ", + paste(time_periods, collapse = " + "), + ifelse("cov_cat_sex" %in% colnames(df), " + cov_cat_sex", ""), + ifelse( + is.null(strata), + "", + paste(" +", paste0("rms::strat(", strata, ")"), collapse = " + ") + ), + ifelse(isTRUE(ipw), " + cluster(patient_id)", "") + ) + + # Add age covariate, specifying knot placement for age spline if applicable -- + + if ("cov_num_age" %in% colnames(df)) { + print("Add age covariate") + + if (age_spline == TRUE) { + print("Specify knot placement for age spline") + + knot_placement <- as.numeric(quantile( + df$cov_num_age, + probs = c(0.1, 0.5, 0.9) + )) + + print(paste0( + "Knots will be placed at: ", + paste0(knot_placement, collapse = ", ") + )) + + surv_formula <- paste0( + surv_formula, + " + rms::rcs(cov_num_age, parms=knot_placement)" + ) + } else { + surv_formula <- paste0(surv_formula, " + cov_num_age + cov_num_age_sq") + } + } + + print(surv_formula) + + # Fit Cox model ---------------------------------------------------------------- + print("Fit Cox model") + + dd <<- rms::datadist(df) + + withr::local_options(list( + datadist = "dd", + contrasts = c("contr.treatment", "contr.treatment") + )) + + if (ipw == TRUE) { + N_obs_in <- nrow(df) + + fit_cox_model <- rms::cph( + formula = as.formula(surv_formula), + data = df, + weight = df$cox_stacked_weight, # weights for stacked dataset (across MI) + method = "breslow", + robust = TRUE, # use robust var -> robust SE + surv = TRUE, + x = TRUE, + y = TRUE + ) + + N_obs_out <- sum(fit_cox_model$n) + } else { + N_obs_in <- nrow(df) + + fit_cox_model <- rms::cph( + formula = as.formula(surv_formula), + data = df, + method = "breslow", + robust = TRUE, # use robust var -> robust SE + surv = TRUE, + x = TRUE, + y = TRUE + ) + + N_obs_out <- sum(fit_cox_model$n) + } + + print(fit_cox_model) + + # Format results --------------------------------------------------------------- + print("Format results") + + results <- data.frame( + term = names(fit_cox_model$coefficients), + lnhr = fit_cox_model$coefficients, + se_lnhr = sqrt(diag(vcov(fit_cox_model))), + model = "mdl_age_sex", + surv_formula = surv_formula, + covariate_removed = "", + covariate_collapsed = "", + obs_warning = ifelse( + N_obs_in == N_obs_out, + "", + paste0( + N_obs_in, + " observations provided. ", + N_obs_out, + " observations used." + ) + ), + stringsAsFactors = FALSE + ) + + row.names(results) <- NULL + + # If covariates are specified, run an additional model including them -------- + + covariates <- setdiff(covariates, covariate_removed) + + if (!is.null(covariates) && length(covariates) > 0) { + # Add covariates to model formula ------------------------------------------ + print("Add covariates to model formula") + + surv_formula_adj <- paste0( + surv_formula, + " + ", + paste(covariates, collapse = " + ") + ) + + print(surv_formula_adj) + + # Fit Cox model ------------------------------------------------------------ + print("Fit Cox model with covariates") + + dd_adj <<- rms::datadist(df) + + withr::local_options(list( + datadist = "dd_adj", + contrasts = c("contr.treatment", "contr.treatment") + )) + + if (ipw == TRUE) { + N_obs_in <- nrow(df) + + fit_cox_model_adj <- rms::cph( + formula = as.formula(surv_formula_adj), + data = df, + weight = df$cox_weight, + method = "breslow", + robust = TRUE, # use robust var -> robust SE + surv = TRUE, + x = TRUE, + y = TRUE + ) + + N_obs_out <- sum(fit_cox_model$n) + } else { + N_obs_in <- nrow(df) + + fit_cox_model_adj <- rms::cph( + formula = as.formula(surv_formula_adj), + data = df, + method = "breslow", + robust = TRUE, # use robust var -> robust SE + surv = TRUE, + x = TRUE, + y = TRUE + ) + + N_obs_out <- sum(fit_cox_model$n) + } + + print(fit_cox_model_adj) + + # Format results --------------------------------------------------------------- + print("Format results") + + results_adj <- data.frame( + term = names(fit_cox_model_adj$coefficients), + lnhr = fit_cox_model_adj$coefficients, + se_lnhr = sqrt(diag(vcov(fit_cox_model_adj))), + model = "mdl_max_adj", + surv_formula = surv_formula_adj, + covariate_removed = paste0(covariate_removed, collapse = ";"), + covariate_collapsed = paste0(covariate_collapsed, collapse = ";"), + obs_warning = ifelse( + N_obs_in == N_obs_out, + "", + paste0( + N_obs_in, + " observations provided. ", + N_obs_out, + " observations used." + ) + ), + stringsAsFactors = FALSE + ) + + row.names(results_adj) <- NULL + + # Bind to other results ---------------------------------------------------- + print("Bind to other results") + + results <- rbind(results, results_adj) + } + + # Return results ------------------------------------------------------------- + print("Return results") + + return(results) + print(summary(results)) +} diff --git a/analysis/cox_ipw/fn-get_episode_info.R b/analysis/cox_ipw/fn-get_episode_info.R new file mode 100644 index 0000000..83e2717 --- /dev/null +++ b/analysis/cox_ipw/fn-get_episode_info.R @@ -0,0 +1,86 @@ +get_episode_info <- function(df, cut_points, episode_labels, ipw) { + library(magrittr) + + # Calculate number of events per episode ------------------------------------- + print("Calculate number of events per episode") + + events <- df[df$outcome_status == 1, c("patient_id", "episode")] + + events <- aggregate(episode ~ patient_id, data = events, FUN = max) + + events <- data.frame(table(events$episode), stringsAsFactors = FALSE) + + events <- dplyr::rename(events, "episode" = "Var1", "N_events" = "Freq") + + print(events) + + # Add number of events to episode info table --------------------------------- + print("Add number of events to episode info table") + + episode_info <- merge(episode_labels, events, by = "episode", all.x = TRUE) + episode_info$N_events <- ifelse( + is.na(episode_info$N_events), + 0, + episode_info$N_events + ) + + print(episode_info) + + # Calculate person-time in each episode -------------------------------------- + print("Calculate person-time in each episode") + + if (ipw == TRUE) { + tmp <- df[, c("episode", "tstart", "tstop", "cox_weight")] + tmp$person_time_total <- (tmp$tstop - tmp$tstart) * tmp$cox_weight + tmp[, c("tstart", "tstop", "cox_weight")] <- NULL + } + + if (ipw == FALSE) { + tmp <- df[, c("episode", "tstart", "tstop")] + tmp$person_time_total <- (tmp$tstop - tmp$tstart) + tmp[, c("tstart", "tstop")] <- NULL + } + + tmp <- aggregate(person_time_total ~ episode, data = tmp, FUN = sum) + + episode_info <- merge(episode_info, tmp, by = "episode", all.x = TRUE) + + print(episode_info) + + # Calculate median person-time ----------------------------------------------- + print("Calculate median person-time") + + tmp <- df[df$outcome_status == 1, c("episode", "tstart", "tstop")] + + tmp$outcome_time <- tmp$tstop - tmp$tstart + + tmp <- tmp %>% + dplyr::group_by(episode) %>% + dplyr::mutate(outcome_time_median = median(outcome_time)) %>% + dplyr::ungroup(episode) + + tmp <- unique(tmp[, c("episode", "outcome_time_median")]) + + episode_info <- merge(episode_info, tmp, by = "episode", all.x = TRUE) + + print(episode_info) + + # Add time from previous episodes to median person-time ----------------------- + print("Add time from previous episodes to median person-time") + + episode_info$add <- as.numeric(gsub( + "days", + "", + gsub("_.*", "", episode_info$time_period) + )) + episode_info$add <- ifelse(is.na(episode_info$add), 0, episode_info$add) + episode_info$outcome_time_median <- episode_info$outcome_time_median + + episode_info$add + episode_info$add <- NULL + + print(episode_info) + + # Return episode_info table -------------------------------------------------- + + return(episode_info) +} diff --git a/analysis/cox_ipw/fn-ipw_sample.R b/analysis/cox_ipw/fn-ipw_sample.R new file mode 100644 index 0000000..4b88e33 --- /dev/null +++ b/analysis/cox_ipw/fn-ipw_sample.R @@ -0,0 +1,101 @@ +ipw_sample <- function(df, controls_per_case, seed = 137, sample_exposed) { + # Set seed ------------------------------------------------------------------- + print("Set seed") + + set.seed(seed) + + # Split cases and controls --------------------------------------------------- + print("Split cases and controls") + + cases <- df[df$outcome_status == TRUE, ] + controls <- df[df$outcome_status == FALSE, ] + + print(paste0("Cases: ", nrow(cases))) + print(summary(cases)) + + print(paste0("Controls: ", nrow(controls))) + print(summary(controls)) + + # Sample controls if more than enough, otherwise retain all controls --------- + + if (sample_exposed == TRUE) { + if (nrow(cases) * controls_per_case < nrow(controls)) { + print("Sample controls, including exposed control individuals") + controls <- controls[ + sample( + seq_len(nrow(controls)), + nrow(cases) * controls_per_case, + replace = FALSE + ), + ] + controls$cox_weight <- (nrow(df) - nrow(cases)) / nrow(controls) + print(paste0( + nrow(controls), + " controls sampled with Cox weight of ", + controls$cox_weight[1] + )) + } else { + print("Retain all controls") + controls$cox_weight <- 1 + } + } + + if (sample_exposed == FALSE) { + print("Separate exposed controls so they are not sampled") + controls_exposed <- controls[!is.na(controls$exposure), ] + controls_exposed$cox_weight <- 1 + print(paste0(nrow(controls_exposed), " exposed controls")) + + print("Exposed controls:") + print(summary(controls_exposed)) + + print("Sample unexposed controls") + controls_unexposed <- controls[is.na(controls$exposure), ] + + if (nrow(cases) * controls_per_case < nrow(controls_unexposed)) { + controls_unexposed <- controls_unexposed[ + sample( + seq_len(nrow(controls_unexposed)), + nrow(cases) * controls_per_case, + replace = FALSE + ), + ] + controls_unexposed$cox_weight <- (nrow(df) - + nrow(cases) - + nrow(controls_exposed)) / + nrow(controls_unexposed) + print(paste0( + nrow(controls_unexposed), + " unexposed controls sampled with Cox weight of ", + controls_unexposed$cox_weight[1] + )) + + print("Unexposed controls:") + print(summary(controls_unexposed)) + + print("Add exposed control individuals back to control dataset") + controls <- NULL + controls <- rbind(controls_unexposed, controls_exposed) + print(paste0("Controls (N=", nrow(controls), "):")) + print(summary(controls)) + } else { + print("Insufficient controls so retain all controls") + rm(controls_exposed, controls_unexposed) + controls$cox_weight <- 1 + } + } + + # Specify cox weight for cases ----------------------------------------------- + print("Specify cox weight for cases") + + cases$cox_weight <- 1 + + # Recombine cases and controls ----------------------------------------------- + print("Recombine cases and controls") + + return_df <- rbind(cases, controls) + + # Return dataset ------------------------------------------------------------- + + return(return_df) +} diff --git a/analysis/cox_ipw/fn-survival_data_setup.R b/analysis/cox_ipw/fn-survival_data_setup.R new file mode 100644 index 0000000..7d7e113 --- /dev/null +++ b/analysis/cox_ipw/fn-survival_data_setup.R @@ -0,0 +1,159 @@ +survival_data_setup <- function(df, cut_points, episode_labels) { + # Calculate where follow-up occurs in study period --------------------------- + print("Calculate where follow-up occurs in study period") + + df$days_to_start <- as.numeric(df$fup_start - df$study_start) + df$days_to_end <- as.numeric(df$fup_stop - df$study_start) + 1 + + # Set survival data for the exposed ------------------------------------------ + print("Set survival data for the exposed") + + ## Filter data to exposed + print("Filter data to exposed") + exposed <- df[!is.na(df$exposure), ] + + ## Calculate days to exposure + print("Calculate days to exposure") + exposed$days_to_exp <- as.numeric(exposed$exposure - exposed$study_start) + + ## Put into survival format + print("Put into survival format (dataset d1)") + d1 <- exposed[, + !(colnames(exposed) %in% + c("days_to_start", "days_to_exp", "days_to_end", "outcome_status")) + ] + print(summary(d1)) + + print("Put into survival format (dataset d2)") + d2 <- exposed[, c( + "patient_id", + "days_to_start", + "days_to_exp", + "days_to_end", + "outcome_status" + )] + print(summary(d2)) + + print("Put into survival format (tmerge)") + exposed <- survival::tmerge( + data1 = d1, + data2 = d2, + id = patient_id, + outcome_status = event(days_to_end, outcome_status), + tstart = days_to_start, + tstop = days_to_end, + exposure_status = tdc(days_to_exp) + ) + print(summary(exposed)) + + # Split post-exposure time for the exposed ----------------------------------- + print("Split post-exposure time for the exposed") + + ## Filter to post-exposure data + print("Filter to post-exposure data") + exposed_post <- exposed[exposed$exposure_status == 1, ] + + ## Format tstart and tstop + print("Format tstart and tstop") + exposed_post <- dplyr::rename(exposed_post, t0 = tstart, t = tstop) + exposed_post$tstart <- 0 + exposed_post$tstop <- exposed_post$t - exposed_post$t0 + + ## Implement post-exposure cut points + print("Implement post-exposure cut points") + + print(summary(exposed_post)) + + exposed_post <- survival::survSplit( + Surv(tstop, outcome_status) ~ ., + exposed_post, + cut = cut_points, + episode = "episode" + ) + + print(summary(exposed_post)) + + ## Account for pre-exposure time + print("Account for pre-exposure time") + + exposed_post$tstart <- exposed_post$tstart + exposed_post$t0 + exposed_post$tstop <- exposed_post$tstop + exposed_post$t0 + exposed_post[, c("t0", "t")] <- NULL + + print(summary(exposed_post)) + + # Combine pre- and post-exposure time for the exposed ------------------------ + print("Combine pre- and post-exposure time for the exposed") + + ## Filter to pre-exposure data + print("Filter to pre-exposure data") + exposed_pre <- exposed[exposed$exposure_status == 0, ] + + if (nrow(exposed_pre) > 0) { + ## Set pre-exposure to be episode 0 + print("Set pre-exposure to be episode 0") + exposed_pre$episode <- 0 + print(summary(exposed_pre)) + + ## Bind pre- and post-exposure time + print("Bind pre- and post-exposure time") + exposed <- plyr::rbind.fill(exposed_pre, exposed_post) + } else { + ## No pre-exposure time so exposed = post-exposure time only + print("No pre-exposure time so exposed = post-exposure time only") + exposed <- exposed_post + } + + print(summary(exposed)) + + # Set survival data for the unexposed ---------------------------------------- + print("Set survival data for the unexposed") + + ## Filter data to unexposed + print("Filter data to unexposed") + unexposed <- df[is.na(df$exposure), ] + + ## Rename variables to survival variable names + print("Rename variables to survival variable names") + + unexposed <- dplyr::rename( + unexposed, + "tstart" = "days_to_start", + "tstop" = "days_to_end" + ) + + unexposed$exposure_status <- c(0) + unexposed$episode <- c(0) + unexposed$outcome_status <- as.numeric(unexposed$outcome_status) + + print(summary(unexposed)) + + # Combine exposed and unexposed individuals ---------------------------------- + print("Combine exposed and unexposed individuals") + + exposed <- exposed[, intersect(colnames(unexposed), colnames(exposed))] + unexposed <- unexposed[, intersect(colnames(unexposed), colnames(exposed))] + df <- rbind(exposed, unexposed) + + print(summary(df)) + + # Add indicators for episode ------------------------------------------------- + print("Add indicators for episode") + + for (i in 1:max(episode_labels$episode)) { + preserve_cols <- colnames(df) + + df$tmp <- as.numeric(df$episode == i) + + colnames(df) <- c( + preserve_cols, + episode_labels[episode_labels$episode == i, ]$time_period + ) + } + + print(summary(df)) + + # Return dataset ------------------------------------------------------------- + + return(df) +} diff --git a/analysis/cox_ipw/move-covariate-other-txt.R b/analysis/cox_ipw/move-covariate-other-txt.R new file mode 100644 index 0000000..e9575f9 --- /dev/null +++ b/analysis/cox_ipw/move-covariate-other-txt.R @@ -0,0 +1 @@ +file.copy("dummy_tables/covariate-other.txt", "output/covariate-other.txt") diff --git a/analysis/cox_ipw/move-test-data-arrow-subdir.R b/analysis/cox_ipw/move-test-data-arrow-subdir.R new file mode 100644 index 0000000..5402a84 --- /dev/null +++ b/analysis/cox_ipw/move-test-data-arrow-subdir.R @@ -0,0 +1,3 @@ +R.utils::gunzip("dummy_tables/input.csv.gz") +if (!dir.exists("output/subdir")) dir.create("output/subdir") +arrow::write_feather(readr::read_csv("dummy_tables/input.csv"), "output/subdir/input-5.arrow") diff --git a/analysis/cox_ipw/move-test-data-arrow.R b/analysis/cox_ipw/move-test-data-arrow.R new file mode 100644 index 0000000..fb9082f --- /dev/null +++ b/analysis/cox_ipw/move-test-data-arrow.R @@ -0,0 +1,2 @@ +R.utils::gunzip("dummy_tables/input.csv.gz") +arrow::write_feather(readr::read_csv("dummy_tables/input.csv"), "output/input-4.arrow") diff --git a/analysis/cox_ipw/move-test-data-csv-gz.R b/analysis/cox_ipw/move-test-data-csv-gz.R new file mode 100644 index 0000000..312793a --- /dev/null +++ b/analysis/cox_ipw/move-test-data-csv-gz.R @@ -0,0 +1 @@ +file.copy("dummy_tables/input.csv.gz", "output/input-3.csv.gz") diff --git a/analysis/cox_ipw/move-test-data.R b/analysis/cox_ipw/move-test-data.R new file mode 100644 index 0000000..2940936 --- /dev/null +++ b/analysis/cox_ipw/move-test-data.R @@ -0,0 +1 @@ +R.utils::gunzip("dummy_tables/input.csv.gz", "output/input.csv") diff --git a/analysis/create_project_actions.R b/analysis/create_project_actions.R index 52bc7dc..b4a1674 100644 --- a/analysis/create_project_actions.R +++ b/analysis/create_project_actions.R @@ -184,21 +184,69 @@ clean_data <- function(cohort, describe = describe) { } -# Create function to define post-hoc variables for clean data ------------------- -post_hoc_vars <- function(cohort) { +# Create function to make missingness table --------------------------- +table_missingness <- function(cohort) { splice( - comment(glue("post_hoc_vars_cohort_{cohort}")), + comment(glue("table_missingness_cohort_{cohort}")), action( - name = glue("post_hoc_vars_cohort_{cohort}"), + name = glue("table_missingness_cohort_{cohort}"), run = glue( - "r:latest analysis/post_hoc_vars/post_hoc_vars.R" + "r:latest analysis/table_missingness/table_missingness.R" + ), + arguments = c(c(cohort)), + needs = list( + glue("generate_input_{cohort}_clean"), + glue("generate_input_{cohort}") + ), + moderately_sensitive = list( + table_missingness = glue("output/table_missingness/table_missingness_cohort_{cohort}.csv") + ) + ) + ) +} + + +# Create function to apply across multiple imputation --------------------------- +apply_across_MI <- function(cohort) { + splice( + comment(glue("apply_across_MI_cohort_{cohort}")), + action( + name = glue("apply_across_MI_cohort_{cohort}"), + run = glue( + "r:latest analysis/apply_across_MI/apply_across_MI.R" ), arguments = c(c(cohort)), needs = list( glue("generate_input_{cohort}_clean") ), highly_sensitive = list( - cohort_clean = glue("output/dataset_clean/input_{cohort}_clean.rds") + cohort_clean_ami = glue("output/dataset_clean/input_{cohort}_clean_ami.rds"), + cohort_clean_sahhs = glue("output/dataset_clean/input_{cohort}_clean_sahhs.rds"), + nelsonaalen_ami = glue("output/dataset_clean/nelson_aalen_{cohort}_ami.rds"), + nelsonaalen_sahhs = glue("output/dataset_clean/nelson_aalen_{cohort}_sahhs.rds") + ) + ) + ) +} + + +nelson_aalen_plots <- function(cohort) { + splice( + comment(glue("nelson_aalen_plots_cohort_{cohort}")), + action( + name = glue("nelson_aalen_plots_cohort_{cohort}"), + run = glue( + "r:latest analysis/nelson_aalen_plots/nelson_aalen_plots.R" + ), + arguments = c(c(cohort)), + needs = list( + glue("apply_across_MI_cohort_{cohort}") + ), + highly_sensitive = list( + nelsonaalen_plot_ami = glue("output/dataset_clean/nelson_aalen_{cohort}_ami.png"), + nelsonaalen_plot_ami_ggplot = glue("output/dataset_clean/nelson_aalen_{cohort}_ami_ggplot.png"), + nelsonaalen_plot_sahhs = glue("output/dataset_clean/nelson_aalen_{cohort}_sahhs.png"), + nelsonaalen_plot_sahhs_ggplot = glue("output/dataset_clean/nelson_aalen_{cohort}_sahhs_ggplot.png") ) ) ) @@ -216,10 +264,11 @@ generate_subsample_cohort <- function(cohort) { ), arguments = c(c(cohort)), needs = list( - glue("post_hoc_vars_cohort_{cohort}") # , glue("make_model_input-{name}") + glue("apply_across_MI_cohort_{cohort}") # , glue("make_model_input-{name}") ), highly_sensitive = list( - cohort_clean_subsample = glue("output/generate_subsample/input_{cohort}_clean_subsample.rds") + cohort_clean_subsample_ami = glue("output/generate_subsample/input_{cohort}_clean_subsample_ami.rds"), + cohort_clean_subsample_sahhs = glue("output/generate_subsample/input_{cohort}_clean_subsample_sahhs.rds") ) ) ) @@ -274,13 +323,19 @@ table1 <- function(cohort, ages = "18;40;60;80", preex = "All") { name = glue("table1-cohort_{cohort}{preex_str}"), run = "r:v2 analysis/table1/table1.R", arguments = c(c(cohort), c(ages), c(preex)), - needs = list(glue("post_hoc_vars_cohort_{cohort}")), + needs = list(glue("apply_across_MI_cohort_{cohort}")), moderately_sensitive = list( - table1 = glue( - "output/table1/table1-cohort_{cohort}{preex_str}.csv" + table1_ami = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_ami.csv" + ), + table1_ami_midpoint6 = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_ami-midpoint6.csv" ), - table1_midpoint6 = glue( - "output/table1/table1-cohort_{cohort}{preex_str}-midpoint6.csv" + table1_sahhs = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_sahhs.csv" + ), + table1_sahhs_midpoint6 = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_sahhs-midpoint6.csv" ) ) ) @@ -302,13 +357,22 @@ table1_subsample <- function(cohort, ages = "18;40;60;80", preex = "All") { name = glue("table1-cohort_{cohort}{preex_str}_subsample"), run = "r:v2 analysis/table1/table1_subsample.R", arguments = c(c(cohort), c(ages), c(preex)), - needs = list(glue("generate_subsample_cohort_{cohort}")), + needs = list( + glue("apply_across_MI_cohort_{cohort}"), + glue("generate_subsample_cohort_{cohort}") + ), moderately_sensitive = list( - table1_subsample = glue( - "output/table1/table1-cohort_{cohort}{preex_str}_subsample.csv" + table1_ami_subsample = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_ami_subsample.csv" + ), + table1_ami_midpoint6_subsample = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_ami-midpoint6_subsample.csv" + ), + table1_sahhs_subsample = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_sahhs_subsample.csv" ), - table1_midpoint6_subsample = glue( - "output/table1/table1-cohort_{cohort}{preex_str}-midpoint6_subsample.csv" + table1_sahhs_midpoint6_subsample = glue( + "output/table1/table1-cohort_{cohort}{preex_str}_sahhs-midpoint6_subsample.csv" ) ) ) @@ -501,7 +565,7 @@ apply_model_function <- function( action( name = glue("make_model_input-{name}"), run = glue("r:latest analysis/model/make_model_input.R {name}"), - needs = as.list(glue("post_hoc_vars_cohort_{cohort}")), + needs = as.list(glue("apply_across_MI_cohort_{cohort}")), highly_sensitive = list( model_input = glue("output/model/model_input-{name}.rds") ) @@ -509,7 +573,7 @@ apply_model_function <- function( action( name = glue("cox_ipw-{name}"), run = glue( - "cox-ipw:v0.0.39 + "r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-{name}.rds --ipw={ipw} --exposure=exp_date @@ -589,7 +653,7 @@ apply_lasso_cox_model_function <- function( action( name = glue("lasso_cox_ipw-{name}"), run = glue( - "cox-ipw:v0.0.39 + "r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-{name}.rds --ipw={ipw} --exposure=exp_date @@ -648,7 +712,7 @@ apply_lasso_X_cox_model_function <- function( action( name = glue("lasso_X_cox_ipw-{name}"), run = glue( - "cox-ipw:v0.0.39 + "r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-{name}.rds --ipw={ipw} --exposure=exp_date @@ -707,7 +771,7 @@ apply_lasso_union_cox_model_function <- function( action( name = glue("lasso_union_cox_ipw-{name}"), run = glue( - "cox-ipw:v0.0.39 + "r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-{name}.rds --ipw={ipw} --exposure=exp_date @@ -949,6 +1013,56 @@ make_lasso_union_model_output <- function(subgroup) { } +# Create funtion for making combined table/venn outputs ------------------------ + +make_table1_output <- function(action_name, cohort, subgroup = "") { + cohort_names <- stringr::str_split(as.vector(cohort), ";")[[1]] + if (subgroup == "All" | subgroup == "") { + sub_str <- "" + } else { + if (grepl("preex", subgroup)) { + sub_str <- paste0("-", subgroup) + } else { + sub_str <- paste0("-sub_", subgroup) + } + } + + splice( + comment(glue("Generate make-{action_name}{sub_str}-output")), + action( + name = glue("make-{action_name}{sub_str}-output"), + run = "r:v2 analysis/make_output/make_table1_output.R", + arguments = unlist(lapply( + list( + c(action_name, cohort, subgroup) + ), + function(x) { + x[x != ""] + } + )), + needs = c( + glue("apply_across_MI_cohort_prevax"), + as.list( + paste0( + action_name, + "-cohort_", + cohort_names, + sub_str + ) + ) + ), + moderately_sensitive = list( + other_output_ami_midpoint6 = glue( + "output/make_output/{action_name}{sub_str}_output_ami_midpoint6.csv" + ), + other_output_sahhs_midpoint6 = glue( + "output/make_output/{action_name}{sub_str}_output_sahhs_midpoint6.csv" + ) + ) + ) + ) +} + # Create funtion for making combined table/venn outputs ------------------------ make_other_output <- function(action_name, cohort, subgroup = "") { @@ -1047,15 +1161,33 @@ actions_list <- splice( ) ), - ## Define post hoc variables ---------------------------------- + ## Make missingness table ---------------------------------------------------- + + splice( + unlist( + lapply(cohorts, function(x) table_missingness(cohort = x)), + recursive = FALSE + ) + ), + ## Apply across multiple imputation ------------------------------------------ + splice( unlist( - lapply(cohorts, function(x) post_hoc_vars(cohort = x)), + lapply(cohorts, function(x) apply_across_MI(cohort = x)), recursive = FALSE ) ), + ## Make Nelson-Aalen plots --------------------------------------------------- + + # splice( + # unlist( + # lapply(cohorts, function(x) nelson_aalen_plots(cohort = x)), + # recursive = FALSE + # ) + # ), + ## Generate 10% subsample study population ---------------------------------- splice( @@ -1113,7 +1245,7 @@ actions_list <- splice( ), splice( - make_other_output( + make_table1_output( action_name = "table1", cohort = paste0(cohorts, collapse = ";"), subgroup = "" diff --git a/analysis/dataset_definition/codelists.py b/analysis/dataset_definition/codelists.py index 2b0caa4..2252156 100644 --- a/analysis/dataset_definition/codelists.py +++ b/analysis/dataset_definition/codelists.py @@ -114,6 +114,12 @@ column="code" ) +## BMI Numerical +bmi_numeric_nhsd = codelist_from_csv( + "codelists/reducehf-body-mass-index-numeric-value.csv", + column="code" +) + ## Severe obesity code recorded sev_obesity_primis = codelist_from_csv( "codelists/primis-covid19-vacc-uptake-sev_obesity.csv", diff --git a/analysis/dataset_definition/variable_helper_functions.py b/analysis/dataset_definition/variable_helper_functions.py index ea1aea6..a6be899 100644 --- a/analysis/dataset_definition/variable_helper_functions.py +++ b/analysis/dataset_definition/variable_helper_functions.py @@ -1,5 +1,6 @@ import operator -from ehrql import case, when +from ehrql import case, when, days +from ehrql.codes import CTV3Code from functools import reduce # for function building, e.g. any_of from ehrql.tables.tpp import ( apcs, @@ -231,3 +232,22 @@ def get_latest_ethnicity( ) return ethnicity_combined + + +### BMI calculation +def most_recent_bmi(*, patient_date_of_birth, minimum_age_at_measurement, where=True): + age_threshold = patient_date_of_birth + days( + # This is obviously inexact but, given that the dates of birth are rounded to + # the first of the month anyway, there's no point trying to be more accurate + int(365.25 * minimum_age_at_measurement) + ) + return ( + # This captures just explicitly recorded BMI observations rather than attempting + # to calculate it from height and weight measurements. Investigation has shown + # this to have no real benefit it terms of coverage or accuracy. + clinical_events.where(where) + .where(clinical_events.ctv3_code == CTV3Code("22K..")) + .where(clinical_events.date >= age_threshold) + .sort_by(clinical_events.date) + .last_for_patient() + ) diff --git a/analysis/dataset_definition/variables_cohorts.py b/analysis/dataset_definition/variables_cohorts.py index 803c414..11d83b5 100644 --- a/analysis/dataset_definition/variables_cohorts.py +++ b/analysis/dataset_definition/variables_cohorts.py @@ -36,6 +36,7 @@ matching_death_before, filter_codes_by_category, get_latest_ethnicity, + most_recent_bmi ) # Define generate variables function @@ -194,6 +195,21 @@ def generate_variables(index_date, end_date_exp, end_date_out): ### Sex cov_cat_sex = patients.sex + ## BMI + # bmi_primis -> 99.95% missing + # bmi_stage_primis -> 99.03% missing + # bmi_numeric_nhsd -> 99.67% missing + cov_num_bmi = last_matching_event_clinical_snomed_before( + bmi_stage_primis, index_date + ).numeric_value + + # # 99.93% missing + # cov_num_bmi = most_recent_bmi( + # where=clinical_events.date.is_on_or_between(index_date - days(2 * 366), index_date), + # patient_date_of_birth = patients.date_of_birth, + # minimum_age_at_measurement = 16, + # ).numeric_value + ### Ethnicity ### Grouping refers to the number of pre-defined categories (6 or 16) (White, Mixed, etc...) ### See https://www.opencodelists.org/codelist/opensafely/ethnicity-snomed-0removed/22911876/ @@ -534,6 +550,7 @@ def generate_variables(index_date, end_date_exp, end_date_out): ### Core covariates cov_num_age = cov_num_age, cov_cat_sex = cov_cat_sex, + cov_num_bmi = cov_num_bmi, cov_cat_ethnicity = cov_cat_ethnicity, cov_cat_imd = cov_cat_imd, cov_cat_smoking = cov_cat_smoking, diff --git a/analysis/generate_subsample/generate_subsample.R b/analysis/generate_subsample/generate_subsample.R index 5520890..ea6766a 100644 --- a/analysis/generate_subsample/generate_subsample.R +++ b/analysis/generate_subsample/generate_subsample.R @@ -74,68 +74,124 @@ if (length(args) == 0) { # Load data -------------------------------------------------------------------- print("Load data") -df <- readr::read_rds(paste0( +df_ami <- readr::read_rds(paste0( "output/dataset_clean/input_", cohort, - "_clean.rds" + "_clean_ami.rds" +)) + +df_sahhs <- readr::read_rds(paste0( + "output/dataset_clean/input_", + cohort, + "_clean_sahhs.rds" )) # Sanity check all covariate data types ---------------------------------------- print("Sanity check all covariate data types") -df$cov_bin_ami <- as.factor(df$cov_bin_ami) # outcome -df$cov_bin_sahhs <- as.factor(df$cov_bin_sahhs) # outcome -df$cov_bin_covid <- as.factor(df$cov_bin_covid) # exposure - -df$cov_num_age <- as.numeric(df$cov_num_age) -df$cov_cat_sex <- as.factor(df$cov_cat_sex) -df$cov_cat_ethnicity <- as.factor(df$cov_cat_ethnicity) -df$cov_cat_imd <- as.factor(df$cov_cat_imd) -df$cov_cat_smoking <- as.factor(df$cov_cat_smoking) - -df$cov_bin_carehome <- as.factor(df$cov_bin_carehome) -df$cov_bin_hcworker <- as.factor(df$cov_bin_hcworker) -df$cov_bin_dementia <- as.factor(df$cov_bin_dementia) -df$cov_bin_liver_disease <- as.factor(df$cov_bin_liver_disease) -df$cov_bin_ckd <- as.factor(df$cov_bin_ckd) - -df$cov_bin_cancer <- as.factor(df$cov_bin_cancer) -df$cov_bin_hypertension <- as.factor(df$cov_bin_hypertension) -df$cov_bin_diabetes <- as.factor(df$cov_bin_diabetes) -df$cov_bin_obesity <- as.factor(df$cov_bin_obesity) -df$cov_bin_copd <- as.factor(df$cov_bin_copd) - -df$cov_bin_depression <- as.factor(df$cov_bin_depression) -df$cov_bin_stroke_all <- as.factor(df$cov_bin_stroke_all) -df$cov_bin_other_ae <- as.factor(df$cov_bin_other_ae) -df$cov_bin_vte <- as.factor(df$cov_bin_vte) -df$cov_bin_hf <- as.factor(df$cov_bin_hf) +df_ami$cov_bin_ami <- as.factor(df_ami$cov_bin_ami) # outcome +df_ami$cov_bin_sahhs <- as.factor(df_ami$cov_bin_sahhs) # outcome +df_ami$cov_bin_covid <- as.factor(df_ami$cov_bin_covid) # exposure + +df_ami$cov_num_age <- as.numeric(df_ami$cov_num_age) +df_ami$cov_cat_sex <- as.factor(df_ami$cov_cat_sex) +df_ami$cov_num_bmi <- as.factor(df_ami$cov_num_bmi) +df_ami$cov_cat_ethnicity <- as.factor(df_ami$cov_cat_ethnicity) +df_ami$cov_cat_imd <- as.factor(df_ami$cov_cat_imd) +df_ami$cov_cat_smoking <- as.factor(df_ami$cov_cat_smoking) + +df_ami$cov_bin_carehome <- as.factor(df_ami$cov_bin_carehome) +df_ami$cov_bin_hcworker <- as.factor(df_ami$cov_bin_hcworker) +df_ami$cov_bin_dementia <- as.factor(df_ami$cov_bin_dementia) +df_ami$cov_bin_liver_disease <- as.factor(df_ami$cov_bin_liver_disease) +df_ami$cov_bin_ckd <- as.factor(df_ami$cov_bin_ckd) + +df_ami$cov_bin_cancer <- as.factor(df_ami$cov_bin_cancer) +df_ami$cov_bin_hypertension <- as.factor(df_ami$cov_bin_hypertension) +df_ami$cov_bin_diabetes <- as.factor(df_ami$cov_bin_diabetes) +df_ami$cov_bin_obesity <- as.factor(df_ami$cov_bin_obesity) +df_ami$cov_bin_copd <- as.factor(df_ami$cov_bin_copd) + +df_ami$cov_bin_depression <- as.factor(df_ami$cov_bin_depression) +df_ami$cov_bin_stroke_all <- as.factor(df_ami$cov_bin_stroke_all) +df_ami$cov_bin_other_ae <- as.factor(df_ami$cov_bin_other_ae) +df_ami$cov_bin_vte <- as.factor(df_ami$cov_bin_vte) +df_ami$cov_bin_hf <- as.factor(df_ami$cov_bin_hf) + +df_ami$cov_bin_angina <- as.factor(df_ami$cov_bin_angina) +df_ami$cov_bin_lipidmed <- as.factor(df_ami$cov_bin_lipidmed) +df_ami$cov_bin_antiplatelet <- as.factor(df_ami$cov_bin_antiplatelet) +df_ami$cov_bin_anticoagulant <- as.factor(df_ami$cov_bin_anticoagulant) +df_ami$cov_bin_cocp <- as.factor(df_ami$cov_bin_cocp) + +df_ami$cov_bin_hrt <- as.factor(df_ami$cov_bin_hrt) +df_ami$strat_cat_region <- as.factor(df_ami$strat_cat_region) + + +df_sahhs$cov_bin_sahhs <- as.factor(df_sahhs$cov_bin_sahhs) # outcome +df_sahhs$cov_bin_sahhs <- as.factor(df_sahhs$cov_bin_sahhs) # outcome +df_sahhs$cov_bin_covid <- as.factor(df_sahhs$cov_bin_covid) # exposure + +df_sahhs$cov_num_age <- as.numeric(df_sahhs$cov_num_age) +df_sahhs$cov_cat_sex <- as.factor(df_sahhs$cov_cat_sex) +df_sahhs$cov_num_bmi <- as.factor(df_sahhs$cov_num_bmi) +df_sahhs$cov_cat_ethnicity <- as.factor(df_sahhs$cov_cat_ethnicity) +df_sahhs$cov_cat_imd <- as.factor(df_sahhs$cov_cat_imd) +df_sahhs$cov_cat_smoking <- as.factor(df_sahhs$cov_cat_smoking) + +df_sahhs$cov_bin_carehome <- as.factor(df_sahhs$cov_bin_carehome) +df_sahhs$cov_bin_hcworker <- as.factor(df_sahhs$cov_bin_hcworker) +df_sahhs$cov_bin_dementia <- as.factor(df_sahhs$cov_bin_dementia) +df_sahhs$cov_bin_liver_disease <- as.factor(df_sahhs$cov_bin_liver_disease) +df_sahhs$cov_bin_ckd <- as.factor(df_sahhs$cov_bin_ckd) + +df_sahhs$cov_bin_cancer <- as.factor(df_sahhs$cov_bin_cancer) +df_sahhs$cov_bin_hypertension <- as.factor(df_sahhs$cov_bin_hypertension) +df_sahhs$cov_bin_diabetes <- as.factor(df_sahhs$cov_bin_diabetes) +df_sahhs$cov_bin_obesity <- as.factor(df_sahhs$cov_bin_obesity) +df_sahhs$cov_bin_copd <- as.factor(df_sahhs$cov_bin_copd) + +df_sahhs$cov_bin_depression <- as.factor(df_sahhs$cov_bin_depression) +df_sahhs$cov_bin_stroke_all <- as.factor(df_sahhs$cov_bin_stroke_all) +df_sahhs$cov_bin_other_ae <- as.factor(df_sahhs$cov_bin_other_ae) +df_sahhs$cov_bin_vte <- as.factor(df_sahhs$cov_bin_vte) +df_sahhs$cov_bin_hf <- as.factor(df_sahhs$cov_bin_hf) + +df_sahhs$cov_bin_angina <- as.factor(df_sahhs$cov_bin_angina) +df_sahhs$cov_bin_lipidmed <- as.factor(df_sahhs$cov_bin_lipidmed) +df_sahhs$cov_bin_antiplatelet <- as.factor(df_sahhs$cov_bin_antiplatelet) +df_sahhs$cov_bin_anticoagulant <- as.factor(df_sahhs$cov_bin_anticoagulant) +df_sahhs$cov_bin_cocp <- as.factor(df_sahhs$cov_bin_cocp) + +df_sahhs$cov_bin_hrt <- as.factor(df_sahhs$cov_bin_hrt) +df_sahhs$strat_cat_region <- as.factor(df_sahhs$strat_cat_region) -df$cov_bin_angina <- as.factor(df$cov_bin_angina) -df$cov_bin_lipidmed <- as.factor(df$cov_bin_lipidmed) -df$cov_bin_antiplatelet <- as.factor(df$cov_bin_antiplatelet) -df$cov_bin_anticoagulant <- as.factor(df$cov_bin_anticoagulant) -df$cov_bin_cocp <- as.factor(df$cov_bin_cocp) - -df$cov_bin_hrt <- as.factor(df$cov_bin_hrt) -df$strat_cat_region <- as.factor(df$strat_cat_region) # Generate 10% subsample ------------------------------------------------------ print("Generate 10% subsample") set.seed(2026) # fixed for reproducibility, no overlapping RNG sequences so fine to handle in this way -sample_size <- nrow(df) -selection <- sample(x = c(1:sample_size), size = ceiling(sample_size/10), replace = FALSE) -subsample_df <- df[selection, ] + +sample_size <- nrow(df_ami) # same for both +selection <- sample(x = c(1:sample_size), size = ceiling(sample_size/10), replace = FALSE) # sxame for both + +subsample_df_ami <- df_ami[selection, ] +subsample_df_sahhs <- df_sahhs[selection, ] # Save subsample -------------------------------------------------------------- print("Save subsample") saveRDS( - subsample_df, - file = paste0(generate_subsample_dir, "input_", cohort, "_clean_subsample.rds"), + subsample_df_ami, + file = paste0(generate_subsample_dir, "input_", cohort, "_clean_subsample_ami.rds"), + compress = TRUE +) + +saveRDS( + subsample_df_sahhs, + file = paste0(generate_subsample_dir, "input_", cohort, "_clean_subsample_sahhs.rds"), compress = TRUE ) diff --git a/analysis/lasso_X_var_selection/lasso_X_var_selection.R b/analysis/lasso_X_var_selection/lasso_X_var_selection.R index 0505b7c..816db67 100644 --- a/analysis/lasso_X_var_selection/lasso_X_var_selection.R +++ b/analysis/lasso_X_var_selection/lasso_X_var_selection.R @@ -78,11 +78,21 @@ preex_string <- "" # Load subsample data ---------------------------------------------------------- print("Load subsample data") -df <- readr::read_rds(paste0( - "output/generate_subsample/input_", - cohort, - "_clean_subsample.rds" -)) +if (grepl("ami", name)) { + # subsample + df <- readr::read_rds(paste0( + "output/generate_subsample/input_", + cohort, + "_clean_subsample_ami.rds" + )) +} else { + # subsample + df <- readr::read_rds(paste0( + "output/generate_subsample/input_", + cohort, + "_clean_subsample_sahhs.rds" + )) +} # Data preparation for fully adjusted logistic model --------------------- @@ -93,6 +103,7 @@ df2 <- (df %>% select(c( cov_num_age, cov_cat_sex, + cov_num_bmi, cov_cat_ethnicity, cov_cat_imd, cov_cat_smoking, @@ -125,6 +136,9 @@ df2 <- (df %>% select(c( strat_cat_region ))) +# prevent conversion to factor +df2$cov_num_bmi <- as.numeric(df2$cov_num_bmi) + # Data preparation for the lasso_X logistic model ------------------------------ message("Data preparation for the lasso_X logistic model") @@ -134,6 +148,7 @@ message("Data preparation for the lasso_X logistic model") lasso_X_conf_matrix <- (df %>% select(c( cov_num_age, cov_cat_sex, + cov_num_bmi, cov_cat_ethnicity, cov_cat_imd, cov_cat_smoking, @@ -166,6 +181,9 @@ lasso_X_conf_matrix <- (df %>% select(c( strat_cat_region ))) +# prevent conversion to factor +lasso_X_conf_matrix$cov_num_bmi <- as.numeric(lasso_X_conf_matrix$cov_num_bmi) + lasso_X_conf_matrix_preserving_factors <- model.matrix( ~ ., # formula meaning take all terms data = lasso_X_conf_matrix @@ -188,6 +206,7 @@ cv_lasso_X_model <- cv.glmnet(x = lasso_X_conf_matrix_preserving_factors, y = lasso_X_exposure_matrix_preserving_factors, nfolds = 20, # number of cv datasets family = "binomial", # logistic regression + weights = generate_weights(sample_size = nrow(df)), alpha = 1) # LASSO penalty # tune regularisation parameter lambda to minimise cross-validated error (cvm) @@ -196,6 +215,7 @@ lambda <- cv_lasso_X_model$lambda.min lasso_X_model <- glmnet(x = lasso_X_conf_matrix_preserving_factors, y = lasso_X_exposure_matrix_preserving_factors, family = "binomial", # logistic regression + weights = generate_weights(sample_size = nrow(df)), alpha = 1, # LASSO penalty lambda = lambda) # optimal lambda @@ -208,7 +228,8 @@ fully_adjusted_formula <- "cov_bin_covid ~ ." fully_adjusted_logistic <- glm( fully_adjusted_formula, family = "binomial", - data = df2 + data = df2, + weights = generate_weights(sample_size = nrow(df2)) ) @@ -241,7 +262,7 @@ write.csv( write.csv( lasso_X_coefs, paste0(lasso_X_var_selection_dir, "lasso_X_var_selection-coefs-", name, preex_string, ".csv"), - row.names = FALSE + row.names = TRUE ) write.csv( diff --git a/analysis/lasso_var_selection/lasso_var_selection.R b/analysis/lasso_var_selection/lasso_var_selection.R index 8b1a3e2..c64b5d8 100644 --- a/analysis/lasso_var_selection/lasso_var_selection.R +++ b/analysis/lasso_var_selection/lasso_var_selection.R @@ -86,12 +86,21 @@ preex_string <- "" # Load subsample data ---------------------------------------------------------- print("Load subsample data") -# subsample -df <- readr::read_rds(paste0( - "output/generate_subsample/input_", - cohort, - "_clean_subsample.rds" -)) +if (grepl("ami", name)) { + # subsample + df <- readr::read_rds(paste0( + "output/generate_subsample/input_", + cohort, + "_clean_subsample_ami.rds" + )) +} else { + # subsample + df <- readr::read_rds(paste0( + "output/generate_subsample/input_", + cohort, + "_clean_subsample_sahhs.rds" + )) +} # subsample model_input_df <- readr::read_rds(paste0( @@ -114,6 +123,7 @@ df2 <- (model_input_df %>% select(c( cov_num_age, cov_cat_sex, + cov_num_bmi, cov_cat_ethnicity, cov_cat_imd, cov_cat_smoking, @@ -146,6 +156,9 @@ df2 <- (model_input_df %>% select(c( strat_cat_region ))) +# prevent conversion to factor +df2$cov_num_bmi <- as.numeric(df2$cov_num_bmi) + # Data preparation for lasso cox model ---------------------------------- message("Data preparation for lasso cox model") @@ -155,6 +168,7 @@ lasso_cox_conf_matrix <- (model_input_df %>% select(c( cov_num_age, cov_cat_sex, + cov_num_bmi, cov_cat_ethnicity, cov_cat_imd, cov_cat_smoking, @@ -187,6 +201,9 @@ lasso_cox_conf_matrix <- (model_input_df %>% select(c( strat_cat_region ))) +# prevent conversion to factor +lasso_cox_conf_matrix$cov_num_bmi <- as.numeric(lasso_cox_conf_matrix$cov_num_bmi) + lasso_cox_conf_matrix_preserving_factors <- model.matrix( ~ ., # formula meaning take all terms data = lasso_cox_conf_matrix @@ -219,11 +236,12 @@ lasso_cox_outcome_survival <- Surv(time = as.numeric(outcome_cox_dates), # Fitting the lasso cox model ---------------------------------------------------- message("Fitting the lasso cox model") -cv_lasso_cox_model <- cv.glmnet(x = lasso_cox_conf_matrix_preserving_factors, - y = lasso_cox_outcome_survival, - nfolds = 20, # number of cv datasets - family = "cox", # cox regression - alpha = 1) # LASSO penalty +cv_lasso_cox_model <- cv.glmnet(x = lasso_cox_conf_matrix_preserving_factors, + y = lasso_cox_outcome_survival, + nfolds = 20, # number of cv datasets + family = "cox", # cox regression + weights = generate_weights(sample_size = nrow(model_input_df)), + alpha = 1) # LASSO penalty # tune regularisation parameter lambda to minimise cross-validated error (cvm) lambda <- cv_lasso_cox_model$lambda.min @@ -231,6 +249,7 @@ lambda <- cv_lasso_cox_model$lambda.min lasso_cox_model <- glmnet(x = lasso_cox_conf_matrix_preserving_factors, y = lasso_cox_outcome_survival, family = "cox", # cox regression + weights = generate_weights(sample_size = nrow(model_input_df)), alpha = 1, # LASSO penalty lambda = lambda) # optimal lambda @@ -258,7 +277,8 @@ fully_adjusted_outcome_regression_formula <- make_outcome_formula( fully_adjusted_cox <- coxph( formula = as.formula(fully_adjusted_outcome_regression_formula), - data = df2 + data = df2, + weights = generate_weights(sample_size = nrow(model_input_df)) ) diff --git a/analysis/make_output/make_other_output.R b/analysis/make_output/make_other_output.R index ba6c1e8..e5576f2 100644 --- a/analysis/make_output/make_other_output.R +++ b/analysis/make_output/make_other_output.R @@ -1,6 +1,6 @@ # ------------------------------------------------------------------------------ # -# make_other_output.R +# make_table2_output.R # # This file handles any other required post-processing of tables 1 and 2 # (patient characteristics and outcomes) @@ -112,39 +112,6 @@ for (i in cohorts) { df <- df[df["cohort"] != TRUE, ] -# table1-specific processing --------------------------------------------------- -if (output == "table1") { - print("table1 processing") - df <- pivot_wider( - df, - names_from = "cohort", - values_from = c( - "N [midpoint6_derived]", - "(%) [midpoint6_derived]", - "COVID-19 diagnoses [midpoint6]", - "percent_exposed_midpoint6" - ), - names_vary = "slowest" - ) - - # ensure same decimal places for n_percent_midpoint6 column - n_percent_midpoint6 <- unname(unlist(df[4])) - for (i in seq_along(n_percent_midpoint6)) { - if ((!is.na(n_percent_midpoint6[i])) && (!grepl(".", n_percent_midpoint6[i], fixed = TRUE))) { - n_percent_midpoint6[i] <- paste0(str_sub(n_percent_midpoint6[i], end=-2), ".0%") - } - } - df[4] <- n_percent_midpoint6 - - # ensure same decimal places for percent_exposed_midpoint6_prevax column - percent_exposed_midpoint6_prevax <- unname(unlist(df[6])) - for (i in seq_along(percent_exposed_midpoint6_prevax)) { - if ((!is.na(percent_exposed_midpoint6_prevax[i])) && (!grepl(".", percent_exposed_midpoint6_prevax[i], fixed = TRUE))) { - percent_exposed_midpoint6_prevax[i] <- paste0(str_sub(percent_exposed_midpoint6_prevax[i], end=-2), ".0%") - } - } - df[6] <- percent_exposed_midpoint6_prevax -} # table2-specific processing --------------------------------------------------- if (output == "table2") { @@ -164,4 +131,4 @@ readr::write_csv( df, paste0(makeout_dir, output, sub_str, "_output_midpoint6.csv"), na = "-" -) +) \ No newline at end of file diff --git a/analysis/make_output/make_table1_output.R b/analysis/make_output/make_table1_output.R new file mode 100644 index 0000000..88e88fb --- /dev/null +++ b/analysis/make_output/make_table1_output.R @@ -0,0 +1,258 @@ +# ------------------------------------------------------------------------------ +# +# make_table1_output.R +# +# This file handles any other required post-processing of tables 1 and 2 +# (patient characteristics and outcomes) +# for the replicated Covid-19 x Cardiovascular paper +# Original text: https://doi.org/10.1038/s41467-024-46497-0 +# +# Arguments: +# - output - string, spoecifies which table is being processed +# (table1, table2, etc) +# - cohort - string, specifies which study cohort is being processed +# (prevax, vax, unvax) +# - subgroup - string, specifies other subgroups by which the study population +# can be divided (e.g. main, covid hospitalisation) +# +# Returns: +# - The post-processed patient characteristics data table +# (output/make_output/table1_output_midpoint6.csv) +# - The post-processed outcomes table +# (output/make_output/table2-sub_covidhospital_output_midpoint.csv) +# +# Authors: Emma Tarmey, Venexia Walker, UoB ehrQL Team +# +# ------------------------------------------------------------------------------ + + +# Load packages ---------------------------------------------------------------- +print('Load packages') + +library(magrittr) +library(data.table) +library(stringr) +library(tidyr) + + +# Source common functions ------------------------------------------------------ +print('Source common functions') + +source("analysis/utility.R") + + +# Define make_output folder ------------------------------------------ +print("Creating output/make_output output folder") + +makeout_dir <- "output/make_output/" +fs::dir_create(here::here(makeout_dir)) + + +# Specify arguments ------------------------------------------------------------ +print('Specify arguments') + +args <- commandArgs(trailingOnly = TRUE) + +if (length(args) == 0) { + output <- "table1" # the action to apply + cohorts <- "prevax-preex_FALSE;vax-preex_FALSE;unvax-preex_FALSE;prevax-preex_TRUE;vax-preex_TRUE;unvax-preex_TRUE" # The iterative label + subgroup <- "" +} else { + output <- args[[1]] + cohorts <- args[[2]] + if (length(args) < 3) { + subgroup <- "" # an optional subgroup label (e.g. preex_FALSE) + } else { + subgroup <- args[[3]] + } +} + + +# Separate cohorts ------------------------------------------------------------- +print('Separate cohorts') + +cohorts <- stringr::str_split(as.vector(cohorts), ";")[[1]] + + +# Generate output/saving string ------------------------------------------------ +print('Generate strings') + +if (subgroup == "All" | subgroup == "") { + sub_str <- "" +} else { + if (grepl("preex", subgroup)) { + sub_str <- paste0("-", subgroup) + } else { + sub_str <- paste0("-sub_", subgroup) + } +} + + +# Create blank table ----------------------------------------------------------- +print('Create blank table') + +df_ami <- NULL +df_sahhs <- NULL + + +# Add output from each cohort -------------------------------------------------- +print('Add output from each cohort') + +for (i in cohorts) { + # load input + tmp_ami <- readr::read_csv(paste0( + "output/", + output, + "/", + output, + "-cohort_", + i, + sub_str, + "_ami-midpoint6.csv" + )) + + # create column for cohort + tmp_ami$cohort <- i + + # combining dataframes + df_ami <- rbind(df_ami, tmp_ami, fill = TRUE) +} + +df_ami <- df_ami[df_ami["cohort"] != TRUE, ] + +for (i in cohorts) { + # load input + tmp_sahhs <- readr::read_csv(paste0( + "output/", + output, + "/", + output, + "-cohort_", + i, + sub_str, + "_sahhs-midpoint6.csv" + )) + + # create column for cohort + tmp_sahhs$cohort <- i + + # combining dataframes + df_sahhs <- rbind(df_sahhs, tmp_sahhs, fill = TRUE) +} + +df_sahhs <- df_sahhs[df_sahhs["cohort"] != TRUE, ] + + +# table1-specific processing --------------------------------------------------- +if (output == "table1") { + print("table1 processing (ami)") + + df_ami <- pivot_wider( + df_ami, + names_from = "cohort", + values_from = c( + "N [midpoint6_derived]", + "(%) [midpoint6_derived]", + "COVID-19 diagnoses [midpoint6]", + "percent_exposed_midpoint6" + ), + names_vary = "slowest" + ) + + # ensure same decimal places for n_percent_midpoint6 column + n_percent_midpoint6 <- unname(unlist(df_ami[4])) + for (i in seq_along(n_percent_midpoint6)) { + if ((!is.na(n_percent_midpoint6[i])) && (!grepl(".", n_percent_midpoint6[i], fixed = TRUE))) { + n_percent_midpoint6[i] <- paste0(str_sub(n_percent_midpoint6[i], end=-2), ".0%") + } + } + df_ami[4] <- n_percent_midpoint6 + + # ensure same decimal places for percent_exposed_midpoint6_prevax column + percent_exposed_midpoint6_prevax <- unname(unlist(df_ami[6])) + for (i in seq_along(percent_exposed_midpoint6_prevax)) { + if ((!is.na(percent_exposed_midpoint6_prevax[i])) && (!grepl(".", percent_exposed_midpoint6_prevax[i], fixed = TRUE))) { + percent_exposed_midpoint6_prevax[i] <- paste0(str_sub(percent_exposed_midpoint6_prevax[i], end=-2), ".0%") + } + } + df_ami[6] <- percent_exposed_midpoint6_prevax +} + +if (output == "table1") { + print("table1 processing (sahhs)") + + df_sahhs <- pivot_wider( + df_sahhs, + names_from = "cohort", + values_from = c( + "N [midpoint6_derived]", + "(%) [midpoint6_derived]", + "COVID-19 diagnoses [midpoint6]", + "percent_exposed_midpoint6" + ), + names_vary = "slowest" + ) + + # ensure same decimal places for n_percent_midpoint6 column + n_percent_midpoint6 <- unname(unlist(df_sahhs[4])) + for (i in seq_along(n_percent_midpoint6)) { + if ((!is.na(n_percent_midpoint6[i])) && (!grepl(".", n_percent_midpoint6[i], fixed = TRUE))) { + n_percent_midpoint6[i] <- paste0(str_sub(n_percent_midpoint6[i], end=-2), ".0%") + } + } + df_sahhs[4] <- n_percent_midpoint6 + + # ensure same decimal places for percent_exposed_midpoint6_prevax column + percent_exposed_midpoint6_prevax <- unname(unlist(df_sahhs[6])) + for (i in seq_along(percent_exposed_midpoint6_prevax)) { + if ((!is.na(percent_exposed_midpoint6_prevax[i])) && (!grepl(".", percent_exposed_midpoint6_prevax[i], fixed = TRUE))) { + percent_exposed_midpoint6_prevax[i] <- paste0(str_sub(percent_exposed_midpoint6_prevax[i], end=-2), ".0%") + } + } + df_sahhs[6] <- percent_exposed_midpoint6_prevax +} + + +# Nelson-Aalen estimates ------------------------------------------------------- +print("Nelson-Aalen estimates") + +# ami estimate +df_ami_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_prevax_ami.rds")) +nelsonaalen_ami_summary <- summary(df_ami_nelsonaalen$H0) + +nelson_aalen_mean_ami_row <- which(df_ami$Characteristic == "Nelson-Aalen Estimator (H0)" & df_ami$Subcharacteristic == "Mean") +nelson_aalen_median_ami_row <- which(df_ami$Characteristic == "Nelson-Aalen Estimator (H0)" & df_ami$Subcharacteristic == "Median") +nelson_aalen_IQR_ami_row <- which(df_ami$Characteristic == "Nelson-Aalen Estimator (H0)" & df_ami$Subcharacteristic == "IQR") + +df_ami[nelson_aalen_mean_ami_row, "N [midpoint6_derived]_prevax" ] <- as.character(nelsonaalen_ami_summary[["Mean"]]) +df_ami[nelson_aalen_median_ami_row, "N [midpoint6_derived]_prevax" ] <- as.character(nelsonaalen_ami_summary[["Median"]]) +df_ami[nelson_aalen_IQR_ami_row, "N [midpoint6_derived]_prevax" ] <- as.character((nelsonaalen_ami_summary[["3rd Qu."]] - nelsonaalen_ami_summary[["1st Qu."]])) + + +# sahhs estimate +df_sahhs_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_prevax_sahhs.rds")) +nelsonaalen_sahhs_summary <- summary(df_sahhs_nelsonaalen$H0) + +nelson_aalen_mean_sahhs_row <- which(df_sahhs$Characteristic == "Nelson-Aalen Estimator (H0)" & df_sahhs$Subcharacteristic == "Mean") +nelson_aalen_median_sahhs_row <- which(df_sahhs$Characteristic == "Nelson-Aalen Estimator (H0)" & df_sahhs$Subcharacteristic == "Median") +nelson_aalen_IQR_sahhs_row <- which(df_sahhs$Characteristic == "Nelson-Aalen Estimator (H0)" & df_sahhs$Subcharacteristic == "IQR") + +df_sahhs[nelson_aalen_mean_sahhs_row, "N [midpoint6_derived]_prevax" ] <- as.character(nelsonaalen_sahhs_summary[["Mean"]]) +df_sahhs[nelson_aalen_median_sahhs_row, "N [midpoint6_derived]_prevax" ] <- as.character(nelsonaalen_sahhs_summary[["Median"]]) +df_sahhs[nelson_aalen_IQR_sahhs_row, "N [midpoint6_derived]_prevax" ] <- as.character((nelsonaalen_sahhs_summary[["3rd Qu."]] - nelsonaalen_sahhs_summary[["1st Qu."]])) + + +# Save output ------------------------------------------------------------------ +print('Save output') + +readr::write_csv( + df_ami, + paste0(makeout_dir, output, sub_str, "_output_ami_midpoint6.csv"), + na = "-" +) + +readr::write_csv( + df_sahhs, + paste0(makeout_dir, output, sub_str, "_output_sahhs_midpoint6.csv"), + na = "-" +) diff --git a/analysis/model/fn-prepare_model_input.R b/analysis/model/fn-prepare_model_input.R index ea6fe83..0b3ce4c 100644 --- a/analysis/model/fn-prepare_model_input.R +++ b/analysis/model/fn-prepare_model_input.R @@ -16,11 +16,19 @@ prepare_model_input <- function(name) { # Load data ------------------------------------------------------------------ print(paste0("Load data for ", active_analyses$name)) - input <- readr::read_rds(paste0( - "output/dataset_clean/input_", - active_analyses$cohort, - "_clean.rds" - )) + if (grepl("ami", name)) { + input <- readr::read_rds(paste0( + "output/dataset_clean/input_", + active_analyses$cohort, + "_clean_ami.rds" + )) + } else { + input <- readr::read_rds(paste0( + "output/dataset_clean/input_", + active_analyses$cohort, + "_clean_sahhs.rds" + )) + } # Restrict to required variables for dataset preparation --------------------- print("Restrict to required variables for dataset preparation") @@ -35,6 +43,7 @@ prepare_model_input <- function(name) { active_analyses$strata, active_analyses$covariate_age, "cov_cat_sex", + "cov_num_bmi", "cov_cat_ethnicity", "cov_cat_smoking", "cov_bin_covid", @@ -114,7 +123,7 @@ prepare_model_input <- function(name) { ) %>% dplyr::ungroup() - keep <- c(keep, "cov_bin_covid", "cov_bin_sahhs") + keep <- c(keep, "cov_num_bmi", "cov_bin_covid", "cov_bin_sahhs") return(list(input = input, keep = keep)) -} +} \ No newline at end of file diff --git a/analysis/model/fn-prepare_model_input_subsample.R b/analysis/model/fn-prepare_model_input_subsample.R index 72bb99c..0b83193 100644 --- a/analysis/model/fn-prepare_model_input_subsample.R +++ b/analysis/model/fn-prepare_model_input_subsample.R @@ -16,12 +16,21 @@ prepare_model_input_subsample <- function(name) { # Load subsample data --------------------------------------------------------- print(paste0("Load subsample data for ", active_analyses$name)) - # subsample in place of study population - input <- readr::read_rds(paste0( - "output/generate_subsample/input_", - active_analyses$cohort, - "_clean_subsample.rds" - )) + if (grepl("ami", name)) { + # subsample + input <- readr::read_rds(paste0( + "output/generate_subsample/input_", + active_analyses$cohort, + "_clean_subsample_ami.rds" + )) + } else { + # subsample + input <- readr::read_rds(paste0( + "output/generate_subsample/input_", + active_analyses$cohort, + "_clean_subsample_sahhs.rds" + )) + } # Restrict to required variables for dataset preparation --------------------- print("Restrict to required variables for dataset preparation") @@ -36,6 +45,7 @@ prepare_model_input_subsample <- function(name) { active_analyses$strata, active_analyses$covariate_age, "cov_cat_sex", + "cov_num_bmi", "cov_cat_ethnicity", "cov_cat_smoking", "cov_bin_covid", @@ -115,7 +125,7 @@ prepare_model_input_subsample <- function(name) { ) %>% dplyr::ungroup() - keep <- c(keep, "cov_bin_covid", "cov_bin_sahhs") + keep <- c(keep, "cov_num_bmi", "cov_bin_covid", "cov_bin_sahhs") return(list(input = input, keep = keep)) -} +} \ No newline at end of file diff --git a/analysis/nelson_aalen_plots/nelson_aalen_plots.R b/analysis/nelson_aalen_plots/nelson_aalen_plots.R new file mode 100644 index 0000000..0f8051d --- /dev/null +++ b/analysis/nelson_aalen_plots/nelson_aalen_plots.R @@ -0,0 +1,213 @@ +# ------------------------------------------------------------------------------ +# +# nelson_aalen_plots.R +# +# This file generates plots of the nelson_aalen estimator for both outcomes +# +# Arguments: +# - cohort - string, defines which of three opensafely cohorts to describe +# (prevax, vax, unvax) +# +# Returns: +# - plots! +# +# Authors: Emma Tarmey +# +# ------------------------------------------------------------------------------ + + +# Load libraries --------------------------------------------------------------- +print("Load libraries") + +library(magrittr) +library(mice) +library(here) +library(dplyr) +library(fs) +library(survival) +library(ggplot2) + + +# Source common functions ------------------------------------------------------ +print("Source common functions") + +source("analysis/utility.R") + + +# Specify arguments ------------------------------------------------------------ +print("Specify arguments") + +args <- commandArgs(trailingOnly = TRUE) +print(length(args)) + +if (length(args) == 0) { + # default argument values + cohort <- "prevax" +} else { + # YAML arguments + cohort <- args[[1]] +} + + +# Load data -------------------------------------------------------------------- +print("Load data") + +df_ami_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_", cohort, "_ami.rds")) +df_sahhs_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_", cohort, "_sahhs.rds")) + + +# Take random 100-point sample ------------------------------------------------- +print("Take random 100-point sample") + +sample_size <- nrow(df_ami_nelsonaalen) +random_sample <- sample(x = c(1:sample_size), size = 100, replace = FALSE) + +df_ami_nelsonaalen <- df_ami_nelsonaalen[random_sample, ] +df_sahhs_nelsonaalen <- df_sahhs_nelsonaalen[random_sample, ] + + +# Sort by x-axis value (time) -------------------------------------------------- +print("Sort by x-axis value (time)") + +df_ami_nelsonaalen <- df_ami_nelsonaalen[order(df_ami_nelsonaalen$time), ] +df_sahhs_nelsonaalen <- df_sahhs_nelsonaalen[order(df_sahhs_nelsonaalen$time), ] + + +# Convert time to date format -------------------------------------------------- +print("Convert time to date format") + +df_ami_nelsonaalen$time <- as.Date(df_ami_nelsonaalen$time, origin = lubridate::origin) +df_sahhs_nelsonaalen$time <- as.Date(df_sahhs_nelsonaalen$time, origin = lubridate::origin) + + +# Calculate 95% Confidence Interval -------------------------------------------- +print("Calculate 95% Confidence Interval") + +# z value is 1.96 for 95% confidence interval +z_value <- qnorm(p=(0.05/2), lower.tail = FALSE) + +df_ami_nelsonaalen$upper_CI <- df_ami_nelsonaalen$H0 + (z_value * df_ami_nelsonaalen$se) +df_ami_nelsonaalen$lower_CI <- df_ami_nelsonaalen$H0 - (z_value * df_ami_nelsonaalen$se) + +df_sahhs_nelsonaalen$upper_CI <- df_sahhs_nelsonaalen$H0 + (z_value * df_sahhs_nelsonaalen$se) +df_sahhs_nelsonaalen$lower_CI <- df_sahhs_nelsonaalen$H0 - (z_value * df_sahhs_nelsonaalen$se) + + +# Generate Nelson-Aalen plots ----- +print("Generate Nelson-Aalen plots") + +png( + paste0("output/dataset_clean/nelson_aalen_", cohort, "_ami_ggplot.png"), + width = 1000, + height = 800 +) +ggplot(df_ami_nelsonaalen, aes(x=time)) + + geom_step(mapping = aes(y=H0, colour = "Nelson Aalen Estimate"), linetype = 1, linewidth = 0.5) + + geom_step(mapping = aes(y=lower_CI, colour = "Lower 95% Confidence Interval"), linetype = 2, linewidth = 0.5) + + geom_step(mapping = aes(y=upper_CI, colour = "Upper 95% Confidence Interval"), linetype = 2, linewidth = 0.5) + + scale_colour_manual( + name = 'Legend', + breaks = c('Nelson Aalen Estimate', 'Upper 95% Confidence Interval', 'Lower 95% Confidence Interval'), + values = c('Nelson Aalen Estimate'='black', 'Upper 95% Confidence Interval'='blue', 'Lower 95% Confidence Interval'='red') + ) + + labs( + title = "Nelson-Aalen Estimator (H0) Plot", + subtitle = "Acute Myocardial Infarction (ami) (sample 100 points)", + x = "Time", + y = "Cumulative hazard (H0)" + ) + + theme(legend.position = "right") +dev.off() + +png( + paste0("output/dataset_clean/nelson_aalen_", cohort, "_ami.png"), + width = 1000, + height = 800 +) +plot( + x = df_ami_nelsonaalen$time, + y = df_ami_nelsonaalen$H0, + type = "s", + main = "Nelson-Aalen Estimator (H0) Plot\nAcute Myocardial Infarction (ami) (sample 100 points)", + xlab = "Time", + ylab = "Cumulative hazard (H0)" +) +lines( + x = df_ami_nelsonaalen$time, + y = df_ami_nelsonaalen$upper_CI, + type = "s", + lty = 2, + col = "blue" +) +lines( + x = df_ami_nelsonaalen$time, + y = df_ami_nelsonaalen$lower_CI, + type = "s", + lty = 2, + col = "red" +) +legend( + "bottomright", + legend = c("Nelson Aalen Estimate", "Upper 95% Confidence Interval", "Lower 95% Confidence Interval"), + col = c("black", "blue", "red"), + lty = c(1, 2, 2) +) +dev.off() + +png( + paste0("output/dataset_clean/nelson_aalen_", cohort, "_sahhs_ggplot.png"), + width = 1000, + height = 800 +) +ggplot(df_sahhs_nelsonaalen, aes(x=time)) + + geom_step(mapping = aes(y=H0, colour = "Nelson Aalen Estimate"), linetype = 1, linewidth = 0.5) + + geom_step(mapping = aes(y=lower_CI, colour = "Lower 95% Confidence Interval"), linetype = 2, linewidth = 0.5) + + geom_step(mapping = aes(y=upper_CI, colour = "Upper 95% Confidence Interval"), linetype = 2, linewidth = 0.5) + + scale_colour_manual( + name = 'Legend', + breaks = c('Nelson Aalen Estimate', 'Upper 95% Confidence Interval', 'Lower 95% Confidence Interval'), + values = c('Nelson Aalen Estimate'='black', 'Upper 95% Confidence Interval'='blue', 'Lower 95% Confidence Interval'='red') + ) + + labs( + title = "Nelson-Aalen Estimator (H0) Plot", + subtitle = "Subarachnoid haemorrhage and haemorrhage stroke (sahhs) (sample 100 points)", + x = "Time", + y = "Cumulative hazard (H0)" + ) + + theme(legend.position = "right") +dev.off() + +png( + paste0("output/dataset_clean/nelson_aalen_", cohort, "_sahhs.png"), + width = 1000, + height = 800 +) +plot( + x = df_sahhs_nelsonaalen$time, + y = df_sahhs_nelsonaalen$H0, + type = "s", + main = "Nelson-Aalen Estimator (H0) Plot\nSubarachnoid haemhorrhage & haemhorrhage stroke (sahhs) (sample 100 points)", + xlab = "Time", + ylab = "Cumulative hazard (H0)" +) +lines( + x = df_sahhs_nelsonaalen$time, + y = df_sahhs_nelsonaalen$upper_CI, + type = "s", + lty = 2, + col = "blue" +) +lines( + x = df_sahhs_nelsonaalen$time, + y = df_sahhs_nelsonaalen$lower_CI, + type = "s", + lty = 2, + col = "red" +) +legend( + "bottomright", + legend = c("Nelson Aalen Estimate", "Upper 95% Confidence Interval", "Lower 95% Confidence Interval"), + col = c("black", "blue", "red"), + lty = c(1, 2, 2) +) +dev.off() diff --git a/analysis/post_hoc_vars/post_hoc_vars.R b/analysis/post_hoc_vars/post_hoc_vars.R deleted file mode 100644 index 3d35955..0000000 --- a/analysis/post_hoc_vars/post_hoc_vars.R +++ /dev/null @@ -1,94 +0,0 @@ -# ------------------------------------------------------------------------------ -# -# post_hoc_vars.R -# -# This file generates any commonly used variables defined after the dataset -# defintion runs, such as indicator variables. -# -# Authors: Emma Tarmey, Venexia Walker, UoB ehrQL Team -# -# ------------------------------------------------------------------------------ - -# Refresh local R session ------------------------------------------------------ -print("Refresh local R session") - -rm(list=ls()) - - -# Specify arguments ------------------------------------------------------------ -print("Specify arguments") - -args <- commandArgs(trailingOnly = TRUE) -if (length(args) == 0) { - cohort <- "prevax" -} else { - cohort <- args[[1]] -} - - -# Load data -------------------------------------------------------------------- -print("Load data") - -df <- readr::read_rds(paste0( - "output/dataset_clean/input_", - cohort, - "_clean_prehoc.rds" -)) - - -# Define variables ------------------------------------------------------------- -print("Define variables") - -df$cov_bin_covid <- !is.na(df$exp_date_covid) -df$cov_bin_sahhs <- !is.na(df$out_date_stroke_sahhs) - - -# Check all covariates types --------------------------------------------------- -print("Check all covariate types") - -df$cov_bin_ami <- as.factor(df$cov_bin_ami) # outcome -df$cov_bin_sahhs <- as.factor(df$cov_bin_sahhs) # outcome -df$cov_bin_covid <- as.factor(df$cov_bin_covid) # exposure - -df$cov_num_age <- as.numeric(df$cov_num_age) -df$cov_cat_sex <- as.factor(df$cov_cat_sex) -df$cov_cat_ethnicity <- as.factor(df$cov_cat_ethnicity) -df$cov_cat_imd <- as.factor(df$cov_cat_imd) -df$cov_cat_smoking <- as.factor(df$cov_cat_smoking) - -df$cov_bin_carehome <- as.factor(df$cov_bin_carehome) -df$cov_bin_hcworker <- as.factor(df$cov_bin_hcworker) -df$cov_bin_dementia <- as.factor(df$cov_bin_dementia) -df$cov_bin_liver_disease <- as.factor(df$cov_bin_liver_disease) -df$cov_bin_ckd <- as.factor(df$cov_bin_ckd) - -df$cov_bin_cancer <- as.factor(df$cov_bin_cancer) -df$cov_bin_hypertension <- as.factor(df$cov_bin_hypertension) -df$cov_bin_diabetes <- as.factor(df$cov_bin_diabetes) -df$cov_bin_obesity <- as.factor(df$cov_bin_obesity) -df$cov_bin_copd <- as.factor(df$cov_bin_copd) - -df$cov_bin_depression <- as.factor(df$cov_bin_depression) -df$cov_bin_stroke_all <- as.factor(df$cov_bin_stroke_all) -df$cov_bin_other_ae <- as.factor(df$cov_bin_other_ae) -df$cov_bin_vte <- as.factor(df$cov_bin_vte) -df$cov_bin_hf <- as.factor(df$cov_bin_hf) - -df$cov_bin_angina <- as.factor(df$cov_bin_angina) -df$cov_bin_lipidmed <- as.factor(df$cov_bin_lipidmed) -df$cov_bin_antiplatelet <- as.factor(df$cov_bin_antiplatelet) -df$cov_bin_anticoagulant <- as.factor(df$cov_bin_anticoagulant) -df$cov_bin_cocp <- as.factor(df$cov_bin_cocp) - -df$cov_bin_hrt <- as.factor(df$cov_bin_hrt) -df$strat_cat_region <- as.factor(df$strat_cat_region) - - -# Save results ----------------------------------------------------------------- -print("Save results") - -saveRDS( - df, - file = paste0("output/dataset_clean/input_", cohort, "_clean.rds"), - compress = TRUE -) diff --git a/analysis/table1/table1.R b/analysis/table1/table1.R index 5a52fc7..23707f2 100644 --- a/analysis/table1/table1.R +++ b/analysis/table1/table1.R @@ -80,52 +80,102 @@ age_bounds <- as.numeric(stringr::str_split(as.vector(age_str), ";")[[1]]) # Load data -------------------------------------------------------------------- print("Load data") -df <- readr::read_rds(paste0( +# full population, ami +df_ami <- readr::read_rds(paste0( "output/dataset_clean/input_", cohort, - "_clean.rds" + "_clean_ami.rds" )) +# full population, sahhs +df_sahhs <- readr::read_rds(paste0( + "output/dataset_clean/input_", + cohort, + "_clean_sahhs.rds" +)) -# Check all covariates -------------------------------------------------------- -message("Check all covariates") - -df$cov_bin_ami <- as.logical(df$cov_bin_ami) # outcome -df$cov_bin_sahhs <- as.logical(df$cov_bin_sahhs) # outcome -df$cov_bin_covid <- as.logical(df$cov_bin_covid) # exposure - -df$cov_num_age <- as.numeric(df$cov_num_age) -df$cov_cat_sex <- as.factor(df$cov_cat_sex) -df$cov_cat_ethnicity <- as.factor(df$cov_cat_ethnicity) -df$cov_cat_imd <- as.factor(df$cov_cat_imd) -df$cov_cat_smoking <- as.factor(df$cov_cat_smoking) - -df$cov_bin_carehome <- as.factor(df$cov_bin_carehome) -df$cov_bin_hcworker <- as.factor(df$cov_bin_hcworker) -df$cov_bin_dementia <- as.factor(df$cov_bin_dementia) -df$cov_bin_liver_disease <- as.factor(df$cov_bin_liver_disease) -df$cov_bin_ckd <- as.factor(df$cov_bin_ckd) - -df$cov_bin_cancer <- as.factor(df$cov_bin_cancer) -df$cov_bin_hypertension <- as.factor(df$cov_bin_hypertension) -df$cov_bin_diabetes <- as.factor(df$cov_bin_diabetes) -df$cov_bin_obesity <- as.factor(df$cov_bin_obesity) -df$cov_bin_copd <- as.factor(df$cov_bin_copd) +# nelson aalen estimator +df_ami_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_", cohort, "_ami.rds")) +df_sahhs_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_", cohort, "_sahhs.rds")) -df$cov_bin_depression <- as.factor(df$cov_bin_depression) -df$cov_bin_stroke_all <- as.factor(df$cov_bin_stroke_all) -df$cov_bin_other_ae <- as.factor(df$cov_bin_other_ae) -df$cov_bin_vte <- as.factor(df$cov_bin_vte) -df$cov_bin_hf <- as.factor(df$cov_bin_hf) -df$cov_bin_angina <- as.factor(df$cov_bin_angina) -df$cov_bin_lipidmed <- as.factor(df$cov_bin_lipidmed) -df$cov_bin_antiplatelet <- as.factor(df$cov_bin_antiplatelet) -df$cov_bin_anticoagulant <- as.factor(df$cov_bin_anticoagulant) -df$cov_bin_cocp <- as.factor(df$cov_bin_cocp) +# Check all covariates -------------------------------------------------------- +message("Check all covariates") -df$cov_bin_hrt <- as.factor(df$cov_bin_hrt) -df$strat_cat_region <- as.factor(df$strat_cat_region) +df_ami$cov_bin_ami <- as.logical(df_ami$cov_bin_ami) # outcome +df_ami$cov_bin_sahhs <- as.logical(df_ami$cov_bin_sahhs) # outcome +df_ami$cov_bin_covid <- as.logical(df_ami$cov_bin_covid) # exposure + +df_ami$cov_num_age <- as.numeric(df_ami$cov_num_age) +df_ami$cov_cat_sex <- as.factor(df_ami$cov_cat_sex) +df_ami$cov_cat_ethnicity <- as.factor(df_ami$cov_cat_ethnicity) +df_ami$cov_cat_imd <- as.factor(df_ami$cov_cat_imd) +df_ami$cov_cat_smoking <- as.factor(df_ami$cov_cat_smoking) + +df_ami$cov_bin_carehome <- as.factor(df_ami$cov_bin_carehome) +df_ami$cov_bin_hcworker <- as.factor(df_ami$cov_bin_hcworker) +df_ami$cov_bin_dementia <- as.factor(df_ami$cov_bin_dementia) +df_ami$cov_bin_liver_disease <- as.factor(df_ami$cov_bin_liver_disease) +df_ami$cov_bin_ckd <- as.factor(df_ami$cov_bin_ckd) + +df_ami$cov_bin_cancer <- as.factor(df_ami$cov_bin_cancer) +df_ami$cov_bin_hypertension <- as.factor(df_ami$cov_bin_hypertension) +df_ami$cov_bin_diabetes <- as.factor(df_ami$cov_bin_diabetes) +df_ami$cov_bin_obesity <- as.factor(df_ami$cov_bin_obesity) +df_ami$cov_bin_copd <- as.factor(df_ami$cov_bin_copd) + +df_ami$cov_bin_depression <- as.factor(df_ami$cov_bin_depression) +df_ami$cov_bin_stroke_all <- as.factor(df_ami$cov_bin_stroke_all) +df_ami$cov_bin_other_ae <- as.factor(df_ami$cov_bin_other_ae) +df_ami$cov_bin_vte <- as.factor(df_ami$cov_bin_vte) +df_ami$cov_bin_hf <- as.factor(df_ami$cov_bin_hf) + +df_ami$cov_bin_angina <- as.factor(df_ami$cov_bin_angina) +df_ami$cov_bin_lipidmed <- as.factor(df_ami$cov_bin_lipidmed) +df_ami$cov_bin_antiplatelet <- as.factor(df_ami$cov_bin_antiplatelet) +df_ami$cov_bin_anticoagulant <- as.factor(df_ami$cov_bin_anticoagulant) +df_ami$cov_bin_cocp <- as.factor(df_ami$cov_bin_cocp) + +df_ami$cov_bin_hrt <- as.factor(df_ami$cov_bin_hrt) +df_ami$strat_cat_region <- as.factor(df_ami$strat_cat_region) + + +df_sahhs$cov_bin_sahhs <- as.logical(df_sahhs$cov_bin_sahhs) # outcome +df_sahhs$cov_bin_sahhs <- as.logical(df_sahhs$cov_bin_sahhs) # outcome +df_sahhs$cov_bin_covid <- as.logical(df_sahhs$cov_bin_covid) # exposure + +df_sahhs$cov_num_age <- as.numeric(df_sahhs$cov_num_age) +df_sahhs$cov_cat_sex <- as.factor(df_sahhs$cov_cat_sex) +df_sahhs$cov_cat_ethnicity <- as.factor(df_sahhs$cov_cat_ethnicity) +df_sahhs$cov_cat_imd <- as.factor(df_sahhs$cov_cat_imd) +df_sahhs$cov_cat_smoking <- as.factor(df_sahhs$cov_cat_smoking) + +df_sahhs$cov_bin_carehome <- as.factor(df_sahhs$cov_bin_carehome) +df_sahhs$cov_bin_hcworker <- as.factor(df_sahhs$cov_bin_hcworker) +df_sahhs$cov_bin_dementia <- as.factor(df_sahhs$cov_bin_dementia) +df_sahhs$cov_bin_liver_disease <- as.factor(df_sahhs$cov_bin_liver_disease) +df_sahhs$cov_bin_ckd <- as.factor(df_sahhs$cov_bin_ckd) + +df_sahhs$cov_bin_cancer <- as.factor(df_sahhs$cov_bin_cancer) +df_sahhs$cov_bin_hypertension <- as.factor(df_sahhs$cov_bin_hypertension) +df_sahhs$cov_bin_diabetes <- as.factor(df_sahhs$cov_bin_diabetes) +df_sahhs$cov_bin_obesity <- as.factor(df_sahhs$cov_bin_obesity) +df_sahhs$cov_bin_copd <- as.factor(df_sahhs$cov_bin_copd) + +df_sahhs$cov_bin_depression <- as.factor(df_sahhs$cov_bin_depression) +df_sahhs$cov_bin_stroke_all <- as.factor(df_sahhs$cov_bin_stroke_all) +df_sahhs$cov_bin_other_ae <- as.factor(df_sahhs$cov_bin_other_ae) +df_sahhs$cov_bin_vte <- as.factor(df_sahhs$cov_bin_vte) +df_sahhs$cov_bin_hf <- as.factor(df_sahhs$cov_bin_hf) + +df_sahhs$cov_bin_angina <- as.factor(df_sahhs$cov_bin_angina) +df_sahhs$cov_bin_lipidmed <- as.factor(df_sahhs$cov_bin_lipidmed) +df_sahhs$cov_bin_antiplatelet <- as.factor(df_sahhs$cov_bin_antiplatelet) +df_sahhs$cov_bin_anticoagulant <- as.factor(df_sahhs$cov_bin_anticoagulant) +df_sahhs$cov_bin_cocp <- as.factor(df_sahhs$cov_bin_cocp) + +df_sahhs$cov_bin_hrt <- as.factor(df_sahhs$cov_bin_hrt) +df_sahhs$strat_cat_region <- as.factor(df_sahhs$strat_cat_region) # Table 1 Processing Start ----------------------------------------------------- @@ -134,12 +184,16 @@ print("Table 1 processing") # Remove people with history of COVID-19 --------------------------------------- print("Remove people with history of COVID-19") -df <- df[df$sub_bin_covidhistory == FALSE, ] +df_ami <- df_ami[df_ami$sub_bin_covidhistory == FALSE, ] + +df_sahhs <- df_sahhs[df_sahhs$sub_bin_covidhistory == FALSE, ] # Create exposure indicator ---------------------------------------------------- print("Create exposure indicator") -df$exposed <- df$cov_bin_covid +df_ami$exposed <- df_ami$cov_bin_covid + +df_sahhs$exposed <- df_sahhs$cov_bin_covid # Select for pre-existing conditions print("Select for pre-existing conditions") @@ -148,76 +202,139 @@ print("Select for pre-existing conditions") # because neither asthma nor copdoutcomes are present preex_string <- "" # if (preex != "All") { -# df <- df[df$sup_bin_preex == preex, ] +# df_ami <- df_ami[df_ami$sup_bin_preex == preex, ] # preex_string <- paste0("-preex_", preex) # } +# preex optional argument deliberately ignored, sup_bin_preex does not exist +# because neither asthma nor copdoutcomes are present +preex_string <- "" +# if (preex != "All") { +# df_sahhs <- df_sahhs[df_sahhs$sup_bin_preex == preex, ] +# preex_string <- paste0("-preex_", preex) +# } + + # Define age groups ------------------------------------------------------------ print("Define age groups") -df$cov_cat_age_group <- numerical_to_categorical(df$cov_num_age, age_bounds) # See utility.R +df_ami$cov_cat_age_group <- numerical_to_categorical(df_ami$cov_num_age, age_bounds) # See utility.R + +# df_ami$cov_cat_consrate2019 <- numerical_to_categorical( +# df_ami$cov_num_consrate2019, +# c(1, 6), +# zero_flag = TRUE +# ) + +median_iqr_age <- create_median_iqr_string(df_ami$cov_num_age) # See utility.R -# df$cov_cat_consrate2019 <- numerical_to_categorical( -# df$cov_num_consrate2019, + +df_sahhs$cov_cat_age_group <- numerical_to_categorical(df_sahhs$cov_num_age, age_bounds) # See utility.R + +# df_sahhs$cov_cat_consrate2019 <- numerical_to_categorical( +# df_sahhs$cov_num_consrate2019, # c(1, 6), # zero_flag = TRUE # ) -median_iqr_age <- create_median_iqr_string(df$cov_num_age) # See utility.R +median_iqr_age <- create_median_iqr_string(df_sahhs$cov_num_age) # See utility.R + # Filter data ------------------------------------------------------------------ print("Filter data") -df <- df[, c( +df_ami <- df_ami[, c( "patient_id", "exposed", - colnames(df)[grepl("cov_cat_", colnames(df))], - colnames(df)[grepl("strat_cat_", colnames(df))], - colnames(df)[grepl("cov_bin_", colnames(df))] + colnames(df_ami)[grepl("cov_cat_", colnames(df_ami))], + colnames(df_ami)[grepl("strat_cat_", colnames(df_ami))], + colnames(df_ami)[grepl("cov_bin_", colnames(df_ami))] )] -df$All <- "All" +df_ami$All <- "All" # Filter binary data -for (colname in colnames(df)[grepl("cov_bin_", colnames(df))]) { - df[[colname]] <- sapply(df[[colname]], as.character) +for (colname in colnames(df_ami)[grepl("cov_bin_", colnames(df_ami))]) { + df_ami[[colname]] <- sapply(df_ami[[colname]], as.character) } -df <- df %>% +df_ami <- df_ami %>% + mutate(across(where(is.factor), as.character)) + +df_sahhs <- df_sahhs[, c( + "patient_id", + "exposed", + colnames(df_sahhs)[grepl("cov_cat_", colnames(df_sahhs))], + colnames(df_sahhs)[grepl("strat_cat_", colnames(df_sahhs))], + colnames(df_sahhs)[grepl("cov_bin_", colnames(df_sahhs))] +)] + +df_sahhs$All <- "All" + +# Filter binary data + +for (colname in colnames(df_sahhs)[grepl("cov_bin_", colnames(df_sahhs))]) { + df_sahhs[[colname]] <- sapply(df_sahhs[[colname]], as.character) +} + +df_sahhs <- df_sahhs %>% mutate(across(where(is.factor), as.character)) # Convert to characteristics and subcharacteristics ---------------------------- print("Convert to characteristics and subcharacteristics") -df <- tidyr::pivot_longer( - df, - cols = setdiff(colnames(df), c("patient_id", "exposed")), +df_ami <- tidyr::pivot_longer( + df_ami, + cols = setdiff(colnames(df_ami), c("patient_id", "exposed")), + names_to = "characteristic", + values_to = "subcharacteristic" +) + +df_ami$total <- 1 + +df_sahhs <- tidyr::pivot_longer( + df_sahhs, + cols = setdiff(colnames(df_sahhs), c("patient_id", "exposed")), names_to = "characteristic", values_to = "subcharacteristic" ) -df$total <- 1 +df_sahhs$total <- 1 # Tidy missing data labels ----------------------------------------------------- print("Tidy missing data labels") -df$subcharacteristic <- ifelse( - df$subcharacteristic == "" | - df$subcharacteristic == "unknown" | - is.na(df$subcharacteristic), +df_ami$subcharacteristic <- ifelse( + df_ami$subcharacteristic == "" | + df_ami$subcharacteristic == "unknown" | + is.na(df_ami$subcharacteristic), "Missing", - df$subcharacteristic + df_ami$subcharacteristic +) + +df_sahhs$subcharacteristic <- ifelse( + df_sahhs$subcharacteristic == "" | + df_sahhs$subcharacteristic == "unknown" | + is.na(df_sahhs$subcharacteristic), + "Missing", + df_sahhs$subcharacteristic ) # Aggregate data --------------------------------------------------------------- print("Aggregate data") -df <- aggregate( +df_ami <- aggregate( + cbind(total, exposed) ~ characteristic + subcharacteristic, + data = df_ami, + sum +) + +df_sahhs <- aggregate( cbind(total, exposed) ~ characteristic + subcharacteristic, - data = df, + data = df_sahhs, sum ) @@ -225,45 +342,133 @@ df <- aggregate( # Sort characteristics --------------------------------------------------------- print("Sort characteristics") -df <- df[order(df$characteristic, df$subcharacteristic), ] +df_ami <- df_ami[order(df_ami$characteristic, df_ami$subcharacteristic), ] + +df_sahhs <- df_sahhs[order(df_sahhs$characteristic, df_sahhs$subcharacteristic), ] + # Add in Median IQR print('Add median (IQR) age') # Pastes: "Mean Age (LQ Age - UQ Age)" as a string for each cohort -df[nrow(df) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, 0) +df_ami[nrow(df_ami) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, 0) +df_sahhs[nrow(df_sahhs) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, 0) + + +# Nelson-Aalen summary statistics --------------------------------------------- +print("Nelson-Aalen summary statistics") + +nelsonaalen_ami_summary <- summary(df_ami_nelsonaalen$H0) +df_ami[nrow(df_ami) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Mean", nelsonaalen_ami_summary[["Mean"]], 0) +df_ami[nrow(df_ami) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Median", nelsonaalen_ami_summary[["Median"]], 0) +df_ami[nrow(df_ami) + 1,] <- c("Nelson-Aalen Estimator (H0)", "IQR", (nelsonaalen_ami_summary[["3rd Qu."]] - nelsonaalen_ami_summary[["1st Qu."]]), 0) + +nelsonaalen_sahhs_summary <- summary(df_sahhs_nelsonaalen$H0) +df_sahhs[nrow(df_sahhs) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Mean", nelsonaalen_sahhs_summary[["Mean"]], 0) +df_sahhs[nrow(df_sahhs) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Median", nelsonaalen_sahhs_summary[["Median"]], 0) +df_sahhs[nrow(df_sahhs) + 1,] <- c("Nelson-Aalen Estimator (H0)", "IQR", (nelsonaalen_sahhs_summary[["3rd Qu."]] - nelsonaalen_sahhs_summary[["1st Qu."]]), 0) + # Save Table 1 ----------------------------------------------------------------- print("Save Table 1") write.csv( - df, - paste0(table1_dir, "table1-cohort_", cohort, preex_string, ".csv"), + df_ami, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_ami.csv"), row.names = FALSE ) +write.csv( + df_sahhs, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_sahhs.csv"), + row.names = FALSE +) + + # Perform redaction ------------------------------------------------------------ print("Perform redaction") -df <- df[df$subcharacteristic != "Median (IQR)", ] # Remove Median IQR row -df <- df[df$subcharacteristic != FALSE, ] # Remove False binary data +df_ami <- df_ami[df_ami$subcharacteristic != "Median (IQR)", ] # Remove Median IQR row +df_ami <- df_ami[df_ami$subcharacteristic != FALSE, ] # Remove False binary data + +df_ami$total_midpoint6 <- roundmid_any(df_ami$total) +df_ami$exposed_midpoint6 <- roundmid_any(df_ami$exposed) + +df_sahhs <- df_sahhs[df_sahhs$subcharacteristic != "Median (IQR)", ] # Remove Median IQR row +df_sahhs <- df_sahhs[df_sahhs$subcharacteristic != FALSE, ] # Remove False binary data + +df_sahhs$total_midpoint6 <- roundmid_any(df_sahhs$total) +df_sahhs$exposed_midpoint6 <- roundmid_any(df_sahhs$exposed) -df$total_midpoint6 <- roundmid_any(df$total) -df$exposed_midpoint6 <- roundmid_any(df$exposed) # Calculate column percentages ------------------------------------------------- +print("Calculate column percentages") + +df_ami$N_midpoint6_derived <- df_ami$total_midpoint6 + +df_ami$percent_midpoint6_derived <- paste0( + ifelse( + df_ami$characteristic == "All", + "", + paste0( + round( + 100 * + (df_ami$total_midpoint6 / + df_ami[df_ami$characteristic == "All", "total_midpoint6"]), + 1 + ), + "%" + ) + ) +) + +df_ami$percent_exposed_midpoint6 <- paste0( + ifelse( + df_ami$characteristic == "All", + "", + paste0( + round( + 100 * + (df_ami$exposed_midpoint6 / + df_ami[df_ami$characteristic == "All", "exposed_midpoint6"]), + 1 + ), + "%" + ) + ) +) + +df_ami <- df_ami[, c( + "characteristic", + "subcharacteristic", + "N_midpoint6_derived", + "percent_midpoint6_derived", + "exposed_midpoint6", + "percent_exposed_midpoint6" +)] + +df_ami[nrow(df_ami) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, "", 0, "") + +df_ami <- dplyr::rename( + df_ami, + "Characteristic" = "characteristic", + "Subcharacteristic" = "subcharacteristic", + "N [midpoint6_derived]" = "N_midpoint6_derived", + "(%) [midpoint6_derived]" = "percent_midpoint6_derived", + "COVID-19 diagnoses [midpoint6]" = "exposed_midpoint6" +) -df$N_midpoint6_derived <- df$total_midpoint6 +df_sahhs$N_midpoint6_derived <- df_sahhs$total_midpoint6 -df$percent_midpoint6_derived <- paste0( +df_sahhs$percent_midpoint6_derived <- paste0( ifelse( - df$characteristic == "All", + df_sahhs$characteristic == "All", "", paste0( round( 100 * - (df$total_midpoint6 / - df[df$characteristic == "All", "total_midpoint6"]), + (df_sahhs$total_midpoint6 / + df_sahhs[df_sahhs$characteristic == "All", "total_midpoint6"]), 1 ), "%" @@ -271,15 +476,15 @@ df$percent_midpoint6_derived <- paste0( ) ) -df$percent_exposed_midpoint6 <- paste0( +df_sahhs$percent_exposed_midpoint6 <- paste0( ifelse( - df$characteristic == "All", + df_sahhs$characteristic == "All", "", paste0( round( 100 * - (df$exposed_midpoint6 / - df[df$characteristic == "All", "exposed_midpoint6"]), + (df_sahhs$exposed_midpoint6 / + df_sahhs[df_sahhs$characteristic == "All", "exposed_midpoint6"]), 1 ), "%" @@ -287,7 +492,7 @@ df$percent_exposed_midpoint6 <- paste0( ) ) -df <- df[, c( +df_sahhs <- df_sahhs[, c( "characteristic", "subcharacteristic", "N_midpoint6_derived", @@ -296,10 +501,10 @@ df <- df[, c( "percent_exposed_midpoint6" )] -df[nrow(df) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, "", 0, "") +df_sahhs[nrow(df_sahhs) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, "", 0, "") -df <- dplyr::rename( - df, +df_sahhs <- dplyr::rename( + df_sahhs, "Characteristic" = "characteristic", "Subcharacteristic" = "subcharacteristic", "N [midpoint6_derived]" = "N_midpoint6_derived", @@ -307,11 +512,18 @@ df <- dplyr::rename( "COVID-19 diagnoses [midpoint6]" = "exposed_midpoint6" ) + # Save Table 1 ----------------------------------------------------------------- print("Save rounded Table 1") write.csv( - df, - paste0(table1_dir, "table1-cohort_", cohort, preex_string, "-midpoint6.csv"), + df_ami, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_ami-midpoint6.csv"), + row.names = FALSE +) + +write.csv( + df_sahhs, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_sahhs-midpoint6.csv"), row.names = FALSE ) diff --git a/analysis/table1/table1_subsample.R b/analysis/table1/table1_subsample.R index ad9d7fc..3f9c36b 100644 --- a/analysis/table1/table1_subsample.R +++ b/analysis/table1/table1_subsample.R @@ -81,52 +81,102 @@ age_bounds <- as.numeric(stringr::str_split(as.vector(age_str), ";")[[1]]) # Load data -------------------------------------------------------------------- print("Load data") -df <- readr::read_rds(paste0( +# subsample, ami +df_ami <- readr::read_rds(paste0( "output/generate_subsample/input_", cohort, - "_clean_subsample.rds" + "_clean_subsample_ami.rds" )) +# subsample, sahhs +df_sahhs <- readr::read_rds(paste0( + "output/generate_subsample/input_", + cohort, + "_clean_subsample_sahhs.rds" +)) -# Check all covariates -------------------------------------------------------- -message("Check all covariates") - -df$cov_bin_ami <- as.logical(df$cov_bin_ami) # outcome -df$cov_bin_sahhs <- as.logical(df$cov_bin_sahhs) # outcome -df$cov_bin_covid <- as.logical(df$cov_bin_covid) # exposure - -df$cov_num_age <- as.numeric(df$cov_num_age) -df$cov_cat_sex <- as.factor(df$cov_cat_sex) -df$cov_cat_ethnicity <- as.factor(df$cov_cat_ethnicity) -df$cov_cat_imd <- as.factor(df$cov_cat_imd) -df$cov_cat_smoking <- as.factor(df$cov_cat_smoking) - -df$cov_bin_carehome <- as.factor(df$cov_bin_carehome) -df$cov_bin_hcworker <- as.factor(df$cov_bin_hcworker) -df$cov_bin_dementia <- as.factor(df$cov_bin_dementia) -df$cov_bin_liver_disease <- as.factor(df$cov_bin_liver_disease) -df$cov_bin_ckd <- as.factor(df$cov_bin_ckd) - -df$cov_bin_cancer <- as.factor(df$cov_bin_cancer) -df$cov_bin_hypertension <- as.factor(df$cov_bin_hypertension) -df$cov_bin_diabetes <- as.factor(df$cov_bin_diabetes) -df$cov_bin_obesity <- as.factor(df$cov_bin_obesity) -df$cov_bin_copd <- as.factor(df$cov_bin_copd) +# nelson aalen estimator +df_ami_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_", cohort, "_ami.rds")) +df_sahhs_nelsonaalen <- readRDS(paste0("output/dataset_clean/nelson_aalen_", cohort, "_sahhs.rds")) -df$cov_bin_depression <- as.factor(df$cov_bin_depression) -df$cov_bin_stroke_all <- as.factor(df$cov_bin_stroke_all) -df$cov_bin_other_ae <- as.factor(df$cov_bin_other_ae) -df$cov_bin_vte <- as.factor(df$cov_bin_vte) -df$cov_bin_hf <- as.factor(df$cov_bin_hf) -df$cov_bin_angina <- as.factor(df$cov_bin_angina) -df$cov_bin_lipidmed <- as.factor(df$cov_bin_lipidmed) -df$cov_bin_antiplatelet <- as.factor(df$cov_bin_antiplatelet) -df$cov_bin_anticoagulant <- as.factor(df$cov_bin_anticoagulant) -df$cov_bin_cocp <- as.factor(df$cov_bin_cocp) +# Check all covariates -------------------------------------------------------- +message("Check all covariates") -df$cov_bin_hrt <- as.factor(df$cov_bin_hrt) -df$strat_cat_region <- as.factor(df$strat_cat_region) +df_ami$cov_bin_ami <- as.logical(df_ami$cov_bin_ami) # outcome +df_ami$cov_bin_sahhs <- as.logical(df_ami$cov_bin_sahhs) # outcome +df_ami$cov_bin_covid <- as.logical(df_ami$cov_bin_covid) # exposure + +df_ami$cov_num_age <- as.numeric(df_ami$cov_num_age) +df_ami$cov_cat_sex <- as.factor(df_ami$cov_cat_sex) +df_ami$cov_cat_ethnicity <- as.factor(df_ami$cov_cat_ethnicity) +df_ami$cov_cat_imd <- as.factor(df_ami$cov_cat_imd) +df_ami$cov_cat_smoking <- as.factor(df_ami$cov_cat_smoking) + +df_ami$cov_bin_carehome <- as.factor(df_ami$cov_bin_carehome) +df_ami$cov_bin_hcworker <- as.factor(df_ami$cov_bin_hcworker) +df_ami$cov_bin_dementia <- as.factor(df_ami$cov_bin_dementia) +df_ami$cov_bin_liver_disease <- as.factor(df_ami$cov_bin_liver_disease) +df_ami$cov_bin_ckd <- as.factor(df_ami$cov_bin_ckd) + +df_ami$cov_bin_cancer <- as.factor(df_ami$cov_bin_cancer) +df_ami$cov_bin_hypertension <- as.factor(df_ami$cov_bin_hypertension) +df_ami$cov_bin_diabetes <- as.factor(df_ami$cov_bin_diabetes) +df_ami$cov_bin_obesity <- as.factor(df_ami$cov_bin_obesity) +df_ami$cov_bin_copd <- as.factor(df_ami$cov_bin_copd) + +df_ami$cov_bin_depression <- as.factor(df_ami$cov_bin_depression) +df_ami$cov_bin_stroke_all <- as.factor(df_ami$cov_bin_stroke_all) +df_ami$cov_bin_other_ae <- as.factor(df_ami$cov_bin_other_ae) +df_ami$cov_bin_vte <- as.factor(df_ami$cov_bin_vte) +df_ami$cov_bin_hf <- as.factor(df_ami$cov_bin_hf) + +df_ami$cov_bin_angina <- as.factor(df_ami$cov_bin_angina) +df_ami$cov_bin_lipidmed <- as.factor(df_ami$cov_bin_lipidmed) +df_ami$cov_bin_antiplatelet <- as.factor(df_ami$cov_bin_antiplatelet) +df_ami$cov_bin_anticoagulant <- as.factor(df_ami$cov_bin_anticoagulant) +df_ami$cov_bin_cocp <- as.factor(df_ami$cov_bin_cocp) + +df_ami$cov_bin_hrt <- as.factor(df_ami$cov_bin_hrt) +df_ami$strat_cat_region <- as.factor(df_ami$strat_cat_region) + + +df_sahhs$cov_bin_sahhs <- as.logical(df_sahhs$cov_bin_sahhs) # outcome +df_sahhs$cov_bin_sahhs <- as.logical(df_sahhs$cov_bin_sahhs) # outcome +df_sahhs$cov_bin_covid <- as.logical(df_sahhs$cov_bin_covid) # exposure + +df_sahhs$cov_num_age <- as.numeric(df_sahhs$cov_num_age) +df_sahhs$cov_cat_sex <- as.factor(df_sahhs$cov_cat_sex) +df_sahhs$cov_cat_ethnicity <- as.factor(df_sahhs$cov_cat_ethnicity) +df_sahhs$cov_cat_imd <- as.factor(df_sahhs$cov_cat_imd) +df_sahhs$cov_cat_smoking <- as.factor(df_sahhs$cov_cat_smoking) + +df_sahhs$cov_bin_carehome <- as.factor(df_sahhs$cov_bin_carehome) +df_sahhs$cov_bin_hcworker <- as.factor(df_sahhs$cov_bin_hcworker) +df_sahhs$cov_bin_dementia <- as.factor(df_sahhs$cov_bin_dementia) +df_sahhs$cov_bin_liver_disease <- as.factor(df_sahhs$cov_bin_liver_disease) +df_sahhs$cov_bin_ckd <- as.factor(df_sahhs$cov_bin_ckd) + +df_sahhs$cov_bin_cancer <- as.factor(df_sahhs$cov_bin_cancer) +df_sahhs$cov_bin_hypertension <- as.factor(df_sahhs$cov_bin_hypertension) +df_sahhs$cov_bin_diabetes <- as.factor(df_sahhs$cov_bin_diabetes) +df_sahhs$cov_bin_obesity <- as.factor(df_sahhs$cov_bin_obesity) +df_sahhs$cov_bin_copd <- as.factor(df_sahhs$cov_bin_copd) + +df_sahhs$cov_bin_depression <- as.factor(df_sahhs$cov_bin_depression) +df_sahhs$cov_bin_stroke_all <- as.factor(df_sahhs$cov_bin_stroke_all) +df_sahhs$cov_bin_other_ae <- as.factor(df_sahhs$cov_bin_other_ae) +df_sahhs$cov_bin_vte <- as.factor(df_sahhs$cov_bin_vte) +df_sahhs$cov_bin_hf <- as.factor(df_sahhs$cov_bin_hf) + +df_sahhs$cov_bin_angina <- as.factor(df_sahhs$cov_bin_angina) +df_sahhs$cov_bin_lipidmed <- as.factor(df_sahhs$cov_bin_lipidmed) +df_sahhs$cov_bin_antiplatelet <- as.factor(df_sahhs$cov_bin_antiplatelet) +df_sahhs$cov_bin_anticoagulant <- as.factor(df_sahhs$cov_bin_anticoagulant) +df_sahhs$cov_bin_cocp <- as.factor(df_sahhs$cov_bin_cocp) + +df_sahhs$cov_bin_hrt <- as.factor(df_sahhs$cov_bin_hrt) +df_sahhs$strat_cat_region <- as.factor(df_sahhs$strat_cat_region) # Table 1 Processing Start ----------------------------------------------------- @@ -135,12 +185,16 @@ print("Table 1 processing") # Remove people with history of COVID-19 --------------------------------------- print("Remove people with history of COVID-19") -df <- df[df$sub_bin_covidhistory == FALSE, ] +df_ami <- df_ami[df_ami$sub_bin_covidhistory == FALSE, ] + +df_sahhs <- df_sahhs[df_sahhs$sub_bin_covidhistory == FALSE, ] # Create exposure indicator ---------------------------------------------------- print("Create exposure indicator") -df$exposed <- df$cov_bin_covid +df_ami$exposed <- df_ami$cov_bin_covid + +df_sahhs$exposed <- df_sahhs$cov_bin_covid # Select for pre-existing conditions print("Select for pre-existing conditions") @@ -149,76 +203,139 @@ print("Select for pre-existing conditions") # because neither asthma nor copdoutcomes are present preex_string <- "" # if (preex != "All") { -# df <- df[df$sup_bin_preex == preex, ] +# df_ami <- df_ami[df_ami$sup_bin_preex == preex, ] +# preex_string <- paste0("-preex_", preex) +# } + +# preex optional argument deliberately ignored, sup_bin_preex does not exist +# because neither asthma nor copdoutcomes are present +preex_string <- "" +# if (preex != "All") { +# df_sahhs <- df_sahhs[df_sahhs$sup_bin_preex == preex, ] # preex_string <- paste0("-preex_", preex) # } + # Define age groups ------------------------------------------------------------ print("Define age groups") -df$cov_cat_age_group <- numerical_to_categorical(df$cov_num_age, age_bounds) # See utility.R +df_ami$cov_cat_age_group <- numerical_to_categorical(df_ami$cov_num_age, age_bounds) # See utility.R + +# df_ami$cov_cat_consrate2019 <- numerical_to_categorical( +# df_ami$cov_num_consrate2019, +# c(1, 6), +# zero_flag = TRUE +# ) + +median_iqr_age <- create_median_iqr_string(df_ami$cov_num_age) # See utility.R + + +df_sahhs$cov_cat_age_group <- numerical_to_categorical(df_sahhs$cov_num_age, age_bounds) # See utility.R -# df$cov_cat_consrate2019 <- numerical_to_categorical( -# df$cov_num_consrate2019, +# df_sahhs$cov_cat_consrate2019 <- numerical_to_categorical( +# df_sahhs$cov_num_consrate2019, # c(1, 6), # zero_flag = TRUE # ) -median_iqr_age <- create_median_iqr_string(df$cov_num_age) # See utility.R +median_iqr_age <- create_median_iqr_string(df_sahhs$cov_num_age) # See utility.R + # Filter data ------------------------------------------------------------------ print("Filter data") -df <- df[, c( +df_ami <- df_ami[, c( + "patient_id", + "exposed", + colnames(df_ami)[grepl("cov_cat_", colnames(df_ami))], + colnames(df_ami)[grepl("strat_cat_", colnames(df_ami))], + colnames(df_ami)[grepl("cov_bin_", colnames(df_ami))] +)] + +df_ami$All <- "All" + +# Filter binary data + +for (colname in colnames(df_ami)[grepl("cov_bin_", colnames(df_ami))]) { + df_ami[[colname]] <- sapply(df_ami[[colname]], as.character) +} + +df_ami <- df_ami %>% + mutate(across(where(is.factor), as.character)) + +df_sahhs <- df_sahhs[, c( "patient_id", "exposed", - colnames(df)[grepl("cov_cat_", colnames(df))], - colnames(df)[grepl("strat_cat_", colnames(df))], - colnames(df)[grepl("cov_bin_", colnames(df))] + colnames(df_sahhs)[grepl("cov_cat_", colnames(df_sahhs))], + colnames(df_sahhs)[grepl("strat_cat_", colnames(df_sahhs))], + colnames(df_sahhs)[grepl("cov_bin_", colnames(df_sahhs))] )] -df$All <- "All" +df_sahhs$All <- "All" # Filter binary data -for (colname in colnames(df)[grepl("cov_bin_", colnames(df))]) { - df[[colname]] <- sapply(df[[colname]], as.character) +for (colname in colnames(df_sahhs)[grepl("cov_bin_", colnames(df_sahhs))]) { + df_sahhs[[colname]] <- sapply(df_sahhs[[colname]], as.character) } -df <- df %>% +df_sahhs <- df_sahhs %>% mutate(across(where(is.factor), as.character)) # Convert to characteristics and subcharacteristics ---------------------------- print("Convert to characteristics and subcharacteristics") -df <- tidyr::pivot_longer( - df, - cols = setdiff(colnames(df), c("patient_id", "exposed")), +df_ami <- tidyr::pivot_longer( + df_ami, + cols = setdiff(colnames(df_ami), c("patient_id", "exposed")), + names_to = "characteristic", + values_to = "subcharacteristic" +) + +df_ami$total <- 1 + +df_sahhs <- tidyr::pivot_longer( + df_sahhs, + cols = setdiff(colnames(df_sahhs), c("patient_id", "exposed")), names_to = "characteristic", values_to = "subcharacteristic" ) -df$total <- 1 +df_sahhs$total <- 1 # Tidy missing data labels ----------------------------------------------------- print("Tidy missing data labels") -df$subcharacteristic <- ifelse( - df$subcharacteristic == "" | - df$subcharacteristic == "unknown" | - is.na(df$subcharacteristic), +df_ami$subcharacteristic <- ifelse( + df_ami$subcharacteristic == "" | + df_ami$subcharacteristic == "unknown" | + is.na(df_ami$subcharacteristic), "Missing", - df$subcharacteristic + df_ami$subcharacteristic +) + +df_sahhs$subcharacteristic <- ifelse( + df_sahhs$subcharacteristic == "" | + df_sahhs$subcharacteristic == "unknown" | + is.na(df_sahhs$subcharacteristic), + "Missing", + df_sahhs$subcharacteristic ) # Aggregate data --------------------------------------------------------------- print("Aggregate data") -df <- aggregate( +df_ami <- aggregate( + cbind(total, exposed) ~ characteristic + subcharacteristic, + data = df_ami, + sum +) + +df_sahhs <- aggregate( cbind(total, exposed) ~ characteristic + subcharacteristic, - data = df, + data = df_sahhs, sum ) @@ -226,45 +343,132 @@ df <- aggregate( # Sort characteristics --------------------------------------------------------- print("Sort characteristics") -df <- df[order(df$characteristic, df$subcharacteristic), ] +df_ami <- df_ami[order(df_ami$characteristic, df_ami$subcharacteristic), ] + +df_sahhs <- df_sahhs[order(df_sahhs$characteristic, df_sahhs$subcharacteristic), ] + # Add in Median IQR print('Add median (IQR) age') # Pastes: "Mean Age (LQ Age - UQ Age)" as a string for each cohort -df[nrow(df) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, 0) +df_sahhs[nrow(df_sahhs) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, 0) + + +# Nelson-Aalen summary statistics --------------------------------------------- +print("Nelson-Aalen summary statistics") + +nelsonaalen_ami_summary <- summary(df_ami_nelsonaalen$H0) +df_ami[nrow(df_ami) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Mean", nelsonaalen_ami_summary[["Mean"]], 0) +df_ami[nrow(df_ami) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Median", nelsonaalen_ami_summary[["Median"]], 0) +df_ami[nrow(df_ami) + 1,] <- c("Nelson-Aalen Estimator (H0)", "IQR", (nelsonaalen_ami_summary[["3rd Qu."]] - nelsonaalen_ami_summary[["1st Qu."]]), 0) + +nelsonaalen_sahhs_summary <- summary(df_sahhs_nelsonaalen$H0) +df_sahhs[nrow(df_sahhs) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Mean", nelsonaalen_sahhs_summary[["Mean"]], 0) +df_sahhs[nrow(df_sahhs) + 1,] <- c("Nelson-Aalen Estimator (H0)", "Median", nelsonaalen_sahhs_summary[["Median"]], 0) +df_sahhs[nrow(df_sahhs) + 1,] <- c("Nelson-Aalen Estimator (H0)", "IQR", (nelsonaalen_sahhs_summary[["3rd Qu."]] - nelsonaalen_sahhs_summary[["1st Qu."]]), 0) + # Save Table 1 ----------------------------------------------------------------- print("Save Table 1") write.csv( - df, - paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_subsample.csv"), + df_ami, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_ami_subsample.csv"), + row.names = FALSE +) + +write.csv( + df_sahhs, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_sahhs_subsample.csv"), row.names = FALSE ) + # Perform redaction ------------------------------------------------------------ print("Perform redaction") -df <- df[df$subcharacteristic != "Median (IQR)", ] # Remove Median IQR row -df <- df[df$subcharacteristic != FALSE, ] # Remove False binary data +df_ami <- df_ami[df_ami$subcharacteristic != "Median (IQR)", ] # Remove Median IQR row +df_ami <- df_ami[df_ami$subcharacteristic != FALSE, ] # Remove False binary data + +df_ami$total_midpoint6 <- roundmid_any(df_ami$total) +df_ami$exposed_midpoint6 <- roundmid_any(df_ami$exposed) + +df_sahhs <- df_sahhs[df_sahhs$subcharacteristic != "Median (IQR)", ] # Remove Median IQR row +df_sahhs <- df_sahhs[df_sahhs$subcharacteristic != FALSE, ] # Remove False binary data + +df_sahhs$total_midpoint6 <- roundmid_any(df_sahhs$total) +df_sahhs$exposed_midpoint6 <- roundmid_any(df_sahhs$exposed) -df$total_midpoint6 <- roundmid_any(df$total) -df$exposed_midpoint6 <- roundmid_any(df$exposed) # Calculate column percentages ------------------------------------------------- +print("Calculate column percentages") + +df_ami$N_midpoint6_derived <- df_ami$total_midpoint6 + +df_ami$percent_midpoint6_derived <- paste0( + ifelse( + df_ami$characteristic == "All", + "", + paste0( + round( + 100 * + (df_ami$total_midpoint6 / + df_ami[df_ami$characteristic == "All", "total_midpoint6"]), + 1 + ), + "%" + ) + ) +) + +df_ami$percent_exposed_midpoint6 <- paste0( + ifelse( + df_ami$characteristic == "All", + "", + paste0( + round( + 100 * + (df_ami$exposed_midpoint6 / + df_ami[df_ami$characteristic == "All", "exposed_midpoint6"]), + 1 + ), + "%" + ) + ) +) + +df_ami <- df_ami[, c( + "characteristic", + "subcharacteristic", + "N_midpoint6_derived", + "percent_midpoint6_derived", + "exposed_midpoint6", + "percent_exposed_midpoint6" +)] + +df_ami[nrow(df_ami) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, "", 0, "") + +df_ami <- dplyr::rename( + df_ami, + "Characteristic" = "characteristic", + "Subcharacteristic" = "subcharacteristic", + "N [midpoint6_derived]" = "N_midpoint6_derived", + "(%) [midpoint6_derived]" = "percent_midpoint6_derived", + "COVID-19 diagnoses [midpoint6]" = "exposed_midpoint6" +) -df$N_midpoint6_derived <- df$total_midpoint6 +df_sahhs$N_midpoint6_derived <- df_sahhs$total_midpoint6 -df$percent_midpoint6_derived <- paste0( +df_sahhs$percent_midpoint6_derived <- paste0( ifelse( - df$characteristic == "All", + df_sahhs$characteristic == "All", "", paste0( round( 100 * - (df$total_midpoint6 / - df[df$characteristic == "All", "total_midpoint6"]), + (df_sahhs$total_midpoint6 / + df_sahhs[df_sahhs$characteristic == "All", "total_midpoint6"]), 1 ), "%" @@ -272,15 +476,15 @@ df$percent_midpoint6_derived <- paste0( ) ) -df$percent_exposed_midpoint6 <- paste0( +df_sahhs$percent_exposed_midpoint6 <- paste0( ifelse( - df$characteristic == "All", + df_sahhs$characteristic == "All", "", paste0( round( 100 * - (df$exposed_midpoint6 / - df[df$characteristic == "All", "exposed_midpoint6"]), + (df_sahhs$exposed_midpoint6 / + df_sahhs[df_sahhs$characteristic == "All", "exposed_midpoint6"]), 1 ), "%" @@ -288,7 +492,7 @@ df$percent_exposed_midpoint6 <- paste0( ) ) -df <- df[, c( +df_sahhs <- df_sahhs[, c( "characteristic", "subcharacteristic", "N_midpoint6_derived", @@ -297,10 +501,10 @@ df <- df[, c( "percent_exposed_midpoint6" )] -df[nrow(df) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, "", 0, "") +df_sahhs[nrow(df_sahhs) + 1, ] <- c("Age, years", "Median (IQR)", median_iqr_age, "", 0, "") -df <- dplyr::rename( - df, +df_sahhs <- dplyr::rename( + df_sahhs, "Characteristic" = "characteristic", "Subcharacteristic" = "subcharacteristic", "N [midpoint6_derived]" = "N_midpoint6_derived", @@ -308,11 +512,18 @@ df <- dplyr::rename( "COVID-19 diagnoses [midpoint6]" = "exposed_midpoint6" ) + # Save Table 1 ----------------------------------------------------------------- print("Save rounded Table 1") write.csv( - df, - paste0(table1_dir, "table1-cohort_", cohort, preex_string, "-midpoint6_subsample.csv"), + df_ami, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_ami-midpoint6_subsample.csv"), row.names = FALSE ) + +write.csv( + df_sahhs, + paste0(table1_dir, "table1-cohort_", cohort, preex_string, "_sahhs-midpoint6_subsample.csv"), + row.names = FALSE +) \ No newline at end of file diff --git a/analysis/table_missingness/table_missingness.R b/analysis/table_missingness/table_missingness.R new file mode 100644 index 0000000..1783bed --- /dev/null +++ b/analysis/table_missingness/table_missingness.R @@ -0,0 +1,116 @@ +# ------------------------------------------------------------------------------ +# +# table_missingness.R +# +# This file generates the missingness table +# +# Arguments: +# - cohort - string, defines which of three opensafely cohorts to describe +# (prevax, vax, unvax) +# +# Returns: +# - table! +# +# Authors: Emma Tarmey, Venexia Walker, UoB ehrQL Team +# +# ------------------------------------------------------------------------------ + + +# Refresh local R session ------------------------------------------------------ +print("Refresh local R session") + +rm(list=ls()) + + +# Load libraries --------------------------------------------------------------- +print("Load libraries") + +library(magrittr) +library(here) +library(dplyr) +library(fs) + + +# Define table_missingness output folder --------------------------------------------------------- +print("Creating output/table_missingness output folder") + +table_missingness_dir <- "output/table_missingness/" +dir_create(here::here(table_missingness_dir)) + + +# Source common functions ------------------------------------------------------ +print("Source common functions") +source("analysis/utility.R") + + +# Specify arguments ------------------------------------------------------------ +print("Specify arguments") + +args <- commandArgs(trailingOnly = TRUE) + +print(length(args)) + +if (length(args) == 0) { + # default argument values + cohort <- "prevax" +} else { + # YAML arguments + cohort <- args[[1]] +} + + +# Load data -------------------------------------------------------------------- +print("Load data") + +df <- readr::read_rds(paste0( + "output/dataset_clean/input_", + cohort, + "_clean_prehoc.rds" +)) + +# df <- readr::read_rds(paste0( +# "output/dataset_definition/input_", +# cohort, +# ".rds" +# )) + +df <- as.data.frame(df) + + +# Convert missingness categories to NA ----------------------------------------- + +df$cov_cat_ethnicity[df$cov_cat_ethnicity == "Missing"] <- NA +df$cov_cat_smoking[df$cov_cat_smoking == "Missing"] <- NA + + +# Make table ------------------------------------------------------------------- + +all_var_names <- c( + "cov_num_age", "cov_cat_sex", "cov_num_bmi", "cov_cat_ethnicity", "cov_cat_imd", + "cov_cat_smoking", "cov_bin_carehome", "cov_bin_hcworker", "cov_bin_dementia", + "cov_bin_liver_disease", "cov_bin_ckd", "cov_bin_cancer", "cov_bin_hypertension", + "cov_bin_diabetes", "cov_bin_obesity", "cov_bin_copd", "cov_bin_depression", "cov_bin_stroke_all", + "cov_bin_other_ae", "cov_bin_vte", "cov_bin_hf", "cov_bin_angina", "cov_bin_lipidmed", + "cov_bin_antiplatelet", "cov_bin_anticoagulant", "cov_bin_cocp", "cov_bin_hrt", "strat_cat_region" +) + +all_var_missingness <- rep(0, length.out = length(all_var_names)) +names(all_var_missingness) <- all_var_names + +all_var_missingness["cov_num_bmi"] <- signif(100 * (sum(is.na(df$cov_num_bmi)) / length(df$cov_num_bmi)), digits = 4) +all_var_missingness["cov_cat_ethnicity"] <- signif(100 * (sum(is.na(df$cov_cat_ethnicity)) / length(df$cov_cat_ethnicity)), digits = 4) +all_var_missingness["cov_cat_smoking"] <- signif(100 * (sum(is.na(df$cov_cat_smoking)) / length(df$cov_cat_smoking)), digits = 4) + +missingness_table <- data.frame( + Covariates = all_var_names, + Missingness = all_var_missingness +) + + +# Save results ----------------------------------------------------------------- + +write.csv( + missingness_table, + paste0(table_missingness_dir, "table_missingness_cohort_", cohort, ".csv"), + row.names = TRUE +) diff --git a/analysis/unconfoundedness_test/unconfoundedness_test.R b/analysis/unconfoundedness_test/unconfoundedness_test.R index 9b2ccb4..95ee669 100644 --- a/analysis/unconfoundedness_test/unconfoundedness_test.R +++ b/analysis/unconfoundedness_test/unconfoundedness_test.R @@ -84,12 +84,21 @@ if (grepl("ami", name)) { # Load data -------------------------------------------------------------------- print("Load data") -# subsample -df <- readr::read_rds(paste0( - "output/generate_subsample/input_", - cohort, - "_clean_subsample.rds" -)) +if (grepl("ami", name)) { + # subsample + df <- readr::read_rds(paste0( + "output/generate_subsample/input_", + cohort, + "_clean_subsample_ami.rds" + )) +} else { + # subsample + df <- readr::read_rds(paste0( + "output/generate_subsample/input_", + cohort, + "_clean_subsample_sahhs.rds" + )) +} # subsample model_input_df <- readr::read_rds(paste0( @@ -113,6 +122,7 @@ logistic_df <- (df %>% select(c( cov_num_age, cov_cat_sex, + cov_num_bmi, cov_cat_ethnicity, cov_cat_imd, cov_cat_smoking, @@ -145,6 +155,9 @@ logistic_df <- (df %>% select(c( strat_cat_region ))) +# prevent conversion to factor +logistic_df$cov_num_bmi <- as.numeric(logistic_df$cov_num_bmi) + # Data preparation for cox model ---------------------------------- message("Data preparation for cox model") @@ -159,6 +172,7 @@ cox_df <- (model_input_df %>% select(c( cov_num_age, cov_cat_sex, + cov_num_bmi, cov_cat_ethnicity, cov_cat_imd, cov_cat_smoking, @@ -191,6 +205,9 @@ cox_df <- (model_input_df %>% select(c( strat_cat_region ))) +# prevent conversion to factor +cox_df$cov_num_bmi <- as.numeric(cox_df$cov_num_bmi) + cox_df$outcome_cox_dates <- rep(as.Date(NA), times = nrow(cox_df)) cox_df$cens_status <- rep(NA, times = nrow(cox_df)) @@ -234,11 +251,13 @@ fully_adjusted_outcome_regression_formula <- make_outcome_formula( fully_adjusted_exposure_regression <- glm( fully_adjusted_exposure_regression_formula, family = "binomial", - data = logistic_df + data = logistic_df, + weights = generate_weights(sample_size = nrow(logistic_df)) ) fully_adjusted_outcome_regression <- coxph( formula = as.formula(fully_adjusted_outcome_regression_formula), - data = cox_df + data = cox_df, + weights = generate_weights(sample_size = nrow(cox_df)) ) # extract regression results @@ -325,11 +344,13 @@ lasso_outcome_regression_formula <- make_outcome_formula( lasso_exposure_regression <- glm( lasso_exposure_regression_formula, family = "binomial", - data = logistic_df + data = logistic_df, + weights = generate_weights(sample_size = nrow(logistic_df)) ) lasso_outcome_regression <- coxph( formula = as.formula(lasso_outcome_regression_formula), - data = cox_df + data = cox_df, + weights = generate_weights(sample_size = nrow(cox_df)) ) # extract regression results @@ -416,11 +437,13 @@ lasso_X_outcome_regression_formula <- make_outcome_formula( lasso_X_exposure_regression <- glm( lasso_X_exposure_regression_formula, family = "binomial", - data = logistic_df + data = logistic_df, + weights = generate_weights(sample_size = nrow(logistic_df)) ) lasso_X_outcome_regression <- coxph( formula = as.formula(lasso_X_outcome_regression_formula), - data = cox_df + data = cox_df, + weights = generate_weights(sample_size = nrow(cox_df)) ) # extract regression results @@ -507,11 +530,13 @@ lasso_union_outcome_regression_formula <- make_outcome_formula( lasso_union_exposure_regression <- glm( lasso_union_exposure_regression_formula, family = "binomial", - data = logistic_df + data = logistic_df, + weights = generate_weights(sample_size = nrow(logistic_df)) ) lasso_union_outcome_regression <- coxph( formula = as.formula(lasso_union_outcome_regression_formula), - data = cox_df + data = cox_df, + weights = generate_weights(sample_size = nrow(cox_df)) ) # extract regression results diff --git a/analysis/utility.R b/analysis/utility.R index 48a5d64..ccaf8bf 100644 --- a/analysis/utility.R +++ b/analysis/utility.R @@ -138,22 +138,49 @@ make_outcome_formula <- function(vars_selected = NULL, outcome = NULL) { } -convert_terms_to_vars <- function(terms = NULL, all_var_names = NULL) { +convert_terms_to_vars <- function(terms = NULL) { + all_var_names <- c( + "cov_bin_covid", + "cov_num_age", "cov_cat_sex", "cov_num_bmi", "cov_cat_ethnicity", "cov_cat_imd", + "cov_cat_smoking", "cov_bin_carehome", "cov_bin_hcworker", "cov_bin_dementia", + "cov_bin_liver_disease", "cov_bin_ckd", "cov_bin_cancer", "cov_bin_hypertension", + "cov_bin_diabetes", "cov_bin_obesity", "cov_bin_copd", "cov_bin_depression", "cov_bin_stroke_all", + "cov_bin_other_ae", "cov_bin_vte", "cov_bin_hf", "cov_bin_angina", "cov_bin_lipidmed", + "cov_bin_antiplatelet", "cov_bin_anticoagulant", "cov_bin_cocp", "cov_bin_hrt", "strat_cat_region" + ) + vars <- c() for (term in terms) { - # split by capital letter, extract first term - # covariate names always all lower case - # Factor Level Names always begin with upper case - new_term <- sapply(strsplit(x = term, split = '([[:upper:]])'), `[`, 1) - vars <- c(vars, new_term) + for (name in all_var_names) { + if (stringr::str_detect(term, name)) { + vars <- append(vars, name) + } + } } # remove duplicates (i.e. two levels are significant) vars <- unique(vars) - # remove any trailing . (leftover from levels) - vars <- gsub('.', '', vars) - return (vars) } + + +generate_weights <- function(initial_weights = NULL, sample_size = NULL, num_imps = NULL) { + if (is.null(num_imps)) { + num_imps <- get_number_of_imputed_datasets() + } + + if (is.null(initial_weights)) { + weights <- rep((1 / num_imps), length.out = sample_size) + } else { + weights <- ((1 / num_imps) * initial_weights) + } + + return (weights) +} + + +get_number_of_imputed_datasets <- function() { + return (10) +} diff --git a/project.yaml b/project.yaml index a7e2fc9..e3d8dad 100644 --- a/project.yaml +++ b/project.yaml @@ -56,25 +56,40 @@ actions: venn: output/dataset_clean/venn-cohort_prevax.rds cohort_clean: output/dataset_clean/input_prevax_clean_prehoc.rds - ## post_hoc_vars_cohort_prevax + ## table_missingness_cohort_prevax - post_hoc_vars_cohort_prevax: - run: r:latest analysis/post_hoc_vars/post_hoc_vars.R prevax + table_missingness_cohort_prevax: + run: r:latest analysis/table_missingness/table_missingness.R prevax + needs: + - generate_input_prevax_clean + - generate_input_prevax + outputs: + moderately_sensitive: + table_missingness: output/table_missingness/table_missingness_cohort_prevax.csv + + ## apply_across_MI_cohort_prevax + + apply_across_MI_cohort_prevax: + run: r:latest analysis/apply_across_MI/apply_across_MI.R prevax needs: - generate_input_prevax_clean outputs: highly_sensitive: - cohort_clean: output/dataset_clean/input_prevax_clean.rds + cohort_clean_ami: output/dataset_clean/input_prevax_clean_ami.rds + cohort_clean_sahhs: output/dataset_clean/input_prevax_clean_sahhs.rds + nelsonaalen_ami: output/dataset_clean/nelson_aalen_prevax_ami.rds + nelsonaalen_sahhs: output/dataset_clean/nelson_aalen_prevax_sahhs.rds ## generate_subsample_cohort_prevax generate_subsample_cohort_prevax: run: r:latest analysis/generate_subsample/generate_subsample.R prevax needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: highly_sensitive: - cohort_clean_subsample: output/generate_subsample/input_prevax_clean_subsample.rds + cohort_clean_subsample_ami: output/generate_subsample/input_prevax_clean_subsample_ami.rds + cohort_clean_subsample_sahhs: output/generate_subsample/input_prevax_clean_subsample_sahhs.rds ## Generate cox model input data for 10% subsample study population ## make_model_input_subsample-cohort_prevax-main-ami @@ -146,32 +161,39 @@ actions: table1-cohort_prevax: run: 'r:v2 analysis/table1/table1.R prevax 18;30;40;50;60;70;80;90 ' needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: moderately_sensitive: - table1: output/table1/table1-cohort_prevax.csv - table1_midpoint6: output/table1/table1-cohort_prevax-midpoint6.csv + table1_ami: output/table1/table1-cohort_prevax_ami.csv + table1_ami_midpoint6: output/table1/table1-cohort_prevax_ami-midpoint6.csv + table1_sahhs: output/table1/table1-cohort_prevax_sahhs.csv + table1_sahhs_midpoint6: output/table1/table1-cohort_prevax_sahhs-midpoint6.csv ## Generate make-table1-output make-table1-output: - run: r:v2 analysis/make_output/make_other_output.R table1 prevax + run: r:v2 analysis/make_output/make_table1_output.R table1 prevax needs: + - apply_across_MI_cohort_prevax - table1-cohort_prevax outputs: moderately_sensitive: - other_output_midpoint6: output/make_output/table1_output_midpoint6.csv + other_output_ami_midpoint6: output/make_output/table1_output_ami_midpoint6.csv + other_output_sahhs_midpoint6: output/make_output/table1_output_sahhs_midpoint6.csv ## Generate table1_cohort_prevax_subsample table1-cohort_prevax_subsample: run: 'r:v2 analysis/table1/table1_subsample.R prevax 18;30;40;50;60;70;80;90 ' needs: + - apply_across_MI_cohort_prevax - generate_subsample_cohort_prevax outputs: moderately_sensitive: - table1_subsample: output/table1/table1-cohort_prevax_subsample.csv - table1_midpoint6_subsample: output/table1/table1-cohort_prevax-midpoint6_subsample.csv + table1_ami_subsample: output/table1/table1-cohort_prevax_ami_subsample.csv + table1_ami_midpoint6_subsample: output/table1/table1-cohort_prevax_ami-midpoint6_subsample.csv + table1_sahhs_subsample: output/table1/table1-cohort_prevax_sahhs_subsample.csv + table1_sahhs_midpoint6_subsample: output/table1/table1-cohort_prevax_sahhs-midpoint6_subsample.csv ## Generate lasso_var_selection-cohort_prevax-main-ami @@ -530,14 +552,14 @@ actions: make_model_input-cohort_prevax-main-ami: run: r:latest analysis/model/make_model_input.R cohort_prevax-main-ami needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: highly_sensitive: model_input: output/model/model_input-cohort_prevax-main-ami.rds cox_ipw-cohort_prevax-main-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-ami.rds --ipw=TRUE --exposure=exp_date @@ -568,14 +590,14 @@ actions: make_model_input-cohort_prevax-main-stroke_sahhs: run: r:latest analysis/model/make_model_input.R cohort_prevax-main-stroke_sahhs needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: highly_sensitive: model_input: output/model/model_input-cohort_prevax-main-stroke_sahhs.rds cox_ipw-cohort_prevax-main-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -606,14 +628,14 @@ actions: make_model_input-cohort_prevax-sub_covidhospital_FALSE-ami: run: r:latest analysis/model/make_model_input.R cohort_prevax-sub_covidhospital_FALSE-ami needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: highly_sensitive: model_input: output/model/model_input-cohort_prevax-sub_covidhospital_FALSE-ami.rds cox_ipw-cohort_prevax-sub_covidhospital_FALSE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-ami.rds --ipw=TRUE --exposure=exp_date @@ -644,14 +666,14 @@ actions: make_model_input-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs: run: r:latest analysis/model/make_model_input.R cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: highly_sensitive: model_input: output/model/model_input-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs.rds cox_ipw-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -682,14 +704,14 @@ actions: make_model_input-cohort_prevax-sub_covidhospital_TRUE-ami: run: r:latest analysis/model/make_model_input.R cohort_prevax-sub_covidhospital_TRUE-ami needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: highly_sensitive: model_input: output/model/model_input-cohort_prevax-sub_covidhospital_TRUE-ami.rds cox_ipw-cohort_prevax-sub_covidhospital_TRUE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-ami.rds --ipw=TRUE --exposure=exp_date @@ -720,14 +742,14 @@ actions: make_model_input-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs: run: r:latest analysis/model/make_model_input.R cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs needs: - - post_hoc_vars_cohort_prevax + - apply_across_MI_cohort_prevax outputs: highly_sensitive: model_input: output/model/model_input-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs.rds cox_ipw-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -851,7 +873,7 @@ actions: lasso_cox_ipw-cohort_prevax-main-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-ami.rds --ipw=TRUE --exposure=exp_date @@ -882,7 +904,7 @@ actions: lasso_cox_ipw-cohort_prevax-main-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -913,7 +935,7 @@ actions: lasso_cox_ipw-cohort_prevax-sub_covidhospital_FALSE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-ami.rds --ipw=TRUE --exposure=exp_date @@ -944,7 +966,7 @@ actions: lasso_cox_ipw-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -975,7 +997,7 @@ actions: lasso_cox_ipw-cohort_prevax-sub_covidhospital_TRUE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-ami.rds --ipw=TRUE --exposure=exp_date @@ -1006,7 +1028,7 @@ actions: lasso_cox_ipw-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -1039,7 +1061,7 @@ actions: lasso_X_cox_ipw-cohort_prevax-main-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-ami.rds --ipw=TRUE --exposure=exp_date @@ -1070,7 +1092,7 @@ actions: lasso_X_cox_ipw-cohort_prevax-main-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -1101,7 +1123,7 @@ actions: lasso_X_cox_ipw-cohort_prevax-sub_covidhospital_FALSE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-ami.rds --ipw=TRUE --exposure=exp_date @@ -1132,7 +1154,7 @@ actions: lasso_X_cox_ipw-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -1163,7 +1185,7 @@ actions: lasso_X_cox_ipw-cohort_prevax-sub_covidhospital_TRUE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-ami.rds --ipw=TRUE --exposure=exp_date @@ -1194,7 +1216,7 @@ actions: lasso_X_cox_ipw-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -1227,7 +1249,7 @@ actions: lasso_union_cox_ipw-cohort_prevax-main-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-ami.rds --ipw=TRUE --exposure=exp_date @@ -1258,7 +1280,7 @@ actions: lasso_union_cox_ipw-cohort_prevax-main-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-main-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -1289,7 +1311,7 @@ actions: lasso_union_cox_ipw-cohort_prevax-sub_covidhospital_FALSE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-ami.rds --ipw=TRUE --exposure=exp_date @@ -1320,7 +1342,7 @@ actions: lasso_union_cox_ipw-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_FALSE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date @@ -1351,7 +1373,7 @@ actions: lasso_union_cox_ipw-cohort_prevax-sub_covidhospital_TRUE-ami: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-ami.rds --ipw=TRUE --exposure=exp_date @@ -1382,7 +1404,7 @@ actions: lasso_union_cox_ipw-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs: run: |- - cox-ipw:v0.0.39 + r:v2 analysis/cox_ipw/cox-ipw.R --df_input=model/model_input-cohort_prevax-sub_covidhospital_TRUE-stroke_sahhs.rds --ipw=TRUE --exposure=exp_date