Multicollinearity

Introduction

The text is motivated by Marcin Rutecki and Manimala, both of the Kaggle project.

Multicollinearity - definition

Multicollinearity means the situation, where the explanatory variables are highly pairwise correlated.

Consequences of Multicollinearity (according to Rutecki)

Multicollinearity could result in significant problems during model fitting. It can reduce the overall performance of regression and classification models:

• does not increase bias, but it can increase variance (overfitting);
• make the estimates very sensitive to minor changes in the model;
• doesn’t affect the predictive power, but individual predictor variable’s impact on the response variable could be calculated wrongly.

The result is that the coefficient estimates are unstable and difficult to interpret. Multicollinearity saps the statistical power of the analysis, can cause the coefficients to switch signs, and makes it more difficult to specify the correct model.

Remedies to Multicollinearity problem (according to Rutecki)

1. Increase the sample size. Enlarging the sample will introduce more variation in the data series, which reduces the effect of sampling error and helps increase precision when estimating various properties of the data. Increased sample sizes can reduce either the presence or the impact of multicollinearity, or both.

2. Make some transformation of the data (for example, replace level data by their differences)

3. Remove some of the highly correlated features.

• Manual Method - Variance Inflation Factor (VIF)
• Automatic Method - Recursive Feature Elimination (RFE) (I’ll skip it in this study)
• Feature Elmination using PCA Decomposition (I’ll skip it in this study)
• Replace highly correlated regressors with a linear combination of them. (I’ll skip it in this study)

3. Use regularization methods such as RIDGE and LASSO or Bayesian regression. (I’ll skip it in this study)

Boston House Prices

Each record in the database describes a Boston suburb or town. The data was drawn from the Boston Standard Metropolitan Statistical Area (SMSA) in 1970. The attributes are defined as follows (taken from the UCI Machine Learning Repository1): CRIM: per capita crime rate by town

• ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
• INDUS: proportion of non-retail business acres per town
• CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
• NOX: nitric oxides concentration (parts per 10 million)
• RM: average number of rooms per dwelling
• AGE: proportion of owner-occupied units built prior to 1940
• DIS: weighted distances to five Boston employment centers
• RAD: index of accessibility to radial highways
• TAX: full-value property-tax rate per $10,000
• PTRATIO: pupil-teacher ratio by town 12. B: 1000(Bk−0.63)2 where Bk is the proportion of blacks by town 13. LSTAT: % lower status of the population
• MEDV: Median value of owner-occupied homes in $1000s

We can see that the input attributes have a mixture of units.

Computations

rm(list=ls())
library(lmtest)
library(car)
library(glmnet)

column_names <- c("CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT", "MEDV")
df <- read.table("adot/housing.data", header = FALSE, sep = "", col.names = column_names)

### Plotting the correlation between various columns
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(df, lower.panel = panel.smooth, upper.panel = panel.cor,
      gap=0, row1attop=FALSE)

Overparametrized model

model_all <- lm(MEDV ~ ., data=df)  # with all the independent variables in the dataframe

summary(model_all)

Call:
lm(formula = MEDV ~ ., data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.595  -2.730  -0.518   1.777  26.199 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
CRIM        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
ZN           4.642e-02  1.373e-02   3.382 0.000778 ***
INDUS        2.056e-02  6.150e-02   0.334 0.738288    
CHAS         2.687e+00  8.616e-01   3.118 0.001925 ** 
NOX         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
RM           3.810e+00  4.179e-01   9.116  < 2e-16 ***
AGE          6.922e-04  1.321e-02   0.052 0.958229    
DIS         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
RAD          3.060e-01  6.635e-02   4.613 5.07e-06 ***
TAX         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
PTRATIO     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
B            9.312e-03  2.686e-03   3.467 0.000573 ***
LSTAT       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406,    Adjusted R-squared:  0.7338 
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

Multicollinearity testing

VIF

\[VIF_i = \frac{1}{1-R^2_i}\] niekedy sa hovorí, že tento ukazovateľ pre danú premennú nesmie byť väčší, ako 5, niekedy sa hovorí o hranici 10.

vif(model_all)

    CRIM       ZN    INDUS     CHAS      NOX       RM      AGE      DIS      RAD 
1.792192 2.298758 3.991596 1.073995 4.393720 1.933744 3.100826 3.955945 7.484496 
     TAX  PTRATIO        B    LSTAT 
9.008554 1.799084 1.348521 2.941491 

Vizualizácia VIF

vif_values <- vif(model_all)           #create vector of VIF values

barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue") #create horizontal bar chart to display each VIF value

abline(v = 5, lwd = 3, lty = 2)    #add vertical line at 5 as after 5 there is severe correlation

NA
NA

df_9_10 <- df[,-c(9,10)]
model_9_10 <- lm(MEDV ~., data=df_9_10)
summary(model_9_10)

Call:
lm(formula = MEDV ~ ., data = df_9_10)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.7734  -2.7928  -0.6504   1.5477  27.9833 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  29.210710   4.893167   5.970 4.54e-09 ***
CRIM         -0.061435   0.030434  -2.019 0.044066 *  
ZN            0.041516   0.013551   3.064 0.002306 ** 
INDUS        -0.045831   0.055961  -0.819 0.413190    
CHAS          3.085016   0.871868   3.538 0.000441 ***
NOX         -14.594752   3.670109  -3.977 8.03e-05 ***
RM            4.134272   0.419341   9.859  < 2e-16 ***
AGE          -0.004295   0.013415  -0.320 0.748983    
DIS          -1.488851   0.203337  -7.322 9.98e-13 ***
PTRATIO      -0.811272   0.121647  -6.669 6.91e-11 ***
B             0.008205   0.002706   3.032 0.002558 ** 
LSTAT        -0.515812   0.051668  -9.983  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.839 on 494 degrees of freedom
Multiple R-squared:  0.7293,    Adjusted R-squared:  0.7232 
F-statistic:   121 on 11 and 494 DF,  p-value: < 2.2e-16

vif(model_9_10)

    CRIM       ZN    INDUS     CHAS      NOX       RM      AGE      DIS  PTRATIO 
1.478206 2.154483 3.179166 1.057805 3.901348 1.872532 3.075755 3.954443 1.496077 
       B    LSTAT 
1.316559 2.936487 

Is heteroskedasticity present in the new model?

df_9_10 <- df[,-c(9,10)]
model_9_10 <- lm(MEDV ~., data=df_9_10)
bptest(model_9_10)

    studentized Breusch-Pagan test

data:  model_9_10
BP = 50.806, df = 11, p-value = 4.481e-07

Yes, it is…