---
author: "Sean Laverty"
title: "Class 8 notes"
date: "`r Sys.Date()`"
format: html
editor: source
---

## Experiments in `ggplot`

Load the package.

```{r}
#| label: packages
library(ggplot2)
```

## iNaturalist data
This is a great source of data, with a somewhat intricate query process.  This is just a subset of observations I pulled from a previous dataset for our use. Even at that, there is a _lot_ to consider.

```{r}
#| label: read-inat-data
dat <- read.delim("./data/observations-526962.csv", header = TRUE, sep = ',')
head(dat, n = 2)
```

## A subset of birds
Below we've subset the data to look at two randomly chosen bird species and their pattern of occurrence in Oklahoma County, OK. In reality it may be of greater interest to choose species after browsing a table of occurrences.

```{r}
#| label: inat-tab
tab <- table(dat$common_name)
tab <- tab[order(tab, decreasing = TRUE)]
c(head(tab), tail(tab))
```

Interestingly, combining `head()` and `tail()` we can determine commonly and less commonly encountered species. By chance alone we chose the Northern Cardinal ($n_\text{obs} = 348$) and the Canada Goose ($n_\text{obs} = 373$), two of top four species listed in the data, each representing about 3\\% of total observations.  Using commands such as `tab[tab >10]` shows the species with more than ten observations (131 total species).

```{r}
#| label: tabled-counts
tab2 <- table(tab)
tab2 <- tab2[order(tab2, decreasing = TRUE)]
tab2
```
The previous table is useful to understand overall patterns.  Of all recordings, 62 different species were recorded only once.  That is 20\\% of the data being represented by at most a single observation. 

With that understanding, we now color the points of two familiar species, which happen to be among the most common, with a color generated from their `common_name` entry. Using the transparency provided by `alpha = 0.5` ensures that you can get a sense of overlap in locations where many observations are made.

```{r}
#| label: fig-cardinal-goose
#| fig.width: 5
#| fig.height: 5
#| fig.cap: "A plot of iNaturalist observations of Norther Cardinal and Canada Goose birds."
p <- ggplot(data = subset(dat, subset = common_name %in% c("Northern Cardinal", "Canada Goose")) ) + geom_point(aes(latitude, longitude, color = common_name), alpha = 0.5)
p
```

Using `labs()` we can quickly stylize the axes labels and provide a title and subtitle. 

```{r}
#| label: fig-cardinal-goose-labs
#| fig.width: 5
#| fig.height: 5
#| fig.cap: "A plot of iNaturalist observations of Norther Cardinal and Canada Goose birds with improved text annotation."
p <- p + labs(x = "Latitude", y = "Longitude", 
              subtitle = "iNat observations", title = "Birds in Oklahoma County, OK, USA")
p
```

As a proof of concept, we show that we can label points by the text of their `common_name` entries.  The code for `geom_label()` "works", but this is _very slow_ to render -- it requires drawing a label "box" for each point and a printing its corresponding label.  These sort of "screen output" tasks can be very slow.

```{r}
#| label: fig-cardinal-goose-labs2
#| fig.width: 5
#| fig.height: 5
#| fig.cap: "A plot of iNaturalist observations of Norther Cardinal and Canada Goose birds."
p <- p + geom_text(mapping = aes(x = latitude, y = longitude, label = common_name), hjust = 0) 
# + geom_label(mapping = aes(x = latitude, y = longitude, label = common_name))
p
```

## Olympic athlete data

Just like labeling all spatial locations of bird observations is possible, but not that useful -- at least for the most common of birds, giving colored symbols for every Olympic sport is a visually-challenging task.  Applying additional knowledge, it is reasonable to group sports by some meaningful characteristics. Below is a sample approach for grouping Olymipic sports by some shared characteristic. Of course if a classification scheme existed already, we could read it in and merge it. But, making and recording decisions is a valuable part of exploration.

```{r}
#| label: data-olympics
dat <- read.delim("./data/athletes.csv", header = TRUE, sep = ',')
head(dat)
dat$group <- dat$sport
dat[dat$group %in% c("aquatics", "sailing", "rowing"), "group"] <- "water"
#dat[dat$group %in% c("aquatics", "sailing", "rowing"), "group"] <- "field"

```

Bar graphs, which aren't great to begin with, are rather hard to produce as we might like, in `ggplot`.  There is interest in representing the accumulation of different medals by `nationality` (or possibly by `sex` and `sport`).  I find this not the easiest to accomplish in `ggplot`.  While "model specification syntax" (try `aggregate(gold ~ nationality + sex, dat, sum)`) is nice, simultaneously aggregating multiple columns is also nice.  Of course, there are tools in `dyplyr` that help with this, but the base R commands are quite clear and powerful.

```{r}
#| label: aggregate-athelete
agg <- aggregate(dat[, c("gold", "silver", "bronze")], by = list(dat$nationality, dat$sex), sum)
head(agg)
names(agg)[1:2] <- c("nationality", "sex")
agg <- agg[order(agg$gold, decreasing = TRUE), ]
```

Given the summarized data, we can consider graphs, likely showing how some combinations of variables explain amounts and combinations of medals. Among you there was interest in making bar graphs showing medals by nationality or scatterplots showing heights or weights, with annotation for type of sport, sex, or other category.  Since there are so many nationalities and relatively many sports, using colors (or especially symbols) for those is not that helpful. Instead, you could create new variables that group nations by continent (or some other metric) or group sports by relevant descriptions (e.g., water, team, solo, field, gym).

```{r}
#| label: recode-sports
sport.cat <- dat$sport
sport.cat[sport.cat %in% c("rowing", "sailing", "aquatics")] <- "water"
```

After you were satisfied, you would probably want to run `table(sport.cat)` to make sure you hadn't missed any sports. Finally you could create this as a variable by adding a column to the original data `dat$sport.cat <- factor(sport.cat)`. The name of our variable doesn't have to match the new column we add to our data.

This requires some work and interpretation.  If you had a separate dataset that listed other classifiers,  you could merge that onto this dataset.  It may even be the case that some `sport` entries contain both team and solo sports, complicating that assignment.  The motivation, though, is a smaller number of distinct entries.
  
## Pharmaceutical spending data

We read in the "drug" data from <https://github.com/datasets/pharmaceutical-drug-spending>.

>Pharmaceutical Drug Spending by countries with indicators such as a share of total health spending, in USD per capita (using economy-wide PPPs) and as a share of GDP. Plus, total spending by each countries in the specific year.

A question was how we could identify which country had the greatest spending each year.  We can aggregate spekding by year using the `max` function to select which spending value was the annual maximum.  That doesn't tell us *where* the spending occurred. We could search for our with-year maximum spending values on a spreadsheet and take notes, or we could merge original data back to this, matching year and spending, and study the years.

```{r}
#label: data-drug
dat <- read.delim("./data/drug.csv", header = TRUE, sep = ',')
head(dat)
agg <- aggregate(dat$TOTAL_SPEND, by = list(dat$TIME), max)
agg2 <- aggregate(TOTAL_SPEND ~ TIME, dat, max)
head(agg2)
```

We merge original data `dat` onto our 'aggregated data' `agg`. Since both contain `TIME` and `TOTAL_SPEND`, we match by those to create a dataset that contains the corresponding country code (and other variable that remain in `dat`). If for some reason this didn't work out, we could have used `all.x` to keep only data present in `agg2`. Or, we could have used `all.y` to keep unique combinations of `TIME` and `TOTAL_SPEND` present in the original data, though this would have resulted in `NA` values being included in the merged result.

```{r}
#label: id-max-country
tmp <- merge(agg2, dat, by = c("TIME", "TOTAL_SPEND"))
head(tmp)
```

From here this dataset could be graphed and annotated.  Becuase the original axes of @fig-total-drug-spend are not that nice, we can use provided code to write powers of 10.

```{r}
#| label: fig-total-drug-spend
#| fig.width: 5
#| fig.height: 5
#| fig.cap: "Maximum total drug spending."
ggplot(tmp) + aes(x = TIME, y = TOTAL_SPEND, label = LOCATION, color = LOCATION) + geom_point() + scale_y_log10(breaks = scales::trans_breaks("log10", function(x) 10^x),
              labels = scales::trans_format("log10", scales::math_format(10^.x)))# + geom_text(hjust = 0)
```
I wouldn't make too much of the most recent, surprisingly low, year being Italy.  Without looking at this too closely, it would be *very* reasonable to hypothesize that result is based on incomplete data. This was visible in some of our "bad graphs" from newspapers and other published sources. It's reasonable to believe that there might be a specific way to do this without pre-calculating the maximum yearly expenses and identifying the corresponding country separately. Additionally, using a legend instead of the commented out `geom_text()` winds up looking quite a bit nicer.

## The mouse data

Just to revisit the ideas from the iNaturalist data using a smaller dataset, we will revisit the "mouse" data.  Using the optional argument `check_overlap = TRUE` will suppress one of any labels that are too close together to print without overlap. Using additional arguments as used in @fig-total-drug-spend will result in "nicer" tick mark labels. In fact, this is pretty bad because the rightmost label reads as `1e+0` which has been truncated from `1e+04`.

```{r}
#| label: fig-mouse
#| fig.width: 5
#| fig.height: 5
#| fig.cap: "A thoroughly annotated mouse-elephant curve."
dat <- read.delim("./data/mouse.csv", header = TRUE, sep = ',')
ggplot(dat, aes(Mass, MetabolicRate, label = Animal)) + geom_point() + scale_x_log10(limits = c(0.01, 50000)) + scale_y_log10() + geom_text(hjust = 0, check_overlap = FALSE, nudge_x = 0.1)
```

## Conclusion
If you have made it this far, be sure to let me know if you've encountered errors replicating this or have other questions or comments.
