week11_tutorial.qmd

---
title: "ETC1010/5510 Tutorial 11"
subtitle: "Introduction to Data Analysis"
author: "Patrick Li"
date: "Oct 7, 2024"
format: 
  html:
    toc: true
    embed-resources: true
---


```{r setup, include = FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  eval = FALSE,
  message = FALSE, 
  warning = FALSE,
  error = FALSE, 
  out.width = "70%",
  fig.width = 8, 
  fig.height = 6,
  fig.retina = 3)
set.seed(6)
filter <- dplyr::filter
```


## `r emo::ji("target")` Workshop Objectives

- Introduction to linear models
- Interpreting linear models
- Linear model diagnostic


## `r emo::ji("wrench")` Instructions

1. Log in to RStudio Cloud using Google and enter your Monash user details. 
2. In each question, you will replace '___' with your answer. Please note that the Rmd will not knit until you've answered all of the questions.
3.  Once you have filled up all the blanks, remember to go to `knitr::opts_chunk` at the top of the document, change `eval = TRUE`, then knit the document.
 
##  🌶️Exercise 11A: Introduction to Linear Model

```{r}
library(tidyverse)
library(broom)
```


### Section A: Practice with Advertising Data

#### 1. Read the data into RStudio. Remove the first column.

```{r pp}
advert <- read_csv("data/Advertising.csv") %>% 
  select(___)

head(advert)
```

Before fitting a model, let's plot the data with `sales` on the y-axis, and the rest of the variables on the x-axis in 3 separate plots. 
  * Use the following incomplete set of code as a guide. 
  * Discuss the relationships among `sales` and the other variables.
```{r gg-height-width}
p1 <- ggplot(advert,
       aes(x = ___,
           y = ___)) +
  geom_point(alpha = 0.2)


gridExtra::grid.arrange(p1, p2, p3, nrow = 1)
```



#### 2. Fit a model with just `radio` as the x variable.

The simple linear regression is of the form:

$$y_i=\beta_0+\beta_1x_{i}+\varepsilon_i, ~~~ i=1, \dots, n$$
We fit such models in R using the `lm(response ~ explanatory,data=df)` function. The formula `response ~ explanatory` specifies the response variable and explanatory variable from the data in `df`.  Note that by default, the `lm()` function will include the intercept term.

  
Suppose we are interested in the amount of variation in `sales` that is explained by the variation in `radio`. 
  * Fit the linear model that takes sales as `response` and radio `exlanatory`. 
```{r lm-fit}
advert_m1 <- lm(___, data = advert)
advert_m1
```
  * extract components of the regression line using the `broom` package, and the function `tidy()`.
```{r coef-tidy}
library(___)
coef_advert_radio <- ___(advert_m1)
coef_advert_radio
```

The fitted model has the following functional form:

$$\widehat{\text{sales}} = 9.311 + 0.202~\text{radio}$$

The coefficients can be interpreted as 

**Slope:** For each additional thousand dollars spending on radio advertising, on average, the sales increases by 202 units.

**Intercept:** When there is no advertising spending on radio, sales increases by on average 9 thousand units.

We can draw our fitted linear model with `geom_smooth(method = "lm")`

```{r measure-smooth}
ggplot(advert,
       aes(x = radio,
           y = sales)) +
  geom_point(alpha = 0.2) +
  ___
```

#### 3. Diagnostics and augment

The residuals from the model are the differences between the observed values and the fitted values. 
The residuals and fitted values can be extracted from using `augment()` function as part of the `broom` package:

```{r residuals}
advert_m1_diagnostics <- ___(advert_m1)
advert_m1_diagnostics
```

We can visualise the residuals directly on the scatter plot:

```{r vis-res}
ggplot(advert,
       aes(x = radio,
           y = sales)) +
  geom_point(alpha = 0.2)+
  geom_smooth(method = "lm", se = FALSE) +
  # overlay fitted values
  geom_point(data = advert_m1_diagnostics, 
             aes(y = ___), 
             color = "blue", 
             alpha = 0.2) +
  # draw a line segment from the fitted value to observed value
  geom_segment(data = advert_m1_diagnostics, 
               aes(xend = radio, y = .fitted, yend = sales),
               color = "blue",
               alpha = 0.2)
```

Another way of visualising residuals is to produce a scatterplot of the residuals against the fitted values:

```{r fitted-residuals}
ggplot(advert_m1_diagnostics,
       aes(x = ___, y = ___)) +
  geom_point()
```

Yet, another and a better way to visualise residuals is using `resid_panel()` which gives you a panel of diagnostic residual plots.

```{r}
library(ggResidpanel)
___(advert_m1, plots = "all")
```

Recall the assumptions of the linear model from the lecture. What are they?

- Are they met from the plots?
- Are there observations with large positive residuals? What do they mean?

#### 4. How good is your fit? 

- $R^2$ is a common measurement of strength of linear model fit.
- $R^2$ tells us % variability in response explained by the model.

We can extract model fit summaries using `glance()`

```{r r2}
advert_m1_summary <- ___(advert_m1)
advert_m1_summary
```

The $R^2$ value here can be interpreted as the proportion of variability in sales explained by spending on radio advertisement. That is approximately 33% of variability in sales is explained by spending on radio advertisement. Hmm, that is not much! 😒

### Section B: Multiple Linear Regression

In reality, we will not only spend all the advertising spending on just one medium. We may choose to have other media of advertisement too such as newspaper and TV.

#### 1.  Run a multiple regression model to show the effect of different advertising media on sales. Obtain a summary of the regression results.

```{r}
advert_m2 <- lm(___, 
                data = advert)
___(advert_m2)
```


#### 2. Is the relationship linear?

```{r}
advert_m2_aug <- ___(advert_m2)
ggplot(data = advert_m2_aug,
       aes(x = ___,
           y = ___)) +
  geom_point() +
  geom_hline(yintercept = 0, colour = "red")
```

#### 3.Residual Plots
Residual plots can be used in order to identify non-linearity or other problems with our regression model. 
  * Plot the residuals using `resid_panel()`.
**Take this out of the exercise, but leave in the solutions**
```{r}
library(ggResidpanel)
resid_panel(advert_m2, plots = ___)
```
  * Do you see any patterns? **[ANS: Looks fine but the residuals tends to be more positive when the fitted values are at the lower or higher values.]**





## 👴 Exercise 11B: Data modelling using the gapminder data

### Gapminder

- Hans Rosling was a Swedish doctor, academic and statistician, Professor of International Health at Karolinska Institute. Sadly he passed away in 2017. 
- He developed a keen interest in health and wealth across the globe, and the relationship with other factors like agriculture, education, and energy. 
- You can play with the gapminder data using animations at https://www.gapminder.org/tools/. 

R package `gapminder` contains subset of the data on five year intervals from 1952 to 2007.

```{r show-gapminder}
library(gapminder)
glimpse(gapminder)
```

### Section A: Australia

#### 1. How has life expectancy changed over the years, for each country?
Use a line plot to answer this question and provide your explanation here. 
```{r gg-gapminder-line, fig.height = 4, echo = FALSE}
gapminder %>% 
  ggplot(aes(___, ___,
             group = ___,
             color = continent)) +
    geom_line(alpha = 1/3)

```



#### 2. How has life expectancy in Australia changed over the years. 
Use a line plot to answer this question and provide your explanation here. 

```{r gapminder-oz}
oz <- gapminder %>% 
  ___(country ___)
oz
```

#### 3. Fit the linear model just for Australia: `lifeExp ~ year`

```{r lm-oz-gapminder}
oz_lm <- lm(___, data = oz)

oz_lm

___(oz_lm)
```

#### 5. Tidy the model.

```{r tidy-oz-gapminder}
___(oz_lm)
```

$$\widehat{lifeExp} = -376.1163 + 0.2277~year$$

#### 6. Center year

Let us treat 1950 as the first year, so for model fitting, we are going to shift the years to begin in 1950. This makes interpretation easier.

```{r center-eyar}
gap <- gapminder %>% 
  mutate(year1950 = year - 1950)

oz <- gap %>% 
  filter(country == "Australia")
```

#### 7. Model for centered year

```{r oz-gapminder-centered-year}
oz_lm <- lm(lifeExp ~ ___, data = oz)

oz_lm
```

#### 8. Tidy the model centered year

```{r tidy-oz-gapminder-centered-year}
___(oz_lm)
```

$$\widehat{lifeExp} = 67.9 + 0.2277~year1950 $$


#### 9. Extract residuals and fitted values using `augment()`

```{r oz-gapminder-augment, fig.height=3, fig.width=8}
oz_aug <- ___(oz_lm, oz)
oz_aug

oz_aug$.fitted
oz_aug$.resid
```

#### 10. Plot the `year` and `lifeExp` with points and add in the line.

```{r oz-gap-aug, fig.height=3, fig.width=8}
ggplot(data = oz_aug, 
       aes(x = ___, 
           y = ___)) + 
  ___(colour = "blue") + 
  geom_point(aes(x = year,
                 y = lifeExp))
```

#### 11. Plot residuals against fitted values to reveal problems with the fit and interpret it.

```{r oz-gap-year-resid, fig.height=3, fig.width=8}
ggplot(oz_aug,
       aes(x = ___, y = ___)) +
  ylim(c(-5,5)) +
  geom_point()

# Another way to look at Residuals: .std.resid with year:
ggplot(data = oz_aug, 
             aes(x = year, 
                 y = ___)) +
  ylim(c(-5,5)) + 
  geom_hline(yintercept = 0,
             colour = "white", 
             size = 2)  +
  geom_line() 
```

Provide explanation here.

Note: We use standardized residuals here is to make the residuals unitless. Therefore, it is easier for us to spot any unusual observations, which are those above the line of 2 and below -2. Remember the rules of 95% confidence interval?

### Section B: Other countries?? 

#### 1. Can we fit a model for New Zealand?

```{r lm-nz}
nz <- gap %>% 
  filter(country == "New Zealand")

nz_lm <- lm(___, data = nz)
nz_lm
```

#### 2. Can we fit a model for Japan?

```{r lm-japan}
japan <- gap %>%
  filter(country == ___)

japan_lm <- lm(lifeExp ~ year1950, data = japan)
japan_lm
```

#### 3. Can we fit a model for Italy?

```{r lm-italy}
italy <- gap %>%
  filter(country == ___)

italy_lm <- lm(lifeExp ~ year1950, data = italy)
italy_lm
```

#### 4. Is there a better way? What if we wanted to fit a model for ALL countries?

YES, let's begin. 🙂


1. Nest country level data (one row = one country)

```{r nest}
by_country <- gap %>% 
  ___(country, year1950, lifeExp, continent) %>%
  ___(country, continent) %>% 
  nest()

by_country
```

2. What is in the `data` column?

```{r show-nest}
by_country$data[[1]]
```

It's a tibble! 🤪

Let's fit a linear model to each one of them.

```{r fit-many-model}
lm_afganistan <- lm(lifeExp ~ year1950, 
                    data = by_country$data[[1]])
lm_albania <- lm(lifeExp ~ year1950, 
                 data = by_country$data[[2]])
lm_algeria <- lm(lifeExp ~ year1950, 
                 data = by_country$data[[3]])
```

🙈 But, we are copying and pasting this code more than twice... Is there a better way?

A case for `map`???

`map(<data object>, <function>)`

#### 1. Write a function to fit lm to all countries and map the function to data column of `by_country`

```{r mapped-lm}
fit_lm <- function(x){
  lm(lifeExp ~ year1950, data = x) 
}

mapped_lm <- map(___, fit_lm)

head(mapped_lm)
```

#### 2. Map inside the data?

```{r map-country}
country_model <- by_country %>% 
                    ___(model = map(data, function(x){
                      lm(lifeExp ~ year1950, data = x)
                      })
                      )

country_model
```

#### 3. Case for map (shorthand function)

```{r map-country-short}
country_model <- by_country %>% 
  ___(model = ___(data, ~lm(lifeExp ~ year1950,
                            data = .)))

country_model
```

#### 4. What's the model?

```{r print-model}
country_model$model[[1]]
```

#### 5. How do we summarise the content in the model?

```{r tidy-printed-model}
___(country_model$model[[1]])
```

#### 6. So should we repeat it for each one?

```{r tidy-many-printed-model}
tidy(country_model$model[[1]])
tidy(country_model$model[[2]])
tidy(country_model$model[[3]])
```

NO!! There's a better way. 😙

```{r}
country_model %>%
  mutate(tidy = ___(model, tidy))
```


Data Wrangle and tidy:

```{r map-tidy-model}
country_coefs <- country_model %>%
                    ___(tidy = map(model, tidy)) %>%
                    unnest(tidy) %>%
                    select(country, continent, term, estimate)

country_coefs
```

```{r tidy-unnested-coefs}
tidy_country_coefs <- country_coefs %>%
                          ___(id_cols = c(country, continent), 
                                      names_from =  term,
                                      values_from = estimate) %>%
                          rename(intercept = `(Intercept)`)

tidy_country_coefs
```

#### 7. Check for Australia

```{r oz-unnested-coef, echo=TRUE}
tidy_country_coefs %>%
  filter(country == "Australia")
```

### Section C: Extension Exercise

- Fit the models to all countries.
- Pick your favourite country (other than Australia), print the coefficients, and make a hand sketch of the the model fit.


#### 1. Plot all the models

```{r augmente-many-countries, fig.height=4}
country_aug <- country_model %>% 
  mutate(augmented = map(model, augment)) %>%
  unnest(augmented)

country_aug
```


```{r plot-gapminder-data, fig.height=4}
p1 <- gapminder %>% 
  ggplot(aes(year, lifeExp, group = country)) +
    geom_line(alpha = 1/3) + ggtitle("Data")

p2 <- ggplot(country_aug) + 
  geom_line(aes(x = year1950 + 1950, 
                y = .fitted, 
                group = country), 
            alpha = 1/3) +
  xlab("year") +
  ggtitle("Fitted models")
```

```{r plot-print-gapminder, fig.height = 4}
library(gridExtra)
grid.arrange(p1, p2, ncol=2)
```


#### 2. Plot all the model coefficients

```{r ggplotly-aug-coef}
p <- ggplot(tidy_country_coefs, 
            aes(x = intercept, 
                y = year1950, 
                colour = continent, 
                label = country)) +
  geom_point(alpha = 0.5, 
             size = 2) +
  scale_color_brewer(palette = "Dark2")
```

```{r ggplotly-aug-data}
library(plotly)
ggplotly(p)
```


Let's summarise the information we learnt from the model coefficients.

- Generally, the relationship is negative: this means that if a country starts with a high intercept, it tends to have lower rate of increase.
- There is a difference across the continents: Countries in Europe and Oceania tend to start with a higher life expectancy and increased; countries in Asia and America tended to start lower but have high rates of improvement; Africa tends to start lower and have a huge range in rate of change.
- Three countries had negative growth in life expectancy: Rwand, Zimbabwe, Zambia


#### Model diagnostics by country

```{r glance-country}
country_glance <- country_model %>% 
  mutate(glance = map(model, glance)) %>%
  unnest(glance)

country_glance
```


#### 3. Plot the $R^2$ values using a histogram. 

```{r country-fit, fig.height = 3.5}
ggplot(country_glance, 
       aes(x = r.squared)) + 
  geom_histogram()
```


#### 4. Countries with worst fit

Examine the countries with the worst fit, countries with $R^2<0.45$, by making scatterplots of the data, with the linear model overlaid. 

```{r show-bad-fit}
badfit <- country_glance %>% filter(r.squared <= 0.45)

gap_bad <- gap %>% filter(country %in% badfit$country)

gg_bad_fit <-
ggplot(data = gap_bad, 
       aes(x = year, 
           y = lifeExp)) + 
         geom_point() +
  facet_wrap(~country) +
  scale_x_continuous(breaks = seq(1950,2000,10), 
                     labels = c("1950", "60","70", "80","90","2000")) +
  geom_smooth(method = "lm", 
              se = FALSE)
```


```{r gg-show-bad-fit, fig.height=4}
gg_bad_fit
```

Each of these countries was moving on a nice trajectory of increasing life expectancy, but later suffered a big dip during the time period.