Linear models form an important starting point for models (more) dedicated to spatial and/or temporal data. The basic linear model can be seen as a multiple linear regression model:
\(y = (y_1, ..., y_n)'\) the vector with \(n\) observations,
\(X = [x_0,x_1,...,x_p]\) the design matrix holding \(p+1\) known predictors,
\(x_i = (x_{1,i},...,x_{n,i})'\) the \(n\) known values of predictor \(i\), \(x_0\) usually containing only ones,
\(\beta=(\beta_0,...,\beta_p)'\) the vector with \(p+1\) unobserved parameters,
\(e = (e_1, ..., e_n)'\) the vector with \(n\) unobserved residuals.
Notation: \((1,2,3)\) is a row vector, \((1,2,3)'\) is a column vector, \('\) refers to transpose (swap rows and columns).
In the linear we assume \(e\) to be a random variable, and \(X\) and \(\beta\) to be fixed (“fixed effects”). This means that \(y\) is also random, and since \(E(e) = 0\) we have \(E(y) = X\beta\),
As an example, we could regress (yearly) average temperature on geographic latitude and elevation, so that
\(y\) contains the average temperatures at sensor locations
\(x_0=1\),
\(x_1\) contains the latitudes corresponding to the rows in \(y\), and
\(x_2\) contains the elevations corresponding to the rows in \(y\)
If we had yearly mean measurements over a number of years, we could add
\(x_3\) the year corresponding to each row in \(y\)
NoteExercise
Find such a dataset and carry out this regression, first without time for a given year, then with year
2.1 OLS, Normal equations
We need a way to estimate \(\beta\), the vector with unknown parameters in the linear model. When we assume that
the residuals \(e\) are independent and
the residuals \(e\) have equal (constant) variance \(\sigma^2\)
then life becomes relatively simple. Ordinary least squares (OLS) estimates of \(\beta\) are obtained by setting the derivative of the scalar \(R = e'e\) (the sum of squared residuals) to \(\beta\) to zero:
(where we use \((AB)'=B'A'\) and the fact that the transpose of a scalar is itself)
Derivative to \(\beta\):
\[\frac{d R}{d \beta} = -2X'y + 2 X'X \beta\]
Set to zero:
\[-2X'y + 2 X'X \beta=0\]
\[ X'X \beta = X'y\]
These are called the normal equations, and solving them involves solving \(p+1\) equations with \(p+1\) unknowns. We can write the solution explicitly:
\[ \hat{\beta} = (X'X)^{-1} X'y\]
with \(()^{-1}\) indicating matrix inverse and \(\hat{\beta}\) indicating that we estimate\(\beta\) from \(y\) and \(X\).
NoteExercise
Not every matrix can be inverted; which property does \(X\) need to have for \(X'X\) to be invertible?
We can estimate the residual variance \(\sigma^2\) by
\[s^2 = \frac{1}{n-(p+1)} \hat{e}'\hat{e}\]
with \(\hat{e}=y-X\hat{\beta}\).
2.2 Inference
If, in addition, we assume that the errors are normally distributed (\(\mathcal{N}(0, \sigma^2 I)\) with \(I\) the identity matrix), we can do tests on \(\beta\) and create prediction intervals for unobserved values \(y\).
Standard errors for individual \(\beta_i\) are then estimated by \(\sqrt{s^2
(X'X)^{-1}_{i,i}}\); confidence intervals can be constructed using a t distribution with \(n-p-1\) degrees of freedom.
NoteExercise
Write down the equation for a 95% confidence interval for \(\beta_1\).
Prediction of \(y\) at a particular vector with predictor values \(\tilde{x}\) can take two forms. For both, the predicted value is \(\hat{y} = \tilde{x} \hat{\beta}\), but they differ in variance:
the mean value of of \(\hat{y}\) has variance \(Var(\tilde{x}\hat{\beta}) = \tilde{x}Var(\hat{\beta})\tilde{x}' = s^2 \tilde{x}(X'X)^{-1}\tilde{x}'\)
an individual value of of \(\hat{y}\) has variance \(s^2 (1 + \tilde{x}(X'X)^{-1}\tilde{x}')\), adding \(s^2\) the variance of measurements around their mean (the regression surface)
2.3 Weighted least squares
If the observations \(y\) are independent but have a non-constant variance, one can carry out weighted least squares (WLS) to estimate the parameters, giving more accurate observations (with smaller variance) more weight then less accurate ones.
If \(Var(y) = \sigma^2 W^2\) with \(W\) a diagonal matrix with weight \(w_{i,i}\) the weight of observation \(i\). It then follows that \(Var(W^{-1} y) = W^{-1} Var(y) W^{-1} = W^{-1} (\sigma^2 W^2) W^{-1}
= \sigma^2 I\) with \(I\) the identity matrix. So for WLS estimation, one can “unweight” the observations and predictors, and use \[W^{-1} y = W^{-1} (X \beta + e)\]
Since \(Var(W^{-1}e)=\sigma^2I\), after unweighting \(y\) and \(X\) we can use OLS to estimate \(\beta\), which is called WLS.
2.4 Generalized least squares
Similar to WLS, if the observations \(y\) are dependent and have a known covariance \(\Sigma\), one can carry out generalized least squares (GLS) to estimate the parameters, giving more weight to observations that carry more independent information, compared to those that are more strongly correlated with others.
If \(Var(y) = \Sigma\) with \(\Sigma_{i,j}=Cov(y_i,y_j)\), and Choleski decomposition \(\Sigma=CC'\), then we can decorrelate \(y\) by pre-multiplying it with the inverse of \(C\):
\(Var(C^{-1} y) = C^{-1} Var(y) (C^{-1})' = C^{-1} \Sigma (C^{-1})' = C^{-1}C C' (C^{-1})' = I\) with \(I\) the identity matrix (so \(y\) now also has variance 1). So for GLS estimation, one can “weigh” the observations and predictors, and use \[C^{-1} y = C^{-1} X \beta + C^{-1} e\]
Since \(Var(\Sigma^{-1/2}e)=I\), after unweighting \(y\) and \(X\) we can use OLS to estimate \(\beta\), which is called generalized least squares (GLS). Substituting into the normal equations, one gets the generalized least squares estimate for \(\beta\), \[\hat{\beta}_{gls} = (X' \Sigma^{-1}X)^{-1}X' \Sigma^{-1}y,\] which has estimation error (co)variances \(Cov(\hat{\beta})=(X'\Sigma^{-1}X)^{-1}\).
2.5 Partial correlations, causality
Correlation does not imply causality. Regression and partial correlation however can give some insights in possible causal networks. If X1 influences X2, and X2 influences Y, Y is correlated with both X1 and X2, but multiple linear regression will only indicate X2 as an influence:
X1 X2 Y
X1 1.0000000 0.8792068 0.8766155
X2 0.8792068 1.0000000 0.9960418
Y 0.8766155 0.9960418 1.0000000
summary(lm(Y~X1))
Call:
lm(formula = Y ~ X1)
Residuals:
Min 1Q Median 3Q Max
-1.3596 -0.3312 -0.0212 0.3440 1.4886
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.006774 0.023009 0.294 0.769
X1 0.969149 0.023839 40.654 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5141 on 498 degrees of freedom
Multiple R-squared: 0.7685, Adjusted R-squared: 0.768
F-statistic: 1653 on 1 and 498 DF, p-value: < 2.2e-16
summary(lm(Y~X2))
Call:
lm(formula = Y ~ X2)
Residuals:
Min 1Q Median 3Q Max
-0.24705 -0.06609 -0.00489 0.06241 0.33890
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.013452 0.004248 -3.167 0.00163 **
X2 0.500865 0.002003 250.069 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.09497 on 498 degrees of freedom
Multiple R-squared: 0.9921, Adjusted R-squared: 0.9921
F-statistic: 6.253e+04 on 1 and 498 DF, p-value: < 2.2e-16
summary(lm(Y~X1+X2))
Call:
lm(formula = Y ~ X1 + X2)
Residuals:
Min 1Q Median 3Q Max
-0.24332 -0.06716 -0.00421 0.06197 0.33740
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.013346 0.004257 -3.135 0.00182 **
X1 0.004329 0.009250 0.468 0.64002
X2 0.499134 0.004207 118.638 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.09504 on 497 degrees of freedom
Multiple R-squared: 0.9921, Adjusted R-squared: 0.9921
F-statistic: 3.122e+04 on 2 and 497 DF, p-value: < 2.2e-16
The direction of the causal graph however is not revealed:
summary(lm(X1 ~ X2+Y))
Call:
lm(formula = X1 ~ X2 + Y)
Residuals:
Min 1Q Median 3Q Max
-1.24182 -0.34292 0.02147 0.32267 1.16818
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02317 0.02082 -1.113 0.26616
X2 0.34894 0.10934 3.191 0.00151 **
Y 0.10175 0.21743 0.468 0.64002
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4608 on 497 degrees of freedom
Multiple R-squared: 0.7731, Adjusted R-squared: 0.7722
F-statistic: 846.7 on 2 and 497 DF, p-value: < 2.2e-16
Partial correlations consider the same thing, but do not have the asymetry of dependent - independent variables that regression has; in addition the coefficients are standardized to range between -1 and 1 (where regression coefficients depend on the units of X and Y):
library(ppcor)
Loading required package: MASS
pcor(cbind(X1,X2,Y))
$estimate
X1 X2 Y
X1 1.00000000 0.1417106 0.02098622
X2 0.14171057 1.0000000 0.98279878
Y 0.02098622 0.9827988 1.00000000
$p.value
X1 X2 Y
X1 0.000000000 0.001505269 0.6400189
X2 0.001505269 0.000000000 0.0000000
Y 0.640018857 0.000000000 0.0000000
$statistic
X1 X2 Y
X1 0.0000000 3.191432 0.4679593
X2 3.1914318 0.000000 118.6380124
Y 0.4679593 118.638012 0.0000000
$n
[1] 500
$gp
[1] 1
$method
[1] "pearson"