Preliminary Readings

• Using R for Introductory Econometrics, Florian Heiss, pp:81-158

• Introduction to Econometrics with R, Christoph Hanck, Martin Arnold, Alexander Gerber, and Martin Schmelzer, Chapters:4-7

• An Introduction to Statistical Learning with Applications in R, Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani, Chapter 3

Learning Objectives

• How to evaluate Linear Model (Regression)

• [Some] Major Assumptions

  • Linear Relationship
    
  • Constant Variance of Error Term, \(V(\varepsilon_i)=s^2\)
    
  • Errors are uncorrelated \(Cov(\varepsilon_i,\varepsilon_j)=0 \; for \; i \neq j\)
    
  • …

• Multicollinearity

• Outliers

Remember:

• Inputs: independent variable, covariates, predictors, features, regressors.

• Output: dependent variable, variates, labels, regressand.

Linear Models (Linear Regression)

Steps involve Regression analysis:

• Model Formulation (Model Selection will be discussed later)

• Model Estimation

• Model Evaluation

• Model Use

Some Questions: Linear Models

• Is there a [linear] relationship between dependent and independent variable

• How strong is this relationships

• How accurately we can predict the dependent variable given the values of independent variable

• …

Are the assumptions valid given data

Linear Relationship: Simple Case

• Which one of the model is a better fit to the data

Linear Relationship: Simple Case

Linear Relationship: Simple Case

Coefficients

Residuals

• Before comparing the models, let’s see the residuals

Some Major Assumptions

• Linear Relationship

• Error term \(\varepsilon\) has zero mean \(E[\varepsilon]=0\)

• Error term \(\varepsilon\) has constant variance \(V(\varepsilon_i)=s^2\)

• Errors are uncorrelated \(Cov(\varepsilon_i,\varepsilon_j)=0 \; for \; i \neq j\)

• Errors are Normally Distributed \(\varepsilon \sim N(0,s^2)\) (or we have enough observations)

\(\varepsilon_i=y_i-\hat y_i\)

\(\hat y_i=X \hat \beta\)

\(E[\varepsilon]=\bar \varepsilon=\frac{1}{n}\sum_{i}\varepsilon_i=0\)

\(V[\varepsilon]=s^2=\frac{1}{n-p-1}\sum_{i}\varepsilon_i^2\)

An Example: Advertising Data

An Example: Advertising Data

Call:
lm(formula = sales ~ ., data = advertising)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8277 -0.8908  0.2418  1.1893  2.8292 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.938889   0.311908   9.422   <2e-16 ***
TV           0.045765   0.001395  32.809   <2e-16 ***
radio        0.188530   0.008611  21.893   <2e-16 ***
newspaper   -0.001037   0.005871  -0.177     0.86    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.686 on 196 degrees of freedom
Multiple R-squared:  0.8972,    Adjusted R-squared:  0.8956 
F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

Residuals

• Linear Relationship: Visual Clues (at least for now, Ramsey Regression Equation Specification Error Test, RESET, may be used)

library(lmtest)
resettest(lm1, power = 2, type = "regressor")

    RESET test

data:  lm1
RESET = 15.902, df1 = 3, df2 = 193, p-value = 2.795e-09

• Error term \(\varepsilon\) has zero mean \(E[\varepsilon]=0\)

\(E[\varepsilon]=\bar \varepsilon=\frac{1}{n}\sum_{i}\varepsilon_i=0\)

Residuals

• Error term \(\varepsilon\) has constant variance \(V(\varepsilon_i)=s^2\)

Click here to read more

• Homoscedasticity is necessary to calculate accurate standard errors for parameter estimates

To diagnose, visually check plots of residuals against fitted values or predictors for constant variance (at least for now, the Breusch-Pagan test against heteroscedaticity can be used)

library(car)
# error variance changes with the level of the response
ncvTest(lm1)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 5.355982, Df = 1, p = 0.020651

Residuals

• Plot the residuals variance against the fitted values and against each of the predictors

• Errors are independent after controlling predictors \(Cov(\varepsilon_i,\varepsilon_j|X)=0 \; for \; i \neq j\)

Omitting variables may result in biased estimation of parameters and standard errors

Statistical tests as well as visual tests of independence are usually insufficient

Data Collection, sample selection, repeated measures, spatial/temporal observations, etc. are important questions

Residuals

• Errors are Normally Distributed \(\varepsilon \sim N(0,s^2)\)
  Both visual (Q-Q plot) and statistical tests (Shapiro test for example for relatively small samples, \(n \leq 5000\)) can be used to diagnose whether error terms are normally distributed

Get Standardized residuals (rstandard()) and perform shapiro test (shapiro.test())

Standardized residuals, \(z_i=\frac{\varepsilon_i}{\sqrt{MSE}}, \; z_i \sim N(0,1)\)

shapiro.test(rstandard(lm1))

    Shapiro-Wilk normality test

data:  rstandard(lm1)
W = 0.91662, p-value = 3.312e-09

Residuals

Q-Q Plot

Skewness of standardized residuals of the advertising model is -1.332

Simulated Random Values and Q-Q Plot

Q-Q plot for simulated random variable that follows normal distribution

The Shapiro-Wilk test

shapiro.test(x)

    Shapiro-Wilk normality test

data:  x
W = 0.9995, p-value = 0.9086

Residuals

Skewed distributions (Click here to read more)

Positively Skewed

Residuals

Negatively Skewed

Residuals

Fat Tails (example, t-distribution with 4 degrees of freedom)

Residuals

With 2 peaks

Residuals

With 3 peaks

What to do

Skew

Fit GLM model (not a part of this course)

Transform the dependent variable

Multiple peaks

Check for omitted categorical variable(s)

Fat-tailed

Transform the dependent variable

No Multicollinearity

• The information contained in one variables should not be contained by other variable(s) significantly

• Multicollinearity results in increased standard errors

• Variance Inflation Factor (VIF) can be used to diagnose (VIF value should be less than 4 or 5 most of the models)

• These variables can be combined (factors), some of them may be dropped

let, \(y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \varepsilon_i\)

If the correlation between \(x_{1i} \; ve \; x_{2i}\) are very high, the estimate of \(\beta_1\) is consistent and unbiased but it has a large variance

\(\sigma^2_{\hat\beta_1} = \frac{1}{n} \left( \frac{1}{1-\rho^2_{X_1,X_2}} \right) \frac{\sigma^2_\varepsilon}{\sigma^2_{X_1}}\)

To see the effect of high correlation, let’s have two simulation examples where the correlations between these variables are 0.25 and 0.85

Multicollinearity

\[x_i = (x_{1i}, x_{2i}) \overset{i.i.d.}{\sim} \mathcal{N} \left[\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 10 & 2.5 \\ 2.5 & 10 \end{pmatrix} \right]\]

\[\rho_{X_1,X_2} = \frac{Cov(X_1,X_2)}{\sqrt{Var(X_1)}\sqrt{Var{(X_2)}}} = \frac{2.5}{10} = 0.25\]
Here is the simulated Covariance Matrix and the joint density estimation of \(\hat{\beta_i}\)

hat_beta_1 hat_beta_2 
0.05674375 0.05712459 

Multicollinearity

\[x_i = (x_{1i}, x_{2i}) \overset{i.i.d.}{\sim} \mathcal{N} \left[\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 10 & 8.5 \\ 8.5 & 10 \end{pmatrix} \right]\]

\[\rho_{X_1,X_2} = \frac{Cov(X_1,X_2)}{\sqrt{Var(X_1)}\sqrt{Var{(X_2)}}} = \frac{8.5}{10} = 0.85\]

hat_beta_1 hat_beta_2 
 0.1904949  0.1909056 

Variance Inflation Factor

\(VIF_j=\frac{1}{1-R_{j}^{2}}\)

where \(R_{j}^{2}\), is obtained from a model where one predictor, \(j^{th}\), is regressed on all of the other predictors

\(VIF_j=1\) implies perfectly uncorrelated predictors where \(VIF_j \to \infty\) for perfectly correlated predictors

Advertising Data Example:

No Outliers

• Outliers have disproportionately large influence on the predicted values and/or model parameter estimates

• They may be classified into three types:

  • Extreme values
    
  • Leverage
    
  • Influence,
    • Influence on predicted values
      
    • Influence on parameter estimates

No Outliers

Extreme Values

They can be identified as points with high residuals

Optional: Studentized Residuals

\(MSE\) is the approximation of the variance of \(\varepsilon_i\). One can improve this by scaling,

\(V(\varepsilon_i)=\sigma^2(1-h_i)\)

where, \(h_i\) is the \(i^{th}\) element on the diagonal of the hat matrix, \(H=X(X^TX)^{-1}X^T\)

Studentized residual, \(r_i=\frac{\varepsilon_i}{\sqrt{MSE(1-h_i)}}\)
\(r_i \sim t-distribution\) with \(n-p-1\) degrees of freedom. One can consider values outside the range \([-2, 2]\) as potential outliers

rstudent() function helps us to obtain studentized residuals

In our advertising model, the potential outliers (extreme values of residuals) are

No Outliers

Leverage

• It is a measure of distance between individual values of a predictor and other values of the predictor

• Points with high leverage have the potential to influence our model estimates

• Values range from 0 to 1 where \(Leverage_i > \frac {2k} {n}\) can be considered large value* (k is the number of parameters, n is the sample size)

*: Belsley, D. A., Kuh, E., and Welsch, R. E. (1980). Regression Diagnostics: Identifying influential data and sources of collinearity. Wiley. https://doi.org/10.1002/0471725153

hatvalues() can be used for leverage

Advertising data model then have a cut-off value approximately \(Leverage_i > \frac {2k} {n}=0.04\)

Here is a table showing these observations

No Outliers

Influence

• Leverage and residuals are abstract outlier metrics

• Substantive impact of an observation in model, Influence, are more important

• Influence is a measure of how much an observation affects our model estimates

Measures of Influences that are used widely:

• Cook’s Distance

• DFFITS: standardized measure of how much the prediction for a given observation would change if it were deleted from the model
  \(DFFITS_i = r_i \times \sqrt{ \frac { h_i } { 1 - h_i } }\)

where \(r_i\) is studentized residual, \(h_i\) is the leverage

Cut-off value \(|DFFITS_i| > 2 \times \sqrt{ \frac {k} {n} }\) where \(k\) is the number of parameters

dffits() function can be used

Here is a table showing ten of these observations

• DFBETAS: standardized differences between regression coefficients in a model with a given observation, and a model without that observation

• DFBETAS are standardized by the standard error of the coefficient. A model, then, has \(n \times k \times DFBETAS\), one for each combination of observations and predictors

Cut-off value \(|DFBETAS_i| > \frac {2} {\sqrt{n}}\)

Here is a table showing some of these observations