---
title: "Plotting data with ggplot2"
author: "Julian Faraway"
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  html_document:
    keep_md: yes
---

```{r echo=FALSE}
library(knitr)
opts_chunk$set(comment=NA, warning=FALSE, fig.path="Figs/")
```

The `ggplot2` package is part of the `tidyverse`. It is useful for plotting data.

A reference to ggplot commands can be found at
<https://ggplot2.tidyverse.org/reference/index.html>
or a collection of examples at
<http://www.cookbook-r.com/Graphs/>

We need to start by loading the `tidyverse` and the data.

```{r}
library(tidyverse)
load("uliral.rda")
```

If you continue immediately from the previous exercise, the data would already be loaded. But is sensible to make
your scripts work individually.


# One dimensional plots


Let's start by looking at plots of one variable at a time. Here is an example:


```{r}
ggplot(uliral, aes(x=Gender)) + geom_bar()
```

The first part of the command specifies the data as the first argument. This
must be a data frame (not vectors or matrices). Our data frame is `uliral`. The
*aesthetic* part of the command is give by `aes()`. This specifies how the
different variables should be represented in the plot - we might have an *x*
direction (horizontal) and a *y* direction. We might also have colors,
linetypes, different plot symbols and more. In this example, there is only one
aesthetic variable which goes on the horizontal axis.

The second part of the command gives the *geometry* of the plot. This tells us
what kind of plot we want - for example a histogram or a scatterplot. In this
example, we choose a bar plot.

In the following example, we add a third element to the plot. We can specify the label for the x-axis:

```{r}
ggplot(uliral, aes(x=City)) +  geom_bar() + xlab("Hospital City")
```

We can see that the labels for the cities overlap. One solution for this is to change the orientation of the plot.

```{r}
ggplot(uliral, aes(x=City)) +  geom_bar() + xlab("Hospital City") + coord_flip()
```

Another option is:

```{r}
ggplot(uliral, aes(x=City)) +  geom_bar() + xlab("Hospital City") + theme(axis.text.x = element_text(angle=30, hjust=1))
```

The `theme` component modifies the appearance of the plot. In this case, it changes the labelling. If you type `?theme` at the prompt, you will get a summary of the many options.


For more continuous variables, a histogram is the standard choice:

```{r}
ggplot(uliral, aes(x=Age)) + geom_histogram() 
```

Sometimes you should adjust the number or width of the bins:

```{r}
ggplot(uliral, aes(x=Age)) + geom_histogram(bins=100)
```

You can choose a different geometry provided you have the same inputs. Here is a density estimate:

```{r}
ggplot(uliral, aes(x=Age)) + geom_density()
```

## Exercises

1. Make a plot of `Education`.

2. Make a plot of `TreatDays`.

# Two variables, both quantitative

Before we go on, let's filter the data to only Ulaanbaatar. This reduces the size of the data and makes the plots render faster:

```{r}
uliral %>% filter(City == "Ulaanbaatar") -> uliral
```


The appropriate choice of plot depends on the variable types. For two
quantitative variables, a scatterplot is the standard choice:

```{r}
ggplot(uliral,aes(x=TreatDays, y=Cost)) + geom_point()
```

In this case, we see there are some high cost patients skewing the plot. We should
investigate these cases more closely.

When a variable is skewed with many small values and few larger values, a log
scale may be a better choice:

```{r}
ggplot(uliral,aes(TreatDays, Cost)) + geom_point() + scale_y_log10()
```

For many patients, the cost is listed as zero. Since the log(0) = - Infinity, this causes a warning
message. Note that these points are displayed at the bottom of the plot.

We can also add a smooth curve fit to the plot:

```{r}
ggplot(uliral,aes(TreatDays, Cost)) + geom_point() + scale_y_log10() + stat_smooth()
```

The curve is computed using a robust local average fitting method. The grey shaded
area represents a 95% confidence band. We can see that cost is related to the number of treatment days. 

Let's try another pair of quantitative variables:

```{r}
p = ggplot(uliral,aes(Age, TreatDays))
p + geom_point()
```

I have saved part of the plot as `p` and then displayed it with the points. We'll see why soon.

The plot has some problems. Both variables are discrete so we cannot
distinguish all the points. One solution to this problem is to *jitter* the
points by adding small random perturbations:

```{r}
p + geom_jitter(width=0.2,alpha=0.25,size=0.5)
```

The `width=0.2` specifies that the random perturbation should be in the horizontal
direction only with 0.2 giving the amount (some experimentation with this may be
needed). I have also added some *alpha transparency*. This makes the points partially
transparent which makes it easier to distinguish them. I also decreased the plot size
with `size=0.5`.

Notice how saving `p` was useful - we did not need to repeat the specification
of much of the plot.

## Exercises

Make a plot of `Age` against `Cost`.  Think how to present this plot well. What is the relationship between these two variables?


# Two variables, one quantitative, one categorical

Let's try to plot `Education` (categorical) against `Cost` (quantitative). 

```{r}
p = ggplot(uliral, aes(Education, Cost))
p + geom_point()
```

As with the previous plot, we have problems in distinguishing the points. We make a similar modification:

```{r}
p + geom_jitter(width=0.1,alpha=0.25, size=0.5)
```

A popular alternative is the *boxplot*. This displays a five number (min, first quartile,
median, third quartile, max) summary of each level of the categorical variable. More
extreme cases are plotted individually. Here are the full details on [boxplots](https://en.wikipedia.org/wiki/Box_plot).

```{r}
p + geom_boxplot()
```

The boxplot makes it easy to compare the distributions within each level. In this case, the cost is very skewed so a log-transformation is helpful.

```{r}
p + geom_boxplot() + scale_y_log10()
```


We can see that the cost varies somewhat among the education levels. The
roughly equally length of the boxes indicates that the variation in the three
groups is similar.

The *violin* plot is an elaboration of a boxplot. It shows the whole distribution
and not just the five number summary:

```{r}
p + geom_violin()+ scale_y_log10()
```

We can see the distributions.

These plots have been excluding those patients who have zero cost. We should analyse how many there are:

```{r}
uliral %>% mutate(costzero = (Cost == 0)) %>% 
  group_by(Education) %>% 
  summarise(n=n(),zc=sum(costzero), fraction=zc/n) -> edcost
edcost
```

We can see the number in each education level, the number with zero cost and then represent this as a fraction.

We can plot this information. We get rid of the NA cases before making the barplot. We also show
how to order the levels by frequency.

```{r}
edcost %>% 
  filter(!is.na(Education)) %>% 
  ggplot(aes(reorder(Education,-fraction),fraction)) + 
  geom_col() +
  xlab("Education Level") + ylab("Fraction with zero cost")
```


## Exercises

1. Make one or more plots that show how `TreatDays` varies with `Education`.

2. In the previous section, we showed how `Cost` varied with `Age`.
Instead of a scatter plot, make a boxplot that shows this relationship. You will need to specify
the aesthetic like this:
```
aes(Age, Cost, group=Age)
```
What is the advantage of this plot?


# Two categorical variables

There are several ways to plot two categorical variables. Here
we use the fill color to show the count of the second variable:

```{r}
ggplot(uliral, aes(Education, fill=Gender)) + geom_bar()
```

We can display the second variable using side-by-side bars:

```{r}
ggplot(uliral, aes(Education, fill=Gender)) + geom_bar(position = "dodge2")
```

This makes it easier to see the slight imbalance.

## Exercise

For this question, we make two new variables expressing the day of the week that the patient was
admitted and departed.

```{r}
library(lubridate) 
uliral %>% mutate(awday = wday(Admitted, label=TRUE),
                  dwday = wday(Departure, label=TRUE)) -> uliral
uliral %>% group_by(awday) %>% summarise(n())
```

Make plots showing the relationship between admission day and departure day.

# More variables

There are a number of ways to display more than 2 variables in a single plot. Consider
the variables `Age`, `TreatDays` and `Color` - two quantitative and one categorical. We
can display all three in the same plot:

```{r}
ggplot(uliral, aes(x=Age,y=TreatDays,color=Gender)) + geom_jitter(alpha=0.25,size=0.5)
```

We have three aesthetics, `TreatDays` in the y-direction, `Age` in x-direction and `Gender` displayed using color.

Some people are color-blind and not all printers have color so we might also want to distinguish
the `area` using the shape of the plotting character:

```{r}
uliral %>% 
  filter(Hospital == 2026) %>% 
  ggplot(aes(x=Age,y=TreatDays,color=Gender,shape=Gender)) + 
  geom_jitter(size=2)
```

In this plot, we have filtered to a single hospital so that it is easier to distinguish the points.

With so many points, it can be difficult to distinguish trends in the two groups. We can
add smooth curves to help us detect the trends:

```{r}
uliral %>% 
  filter(Hospital == 2026) %>% 
  ggplot(aes(x=Age,y=TreatDays,color=Gender,shape=Gender)) + 
  geom_jitter(size=2) +
  geom_smooth()
```

We can also display separate plots for each level of the categorical variable `area`. This
idea is called *facetting*.

```{r}
ggplot(uliral, aes(x=Age,y=TreatDays)) + 
  geom_jitter(alpha=0.25,size=0.5) +
  facet_grid(~ Gender)
```

We might want to add the smooth curve fits again:

```{r}
ggplot(uliral, aes(x=Age,y=TreatDays)) + 
  geom_jitter(alpha=0.25,size=0.5) +
  facet_grid(~ Gender) +
  geom_smooth()
```

The facetting idea can be applied whenever we have a categorical variable with a small number of levels. Let's redo
one of our previous examples:

```{r}
ggplot(uliral, aes(Education)) + geom_bar() + facet_grid(~ Gender)
```

When the categorical variable has more levels, it is better to wrap than grid:

```{r}
ggplot(uliral, aes(Gender)) + geom_bar() + facet_wrap(~ Education)
```

## Exercises

1. Make plots to show how the relationship between `Age` and `Cost` depends on the day of the week
when the patient was admitted.

2. Show how the age distribution varies by `Age` and `awday`.

