## Overview
Last week we learned about distributions and their visualizations. To wrap up, we will also review boxplots, which are a useful way to inspect distributions, especially from a comparative perspective. 

Next we will we move on to explore the concept of correlation and the different ways of visualizing and exploring correlations in our data. We will take an intuitive approach for how to sensibly interpret correlations in your data. 

To begin, we will use data on airquality in New York City. You can visit: https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/airquality.html for more information about this data set. Note: Ozone concentrations refer to the quantity of ozone (O3) molecules in the air. 

*Motivating question: Are Heat Waves Associated with Lower Air Quality?*

- https://www.thoughtco.com/heat-waves-make-air-quality-worse-1204013

From *thoughtCo:* "Air quality decreases during times of hot temperatures because the heat and sunlight essentially cook the air along with all the chemical compounds lingering within it. This chemical soup combines with the nitrogen oxide emissions present in the air, creating a “smog” of ground-level ozone gas."

## Explore the data
```{r}
library(datasets)
library(ggplot2)

# together: check the data
data(airquality)
head(airquality)
summary(airquality)

# together: how would you check top 5 cases in data (this is very similar to head())
airquality[1:5,]

# together: check information about the distribution of a variable you are interested in 
# (hint: look into the mean, summary, max, etc)
summary(airquality)
mean(airquality$Temp)
max()
quantile(airquality$Temp)
summary(airquality$Ozone)
```

Notice that “Month” and “Day” are coded as numeric variables even though they are clearly categorical. This can be mended as follows, e.g.:

```{r}
# together: change month to a factor with labels
airquality$Month<-factor(airquality$Month,
                         labels=c("May", "Jun", "Jul","Aug", "Sep"))
```

From R documentation: The function "factor" is used to encode a vector as a factor (the terms ‘category’ and ‘enumerated type’ are also used for factors). This means that you can turn a different kind of variable into a categorical variable. Also, you can label the "levels" of this new variable, as shown above. 

## Basic boxplot 

R has a pretty easy way to use base R function for boxplots. We can examine the numeric data using this built-in function (Notice that missing data causes no problems to the boxplot function (similar to summary)):

```{r}
# together, let's use boxplot()
boxplot(airquality[,1:4])
```

The figure we just created is not always meaningful, as the variables may not be on comparable scales. The real power of box plots is to do comparisons of variables by sub-grouping. For example, we may be interested in comparing the fluctuations in temperature across months. To create boxplots of temperature data grouped by the factor “month”, we use the command: 

```{r}
# together, let's compare temperatures across months
boxplot(Temp ~ Month, data=airquality)
```

To learn more about the arguments used inside the boxplot() function check out the helpfile or this tutorial: http://www.sthda.com/english/wiki/box-plots-r-base-graphs

## Boxplots in ggplot

Now let's toy with the same information, but in ggplot. And, recall our motivating question: *Motivating question: Are Heat Waves Associated with Lower Air Quality?* In order to begin we will create a plot and tell `ggplot` that `airquality` is our data, and specify that our x-axis plots the `Month` variable and our y-axis plots the `Ozone` variable. We then instruct ggplot to graph this as a boxplot by adding the geom_boxplot() option. 

```{r}
# together, let's create our first boxplot using ggplot
plot1 <- ggplot(airquality, aes(x=Month, y=Ozone)) + 
      geom_boxplot()
```

Why did we get this warning message? If we were to examine our variables using the summary() function, we might have a pretty clear answer:

```{r include=FALSE}
summary(airquality$Ozone)
```

We can see here that the number of NAs = 37, which is the same number given in the warning. If we remove NAs (we will talk about missing data later in the class), we will no longer get this warning. Below is one way you could subset NA's here. 

Here is a nice bit of info on pipes: https://www.datacamp.com/community/tutorials/pipe-r-tutorial

```{r include=FALSE}
# together: library dplyr, subset the data for NAs, plot again
library(dplyr)
airSub = airquality %>% subset(!is.na(Ozone)) 

plot1 <- ggplot(airSub, aes(x = Month, y = Ozone)) +
        geom_boxplot()
plot1
```

Let's fix up those axes a bit. Let’s make the y-axis ticks appear at every 25 units rather than 50 using the breaks = seq(0, 175, 25) argument in scale_y_continuous. (The seq function is a base R function that indicates the start and endpoints and the units to increment by respectively. See help(seq) for more information.) We ensure that the y-axis begins and ends where we want by also adding the argument limits = c(0, 175) to scale_y_continuous.

```{r}
plot1 <- plot1 + scale_x_discrete(name = "Month") +
       scale_y_continuous(name = "Mean ozone in\nparts per billion",
                              breaks = seq(0, 175, 25),
                              limits=c(0, 175)) + 
        ggtitle("Boxplot of mean ozone by month")
plot1
```

As mentioned earlier, the visual power of a boxplot is often found in its ability to quickly relay differences across variables (or in the above example, the same variable but for different time periods). To go farther into comparisons, we first need to do a little data wrangling. In order to make the graphs a bit clearer, we’ve kept only months “July”, “Aug” and “Sep” in a new dataset called `airquality_trimmed`. We’ve also used the mean as way to generate categories (high and low).

Helpful info:

- Recall: https://www.statmethods.net/management/operators.html
- "|" logical operator in R refers to "OR" statements. 

_In pairs:_ For the code below, get with a partner and walk through the code to try to understand what it does. Write a comment for each line within the R chunk.

# Boxplot across groups
```{r}
# review the code below:
# trim the data for the months you want
airquality_trimmed <- airSub[which(airSub$Month == "Jul" |
                                       airSub$Month == "Aug" |
                                       airSub$Month == "Sep"), ]

# create a new variable based on mean temperature and "low high" groupings
airquality_trimmed$Temp.f <- factor(ifelse(airquality_trimmed$Temp > mean(airquality_trimmed$Temp), 1, 0),labels = c("Low temp", "High temp"))

# now create a boxplot that compared low and high averages. What is new here? (what does "fill" do?)
# challenge yourself to write a comment to explain each line of code
finalPlot <- ggplot(airquality_trimmed, aes(x = Month, y = Ozone, fill = Temp.f)) +
        geom_boxplot(alpha=0.7) +
        scale_y_continuous(name = "Mean ozone in\nparts per billion",
                           breaks = seq(0, 175, 25),
                           limits=c(0, 175)) +
        scale_x_discrete(name = "Month") +
        ggtitle("Boxplot of mean ozone by month") +
        theme_bw() +
        theme(plot.title = element_text(size = 14, face = "bold"),
              text = element_text(size = 12),
              axis.title = element_text(face="bold"),
              axis.text.x=element_text(size = 11)) +
        scale_fill_brewer(palette = "Accent")
finalPlot
```

There are a lot of things you might want to toy with here, such as color, legend position, or legend label. Spend time outside of class with the above code chunk and get familiar with the different ways we could have examined this data. 

Summary: our motivating question was *Are Heat Waves Associated with Lower Air Quality?* Is the evidence suggestive here? 

# Correlations

A *correlation coefficient*, often referred to as "r", is a number between -1 and 1 that captures the strength of the linear relationship between two numeric variables. For instance, say you’ve asked 30 people about their resting heart rate and how often they work out. 

- If all 30 data points fall perfectly on an increasing line, then the correlation between these two variables will be r = 1.
- If the relationship is a decreasing one, then r will lie between 0 and -1.

We are going to explore a couple of simplified examples and then for your homework you will dig deeper into data that you have not yet worked with. 

## Correlation basics 

Pearson correlation (r), which measures:

- a linear dependence between two continuous variables (x and y)
- or the strength of the linear relationship between two variables. 
- This approach is so widely used that when most people refer to correlation they are referring to the Pearson correlation approach.

## R functions for Correlation

Correlation coefficients can be computed in R by using the functions `cor()` and `cor.test()`. We will focus on `cor()` for now (though you are welcome to explore `cor.test` on your own!)

```{r}
?cor
```

The simplified format of this function is `cor(x, y, method = c("pearson", "kendall", "spearman"))`. Where:

  - x, y: numeric vectors with the same length
  - method: correlation method (we will use "pearson")
  - if the data contain missing values, the following R code can be used to handle missing values by case-wise deletion:

```{r}
#example code:
#cor(x, y,  method = "pearson", use = "complete.obs")
```

Let's try this out from the data you used for your QUICKFIRE competition last time.

```{r}
library(readr)
metro_nash_schools <- read_csv("/Users/dorffc/ProjectsGit/instructor-copy-eda/classes/week4/metro-nash-schools.csv")
```

Let's specifically examine a relationship that many students tried to study with the data on their own: the relationship between Disability and Economic disadvantages in the nashville metro schools data. 

```{r}
library(tidyverse)
head(metro_nash_schools)

# let's get two pesky things out of the way 
# rename the data to something easier
# get rid of spaces in variable names, which R hates
# together:
schools<-metro_nash_schools
names(schools) <- str_replace_all(names(schools), c(" " = "."))

# now let's try a quick scatterplot
ggplot(schools, aes(x=Economically.Disadvantaged, y=Disability))+
    geom_point()
```

Before we 'beautify' the plot or increase its interpretability, we might want to know what the correlation coefficient is! 

```{r include=FALSE}
cor(schools$Economically.Disadvantaged, schools$Disability)

#oops! NA return. Why?
summary(schools)

# We have an NA in the Economically Disadvantaged variable we need to remove
# let's go ahead and do that so that we know it is not showing up in any other analyses
# or plots in the data
schools<-schools[complete.cases(schools[,"Economically.Disadvantaged"]),]

# NOTE: we could also use the argument within the function:
# cor(schools$Economically.Disadvantaged, schools$Disability, use="complete.obs")
```

Now we can examine the correlation statistic again. If you would like to read more about how to calculate correlation, you will learn more in your stats class! Or you can review online materials such as: https://www.toppr.com/guides/business-mathematics-and-statistics/correlation-and-regression/karl-pearsons-coefficient-correlation/

```{r}
# correlation of e# of conomically disadvantaged students and students with disability
cor(schools$Economically.Disadvantaged, schools$Disability)
```

The correlation coefficient is about 0.68 which reflects a high, positive correlation. Let's plot this relationship below.

```{r}
ggplot(schools, aes(x = Economically.Disadvantaged, y = Disability)) +
  geom_point(na.rm = TRUE) +
  labs(title="Number of students with Disability \n vs. Number of Economically Disadvantaged Students in Nashville Schools", y="Disability", x="Number of Economically Disadvantaged Students") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5))
```

Now, can we add the correlation coefficient to the graph:

```{r}
ggplot(schools, aes(x = Economically.Disadvantaged, y = Disability)) +
  geom_point(na.rm = TRUE) +
  labs(title="Number of students with Disability \n vs. Number of Economically Disadvantaged Students in Nashville Schools \n (Correlation=.68)", y="Disability", x="Number of Economically Disadvantaged Students") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5))
```

That was a cheap an easy way to do it. How else, or where else, might we add it?

```{r}
# on your own, edit the graph above
# specifically, plot the correlation coefficient onto the graph itself (not in a title)

ggplot(schools, aes(x = Economically.Disadvantaged, y = Disability)) +
  geom_point(na.rm = TRUE) +
  labs(title="Number of students with Disability \n vs. Number of Economically Disadvantaged Students in Nashville Schools", y="Disability", x="Number of Economically Disadvantaged Students") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) + 
  geom_text(x=750, y=350, label="Correlation=.68", size=5)


```

Finally, for fun and in pairs, we can explore a similar trend to Yasi's plot in class! Below can you add an aesthetic mapping to map `Limited English Proficiency` onto the geom points? Edit your scale color gradient to be something other than the default.
 
```{r}
ggplot(schools, aes(x = Economically.Disadvantaged, y = Disability)) +
  geom_point(na.rm = TRUE, aes(col=Limited.English.Proficiency)) +
  labs(title="Number of students with Disability \n vs. Number of Economically Disadvantaged Students in Nashville Schools", y="Disability", x="Number of Economically Disadvantaged Students") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) + 
  geom_text(x=750, y=350, label="Correlation=.68", size=5) +
  scale_color_gradient(low="yellow", high="brown")

```

## On your own, some more! 

Here is a different set of data for you to examine on your own with your partner:

```{r}
gapminder <- read.csv("https://raw.githubusercontent.com/swcarpentry/r-novice-gapminder/gh-pages/_episodes_rmd/data/gapminder-FiveYearData.csv")

head(gapminder)

# code on your own:
ggplot(gapminder, aes(x = continent, y = lifeExp, fill = continent)) +
  geom_boxplot()

```

*HW:* Your homework will be posted online

## References
- http://t-redactyl.io/blog/2016/04/creating-plots-in-r-using-ggplot2-part-10-boxplots.html
- http://people.bu.edu/aimcinto/720/lec2.pdf
- http://www.rebeccabarter.com/blog/2017-11-17-ggplot2_tutorial/
- https://janhove.github.io/teaching/2016/11/21/what-correlations-look-like
