13-model.Rmd

# Linear Model {#model}

### Setup {-}

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

* broom: https://cran.r-project.org/web/packages/broom/index.html
* Introduction to broom: https://cran.r-project.org/web/packages/broom/vignettes/broom.html


## What is modeling in EDA

Model is a simple summary of data

Goal: A simple low-dimensional summary of a dataset. Ideally, the model will capture true “signals” (i.e. patterns generated by the phenomenon of interest), and ignore “noise” (i.e. random variation that you’re not interested in).

![image](data/data-science.png)



```{r}
df0 <- as_tibble(iris) %>% filter(Species != "setosa")
```

Let us look at two charts.


```{r}
df0 %>% ggplot(aes(Petal.Width, Petal.Length)) + geom_point()
```



```{r}
df0 %>% ggplot(aes(Petal.Width, Petal.Length)) + geom_point() + geom_smooth(method="lm",formula=y~x, se=FALSE)
```



```{r}
df0 %>% ggplot(aes(Sepal.Width, Sepal.Length)) + geom_point()
```



```{r}
df0 %>% ggplot(aes(Sepal.Width, Sepal.Length)) + geom_point() + geom_smooth(method="lm",formula=y~x, se=FALSE)
```



## Linear Model: Petal.Length ~ Petal.Width

```{r}
df0 %>% lm(Petal.Length ~ Petal.Width, .)
```



## Formula: $\text{Petal.Length} = 2.224 + 1.600\cdot \text{Petal.Width}$

```{r echo=FALSE}
df0 %>% ggplot(aes(Petal.Width, Petal.Length)) + geom_point() + geom_smooth(method="lm",formula=y~x, se=FALSE)
```



## Linear Model: Sepal.Length ~ Sepal.Width

```{r}
df0 %>% lm(Sepal.Length ~ Sepal.Width, .)
```



## Formula: $\text{Sepal.Length} = 3.093 + 1.103\cdot \text{Sepal.Width}$

```{r echo=FALSE}
df0 %>% ggplot(aes(Sepal.Width, Sepal.Length)) + geom_point() + geom_smooth(method="lm",formula=y~x, se=FALSE)
```



## Petal.Length ~ Petal.Width: R squared = 0.6779 - 68%

```{r}
df0 %>% lm(Petal.Length ~ Petal.Width, .) %>% summary()
```



## Sepal.Length ~ Sepal.Width: R squared = 0.3068 - 31%

```{r}
df0 %>% lm(Sepal.Length ~ Sepal.Width, .) %>% summary()
```



## Linear Model Basics: y ~ x

```
lm(y~x, data)
```
```
data %>% lm(y~x, .)
```

y-intercept, and slope: rate of increase or decrease


```
summary(lm(y~x, data))
```
```
data %>% lm(y~x, .) %>% summary()
```

(Multiple) R Squared: a value between 0 and 1, the model's strength. It is a measurement of the model quality. If the value is close to 1, the model quality is high. If it is close to 0, the model quality is low.



```{r echo=FALSE}
df0 %>% ggplot(aes(Petal.Width, Petal.Length)) + geom_point() + geom_smooth(method="lm",formula=y~x, se=FALSE) + geom_hline(yintercept = mean(df0$Petal.Length), col = "red")
```



```{r eval = FALSE}
df1 <- data.frame(x = c(1,2,3,4), y = c(1,0.5,2, 1.5))
ybar <- mean(df1$y)
mod1 <- lm(y~x, df1)
augment(mod1) %>% ggplot() + geom_point(aes(x,y)) + geom_smooth(aes(x,y), formula = y~x, method = "lm", se = FALSE) + geom_hline(yintercept = ybar, linetype="longdash", col = "red") + geom_point(aes(x, ybar), shape=4) + geom_point(aes(x, .fitted), shape =9, size=2) + geom_text(aes(x, ybar, label = paste0("(",x,",",1.25,")")), nudge_y = -0.1, col = "red") + geom_text(aes(x, .fitted, label = paste0("(",x,",",.fitted,")")), nudge_y = -0.1, col = "blue") + geom_text(aes(x, y, label = paste0("(",x,",",y,")")), nudge_y = -0.1) 
```

* $(x_1, y_1)$, $(x_2,y_2)$, $(x_3, y_3)$, $(x_4, y_4)$: Data points
* $\bar{y}$: mean of y = $(y_1 + y_2 + y_3 + y_4)/4$.
* $\hat{y}_i$: prediction at $x_i$, 
  - $(x_1, \hat{y}_1)$, $(x_2, \hat{y}_2)$, $(x_3, \hat{y}_3)$, $(x_4, \hat{y}_4)$ are on the regression line.
* $y_1-\hat{y}_1$, $y_2-\hat{y}_2$, $y_2-\hat{y}_2$, $y_2-\hat{y}_2$ are called residues.



```{r echo=FALSE}
df1 <- data.frame(x = c(1,2,3,4), y = c(1,0.5,2, 1.5))
ybar <- mean(df1$y)
mod1 <- lm(y~x, df1)
augment(mod1) %>% ggplot() + geom_point(aes(x,y)) + geom_smooth(aes(x,y), formula = y~x, method = "lm", se = FALSE) + geom_hline(yintercept = ybar, linetype="longdash", col = "red") + geom_point(aes(x, ybar), shape=4) + geom_point(aes(x, .fitted), shape =9, size=2) + geom_text(aes(x, ybar, label = paste0("(",x,",",1.25,")")), nudge_y = -0.1, col = "red") + geom_text(aes(x, .fitted, label = paste0("(",x,",",.fitted,")")), nudge_y = -0.1, col = "blue") + geom_text(aes(x, y, label = paste0("(",x,",",y,")")), nudge_y = -0.1) 
```



## R Squared

$$SS_{tot} = (1-1.25)^2 + (0.5-1.25)^2 + (2-1.25)^2 + (1.5-1.25)^2 = 1.25$$
$$SS_{res} = (1-0.8)^2 + (0.5-1.1)^2 + (2-1.4)^2 + (1.5-1.7)^2 = 0.8$$
$$R^2 = 1 - \frac{SS_{res}}{SS_{tot}} = 1- \frac{0.8}{1.25} = 0.36.$$



```{r}
summary(mod1)$r.squared
```

```{r}
mod1 %>% glance() %>% pull(r.squared)
```

```{r}
mod1 %>% glance() %>% select(`R Squared` = r.squared)
```

```{r}
mod1 %>% summary() %>% glimpse()
```




## Useful Mathematical Formula

* Let $x = c(x_1, x_2, \ldots, x_n)$ be the independent variable, i.e., Sepal.L
* Let $y = c(y_1, y_2, \ldots, y_n)$ be the dependent variable, i.e., Sepal.W
* Let $\mbox{pred} = c(\hat{y}_1, \hat{y}_2, \ldots, \hat{y}_n)$ be the predicted values by linear regression.

$$
\begin{aligned}
\mbox{slope of the regression line}  &= \frac{cov(x,y)}{var(x)} = \frac{cor(x,y)\sqrt{var(y)}}{\sqrt{var(x)}}\\
\mbox{total sum of squares} &= SS_{tot} = \sum_{i}(y_i-mean(y))^2\\
\mbox{residual sum of squares} &= SS_{res} = \sum_{i}(y_i-\mbox{pred}_i)^2 = \sum_{i}(y_i-\hat{y}_i)^2\\
\mbox{R squared} = R^2 & = 1 - \frac{SS_{res}}{SS_{tot}} = cor(x,y)^2
\end{aligned}
$$



### Adjusted R Squared

$$\text{Adjusted }R^2 = 1- \frac{(1-R^2)(n-1)}{n-k-1}$$
$n$: number of observations, the number of rows

$k$: number of variables used for prediction



```{r}
df0 %>% select(1:4) %>% cor()
```
```{r}
cormat <- df0 %>% select(1:4) %>% cor()
cormat*cormat
```



```{r}
as_tibble(iris) %>% filter(Species == "setosa") %>% select(-5) %>% cor()
```
```{r}
as_tibble(iris) %>% filter(Species == "virginica") %>% select(-5) %>% cor()
```
```{r}
as_tibble(iris) %>% filter(Species == "versicolor") %>% select(-5) %>% cor()
```



```{r}
as_tibble(iris) %>% filter(Species == "virginica") %>% ggplot(aes(Sepal.Length, Petal.Length)) + geom_point() + geom_smooth(method = "lm", formula = y~x, se = FALSE)
```



```{r}
as_tibble(iris) %>% filter(Species == "virginica") %>% lm(Petal.Length ~ Sepal.Length, .) %>% glance() %>% pull(r.squared) %>% sqrt()
```


Correlations of the data suggest the possible strength of linear model y ~ x.

```{r}
iris %>% select(-5) %>% cor()
```



## Examples: WDI

* SP.DYN.LE00.IN: Life expectancy at birth, total (years)

```{r cash = TRUE, eval=FALSE}
wdi_lifeExp <- WDI(indicator = c(lifeExp = "SP.DYN.LE00.IN"))
```

```{r include=FALSE, eval=FALSE}
write_csv(wdi_lifeExp, "./data/wdi_lifeExp.csv")
```

```{r echo=FALSE}
wdi_lifeExp <- read_csv("./data/wdi_lifeExp.csv")
```



```{r}
wdi_lifeExp %>% filter(country == "World") %>% drop_na(lifeExp) %>%
  ggplot(aes(year, lifeExp)) + geom_point() + geom_smooth(method = "lm", se = FALSE)
```



```{r}
wdi_lifeExp %>% lm(lifeExp ~ year, .) %>% summary()
```

$$lifeExp \sim -557.4 + 0.3123 \cdot year$$
Each year, life expectancy at birth increases approximately 0.3123 years. R-squared of this model is 0.2392, and the model explains 24%.



```{r}
wdi_lifeExp %>% filter(country == "World", year >= 1962, year <= 2019) %>% drop_na(lifeExp) %>% lm(lifeExp ~ year, .) %>% summary()
```



## BRICs

```{r}
mod_brics <- wdi_lifeExp %>% filter(country %in% c("Brazil", "Russian Federation", "India", "China")) %>% drop_na(lifeExp) %>% lm(lifeExp ~ year, .) %>% summary()
mod_brics$r.squared
```


```{r eval=FALSE}
wdi_lifeExp %>% filter(country %in% c("Brazil", "Russian Federation", "India", "China")) %>% drop_na(lifeExp) %>%
  ggplot(aes(year, lifeExp)) + geom_point() + geom_smooth(formula = y~x, method = "lm", se = FALSE)
```



```{r echo=FALSE}
wdi_lifeExp %>% filter(country %in% c("Brazil", "Russian Federation", "India", "China")) %>% drop_na(lifeExp) %>%
  ggplot(aes(year, lifeExp)) + geom_point() + geom_smooth(formula = y~x, method = "lm", se = FALSE)
```





```{r}
wdi_lifeExp %>% filter(country %in% c("Brazil", "Russian Federation", "India", "China")) %>% drop_na(lifeExp) %>%
  ggplot(aes(year, lifeExp, color = country)) + geom_point(aes(shape = country)) + geom_smooth(formula = y~x, method = "lm", se = FALSE)
```


Need to work

```{r eval=FALSE}
country_model <- function(df) {
  lm(lifeExp ~ year, data = df)
}

by_country <- wdi_lifeExp %>% filter(country %in% c("Brazil", "Russian Federation", "India", "China")) %>% drop_na(lifeExp) %>% group_by(country) %>% nest()

by_country <- by_country %>% 
  mutate(model = map(data, country_model))

by_country %>% 
  mutate(tidy = map(model, broom::tidy)) %>% 
  unnest(tidy)

by_country %>% 
  mutate(glance = map(model, broom::glance)) %>% 
  unnest(glance)
```



## Government Expenditure, (% of GDP)

```{r}
wdi_cache <- read_rds("./data/wdi_cache.RData")
```

```{r}
WDIsearch("expenditure", "name", cache = wdi_cache) %>% 
  inner_join(WDIsearch("% of GDP", "name", cache = wdi_cache))
```



```{r}
wdi_cache$series %>% filter(grepl("expenditure", name), grepl("% of GDP", name))
```



* NY.GDP.PCAP.KD: GDP per capita (constant 2015 US$)
* SP.DYN.LE00.IN: Life expectancy at birth, total (years)
* SP.POP.TOTL: Population, total

* GB.XPD.RSDV.GD.ZS: Research and development expenditure (% of GDP) - 2
* MS.MIL.XPND.GD.ZS: Military expenditure (% of GDP) - 6
* SE.XPD.TOTL.GD.ZS: Government expenditure on education, total (% of GDP)


```{r cash = TRUE, eval=FALSE}
wdi_world <- WDI(country = "all", indicator = c(gdpPcap = "NY.GDP.PCAP.KD", lifeExp = "SP.DYN.LE00.IN", pop = "SP.POP.TOTL", research = "GB.XPD.RSDV.GD.ZS", military = "MS.MIL.XPND.GD.ZS", education = "SE.XPD.TOTL.GD.ZS"), 1990, extra = TRUE, cache = wdi_cache)
```

```{r include=FALSE, eval=FALSE}
write_csv(wdi_world, "./data/wdi_world")
```

```{r echo=FALSE}
wdi_world <- read_csv("./data/wdi_world")
```

```{r}
wdi_world
```

SE.XPD.TOTL.GB.ZS: Government expenditure on education, total (% of government expenditure)
SE.XPD.TOTL.GD.ZS: Government expenditure on education, total (% of GDP)
SE.XPD.PRIM.PC.ZS: Government expenditure per student, primary (% of GDP per capita)
MS.MIL.XPND.ZS: Military expenditure (% of general government expenditure)
SE.XPD.TERT.ZS: Expenditure on tertiary education (% of government expenditure on education)


```{r}
mod_e <- lm(lifeExp ~ education, wdi_world); mod_e
```
```{r}
wdi_world %>% ggplot(aes(education, lifeExp)) + geom_point() + geom_smooth(formula = y ~ x, method = "lm", se=FALSE)
```
```{r}
wdi_world %>% filter(income != "Aggregates") %>% drop_na(education, lifeExp) %>% ggplot(aes(education, lifeExp)) + geom_point() + geom_smooth(formula = y ~ x, method = "lm", se=FALSE)
```

```{r}
wdi_world_el <- wdi_world %>% select(country, year, education, lifeExp, gdpPcap, pop, research, military, region, income) %>% filter(income != "Aggregates") %>% drop_na(education, lifeExp)
```

```{r}
wdi_world_el %>% ggplot(aes(education)) + geom_histogram()
```

```{r}
wdi_world_el %>% filter(year==2020) %>% ggplot(aes(x = income, y = education, fill = income)) + geom_boxplot()
```

```{r}
wdi_world_el %>% filter(year==2020) %>% arrange(desc(education))
```


```{r}
wdi_world_el %>% filter(year==2020) %>% arrange(desc(education))
```

```{r}
wdi_world_el %>% filter(year==2020) %>% lm(gdpPcap ~ education, .)
```
```{r}
wdi_world_el %>% filter(year==2020) %>% lm(gdpPcap ~ education, .) %>% glance()
```
```{r}
wdi_world_el %>% lm(lifeExp ~ education + research + military, .) %>% glance()
```

```{r}
wdi_world_el %>% lm(lifeExp ~ education + research + military, .) %>% tidy()
```

$$lifeExp \sim 70.22 + 0.08\cdot education + 3.84 \cdot research - 0.07 \cdot military$$

```{r}
wdi_world_el %>% lm(gdpPcap ~ education + research + military, .) %>% tidy()
```

```{r}
wdi_world_el %>% lm(gdpPcap ~ education + research + military, .) %>% glance()
```

$$gdpPcap \sim 1077 + 1024\cdot education + 12792 \cdot research - 967 \cdot military$$

```{r}
mod_r <- lm(lifeExp ~ research, wdi_world); mod_e
```

## model and Linear Regression Quick Reference

* R4DS: Model basics
  - https://r4ds.had.co.nz/model-basics.html

For explanation of other indices, please see.

* r-statistics.co by Selva Prabhakaran: 
  - http://r-statistics.co/Linear-Regression.html


* Meaning Behind Each Section of Summary()
  - https://www.learnbymarketing.com/tutorials/explaining-the-lm-summary-in-r/