```{r setup, include=FALSE}
library("knitr")
library("dplyr")
library("plotly")
library("ggplot2")
library("tidyverse")
```

## 1.1 Linear regression: a basic model for parameter estimation

Here we look more in detail at how the estimation of model parameters works. We have a generating model of the following form:

$$
y = 1.25 \cdot x
$$

The **true** $\theta$ (model parameter: slope) is therefore $1.25$

```{r lin_reg}
lin_reg <- function(x) 1.25*x
```

```{r}
x = 1
y = lin_reg(x)
print(paste("y as function of x=1 in the model above:",y))
```

Now apply the function to a bunch of data

```{r}
x = seq(-5,+5,0.5) ## independent variable
y = lin_reg(x) ## dependent variable

data.frame("x" = x, "y" = y)
```

```{r echo=FALSE}
p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x))
p <- p + stat_function(fun = lin_reg) + xlim(-5,6)
p <- p + theme(axis.title.y = element_text(angle=0, vjust = 0.5))
p
```

### The loss function

For linear regression (simple or multiple, as long as $n > p$), the least squares method can be used, where the residual sum of squares is minimised through differentiation of vector and matrix expressions (linear algebra $\rightarrow$ *normal equations*).

$$
\mathbf{y} = \mathbf{Xb}+\mathbf{e} \rightarrow \mathbf{X’y} = \mathbf{X’Xb} \rightarrow \mathbf{X’Xb} = \mathbf{X’y} \rightarrow \mathbf{b} = (\mathbf{X’X})^{-1} \cdot \mathbf{X’y}
$$

$$
b = \frac{cov(x,y)}{var(x)}
$$

$$
\mu = avg(y) - b \cdot avg(x)
$$

However, from the **perspective of machine learning** a different approach is taken. First, a **loss function** is chosen: a common choice for (multiple) linear regression is the **normalised squared error function**:

$$
J(\beta) = \frac{1}{2n} \sum_{i=1}^n \left( \beta_i X_i -y_i \right)^2
$$

```{r}
loss_function <- function(x,beta) {
  
  n = length(x)
  y = lin_reg(x)
  normalised_squared_error = sum((y - beta*x)^2)/(2*n)
  
  return(normalised_squared_error)
}
```

We then calculate the loss function for different values of the parameter(s) to estimate. We take 11 datapoints (from 0 to 10) for our linear regression model and try different values for beta:

```{r}
x <- seq(0,10,0.25)
beta <- seq(0.25,2.25,0.05)

print("The data points:")
print(x)

cat("\nand the values for the parameter beta:\n")
print(beta)
```

Now we calculate the loss function for all different values of the parameter $\beta$ to be estimated:

```{r}
cost <- sapply(beta, function(z) loss_function(x,z))
res <- data.frame("beta" = beta, "loss" = cost)
beta_min = res[which.min(res$loss),"beta"]
print(res)
```

```{r}
print(paste("parameter value for which the cost function is minimised:", beta_min))
```

```{r pressure, echo=FALSE}
plot(beta,cost, type="l",xlab = "Values for parameter beta", ylab = "Value for loss function")
```

Note: sometimes the term **loss function** is used to define the individual loss (each single record, i.e. $(y_i-\hat{y_i})^2$) and the term **cost function** is used for the sum of individual losses

### Random dataset

```{r}
# Construct sample data set
set.seed(132)
x <- sort(rnorm(30))
y <- sort(rnorm(30)/sample(3,1))

plot(x,y)
```

We can model this data using a **simple linear regression model** of the following form:

$$
y = \beta_0 + \beta_1 x + e
$$

To find the linear regression that "best" fits the data, we now try a range of values for the parameters $\beta_0$ and $\beta_1$, and modify the **loss function** accordingly:

```{r}
beta0 = seq(-1,1,0.1)
beta1 = seq(0.1,3,0.1)
parameters <- expand_grid(beta0,beta1)
```

```{r}
loss_function <- function(x,y,beta0,beta1) {
  
  n = length(x)
  y_hat = beta0 + beta1*x
  normalised_squared_error = sum((y - y_hat)^2)/(2*n)
  
  return(normalised_squared_error)
}
```

```{r}
res <- parameters |> group_by(beta0,beta1) |> 
  mutate(cost = loss_function(x=x, y=y, beta0=beta0, beta1=beta1))
print(res)
```

```{r}
print(paste("The minimum value for the cost function is", min(res$cost)))
```

We can plot the cost function to appraise it visually:

```{r}
library("plotly")

plot_ly(x = res$beta0, y=res$beta1, z=res$cost, type = "scatter3d", mode = "markers", color = res$cost) %>%
 layout(
     title = "cost function with 2 parameters",
      scene = list(
          xaxis = list(title = "beta0"),
          yaxis = list(title = "beta1"),
          zaxis = list(title = "cost")
      )
 )
```

```{r}
cat("Estimated coefficents 
    (values of parameters that minimise the cost function:\n")

beta_min = res[which.min(res$cost),c("beta0","beta1")]
beta_min
```

```{r}
plot(x,y)
abline(beta_min$beta0, beta_min$beta1, col="red")
```

#### Compare with least squares

```{r}
X = matrix(c(rep(1,length(x)), x), ncol = 2)
beta_coeffs = solve(t(X)%*%X) %*% t(X)%*%y
print(beta_coeffs)
```

```{r}
plot(x,y)
abline(beta_min$beta0, beta_min$beta1, col="red")
abline(beta_coeffs[1], beta_coeffs[2], col="blue")
```

### Exercise 1.1

Try to estimate your own model coefficient:

1.  Create your random data (**dare try with more than two parameters?**):

```{r results="asis"}
x1 <- NULL
# x2 <- NULL
y <- NULL
```

2.  Define your own loss function:

```{r results="asis"}
loss_function.1 <- function() {}
```

3.  Choose a set of values for the $\beta$'s to be tested (and combine them in a grid of unique value combinations):

```{r results="asis"}
beta0 = c()
beta1 = c()
beta2 = c()
## etc.
```

4.  Calculate the values for the loss function and get the coefficients corresponding to the minimum loss:

```{r results="asis"}
cost = c()
```

5.  [optional] Plot the cost function vs values of the parameter(s):

```{r results="asis"}
# plot()
```

6.  [optional] Compare your results with those from least squares:

```{r results="asis"}
# plot()
```

## 1.2 Linear regression: measuring performance

In practice, we are not going to manually minimise the loss function to estimate model parameters for our predictive machine: instead, higher-level *R* functions are used, like `lm()`.

An important aspect of predictive statistics is to measure the performance of the developed predictive model (predictive machine).

Let's start by using an example dataset from base R: the ChickWeight dataset, with weight and age of chicks.

```{r}
data(ChickWeight) ##
dataset <- rename(ChickWeight, y = weight, x = Time) %>% select(y,x)
head(dataset)
```

y: weight (grams) x: time (days)

```{r}
summary(dataset)
```

We now fit a simple linear regression model:

$$
y = \mu + \beta \cdot x + e
$$

```{r}
fit <- lm(y ~ x, data = dataset)
coef(fit)
```

[Question: how do we interpret the estimated model coefficients?]{style="color:red"}

```{r}
ggplot(dataset, aes(x = x, y = y)) + geom_jitter() + 
  geom_smooth(method = "lm", se = FALSE)
```

We now have all the ingredients to obtain predictions: either by explicitly using the estimated coefficients:

```{r}
predictions <- dataset$x*coef(fit)[2] + coef(fit)[1]
head(predictions)
```

or by using the *R* `predict()` function:

```{r}
?predict
predictions <- predict(fit, newdata = dataset)
head(predictions)
```

The two approaches are obviously equivalent

```{r}
concordance <- predict(fit, newdata = dataset) == dataset$x*coef(fit)[2] + 
  coef(fit)[1]
print(sum(concordance)/length(predictions))
```

The predict function is more flexible and can for instance also give us a confidence interval for predictions:

```{r}
predict(fit, newdata = dataset, interval = "confidence") %>%
  head()
```

Finally, we can plot predictions against observations:

```{r}
dataset$predictions <- predictions
dataset <- dataset |>
  mutate(error = abs(predictions-y)/y)

ggplot(data = dataset, aes(y,predictions)) + geom_jitter(aes(color=error)) + xlab("observed values") +
  scale_color_gradient(low="black", high="red")
#abline(fit)
```

Besides visualizing how predictions relate to observations, we need also to measure (quantify) the predictive performance of the model.

Several metrics exist for regression problems. Here we list a few of the most commonly used.

1.  **MSE** (mean squared error)

```{r}
mse <- function(y,y_hat) {
  
  n = length(y)
  se = sum((y-y_hat)^2)
  mse = se/n
  
  return(mse)
}

error = mse(y = dataset$y, y_hat = predictions)
error
```

The MSE is **`r round(error,3)`**.

2.  **RMSE** (root mean squared error): this is on the same scale as the target variable

```{r}
rmse = sqrt(error)
rmse
```

The RMSE is **`r round(rmse,3)`**.

3.  **MAE** (mean absolute error)

```{r}
mae <- function(y,y_hat) {
  
  n = length(y)
  se = sum(abs(y-y_hat))
  mae = se/n
  
  return(mae)
}

error = mae(y = dataset$y, y_hat = predictions)
error
```

The MAE is **`r round(error,3)`**.

**RMSE** and **MAE** can also be expressed as percentage of the response mean value:

```{r}
avg = mean(dataset$y)

print(paste("RMSE is", round((rmse/avg)*100,2), "% of the average chick weight"))
print(paste("MAE is", round((error/avg)*100,2), "% of the average chick weight"))
```

------------------------------------------------------------------------

The we have correlations:

4.  **Pearson's linear** correlation coefficient
5.  **Spearman's rank** correlation coefficient

```{r}
r_pearson = cor(dataset$y, predictions, method = "pearson")
r_spearman = cor(dataset$y, predictions, method = "spearman")
```

```{r}
print(r_pearson)
print(r_spearman)
```

------------------------------------------------------------------------

[**Food for thought**: what if the relationship between the response variable and the predictors is not linear?]{style="color:green"}

```{r}
ggplot(dataset, aes(x = x, y = y)) + geom_jitter() + geom_smooth(method = "lm", formula = y ~ poly(x,2), se = FALSE)
```

------------------------------------------------------------------------

## Exercise 1.2

Choose a built-in dataset, fit a linear model and measure the accuracy of predictions:

1.  you can choose one of the many built-in R datasets, using the function `data()`:
2.  or you can generate a dataset (e.g. sampling from a Gaussian distribution)

#### Part 1

```{r}
data()
```

```{r results="asis"}
# data(NULL) ## choose a dataset
# dataset <- rename() %>% select()
# head(dataset)
```

2.  Fit a linear model

```{r results="asis"}
# fit <- lm()
```

#### Part 2

3.  Obtain predictions

```{r results="asis"}
# predictions <- predict()
```

4.  Plot observations vs predictions

```{r}
# p <- ggplot()
```

5.  Choose a metric to measure the accuracy of predictions (performance)

```{r results="asis"}
# metric = 
```
