01-regression-simple.Rmd

<!--
kate: indent-width 4; word-wrap-column 74; default-dictionary en_AU
Copyright (C) 2020-2021, Marek Gagolewski <https://www.gagolewski.com>
This material is licensed under the Creative Commons BY-NC-ND 4.0 License.
-->

# Simple Linear Regression


{ LATEX \renewcommand{\mathbf}[1]{\mathbfup{#1}} }
{ LATEX \renewcommand{\boldsymbol}[1]{\mathbfit{#1}} }


```{r chapter-header-motd,echo=FALSE,results="asis"}
cat(readLines("chapter-header-motd.md"), sep="\n")
```



<!-- TODO


"Good" data sets are hard to find, most textbooks fool you



move some quality metrics from Ch.2 to Ch.1

what is high correlation? 0.6 is not! time to stop tolerating weak
associations in the social sciences, you can do better, get
more data or don't do what you do... :/

Social sciences is built on  sand

If your data is crap or your models are weak, admit it and give up.




different coefficients of correlation (Spearman, Kendall) → Chapter 2



We need more than just ML:

Mathematical modelling, differential equations, discrete maths

Statistics for small samples but not needed for big data




It's not all in the data, real life processes are not static

Machine learning models are inherently static, they describe today. Tomorrow may be different, especially if people start to trick the system. Think spam detection , spammers will change their behaviour as new detection algorithms arise

People are smarter than your think, they will challenge the system, by choosing a quality metric you change the direction, but can you really predict what other factors gonna be like when we get there?


Examples : covid prediction not taking into account country policies and other things going on

Garmin vo2max estimation and the risk of going mad/overtraining


-->


<!-- TODO


# What is AI?

**Artificial Intelligence** (AI) aims to construct algorithms and devices
to *make decisions that maximise their chance of successfully achieving
their goals*.



We build AI systems to **solve various problems** under
the assumption of **availability of data**.



Crucial features of AI systems include the ability to:

* perceive the environment

    * input data can come from various sources, e.g., physical sensors,
    files, databases, (pseudo)random number generators
    * data can be of different forms, e.g., vectors, matrices
    and other tensors, graphs, audio/video streams, text

* learn from past experience

    * information on how to correctly behave (what is expected of the system)
    in certain situations might be available
    * sometimes feedback can be given after making a decision ("that was nice!
    well done! look at you go!")
    * some systems can adapt based on current experience as well
    (learn to modify its behaviour *on-line*)

* make decisions, predictions, explanations

    * outputs can be interpreted by a human being and used (critically) for whatever  purpose
    (e.g., a GP for diagnosing illnesses)
    * outputs can be provided to other algorithms to take immediate actions
    (e.g., buy those market shares immediately!)



# Common Applications


Application domains of AI include, but are not limited to:


* Life science and medical data analysis
* Financial services (banking, insurance, investment funds)
* Oil and gas and other energy (solar, mining)
* Real estate
* Pharmaceuticals
* Scientific simulations
* Genomics
* Advertising
* Transportation
* Retail
* Healthcare
* Food production

> Discussion (exercise):
Think of different ways in which
these sectors could benefit from AI solutions.

Basically, wherever we have data and there is a need to improve
things, there is a place for AI solutions.

In the course of this book,
we're going to take a deep dive into what AI really is, what
are its limitations but also how we can use it for improving different processes.

 *AI is a tool*.
Yes, AI is really useful, but it's not a remedy for all our problems.


Behind every success story related to  AI, there are always people involved. Sometimes large
teams of:

* those who deal with hardware, data storage, high performance computing
* those who wrangle data, fit models to data and develop prototypes
* those who are responsible for making the algorithms
work well in production environments
* those who are responsible for providing you
with software/algorithms/models that enable everyone to perform the above,
e.g., languages such as R and Python, libraries, packages

You can call them data engineers, data analysts, machine learning engineers,
data scientists, whatever.

We would like to appreciate the role of open source software.
Most started as non-for-profit and it's beautiful that some of them
continue this way.
Not everything is for sale. Not everything has its price.
Let's not forget it.

We can do a lot of great work for greater good,
with the increased availability of open data,
everyone can be a reporter, an engaged citizen that seeks for truth.
There are NGOs.
Finally, there are researchers (remember that in many countries
still the university's main role is to spread knowledge and not make money!)
that need these methods to make new discoveries, e.g., in psychology,
sociology, agriculture, engineering, biotechnology, pharmacy, medicine, you name it.
-->

## Machine Learning

### What is Machine Learning?




An **algorithm** is a well-defined sequence of instructions that,
for a given sequence of input arguments,
yields some desired output.

In other words, it is a specific recipe for a **function**.

Developing algorithms is a tedious task.





In **machine learning**, we build and study computer algorithms
that make *predictions* or *decisions* but which are not
manually programmed.

**Learning** needs some material based upon which new knowledge is to be acquired.

In other words, we need **data**.



### Main Types of Machine Learning Problems




Machine Learning Problems include, but are not limited to:

- **Supervised learning** -- for every input point (e.g., a photo)
there is an associated desired output (e.g., whether it depicts a crosswalk
or how many cars can be seen on it)

- **Unsupervised learning** -- inputs are unlabelled, the aim is to discover
the underlying structure in the data (e.g., automatically group customers
w.r.t. common behavioural patterns)

- **Semi-supervised learning** -- some inputs are labelled, the others
are not (definitely a cheaper scenario)


- **Reinforcement learning** -- learn to act based on a
feedback given after the actual decision was made
(e.g., learn to play The Witcher 7 by testing different hypotheses
what to do to survive as long as possible)




## Supervised Learning

### Formalism




Let $\mathbf{X}=\{\mathfrak{X}_1,\dots,\mathfrak{X}_n\}$
be an input sample ("a database")
that consists of $n$ objects.



Most often we assume that each object $\mathfrak{X}_i$
is represented using $p$ numbers for some $p$.

We denote this fact as $\mathfrak{X}_i\in \mathbb{R}^p$
(it is *a $p$-dimensional real vector* or
*a sequence of $p$ numbers* or
*a point in a $p$-dimensional real space*
or *an element of a real $p$-space* etc.).

If we have "complex" objects on input,
we can always try representing them as **feature vectors** (e.g.,
come up with numeric attributes that best describe them in a task at hand).

{ BEGIN exercise }
Consider the following problems:

1. How would you represent a patient in a clinic?

2. How would you represent a car in an insurance company's database?

3. How would you represent a student in an university?
{ END exercise }





Of course, our setting is *abstract* in the sense that
there might be different realities *hidden* behind these symbols.

This is what maths is for -- creating *abstractions* or *models*
of complex entities/phenomena so that they can be much more easily manipulated
or understood.
This is very powerful -- spend a moment
contemplating how many real-world situations fit into our framework.



This also includes image/video data, e.g., a 1920×1080 pixel image
can be "unwound" to a "flat" vector of length 2,073,600.

(\*) There are some algorithms such as Multidimensional Scaling,
Locally Linear Embedding, IsoMap etc.
that can do that automagically.

. . .



In cases such as this we say that we deal with *structured (tabular) data*\
-- $\mathbf{X}$ can be written as an ($n\times p$)-matrix:
\[
\mathbf{X}=
\left[
\begin{array}{cccc}
x_{1,1} & x_{1,2} & \cdots & x_{1,p} \\
x_{2,1} & x_{2,2} & \cdots & x_{2,p} \\
\vdots & \vdots & \ddots & \vdots \\
x_{n,1} & x_{n,2} & \cdots & x_{n,p} \\
\end{array}
\right]
\]
Mathematically, we denote this as $\mathbf{X}\in\mathbb{R}^{n\times p}$.





Remark.

: Structured data == think: Excel/Calc spreadsheets, SQL tables etc.








For an example, consider the famous Fisher's Iris flower dataset,
see `?iris` in R
and https://en.wikipedia.org/wiki/Iris_flower_data_set.

```{r show_iris}
X <- iris[1:6, 1:4] # first 6 rows and 4 columns
X         # or: print(X)
dim(X)    # gives n and p
dim(iris) # for the full dataset
```




$x_{i,j}\in\mathbb{R}$
represents the $j$-th feature of the $i$-th observation,
$j=1,\dots,p$, $i=1,\dots,n$.

For instance:

```{r show_iris_elem}
X[3, 2] # 3rd row, 2nd column
```





The third observation (data point, row in $\mathbf{X}$)
consists of items $(x_{3,1}, \dots, x_{3,p})$ that can be extracted by calling:

```{r show_iris_row}
X[3,]
as.numeric(X[3,]) # drops names
```

```{r show_iris_row_length}
length(X[3,])
```




Moreover, the second feature/variable/column
is comprised of
$(x_{1,2}, x_{2,2}, \dots, x_{n,2})$:

```{r show_iris_col}
X[,2]
length(X[,2])
```





We will sometimes use the following notation to emphasise that
the $\mathbf{X}$ matrix  consists of $n$ rows
or $p$ columns:

\[
\mathbf{X}=\left[
\begin{array}{c}
\mathbf{x}_{1,\cdot} \\
\mathbf{x}_{2,\cdot} \\
\vdots\\
\mathbf{x}_{n,\cdot} \\
\end{array}
\right]
=
\left[
\begin{array}{cccc}
\mathbf{x}_{\cdot,1} &
\mathbf{x}_{\cdot,2} &
\cdots &
\mathbf{x}_{\cdot,p} \\
\end{array}
\right].
\]




Here, $\mathbf{x}_{i,\cdot}$ is a *row vector* of length $p$,
i.e., a $(1\times p)$-matrix:

\[
\mathbf{x}_{i,\cdot} = \left[
\begin{array}{cccc}
x_{i,1} &
x_{i,2} &
\cdots &
x_{i,p} \\
\end{array}
\right].
\]

Moreover, $\mathbf{x}_{\cdot,j}$ is a *column vector* of length $n$,
i.e., an $(n\times 1)$-matrix:

\[
\mathbf{x}_{\cdot,j} = \left[
\begin{array}{cccc}
x_{1,j} &
x_{2,j} &
\cdots &
x_{n,j} \\
\end{array}
\right]^T=\left[
\begin{array}{c}
{x}_{1,j} \\
{x}_{2,j} \\
\vdots\\
{x}_{n,j} \\
\end{array}
\right],
\]

where $\cdot^T$ denotes the *transpose* of a given matrix --
think of this as a kind of rotation; it allows us to introduce a set of
"vertically stacked" objects using a single inline formula.



### Desired Outputs




In supervised learning,
apart from the inputs we are also given the corresponding
reference/desired outputs.

The aim of supervised learning is to try to create an "algorithm" that,
given an input point, generates an output that is as *close* as possible
to the desired one. The given data sample will be used to "train" this "model".

Usually the reference outputs are encoded as individual  numbers (scalars)
or textual labels.

In other words, with each input $\mathbf{x}_{i,\cdot}$ we associate
the desired output $y_i$:

```{r show_iris_sample,echo=-1}
set.seed(123)
# in iris, iris[, 5] gives the Ys
iris[sample(nrow(iris), 3), ]  # three random rows
```





Hence, our dataset is $[\mathbf{X}\ \mathbf{y}]$ --
each object is represented as a row vector
$[\mathbf{x}_{i,\cdot}\ y_i]$, $i=1,\dots,n$:

\[
[\mathbf{X}\ \mathbf{y}]=
\left[
\begin{array}{ccccc}
x_{1,1} & x_{1,2} & \cdots & x_{1,p} & y_1\\
x_{2,1} & x_{2,2} & \cdots & x_{2,p} & y_2\\
\vdots & \vdots & \ddots & \vdots    & \vdots\\
x_{n,1} & x_{n,2} & \cdots & x_{n,p} & y_n\\
\end{array}
\right],
\]

where:

\[
\mathbf{y} = \left[
\begin{array}{cccc}
y_{1} &
y_{2} &
\cdots &
y_{n} \\
\end{array}
\right]^T=\left[
\begin{array}{c}
{y}_{1} \\
{y}_{2} \\
\vdots\\
{y}_{n} \\
\end{array}
\right].
\]



### Types of Supervised Learning Problems




Depending on the type of the elements in $\mathbf{y}$
(the domain of $\mathbf{y}$),
supervised learning problems are usually
classified as:

- **regression** -- each $y_i$ is a real number

    e.g., $y_i=$ future market stock price
    with $\mathbf{x}_{i,\cdot}=$ prices from $p$ previous days

- **classification** -- each $y_i$ is a discrete label

    e.g., $y_i=$ healthy (0) or ill (1)
    with $\mathbf{x}_{i,\cdot}=$ a patient's health record


- **ordinal regression** (a.k.a. ordinal classification) -- each $y_i$ is a rank

    e.g., $y_i=$ rating of a product on the scale 1--5
    with $\mathbf{x}_{i,\cdot}=$ ratings of $p$ most similar products





{ BEGIN exercise }
Example Problems -- Discussion:

Which of the following are instances of classification problems? Which of them are regression tasks?

What kind of data should you gather in order to tackle them?

- Detect email spam
- Predict a market stock price
- Predict the likeability of a new ad
- Assess credit risk
- Detect tumour tissues in medical images
- Predict time-to-recovery of cancer patients
- Recognise smiling faces on photographs
- Detect unattended luggage in airport security camera footage
- Turn on emergency braking to avoid a collision with pedestrians
{ END exercise }




A single dataset can become an instance of many different ML problems.

Examples -- the `wines` dataset:

```{r wines_load,echo=-1}
set.seed(123)
wines <- read.csv("datasets/winequality-all.csv", comment="#")
wines[1,]
```


`alcohol` is a numeric (quantitative) variable (see Figure \@ref(fig:wines-regression) for a histogram depicting its empirical distribution):


```{r wines-regression,fig.cap="Quantitative (numeric) outputs lead to regression problems"}
summary(wines$alcohol) # continuous variable
hist(wines$alcohol, main="", col="white"); box()
```

`color` is a quantitative variable with two possible outcomes (see Figure \@ref(fig:wines-binary) for a bar plot):



```{r wines-binary,fig.cap="Quantitative outputs lead to classification tasks"}
table(wines$color) # binary variable
barplot(table(wines$color), col="white", ylim=c(0, 6000))
```


Moreover, `response` is an ordinal variable, representing
a wine's rating as assigned by a wine expert
(see  Figure \@ref(fig:wines-ordinal) for a barplot).
Note that although the ranks are represented with numbers,
we they are not continuous variables. Moreover,
these ranks are something more than just labels -- they are linearly
ordered, we know what's the smallest rank and whats the greatest one.

```{r wines-ordinal,fig.cap="Ordinal variables constitute ordinal regression tasks"}
table(wines$response) # ordinal variable
barplot(table(wines$response), col="white", ylim=c(0, 3000))
```




## Simple Regression

### Introduction




**Simple regression** is the easiest setting to start with -- let's assume
$p=1$, i.e., all inputs are 1-dimensional.
Denote $x_i=x_{i,1}$.

We will use it to build many intuitions, for example, it'll be easy
to illustrate all the concepts graphically.

```{r credit-scatter,fig.cap="A scatter plot of  Rating vs. Balance"}
library("ISLR") # Credit dataset
plot(Credit$Balance, Credit$Rating) # scatter plot
```


In what follows we will be modelling the Credit Rating ($Y$)
as a function of the average Credit Card Balance ($X$) in USD
for customers *with positive Balance only*.
It is because it is evident from Figure \@ref(fig:credit-scatter)
that some customers with zero balance obtained a credit rating
based on some external data source that we don't have access to in
our very setting.

```{r credit-XY}
X <- as.matrix(as.numeric(Credit$Balance[Credit$Balance>0]))
Y <- as.matrix(as.numeric(Credit$Rating[Credit$Balance>0]))
```

Figure \@ref(fig:credit-XY-plot) gives the updated scatter plot
with the zero-balance clients "taken care of".

```{r credit-XY-plot,fig.cap="A scatter plot of  Rating vs. Balance with clients of Balance=0 removed"}
plot(X, Y, xlab="X (Balance)", ylab="Y (Rating)")
```

Our aim is to construct a function $f$ that
**models** Rating as a function of Balance,
$f(X)=Y$.

We are equipped with $n=`r length(X)`$ reference (observed) Ratings
$\mathbf{y}=[y_1\ \cdots\ y_n]^T$
for particular Balances $\mathbf{x}=[x_1\ \cdots\ x_n]^T$.





Note the following naming conventions:

* Variable types:

    - $X$ -- independent/explanatory/predictor variable

    - $Y$ -- dependent/response/predicted variable

* Also note that:

    - $Y$ -- idealisation (any possible Rating)

    - $\mathbf{y}=[y_1\ \cdots\ y_n]^T$ -- values actually observed







The model will not be ideal, but it might be usable:

- We will be able to **predict** the rating of any new client.


    What should be the rating of a client with Balance of \$1500?

    What should be the rating of a client with Balance of \$2500?

- We will be able to **describe** (understand) this reality using a single mathematical formula
so as to infer that there is an association between $X$ and $Y$

    Think of "data compression" and laws of physics, e.g., $E=mc^2$.








(\*) Mathematically, we will assume that there is some "true" function that models the data
(true relationship between $Y$ and $X$),
but the observed outputs are subject to **additive error**:
\[Y=f(X)+\varepsilon.\]


$\varepsilon$ is a random term, classically we assume that
errors are independent of each other,
have expected value of $0$ (there is no systematic error = unbiased)
and that they follow a normal distribution.




(\*) We denote this as $\varepsilon\sim\mathcal{N}(0, \sigma)$
(read: random variable $\varepsilon$ follows a normal distribution
with expected value of $0$ and standard deviation of $\sigma$ for some $\sigma\ge 0$).

$\sigma$ controls the amount of noise (and hence, uncertainty).
Figure \@ref(fig:normal-distribs) gives the plot of the probability
distribution function (PDFs, densities)
of $\mathcal{N}(0, \sigma)$ for different $\sigma$s:



```{r normal-distribs, echo=FALSE, message=FALSE, fig.cap="Probability density functions of normal distributions with different standard deviations $\\sigma$."}
y <- seq(-3,3,length.out=101)
plot(y, dnorm(y, mean=0), type='l', ann=FALSE, ylim=c(0,dnorm(0,0,0.5)))
lines(y, dnorm(y, mean=0, sd=0.5), type='l', col=2, lty=2)
legend("topleft", lty=c(1,2), col=c(1,2), legend=expression(sigma*"="*1, sigma*"="*0.5), bg="white")
abline(v=0, lty=3, col="gray")
```



### Search Space and Objective




There are many different functions that can be **fitted** into
the observed $(\mathbf{x},\mathbf{y})$,
compare Figure \@ref(fig:credit-different-models).
Some of them are better than the other (with respect to
different aspects, such as fit quality, simplicity etc.).


```{r credit-different-models, echo=FALSE, message=FALSE,fig.cap="Different polynomial models fitted to data"}
library("RColorBrewer")
library("Cairo")
pal <- RColorBrewer::brewer.pal(6, "Dark2")
plot(X, Y, col="#000000ff")

#library("mda")
#f <- mars(as.matrix(X), Y, degree=5)
#x <- seq(min(X), max(X), length.out=101)
#y <- predict(f, x)
#lines(x,y,col=4,lty=4,lwd=2)

f <- lm(Y~poly(X, 2))
x <- seq(min(X), max(X), length.out=101)
y <- predict(f, data.frame(X=x))
lines(x,y,col=pal[1],lty=1,lwd=3)

f <- lm(Y~poly(X, 15))
x <- seq(min(X), max(X), length.out=101)
y <- predict(f, data.frame(X=x))
lines(x,y,col=pal[2],lty=2,lwd=3)

lines(X[order(X)], Y[order(Y)],col=pal[3],lty=3,lwd=3)

abline(lm(Y~X), col=pal[4],lwd=3,lty=4)
```

Thus, we need a formal **model selection criterion**
that could enable as to tackle the model fitting
task on a computer.







Usually, we will be interested in a model
that minimises appropriately aggregated **residuals**
$f(x_i)-y_i$, i.e.,
**predicted outputs minus observed outputs**,
often denoted with $\hat{y}_i-y_i$,
for $i=1,\dots,n$.

```{r credit-residuals, echo=FALSE, message=FALSE,fig.cap="Residuals are defined as the differences between the predicted and observed outputs $\\hat{y}_i-y_i$"}
plot(X, Y, col=1, xlim=c(1500,2000), ylim=c(400, 1000))
f <- lm(Y~poly(X, 2))
x <- seq(min(X), max(X), length.out=101)
y <- predict(f, data.frame(X=x))
pal <- RColorBrewer::brewer.pal(6, "Pastel1")
lines(x,y,col=pal[1],lty=1,lwd=2)

Y2 <- predict(f, data.frame(X=X))

segments(X, Y, X, Y2, lty=2)
points(X, Y2, pch=7, col="red")


legend("bottomright",
    legend=c("observed output", "predicted output", "fitted model"),
    pch=c(1,7, NA), lty=c(NA, NA, 1), lwd=c(NA, NA, 2),
    col=c(1, "red", pal[1]), bg="white")
```


In Figure \@ref(fig:credit-residuals), the residuals correspond to the
lengths of the dashed line segments -- they measure the discrepancy between
the outputs generated by the model (what we get)
and the true outputs (what we want).


> Note that some sources define residuals as $y_i-\hat{y}_i=y_i-f(x_i)$.



Top choice: sum of squared residuals:
\[
\begin{array}{rl}
\mathrm{SSR}(f|\mathbf{x},\mathbf{y})
& = \left( f(x_1)-y_1 \right)^2 + \dots + \left( f(x_n)-y_n \right)^2 \\
& =
\displaystyle\sum_{i=1}^n \left( f(x_i)-y_i \right)^2
\end{array}
\]


Remark.

: Read "$\sum_{i=1}^n z_i$" as "the sum of $z_i$ for $i$ from $1$ to $n$";
this is just a shorthand for $z_1+z_2+\dots+z_n$.



Remark.

: The notation $\mathrm{SSR}(f|\mathbf{x},\mathbf{y})$ means that
it is the error measure
corresponding to the model $(f)$ *given* our data.\
We could've  denoted it  with $\mathrm{SSR}_{\mathbf{x},\mathbf{y}}(f)$
or even $\mathrm{SSR}(f)$ to emphasise that $\mathbf{x},\mathbf{y}$
are just fixed  values and we are not  interested
in changing them at all (they are "global variables").






We enjoy SSR because (amongst others):

- larger errors are penalised much more than smaller ones

    > (this can be considered a drawback as well)



- (\*\*) statistically speaking, this has a clear underlying interpretation

    > (assuming errors are normally distributed,
    finding a model minimising the SSR is equivalent
    to maximum likelihood estimation)

- the models minimising the SSR can often be found easily

    > (corresponding optimisation tasks have an analytic solution --
    studied already by Gauss in the late 18th century)


(\*\*) Other choices:

- regularised SSR, e.g., lasso or ridge regression (in the case of multiple input variables)
- sum or median of absolute values (robust regression)








Fitting a model to data can be written as an optimisation problem:

\[
\min_{f\in\mathcal{F}} \mathrm{SSR}(f|\mathbf{x},\mathbf{y}),
\]

i.e., find $f$ minimising the SSR **(seek "best" $f$)**\
amongst the set of admissible models $\mathcal{F}$.



Example $\mathcal{F}$s:

- $\mathcal{F}=\{\text{All possible functions of one variable}\}$  --  if there are no repeated
$x_i$'s, this corresponds to data *interpolation*; note that there
are many functions that give SSR of $0$.

- $\mathcal{F}=\{ x\mapsto x^2, x\mapsto \cos(x), x\mapsto \exp(2x+7)-9 \}$ -- obviously
an ad-hoc choice but you can easily choose the best amongst the 3 by computing 3 sums of squares.

- $\mathcal{F}=\{ x\mapsto a+bx\}$ -- the space of linear functions of one variable

- etc.


(e.g., $x\mapsto x^2$ is read "$x$ maps to $x^2$" and is
an elegant way to define an inline function $f$ such that $f(x)=x^2$)




## Simple Linear Regression

### Introduction




If the family of admissible models $\mathcal{F}$ consists only of all linear functions of one variable,
we deal with a **simple linear regression**.

Our problem becomes:


\[
\min_{a,b\in\mathbb{R}} \sum_{i=1}^n \left(
a+bx_i-y_i
\right)^2
\]


In other words, we seek best fitting line in terms of the squared residuals.

This is the **method of least squares**.

This is particularly nice, because our search space
is just $\mathbb{R}^2$ -- easy to handle both analytically and numerically.




```{r credit-different-lines-ssr,echo=FALSE,fig.height=4,fig.cap="Three simple linear models together with the corresponding SSRs"}
plot(X, Y, col="#00000044")
m1 <- lm(Y~X)
abline(m1, col=2,lwd=3)
m2 <- lm(Y~X+0)
abline(m2, col=3,lwd=3)
m3 <- L1pack::lad(Y~X)
abline(m3, col=4,lwd=3)
legend("topleft",
    lty=1,
    col=c(2,3,4),
    legend=c(
        sprintf("y=%03.0f+%.3fx; SSR=%.0f",m1$coefficients[1],m1$coefficients[2], sum(m1$residuals^2)),
        sprintf("y=%03.0f+%.3fx; SSR=%.0f",0,                 m2$coefficients[1], sum(m2$residuals^2)),
        sprintf("y=%03.0f+%.3fx; SSR=%.0f",m3$coefficients[1],m3$coefficients[2], sum(m3$residuals^2))
    ), bg="white")
```


{ BEGIN exercise }
Which of the lines in Figure \@ref(fig:credit-different-lines-ssr) is the least squares solution?
{ END exercise }


<!-- TODO:
heatmap of SSR(a,b) for a = [0,1], b=[0,300]
-->





### Solution in R




Let's fit the linear model minimising the SSR in R.
The `lm()` function (`l`inear `m`odels) has a convenient *formula*-based interface.


```{r credit-fi-lm}
f <- lm(Y~X)
```

In R, the expression "`Y~X`" denotes a formula, which we read as:
variable `Y` is a function of `X`.
Note that the dependent variable is on the left side of the formula.
Here, `X` and `Y` are two R numeric vectors of identical lengths.


Let's print the fitted model:

```{r credit-print-lm}
print(f)
```

Hence, the fitted model is:
\[
Y = f(X) = `r f$coefficient[1]`+`r f$coefficient[2]` X
\qquad (+ \varepsilon)
\]


Coefficient $a$ (intercept):

```{r credit_print_lm_2}
f$coefficient[1]
```

Coefficient $b$ (slope):

```{r credit_print_lm_1}
f$coefficient[2]
```

Plotting, see Figure \@ref(fig:credit-plot-lm):

```{r credit-plot-lm, fig.cap="Fitted regression line"}
plot(X, Y, col="#000000aa")
abline(f, col=2, lwd=3)
```


SSR:

```{r credit_ssr_lm}
sum(f$residuals^2)
sum((f$coefficient[1]+f$coefficient[2]*X-Y)^2) # equivalent
```

We can predict the model's output for yet-unobserved inputs by
writing:

```{r credit_predict_lm}
X_new <- c(1500, 2000, 2500) # example inputs
f$coefficient[1] + f$coefficient[2]*X_new
```


Note that linear models can also be fitted based on formulas that refer
to a data frame's columns. For example, let us wrap
both $\mathbf{x}$ and $\mathbf{y}$ inside a data frame:

```{r credit_XY_df}
XY <- data.frame(Balance=X, Rating=Y)
head(XY, 3)
```

By writing:

```{r credit_XY_df_lm}
f <- lm(Rating~Balance, data=XY)
```

now `Balance` and `Rating` refer to column names in the `XY` data frame,
and not the objects in R's "workspace".

Based on the above, we can make a prediction
using the `predict()` function"

```{r credit_XY_predict}
X_new <- data.frame(Balance=c(1500, 2000, 2500))
predict(f, X_new)
```

Interestingly:

```{r credit_XY_predict_extrapolate}
predict(f, data.frame(Balance=c(5000)))
```

This is more than the highest possible rating -- we have extrapolated
way beyond the observable data range.




Note that our $Y=a+bX$ model is **interpretable**
and **well-behaving**
(not all machine learning models will have this feature,
think: deep neural networks,
which we rather conceive as *black boxes*):

* we know that by increasing $X$ by a small amount,
$Y$ will also increase (positive correlation),

* the model is continuous -- small change in $X$
doesn't yield any drastic change in $Y$,

* we know what will happen if we increase or decrease $X$ by, say, $100$,

* the function is invertible -- if we want Rating of $500$,
we can compute the associated preferred Balance that should yield it
(provided that the model is valid).


### Analytic Solution


It may be shown (which we actually do below)
that the solution is:

\[
\left\{
\begin{array}{rl}
b  = & \dfrac{
n \displaystyle\sum_{i=1}^n x_i y_i - \displaystyle\sum_{i=1}^n  y_i \displaystyle\sum_{i=1}^n x_i
}{
n \displaystyle\sum_{i=1}^n x_i x_i -   \displaystyle\sum_{i=1}^n x_i\displaystyle\sum_{i=1}^n x_i
}\\
a = & \dfrac{1}{n}\displaystyle\sum_{i=1}^n  y_i - b  \dfrac{1}{n} \displaystyle\sum_{i=1}^n x_i  \\
\end{array}
\right.
\]

Which can be implemented in R as follows:

```{r lm_estimate_manual}
n <- length(X)
b <- (n*sum(X*Y)-sum(X)*sum(Y))/(n*sum(X*X)-sum(X)^2)
a <- mean(Y)-b*mean(X)
c(a, b) # the same as f$coefficients
```

### Derivation of the Solution (\*\*)




Remark.

: You can safely skip this part if you are yet to know
how to search for a minimum of a function of many variables
and what are partial derivatives.



Denote with:

\[
E(a,b)=\mathrm{SSR}(x\mapsto a+bx|\mathbf{x},\mathbf{y})
=\sum_{i=1}^n \left( a+bx_i - y_i \right) ^2.
\]

We seek the minimum of $E$ w.r.t. both $a,b$.



Theorem.

: If $E$ has a (local) minimum at $(a^*,b^*)$,
then its partial derivatives vanish therein,
i.e., $\partial E/\partial a(a^*, b^*) = 0$ and $\partial E/\partial b(a^*, b^*)=0$.





We have:

\[
E(a,b)   = \displaystyle\sum_{i=1}^n \left( a+bx_i - y_i \right) ^2.
\]

We need to compute the partial derivatives $\partial E/\partial a$ (derivative of $E$
w.r.t. variable $a$ -- all other terms treated as constants)
and $\partial E/\partial b$ (w.r.t. $b$).



Useful rules -- derivatives w.r.t. $a$ (denote $f'(a)=(f(a))'$):

- $(f(a)+g(a))'=f'(a)+g'(a)$ (derivative of sum is sum of derivatives)
- $(f(a) g(a))' = f'(a)g(a) + f(a)g'(a)$ (derivative of product)
- $(f(g(a)))' = f'(g(a)) g'(a)$ (chain rule)
- $(c)' = 0$ for any constant $c$ (expression not involving $a$)
- $(a^p)' = pa^{p-1}$ for any $p$
- in particular: $(c a^2+d)'=2ca$, $(ca)'=c$, $((ca+d)^2)'=2(ca+d)c$ (application of the above rules)


<!--\[
\begin{array}{rl}
E_a'(a,b)=&2\displaystyle\sum_{i=1}^n \left( a+bx_i - y_i \right) x_i \\
E_b'(a,b)=&2\displaystyle\sum_{i=1}^n \left( a+bx_i - y_i \right) \\
\end{array}
\]-->

We seek $a,b$ such that $\frac{\partial E}{\partial a}(a,b) = 0$
and $\frac{\partial E}{\partial b}(a,b)=0$.




\[
\left\{
\begin{array}{rl}
\frac{\partial E}{\partial a}(a,b)=&2\displaystyle\sum_{i=1}^n \left( a+bx_i - y_i \right) = 0\\
\frac{\partial E}{\partial b}(a,b)=&2\displaystyle\sum_{i=1}^n \left( a+bx_i - y_i \right) x_i = 0 \\
\end{array}
\right.
\]

This is a system of 2 linear equations. Easy.

Rearranging like back in the school days:

\[
\left\{
\begin{array}{rl}
b \displaystyle\sum_{i=1}^n x_i+ a n = & \displaystyle\sum_{i=1}^n  y_i  \\
b \displaystyle\sum_{i=1}^n x_i x_i + a \displaystyle\sum_{i=1}^n x_i = & \displaystyle\sum_{i=1}^n x_i y_i \\
\end{array}
\right.
\]




It is left as an exercise to show that the solution is:

\[
\left\{
\begin{array}{rl}
b^*  = & \dfrac{
n \displaystyle\sum_{i=1}^n x_i y_i - \displaystyle\sum_{i=1}^n  y_i \displaystyle\sum_{i=1}^n x_i
}{
n \displaystyle\sum_{i=1}^n x_i x_i -   \displaystyle\sum_{i=1}^n x_i\displaystyle\sum_{i=1}^n x_i
}\\
a^* = & \dfrac{1}{n}\displaystyle\sum_{i=1}^n  y_i - b^*  \dfrac{1}{n} \displaystyle\sum_{i=1}^n x_i  \\
\end{array}
\right.
\]



(we should additionally perform the second derivative test
to assure that this is the minimum of $E$ -- which is exactly the case though)

<!-- TODO rewrite nicely with means and correlation coefficients etc.

Pearson's r introduced in the next chapter though

(a <- cor(x,y)*sd(y)/sd(x))
(b <- mean(y)-a*mean(x))
-->








(\*\*) In the next chapter, we will introduce the notion of Pearson's
linear coefficient, $r$ (see `cor()` in R). It might be shown that
$a$ and $b$ can also be rewritten as:

```{r ab_coef_cor}
(b <- cor(X,Y)*sd(Y)/sd(X))
(a <- mean(Y)-b*mean(X))
```

<!-- TODO show optim() -->



## Exercises in R

### The Anscombe Quartet

Here is a famous illustrative example proposed by
the statistician Francis Anscombe in the early 1970s.

```{r anscombe}
print(anscombe) # `anscombe` is a built-in object
```

What we see above is a single data frame
that encodes four separate datasets:
`anscombe$x1` and `anscombe$y1` define the first pair of variables,
`anscombe$x2` and `anscombe$y2` define the second pair and so forth.





{ BEGIN exercise }
Split the above data (manually) into four data frames
`ans1`, ..., `ans4` with columns `x` and `y`.

For example, `ans1` should look like:

```{r anscombe1,echo=-(1:4)}
ans1 <- data.frame(x=anscombe$x1, y=anscombe$y1)
ans2 <- data.frame(x=anscombe$x2, y=anscombe$y2)
ans3 <- data.frame(x=anscombe$x3, y=anscombe$y3)
ans4 <- data.frame(x=anscombe$x4, y=anscombe$y4)
print(ans1)
```
{ END exercise }





{ BEGIN solution }

```{r anscombe1apply, eval=FALSE}
<<anscombe1>>
```

{ END solution }



{ BEGIN exercise }
Compute the mean of each `x` and `y` variable.
{ END exercise }

{ BEGIN solution }

```{r anscombe_mean}
mean(ans1$x) # individual column
mean(ans1$y) # individual column
sapply(ans2, mean) # all columns in ans2
sapply(anscombe, mean) # all columns in the full anscombe dataset
```


*Comment: This is really interesting, all the means of `x` columns
as well as the means of `y`s are almost identical.*
{ END solution }




{ BEGIN exercise }
Compute the standard deviation of each `x` and `y` variable.
{ END exercise }

{ BEGIN solution }

The solution is similar to the previous one, just replace `mean` with `sd`.
Here, to learn something new, we will use the `knitr::kable()` function
that pretty-prints a given matrix or data frame:

```{r anscombe_sd}
results <- sapply(anscombe, sd)
knitr::kable(results, col.names="standard deviation")
```


*Comment: This is even more interesting, because the numbers agree up to 2 decimal digits.*

{ END solution }

<!--
TODO: cor() covered later
Compute the Pearson linear correlation coefficient
for each of the four pairs of `x` and `y`.

To recall, this can be done with the `cor()` function.
What is the Pearson $r$ coefficient is explained in Lecture 3.

```{r}
cor(ans1$x, ans1$y)
cor(ans2$x, ans2$y)
cor(ans3$x, ans3$y)
cor(ans4$x, ans4$y)
```


*Comment: We get a feeling that we're being tricked by Anscombe...
All the variables are highly correlated, and the correlation
coefficients are more or less the same.*
-->



{ BEGIN exercise }
Fit a simple linear regression model for each data set.
Draw the scatter plots again (`plot()`)
and add the regression lines (`lines()` or `abline()`).
{ END exercise }

{ BEGIN solution }

To recall, this can be done with the `lm()` function
explained in Lecture 2.

At this point we should already have become lazy -- the tasks are very
repetitious. Let's automate them by writing a single function
that does all the above for any data set:

```{r anscombe_fit}
fit_models <- function(ans) {
    # ans is a data frame with columns x and y
    f <- lm(y~x, data=ans) # fit linear model
    print(f$coefficients) # estimated coefficients
    plot(ans$x, ans$y) # scatter plot
    abline(f, col="red") # regression line
    return(f)
}
```

Now we can apply it on the four particular examples.

```{r anscombe-fit-apply,fig.height=6,fig.cap="Fitted regression lines for the Anscombe quartet"}
par(mfrow=c(2, 2)) # four plots on 1 figure (2x2 grid)
f1 <- fit_models(ans1)
f2 <- fit_models(ans2)
f3 <- fit_models(ans3)
f4 <- fit_models(ans4)
```

*Comment: All the estimated models are virtually the same,
the regression lines are
$y=0.5x+3$, compare Figure \@ref(fig:anscombe-fit-apply).*

{ END solution }





{ BEGIN exercise }
Create scatter plots of the residuals (predicted $\hat{y}_i$ minus
true $y_i$) as a function of the predicted $\hat{y}_i=f(x_i)$ for every
$i=1,\dots,11$.
{ END exercise }

{ BEGIN solution }

To recall, the model predictions can be generated by (amongst others)
calling the `predict()` function.


```{r anscombe_resid}
y_pred1 <- f1$fitted.values # predict(f1, ans1)
y_pred2 <- f2$fitted.values # predict(f2, ans2)
y_pred3 <- f3$fitted.values # predict(f3, ans3)
y_pred4 <- f4$fitted.values # predict(f4, ans4)
```

Plots of residuals as a function of the predicted (fitted) values
are given in Figure \@ref(fig:anscombe-resid-plot).

```{r anscombe-resid-plot,fig.height=6,fig.cap="Residuals vs. fitted values for the regression lines fitted to the Anscombe quartet"}
par(mfrow=c(2, 2)) # four plots on 1 figure (2x2 grid)
plot(y_pred1, y_pred1-ans1$y)
plot(y_pred2, y_pred2-ans2$y)
plot(y_pred3, y_pred3-ans3$y)
plot(y_pred4, y_pred4-ans4$y)
```

*Comment: Ideally, the residuals shouldn't be correlated
with the predicted values -- they should "oscillate" randomly
around 0.     This is only the case of the first dataset.
All the other cases are "alarming" in the sense that
they suggest that the obtained models are "suspicious"
(perhaps data cleansing is needed or a linear model is not at all appropriate).*

{ END solution }



{ BEGIN exercise }
Draw conclusions (in your own words).
{ END exercise }

{ BEGIN solution }
We're being taught a lesson here: don't perform data analysis tasks
automatically, don't look at bare numbers only, visualise your data first!
{ END solution }




{ BEGIN exercise }
Read more about Anscombe's quartet at https://en.wikipedia.org/wiki/Anscombe%27s_quartet
{ END exercise }





## Outro

### Remarks


In supervised learning, with each input point,
there's an associated reference output value.

Learning a model = constructing a function that approximates
(minimising some error measure) the given data.

Regression = the output variable $Y$ is continuous.

We studied linear models with a single independent variable based on
the least squares (SSR) fit.

In the next part we will extend this setting to the case
of many variables, i.e., $p>1$, called multiple regression.




### Further Reading


Recommended further reading:  [@islr: Chapters 1, 2 and 3]

Other: [@esl: Chapter 1, Sections 3.2 and 3.3]