---
title: "Introduction to R studio: Part II"
author: "DATASCI 120 Teaching Team"
date: "`r Sys.Date()`"
output:
  tufte::tufte_html: 
    citation_package: natbib
    latex_engine: xelatex
  tufte::tufte_handout:
    latex_engine: xelatex
    keep_tex: true
    pandoc_args: ["--variable=geometry:portrait"]
  tufte::tufte_book:
    citation_package: natbib
    latex_engine: xelatex
---

```{r setup, include=FALSE}
library(tufte)
library(ggplot2)
library(maps)
library(tidyverse)
library(reshape2)
# invalidate cache when the tufte version changes
knitr::opts_chunk$set(cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
```


# Displaying maps

We now consider in detail one particular type of graphics that is often not covered in standard courses: maps.

Each region on a map is represented as a polygon, which is essentially a shape made up of multiple points connected by lines to form a closed boundary. These polygons are crucial for mapping as they define the spatial boundaries of different geographic areas.

One popular package in R for working with maps is `maps`. The `maps` package provides a simple interface to access a variety of map databases, offering detailed map outlines for countries, states, and world regions. When these map databases are loaded, each region (like a country) is broken down into its polygonal components, each described by a series of longitude and latitude points. The `ggplot2` package can then take these points and visually render them as polygons, allowing for the addition of aesthetic mappings such as colors to represent different data variables (like population, GDP, or health statistics).^[Before running the code, make sure that the required packages `ggplot2` and `maps` are installed.
]

We will explore examples of two distinct types of maps: dot maps and choropleths.

- A dot map is a map overlaid with dots colored according to some statistic.
- A choropleth is a map where each region is colored according to some statistic.

## Dot maps

There is a built-in dataset called quakes that contains information on earthquake occurrences, detailing attributes like magnitude, depth, and location coordinates for seismic events in the Fiji area. Our goal is to produce a dot mat to visualize the magnitude and depth of seismic events. We shall approach this in three steps.

```{r quakes, echo=TRUE}
data("quakes")  # Load the earthquakes data
head(quakes)    # Display the first few rows
```

Step 1: World Map

First, let's use the `map_data()` function to retrieve data for the world map. The columns of the retrieved data set store the coordinates (longitude and latitude), along with identifiers that help in grouping and ordering these points to form the shapes of countries or regions on a map.

```{r retrieve,echo=TRUE}
world_map <- map_data("world")
head(world_map)
```

Now, we are ready to use `ggplot` to produce the map.

```{r maps1, echo=TRUE, fig.cap="Simple world map.", message=FALSE, cache=TRUE}
map1= ggplot() +
  geom_polygon(data = world_map, 
               aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "black") +
  theme_minimal()
map1+labs(x = "Longitude", y = "Latitude")
```

Step 2: Zooming in on a region.

Next, let's use `coord_fixed()` to zoom in on the Fiji area.

```{r maps1_2, fig.cap="World Map with Fiji earthquakes data.", fig.height=3.5, fig.margin=TRUE}
map2 = map1 +
  coord_fixed(ratio = 1, xlim = c(160, 190), ylim = c(-50, -10))
map2 + labs(x = "Longitude", y = "Latitude", 
             color = "Magnitude")
```

Step 3: Finally. we are ready to overlay the map with dots such that a deeper color indicates greater depth and a larger size indicates greater magnitude.

```{r map1_3,fig.cap="Earthquakes around Fiji.",fig.height=5,fig.width=4}
map2+
  geom_point(data = quakes, 
             aes(x = long, y = lat, size =depth , color = mag),
             alpha = 0.5) +
  scale_size(range = c(0.1, 2)) +
  scale_color_gradient(low = "blue", high = "red") +
  labs( x = "Longitude", y = "Latitude",
        color = "Magnitude",size="Depth") +
  theme_minimal()
```

## Choropleths

Next, let's make a choropleth of the United States with each state colored according the the percentage of the population living in urban areas in 1973. We will work with another built-in dataset called `USArrests`, which contains statistics on arrests per 100,000 residents for assault, murder, and rape in each of the U.S. states in 1973, along with the percentage of the population living in urban areas.

```{r usarrest,echo=TRUE}
data("USArrests")
head(USArrests)
```

```{r mapsusarrests, fig.cap="Urban population by state"}
# Get map data for US states (this is in the ggplot package, but you can also download your own lat, long coordinates)
states_map <- map_data("state")

# Merge map data with USArrests data
USArrests$state <- tolower(rownames(USArrests))
arrest_map <- states_map %>%
  left_join(USArrests, by = c("region" = "state"))

# Plot with ggplot2
ggplot(arrest_map, aes(long, lat, group = group, fill = UrbanPop)) +
  geom_polygon(color = "white", linewidth = 0.3) +
  coord_fixed(1.3) +
  scale_fill_gradientn(
    colors = colorRampPalette(c("lightblue", "darkblue"))(100),
    name = "Urban Population (%)"
  ) +
  labs(title = "Urban Population by US State") +
  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5)
  )

```

An alternative to  making choropleths is using the `maps()` package, instead of using `ggplot()`.

# Working with Categorical Variables
Categorical variables in R are stored in a data type called a factor. It stores both the values of the variable and its unique categories, known as the levels of a factor. Even if they are not explicitly specified, R can determine a factor's levels. We will discuss how to work with qualitative and ordinal factors through two example datasets.

## Renaming Variables and Levels
Let's first take a look at the `iris` dataset, which contains 150 total measurements of sepal and petal lengths and widths for three species of Iris flowers. `Species` is a qualitative variable stored as a factor, so we can print the levels and see that the species included are Setosa, Versicolor, and Virginica. We can make plots to visually observe differences between the species.

```{r}
data("iris")
head(iris)
```
```{r}
levels(iris$Species)
```

```{r fig.cap="Petal width and length for Iris species."}
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, col = Species)) + 
  geom_point() + 
  labs(x = "Petal Length", y = "Petal Width") + 
  theme_minimal()
```

Even though the variable name `Species` is informative, someone who is unfamiliar with the dataset or flowers may not know that these are referring to Iris species. We can either rename the variable or the levels to clarify exactly what flowers these correspond to. To rename the variable, we modify the vector returned by `names()` which contains all the column names of `iris`. Alternatively, `tidyverse` provides a function `rename` that gives the same result.

```{r}
# Rename variable
iris_copy <- data.frame(iris)
names(iris_copy)[names(iris_copy) == "Species"] <- "Iris Species"
head(iris_copy)
```

```{r message=FALSE}
# Rename variable using tidyverse
library(tidyverse)
iris_copy <- iris %>% rename("Iris Species" = Species)
```

We could also rename the levels instead of the variable name.

```{r fig.cap="Petal width and length for Iris species."}
iris_copy <- data.frame(iris)
levels(iris_copy$Species)

# To rename specific levels
levels(iris_copy$Species)[levels(iris_copy$Species) == "versicolor"] <- "Iris versicolor"
levels(iris_copy$Species)

# To rename all levels
levels(iris_copy$Species) <- c("Iris setosa", "Iris versicolor", "Iris virginica")
levels(iris_copy$Species)

ggplot(iris_copy, aes(x = Petal.Length, y = Petal.Width, col = Species)) + 
  geom_point() + 
  labs(x = "Petal Length", y = "Petal Width") + 
  theme_minimal()
```

## Creating factors
The next example will show how to create a factor, whether for qualitative or ordinal variables. In some datasets, categorical variables may be stored as another data type and should be converted to a factor first.

Ordinal variables are categorical variables where the categories have an inherent order. One example is the months of a year (January, February, etc.), which are ordered and often represented numerically (1, 2, etc.). In the `airquality` dataset, there are daily air quality measurements from May to September 1973 in New York. We demonstrate how to encode the ordinal variable `Month` as a factor.

```{r}
data("airquality")
head(airquality)
```

Currently, `Month` is stored as integers in the dataframe. To encode a qualitative variable as a factor, we can use the function `factor` or `as.factor`.

```{r}
airquality$Month <- as.factor(airquality$Month)
airquality$Month
```
For an ordinal variable, we want to also specify that the levels are ordered. This allows us to use comparison operators which would not be defined for the unordered factor above.

```{r}
airquality$Month <- factor(airquality$Month, ordered = TRUE,
                           levels = c("5", "6", "7", "8", "9"))
airquality$Month
```
Finally, let's rename the levels to be more informative than its original numerical values. Now we can make plots using `Month` as a categorical variable.

```{r}
levels(airquality$Month) <- c("May", "June", "July", "August", "September")
```

```{r fig.cap="Temperature in New York 1973."}
ggplot(airquality, aes(x = Month, y = Temp)) +
  geom_boxplot() +
  labs(y = "Temperature") +
  theme_minimal()
```

# Facets

Facets are a powerful tool for visualizing how the relationship between two or 
more variables depends on other variables. For instance, let's say we're
interested in understanding the relationship between petal length and petal
width for flowers. Depending on the species, the nature of this relationship 
might differ, so it can be helpful to have separate plots for each flower type
in our data. The `facet_grid` method makes this easy:

```{r fig.cap="Petal length vs width by species"}
plot <- ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) + 
  geom_point() + 
  labs(x = "Petal Length", y = "Petal Width") + 
  theme_minimal()
plot + facet_grid(. ~ Species)
```

Note that the above generated horizontally organized subplots. If vertically 
organized subplots are desired, then instead of `. ~ [variable]`, you should 
use `[variable] ~ .` as an input to `facet_grid`. 

```{r fig.cap="Vertical Subplots"}
plot + facet_grid(Species ~ .)
```

The `facet_grid` method can also be used to generate subplots for each 
combination of two categorical variables. To see this in action, let's start by
creating a new variable that indicates whether a flower has sepal width greater
than half its sepal length. 

```{r}
iris <- mutate(iris, wide_sepals = ifelse(Sepal.Width > 0.5 * Sepal.Length,
               "Wide", "Thin"))
head(iris)
```

Next, using the input `wide_sepals ~ Species`, we generate subplots that are
vertically divided by sepal width and horizontally divided by species.

```{r fig.cap="Example of Two-Variable Facet Grid"}
plot <- ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) + 
  geom_point() + 
  labs(x = "Petal Length", y = "Petal Width") + 
  theme_minimal()
plot + facet_grid(wide_sepals ~ Species)
```

The above plot makes it visually clear that the ratio between sepal width and 
sepal length doesn't seem to affect the relationship between petal width and 
petal length for any of the three species in the `iris` data. Depending on your
data, the best method for visualizing specific relationships might involve a 
combination of color-coding and facets. An example is below:

```{r fig.cap="Example of Combining Facets and Color-Coding"}
plot <- ggplot(iris, aes(x = Petal.Length, y = Petal.Width, col = Species)) + 
  geom_point() + 
  labs(x = "Petal Length", y = "Petal Width") + 
  theme_minimal()
plot + facet_grid(. ~ wide_sepals)
```

## Changing colors 

We can also easily change the colors if we wanted to using `scale_color_manual`. There are a variety of colors, some of which you can see here: https://sape.inf.usi.ch/quick-reference/ggplot2/colour. 

Generally, it is a good idea to use colors consistently throughout your document. For example, if you always display points or lines for setosa in orange, then in any figure, the data related to setosa should be colored orange. 

```{r fig.cap="Example of Changing Colors"}
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, col = Species)) + 
  geom_point() + 
  labs(x = "Petal Length", y = "Petal Width") + 
  facet_grid(. ~ wide_sepals) + 
  theme_minimal() + 
  scale_color_manual(values = c("darkorange", "cyan4", "lightcoral"), 
                     labels = c("Setosa", "Versicolor", "Virgincia"))

```


## Facet Wrap

By default, calling `facet_grid([Variable 1] ~ [Variable 2])` will create an
 grid of subplots, where are the number of categories in
`Variable 1` and `Variable 2` respectively. To manually control the dimensions
of your facet grid, you can use the `facet_wrap()` method, which takes a number
of columns `ncol` as an argument. In the following example, we set the number
of columns to 2:

```{r fig.cap="Example of Facet Wrap"}
plot <- ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) + 
  geom_point() + 
  labs(x = "Petal Length", y = "Petal Width") + 
  theme_minimal()
plot + facet_wrap( ~ Species, ncol = 2)
```


## Pie charts 

Below we make a pie chart for species in the Iris dataset. Note that here there is the same number observations for each species, but we do include a pie chart below for illustration purposes. We also include labeling by percentages, just to illustrate the approach. You can modify the plot below for your own purposes. 

```{r fig.cap="Example of Pie Charts"}
# Count species
species_counts <- as.data.frame(table(iris$Species))
colnames(species_counts) <- c("Species", "Count")

# Calculate percentages and labels
species_counts$Percent <- round(100 * species_counts$Count / sum(species_counts$Count), 1)
species_counts$Label <- paste0(species_counts$Species, "\n", species_counts$Percent, "%")

# Plot with ggplot
ggplot(species_counts, aes(x = "", y = Count, fill = Species)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y") +  # turns bar into pie
  theme_void() +  # remove background
  geom_text(aes(label = Label), 
            position = position_stack(vjust = 0.5)) +
  labs(title = "Distribution of Iris Species")
```


## Timeseries 

The `lubridate` package, which is part of `tidyverse` is designed to handle dates in R. Here is an example: 

```{r}
ymd(c("20160927", "20160230"))
```

```{r}
mdy(c("6/12/16", "2/9/16"))
dmy(c("1/9/2016", "26/9/16"))
```


We will use data from New York's Department of Sanitation below. The data was downloaded from here https://data.cityofnewyork.us/City-Government/DSNY-Monthly-Tonnage-Data/ebb7-mvp5.

You can find the dataset on Canvas in the Rintro folder: https://canvas.stanford.edu/courses/206498/files/folder/Rintro. Download the file and put it into your data directory. 


```{r}
NY_waste <- read.csv("NY waste.csv")
head(NY_waste)
```


This data set contains information about the amount of waste collected in different locations in New York city. The first column specifies the date of collection. The next two columns specify the borough and district (location). The next column specifies the type of rubbish collection ranging from refuse (general waste) to Christmas trees. Let's use the waste dataset to look at trends over time across all of New York city. To do this we should group by date and type and use `summarize()` to add up the amount of waste across all districts.

```{r}
aggregate_waste <- NY_waste %>% 
  group_by(date, type) %>% 
  summarize(amount = sum(value))

aggregate_waste
```

Let's first look at the overall proportions of waste across time. Let's consider a pie-chart for this:

```{r piechartwaste, fig.cap="Pie chart of waste by type"}
# Summarize total waste by type
waste_summary <- aggregate_waste %>%
  group_by(type) %>%
  summarise(total_amount = sum(amount)) %>%
  ungroup()

# Create pie chart
ggplot(waste_summary, aes(x = "", y = total_amount, fill = type)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y") +
  labs(title = "Total Waste by Type", fill = "Waste Type") +
  theme_void()
```

We see that most waste comes from paper, mgp and refuse. Next, we are interested in exploring trends over time, for different types: 

```{r wastetrends, fig.cap="Time series waste by type"}
# Make sure that date is date type 
aggregate_waste <- aggregate_waste %>%
  mutate(date = as.Date(date))

# Plot
ggplot(aggregate_waste, aes(x = date, y = amount)) +
  geom_line(color = "steelblue") +
  facet_wrap(~ type, scales = "free") + 
  labs(title = "Time Series of Waste by Type",
       x = "Date",
       y = "Amount of waste") +
  theme_minimal()
```

Based on the above, we see some interesting patterns: For example, christmas trees spike at certain times, and school organics are zero during certain times. Let's calculate the total by year and month, and filter out refuse, mgp and paper: 

```{r}
monthly_waste <- aggregate_waste %>% 
  mutate(month = month(date),
         year = year(date)) %>%
  group_by(month, type) %>% 
  summarize(amount = mean(amount)) %>% 
  ungroup()

monthly_waste_filtered = monthly_waste %>% 
  filter(!(type %in% c("refuse", "mgp", "paper")))
```


```{r wastetrendsmonthly, fig.cap="Waste by month and type"}
monthly_waste_filtered %>%
  ggplot(aes(x = month, y = amount, color = type)) +
  geom_line() +
  scale_x_continuous(breaks = (1:12))  + 
  labs(title = "Waste by Month",
       x = "Month",
       y = "Amount of waste", 
       col = "Waste Type")
```

We can see a lot of seasonal effects in the last two graphics. Three that stand out are:

- A large number of Christmas trees are collected in January but not in any other month.
- Leaves are collected in the fall. There's a little bit in October, a lot in November and December and again a little bit in January and February.
- There's no organic waste from the schools in July and August! This makes sense as schools are closed for summer.


Let's say we are interested in seeing whether there are any patter by whether it was the Covid year (2020) or not (previous).

```{r comparingcovid, fig.cap="Comparing waste in covid versus non-covid year"}
refuse_2020_b <- NY_waste %>% 
  mutate(year = year(date),
         month = month(date),
         is_2020 = (year == 2020)) %>% 
  filter(year <= 2020,
         type == "refuse") %>% 
  group_by(borough, month, is_2020) %>% 
  summarize(amount = mean(value)) 

ggplot(refuse_2020_b) +
  geom_line(aes(x = month, y = amount, color = is_2020)) +
  facet_wrap(~borough) + 
  labs(x = "Month", y = "Amount", color = "2020?") + 
  scale_x_continuous(breaks = seq(1, 12))
```
We see that during Covid, a business district such as Manhattan had lower waste. However, Bronx, Brooklyn, Queens and Staten Island had higher waste in the second half of the year. 

# Plots - a gallery


The  codes and examples in this section are done using a standard demo package, with its advantages and disadvantages. Take them simply as a preview of what is possible and the opportunity to learn a bit more about any of these if you are interested, including how to customize them to maximize effectiveness.

```{r plot}
demo(graphics)
```



