Using R for Reproducible and Transparent Health Research

bas_bp_summary <- bas_bp %>%
  dplyr::select(
    `Treatment Protocol`,
    `Age, Mean (SD), years`,
    `Sex, No. (%)`,
    `Body mass index, Mean (SD), kg/m2`,
    #`Body mass index class, No. (%)`, 
    `Education, No. (%)`,
    `History of hypertension, No. (%)`,
    `History of diabetes, No. (%)`,
    `History of chronic kidney disease, No. (%)`,
    `History of stroke, No. (%)`,
    `History of heart attack, No. (%)`,
    `Smoker/Tobacco user, No. (%)`, 
    `Alcohol user, No. (%)`,
    `Systolic blood pressure, Mean (SD), mmHg`,
    `Diastolic blood pressure, Mean (SD), mmHg`,
    `Hypertension Control, No. (%)`,
    `Heart rate, Mean (SD), bpm`,
    `Baseline treatment rate, No. (%)`
  )

table_1 <- tbl_summary(bas_bp_summary,
                       by = `Treatment Protocol`,
                       missing = "no",
                       digits = list(all_categorical() ~ c(0, 1), all_continuous() ~ c(0, 1)),
                       statistic = list(all_continuous() ~ "{mean} ({sd})",
                                        all_categorical() ~ "{n} ({p})")) %>% 
  # add_n() %>% 
  add_p(pvalue_fun = ~ style_pvalue(.x, digits = 2)) %>% 
  add_overall() %>% 
  modify_header(label = "**Variables**") %>% 
  bold_p() %>% 
  modify_caption("**Table 1**. Patient baseline demographics, clinical characteristics and blood pressure overall and by study arm.") %>% 
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Treatment Protocol Received**") %>%   modify_footnote(
    update = everything() ~ NA
  )

table_1 

# summary data on htn control
htnc_b <- dat_final %>% 
  group_by(id) %>%
  #filter(time > 2) %>%
  filter(time == min(time)) %>% 
  group_by(study_arm) %>%
  count(htn_c) %>%
  mutate(total = sum(n)) %>% 
  filter(htn_c == "1") %>%
   mutate( #total = sum(n),
           rate = map2(n, total, ~ prop.test(.x, .y, conf.level = 0.95) %>%
                   broom::tidy()),
           htn_c = factor(htn_c, levels = c("0", "1"), labels = c("Not controlled", "Controlled"))) %>%
  unnest(rate) %>%
  filter(htn_c == "Controlled") %>% 
  mutate(
    n_base = n,
    total_base = total,
    estimate_base = estimate
  )

htn_control_baseline_itt <- htnc_b %>% 
  dplyr::select(
    study_arm,
    n_base,
    total_base,
    estimate_base
  )


htn_control_finale_itt <- htnc_b %>% 
  dplyr::select(
    study_arm,
    n_finale = n,
    total_finale = total,
    estimate_finale = estimate
  )



## selecting people who attended the PHC twice and are controlled at the last visit
htnc_new <- dat_final %>% 
  group_by(id) %>%
  filter(time > 2) %>%
  filter(time == max(time)) %>%
  group_by(study_arm) %>%
  count(htn_c) %>% 
  mutate(total = sum(n)) %>% 
  filter(htn_c == "1") %>% 
   mutate( #total = sum(n),
           rate = map2(n, total, ~ prop.test(.x, .y, conf.level = 0.95) %>%
                   broom::tidy()),
           htn_c = factor(htn_c, levels = c("0", "1"), labels = c("Not controlled", "Controlled"))) %>%
  unnest(rate) %>%
  filter(htn_c == "Controlled") %>% 
   mutate(
    n_finale = n,
    total_finale = total,
    estimate_finale = estimate
  )


htn_control_baseline <- htnc_b %>% 
  dplyr::select(
    study_arm,
    n_base,
    total_base,
    estimate_base
  )

htn_control_finale <- htnc_new %>% 
  dplyr::select(
    study_arm,
    n_finale = n,
    total_finale = total,
    estimate_finale = estimate
  )



## difference between hypertension control at baseline and at last visit
## let's join the database of htn control at baseline and at the last visit. 

htn_control_diff <- htn_control_baseline %>% 
  left_join(htn_control_finale, by = "study_arm") %>% 
  mutate(
    htn_diff = estimate_finale - estimate_base
  )

## let's compare the proportion of htn control in baseline and at last visit in each study arm









controlled_least_twice <- ggplot(data = htnc_new, aes(x = study_arm, y = estimate)) + 
  geom_col(aes(fill = htn_c), position = "dodge", show.legend = FALSE) + 
    geom_errorbar(
    aes(ymin = conf.low,
        ymax = conf.high),
     position = position_dodge(0.9),
     width = 0.2
  ) + 
  geom_text(aes(label = scales::percent(estimate), 
                y = estimate,
                group = htn_c),
            position = position_dodge(width = 0.9),
            vjust = -2.1) +
  scale_color_jama() +
  #scale_fill_manual(values=c("grey70", "grey70")) + 
  labs(
    x = "Treatment Protocol",
    y = "Percentage (%) of HTN Controlled at last visit",
    fill = "Hypertension"
  ) +  
   scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
   theme_bw() +
  geom_signif(comparisons=list(c("1", "2")), annotations="Cluster adjusted p-value = 0.29",
              y_position = 0.68,
              tip_length = 0.03, 
              vjust = -0.005
              )
  
  

controlled_least_twice