diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile new file mode 100644 index 0000000..506a40d --- /dev/null +++ b/.devcontainer/Dockerfile @@ -0,0 +1,31 @@ +FROM remlapmot/r-docker:2024-04-02-rstudio + +LABEL org.opencontainers.image.source https://github.com/opensafely/research-template + +# we are going to use an apt cache on the host, so disable the default debian +# docker clean up that deletes that cache on every apt install +RUN rm -f /etc/apt/apt.conf.d/docker-clean + +# Install python 3.10. This is the version used by the python-docker +# image, used for analyses using the OpenSAFELY pipeline. +RUN --mount=type=cache,target=/var/cache/apt \ + echo "deb http://ppa.launchpad.net/deadsnakes/ppa/ubuntu focal main" > /etc/apt/sources.list.d/deadsnakes-ppa.list &&\ + /usr/lib/apt/apt-helper download-file 'https://keyserver.ubuntu.com/pks/lookup?op=get&search=0xf23c5a6cf475977595c89f51ba6932366a755776' /etc/apt/trusted.gpg.d/deadsnakes.asc &&\ + apt update &&\ + apt install -y curl python3.10 python3.10-distutils python3.10-venv &&\ + # Pip for Python 3.10 isn't included in deadsnakes, so install separately + curl https://bootstrap.pypa.io/get-pip.py | python3.10 &&\ + # Set default python, so that the Python virtualenv works as expected + rm /usr/bin/python3 && ln -s /usr/bin/python3.10 /usr/bin/python3 + +# Copy the Python virtualenv from OpenSAFELY Python action image +COPY --from=ghcr.io/opensafely-core/python:v2 /opt/venv /opt/venv + +# Create a local user and give it sudo (aka root) permissions +RUN usermod -aG sudo rstudio +RUN echo '%sudo ALL=(ALL) NOPASSWD:ALL' >> /etc/sudoers + +# Required for installing opensafely cli +ENV PATH="/home/rstudio/.local/bin:${PATH}" + +USER rstudio diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 0000000..8baae2b --- /dev/null +++ b/.devcontainer/devcontainer.json @@ -0,0 +1,43 @@ +// For format details, see https://aka.ms/devcontainer.json. For config options, see the +// README at: https://github.com/devcontainers/templates/tree/main/src/python +{ + "name": "OpenSAFELY", + "image": "ghcr.io/opensafely/research-template:latest", + // Features to add to the dev container. More info: https://containers.dev/features. + "features": { + "ghcr.io/devcontainers/features/docker-in-docker:2": { + "moby": true, + "azureDnsAutoDetection": true, + "installDockerBuildx": true, + "version": "latest", + "dockerDashComposeVersion": "v2" + } + }, + "postCreateCommand": "/bin/bash .devcontainer/postCreate.sh", + "postAttachCommand": { + "rstudio-start": "sudo rstudio-server start" + }, + "forwardPorts": [ + 8787 + ], + "portsAttributes": { + "8787": { + "label": "RStudio IDE" + } + }, + // Configure tool-specific properties. + "customizations": { + "vscode": { + "extensions": [ + "ms-python.python", + "ms-toolsai.jupyter", + "ms-toolsai.jupyter-renderers" + ] + } + }, + // Uncomment to connect as root instead. More info: https://aka.ms/dev-containers-non-root. + // "remoteUser": "root" + "remoteEnv": { + "MAX_WORKERS": "2" + } +} diff --git a/.devcontainer/postCreate.sh b/.devcontainer/postCreate.sh new file mode 100644 index 0000000..cac58e6 --- /dev/null +++ b/.devcontainer/postCreate.sh @@ -0,0 +1,15 @@ +#!/bin/bash + +set -euo pipefail + +pip3 install --user -r .devcontainer/requirements.in + +#set R working directory +! grep -q `pwd` $R_HOME/etc/Rprofile.site && sudo tee -a $R_HOME/etc/Rprofile.site <<< "setwd(\"`pwd`\")" +#set RStudio working directory +! grep -q `pwd` ~/.config/rstudio/rstudio-prefs.json && cat ~/.config/rstudio/rstudio-prefs.json | jq ". + {\"initial_working_directory\":\"`pwd`\"}" > ~/.config/rstudio/rstudio-prefs.json + +#download and extract latest ehrql source +wget https://github.com/opensafely-core/ehrql/archive/main.zip -P .devcontainer +unzip -o .devcontainer/main.zip -d .devcontainer/ +rm .devcontainer/main.zip diff --git a/.devcontainer/requirements.in b/.devcontainer/requirements.in new file mode 100644 index 0000000..8136f40 --- /dev/null +++ b/.devcontainer/requirements.in @@ -0,0 +1 @@ +opensafely diff --git a/.gitignore b/.gitignore index d8eaede..d5ce6cb 100644 --- a/.gitignore +++ b/.gitignore @@ -9,10 +9,12 @@ venv/ .DS_Store .Rproj.user .Rhistory +.devcontainer/ehrql-main reports/variation/* reports/coverage/table_prop_eligible_clinc_demo.csv reports/coverage/mabs_and_antivirals_coverage_report.html reports/coverage/figures/* reports/coverage/tables/* -released_outputs/reports_wip/coverage/* -reports/coverage_wip/* \ No newline at end of file +released_outputs/* +reports/coverage_wip/* +analysis/descriptive/mabs-and-avs-by-stp.html \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json index 4a01863..22c1591 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,7 +1,7 @@ { - "python.linting.pylintEnabled": false, - "python.linting.flake8Enabled": true, - "python.linting.enabled": true, + "python.analysis.extraPaths": [".devcontainer/ehrql-main/"], + "python.defaultInterpreterPath": "/opt/venv/bin/python3.10", + "python.terminal.activateEnvironment": true, "data.preview.create.json.schema": false, "files.associations": { "*.feather": "arrow", diff --git a/analysis/descriptive/coverage_report_data.R b/analysis/descriptive/coverage_report_data.R index af0ef3d..a2b1a23 100644 --- a/analysis/descriptive/coverage_report_data.R +++ b/analysis/descriptive/coverage_report_data.R @@ -751,3 +751,44 @@ groups2 <- groups2 %>% write_csv(rbind(all, groups) %>% filter(!is.na(tb)), fs::path(output_dir, "table_time_to_treat_redacted.csv")) write_csv(groups2, fs::path(output_dir, "table_time_to_treat_groups_redacted.csv")) + + +# Treatment recording over time +plot_data_treatment_codes <- data_processed_clean %>% + filter(!is.na(last_treatment_date)) %>% + select(last_treatment_date, last_treatment_type) %>% + rbind(data_processed_clean %>% + filter(!is.na(last_treatment_date)) %>% + select(last_treatment_date, last_treatment_type) %>% + mutate(last_treatment_type = "All")) %>% + group_by(last_treatment_date, last_treatment_type) %>% + tally() %>% + group_by(last_treatment_type) %>% + arrange(last_treatment_type, last_treatment_date) %>% + complete(last_treatment_date = seq.Date(min(last_treatment_date, na.rm = T), max(last_treatment_date, na.rm = T), by="day")) %>% + mutate(count = ifelse(is.na(n), 0, n), + count_redacted = plyr::round_any(count, 10), + count_redacted = ifelse(count < threshold, NA, count_redacted)) %>% + # cum_count = cumsum(count), + # cum_count_redacted = plyr::round_any(cum_count, 10), + # cum_count_redacted = ifelse(cum_count < threshold, NA, cum_count_redacted) + select(-n) %>% + arrange(last_treatment_type, last_treatment_date) + +plot_order <- plot_data_treatment_codes %>% + group_by(last_treatment_type) %>% + mutate(order = max(count_redacted, na.rm = T)) %>% + arrange(desc(order)) %>% + filter(count_redacted == order) %>% + select(last_treatment_type, order) %>% + distinct() + +treatment_codes_plot_data <- plot_data_treatment_codes %>% + mutate(last_treatment_type = factor(last_treatment_type, levels = plot_order$last_treatment_type)) + +write_csv(treatment_codes_plot_data %>% + select(last_treatment_date, count_redacted, last_treatment_type), + fs::path(output_dir, "table_last_treatment_codes_redacted.csv")) +write_csv(treatment_plot_data, fs::path(output_dir2, "table_last_treatment_codes.csv")) + +print("treatment_codes_plot_data saved") diff --git a/analysis/descriptive/crude_outcomes.Rmd b/analysis/descriptive/crude_outcomes.Rmd index eedb272..7d7786e 100644 --- a/analysis/descriptive/crude_outcomes.Rmd +++ b/analysis/descriptive/crude_outcomes.Rmd @@ -685,7 +685,7 @@ sgtf_outcomes <- data_processed_clean %>% # Redact values < threshold count = ifelse(count < threshold, NA, as.numeric(count)), covid_positive_test = ifelse(covid_positive_test < threshold, NA, covid_positive_test), - covid_hospital_admission = ifelse(covid_hospital_admission < threshold, NA, as.numeric(covid_hospital_admission)), + covid_hospital_admission = ifelse(covid_hospital_admission < threshold, NA_integer_, as.numeric(covid_hospital_admission)), covid_hospitalisation_critical_care = ifelse(covid_hospitalisation_critical_care < threshold, NA, as.numeric(covid_hospitalisation_critical_care)), covid_death = ifelse(covid_death < threshold, NA, as.numeric(covid_death)), any_death = ifelse(any_death < threshold, NA, as.numeric(any_death)), diff --git a/analysis/descriptive/mabs-and-avs-by-stp.Rmd b/analysis/descriptive/mabs-and-avs-by-stp.Rmd index 8b63f24..c1c2890 100644 --- a/analysis/descriptive/mabs-and-avs-by-stp.Rmd +++ b/analysis/descriptive/mabs-and-avs-by-stp.Rmd @@ -498,8 +498,8 @@ table_elig_treat_redacted <- table_elig_treat %>% Paxlovid = plyr::round_any(as.numeric(Paxlovid), 10), Sotrovimab = plyr::round_any(as.numeric(Sotrovimab), 10), Remdesivir = plyr::round_any(as.numeric(Remdesivir), 10), - Molnupiravir = plyr::round_any(Molnupiravir, 10), - Casirivimab = plyr::round_any(Casirivimab, 10)) %>% + Molnupiravir = plyr::round_any(as.numeric(Molnupiravir), 10), + Casirivimab = plyr::round_any(as.numeric(Casirivimab), 10)) %>% arrange(desc(decile)) %>% mutate(Paxlovid_perc = paste(round(Paxlovid/Treated*100, digits = 0), " (", round((Paxlovid/Treated - 1.96*sqrt((Paxlovid/Treated)*(1-Paxlovid/Treated)/Treated))*100, digits = 0), diff --git a/analysis/process/process_data.R b/analysis/process/process_data.R index ebef981..53c08d4 100644 --- a/analysis/process/process_data.R +++ b/analysis/process/process_data.R @@ -59,6 +59,12 @@ data_extract0 <- read_csv( remdesivir_covid_therapeutics = col_date(format = "%Y-%m-%d"), molnupiravir_covid_therapeutics = col_date(format = "%Y-%m-%d"), casirivimab_covid_therapeutics = col_date(format = "%Y-%m-%d"), + + paxlovid_covid_therapeutics_last = col_date(format = "%Y-%m-%d"), + sotrovimab_covid_therapeutics_last = col_date(format = "%Y-%m-%d"), + remdesivir_covid_therapeutics_last = col_date(format = "%Y-%m-%d"), + molnupiravir_covid_therapeutics_last = col_date(format = "%Y-%m-%d"), + casirivimab_covid_therapeutics_last = col_date(format = "%Y-%m-%d"), # ELIGIBILITY CRITERIA VARIABLES ---- covid_test_positive = col_logical(), @@ -120,7 +126,7 @@ data_extract0 <- read_csv( serious_mental_illness_nhsd = col_logical(), sickle_cell_disease_nhsd = col_date(format = "%Y-%m-%d"), vaccination_status = col_character(), - first_lc_code_date = col_date(format = "%Y-%m-%d"), + first_long_covid = col_date(format = "%Y-%m-%d"), # COVID VARIENT sgtf = col_character(), @@ -151,7 +157,13 @@ if(Sys.getenv("OPENSAFELY_BACKEND") %in% c("", "expectations")){ remdesivir_covid_therapeutics = as.Date(ifelse(!is.na(remdesivir_covid_therapeutics), date, NA), origin = "1970-01-01"), molnupiravir_covid_therapeutics = as.Date(ifelse(!is.na(molnupiravir_covid_therapeutics), date, NA), origin = "1970-01-01"), casirivimab_covid_therapeutics = as.Date(ifelse(!is.na(casirivimab_covid_therapeutics), date, NA), origin = "1970-01-01"), - + + paxlovid_covid_therapeutics_last = as.Date(ifelse(!is.na(paxlovid_covid_therapeutics_last), date, NA), origin = "1970-01-01"), + sotrovimab_covid_therapeutics_last = as.Date(ifelse(!is.na(sotrovimab_covid_therapeutics_last), date, NA), origin = "1970-01-01"), + remdesivir_covid_therapeutics_last = as.Date(ifelse(!is.na(remdesivir_covid_therapeutics_last), date, NA), origin = "1970-01-01"), + molnupiravir_covid_therapeutics_last = as.Date(ifelse(!is.na(molnupiravir_covid_therapeutics_last), date, NA), origin = "1970-01-01"), + casirivimab_covid_therapeutics_last = as.Date(ifelse(!is.na(casirivimab_covid_therapeutics_last), date, NA), origin = "1970-01-01"), + covid_positive_test_30_days_post_elig_or_treat = as.Date(ifelse(covid_positive_test_30_days_post_elig_or_treat > covid_test_positive_date + 30, covid_positive_test_30_days_post_elig_or_treat, NA), origin = "1970-01-01")) @@ -192,6 +204,19 @@ data_processed <- data_extract %>% treatment_date == casirivimab_covid_therapeutics ~ "Casirivimab", TRUE ~ NA_character_), + # last MAB or Antiviral recorded + last_treatment_date = as.Date(pmin(paxlovid_covid_therapeutics_last, sotrovimab_covid_therapeutics_last, + remdesivir_covid_therapeutics_last, molnupiravir_covid_therapeutics_last, + casirivimab_covid_therapeutics_last, na.rm = TRUE), origin = "1970-01-01"), + last_treatment_type = case_when( + last_treatment_date == paxlovid_covid_therapeutics_last ~ "Paxlovid", + last_treatment_date == sotrovimab_covid_therapeutics_last ~ "Sotrovimab", + last_treatment_date == remdesivir_covid_therapeutics_last ~ "Remdesivir", + last_treatment_date == molnupiravir_covid_therapeutics_last ~ "Molnupiravir", + last_treatment_date == casirivimab_covid_therapeutics_last ~ "Casirivimab", + TRUE ~ NA_character_), + + # ELIGIBILITY VARIABLES ---- @@ -212,7 +237,7 @@ data_processed <- data_extract %>% sickle_cell_disease_nhsd = ifelse(!is.na(sickle_cell_disease_nhsd), 1, 0), ## Long COVID - long_covid = ifelse(!is.na(first_lc_code_date), 1, 0), + long_covid = ifelse(!is.na(first_long_covid), 1, 0), # Combine subgoups of rare neurological conditions cohort rare_neurological_conditions_nhsd = pmin(multiple_sclerosis_nhsd, motor_neurone_disease_nhsd, myasthenia_gravis_nhsd, @@ -561,6 +586,9 @@ data_processed_clean <- data_processed_combined %>% # Treatment paxlovid_covid_therapeutics, sotrovimab_covid_therapeutics, remdesivir_covid_therapeutics, molnupiravir_covid_therapeutics, casirivimab_covid_therapeutics, treatment_date, treatment_type, + + # Last treatment code recorded + last_treatment_date, last_treatment_type, # High risk cohort downs_syndrome, solid_cancer, haematological_disease, renal_disease, liver_disease, imid, immunosupression, diff --git a/analysis/study_definition.py b/analysis/study_definition.py index 93ebb36..c204538 100644 --- a/analysis/study_definition.py +++ b/analysis/study_definition.py @@ -97,6 +97,20 @@ "incidence": 0.05 }, ), + + paxlovid_covid_therapeutics_last = patients.with_covid_therapeutics( + #with_these_statuses = ["Approved", "Treatment Complete"], + with_these_therapeutics = "Paxlovid", + with_these_indications = "non_hospitalised", + on_or_after = "index_date", + find_last_match_in_period = True, + returning = "date", + date_format = "YYYY-MM-DD", + return_expectations = { + "date": {"earliest": "2022-02-10"}, + "incidence": 0.05 + }, + ), ## Sotrovimab sotrovimab_covid_therapeutics = patients.with_covid_therapeutics( @@ -112,6 +126,20 @@ "incidence": 0.2 }, ), + + sotrovimab_covid_therapeutics_last = patients.with_covid_therapeutics( + #with_these_statuses = ["Approved", "Treatment Complete"], + with_these_therapeutics = "Sotrovimab", + with_these_indications = "non_hospitalised", + on_or_after = "index_date", + find_last_match_in_period = True, + returning = "date", + date_format = "YYYY-MM-DD", + return_expectations = { + "date": {"earliest": "2021-12-16"}, + "incidence": 0.2 + }, + ), ## Remdesivir remdesivir_covid_therapeutics = patients.with_covid_therapeutics( @@ -127,6 +155,20 @@ "incidence": 0.2 }, ), + + remdesivir_covid_therapeutics_last = patients.with_covid_therapeutics( + #with_these_statuses = ["Approved", "Treatment Complete"], + with_these_therapeutics = "Remdesivir", + with_these_indications = "non_hospitalised", + on_or_after = "index_date", + find_last_match_in_period = True, + returning = "date", + date_format = "YYYY-MM-DD", + return_expectations = { + "date": {"earliest": "2021-12-16"}, + "incidence": 0.2 + }, + ), ### Molnupiravir molnupiravir_covid_therapeutics = patients.with_covid_therapeutics( @@ -143,6 +185,20 @@ }, ), + molnupiravir_covid_therapeutics_last = patients.with_covid_therapeutics( + #with_these_statuses = ["Approved", "Treatment Complete"], + with_these_therapeutics = "Molnupiravir", + with_these_indications = "non_hospitalised", + on_or_after = "index_date", + find_last_match_in_period = True, + returning = "date", + date_format = "YYYY-MM-DD", + return_expectations = { + "date": {"earliest": "2021-12-16"}, + "incidence": 0.2 + }, + ), + ### Casirivimab and imdevimab casirivimab_covid_therapeutics = patients.with_covid_therapeutics( #with_these_statuses = ["Approved", "Treatment Complete"], @@ -158,6 +214,19 @@ }, ), + casirivimab_covid_therapeutics_last = patients.with_covid_therapeutics( + #with_these_statuses = ["Approved", "Treatment Complete"], + with_these_therapeutics = "Casirivimab and imdevimab", + with_these_indications = "non_hospitalised", + on_or_after = "index_date", + find_last_match_in_period = True, + returning = "date", + date_format = "YYYY-MM-DD", + return_expectations = { + "date": {"earliest": "2021-12-16"}, + "incidence": 0.05 + }, + ), ## Date treated date_treated = patients.minimum_of( @@ -1394,30 +1463,14 @@ # # Long COVID - first_lc_dx_flag = patients.with_these_clinical_events( - long_covid_nice_dx, - on_or_before = "start_date", - returning = "code", - include_date_of_match = True, - return_first_date_in_period = True, - date_format = "YYYY-MM-DD", - return_expectations = {"incidence": 0.1}, - ), - - first_lc_dxrx_flag = patients.with_these_clinical_events( + first_long_covid = patients.with_these_clinical_events( long_covid_combine, on_or_before = "start_date", - returning = "binary _flag", - include_date_of_match = True, - return_first_date_in_period = True, + returning = "date", + find_first_match_in_period = True, date_format = "YYYY-MM-DD", return_expectations = {"incidence": 0.1}, ), - - first_lc_code_date = patients.minimum_of( - "first_lc_dx_flag_date","first_lc_dxrx_flag_date" - ), - # CLINICAL CO-MORBIDITIES TBC ---- diff --git a/project.yaml b/project.yaml index 81db89f..ed8e4cf 100644 --- a/project.yaml +++ b/project.yaml @@ -12,7 +12,7 @@ version: '3.0' expectations: - population_size: 500000 + population_size: 100000 actions: diff --git a/reports/coverage/mabs_and_antivirals_coverage_report.Rmd b/reports/coverage/mabs_and_antivirals_coverage_report.Rmd index 701c4f9..9ade705 100644 --- a/reports/coverage/mabs_and_antivirals_coverage_report.Rmd +++ b/reports/coverage/mabs_and_antivirals_coverage_report.Rmd @@ -39,7 +39,7 @@ library(ggplot2) library(png) library(readr) library(htmltools) -library(gt) +# library(gt) library(scales) library(lubridate) library(webshot) @@ -477,8 +477,7 @@ weekly_treatment_proportion_plot <- weekly_treatment_proportion %>% geom_line(size = 1) + facet_wrap(~high_risk_cohort, scales = "free_y") + theme_classic(base_size = 8) + - scale_x_date(limit = c(ymd("20211211"), ymd("20230701")), #end date updated to reflect most recent complete data - date_breaks = "2 months", date_labels = "%d %b %Y") + + scale_x_date(date_breaks = "2 months", date_labels = "%d %b %Y") + labs( x = "", y = "Proportion of eligble patients receiving treatment", @@ -862,7 +861,7 @@ proportion_stp_decile <- scale_color_manual(values = c("Decile" = "blue", "Median" = "blue")) + scale_linetype_manual(values = c("Decile" = "dashed", "Median" = "solid")) + theme_classic(base_size = 8) + - scale_x_date(limit = c(ymd("20211211"), ymd("20230701")), #end date updated to reflect most recent complete data + scale_x_date( date_breaks = "2 week", date_labels = "%d %b %Y") + labs( x = "Week", @@ -888,7 +887,62 @@ include_graphics(fs::path(output_dir_rmd, "figures", "figure_proportion_stp_deci ``` +
+ +```{r, treatment type count plot, echo=FALSE, message=FALSE, warning=FALSE, results='asis', fig.width=10, out.width = '100%', fig.cap = "**Figure 6 Total number of patients who received an antiviral or nMAB for treating COVID-19 since 11th December 2021, stratified by treatment type based on the most recent recorded treatment code recorded ** Shorter lines for Paxlovid and casirivimab reflect availability and guidance.", fig.topcaption=TRUE} + +treatment_codes_data <- + read_csv(here::here(input_dir_os, "table_last_treatment_codes_redacted.csv")) %>% + filter(!is.na(last_treatment_date), + last_treatment_date <= end_date, + count_redacted > 0) %>% + mutate(last_treatment_type = case_when(last_treatment_type == "Casirivimab" ~ "Casirivimab/imdevimab", + TRUE ~ last_treatment_type)) + +plot_order <- treatment_codes_data %>% + group_by(last_treatment_type) %>% + mutate(order = max(count_redacted, na.rm = T)) %>% + arrange(desc(order)) %>% + filter(count_redacted == order) %>% + select(last_treatment_type, order) %>% + distinct() +treatment_codes_data <- treatment_codes_data %>% + mutate(last_treatment_type = factor(last_treatment_type, levels = plot_order$last_treatment_type)) + +treatment_codes_plot <- treatment_codes_data %>% + ggplot(aes(x = last_treatment_date, y = count_redacted, colour = last_treatment_type, group = last_treatment_type)) + + geom_step(size = 1) + + theme_classic(base_size = 8) + + scale_x_date(date_breaks = "2 week", date_labels = "%d %b %Y") + + scale_y_continuous(limits = c(0, 200), labels = comma) + + labs( + x = "Most recent treatment date", + y = "Number of patients receiving treatment by most recent treatment date", + colour = "Treatment type",) + + theme( + axis.text = element_text(size = 12), + axis.text.x = element_text(angle = 60, hjust = 1), + axis.title = element_text(size = 12), + legend.text = element_text(size = 12), + legend.title=element_text(size=12, face = "bold"), + legend.position = c(0.1,0.75), + legend.background = element_rect(colour = "white"), + legend.box.background = element_rect(colour = "black"), + axis.line.x = element_line(colour = "black"), + panel.grid.minor.x = element_blank(), + legend.box.margin = margin(t = 1, l = 1, b = 1, r = 1), + plot.title = element_text(size=20, face="bold")) + +ggsave( + here::here(output_dir_rmd, "figures", "figure_last_treatment_code.png"), + treatment_codes_plot, + units = "cm", width = 35, height = 20 +) + +include_graphics(fs::path(output_dir_rmd, "figures", "figure_last_treatment_code.png")) + +```