analysis_2.Rmd

---
title: "analysis_1"
author: "vaishali"
date: "October 31, 2018"
output: word_document
---

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

## Analysis 1

```{r}
library(haven);library(survey);library(tidyverse);library(tableone)
```

Get the dataframes loaded into the document.

```{r get df for each year seperately}
d07 <- read_csv("H:/bita_nis/df/m2007.csv")
d08 <- read_csv("H:/bita_nis/df/m2008.csv")
d09 <- read_csv("H:/bita_nis/df/m2009.csv")
d10 <- read_csv("H:/bita_nis/df/m2010.csv")
d12 <- read_csv("H:/bita_nis/df/m2012.csv")
d11 <- read_csv("H:/bita_nis/df/m2011.csv")
d13 <- read_csv("H:/bita_nis/df/m2013.csv")
d14 <- read_csv("H:/bita_nis/df/m2014.csv")
d15.1 <- read_csv("H:/bita_nis/df/m2015_1.csv")


# would want to add the amonth variable to this df so that we can seperate into quarters.

a07 <- read_csv("H:/bita_nis/df/a07.csv")
a08 <- read_csv("H:/bita_nis/df/a08.csv")
a09 <- read_csv("H:/bita_nis/df/a09.csv")
a10 <- read_csv("H:/bita_nis/df/a10.csv")
a11 <- read_csv("H:/bita_nis/df/a11.csv")
a12 <- read_csv("H:/bita_nis/df/a12.csv")
a13 <- read_csv("H:/bita_nis/df/a13.csv")
a14 <- read_csv("H:/bita_nis/df/a14.csv")
a15 <- read_csv("H:/bita_nis/df/a15.csv")

# reduce to lower case

names(a07) <- tolower(names(a07))
names(a08) <- tolower(names(a08))
names(a09) <- tolower(names(a09))
names(a10) <- tolower(names(a10))
names(a11) <- tolower(names(a11))
names(a12) <- tolower(names(a12))
names(a13) <- tolower(names(a13))
names(a14) <- tolower(names(a14))
names(a15) <- tolower(names(a15))


# join to add amonth to the df

df07 <- left_join(d07, a07, by = "key")
df08 <- left_join(d08, a08, by = "key")
df09 <- left_join(d09, a09, by = "key")
df10 <- left_join(d10, a10, by = "key")
df11 <- left_join(d11, a11, by = "key")
df12 <- left_join(d12, a12, by = "key_nis")
df13 <- left_join(d13, a13, by = "key_nis")
df14 <- left_join(d14, a14, by = "key_nis")
df15.1 <- left_join(d15.1, a15, by = "key_nis")





# keep hosp_locteach as that is a more accurate variable of teaching status. The hosp_teach variable is not included in the df after 2011. hosp_locteach is available from 1998 - 2016; so we should use hosp_locteach



df12$hospid <- df12$hosp_nis
df12$key <- df12$key_nis

df12$trendwt <- df12$discwt.x

df13$hospid <- df13$hosp_nis
df13$key <- df13$key_nis

df13$trendwt <- df13$discwt.x

df14$hospid <- df14$hosp_nis
df14$key <- df14$key_nis

df14$trendwt <- df14$discwt.x

df15.1$hospid <- df15.1$hosp_nis
df15.1$key <- df15.1$key_nis

df15.1$trendwt <- df15.1$discwt.x

# load the 15_4 quarter df here

df15.4 <- read_csv("H:/bita_nis/df/m2015_4.csv")

df16 <- read_csv("H:/bita_nis/df/m2016.csv")

```

Am going to rbind columns according to need to do the survey analysis. 
Set lonely.psu = certainty to keep 1 PSU units

```{r set lonely psu to average}
options("survey.lonely.psu" = 'average')
```


```{r}
names(df07)
```

# Plan to present data as quarters. 

# create another Rmarkdown document to prepare the `amonth` and `key` or `key_nis` dataframes.


```{r}

b07 <- df07 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, trendwt, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1,dx1:dx10, pr1:pr10,discwt.x,amonth)

b08 <- df08 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, trendwt, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1,dx1:dx10, pr1:pr10,discwt.x, amonth)

b09 <- df09 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, trendwt, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1, dx1:dx10, pr1:pr10,discwt.x, amonth)

b10 <- df10 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, trendwt, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1, dx1:dx10, pr1:pr10,discwt.x, amonth)

b11 <- df11 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, trendwt, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1, dx1:dx10, pr1:pr10,discwt.x, amonth)

b12 <- df12 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race, trendwt,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1, dx1:dx10, pr1:pr10,discwt.x, amonth)


b13 <- df13 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race, trendwt,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1,dx1:dx10, pr1:pr10,discwt.x, amonth)

b14 <- df14 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race, trendwt,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1, dx1:dx10, pr1:pr10,discwt.x, amonth)

b15.1 <- df15.1 %>% select(age, bita, died, year.x, dispuniform, female, key, los, nis_stratum.x, valve, priormi, priorpci, chf,  stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, race, trendwt,cm_dm, cm_dmcx, cm_obese, cm_renlfail, pay1, dx1:dx10, pr1:pr10,discwt.x, amonth)
```

Now, we need to create var in years 2012 onwards with same names that are present in 2011 so that we can rbind all data together.

```{r cocactenate all df together}


df <- rbind(b07,b08,b09,b10,b12,b11,b13,b14,b15.1)

# remove patients who had possible single vessel CABG

cabg_n <- as.character(c('3610','3612','3613','3614','3615','3617','3619'))

df$cabg_n <- with(df,ifelse((pr1 %in% cabg_n | pr2 %in% cabg_n | pr3 %in% cabg_n |
                                   pr4 %in% cabg_n | pr5 %in% cabg_n),"yes","no"))



df %>% count(cabg_n)


df2 <- df %>% filter(cabg_n == "yes")

df <- df2


# Insulin dependent diabetics

iddm <- as.character(c('V5867'))


df$iddm_n <- with(df,ifelse((dx1 %in% iddm | dx2 %in% iddm | dx3 %in% iddm | dx4 %in% iddm | dx5 %in% iddm|
                               dx6 %in% iddm| dx7 %in% iddm| dx8 %in% iddm| dx9 %in% iddm| dx10 %in% iddm),"yes","no"))

df %>% count(iddm_n)


media <- as.character("5192")

df$media <- with(df,ifelse((dx1 %in% media  | dx2 %in% media | dx3 %in% media |
                                   dx4 %in% media | dx5 %in% media | dx6 %in% media | dx7 %in% media | dx8 %in% media | dx9 %in% media | dx10 %in% media),"yes","no"))

df %>% count(media)


inf <- as.character("99859")

df$inf <- with(df,ifelse((dx1 %in% inf  | dx2 %in% inf | dx3 %in% inf |
                                   dx4 %in% inf | dx5 %in% inf | dx6 %in% inf | dx7 %in% inf | dx8 %in% inf | dx9 %in% inf | dx10 %in% inf),"yes","no"))

df %>% count(inf)



cs <- as.character(c(78551))

df$cs <- with(df,ifelse((dx1 %in% cs  | dx2 %in% cs | dx3 %in% cs |
                                   dx4 %in% cs | dx5 %in% cs | dx6 %in% cs | dx7 %in% cs | dx8 %in% cs | dx9 %in% cs | dx10 %in% cs),"yes","no"))


df %>% count(cs)

cpb <- as.character("3961")

df$cpb <- with(df,ifelse((pr1 %in% cpb | pr2 %in% cpb | pr3 %in% cpb |
                                   pr4 %in% cpb | pr5 %in% cpb | pr6 %in% cpb | pr7 %in% cpb | pr8 %in% cpb | pr9 %in% cpb | pr10 %in% cpb),"yes","no"))

df %>% count(cpb)


# use of iabp

iabp <- as.character(c(3761))

df$iabp_n <- with(df,ifelse((pr1 %in% iabp | pr2 %in% iabp | pr3 %in% iabp |
                                   pr4 %in% iabp | pr5 %in% iabp),"yes","no"))



df %>% count(iabp_n)


# use of pvad

pvad <- as.character(c(3768))

df$pvad_n <- with(df,ifelse((pr1 %in% pvad | pr2 %in% pvad | pr3 %in% pvad |
                                   pr4 %in% pvad | pr5 %in% pvad),"yes","no"))



df %>% count(pvad_n)

## create quarters for each year in the df dataframe

df$q_n[df$amonth < 4]<- "Q1"
df$q_n[df$amonth>3 & df$amonth <7]<- "Q2"
df$q_n[df$amonth >6 & df$amonth < 10]<- "Q3"
df$q_n[df$amonth >9] <- "Q4"

df %>% count(q_n)


```

Make some changes and add the variables needed to the ICD 10 containing df:


```{r}

# cpb 

cpb <- as.character(c("5A1221Z"))

a <- cpb

df15.4$cpb <- with(df15.4, ifelse((i10_pr1 %in% a | i10_pr2 %in% a | i10_pr3 %in% a | i10_pr4 %in% a  | i10_pr5 %in% a | i10_pr6 %in% a | i10_pr7 %in% a | i10_pr8 %in% a |i10_pr9 %in% a| i10_pr10 %in% a | i10_pr11 %in% a | i10_pr12 %in% a | i10_pr13 %in% a | i10_pr14 %in% a | i10_pr15 %in% a ), "yes","no"))

df15.4 %>% count(cpb)


# cpb 

df16$cpb <- with(df16, ifelse((i10_pr1 %in% a | i10_pr2 %in% a | i10_pr3 %in% a | i10_pr4 %in% a  | i10_pr5 %in% a | i10_pr6 %in% a | i10_pr7 %in% a | i10_pr8 %in% a |i10_pr9 %in% a| i10_pr10 %in% a | i10_pr11 %in% a | i10_pr12 %in% a | i10_pr13 %in% a | i10_pr14 %in% a | i10_pr15 %in% a ), "yes","no"))

df16 %>% count(cpb)



# age groups

df$young[df$age < 51]<- "yes"
df$young[df$age >50]<- "no"


# octogenerians

df$old[df$age < 80]<- "no"
df$old[df$age > 79]<- "yes"



# select smaller df with some variables 

# 2015 last quarter


b15.4 <- df15.4 %>% dplyr::select(drg,hosp_nis.x,i10_dx1,i10_dx2,i10_dx3,i10_dx4,i10_dx5,i10_dx6         ,i10_dx7,i10_dx8 ,i10_dx9,i10_dx10,i10_dx11,i10_dx12,i10_dx13,i10_dx14,i10_dx15,i10_dx16,i10_dx17        ,i10_dx18 ,i10_dx19, i10_dx20,i10_dx21,i10_dx22,i10_dx23,i10_dx24,i10_dx25,i10_dx26,i10_dx27        ,i10_dx28,i10_dx29,i10_dx30,i10_pr1,i10_pr2,i10_pr3,i10_pr4,i10_pr5,i10_pr6,i10_pr7,i10_pr8,i10_pr9, i10_pr10, i10_pr11, i10_pr12,i10_pr13,i10_pr14,i10_pr15,key_nis,discwt.x,hosp_bedsize, hosp_division.x, hosp_locteach, hosp_region,nis_stratum.x ,age,died,discwt.y,dispuniform,elective,female         ,hosp_division.y,los,pay1,race,totchg,lita,rita,bita,cm_chf,cm_aids,cm_anemdef,priormi,priorpci,stemi,
                                  carotid,icd,cpb, year.x, amonth)


# 2016 

b16 <- df16 %>% dplyr::select(drg,hosp_nis.x,i10_dx1,i10_dx2,i10_dx3,i10_dx4,i10_dx5,i10_dx6         ,i10_dx7,i10_dx8 ,i10_dx9,i10_dx10,i10_dx11,i10_dx12,i10_dx13,i10_dx14,i10_dx15,i10_dx16,i10_dx17        ,i10_dx18 ,i10_dx19, i10_dx20,i10_dx21,i10_dx22,i10_dx23,i10_dx24,i10_dx25,i10_dx26,i10_dx27        ,i10_dx28,i10_dx29,i10_dx30,i10_pr1,i10_pr2,i10_pr3,i10_pr4,i10_pr5,i10_pr6,i10_pr7,i10_pr8,i10_pr9, i10_pr10, i10_pr11, i10_pr12,i10_pr13,i10_pr14,i10_pr15,key_nis,discwt.x,hosp_bedsize, hosp_division.x, hosp_locteach, hosp_region,nis_stratum.x ,age,died,discwt.y,dispuniform,elective,female         ,hosp_division.y,los,pay1,race,totchg,lita,rita,bita,cm_chf,cm_aids,cm_anemdef,priormi,priorpci,stemi,
                                  carotid,icd,cpb, year.x, amonth)



# need to do the same as below for the `var` of interest in this df
# 2015 Q4

b15.4 <- b15.4 %>% 
  mutate(race_n = fct_recode(factor(race),
                                  "white" = "1",
                                   "black" = "2",
                                   "hispanic" = "3",
                              "others" = "4",
                              "others" = "5",
                             NULL = "6" ))


b15.4 <- b15.4 %>% 
  mutate(pay_n = fct_recode(factor(pay1),
                            "medicare" = "1",
                            "medicaid" = "2",
                            "private_insurance" = "3",
                            "self_pay" = "4",
                            "other" = "5",
                            "other" = "6"))



b15.4 %>% count(race_n)

b15.4 %>% count(pay_n)

# 2016

b16 <- b16 %>% 
  mutate(race_n = fct_recode(factor(race),
                                  "white" = "1",
                                   "black" = "2",
                                   "hispanic" = "3",
                              "others" = "4",
                              "others" = "5",
                             NULL = "6" ))


b16 <- b16 %>% 
  mutate(pay_n = fct_recode(factor(pay1),
                            "medicare" = "1",
                            "medicaid" = "2",
                            "private_insurance" = "3",
                            "self_pay" = "4",
                            "other" = "5",
                            "other" = "6"))



b16 %>% count(race_n)

b16 %>% count(pay_n)


# rbind both years together

b2 <- rbind(b15.4, b16)




b2bita <- b15.4 %>% filter(bita == "yes")

# need to create surveydesign object for this df 
# Would now combine 2015 Q4 and 2016 to create the dataset and survey object for them together.

s2 <- survey::svydesign(data = b2, weights = ~discwt.y, strata = ~nis_stratum.x, id = ~hosp_nis.x, nest = TRUE)

summary(s2)


```




```{r some more data munging...}
# convert some variables to better form for analysis

# race 
df <- df %>% 
  mutate(race_n = fct_recode(factor(race),
                                  "white" = "1",
                                   "black" = "2",
                                   "hispanic" = "3",
                              "others" = "4",
                              "others" = "5",
                             NULL = "6" ))


# primary payer

df <- df %>% 
  mutate(pay_n = fct_recode(factor(pay1),
                            "medicare" = "1",
                            "medicaid" = "2",
                            "private_insurance" = "3",
                            "self_pay" = "4",
                            "other" = "5",
                            "other" = "6"))

dfbita <- df %>% filter(bita == "yes")


# disposition after dismissal

df <- df %>% 
  mutate(disp_n = fct_recode(factor(dispuniform),
                            "home" = "1",
                            "short term hospital" = "2",
                            "SNF" = "5",
                            "home health care" = "6",
                            "AMA" = "7",
                            "died" = "20"))

df$disp_n[df$disp_n == "21"]<- NA
df$disp_n[df$disp_n == "99"]<- NA

df %>% count(disp_n)

df$bita_i[df$bita == "yes"] <- 1
df$bita_i[df$bita == "no"]<- 0

svydf <- survey::svydesign(data = df, weights = ~discwt.x, strata = ~nis_stratum.x, id = ~hospid, nest = T)



summary(svydf)

svydfbita <- survey::svydesign(data = dfbita, weights = ~discwt.x, strata = ~nis_stratum.x, id = ~hospid, nes = T)

summary(svydfbita)

```


After the data is obtained and cleaned, here are the trends that we want to see:

We want to divide this paper into 2 parts:

1. Trends overall as a country
2. Differences between regions and trends between regions.

1. Trends overall as a country ---- 

a. % done in large hospitals vs rural hospitals == any change ?
b. % done in teaching vs non-teaching hospitals == any change ?
c. BITA use/ year == any change ?
d. OPCABG/ year == any change ?
e. mortality == any change ?
f. LOS == any change ?
g. disposition == any change ?


If we want to focus on only BITA use ?

a. Use of BITA/ year in whole US
b. Use of BITA according to hospital type 
c. Use of BITA according to teaching status 
d. Use of BITA according to CABG volume from 2007 - 2011. But even after, we can see total volume reported and % of BITA used according to that. So, we may be able to determine if there is a relationship between hospital volume and BITA use according. 


```{r add some cm_ to get all the data according to earlier years}

# liver

liver <- as.character(c(
  
  "B180",
  "B181",
  "B182",
  "I8500",
  "I8501",
  "I8510",
  "I8511",
  "K700",
  "K702",
  "K7030",
  "K7031",
  "K7040",
  "K7041",
  "K709",
  "K7210",
  "K7211",
  "K7290",
  "K7291",
  "K730",
  "K731",
  "K732",
  "K738",
  "K739",
  "K740",
  "K741",
  "K742",
  "K743",
  "K744",
  "K745",
  "K7460",
  "K7469",
  "K754",
  "K7581",
  "K760",
  "K766",
  "K7689",
  "K769",
  "Z944"))


a <- liver

b2$liver <- with(b2, ifelse((i10_dx1 %in% a | i10_dx2 %in% a | i10_dx3 %in% a | i10_dx4 %in% a  | i10_dx5 %in% a | i10_dx6 %in% a | i10_dx7 %in% a | i10_dx8 %in% a |i10_dx9 %in% a| i10_dx10 %in% a | i10_dx11 %in% a | i10_dx12 %in% a | i10_dx13 %in% a | i10_dx14 %in% a | i10_dx15 %in% a ), "yes","no"))

b2 %>% count(liver)

```
```{r}
dm <- as.character(c(
  "E0800",
  "E0801",
  "E0810",
  "E0811",
  "E089",
  "E0900",
  "E0901",
  "E0910",
  "E0911",
  "E099",
  "E1010",
  "E1011",
  "E109",
  "E1100",
  "E1101",
  "E1110",
  "E1111",
  "E119",
  "E1300",
  "E1301",
  "E1310",
  "E1311",
  "E139",
  "O24011",
  "O24012",
  "O24013",
  "O24019",
  "O2402",
  "O2403",
  "O24111",
  "O24112",
  "O24113",
  "O24119",
  "O2412",
  "O2413",
  "O24311",
  "O24312",
  "O24313",
  "O24319",
  "O2432",
  "O2433",
  "O24811",
  "O24812",
  "O24813",
  "O24819",
  "O2482",
  "O2483",
  "O24911",
  "O24912",
  "O24913",
  "O24919",
  "O2492",
  "O2493"))

dmcx <- as.character(c(
  
  "E0821",
  "E0822",
  "E0829",
  "E08311",
  "E08319",
  "E08321",
  "E083211",
  "E083212",
  "E083213",
  "E083219",
  "E08329",
  "E083291",
  "E083292",
  "E083293",
  "E083299",
  "E08331",
  "E083311",
  "E083312",
  "E083313",
  "E083319",
  "E08339",
  "E083391",
  "E083392",
  "E083393",
  "E083399",
  "E08341",
  "E083411",
  "E083412",
  "E083413",
  "E083419",
  "E08349",
  "E083491",
  "E083492",
  "E083493",
  "E083499",
  "E08351",
  "E083511",
  "E083512",
  "E083513",
  "E083519",
  "E083521",
  "E083522",
  "E083523",
  "E083529",
  "E083531",
  "E083532",
  "E083533",
  "E083539",
  "E083541",
  "E083542",
  "E083543",
  "E083549",
  "E083551",
  "E083552",
  "E083553",
  "E083559",
  "E08359",
  "E083591",
  "E083592",
  "E083593",
  "E083599",
  "E0836",
  "E0837X1",
  "E0837X2",
  "E0837X3",
  "E0837X9",
  "E0839",
  "E0840",
  "E0841",
  "E0842",
  "E0843",
  "E0844",
  "E0849",
  "E0851",
  "E0852",
  "E0859",
  "E08610",
  "E08618",
  "E08620",
  "E08621",
  "E08622",
  "E08628",
  "E08630",
  "E08638",
  "E08641",
  "E08649",
  "E0865",
  "E0869",
  "E088",
  "E0921",
  "E0922",
  "E0929",
  "E09311",
  "E09319",
  "E09321",
  "E093211",
  "E093212",
  "E093213",
  "E093219",
  "E09329",
  "E093291",
  "E093292",
  "E093293",
  "E093299",
  "E09331",
  "E093311",
  "E093312",
  "E093313",
  "E093319",
  "E09339",
  "E093391",
  "E093392",
  "E093393",
  "E093399",
  "E09341",
  "E093411",
  "E093412",
  "E093413",
  "E093419",
  "E09349",
  "E093491",
  "E093492",
  "E093493",
  "E093499",
  "E09351",
  "E093511",
  "E093512",
  "E093513",
  "E093519",
  "E093521",
  "E093522",
  "E093523",
  "E093529",
  "E093531",
  "E093532",
  "E093533",
  "E093539",
  "E093541",
  "E093542",
  "E093543",
  "E093549",
  "E093551",
  "E093552",
  "E093553",
  "E093559",
  "E09359",
  "E093591",
  "E093592",
  "E093593",
  "E093599",
  "E0936",
  "E0937X1",
  "E0937X2",
  "E0937X3",
  "E0937X9",
  "E0939",
  "E0940",
  "E0941",
  "E0942",
  "E0943",
  "E0944",
  "E0949",
  "E0951",
  "E0952",
  "E0959",
  "E09610",
  "E09618",
  "E09620",
  "E09621",
  "E09622",
  "E09628",
  "E09630",
  "E09638",
  "E09641",
  "E09649",
  "E0965",
  "E0969",
  "E098",
  "E1021",
  "E1022",
  "E1029",
  "E10311",
  "E10319",
  "E10321",
  "E103211",
  "E103212",
  "E103213",
  "E103219",
  "E10329",
  "E103291",
  "E103292",
  "E103293",
  "E103299",
  "E10331",
  "E103311",
  "E103312",
  "E103313",
  "E103319",
  "E10339",
  "E103391",
  "E103392",
  "E103393",
  "E103399",
  "E10341",
  "E103411",
  "E103412",
  "E103413",
  "E103419",
  "E10349",
  "E103491",
  "E103492",
  "E103493",
  "E103499",
  "E10351",
  "E103511",
  "E103512",
  "E103513",
  "E103519",
  "E103521",
  "E103522",
  "E103523",
  "E103529",
  "E103531",
  "E103532",
  "E103533",
  "E103539",
  "E103541",
  "E103542",
  "E103543",
  "E103549",
  "E103551",
  "E103552",
  "E103553",
  "E103559",
  "E10359",
  "E103591",
  "E103592",
  "E103593",
  "E103599",
  "E1036",
  "E1037X1",
  "E1037X2",
  "E1037X3",
  "E1037X9",
  "E1039",
  "E1040",
  "E1041",
  "E1042",
  "E1043",
  "E1044",
  "E1049",
  "E1051",
  "E1052",
  "E1059",
  "E10610",
  "E10618",
  "E10620",
  "E10621",
  "E10622",
  "E10628",
  "E10630",
  "E10638",
  "E10641",
  "E10649",
  "E1065",
  "E1069",
  "E108",
  "E1121",
  "E1122",
  "E1129",
  "E11311",
  "E11319",
  "E11321",
  "E113211",
  "E113212",
  "E113213",
  "E113219",
  "E11329",
  "E113291",
  "E113292",
  "E113293",
  "E113299",
  "E11331",
  "E113311",
  "E113312",
  "E113313",
  "E113319",
  "E11339",
  "E113391",
  "E113392",
  "E113393",
  "E113399",
  "E11341",
  "E113411",
  "E113412",
  "E113413",
  "E113419",
  "E11349",
  "E113491",
  "E113492",
  "E113493",
  "E113499",
  "E11351",
  "E113511",
  "E113512",
  "E113513",
  "E113519",
  "E113521",
  "E113522",
  "E113523",
  "E113529",
  "E113531",
  "E113532",
  "E113533",
  "E113539",
  "E113541",
  "E113542",
  "E113543",
  "E113549",
  "E113551",
  "E113552",
  "E113553",
  "E113559",
  "E11359",
  "E113591",
  "E113592",
  "E113593",
  "E113599",
  "E1136",
  "E1137X1",
  "E1137X2",
  "E1137X3",
  "E1137X9",
  "E1139",
  "E1140",
  "E1141",
  "E1142",
  "E1143",
  "E1144",
  "E1149",
  "E1151",
  "E1152",
  "E1159",
  "E11610",
  "E11618",
  "E11620",
  "E11621",
  "E11622",
  "E11628",
  "E11630",
  "E11638",
  "E11641",
  "E11649",
  "E1165",
  "E1169",
  "E118",
  "E1321",
  "E1322",
  "E1329",
  "E13311",
  "E13319",
  "E13321",
  "E133211",
  "E133212",
  "E133213",
  "E133219",
  "E13329",
  "E133291",
  "E133292",
  "E133293",
  "E133299",
  "E13331",
  "E133311",
  "E133312",
  "E133313",
  "E133319",
  "E13339",
  "E133391",
  "E133392",
  "E133393",
  "E133399",
  "E13341",
  "E133411",
  "E133412",
  "E133413",
  "E133419",
  "E13349",
  "E133491",
  "E133492",
  "E133493",
  "E133499",
  "E13351",
  "E133511",
  "E133512",
  "E133513",
  "E133519",
  "E133521",
  "E133522",
  "E133523",
  "E133529",
  "E133531",
  "E133532",
  "E133533",
  "E133539",
  "E133541",
  "E133542",
  "E133543",
  "E133549",
  "E133551",
  "E133552",
  "E133553",
  "E133559",
  "E13359",
  "E133591",
  "E133592",
  "E133593",
  "E133599",
  "E1336",
  "E1337X1",
  "E1337X2",
  "E1337X3",
  "E1337X9",
  "E1339",
  "E1340",
  "E1341",
  "E1342",
  "E1343",
  "E1344",
  "E1349",
  "E1351",
  "E1352",
  "E1359",
  "E13610",
  "E13618",
  "E13620",
  "E13621",
  "E13622",
  "E13628",
  "E13630",
  "E13638",
  "E13641",
  "E13649",
  "E1365",
  "E1369",
  "E138",
  "P702"))


a <- dm

b2$cm_dm <- with(b2, ifelse((i10_dx1 %in% a | i10_dx2 %in% a | i10_dx3 %in% a | i10_dx4 %in% a  | i10_dx5 %in% a | i10_dx6 %in% a | i10_dx7 %in% a | i10_dx8 %in% a |i10_dx9 %in% a| i10_dx10 %in% a | i10_dx11 %in% a | i10_dx12 %in% a | i10_dx13 %in% a | i10_dx14 %in% a | i10_dx15 %in% a ), "yes","no"))

b2 %>% count(cm_dm)



a <- dmcx

b2$cm_dmcx <- with(b2, ifelse((i10_dx1 %in% a | i10_dx2 %in% a | i10_dx3 %in% a | i10_dx4 %in% a  | i10_dx5 %in% a | i10_dx6 %in% a | i10_dx7 %in% a | i10_dx8 %in% a |i10_dx9 %in% a| i10_dx10 %in% a | i10_dx11 %in% a | i10_dx12 %in% a | i10_dx13 %in% a | i10_dx14 %in% a | i10_dx15 %in% a ), "yes","no"))

b2 %>% count(cm_dmcx)
```

```{r insulin use dm}


ins <- as.character(c('Z794',"E10"))

a <- ins

a <- liver

b2$iddm_n <- with(b2, ifelse((i10_dx1 %in% a | i10_dx2 %in% a | i10_dx3 %in% a | i10_dx4 %in% a  | i10_dx5 %in% a | i10_dx6 %in% a | i10_dx7 %in% a | i10_dx8 %in% a |i10_dx9 %in% a| i10_dx10 %in% a | i10_dx11 %in% a | i10_dx12 %in% a | i10_dx13 %in% a | i10_dx14 %in% a | i10_dx15 %in% a ), "yes","no"))

b2 %>% count(iddm_n)
```








```{r save df so that we do not have to above everytime}


write_csv(df,"H:/bita_nis/df/svydf.csv")

write_csv(b2, "H:/bita_nis/df/b2.csv")

```


##### till here done and df saved. Not actually knitted as markdown document yet...



Prepare table 1:

```{r}
# table 1 with data from 2007 - 2015 3 quarters:

tdf <- df %>% 
  select(age, bita, died, year.x, dispuniform, female, key, 
los, nis_stratum.x, trendwt, priormi, priorpci, 
chf, stemi, hospid, hosp_region, hosp_locteach, hosp_bedsize, 
race, cm_dm, cm_dmcx, cm_obese, cm_renlfail,discwt.x, cabg_n, media, inf, cpb, race_n, pay_n, bita_i,discwt.x, iddm_n, disp_n, iddm_n, iabp_n, pvad_n)

myvars <- c('age', 'bita', 'died', 'year.x', 'dispuniform', 'female',
'los',  'priormi', 'priorpci', 
'chf', 'stemi', 'hosp_region', 'hosp_locteach', 'hosp_bedsize',  'cm_dm', 'cm_dmcx', 'cm_obese', 'cm_renlfail', 'media', 'inf', 'cpb', 'race_n', 'pay_n', 'bita_i', 'disp_n', "iddm_n", "iabp_n", "pvad_n")



## Vector of categorical variables that need transformation

catVars <- c( 'bita', 'died', 'year.x', 'dispuniform', 'female',
  'priormi', 'priorpci', 'chf', 'stemi', 'hosp_region', 'hosp_locteach', 'hosp_bedsize',  'cm_dm', 'cm_dmcx', 'cm_obese', 'cm_renlfail', 'media', 'inf', 'cpb', 'race_n', 'pay_n', 'bita_i', "disp_n", "iddm_n", "iabp_n", "pvad_n")

t1d <- svydesign(data = tdf, strata = ~nis_stratum.x, id = ~hospid, weights = ~discwt.x, nest = TRUE)

tab1 <- svyCreateTableOne(vars = myvars , data = t1d, factorVars = catVars, includeNA = TRUE, smd = TRUE)


table1 <- print(tab1)
```


Doing some initial survey analysis to get basic #:

```{r some survey analysis}

# 1. type of admission: --- we need to work on the 2012 - 2015 to include type of admission

# type of admission NA 2012 onwards. Need to do seperate analysis for atype



```


Survey analysis by year to determine trends per year:

```{r}

# mean age per year for all CABG patients

svyby(~age, by = ~year.x, FUN = svymean, design = svydf)

# mean age for last period

age <- svyby(~age, by = ~year.x, FUN = svymean, design = b2)

#################################################################

##specific age groups

young <- svytable(~year.x + young, svydf)


y1 <- as.data.frame(young)

y2 <- y1 %>% spread(key = young,value = Freq )

y2 %>% tbl_df

y3 <- y2 %>% mutate(total = yes + no)

y4 <- y3 %>% mutate(perc = (yes/total)*100)

y4


## octogenarians

o <- svytable(~year.x + old, svydf)

o1 <- as.data.frame(o)

o2 <- o1 %>% spread(key = old,value = Freq )

o2 %>% tbl_df

o3 <- o2 %>% mutate(total = yes + no)

o4 <- o3 %>% mutate(perc = (yes/total)*100)

o4


# gender per year

f <- svytable(~year.x + female, svydf)

f1 <- as.data.frame(f)

tbl_df(f1)

f3 <- f1 %>% spread(key = female, value = Freq)

write_csv(f3, "H:/bita_nis/st/female_prop.csv")

# analysis for Q4

f_l <- svytable(~year.x + female, s15.4)

f_l




######################################################################
# BITA procedures per year of study

b <- svytable(~year.x + bita, svydf)

b1 <- as.data.frame(b)

tbl_df(b1)

b3 <- b1 %>% spread(key = bita, value = Freq)

b4 <- mutate(b3, total = no + yes)

b5 <- b4 %>% mutate(perc = ((yes/total)*100))

write_csv(b5, "H:/bita_nis/st/bita_prop.csv")

## Analysis for the last Q of 2015

b_l <- svytable(~year.x + bita, s15.4)

b_l

###################################################################


# prior mi per year

pmi <- svytable(~year.x + priormi, svydf)

as.data.frame(pmi)

pm1 <- pmi %>% tbl_df

pm2 <- pm1 %>% spread(key = priormi, value = n)

pm2

pm3 <- pm2 %>% mutate(total = yes + no)

pm4 <- pm3 %>% mutate(perc = ((yes/total)*100))

pm4


###########################################################################

# cardiogenic shock at presentation 

cy <- svytable(~year.x + cs, svydf)

as.data.frame(cy)

cy1 <- cy %>% tbl_df

cy2 <- cy1 %>% spread(key = cs, value = n)

cy2

cy3 <- cy2 %>% mutate(total = yes + no)

cy4 <- cy3 %>% mutate(perc = ((yes/total)*100))

cy4

###########################################################

# important to determine use of IABP or PVAD prior to surgery; so see if PVAD/IABP occurred same day or prior day.

# iabp

ia <- svytable(~year.x + iabp_n, svydf)

as.data.frame(ia)

ia1 <- ia %>% tbl_df

ia2 <- ia1 %>% spread(key = iabp_n, value = n)

ia2

ia3 <- ia2 %>% mutate(total = yes + no)

ia4 <- ia3 %>% mutate(perc = ((yes/total)*100))

ia4

############################################################
# pvad use during the study period

pv <-  svytable(~year.x + pvad_n, svydf)

as.data.frame(pv)

pv1 <- pv %>% tbl_df

pv2 <- pv1 %>% spread(key = pvad_n, value = n)

pv2

pv3 <- pv2 %>% mutate(total = yes + no)

pv4 <- pv3 %>% mutate(perc = ((yes/total)*100))

pv4

```




```{r doing some survey analysis}



# cabg done by region

svytable(~cabg_n + hosp_region, design = svydf)

svytable(~year.x + hosp_region, design = s15.4)

# gender distribution by region

loc_cabg <- svytable(~cabg_n + hosp_locteach, design = svydf)

loc_cabg


svytable(~year.x + hosp_locteach, design = s15.4)



# mean age per year for only patients undergoing BITA use

age_l <- svyby(~age, by = ~year.x, FUN = svymean, design = subset(svydf, bita == "yes"))

age
age_l

# gender per year for only patients with BITA used

fb <- svytable(~year.x + female, design = subset(svydf, bita == "yes"))

fb1 <- as.data.frame(fb)

write_csv(fb1, "H:/bita_nis/st/bita_female.csv")

fb_l <- svytable(~year.x + female, design = subset(s15.4, bita == "yes"))

fb_l


# hospital region

df %>% count(hosp_region)

region <- svytable(~hosp_region + bita, svydf)


r1 <- as.data.frame(region)

r3 <- r1 %>% spread(key = bita, value = Freq)

r3 %>% tbl_df

r4 <- mutate(r3, total = yes + no)

r5 <- mutate(r4, perc = ((yes/total)*100))

r5

write_csv(r5, "H:/bita_nis/st/region_bita.csv")


# teaching nature of hospital for BITA use 

df %>% count(hosp_locteach)


teach <- svytable(~hosp_locteach + bita, svydf)

t1 <- as.data.frame(teach)

t1 

t3 <- t1 %>% spread(key= bita, value = Freq)

as_tibble(t3)

t4 <- t3 %>% mutate(total = yes + no)

t5 <- t4 %>% mutate(perc = ((yes/total)*100))

t5

write_csv(t5, "H:/bita_nis/st/teach_bita.csv")
```




```{r}

# BITA use according to race

df %>% count(race_n)


r <- svytable(~race_n + bita, design = svydf)

r3 <- as.data.frame(r)

r4 <- r3 %>% spread(key = bita, value = Freq)

r5 <- r4 %>% mutate(total = yes + no)

r5 <- r5 %>% mutate(perc = ((yes/total)*100))

r5
```



```{r payment type and BITA use}

p <- svytable(~pay_n + bita, design = svydf)

svychisq(~pay_n + bita, design = svydf)

p3 <- as.data.frame(p)

p4 <- p3 %>% spread(key = bita, value = Freq)

p5 <- p4 %>% mutate(total = yes + no)

p5 <- p5 %>% mutate(perc = ((yes/total)*100))
p5


po <- svytable(~pay_n + bita, design = subset(svydf, age > 64))

svychisq(~pay_n + bita, design = subset(svydf, age > 64))

po3 <- as.data.frame(po)

po4 <- po3 %>% spread(key = bita, value = Freq)

po5 <- po4 %>% mutate(total = yes + no)

po5 <- po5 %>% mutate(perc = ((yes/total)*100))
po5

```




```{r}





mod <- svyglm( bita_i ~ age + factor(hosp_region) + cm_dm + cm_dmcx + cm_obese + race_n + pay_n, design = svydf)

or <- exp(coef(mod)) %>% tbl_df
ci <- exp(confint(mod)) %>% tbl_df

tab <- cbind(or,ci)

tab
```




```{r}

dm <- svytable(~cm_dm + bita, design = svydf)

dmonly <- svytable(~year.x + bita, design = subset(svydf, cm_dm == 1))

dm2 <- as.data.frame(dmonly) 

dmon <- dm2 %>% spread(key = bita, value = Freq)

dmo1 <- dmon %>% mutate(total = yes + no)

dmon2 <- dmo1 %>% mutate(perc = ((yes/total)*100))

dmon2


# nondiabetic patients 

nodm <- svytable(~year.x + bita, design = subset(svydf, cm_dm == 0  & cm_dmcx == 0))

nodm2 <- as.data.frame(nodm) 

nodm3 <- nodm2 %>% spread(key = bita, value = Freq)

nodm4 <- nodm3 %>% mutate(total = yes + no)

nodm5 <- nodm4 %>% mutate(perc = ((yes/total)*100))

nodm5
```




```{r for NE region}

b_ne <- svytable(~year.x + bita, design = subset(svydf, hosp_region == 1))

b_ne1 <- as.data.frame(b_ne)

tbl_df(b_ne1)

b_ne3 <- b_ne1 %>% spread(key = bita, value = Freq)

b_ne4 <- mutate(b_ne3, total = no + yes)

b_ne5 <- b_ne4 %>% mutate(perc = ((yes/total)*100))
```



```{r for young adults}

b_y <- svytable(~year.x + bita, design = subset(svydf, age < 51))

b_y1 <- as.data.frame(b_y)

tbl_df(b_y1)

b_y3 <- b_y1 %>% spread(key = bita, value = Freq)

b_y4 <- mutate(b_y3, total = no + yes)

b_y5 <- b_y4 %>% mutate(perc = ((yes/total)*100))

b_y5

b_yno <- svytable(~year.x + bita, design = subset(svydf, age < 51 & cm_dm == 0))



b_yno1 <- as.data.frame(b_yno)

tbl_df(b_yno1)

b_yno3 <- b_yno1 %>% spread(key = bita, value = Freq)

b_yno4 <- mutate(b_yno3, total = no + yes)

b_yno5 <- b_yno4 %>% mutate(perc = ((yes/total)*100))

b_yno5

```


Use of CPB during surgery:
 
```{r}
cpb_n <- svytable(~year.x + cpb, design = svydf)

cpb_n1 <- as.data.frame(cpb_n)

tbl_df(cpb_n1)

cpb_n3 <- cpb_n1 %>% spread(key = cpb, value = Freq)

cpb_n4 <- cpb_n3 %>% mutate(total = yes + no)

cpb_n5 <- cpb_n4 %>% mutate(perc = ((no/total)*100))

cpb_n5



reg_c <- svytable(~hosp_region + cpb, design = svydf)

as.data.frame(reg_c) 

reg_c1 <- reg_c %>% tbl_df

reg_c2 <- reg_c1 %>% spread(key = cpb, value = n)

reg_c2
```