94-appendixD.Rmd

# Learning Check Solutions {#appendixD}

```{r setup_lc_solutions, include=FALSE, purl=FALSE}
opts_chunk$set(tidy = FALSE, out.width = "\\textwidth")
# This bit of code is a bug fix on asis blocks, which we use to show/not show LC solutions, which are written like markdown text. In theory, it shouldn't be necessary for knitr versions <=1.11.6, but I've found I still need to for everything to knit properly in asis blocks. More info here:
# https://stackoverflow.com/questions/32944715/conditionally-display-block-of-markdown-text-using-knitr
knit_engines$set(asis = function(options) {
  if (options$echo && options$eval) knit_child(text = options$code)
})
```

Here are the solutions to all Learning Checks throughout the book.



## Chapter 1 Solutions 

```{r, include=FALSE, purl=FALSE}
chap <- 1
lc <- 0
```

```{r message=FALSE}
library(dplyr)
library(ggplot2)
library(nycflights13)
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  Repeat the above installing steps, but for the `dplyr`, `nycflights13`, and `knitr` packages. This will install the earlier mentioned `dplyr` package, the `nycflights13` package containing data on all domestic flights leaving a NYC airport in 2013, and the `knitr` package for writing reports in R. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** "Load" the `dplyr`, `nycflights13`, and `knitr` packages as well by repeating the above steps.

**Solution**: If the following code runs with no errors, you've succeeded!

```{r, eval=FALSE}
library(dplyr)
library(nycflights13)
library(knitr)
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What does any *ONE* row in this `flights` dataset refer to?

- A. Data on an airline 
- B. Data on a flight
- C. Data on an airport
- D. Data on multiple flights

**Solution**: This is data on a flight. Not a flight path! Example: 

* a flight path would be United 1545 to Houston
* a flight would be United 1545 to Houston at a specific date/time. For example: 2013/1/1 at 5:15am.


**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What are some examples in this dataset of **categorical** variables?  What makes them different than **quantitative** variables?

**Solution**: Hint: Type `?flights` in the console to see what all the variables mean!

* Categorical:
    + `carrier` the company
    + `dest` the destination
    + `flight` the flight number. Even though this is a number, its simply a label. Example United 1545 is not  less than United 1714
* Quantitative:
    + `distance` the distance in miles
    + `time_hour` time

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What properties of the observational unit do each of `lat`, `lon`, `alt`, `tz`, `dst`, and `tzone` describe for the `airports` data frame?  Note that you may want to use `?airports` to get more information.

**Solution**: `lat` and `long` represent the airport geographic coordinates, `alt` is the altitude above sea level of the airport (Run `airports %>% filter(faa == "DEN")` to see the altitude of Denver International Airport), `tz` is the time zone difference with respect to GMT in London UK, `dst` is the daylight savings time zone, and `tzone` is the time zone label.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Provide the names of variables in a data frame with at least three variables in which one of them is an identification variable and the other two are not.  In other words, create your own tidy dataset that matches these conditions.

**Solution**: 

* In the `weather` example in LC2.8, the combination of `origin`, `year`, `month`, `day`, `hour`  are identification variables as they identify the observation in question.
* Anything else pertains to observations: `temp`, `humid`, `wind_speed`, etc.



***



## Chapter 2 Solutions 

```{r, include=FALSE, purl=FALSE}
chap <- 2
lc <- 0
```

```{r message=FALSE}
library(nycflights13)
library(ggplot2)
library(dplyr)
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  Take a look at both the `flights` data frame from the `nycflights13` package and the `alaska_flights` data frame from the `moderndive` package by running `View(flights)` and `View(alaska_flights)`. In what respect do these data frames differ? For example, think about the number of rows in each dataset.

**Solution**: `flights` contains all flight data, while `alaska_flights` contains only data from Alaskan carrier "AS". We can see that flights has `r nrow(flights)` rows while `alaska_flights` has only `r nrow(alaska_flights)`.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  What are some practical reasons why `dep_delay` and `arr_delay` have a positive relationship?

**Solution**: The later a plane departs, typically the later it will arrive.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  What variables in the `weather` data frame would you expect to have a negative correlation (i.e. a negative relationship) with `dep_delay`? Why? Remember that we are focusing on numerical variables here. Hint: Explore the `weather` dataset by using the `View()` function. 

**Solution**: An example in the `weather` dataset is `visibility`, which measure visibility in miles. As visibility increases, we would expect departure delays to decrease. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why do you believe there is a cluster of points near (0, 0)? What does (0, 0) correspond to in terms of the Alaskan flights?

**Solution**: The point (0,0) means no delay in departure nor arrival. From the point of view of Alaska airlines, this means the flight was on time. It seems most flights are at least close to being on time.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What are some other features of the plot that stand out to you?

**Solution**: Different people will answer this one differently. One answer is most flights depart and arrive less than an hour late. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Create a new scatterplot using different variables in the `alaska_flights` data frame by modifying the example above.

**Solution**: Many possibilities for this one, see the plot below. Is there a pattern in departure delay depending on when the flight is scheduled to depart? Interestingly, there seems to be only two blocks of time where flights depart. 

```{r}
ggplot(data = alaska_flights, mapping = aes(x = dep_time, y = dep_delay)) +
  geom_point()
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  Why is setting the `alpha` argument value useful with scatterplots? What further information does it give you that a regular scatterplot cannot?

**Solution**: It thins out the points so we address overplotting. But more importantly it hints at the (statistical) **density** and **distribution** of the points: where are the points concentrated, where do they occur.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** After viewing the Figure \@ref(fig:alpha) above, give an approximate range of arrival delays and departure delays that occur the most frequently.  How has that region changed compared to when you observed the same plot without the `alpha = 0.2` set in Figure \@ref(fig:noalpha)?

**Solution**: The lower plot suggests that most Alaska flights from NYC depart between 12 minutes early and on time and arrive between 50 minutes early and on time.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Take a look at both the `weather` data frame from the `nycflights13` package and the `early_january_weather` data frame from the `moderndive` package by running `View(weather)` and `View(early_january_weather)`. In what respect do these data frames differ?

**Solution**: The rows of `early_january_weather` are a subset of `weather`.  

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** `View()` the `flights` data frame again. Why does the `time_hour` variable uniquely identify the hour of the measurement whereas the `hour` variable does not? 

**Solution**: Because to uniquely identify an hour, we need the `year`/`month`/`day`/`hour` sequence, whereas there are only 24 possible `hour`'s. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  Why should linegraphs be avoided when there is not a clear ordering of the horizontal axis?

**Solution**: Because lines suggest connectedness and ordering.  

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why are linegraphs frequently used when time is the explanatory variable?

**Solution**: Because time is sequential: subsequent observations are closely related to each other.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Plot a time series of a variable other than `temp` for Newark Airport in the first 15 days of January 2013.

**Solution**: Humidity is a good one to look at, since this very closely related to the cycles of a day.
```{r}
ggplot(data = early_january_weather, mapping = aes(x = time_hour, y = humid)) +
  geom_line()
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What does changing the number of bins from 30 to 40 tell us about the distribution of temperatures?

**Solution**: The distribution doesn't change much. But by refining the bin width, we see that the temperature data has a high degree of accuracy. What do I mean by accuracy? Looking at the `temp` variable by `View(weather)`, we see that the precision of each temperature recording is 2 decimal places.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Would you classify the distribution of temperatures as symmetric or skewed?

**Solution**: It is rather symmetric, i.e. there are no __long tails__ on only one side of the distribution.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What would you guess is the "center" value in this distribution?  Why did you make that choice?

**Solution**: The center is around `r mean(weather$temp, na.rm=TRUE)`&deg;F. By running the `summary()` command, we see that the mean and median are very similar. In fact, when the distribution is symmetric the mean equals the median.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Is this data spread out greatly from the center or is it close?  Why?

**Solution**: This can only be answered relatively speaking! Let's pick things to be relative to Seattle, WA temperatures:

```{r SEATAC, out.width='40%', echo=FALSE, fig.cap="Annual temperatures at SEATAC Airport.", purl=FALSE}
include_graphics("images/SEATAC_temperatures.png")
```

While, it appears that Seattle weather has a similar center of 55&deg;F, its
temperatures are almost entirely between 35&deg;F and 75&deg;F for a range of
about 40&deg;F. Seattle temperatures are much less spread out than New York i.e.
much more consistent over the year. New York on the other hand has much colder
days in the winter and much hotter days in the summer. Expressed differently,
the middle 50% of values, as delineated by the interquartile range is 30&deg;F.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What other things do you notice about the faceted plot above?  How does a faceted plot help us see relationships between two variables?

**Solution**: 

* Certain months have much more consistent weather (August in particular), while others have crazy variability like January and October, representing changes in the seasons.
* Because we see `temp` recordings split by `month`, we are considering the relationship between these two variables. For example, for summer months, temperatures tend to be higher. 
**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What do the numbers 1-12 correspond to in the plot above?  What about 25, 50, 75, 100?

**Solution**: 

* They correspond to the month of the flight. While month is technically a number between 1-12, we're viewing it as a categorical variable here. Specifically, this is an **ordinal categorical** variable since there is an ordering to the categories.
* 25, 50, 75, 100 are temperatures.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** For which types of datasets would these types of faceted plots not work well in comparing relationships between variables? Give an example describing the nature of these variables and other important characteristics.

**Solution**: 

* It would not work if we had a very large number of facets. For example, if we faceted by individual days rather than months, as we would have 365 facets to look at. When considering all days in 2013, it could be argued that we shouldn't care about day-to-day fluctuation in weather so much, but rather month-to-month fluctuations, allowing us to focus on seasonal trends.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Does the `temp` variable in the `weather` dataset have a lot of variability?  Why do you say that?

**Solution**: Again, like in LC `r paste0("(LC", chap, ".", (lc-4), ")")`, this is a relative question. I would say yes, because in New York City, you have 4 clear seasons with different weather. Whereas in Seattle WA and Portland OR, you have two seasons: summer and rain!

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What does the dot at the bottom of the plot for May correspond to?  Explain what might have occurred in May to produce this point.

**Solution**: It appears to be an outlier. Let's revisit the use of the `filter` command to hone in on it. We want all data points where the `month` is 5 and `temp<25`.

```{r}
weather %>%
  filter(month == 5 & temp < 25)
```

There appears to be only one hour and only at JFK that recorded 13.1 F (-10.5 C) in the month of May. This is probably a data entry mistake! Why wasn't the weather at least similar at EWR (Newark) and LGA (LaGuardia)?

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Which months have the highest variability in temperature?  What reasons do you think this is?

**Solution**: We are now interested in the **spread** of the data. One measure some of you may have seen previously is the standard deviation. But in this plot we can read off the Interquartile Range (IQR):

* The distance from the 1st to the 3rd quartiles i.e. the length of the boxes
* You can also think of this as the spread of the **middle 50%** of the data

Just from eyeballing it, it seems

* November has the biggest IQR, i.e. the widest box, so has the most variation in temperature
* August has the smallest IQR, i.e. the narrowest box, so is the most consistent temperature-wise

Here's how we compute the exact IQR values for each month (we'll see this more in depth Chapter \@ref(wrangling) of the text):

1. `group` the observations by `month` then
1. for each `group`, i.e. `month`, `summarize` it by applying the summary statistic function `IQR()`, while making sure to skip over missing data via `na.rm=TRUE` then
1. `arrange` the table in `desc`ending order of `IQR`

```{r, eval=FALSE}
weather %>%
  group_by(month) %>%
  summarize(IQR = IQR(temp, na.rm = TRUE)) %>%
  arrange(desc(IQR))
```
```{r, echo=FALSE}
weather %>%
  group_by(month) %>%
  summarize(IQR = IQR(temp, na.rm = TRUE)) %>%
  arrange(desc(IQR)) %>%
  kable()
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** We looked at the distribution of the numerical variable `temp` split by the numerical variable `month` that we converted to a categorical variable using the `factor()` function. Why would a boxplot of `temp` split by the numerical variable `pressure` similarly converted to a categorical variable using the `factor()` not be informative?

**Solution**: Because there are 12 unique values of `month` yielding only 12 boxes in our boxplot. There are many more unique values of `pressure` (`r weather$pressure %>% unique() %>% length()` unique values in fact), because values are to the first decimal place. This would lead to `r weather$pressure %>% unique() %>% length()` boxes, which is too many for people to digest. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Boxplots provide a simple way to identify outliers.  Why may outliers be easier to identify when looking at a boxplot instead of a faceted histogram?

**Solution**: In a histogram, the bin corresponding to where an outlier lies may not by high enough for us to see. In a boxplot, they are explicitly labelled separately. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why are histograms inappropriate for visualizing categorical variables?

**Solution**: Histograms are for numerical variables i.e. the horizontal part of each histogram bar represents an interval, whereas for a categorical variable each bar represents only one level of the categorical variable.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What is the difference between histograms and barplots?

**Solution**: See above.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** How many Envoy Air flights departed NYC in 2013?

**Solution**: Envoy Air is carrier code `MQ` and thus 26397 flights departed NYC in 2013.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What was the seventh highest airline in terms of departed flights from NYC in 2013? How could we better present the table to get this answer quickly?

**Solution**: The answer is US, AKA U.S. Airways, with 20536 flights. However, picking out the seventh highest airline when the rows are sorted alphabetically by carrier code is difficult. This would be easier to do if the rows were sorted by number. We'll learn how to do this in Chapter \@ref(wrangling) on data wrangling.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why should pie charts be avoided and replaced by barplots?

**Solution**: In our **opinion**, comparisons using horizontal lines are easier than comparing angles and areas of circles.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What is your opinion as to why pie charts continue to be used?

**Solution**: In our **opinion**, pie charts are generally considered as a poorer method for communicating data than bar charts. People's brains are not as good at comparing the size of angles because there is no scale, and in comparison, it is much easier to compare the heights of bars in a bar charts. However, in some circumstances, for example, when representing 25% and 75% of a sample size, if we have 2 bars, in which the higher one is three times in height of the other one, it is difficult to tell the scale of their comparison without labels. But in a bar chart, it would be easy to compare if a circle is divided by 75% and 25%. (Read more at: https://www.displayr.com/why-pie-charts-are-better-than-bar-charts/)

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What kinds of questions are not easily answered by looking at the above figure?

**Solution**: Because the red, green, and blue bars don't all start at 0 (only red does), it makes comparing counts hard.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What can you say, if anything, about the relationship between airline and airport in NYC in 2013 in regards to the number of departing flights?

**Solution**: The different airlines prefer different airports. For example, United is mostly a Newark carrier and JetBlue is a JFK carrier. If airlines didn't prefer airports, each color would be roughly one third of each bar.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why might the side-by-side (AKA dodged) barplot be preferable to a stacked barplot in this case?

**Solution**: We can easily compare the different airports for a given carrier using a single comparison line i.e. things are lined up

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What are the disadvantages of using a side-by-side (AKA dodged) barplot, in general?

**Solution**: It is hard to get totals for each airline.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why is the faceted barplot preferred to the side-by-side and stacked barplots in this case?

**Solution**: Not that different than using side-by-side; depends on how you want to organize your presentation.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What information about the different carriers at different airports is more easily seen in the faceted barplot?

**Solution**: Now we can also compare the different carriers **within** a particular airport easily too. For example, we can read off who the top carrier for each airport is easily using a single horizontal line.



***



## Chapter 3 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 3
lc <- 0
```

```{r message=FALSE}
library(dplyr)
library(ggplot2)
library(nycflights13)
``` 


**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What's another way using the "not" operator `!` to filter only the rows that are not going to Burlington, VT nor Seattle, WA in the `flights` data frame? Test this out using the code above.

**Solution**: 

```{r, eval=FALSE, purl=FALSE}
# Original in book
not_BTV_SEA <- flights %>%
  filter(!(dest == "BTV" | dest == "SEA"))

# Alternative way
not_BTV_SEA <- flights %>%
  filter(!dest == "BTV" & !dest == "SEA")

# Yet another way
not_BTV_SEA <- flights %>%
  filter(dest != "BTV" & dest != "SEA")
```


**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Say a doctor is studying the effect of smoking on lung cancer for a large number of patients who have records measured at five year intervals. She notices that a large number of patients have missing data points because the patient has died, so she chooses to ignore these patients in her analysis. What is wrong with this doctor's approach?

**Solution**: The missing patients may have died of lung cancer! So to ignore them might seriously **bias** your results! It is very important to think of what the consequences on your analysis are of ignoring missing data! Ask yourself:

* There is a systematic reasons why certain values are missing? If so, you might be biasing your results!
* If there isn't, then it might be ok to "sweep missing values under the rug."

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Modify the above `summarize` function to create `summary_temp` to also use the `n()` summary function: `summarize(count = n())`. What does the returned value correspond to?

**Solution**: It corresponds to a count of the number of observations/rows:

```{r, purl=FALSE}
weather %>%
  summarize(count = n())
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why doesn't the following code work?  Run the code line by line instead of all at once, and then look at the data.  In other words, run `summary_temp <- weather %>% summarize(mean = mean(temp, na.rm = TRUE))` first.

```{r eval=FALSE}
summary_temp <- weather %>%
  summarize(mean = mean(temp, na.rm = TRUE)) %>%
  summarize(std_dev = sd(temp, na.rm = TRUE))
```

**Solution**: Consider the output of only running the first two lines:

```{r, purl=FALSE}
weather %>%
  summarize(mean = mean(temp, na.rm = TRUE))
```

Because after the first `summarize()`, the variable `temp` disappears as it has been collapsed to the value `mean`. So when we try to run the second `summarize()`, it can't find the variable `temp` to compute the standard deviation of. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Recall from Chapter \@ref(viz) when we looked at plots of temperatures by months in NYC. What does the standard deviation column in the `summary_monthly_temp` data frame tell us about temperatures in New York City throughout the year?

**Solution**: 

```{r, echo=FALSE, purl=FALSE}
summary_temp_by_month <- weather %>%
  group_by(month) %>%
  summarize(
    mean = mean(temp, na.rm = TRUE),
    std_dev = sd(temp, na.rm = TRUE)
  )

kable(summary_temp_by_month) %>%
  kable_styling(
    font_size = ifelse(is_latex_output(), 10, 16),
    latex_options = c("hold_position")
  )
```

The standard deviation is a quantification of **spread** and **variability**. We
see that the period in November, December, and January has the most variation in
weather, so you can expect very different temperatures on different days.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What code would be required to get the mean and standard deviation temperature for each day in 2013 for NYC?

**Solution**: 

```{r, purl=FALSE}
summary_temp_by_day <- weather %>%
  group_by(year, month, day) %>%
  summarize(
    mean = mean(temp, na.rm = TRUE),
    std_dev = sd(temp, na.rm = TRUE)
  )
summary_temp_by_day
```

Note: `group_by(day)` is not enough, because `day` is a value between 1-31. We need to `group_by(year, month, day)`. 

```{r, purl=FALSE}
library(dplyr)
library(nycflights13)

summary_temp_by_month <- weather %>%
  group_by(month) %>%
  summarize(
    mean = mean(temp, na.rm = TRUE),
    std_dev = sd(temp, na.rm = TRUE)
  )
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Recreate `by_monthly_origin`, but instead of grouping via `group_by(origin, month)`, group variables in a different order `group_by(month, origin)`. What differs in the resulting dataset?

**Solution**: 

```{r, purl=FALSE}
by_monthly_origin <- flights %>%
  group_by(month, origin) %>%
  summarize(count = n())
```

```{r, eval=FALSE, purl=FALSE}
by_monthly_origin
```

```{r, echo=FALSE, purl=FALSE}
kable(by_monthly_origin) %>%
  kable_styling(
    font_size = ifelse(is_latex_output(), 10, 16),
    latex_options = c("hold_position")
  )
```

In `by_monthly_origin` the `month` column is now first and the rows are sorted by `month` instead of origin. If you compare the values of `count` in `by_origin_monthly` and `by_monthly_origin` using the `View()` function, you'll see that the values are actually the same, just presented in a different order. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** How could we identify how many flights left each of the three airports for each `carrier`?

**Solution**: We could summarize the count from each airport using the `n()` function, which *counts rows*.

```{r, purl=FALSE}
count_flights_by_airport <- flights %>%
  group_by(origin, carrier) %>%
  summarize(count = n())
```

```{r, eval=FALSE, purl=FALSE}
count_flights_by_airport
```

```{r, echo=FALSE, purl=FALSE}
kable(count_flights_by_airport) %>%
  kable_styling(
    font_size = ifelse(is_latex_output(), 10, 16),
    latex_options = c("hold_position")
  )
```

All remarkably similar! Note: the `n()` function counts rows, whereas the `sum(VARIABLE_NAME)` function sums all values of a certain numerical variable `VARIABLE_NAME`.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** How does the `filter` operation differ from a `group_by` followed by a `summarize`?

**Solution**: 

* `filter` picks out rows from the original dataset without modifying them, whereas
* `group_by %>% summarize` computes summaries of numerical variables, and hence
reports new values.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What do positive values of the `gain` variable in `flights` correspond to?  What about negative values?  And what about a zero value?

**Solution**: 

* Say a flight departed 20 minutes late, i.e. `dep_delay = 20`.
* Then arrived 10 minutes late, i.e. `arr_delay = 10`.
* Then `gain = dep_delay - arr_delay = 20 - 10  = 10` is positive, so it "made up/gained time in the air."
* 0 means the departure and arrival time were the same, so no time was made up in  the air. We see in most cases that the `gain` is near 0 minutes.
* I never understood this. If the pilot says "we're going make up time in the air"
because of delay by flying faster, why don't you always just fly faster to begin
with?

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Could we create the `dep_delay` and `arr_delay` columns by simply subtracting `dep_time` from `sched_dep_time` and similarly for arrivals?  Try the code out and explain any differences between the result and what actually appears in `flights`.

**Solution**: No because you can't do direct arithmetic on times. The difference in time between 12:03 and 11:59 is 4 minutes, but `1203-1159 = 44`

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What can we say about the distribution of `gain`?  Describe it in a few sentences using the plot and the `gain_summary` data frame values.

**Solution**: Most of the time the gain is a little above zero (the median is 7, meaning gain is above 0 at least 50% of the time) and between -50 and 50 minutes. There are some extreme cases however!

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Looking at Figure \@ref(fig:reldiagram), when joining `flights` and `weather` (or, in other words, matching the hourly weather values with each flight), why do we need to join by all of `year`, `month`, `day`, `hour`, and `origin`, and not just `hour`?

**Solution**: Because `hour` is simply a value between 0 and 23; to identify a *specific* hour, we need to know which year, month, day and at which airport.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What surprises you about the top 10 destinations from NYC in 2013?

**Solution**: This question is subjective! What surprises me is the high number of flights to Boston. Wouldn't it be easier and quicker to take the train?

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What are some advantages of data in normal forms?  What are some disadvantages?

**Solution**: When datasets are in normal form, we can easily `_join` them with other datasets! For example, we can join the `flights` data with the `planes` data.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What are some ways to select all three of the `dest`, `air_time`, and `distance` variables from `flights`?  Give the code showing how to do this in at least three different ways.

**Solution**: 

```{r, purl=FALSE}
# The regular way:
flights %>%
  select(dest, air_time, distance)

# Since they are sequential columns in the dataset
flights %>%
  select(dest:distance)

# Not as effective, by removing everything else
flights %>%
  select(
    -year, -month, -day, -dep_time, -sched_dep_time, -dep_delay, -arr_time,
    -sched_arr_time, -arr_delay, -carrier, -flight, -tailnum, -origin,
    -hour, -minute, -time_hour
  )
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** How could one use `starts_with`, `ends_with`, and `contains` to select columns from the `flights` data frame?  Provide three different examples in total: one for `starts_with`, one for `ends_with`, and one for `contains`.

**Solution**: 

```{r, purl=FALSE}
# Anything that starts with "d"
flights %>%
  select(starts_with("d"))
# Anything related to delays:
flights %>%
  select(ends_with("delay"))
# Anything related to departures:
flights %>%
  select(contains("dep"))
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why might we want to use the `select()` function on a data frame?

**Solution**: To narrow down the data frame, to make it easier to look at. Using `View()` for example.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Create a new data frame that shows the top 5 airports with the largest arrival delays from NYC in 2013.

**Solution**: 

```{r, purl=FALSE}
top_five <- flights %>%
  group_by(dest) %>%
  summarize(avg_delay = mean(arr_delay, na.rm = TRUE)) %>%
  arrange(desc(avg_delay)) %>%
  top_n(n = 5)
top_five
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Using the datasets included in the `nycflights13` package, compute the available seat miles for each airline sorted in descending order. After completing all the necessary data wrangling steps, the resulting data frame should have 16 rows (one for each airline) and 2 columns (airline name and available seat miles). Here are some hints:

1. **Crucial**: Unless you are very confident in what you are doing, it is worthwhile to not starting coding right away, but rather first sketch out on paper all the necessary data wrangling steps not using exact code, but rather high-level *pseudocode* that is informal yet detailed enough to articulate what you are doing. This way you won't confuse *what* you are trying to do (the algorithm) with *how* you are going to do it (writing `dplyr` code). 
1. Take a close look at all the datasets using the `View()` function: `flights`, `weather`, `planes`, `airports`, and `airlines` to identify which variables are necessary to compute available seat miles.
1. Figure \@ref(fig:reldiagram) above showing how the various datasets can be joined will also be useful. 
1. Consider the data wrangling verbs in Table \@ref(tab:wrangle-summary-table) as your toolbox!

**Solution**: Here are some examples of student-written [pseudocode](https://twitter.com/rudeboybert/status/964181298691629056). Based on our own pseudocode, let's first display the entire solution.

```{r, purl=FALSE}
flights %>%
  inner_join(planes, by = "tailnum") %>%
  select(carrier, seats, distance) %>%
  mutate(ASM = seats * distance) %>%
  group_by(carrier) %>%
  summarize(ASM = sum(ASM, na.rm = TRUE)) %>%
  arrange(desc(ASM))
```

Let's now break this down step-by-step. To compute the available seat miles for a given flight, we need the `distance` variable from the `flights` data frame and the `seats` variable from the `planes` data frame, necessitating a join by the key variable `tailnum` as illustrated in Figure \@ref(fig:reldiagram). To keep the resulting data frame easy to view, we'll `select()` only these two variables and `carrier`:

```{r, purl=FALSE}
flights %>%
  inner_join(planes, by = "tailnum") %>%
  select(carrier, seats, distance)
```

Now for each flight we can compute the available seat miles `ASM` by multiplying the number of seats by the distance via a `mutate()`:

```{r, purl=FALSE}
flights %>%
  inner_join(planes, by = "tailnum") %>%
  select(carrier, seats, distance) %>%
  # Added:
  mutate(ASM = seats * distance)
```

Next we want to sum the `ASM` for each carrier. We achieve this by first grouping by `carrier` and then summarizing using the `sum()` function:

```{r, purl=FALSE}
flights %>%
  inner_join(planes, by = "tailnum") %>%
  select(carrier, seats, distance) %>%
  mutate(ASM = seats * distance) %>%
  # Added:
  group_by(carrier) %>%
  summarize(ASM = sum(ASM))
```

However, because for certain carriers certain flights have missing `NA` values, the resulting table also returns `NA`'s. We can eliminate these by adding a `na.rm = TRUE` argument to `sum()`, telling R that we want to remove the `NA`'s in the sum. We saw this in Section \@ref(summarize):

```{r, purl=FALSE}
flights %>%
  inner_join(planes, by = "tailnum") %>%
  select(carrier, seats, distance) %>%
  mutate(ASM = seats * distance) %>%
  group_by(carrier) %>%
  # Modified:
  summarize(ASM = sum(ASM, na.rm = TRUE))
```

Finally, we `arrange()` the data in `desc()`ending order of `ASM`.

```{r, purl=FALSE}
flights %>%
  inner_join(planes, by = "tailnum") %>%
  select(carrier, seats, distance) %>%
  mutate(ASM = seats * distance) %>%
  group_by(carrier) %>%
  summarize(ASM = sum(ASM, na.rm = TRUE)) %>%
  # Added:
  arrange(desc(ASM))
```

While the above data frame is correct, the IATA `carrier` code is not always useful. For example, what carrier is `WN`? We can address this by joining with the `airlines` dataset using `carrier` is the key variable. While this step is not absolutely required, it goes a long way to making the table easier to make sense of. It is important to be empathetic with the ultimate consumers of your presented data!

```{r, purl=FALSE}
flights %>%
  inner_join(planes, by = "tailnum") %>%
  select(carrier, seats, distance) %>%
  mutate(ASM = seats * distance) %>%
  group_by(carrier) %>%
  summarize(ASM = sum(ASM, na.rm = TRUE)) %>%
  arrange(desc(ASM)) %>%
  # Added:
  inner_join(airlines, by = "carrier")
```



***



## Chapter 4 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 4
lc <- 0
```

```{r message=FALSE}
library(dplyr)
library(ggplot2)
library(readr)
library(tidyr)
library(nycflights13)
library(fivethirtyeight)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What are common characteristics of "tidy" datasets?

**Solution**: Rows correspond to observations, while columns correspond to variables.  

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What makes "tidy" datasets useful for organizing data?

**Solution**: Tidy datasets are an organized way of viewing data. This format is required for the `ggplot2` and `dplyr` packages for data visualization and wrangling. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Take a look the `airline_safety` data frame included in the `fivethirtyeight` data. Run the following:

```{r, eval=FALSE}
airline_safety
```

After reading the help file by running `?airline_safety`, we see that `airline_safety` is a data frame containing information on different airlines companies' safety records. This data was originally reported on the data journalism website FiveThirtyEight.com in Nate Silver's article ["Should Travelers Avoid Flying Airlines That Have Had Crashes in the Past?"](https://fivethirtyeight.com/features/should-travelers-avoid-flying-airlines-that-have-had-crashes-in-the-past/). Let's ignore the `incl_reg_subsidiaries` and `avail_seat_km_per_week` variables for simplicity:

```{r}
airline_safety_smaller <- airline_safety %>%
  select(-c(incl_reg_subsidiaries, avail_seat_km_per_week))
airline_safety_smaller
```

This data frame is not in "tidy" format. How would you convert this data frame to be in "tidy" format, in particular so that it has a variable `incident_type_years` indicating the incident type/year and a variable `count` of the counts?

**Solution**: 

This can been done using the `pivot_longer()` function from the `tidyr` package:

```{r}
airline_safety_smaller_tidy <- airline_safety_smaller %>%
  pivot_longer(
    names_to = "incident_type_years",
    values_to = "count",
    cols = -airline
  )
airline_safety_smaller_tidy
```

If you look at the resulting `airline_safety_smaller_tidy` data frame in the spreadsheet viewer, you'll see that the variable `incident_type_years` has 6 possible values: `"incidents_85_99", "fatal_accidents_85_99", "fatalities_85_99", 
"incidents_00_14", "fatal_accidents_00_14", "fatalities_00_14"` corresponding to the 6 columns of `airline_safety_smaller` we tidied. 

Note that prior to `tidyr` version 1.0.0 released to CRAN in September 2019, this could also have been done using the `gather()` function from the `tidyr` package:

```{r}
airline_safety_smaller_tidy <- airline_safety_smaller %>%
  gather(key = incident_type_years, value = count, -airline)
airline_safety_smaller_tidy
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  Convert the `dem_score` data frame into
a tidy data frame and assign the name of `dem_score_tidy` to the resulting long-formatted data frame.

**Solution**: Running the following in the console:

```{r}
dem_score_tidy <- dem_score %>%
  pivot_longer(
    names_to = "year", values_to = "democracy_score",
    cols = -country
  )
#  gather(key = year, value = democracy_score, - country)
```

Let's now compare the `dem_score` and `dem_score_tidy`. `dem_score` has democracy score information for each year in columns, whereas in `dem_score_tidy` there are explicit variables `year` and `democracy_score`. While both representations of the data contain the same information, we can only use `ggplot()` to create plots using the `dem_score_tidy` data frame.

```{r}
dem_score
dem_score_tidy
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  Read in the life expectancy data stored at <https://moderndive.com/data/le_mess.csv> and convert it to a tidy data frame. 

**Solution**: The code is similar:

```{r, eval=FALSE}
life_expectancy <- read_csv("https://moderndive.com/data/le_mess.csv")
life_expectancy_tidy <- life_expectancy %>%
  pivot_longer(
    names_to = "year",
    values_to = "life_expectancy",
    cols = -country
  )
#  gather(key = year, value = life_expectancy, -country)
```
```{r, echo=FALSE, purl=FALSE, message=FALSE, warning=FALSE}
life_expectancy <- read_csv("data/le_mess.csv")
life_expectancy_tidy <- life_expectancy %>%
  pivot_longer(
    names_to = "year",
    values_to = "life_expectancy",
    cols = -country
  )
#  gather(key = year, value = life_expectancy, -country)
```

We observe the same construct structure with respect to `year` in `life_expectancy` vs `life_expectancy_tidy` as we did in `dem_score` vs `dem_score_tidy`:

```{r}
life_expectancy
life_expectancy_tidy
```



***



## Chapter 5 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 5
lc <- 0
```

```{r message=FALSE}
library(tidyverse)
library(moderndive)
library(skimr)
library(gapminder)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Conduct a new exploratory data analysis with the same outcome variable $y$ being `score` but with `age` as the new explanatory variable $x$. Remember, this involves three things:

a) Looking at the raw data values.
a) Computing summary statistics.
a) Creating data visualizations.

What can you say about the relationship between age and teaching scores based on this exploration?

**Solution**: 

- Looking at the raw data values: 
```{r}
glimpse(evals_ch5)
```

- Computing summary statistics:
```{r, eval=FALSE}
skim_with(numeric = list(hist = NULL), integer = list(hist = NULL))
```
```{r, include=FALSE}
skim_with(numeric = list(hist = NULL), integer = list(hist = NULL))
```

```{r, eval=FALSE}
evals_ch5 %>%
  select(score, age) %>%
  skim()
```

(Note that for formatting purposes, the inline histogram that is usually `printed with skim()` has been removed. This can be done by running `skim_with(numeric = list(hist = NULL), integer = list(hist = NULL))` prior to using the `skim()` function as well.)

- Creating data visualizations:
```{r}
ggplot(evals_ch5, aes(x = age, y = score)) +
  geom_point() +
  labs(
    x = "Age", y = "Teaching Score",
    title = "Scatterplot of relationship of teaching score and age"
  )
```

<!--
TODO: Albert needs to double check interpretation:
-->

Based on the scatterplot visualization, there seem to have a weak negative relationship between age and teaching score. As age increases, the teaching score see, to decrease slightly. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Fit a new simple linear regression using `lm(score ~ age, data = evals_ch5)` where `age` is the new explanatory variable $x$. Get information about the "best-fitting" line from the regression table by applying the `get_regression_table()` function. How do the regression results match up with the results from your earlier exploratory data analysis?

**Solution**: 
```{r}
# Fit regression model:
score_age_model <- lm(score ~ age, data = evals_ch5)
# Get regression table:
get_regression_table(score_age_model)
```
$$
\begin{aligned}
\widehat{y} &= b_0 + b_1 \cdot x\\
\widehat{\text{score}} &= b_0 + b_{\text{age}} \cdot\text{age}\\
&= 4.462 - 0.006\cdot\text{age}
\end{aligned}
$$

<!--
TODO: Albert will verify Starry's interpretation:
-->

For every increase of 1 unit in `age`, there is an *associated* decrease of, *on average*, 0.006 units of `score`. It matches with the results from our earlier exploratory data analysis. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Generate a data frame of the residuals of the model where you used `age` as the explanatory $x$ variable.

**Solution**: 
```{r}
score_age_regression_points <- get_regression_points(score_age_model)
score_age_regression_points
```
**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Conduct a new exploratory data analysis with the same explanatory variable $x$ being `continent` but with `gdpPercap` as the new outcome variable $y$. Remember, this involves three things:

1. Most crucially: Looking at the raw data values.
1. Computing summary statistics, such as means, medians, and interquartile ranges.
1. Creating data visualizations.

What can you say about the differences in GDP per capita between continents based on this exploration?

**Solution**: 

- Looking at the raw data values: 
```{r}
glimpse(gapminder2007)
```

- Computing summary statistics, such as means, medians, and interquartile ranges:
```{r eval=FALSE}
gapminder2007 %>%
  select(gdpPercap, continent) %>%
  skim()
```

- Creating data visualizations:
```{r}
ggplot(gapminder2007, aes(x = continent, y = gdpPercap)) +
  geom_boxplot() +
  labs(
    x = "Continent", y = "GPD per capita",
    title = "GDP by continent"
  )
```

<!--
TODO: Albert needs to double check interpretation:
-->

Based on this exploration, it seems that GDP's are very different among different continents, which means that continent might be a statistically significant predictor for an area's GDP. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Fit a new linear regression using `lm(gdpPercap ~ continent, data = gapminder2007)` where `gdpPercap` is the new outcome variable $y$. Get information about the "best-fitting" line from the regression table by applying the `get_regression_table()` function. How do the regression results match up with the results from your previous exploratory data analysis?

**Solution**: 
```{r}
# Fit regression model:
gdp_model <- lm(gdpPercap ~ continent, data = gapminder2007)
# Get regression table:
get_regression_table(gdp_model)
```
$$
\begin{aligned}
\widehat{y} = \widehat{\text{gdpPercap}} &= b_0 + b_{\text{Amer}}\cdot\mathbb{1}_{\mbox{Amer}}(x) + b_{\text{Asia}}\cdot\mathbb{1}_{\mbox{Asia}}(x) + \\
& \qquad b_{\text{Euro}}\cdot\mathbb{1}_{\mbox{Euro}}(x) + b_{\text{Ocean}}\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\
&= 3089 + 7914\cdot\mathbb{1}_{\mbox{Amer}}(x) + 9384\cdot\mathbb{1}_{\mbox{Asia}}(x) + \\
& \qquad 21965\cdot\mathbb{1}_{\mbox{Euro}}(x) + 26721\cdot\mathbb{1}_{\mbox{Ocean}}(x)
\end{aligned}
$$

<!--
TODO: Albert will double check Starry's interpretation:
-->

In our previous exploratory data analysis, it seemed that continent is a statistically significant predictor for an area's GDP. Here, by fit a new linear regression using `lm(gdpPercap ~ continent, data = gapminder2007)` where `gdpPercap` is the new outcome variable $y$, we are able to write an equation to predict `gdpPercap` using the continent as statistically significant predictors. Therefore, the regression results matches with the results from your previous exploratory data analysis. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Using either the sorting functionality of RStudio's spreadsheet viewer or using the data wrangling tools you learned in Chapter \@ref(wrangling), identify the five countries with the five smallest (most negative) residuals? What do these negative residuals say about their life expectancy relative to their continents?

**Solution**: 
Using the sorting functionality of RStudio's spreadsheet viewer, we can identify that the five countries with the five smallest (most negative) residuals are: Afghanistan, Swaziland, Mozambique, Haiti, and Zambia. 

These negative residuals indicate that these data points have the biggest negative deviations from their group means. This means that these five countries' average life expectancies are the lowest comparing to their respective continents' average life expectancies. For example, the residual for Afghanistan is $-26.900$ and it is the smallest residual. This means that the average life expectancy of Afghanistan is $26.900$ years lower than the average life expectancy of its continent, Asia. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Repeat this process, but identify the five countries with the five largest (most positive) residuals. What do these positive residuals say about their life expectancy relative to their continents?

**Solution**: 
Using either the sorting functionality of RStudio's spreadsheet viewer, we can identify that the five countries with the five largest (most positive) residuals are: Reunion, Libya, Tunisia, Mauritius, and Algeria. 

These positive residuals indicate that the data points are above the regression line with the longest distance. This means that these five countries' average life expectancies are the highest comparing to their respective continents' average life expectancies. For example, the residual for Reunion is $21.636$ and it is the largest residual. This means that the average life expectancy of Reunion is $21.636$ years higher than the average life expectancy of its continent, Africa. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Note in the following plot there are 3 points marked with dots along with:

* The "best" fitting solid regression line `r if_else(is_latex_output(), "", "in blue")`
* An arbitrarily chosen dotted `r if_else(is_latex_output(), "", "red")` line 
* Another arbitrarily chosen dashed `r if_else(is_latex_output(), "", "green")` line

```{r, fig.cap="Regression line and two others.", out.width="80%", echo=FALSE}
example <- tibble(
  x = c(0, 0.5, 1),
  y = c(2, 1, 3)
)

ggplot(example, aes(x = x, y = y)) +
  geom_smooth(method = "lm", se = FALSE, fullrange = TRUE) +
  geom_hline(yintercept = 2.5, col = "red", linetype = "dotted", size = 1) +
  geom_abline(
    intercept = 2, slope = -1, col = "forestgreen",
    linetype = "dashed", size = 1
  ) +
  geom_point(size = 4)
```
Compute the sum of squared residuals by hand for each line and show that of these three lines, the regression line `r if_else(is_latex_output(), "", "in blue")` has the smallest value.

**Solution**:  

* The "best" fitting solid regression line `r if_else(is_latex_output(), "", "in blue")`: 

$$
\sum_{i=1}^{n}(y_i - \widehat{y}_i)^2 = (2.0-1.5)^2+(1.0-2.0)^2+(3.0-2.5)^2 = 1.5
$$

* An arbitrarily chosen dotted `r if_else(is_latex_output(), "", "red")` line: 

$$
\sum_{i=1}^{n}(y_i - \widehat{y}_i)^2 = (2.0-2.5)^2+(1.00-2.5)^2+(3.0-2.5)^2 = 2.75
$$

* Another arbitrarily chosen dashed `r if_else(is_latex_output(), "", "green")` line:

$$
\sum_{i=1}^{n}(y_i - \widehat{y}_i)^2 = (2.0-2.0)^2+(1.0-1.5)^2+(3.0-1.0)^2 = 4.25
$$

As calculated, $1.5 < 2.75 < 4.25$. Therefore, we show that the regression line in blue has the smallest value of the residual sum of squares.



***



## Chapter 6 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 6
lc <- 0
```

```{r message=FALSE}
library(tidyverse)
library(moderndive)
library(skimr)
library(ISLR)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Compute the observed values, fitted values, and residuals not for the interaction model as we just did, but rather for the parallel slopes model we saved in `score_model_parallel_slopes`.

**Solution**: 
```{r}
regression_points_parallel <- get_regression_points(score_model_parallel_slopes)
regression_points_parallel
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Conduct a new exploratory data analysis with the same outcome variable $y$ being `debt` but with `credit_rating` and `age` as the new explanatory variables $x_1$ and $x_2$. Remember, this involves three things:

1. Most crucially: Looking at the raw data values.
1. Computing summary statistics, such as means, medians, and interquartile ranges.
1. Creating data visualizations.

What can you say about the relationship between a credit card holder's debt and their credit rating and age?

**Solution**: 

* Most crucially: Looking at the raw data values.
```{r}
credit_ch6 %>%
  select(debt, credit_rating, age) %>%
  head()
```

* Computing summary statistics, such as means, medians, and interquartile ranges.
```{r, eval=FALSE}
skim_with(numeric = list(hist = NULL), integer = list(hist = NULL))
```
```{r, include=FALSE}
skim_with(numeric = list(hist = NULL), integer = list(hist = NULL))
```
```{r, eval=FALSE}
credit_ch6 %>%
  select(debt, credit_rating, age) %>%
  skim()
```
* Creating data visualizations.
```{r}
ggplot(credit_ch6, aes(x = credit_rating, y = debt)) +
  geom_point() +
  labs(
    x = "Credit rating", y = "Credit card debt (in $)",
    title = "Debt and credit rating"
  ) +
  geom_smooth(method = "lm", se = FALSE)

ggplot(credit_ch6, aes(x = age, y = debt)) +
  geom_point() +
  labs(
    x = "Age (in year)", y = "Credit card debt (in $)",
    title = "Debt and age"
  ) +
  geom_smooth(method = "lm", se = FALSE)
```
It seems that there is a positive relationship between one's credit rating and their debt, and very little relationship between one's age and their debt. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Fit a new simple linear regression using `lm(debt ~ credit_rating + age, data = credit_ch6)` where `credit_rating` and `age` are the new numerical explanatory variables $x_1$ and $x_2$. Get information about the "best-fitting" regression plane from the regression table by applying the `get_regression_table()` function. How do the regression results match up with the results from your previous exploratory data analysis? 
```{r}
# Fit regression model:
debt_model_2 <- lm(debt ~ credit_rating + age, data = credit_ch6)
# Get regression table:
get_regression_table(debt_model_2)
```
The coefficients for both new numerical explanatory variables $x_1$ and $x_2$, `credit_rating` and `age`, are $2.59$ and $-2.35$ respectively, which means that `debt` and `credit_rating` are positively correlated, and `debt` and `age` are negatively correlated. This matches up with the results from your previous exploratory data analysis. 



***



## Chapter 7 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 7
lc <- 0
```
```{r message=FALSE}
library(ggplot2)
library(dplyr)
library(moderndive)
library(gapminder)
library(skimr)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why was it important to mix the bowl before we sampled the balls?
  
  **Solution**:  
  
  So that we make sure the sampled balls are randomized. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why is it that our 33 groups of friends did not all have the same numbers of balls that were red out of 50, and hence different proportions red?
  
  **Solution**:  
  
  Because not all pairs have the same portion of the population of the balls, so each pair has a different sampled balls with different color compositions. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why couldn't we study the effects of sampling variation when we used the virtual shovel only once? Why did we need to take more than one virtual sample (in our case 33 virtual samples)?

**Solution**: 

If we use the virtual shovel only once, we only get one sample of the population. We need to take more than one virtual sample to get a range of proportions.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why did we not take 1000 "tactile" samples of 50 balls by hand?

**Solution**: 

That would be way too much repeated work. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Looking at Figure \@ref(fig:samplingdistribution-virtual-1000), would you say that sampling 50 balls where 30% of them were red is likely or not? What about sampling 50 balls where 10% of them were red?

**Solution**: 

According to the Figure, less than 150 out of the 1000 counts were 30% red. So I would say that sampling 50 balls where 30% of them were red is not very likely. Almost no count was only 10% red, so sampling 50 balls where 10% of them were red is extremely unlikely. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** In Figure \@ref(fig:comparing-sampling-distributions), we used shovels to take 1000 samples each, computed the resulting 1000 proportions of the shovel's balls that were red, and then visualized the distribution of these 1000 proportions in a histogram. We did this for shovels with 25, 50, and 100 slots in them. As the size of the shovels increased, the histograms got narrower. In other words, as the size of the shovels increased from 25 to 50 to 100, did the 1000 proportions

- A. vary less,
- B. vary by the same amount, or
- C. vary more?
  
  **Solution**: 
  
  A. As the histograms got narrower, the 1000 proportions varied less. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What summary statistic did we use to quantify how much the 1000 proportions red varied?
  
  - A. The inter-quartile range
- B. The standard deviation
- C. The range: the largest value minus the smallest.

**Solution**: 
  
  B. The standard deviation is used to quantify how much a set of data varies. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** In the case of our bowl activity, what is the *population parameter*? Do we know its value?
  
  **Solution**: 
  
  The *population parameter* in the case of our bowl activity is the population proportion of the red balls in the bowl. Unless we know the exact number of red balls in the bowl, we won't know the value of this population proportion. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What would performing a census in our bowl activity correspond to? Why did we not perform a census?
  
  **Solution**: 
  
  Performing a census in our bowl activity correspond to counting the total number of red balls in all balls, We did not perform a census because it would be too much repetitive work and it is unnecessary. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What purpose do *point estimates* serve in general? What is the name of the point estimate specific to our bowl activity? What is its mathematical notation?
  
  **Solution**: 
  
  *Point estimates* serve to *estimate* an unknown population parameter in the sample. In our bowl activity, our point estimate is the *sample proportion*: the proportion of the shovel's balls that are red. We mathematically denote the sample proportion using $\widehat{p}$.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** How did we ensure that our tactile samples using the shovel were random?

**Solution**: 

We virtually shuffle the sample each time. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why is it important that sampling be done *at random*?

**Solution**: 

So that we get different samples each time to estimate the total population. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What are we *inferring* about the bowl based on the samples using the shovel?

**Solution**: 

We are *inferring* that the samples are representing the total population in the ball. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What purpose did the *sampling distributions* serve?

**Solution**: 

Using the sampling distributions, for a given sample size $n$, we can make statements about what values we can typically expect. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What does the *standard error* of the sample proportion $\widehat{p}$ quantify? 

**Solution**: 

Standard errors quantify the effect of sampling variation induced on our estimates. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** The table that follows is a version of Table \@ref(tab:comparing-n-2) matching sample sizes $n$ to different *standard errors* of the sample proportion $\widehat{p}$, but with the rows randomly re-ordered and the sample sizes removed. Fill in the table by matching the correct sample sizes to the correct standard errors.

```{r comparing-n-3-solution, eval=TRUE, echo=FALSE}
set.seed(76)
comparing_n_table <- virtual_prop %>%
  group_by(n) %>%
  summarize(sd = sd(prop_red)) %>%
  mutate(
    n = str_c("n = ")
  ) %>%
  rename(`Sample size` = n, `Standard error of $\\widehat{p}$` = sd) %>%
  sample_frac(1)

comparing_n_table %>%
  kable(
    digits = 3,
    caption = "Three standard errors of the sample proportion based on n = 25, 50, 100",
    booktabs = TRUE,
    escape = FALSE,
    linesep = ""
  ) %>%
  kable_styling(
    font_size = ifelse(is_latex_output(), 10, 16),
    latex_options = c("hold_position")
  )
```

**Solution**:

$n$ = $25$, $100$, $50$ respectively. 

For the following four learning checks, let the *estimate* be the sample proportion $\widehat{p}$: the proportion of a shovel's balls that were red. It estimates the population proportion $p$: the proportion of the bowl's balls that were red.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What is the difference between an *accurate* estimate and a *precise* estimate? 

**Solution**:

An *accurate* estimate gives an estimate that is close to, but not necessary the exact, actual value. A *precise* estimate gives the exact actual value. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** How do we ensure that an estimate is *accurate*? How do we ensure that an estimate is *precise*?

To ensure that an estimate is *accurate*, we need to have a reasonable range of estimate, and make sure that the estimate is reasonably close to the actual value. To ensure that an estimate is *precise*, we need to make sure the estimate is equivalent to the actual value. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** In a real-life situation, we would not take 1000 different samples to infer about a population, but rather only one. Then, what was the purpose of our exercises where we took 1000 different samples?

**Solution**:

To get a narrower range of the estimates. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Figure \@ref(fig:accuracy-vs-precision) with the targets shows four combinations of "accurate versus precise" estimates. Draw four corresponding *sampling distributions* of the sample proportion $\widehat{p}$, like the one in the left-most plot in Figure \@ref(fig:comparing-sampling-distributions-3).

**Solution**:
```{r, echo=FALSE}
xlim <- c(-6, 6)
p1 <- ggplot(data = data.frame(x = xlim), aes(x)) +
  stat_function(fun = dnorm, n = 101, args = list(mean = 0, sd = .5)) +
  ylab("") +
  # scale_y_continuous(breaks = NULL)+
  geom_vline(xintercept = 0, linetype = "dotted") +
  coord_cartesian(xlim = xlim, ylim = c(0, 0.8)) +
  ggtitle("High precision, high accuracy")

p2 <- ggplot(data = data.frame(x = xlim), aes(x)) +
  stat_function(fun = dnorm, n = 101, args = list(mean = 2.5, sd = .5)) +
  ylab("") +
  # scale_y_continuous(breaks = NULL)+
  geom_vline(xintercept = 0, linetype = "dotted") +
  coord_cartesian(xlim = xlim, ylim = c(0, 0.8)) +
  ggtitle("High precision, low accuracy")

p3 <- ggplot(data = data.frame(x = xlim), aes(x)) +
  stat_function(fun = dnorm, n = 101, args = list(mean = 0, sd = 1.5)) +
  ylab("") +
  # scale_y_continuous(breaks = NULL)+
  geom_vline(xintercept = 0, linetype = "dotted") +
  coord_cartesian(xlim = xlim, ylim = c(0, 0.8)) +
  ggtitle("Low precision, high accuracy")

p4 <- ggplot(data = data.frame(x = xlim), aes(x)) +
  stat_function(fun = dnorm, n = 101, args = list(mean = 2.5, sd = 1.5)) +
  ylab("") +
  # scale_y_continuous(breaks = NULL)+
  geom_vline(xintercept = 0, linetype = "dotted") +
  coord_cartesian(xlim = xlim, ylim = c(0, 0.8)) +
  ggtitle("Low precision, low accuracy")

library(patchwork)
(p1 | p2) / (p3 | p4)
```

Comment on the representativeness of the following *sampling methodologies*:

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** The Royal Air Force wants to study how resistant all their airplanes are to bullets. They study the bullet holes on all the airplanes on the tarmac after an air battle against the Luftwaffe (German Air Force).

**Solution**:

The airplanes on the tarmac after an air battle against the Luftwaffe is not a good representation of all airplanes, because the airplanes which were attacked in less resistant areas did not make it back to the tarmac. This is called *survival bias*. Survivor's bias or survival bias is the logical error of concentrating on the people or things that made it past some selection process and overlooking those that did not, typically because of their lack of visibility. This can lead to false conclusions in several different ways. It is a form of selection bias.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Imagine it is 1993, a time when almost all households had landlines. You want to know the average number of people in each household in your city. You randomly pick out 500 phone numbers from the phone book and conduct a phone survey.

**Solution**:

This is not a good representation, because: (1) adults are more likely to pickup phone calls; (2) households with more people are more likely to have people to be available to pickup phone calls; (3) we are not certain whether all households are in the phone book. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** You want to know the prevalence of illegal downloading of TV shows among students at a local college.  You get the emails of 100 randomly chosen students and ask them, "How many times did you download a pirated TV show last week?".

**Solution**:

This is not a good representation, because it is very likely that students will lie in this survey to stay out of trouble. So we may not get honest data. This is called *volunteer bias*: systematic error due to differences between those who choose to participate in studies and those who do not.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** A local college administrator wants to know the average income of all graduates in the last 10 years. So they get the records of five randomly chosen graduates, contact them, and obtain their answers. 

**Solution**:

This is not a good representation, because the sample size is too small. The sample is representative but not precise. 



***



## Chapter 8 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 8
lc <- 0
```

```{r message=FALSE}
library(tidyverse)
library(moderndive)
library(infer)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What is the chief difference between a bootstrap distribution and a sampling distribution?

**Solution**:

A bootstrap sample is a smaller sample that is "bootstrapped" from a larger sample. Bootstrapping is a type of resampling where large numbers of smaller samples of the same size are repeatedly drawn, with replacement, from a single original sample.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Looking at the bootstrap distribution for the sample mean in Figure \@ref(fig:one-thousand-sample-means), between what two values would you say *most* values lie?

**Solution**:

*Most* values lie in 1990 and 2000. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What condition about the bootstrap distribution must be met for us to be able to construct confidence intervals using the standard error method?

**Solution**: 

We can only use the standard error rule when the bootstrap distribution is roughly normally distributed. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Say we wanted to construct a 68% confidence interval instead of a 95% confidence interval for $\mu$. Describe what changes are needed to make this happen. Hint: we suggest you look at Appendix \@ref(appendix-normal-curve) on the normal distribution.

**Solution**: 

Thus, using our 68% rule of thumb about normal distributions from Appendix \@ref(appendix-normal-curve), we can use the following formula to determine the lower and upper endpoints of a 95% confidence interval for $\mu$:

$$\overline{x} \pm 1 \cdot SE = (\overline{x} - 1 \cdot SE, \overline{x} + 1 \cdot SE)$$

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Construct a 95% confidence interval for the *median* year of minting of *all* US pennies? Use the percentile method and, if appropriate, then use the standard-error method.

**Solution**: 

Using the percentile method:

```{r}
bootstrap_distribution <- pennies_sample %>%
  specify(response = year) %>%
  generate(reps = 1000) %>%
  calculate(stat = "median")
percentile_ci <- bootstrap_distribution %>%
  get_confidence_interval(level = 0.95, type = "percentile")
percentile_ci
```

The standard-error method is not appropriate, because the bootstrap distribution is not bell-shaped:

```{r}
visualize(bootstrap_distribution)
```



***



## Chapter 9 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 9
lc <- 0
```

```{r message=FALSE}
library(tidyverse)
library(infer)
library(moderndive)
library(nycflights13)
library(ggplot2movies)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why does the following code produce an error? In other words, what about the response and predictor variables make this not a possible computation with the `infer` package?

```{r, eval=FALSE}
library(moderndive)
library(infer)
null_distribution_mean <- promotions %>%
  specify(formula = decision ~ gender, success = "promoted") %>% 
  hypothesize(null = "independence") %>% 
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "diff in means", order = c("male", "female"))
```

**Solution**: Because the `decision` outcome variable is categorical and we can't compute a difference in means of a categorical variable. We can however compute the `diff in props` instead. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why are we relatively confident that the distributions of the sample proportions will be good approximations of the population distributions of promotion proportions for the two genders?

**Solution**:  

Because the sample is representative of the population. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Using the definition of _p-value_, write in words what the $p$-value represents for the hypothesis test comparing the promotion rates for males and females.

**Solution**:  

The $p$-value represents for the likelihood that the true mean for the promotion rates for males and females in the population is the same. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Describe in a paragraph how we used Allen Downey's diagram to conclude if a statistical difference existed between the promotion rate of males and females using this study.

**Solution**:  

We use the `promotions` dataset as the input for test statistic. The $H_0$ model is "there is no difference between promotion rates of males and females", and with the p-value from `infer` commands, we reject the $H_0$ model and conclude that there is a statistical difference existed between the promotion rate of males and females. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  What is wrong about saying, "The defendant is innocent." based on the US system of criminal trials?

**Solution**:  

Failing to prove the defendant is guilty does not equal to proving that the defendant is innocent. There will always be the possibility of making errors in the trial. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  What is the purpose of hypothesis testing?

**Solution**:  

The purpose of hypothesis testing is to determine whether there is enough statistical evidence in favor of a certain belief, or hypothesis, about a parameter. (source: https://personal.utdallas.edu/~scniu/OPRE-6301/documents/Hypothesis_Testing.pdf)

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  What are some flaws with hypothesis testing?  How could we alleviate them?

**Solution**:  

The p-value's 0.05 threshold can be misleading researchers to conduct multiple bootstrap tests to get a smaller p-value, therefore validating their statistical results. This threshold is relatively arbitrary (if a p-value is 0.051, does it mean there is no statistical significance?), and trusting it too much may lead to imprecise conclusions. To alleviate this problem, keep in mind that having a smaller p-value can be the result of a "lucky" sampling that is not truly representative, and do multiple bootstrap samplings for hypothesis testing before concluding. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`**  Consider two $\alpha$ significance levels of 0.1 and 0.01. Of the two, which would lead to a more *liberal* hypothesis testing procedure? In other words, one that will, all things being equal, lead to more rejections of the null hypothesis $H_0$.

**Solution**:  

The greater $\alpha$ of 0.1 will lead to a more liberal hypothesis testing procedure, because the required p-value to reject the null hypothesis $H_0$ can be greater. There is a higher chance that the p-value will be less than $\alpha$.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Conduct the same analysis comparing action movies versus romantic movies using the median rating instead of the mean rating. What was different and what was the same? 

**Solution**:  

```{r, eval=FALSE}
library(moderndive)
library(infer)

# In calculate() step replace "diff in means" with "diff in medians"
null_distribution_movies_median <- movies_sample %>% 
  specify(formula = rating ~ genre) %>% 
  hypothesize(null = "independence") %>% 
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "diff in medians", order = c("Action", "Romance"))

# compute observed "diff in medians"
obs_diff_medians <- movies_sample %>% 
  specify(formula = rating ~ genre) %>% 
  calculate(stat = "diff in medians", order = c("Action", "Romance"))
obs_diff_medians

# Visualize p-value. Observing this difference in medians under H0
# is very unlikely! Suggesting H0 is false, similarly to when we used
# "diff in means" as the test statistic. 
visualize(null_distribution_movies_median, bins = 10) + 
  shade_p_value(obs_stat = obs_diff_medians, direction = "both")

# p-value is very small, just like when we used "diff in means"
# as the test statistic. 
null_distribution_movies_median %>% 
  get_p_value(obs_stat = obs_diff_medians, direction = "both")
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What conclusions can you make from viewing the faceted histogram looking at `rating` versus `genre` that you couldn't see when looking at the boxplot?

**Solution**:  

From the faceted histogram, we can also see the comparison of `rating` versus `genre` over each rating range, but we cannot conclude them from the boxplot. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Describe in a paragraph how we used Allen Downey's diagram to conclude if a statistical difference existed between mean movie ratings for action and romance movies.

**Solution**:  

We use the `movies_sample` dataset as the input for test statistic. The $H_0$ model is "there is no statistical difference existed between mean movie ratings for action and romance movies", and with the p-value from `infer` commands, we fail to reject the $H_0$ model and conclude that there is insufficient evidence to conclude that a statistical difference existed between mean movie ratings for action and romance movies.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Why are we relatively confident that the distributions of the sample ratings will be good approximations of the population distributions of ratings for the two genres?

**Solution**:  

Because the sample is representative of the population. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Using the definition of $p$-value, write in words what the $p$-value represents for the hypothesis test comparing the mean rating of romance to action movies.

**Solution**:  

The $p$-value represent the probability that the difference  between mean movie ratings for action and romance movies in the sample is natural, i.e., the probability that there is no statistical difference between mean movie ratings for action and romance movies in the population. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What is the value of the $p$-value for the hypothesis test comparing the mean rating of romance to action movies?

**Solution**: 

The $p$-value here is $0.004$.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Test your data wrangling knowledge and EDA skills:

- Use `dplyr` and `tidyr` to create the necessary data frame focused on only action and romance movies (but not both) from the `movies` data frame in the `ggplot2movies` package.
- Make a boxplot and a faceted histogram of this population data comparing ratings of action and romance movies from IMDb.
- Discuss how these plots compare to the similar plots produced for the `movies_sample` data.

**Solution**: 

- Use `dplyr` and `tidyr` to create the necessary data frame focused on only action and romance movies (but not both) from the `movies` data frame in the `ggplot2movies` package. 

```{r}
action_romance <- movies %>%
  select(title, year, rating, Action, Romance) %>%
  filter((Action == 1 | Romance == 1) & !(Action == 1 & Romance == 1)) %>% 
  # filter(Action + Romance == 1) %>%
  mutate(genre = ifelse(Action == 1, "Action", "Romance"))
```

Note above:

1. In the `filter()` command that `(Action == 1 | Romance == 1)` picks out the movies that are either action or romance while the `!(Action == 1 & Romance == 1)` leaves out the movies that are classified as both. Alternatively, as only one of `Action` and `Romance` can be 1, the sum of these two variables for those movies of either genre will be 1. Therefore, you can also use the commented-out filter to pick out the movies.
2. In the `mutate()` we use `ifelse()` to create a new variable `genre` whose values are either `"Action"` or `"Romance"`.

- Make a boxplot and a faceted histogram of this population data comparing ratings of action and romance movies from IMDb.

```{r}
ggplot(action_romance, aes(rating)) +
  geom_histogram() +
  facet_wrap(~genre)
```



***



## Chapter 10 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 10
lc <- 0
```

```{r message=FALSE}
library(tidyverse)
library(moderndive)
library(infer)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Continuing with our regression using `age` as the explanatory variable and teaching `score` as the outcome variable.

- Use the `get_regression_points()` function to get the observed values, fitted values, and residuals for all `r nrow(evals_ch5)` instructors. 

```{r, eval=TRUE, echo=TRUE}
evals_ch5 <- evals %>%
  select(ID, score, bty_avg, age)

# Fit regression model:
score_age_model <- lm(score ~ age, data = evals_ch5)
# Get regression points:
regression_points <- get_regression_points(score_age_model)
regression_points
```

- Perform a residual analysis and look for any systematic patterns in the residuals. Ideally, there should be little to no pattern but comment on what you find here.

The first condition is that the relationship between the outcome variable $y$ and the explanatory variable $x$ must be **L**inear. 

```{r non-linear2, fig.cap="Example of a clearly non-linear relationship.", echo=FALSE, fig.height=3.3}
set.seed(76)
evals_ch5 %>%
  ggplot(aes(x = age, y = score)) +
  geom_point() +
  labs(x = "Age", y = "Teaching Score") +
  geom_smooth(method = "lm", se = FALSE) +
  expand_limits(y = 10)
```
The relationship between score and age does not seem to be linear. 

The second condition is that the residuals must be **I**ndependent. In other words, the different observations in our data must be independent of one another. As explained in \@ref(independence-of-residuals),  "we say there exists *dependence* between observations".

The third condition is that the residuals should follow a **N**ormal distribution. 
```{r model1residualshist2, echo=FALSE, warning=FALSE, fig.cap="Histogram of residuals."}
ggplot(regression_points, aes(x = residual)) +
  geom_histogram(binwidth = 0.25, color = "white") +
  labs(x = "Residual")
```

The fourth and final condition is that the residuals should exhibit **E**qual variance across all values of the explanatory variable $x$. In other words, the value and spread of the residuals should not depend on the value of the explanatory variable $x$.  
```{r numxplot62, echo=FALSE, warning=FALSE, fig.cap="Plot of residuals over beauty score."}
ggplot(regression_points, aes(x = age, y = residual)) +
  geom_point() +
  labs(x = "Beauty Score", y = "Residual") +
  geom_hline(yintercept = 0, col = "blue", size = 1)
```
This plot seems to fit equality of variance. 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Repeat the inference but this time for the correlation coefficient instead of the slope. Note the implementation of `stat = "correlation"` in the `calculate()` function of the `infer` package.
<!--
TODO: Clarify what the deliverable for this LC should be.
-->
```{r eval=FALSE}
bootstrap_distn_slope <- evals_ch5 %>%
  specify(formula = score ~ bty_avg) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "correlation")

if (!file.exists("rds/bootstrap_distn_slope.rds")) {
  set.seed(76)
  bootstrap_distn_slope <- evals %>%
    specify(score ~ bty_avg) %>%
    generate(reps = 1000, type = "bootstrap") %>%
    calculate(stat = "slope")
  saveRDS(
    object = bootstrap_distn_slope,
    "rds/bootstrap_distn_slope.rds"
  )
} else {
  bootstrap_distn_slope <- readRDS("rds/bootstrap_distn_slope.rds")
}

observed_slope <- evals %>%
  specify(score ~ bty_avg) %>%
  calculate(stat = "correlation")
observed_slope
```



***



## Chapter 11 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- 11
lc <- 0
```

```{r message=FALSE}
library(tidyverse)
library(moderndive)
library(skimr)
library(fivethirtyeight)
``` 

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** Repeat the regression modeling in Subsection \@ref(house-prices-regression) and the prediction making you just did on the house of condition 5 and size 1900 square feet in Subsection \@ref(house-prices-making-predictions), but using the parallel slopes model you visualized in Figure \@ref(fig:house-price-parallel-slopes). Show that it's `r 10^5.720317 %>% dollar()`!

```{r, eval=FALSE}
house_prices <- house_prices %>%
  mutate(
    log10_price = log10(price),
    log10_size = log10(sqft_living)
  )
# Fit regression model:
price_interaction <- lm(log10_price ~ log10_size + condition,
  data = house_prices
)

# Get regression table:
get_regression_table(price_interaction)

10^(2.88 + 0.096 + 0.837 * log10(1900))
```

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What date between 1994 and 2003 has the fewest number of births in the US? What story could you tell about why this is the case?

```{r}
US_births_1994_2003 %>%
  arrange(births)
```

The dates with the fewest number of births in the US was 12/25 of the years of 2001, 2000, 2003, 2002, and 1999. Because it is Christmas Day and hospitals don't generally induce labor on that day. 



***



## Appendix A2 Solutions

```{r, include=FALSE, purl=FALSE}
chap <- "A2"
lc <- 0
```


**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What proportion of the area under the normal curve is less than 3? Greater than 12? Between 0 and 12?

* Less than 3: 3 is one standard deviation less than the mean of 6, since $\frac{3 - \mu}{\sigma} = \frac{3-6}{3} = -1$. Thus we compute the area to the left of z = -1 in Figure \@ref(fig:normal-rule-of-thumb): 0.15% + 2.35% + 13.5% = 16%.
* Greater than 12: 12 is two standard deviations greater than the mean of 6, since $\frac{12 - \mu}{\sigma} = \frac{12-6}{3} = +2$, we compute the area to the right of z = +2 in Figure \@ref(fig:normal-rule-of-thumb): 2.35% + 0.15% = 2.5%.
* Between 0 and 12: 0 is two standard deviations less than the mean of 6, since $\frac{0 - \mu}{\sigma} = \frac{0-6}{3} = -2$. Thus we compute the area between z = -2 and z = +2 in Figure \@ref(fig:normal-rule-of-thumb): 13.5% + 34% + 34% + 13.5% = 95%.

**`r paste0("(LC", chap, ".", (lc <- lc + 1), ")")`** What is the 2.5th percentile of the area under the normal curve? The 97.5th percentile? The 100th percentile?

* 2.5th percentile: Starting from the left of Figure \@ref(fig:normal-rule-of-thumb), since 0.15% + 2.35% = 2.5%, the 2.5th percentile is z = - 2. However, this is in standard units. Thus we need the value in the normal distribution that is two standard deviations lower than the mean: $\mu - 2 \cdot\sigma = 6 - 2 \cdot 3 = 0$.
* 97.5th percentile: Starting from the left of Figure \@ref(fig:normal-rule-of-thumb), since 0.15% + 2.35% + 13.5% + 34% + 34% + 13.5% = 97.5%, the 97.5th percentile is z = +2. However, this is in standard units. Thus we need the value in the normal distribution that is two standard deviations higher than the mean: $\mu + 2 \cdot\sigma = 6 + 2 \cdot 3 = 12$.
* 100th percentile: $+\infty$. In other words, 100% of values will be less than positive infinity.