WA_Can_HCA_Tbls.Rmd

---
title: "WA_Can_HCA_Tbls"
author: "Daniel Cockson"
date: "5/19/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

# Set working directory
setwd("C:/wa-cannabis-hca")

# Clear workspace 
rm(list = ls())

# Load required libraries 
library(survey)
library(epiR)
library(tidyverse)
library(foreign)
library(haven)
library(table1)
library(kableExtra)
library(gtsummary)
library(gt)
library(broom)

# Read in the cleaned data
load("data/WA_2015_2019_BRFSS_Clean.Rds")

```

```{r}
#### Creating Table 1 and other exploratory tables ####

# Change to row percent!

db <- svydesign(
  id = ~1,
  strata = ~X_ststr,
  weights = ~X_llcpwt,
  data = df, 
  nest = TRUE) 

# main table 1 by current vs non-current cannabis use 
t1 <- tbl_svysummary(db, by = mj30day, percent = "row", digits = everything() ~ 0,
               include = c(X_age_g, sex, race, edu, income_wa, dHCA, mentComp, curAlc, curSmk),
               statistic = list(everything() ~ "{n_unweighted} ({p}%)")) %>%
  modify_header(all_stat_cols(FALSE) ~ "**{level}**, N={n_unweighted}") %>%
  modify_caption("**Table 1. Past Month Cannabis User Characteristics**") %>%
  bold_labels()

t1 %>%
  as_gt() %>%
  gtsave("Table1_CA_WA.html")

# alternative table 1's for demographic exploration
t1.alt.medcost <- tbl_svysummary(db, by = dHCA, percent = "row", digits = everything() ~ 0,
               include = c(X_age_g, sex, race, edu, income_wa, mj30day, mentComp, curAlc, curSmk),
               statistic = list(everything() ~ "{n_unweighted} ({p}%)")) %>%
  modify_header(stat_by = "**{level}**, N={n_unweighted}") %>%
  modify_caption("**Table 1Alt. dHCA Characteristics**") %>%
  bold_labels()

t1.alt.medcost %>%
  as_gt() %>%
  gtsave("Table1_Med_CA_WA.html")


t1.alt.ment14d <- tbl_svysummary(db, by = mentComp, percent = "row", digits = everything() ~ 0,
               include = c(X_age_g, sex, race, edu, income_wa, dHCA, mj30day, curAlc, curSmk),
               statistic = list(everything() ~ "{n_unweighted} ({p}%)")) %>%
  modify_header(stat_by = "**{level}**, N={n_unweighted}") %>%
  modify_caption("**Table 1Alt. Mental Health Characteristics**") %>%
  bold_labels()

t1.alt.ment14d %>%
  as_gt() %>%
  gtsave("Table1_Ment_CA_WA.html")

```

```{r}
#### Creating Table 2 ####

# Crude state wide PR
# create an empty dataframe for the adjusted PR's
adj.PR <- NULL

# Find the crude PR
adj.PR <- rbind(adj.PR,
                tidy(epi.2by2(table(df$mentComp, df$mj30day), method = "cross.sectional")) %>%
                  filter(term == "PR.strata.wald")) 

adj.PR[1,1] <- "Crude" # Assign variable label to the initialized row

# Creating confounding variable
df$inc_sex_smk <- case_when(!is.na(df$sex) & 
                              !is.na(df$curSmk) &
                              !is.na(df$income_wa) ~
                               paste0(df$sex, "_", df$curSmk, "_", df$income_wa)) 

# Create a table with the associated confounding variable
combined.tbl <- with(df,
                     table(mentComp, mj30day, inc_sex_smk))

# Find the adjusted PR
adj.PR <- rbind(adj.PR,
                tidy(epi.2by2(combined.tbl, method = "cross.sectional")) %>%
                  filter(term == "PR.mh.wald")) 

adj.PR[2,1] <- "Adjusted" # Assign variable label to the second row

```

```{r}
#### Creating the table 2 ####

adj.PR <- adj.PR %>%
  gt(rowname_col = "term") %>% 
  tab_footnote(
    footnote = "Prevalance ratio adjusted for gender, income, and current tobacco smoking status.",
    locations = cells_body(
      columns = estimate,
      rows = 2)
  ) %>%
  fmt_number(
    columns = 2:4,
    decimals = 2
  ) %>%
  tab_header(
    title = md("Table 2: Prevalance Ratios of Cannabis Usage"),
    subtitle = md("14+ Poor Mental Health Days vs. 0-13 Days of Poor Mental Health Days")) %>%
  tab_spanner(label = "95% Confidence Interval", columns = c(conf.low, conf.high)) %>%
  cols_label(
    estimate = "Prevalance Ratio",
    conf.low = "Lower",
    conf.high = "Upper") 

# Save the table
adj.PR %>%
  gtsave("Table2_CA_WA.html")



```

```{r}
#### Table 3: Supplemental R/E analysis ####

# create empty dataframe to store variables
cost.PR <- NULL

# Get the dHCA strata specific PR's
ovr.strat.y <- with(subset(df, dHCA == "Yes"),
          table(mentComp, mj30day, income_wa))
ovr.strat.n <- with(subset(df, dHCA == "No"),
          table(mentComp, mj30day, income_wa))

cost.PR <- rbind(cost.PR,
                tidy(epi.2by2(ovr.strat.y, method = "cross.sectional")) %>%
                  filter(term == "PR.mh.wald")) 
cost.PR <- rbind(cost.PR,
                tidy(epi.2by2(ovr.strat.n, method = "cross.sectional")) %>%
                  filter(term == "PR.mh.wald")) 

# # get the dHCA strata specific PR's by wt.
# asian.strat <- with(subset(df, race == "Non-Hispanic Asian"),
#                     table(mentComp, mj30day, dHCA))
# cost.PR <- rbind(cost.PR,
#                 tidy(epi.2by2(asian.strat, method = "cross.sectional")) %>%
#                   filter(term == "PR.strata.wald")) 
# 
# # get the dHCA strata specific PR's by wt.
# blk.strat <- with(subset(df, race == "Non-Hispanic Black or African American"),
#                     table(mentComp, mj30day, dHCA))
# cost.PR <- rbind(cost.PR,
#                 tidy(epi.2by2(blk.strat, method = "cross.sectional")) %>%
#                   filter(term == "PR.strata.wald")) 

# get the dHCA strata specific PR's by hisp.
hisp.strat.y <- with(subset(df, dHCA == "Yes" & race == "Hispanic, Latino/a, or Spanish origin"),
          table(mentComp, mj30day, income_wa))
hisp.strat.n <- with(subset(df, dHCA == "No" & race == "Hispanic, Latino/a, or Spanish origin"),
          table(mentComp, mj30day, income_wa))

cost.PR <- rbind(cost.PR,
                tidy(epi.2by2(hisp.strat.y, method = "cross.sectional")) %>%
                  filter(term == "PR.mh.wald")) 
cost.PR <- rbind(cost.PR,
                tidy(epi.2by2(hisp.strat.n, method = "cross.sectional")) %>%
                  filter(term == "PR.mh.wald")) 

# get the dHCA strata specific PR's by wt.
wt.strat.y <- with(subset(df, dHCA == "Yes" & race == "Non-Hispanic White"),
          table(mentComp, mj30day, income_wa))
wt.strat.n <- with(subset(df, dHCA == "No" & race == "Non-Hispanic White"),
          table(mentComp, mj30day, income_wa))

cost.PR <- rbind(cost.PR,
                tidy(epi.2by2(wt.strat.y, method = "cross.sectional")) %>%
                  filter(term == "PR.mh.wald")) 
cost.PR <- rbind(cost.PR,
                tidy(epi.2by2(wt.strat.n, method = "cross.sectional")) %>%
                  filter(term == "PR.mh.wald")) 

# creating headers and lables
cost.PR[c("1","3","5"),1] <- "Difficulty accessing healthcare"
cost.PR[c("2","4","6"),1] <- "Did not have difficulty accessing healthcare"

cost.PR$new <- NA # Add blank column
cost.PR <- cost.PR[,c(5,1:4)] # Reorder columns

cost.PR$new[1:2] <- "All adults"
# cost.PR$new[3:4] <- "Non-Hispanic Asian"
# cost.PR$new[5:6] <- "Non-Hispanic Black or African American"
cost.PR$new[3:4] <- "Hispanic, Latino/a, or Spanish origin"
cost.PR$new[5:6] <- "Non-Hispanic White"

```

```{r}
#### Creating the table 3 ####

cost.PR <- cost.PR %>%
  gt(rowname_col = "term", groupname_col = "new") %>%  
  fmt_number(
    columns = 3:5,
    decimals = 2
  ) %>%
  tab_header(
    title = md("Table 3: Prevalance Ratios of Cannabis Usage by Race/Ethnicity Adjusted for Income"),
    subtitle = md("14+ Poor Mental Health Days vs. 0-13 Days of Poor Mental Health Days")) %>%
  tab_spanner(label = "95% Confidence Interval", columns = c(conf.low, conf.high)) %>%
  cols_label(
    estimate = "Prevalance Ratio",
    conf.low = "Lower",
    conf.high = "Upper") %>%
  tab_style(
    style = list(cell_text(weight = "bold")),
    locations = cells_row_groups()
  )

# Save the table
cost.PR %>%
  gtsave("Table3_CA_WA.html")

```

```{r}

d <- as.data.frame(df$mjpast30[df$mjpast30 != 88])
d <- na.omit(d)

ci(d$`df$mjpast30[df$mjpast30 != 88]`, confidence = 0.95)


ovr.strat.y <- with(subset(df, dHCA == "Yes"),
          table(curSmk))
ovr.strat.n <- with(subset(df, dHCA == "No"),
          table(curSmk))

```