[This article was first published on R Views, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't. Vidisha Vachharajani works in the EdTech industry, where she enjoys developing data-driven strategy solutions for learners. She has been an R user for over 15 years.
As a data professional, I have enjoyed learning and using multiple tools for my workflows. For me, everything used to begin and end with R. Today, SQL is a must-know. Not being able to pull your own custom tables from a warehouse can make things tricky. Then there is tidyverse, the master collection of packages for data science & analytics. As an OG R user, I cannot envision data work without tidyverse.
In this first part of a 2-part article, I want to demonstrate how a data analyst can use one OR the other for the initial stages of data exploration, and then double down on tidyverse, leveraging ggplot2 for a deeper exploration. By no means does this preclude the extensive use of SQL for data wrangling. Rather, this post showcases the wonders of tidyverse (a collection of R packages designed for data science, sharing an underlying design philosophy, grammar, and data structures) and specifically, ggplot2 (the language of elegant graphics) for a SQL user’s benefit.
tidyverse language in tandem, I will split it up into 5 parts, and we will assume that the data is actually available to us in these 5 different pieces, rather than as the whole, cleaned data, since this is typically the case in real life.I will skip the portion about dbplyr, referring readers to the hyperlinked article that will show you how to actually pull data from a remote database using tidyverse’s dbplyr. Typically, this is done using SQL, butdbplyr allows you to do this within R. Rather, I will focus on the initial stages of data exploration, using both SQL and tidyverse for the same output, while extending the tidyverse portion to include ggplot2 visualization examples, using different plot types for each use case. Note that in each case, you can use SQL first, and then use the SQL output as an input for the ggplot2 visualization.
```
``
3. Early explorationsLet’s begin using SQL andtidyverseto answer some initial questions related to the dataset. The primary hypothesis for this data is the **impact of HbA1c measurement on readmission rates**, where “readmission” is our response. We will also answer a number of other questions along the way to understand the data better, usingggplot2` when we can.
3.1 Look at the data3.1.1 Get some countsLet’s take a look at medications and get a sample size for it, first using SQL and then R.
sqldf('SELECT * FROM meds where 1=0') # SQL see col names## [1] uid metformin repaglinide ## [4] nateglinide chlorpropamide glimepiride ## [7] acetohexamide glipizide glyburide ## [10] tolbutamide pioglitazone rosiglitazone ## [13] acarbose miglitol troglitazone ## [16] tolazamide examide citoglipton ## [19] insulin glyburide-metformin glipizide-metformin ## [22] glimepiride-pioglitazone metformin-rosiglitazone metformin-pioglitazone ## [25] change diabetesMed ## <0 rows> (or 0-length row.names)sqldf('SELECT uid, metformin, repaglinide, nateglinide, chlorpropamide FROM meds LIMIT 5') # SQL## uid metformin repaglinide nateglinide chlorpropamide## 1 2278392-8222157 No No No No## 2 149190-55629189 No No No No## 3 64410-86047875 No No No No## 4 500364-82442376 No No No No## 5 16680-42519267 No No No Nohead(meds, n=5) # dplyr## # A tibble: 5 × 26## uid metfo…¹ repag…² nateg…³ chlor…⁴ glime…⁵ aceto…⁶ glipi…⁷ glybu…⁸ tolbu…⁹## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> ## 1 22783… No No No No No No No No No ## 2 14919… No No No No No No No No No ## 3 64410… No No No No No No Steady No No ## 4 50036… No No No No No No No No No ## 5 16680… No No No No No No Steady No No ## # … with 16 more variables: pioglitazone <chr>, rosiglitazone <chr>,## # acarbose <chr>, miglitol <chr>, troglitazone <chr>, tolazamide <chr>,## # examide <chr>, citoglipton <chr>, insulin <chr>,## # `glyburide-metformin` <chr>, `glipizide-metformin` <chr>,## # `glimepiride-pioglitazone` <chr>, `metformin-rosiglitazone` <chr>,## # `metformin-pioglitazone` <chr>, change <chr>, diabetesMed <chr>, and## # abbreviated variable names ¹metformin, ²repaglinide, ³nateglinide, …sqldf('SELECT COUNT(uid) FROM meds') # SQL## COUNT(uid)## 1 101766nrow(meds) # R## [1] 101766
How many patients with a diabetes diagnosis, vs respiratory, circulatory, etc.?
sqldf('SELECT primary\_diag, COUNT(*) FROM results GROUP BY primary\_diag') # SQL## primary\_diag COUNT(*)## 1 circulatory 30437## 2 diabetes 8757## 3 other 48149## 4 respiratory 14423results %>% group\_by(primary\_diag) %>% count(primary\_diag) # R## # A tibble: 4 × 2## # Groups: primary\_diag [4]## primary\_diag n## <chr> <int>## 1 circulatory 30437## 2 diabetes 8757## 3 other 48149## 4 respiratory 14423
How many women came in through an emergency admission type?
sqldf('SELECT gender, admission\_type\_id, COUNT(*) AS n FROM dem LEFT JOIN visits USING(uid) WHERE admission\_type\_id=1 GROUP BY gender') # SQL## gender admission\_type\_id n## 1 Female 1 29448## 2 Male 1 24540## 3 Unknown/Invalid 1 2visits %>% left\_join(dem, by=join\_by(uid)) %>% subset(admission\_type\_id==1) %>% count(gender, admission\_type\_id) # dplyr## # A tibble: 3 × 3## gender admission\_type\_id n## <chr> <dbl> <int>## 1 Female 1 29448## 2 Male 1 24540## 3 Unknown/Invalid 1 2
3.1.2 A mosaic plotInstead of extracting counts manually, let’s use a mosaic plot to get a sense of how 2 count variables are distributed relative to each other. In this case, age and admission type. This plot sheds light into data availability and asymmetric distributions. For example, here, we see that most patients come from emergency, urgent care, or as an elective, and that there is missing or “not available” admission type data. It is important to retain these 2 categories separately, since they mean different things. Note that in the ggplot parameters, I have not yet introduced axes label cleanup, etc.
p0 <- dem %>% left\_join(visits, by=join\_by(uid)) %>% mutate(admission\_type=ifelse(admission\_type\_id==1, "1:Emergency", ifelse(admission\_type\_id==2, "2:Urgent", ifelse(admission\_type\_id==3, "3:Elective", ifelse(admission\_type\_id==4, "4:Newborn", ifelse(admission\_type\_id==5, "5:Not Available", ifelse(admission\_type\_id==6, "6:NULL", ifelse(admission\_type\_id==7, "7:Trauma Center", "8:Not Mapped")))))))) %>% group\_by(admission\_type, age) %>% summarise(n=n()) %>% mutate(freq = n / sum(n)) ggplot(p0, aes(x=age, y=admission\_type)) + geom\_tile(aes(fill=n)) + scale\_fill\_gradient(low="white", high="blue")
3.1.3 A simple joinLet’s join all 5 datasets and look at it. Note that in SQL, in order to look only at the first few columns, we need to know the column names, which is what we first do here.
```
``` 3.2 Explore the response: readmissions3.2.1 Lab proceduresLet’s start with the simplest question – for the primary response variable, “readmitted”, how many lab procedures were done by each category of the response? Note here that “number of lab procedures” is one of a handful of continuous design covariate – rest of the ~45 covariates are all categorical/discrete.
```
``
Since the above doesn’t really tell us much, other than actual counts, proportions by response categories, let’s useggplot2to explore the distribution of “number of lab procedures”, using a barplot/histogram approach, with “readmitted” as thefill` element. This helps us get a better picture of their relationship; we see here how, for a strikingly normally distributed “number of lab procedures” (other than 1 outlier), on average, the higher the volume of procedures, the more the proportion of readmitted.
```
``
Let’s also do this usingggplot’s beautiful density plots. It is a slightly different type of visual, and tells us how the distribution of X shifts left or right by the response orfill`.
ggplot(p1, aes(num\_lab\_procedures)) + geom\_density(aes(fill=factor(readmitted)), alpha=0.8) + labs(x="Number of lab procedures")
3.2.2 DemographicsNext, we ask how readmissions differ across age groups and gender. Let’s also plot this to understand the output better. We first use a population pyramid approach to get the counts and then barplot the proportions to get a better understanding of the variance in readmissions across these groups.
```
```
p22 <- dem %>% left\_join(y, by=join\_by(uid)) %>% group\_by(gender, age, readmitted) %>% summarise(n=n()) %>% mutate(freq = n / sum(n)) %>% subset(gender=="Male"|gender=="Female")ggplot(data=p22, aes(x=age, y=freq, fill=readmitted)) + geom\_col() + facet\_wrap(~ gender) + labs(y="proportions") + geom\_text(aes(label = paste0(round(freq, 4) * 100, "%")), position = position\_stack(vjust = 0.5), size=2.5, angle=90) + theme(axis.text.x = element\_text(angle=90, vjust=.5, hjust=1))
The population pyramid is an intriguing plot type, and already tells us that for most age groups, more women are readmitted. But this could be solely because there are more women than men in the sample. However, from the proportion barchart, we see here that proportion of readmitted women is greater than men, particularly for the 20-30 age group.
3.2.3 Patient diagnosesFinally, how are readmission rates distributed by patient and patient care features. For example, how is it distributed by patient primary diagnosis? In the final section of this post, we will leverage ggplot2’s visualization power to triangulate patient diagnoses with the key covariate and the response. Like in the previous section, we use proportions, adding the relevant labels to more easily infer that we see higher readmission rates for a diabetes diagnosis.
We change around quite a few of the plotting parameters in ggplot2 to make it look much more eye-catching.
```
``
3.2.4 HbA1c measurementOne of the key questions this dataset seeks to answer is the *impact of the A1C test (decision to test) on readmission rates*, in the presence of covariates (especially the primary diagnosis). Output in its raw form (i.e. untransformed) doesn’t always give us the answer clearly. To get around this, we will useCASE WHENin SQL andmutateintidyverse`.
Let’s plot this in 2 ways – a barplot with labels, and a spineplot. The latter allows us to see the “weight” of the underlying categories.
```
```
```
``
We observe a lower readmission rate (<30 days) when there is an A1C measurement taken, vs when it is not measured at all. In the 2nd/spineplot, we see this without actually calculating the percentages, while also inferring that number of patients not measured is much higher than those measured. We do however, manually add in the percentages to the spineplot to get a more complete picture on the relationship between HbA1c measurement and readmission rates.These are key findings which we will explore in greater detail, usingtidyverseandggplot2` more extensively, in the next part of this blog series, including cutting these plots across multiple covariates to explore how HbA1c affects readmissions in the presence of other patient groupings. Stay tuned!
To leave a comment for the author, please follow the link and comment on their blog: R Views.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.Continue reading: A data analyst workflow, part 1: SQL & tidyverse