4-SimpleRegression.qmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r, include=FALSE}
set.seed(42)
```

# Simple linear regression

In this chapter, you will:

-   perform simple linear regression for numeric and categorical predictors
-   interpret regression outputs
-   check the residuals of regression models

## Maximum Likelihood Estimator

$likelihood = P(D|model, parameter)$

The likelihood is the probability to observe the Data given a certain model (which is described by its parameter).

It is an approach to optimize a model/parameter to find the set of parameters that describes best the observed data.

A simple example: we want to estimate the average of a random vector and we assume that our model is a normal distribution (so we assume that the data originated from a normal distribution). We want to optimize the two parameters that describe the normal distribution: the mean, and the standard deviation:

Let's create the "oberved" vector:

```{r}
Xobs = rnorm(n = 100, mean = 0, sd = 1.0)
```

Now, we assume that the mean and the standard deviation of the vector is unknown and we want to find them. Let's write the likelihood function (in reality, we use logLikelihood for numerical reasons):

```{r}
likelihood = function(par) { # we give two parameters, mean and sd
  # calculate for each observation to observe the data given our model:
  lls = dnorm(Xobs, mean = par[1], sd = par[2], log = TRUE) 
  # we use the logLikilihood for numerical reasons
  return(sum(lls))
}
```

Let's guess some values for the mean and sd and see what the function returns

```{r}
likelihood(c(0, 0.2))
```

This number alone doens't tell us much, we have com compare many different guesses to see which combination of parameters gives is the most plausible. To make it easier to understand, we will fix the sd in 1 and plot how the likelihood changes with changing the mean:

```{r}
likelihood_mean = function(mean.par){
  likelihood(c(mean.par, 1.0))}

plot(seq(-5, 5.0, length.out = 100), 
     sapply(seq(-5, 5.0, length.out = 100), likelihood_mean),
     xlab = 'Different mean values', 
     ylab = "negative logLikelihood")
```

The optimum (highest peak) is at 0, which corresponds to the mean we used to sample Xobs. We could do the same with the sd - fix the mean and vary the sd to find it.

However it is tedious to try all values manually to find the best value, especially if we have to optimize several parameters. For that, we can use an optimizer in R which finds for us the best set of parameters:

```{r}
opt = optim(c(0.0, 1.0), fn = function(par) -likelihood(par), hessian = TRUE )
```

We can use the shape of the likelihood to calculate standard errors for our estimates:

```{r}
st_errors = sqrt(diag(solve(opt$hessian)))
```

With that we can calculate the confidence interval for our estimates. When the estimator is repeatedly used, 95% of the calculated confidence intervals will include the true value!

```{r}
cbind(opt$par-1.96*st_errors, opt$par+1.96*st_errors)
```

In short, if we would do a t.test for our Xobs (to test whether the mean is stat. significant different from zero), the test would be non significant, and a strong indicator for that is when the 0 is within the confidence interval. Let's compare our CI to the one calculated by the t-test:

```{r}
t.test(Xobs)
```

The results are lmost the same! The t-test also calculates the MLE to get the standard error and the confidence interval.

## The theory of linear regression

If we want to test for an association between two continuous variables, we can calculate the correlation between the two - and with the cor.test function we can test even for significance. But, the correlation doesn't report the magnitude, the strength, of the effect:

```{r}
X = runif(100)
par(mfrow = c(1,1))
plot(X, 0.5*X, ylim = c(0, 1), type = "p", pch = 15, col = "red", xlab = "X", ylab = "Y")
points(X, 1.0*X, ylim = c(0, 1), type = "p", pch = 15, col = "blue", xlab = "X", ylab = "Y")
cor(X, 0.5*X)
cor(X, 1.0*X)
```

Both have a correlation factor of 1.0! But we see clearly that the blue line has an stronger effect on Y then the red line.

**Solution: Linear regression models**

They describe the relationship between a dependent variable and one or more explanatory variables:

$$
y = a_0 +a_1*x
$$

(x = explanatory variable; y = dependent variable; a~0~=intercept; a~1~= slope)

Examples:

```{r}
plot(X, 0.5*X, ylim = c(0, 1), type = "p", pch = 16, col = "black", xlab = "X", ylab = "Y", lwd = 1.5)
points(X, 0.5*X, col = "red", type = "l", lwd = 1.5)
points(X, 1.0*X, ylim = c(0, 1), type = "p", pch = 16, col = "black", xlab = "X", ylab = "Y", lwd = 1.5)
points(X, 1.0*X, ylim = c(0, 1), type = "l", pch = 16, col = "blue", xlab = "X", ylab = "Y", lwd = 1.5)
legend("topleft", col = c("red", "blue"), lty = 1,legend = c('Y = 0.5*X+0', 'Y = 1.0**X+0'))
```

We can get the parameters (slope and intercept) with the MLE. However, we need first to make another assumptions, usually the model line doesn't perfectly the data because there is an observational error on Y, so the points scatter around the line:

```{r}
plot(X, 0.5*X+rnorm(100, sd = 0.05), ylim = c(0, 1), type = "p", pch = 16, col = "black", xlab = "X", ylab = "Y", lwd = 1.5)
points(X, 0.5*X, col = "red", type = "l", lwd = 1.5)
points(X, 1.0*X+rnorm(100, sd = 0.05), ylim = c(0, 1), type = "p", pch = 16, col = "black", xlab = "X", ylab = "Y", lwd = 1.5)
points(X, 1.0*X, ylim = c(0, 1), type = "l", pch = 16, col = "blue", xlab = "X", ylab = "Y", lwd = 1.5)
legend("topleft", col = c("red", "blue"), lty = 1,legend = c('Y = 0.5*X+0', 'Y = 1.0**X+0'))
```

And the model extends to:

$$
y = a_0 + a_1*x +\epsilon, \epsilon \sim N(0, sd)
$$

Which we can also rewrite into:

$$
y = N(a_0 + a_1*x, sd)
$$

Which is very similar to our previous MLE, right? The only difference is now that the mean depends now on x, let's optimize it again:

```{r}
Xobs = rnorm(100, sd = 1.0)
Y = Xobs + rnorm(100,sd = 0.2)
likelihood = function(par) { # three parameters now
  lls = dnorm(Y, mean = Xobs*par[2]+par[1], sd = par[3], log = TRUE) # calculate for each observation the probability to observe the data given our model
  # we use the logLikilihood because of numerical reasons
  return(sum(lls))
}

likelihood(c(0, 0, 0.2))
opt = optim(c(0.0, 0.0, 1.0), fn = function(par) -likelihood(par), hessian = TRUE )

opt$par
```

Our true parameters are 0.0 for the intercept, 1.0 for the slope, and 0.22 for the sd of the observational error.

Now, we want to test whether the effect (slope) is statistically significant different from 0:

1.  calculate standard error
2.  calculate t-statistic
3.  calculate p-value

```{r}
st_errors = sqrt(diag(solve(opt$hessian)))
st_errors[2]
t_statistic = opt$par[2] / st_errors[2]
pt(t_statistic, df = 100-3, lower.tail = FALSE)*2
```

The p-value is smaller than $\alpha$, so the effect is significant! However, it would be tedious to do that always by hand, and because it is probably one of the most used analysis, there's a function for it in R:

```{r}
model = lm(Y~Xobs) # 1. Get estimates, MLE
model
summary(model) # 2. Calculate standard errors, CI, and p-values
```

## Understanding the linear regression {#sec-lm}

![](resources/lm.png)

Besides the MLE, there are also several tests in a regression. The most important are

-   significance of each parameter. t-test: H~0~ = variable has no effect, that means the estimator for the parameter is 0

-   significance of the model. F-test: H~0~ = none of the explanatory variables has an effect, that means all estimators are 0

**Example:**

```{r}
pairs(airquality)
# first think about what is explanatory / predictor 
# and what is the dependent variable (e.g. in Ozone and Temp)

# par(mfrow = c(1, 1))
plot(Ozone ~ Temp, data = airquality)
fit1 = lm(Ozone ~ Temp, data = airquality)
summary(fit1)
# gives a negative values for the intercept = negative Ozone levels when Temp = 0
# this does not make sense (>extrapolation)

# we can also fit a model without intercept, 
# without means: intercept = 0; y = a*x 
# although this doesn't make much sense here
fit2 = lm(Ozone ~ Temp - 1, data = airquality)
summary(fit2)

plot(Ozone ~ Temp, data = airquality, xlim = c(0,100), ylim = c(-150, 150))
abline(fit1, col = "green")
abline(fit2, col = "red", lty = 2)

# there is no need to check normality of Ozone
hist(airquality$Ozone) # this is not normal, and that's no problem !
```

### Model diagnostics

The regression optimizes the parameters under the condition that the model is correct (e.g. there is really a linear relationship). If the model is not specified correctly, the parameter values are still estimated - to the best of the model's ability, but the result will be misleading, e.g. p-values and effect sizes

What could be wrong:

-   the distribution (e.g. error not normal)
-   the shape of the relationship between explanatory variable and dependent variable (e.g., could be non-linear)

The model's assumptions must always be checked!

We can check the model by looking at the residuals (which are observed - predicted values) which should be normally distributed and should show no patterns:

```{r}
X = runif(50)
Y = X + rnorm(50, sd = 0.2)
fit = lm(Y~X)
par(mfrow = c(1, 3))
plot(X, Y)
abline(fit, col = "red")
plot(X, Y - predict(fit), ylab = "Residuals")
abline(h = 0, col = "red")
hist(Y - predict(fit), main = "", xlab = "Residuals")
par(mfrow = c(1,1))
```

The residuals should match the model assumptions. For linear regression:

-   normal distribution
-   constant variance
-   independence of the data points

Example:

```{r}
fit1 = lm(Ozone~Temp, data = airquality)
residuals(fit1)
hist(residuals(fit1))
# residuals are not normally distributed
# we do not use a test for this, but instead look at the residuals visually

# let's plot residuals versus predictor
plot(airquality$Temp[!is.na(airquality$Ozone)], residuals(fit1))

# model checking plots
oldpar= par(mfrow = c(2,2))
plot(fit1)
par(oldpar)
#> there's a pattern in the residuals > the model does not fit very well!
```

## Linear regression with a categorical variable

We can also use categorical variables as an explanatory variable:

```{r}
m = lm(weight~group, data = PlantGrowth)
summary(m)
```

The lm estimates an effect/intercept for each level in the categorical variable. The first level of the categorical variable is used as a reference, i.e. the true effect for grouptrt1 is Intercept+grouptrt1 = `r 5.0320+(-0.3710)` and grouptrt2 is `r 5.0302+0.4940`. Moreover, the lm tests for a difference of the reference to the other levels. So with this model we know whether the control is significant different from treatment 1 and 2 but we cannot say anything about the difference between trt1 and trt2.

If we are interested in testing trt1 vs trt2 we can, for example, change the reference level of our variable:

```{r}
tmp = PlantGrowth
tmp$group = relevel(tmp$group, ref = "trt1")
m = lm(weight~group, data = tmp)
summary(m)
```

Another example:

```{r}
library(effects)
library(jtools)

summary(chickwts)

plot(weight ~ feed, chickwts)
fit4 = lm(weight ~ feed, chickwts)

summary(fit4)
anova(fit4) #get overall effect of feeding treatment

plot(allEffects(fit4))
plot(allEffects(fit4, partial.residuals = T))
effect_plot(fit4, pred = feed, interval = TRUE, plot.points = F)

old.par = par(mfrow = c(2, 2))
plot(fit4)
par(old.par)

boxplot(residuals(fit4) ~ chickwts$feed)
```

### Linear regression with a quadratic term

What does simple linear regression mean?

-   simple = one predictor!

-   linear = linear in the parameters

If we add a quadratic term, this is still a linear combination: this is called polynomial.

$$ a0 + a1 * x + a2 * x^2 $$

Le'ts model Ozone as a function of temperature and it's quadratic term:

```{r}
fit3 = lm(Ozone ~ Temp + I(Temp^2), data = airquality)
summary(fit3)

oldpar= par(mfrow = c(2,2))
plot(fit3)
par(oldpar)
```

Residual vs. fitted looks okay, but Outliers are still there, and additionally too wide. But for now, let's plot prediction with uncertainty (plot line plus confidence interval):

**OBS:** If the relationship between x and y is not linear, we cannot use abline instead we predict values of x for different values of y based on the model:

```{r}
plot(Ozone ~ Temp, data = airquality)
newDat = data.frame(Temp = 55:100)
predictions = predict(fit3, newdata = newDat, se.fit = T)
# and plot these into our figure:
lines(newDat$Temp, predictions$fit, col= "red")
# let's also plot the confidence intervals:
lines(newDat$Temp, predictions$fit + 1.96*predictions$se.fit, col= "red", lty = 2)
lines(newDat$Temp, predictions$fit - 1.96*predictions$se.fit, col= "red", lty = 2)

# add a polygon (shading for confidence interval)
x = c(newDat$Temp, rev(newDat$Temp))
y = c(predictions$fit - 1.96*predictions$se.fit, 
      rev(predictions$fit + 1.96*predictions$se.fit))

polygon(x,y, col="#99009922", border = F )
```

Alternative: use package effects

```{r}
library(effects)
plot(effect("Temp",fit3))
plot(effect("Temp",fit3,partial.residuals = T) ) 
#to check patterns in residuals (plots measurements and partial residuals)
```

Another alternative:

```{r}
# or jtools package
library(jtools)
effect_plot(fit3, pred = Temp, interval = TRUE, plot.points = TRUE)
```