---
title: "SML201 Week 4 Lecture 1: Predictive Modelling with Linear Regression"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
```

### Linear regression

Let's run our first linear regression on the `gapminder` dataset for 1982.

```{r}
library(gapminder)
gap.1982 <- gapminder %>% filter(year == 1982)
p <- ggplot(data = gap.1982,
            mapping = aes(x = gdpPercap,
                          y = lifeExp))
p + geom_point() +
    geom_smooth(method = "lm")
```


```{r}

fit.gap.1 <- lm(lifeExp ~ gdpPercap, data=gap.1982)
```

We are trying to predict `lifeExp` (so the y here) using `gdpPercap` (the x). We can see the coefficients (the $a_0$ and $a_1$ using)

We can extract the individual coefficients using 

```{r}
fit.gap.1$coefficients[1]
fit.gap.1$coefficients[2]
```

So we can get predictions using

```{r}
my.pred <- fit.gap.1$coefficients[1] + fit.gap.1$coefficients[2] * gapminder$gdpPercap
my.pred[1:10]
```

We can also get predictions using

```{r}
predict(fit.gap.1, newdata=gap.1982)[1:10]
```

(`newdata` needs to be a data frame with the column `gdpPercap`)

Let's now compute the cost function $\sum_i (y^{(i)}-pred(x^{(i)}))^2$:

```{r}
sum((predict(fit.gap.1, gap.1982) - gap.1982$lifeExp)**2)
```

This looks large, but partly that's because we are summing up all the squared errors. Let's compute the *mean* squared error

```{r}
sum((predict(fit.gap.1, gap.1982) - gap.1982$lifeExp)**2)/nrow(gap.1982)
```

(Note: `nrow` computes the number of rows in a data frame).

This corresponds to the mean square error, $\frac{\sum_i (y^{(i)}-pred(x^{(i)}))^2}{m}$.

Now if we want to get some idea of what out actual errors might be, we'd take the square root: $\sqrt{\frac{\sum_i (y^{(i)}-pred(x^{(i)}))^2}{m}}$


```{r}
sqrt(sum((predict(fit.gap.1, gap.1982) - gap.1982$lifeExp)**2)/nrow(gap.1982))
```

So in some sense, we're can expect to be off by about 8 years. Pretty bad!

Recall that before we had the idea of plotting the graph with the x-axis scaled logarithmically.

```{r}
p <- ggplot(data = gap.1982,
            mapping = aes(x = gdpPercap,
                          y = lifeExp))
p + geom_point() +
    scale_x_log10() +  
    geom_smooth(method = "lm")  
```

This seems like the line fits much better. Will this be reflected in the root mean squared prediction error?

```{r}
fit.xlog.1982 <- lm(lifeExp ~ log(gdpPercap), data = gap.1982)
sqrt(sum((predict(fit.xlog.1982, newdata = gap.1982) - gap.1982$lifeExp)**2)/nrow(gap.1982))
```

Yes!

Now, let's try to use the continent as a predictor variable.

```{r}
fit.xlog.cont.1982 <- lm(lifeExp ~ log(gdpPercap) + continent, data = gap.1982)
sqrt(sum((predict(fit.xlog.cont.1982, newdata = gap.1982) - gap.1982$lifeExp)**2)/nrow(gap.1982))
```

We are now below a RMSE of 5 years.

Let's see what the coefficients for the categorical variables are
```{r}
fit.xlog.cont.1982
```

So, for example, if the continent is Asia, we'd add `6.473` to the predicted life expectancy.