Linear models

Example 1: Linear regression with one explanatory variable

Definitions

Consider $n = 7$ farm fields. Each one has two measurements:
- $y_{i}$ - the crop yield
- $x_{i}$ - amount of fertilizer used (a covariate)

# Amount of fertilizer - increases by 100 each field
x=seq(100,700,by=100) 
# Crop yield
y=c(40,45,50,65,70,70,80)

Simple model

We can analyse these with a linear model:$$E(y)=\beta_0 + \beta_1x$$
In this model, we are assuming that:

$β_{0}$ is intercept
- (amount of yield in theoretical field with zero fertilizer)
$β_{1}$ is the slope
- (unit of crop yield per unit of fertilizer).

More complex model

Let's add the assumption that observation $y_{i}$ has independent and identically distributed (iid) Gaussian noise ^[1] ( $ϵ_{i}$ ) in addition to it's mean value, representing Aleatory uncertainty:

y_{i} = E [y_{i}] + ϵ_{i}

y_{i} = β_{0} + β_{1} x + ϵ_{i}

where $ϵ_{i} \sim N (0, σ^{2})$

Implementation in R

model=lm(y~x)
summary(model)

!lecture_notes1_week1, p.20
The sections of the output from R are:

Call

A reminder of the call that you used to generate the object - i.e. the code we used to generate the model

Residuals

This gives a summary of the distribution of residuals, i.e. the difference between our model line and each point. Since this only has 7 points it gives one for each, but for larger datasets it will print the quantiles. Residuals formula:

r_{i} = y_{i} - {\hat{b}}_{0} - {\hat{b}}_{1} x_{1} = y - {\hat{y}}_{i}

I think the hats here represent the model values, so this is another way of saying (real value - predicted value).

Coefficients

This summarizes the estimated parameters of the model in a table for each parameter. In our model there are only two, (Intercept) ( $β_{0}$ ) and x (slope for $x$ , $β_{1}$ ).

Estimate: The least squares estimate for model parameters
Std. Error: The estimated standard error associated with parameter. This is a measure of uncertainty related to the parameter estimates.
t value: Lists the t-values for parameter estimates, given by $t_{j} = {\hat{β}}_{j} / s_{{\hat{β}}_{j}}$ . A standardized estimate of how far each parameter is from zero. This is not the same as the t-value given by the t-test but it's a bit similar.
Pr(>|t|): the p-value corresponding to the given t-value, can be used for hypothesis testing. The probability of getting a t-value this extreme.
- The null hypothesis in linear models is that a particular covariate is not important for explaining the data
- If the null hypothesis is true, the $t$ value should be consistent with what would be drawn from a t-distribution.
- If we have a very low p-value, it means that this t-value is not consistent with what is expected, and that there is a significant effect that this variable has on the model. They add *s and .s after these showing how significant they are.

Signif. codes

A legend explaining the * and . symbols associated with the p-values in the coefficients table.

Residual standard error

The estimate for the standard deviation of the residuals. This is denoted with $\hat{σ}$ in maths.

\hat{σ} = \sqrt{\frac{\sum_{i = 1}^{n} r_{i}^{2}}{n - p}} = \sqrt{\frac{\sum_{i = 1}^{n} (y_{i} - {\hat{b}}_{0} - {\hat{b}}_{1} x_{1})^{2}}{n - p}}

Multiple R-squared

Estimate of the proportion of the variance in the data explained by the regression. It's defined as:

r^{2} = 1 - \frac{\sum_{i = 1}^{n} r_{i}^{2}}{\sum_{i = 1}^{n} (y_{i} - \bar{y})^{2}}

$\bar{y}$ is the mean of all $y_{i}$ (sample mean).
I assume that the $r_{i}$ in the numerator is the residual at value $i$ . Don't think we need to really understand this.

Note that this is a biased Estimator, it increases every time a new parameter is added to the data.

Adjusted R-squared

ANother estimator for the proportion of variance in the data explained by the regression. The problem with R-squared is that it always increases when you add a new predictor to the model, no matter if the new predictor is useful or not.

r_{adj}^{2} = 1 - \frac{\sum_{i = 1}^{n} r_{i}^{2} / (n - p)}{\sum_{i = 1}^{n} (y_{i} - \bar{y})^{2} / (n - 1)}

Where $n$ is the number of observations and $p$ is the number of parameters.

This is an unbiased estimator.

F-statistic

The final line gives the F-statistic and its p-value.

H0: The data was generated by a model with an intercept term only
H1: The fitted model generated the data
By our null hypothesis, the model is doing nothing and the data is totally random, and our alternative hypothesis is the likelihood that our data was generated by this model.

So this statistic tells us whether the model as a whole is useful or not.
It's based on the residual sum of squares or RSS:

The sum of all the residuals from each data point squared.

RSS = \sum_{i = 1}^{n} r_{i}^{2}

This is used for calculating the F-statistic in a linear model.

F = \frac{{RSS}_{0} - {RSS}_{1} / (p - 1)}{{RSS}_{1} / (n - p)}

${RSS}_{0}$ : RSS calculated from the null model
${RSS}_{1}$ : RSS calculated from the given linear model
$p$ : Number of parameters in linear model.

The F-statistic follows a given distribution so we can get a p-value of how likely our F statistic is to be in the distribution.

Confidence Intervals

We can look at the Confidence intervals with the confint() function in R - by default the 95% ones.

!lecture_notes1_week1, p.21

Conclusion

Looks like our model is pretty good and amount of fertilizer is a significant explanatory variable.

Yield increases 0.069 units per unit of fertilizer (slope)
Model has a high adj. R-squared (0.946)

Visualizing the regression

Pretty straightforward, you can plot the regression line with abline().
!lecture_notes1_week1, p.22

Confidence intervals and predictions

We can also use the predict() function to get confidence intervals, then visualize these upper and lower confidence intervals.

y_ci=predict(
	model,
	se.fit=TRUE,
	interval="confidence",
	level=0.95
)
{plot(x,y,xlab="Fertilizer ",ylab="Yield") 
abline(model,col="red")
points(x,y_ci$fit[,"fit"],col="red") 
lines(x,y_ci$fit[,"lwr"], col="green") 
lines(x,y_ci$fit[,"upr"], col="green")}

!lecture_notes1_week1, p.23
Then you can also add uncertainty estimate for new data that could be added (prediction) to the model with this line:

y_pi = predict(
	model, 
	interval="prediction", 
	level=0.95
)
matlines(x, y_pi[,c("lwr","upr")], col="black",lty=2)

!lecture_notes1_week1, p.23
Prediction intervals are wider than confidence intervals since uncertainty is higher with a single value, while confidence intervals work around the mean.

Residuals

You can also plot the residuals as histogram or scatterplot:

resid=residuals(model) 
{par(mfrow=c(1,2)) 
hist(resid,10) 
plot(fitted.values(model),resid) 
abline(h=0)}

A good model will have a random distribution of residuals centered around 0 with no discernible pattern. It seems it is debated whether them being normally distributed is important (?)

Example 2: Linear regression with multiple explanatory variables

This example looks at photosynthesis in small ponds treated with different amounts of nutients and amounts of plankton.

The lecture nots are pretty clear for this example. You start by plotting scatterplots of how all variables effect each other.

Next, looking at the summary() function can let you filter out non-significant explanatory variables. Removing them makes the model better (increases adj r squared)

An anova() analysis (ANOVA (Maths)) comparing the two models shows that keeping the extra variables didn't add anything else, but I'm not sure why exactly (we will learn about this next week I think).

Then looking at some summary plots, we see some problems (explained in the lecture notes). They address these by using a 2nd order polynomial to fit the sum of the two plankton variables (instead of first order).

Categorical covariates in a linear model

This was discussed in week 2:
Categorical covariates

See Central Limit Theorem for an explanation of this and why we assume Gaussian distributions. ↩︎