Multikolinearnost i drugi problemi kod linearne regresije

Kažemo da u linearnom modelu postoji multikolinearnost ako postoje dve prediktora koja su srednje do jako korelisana. Kada imao dva korelisana prediktora u modelu, to loše utiče na regresiju. Postoje dva tipa multikolinearnosti

Strukturna multikolinearnost - na primer imamo prediktore \(x\) i \(x^2\)
Multikolinearnost iz podataka - ako se samo oslanjamo na podatke koje imamo i ne možemo da kontrolišemo prediktore tako da osiguramo da nema multikolinearnosti.

Primer 1

Posmatrajmo bazu bloodpress.txt u kojoj su dati podaci o 20 osoba koje imaju visok krvni pritisak.

\(y\) - BP - krvni pritisak (u mmHg)
\(x_1\) - Age - starost
\(x_2\) - Weight - težina u kg
\(x_3\) - BSA - body surface area (u \(m^2\))
\(x_4\) - Dur - period trajanja hipertenzije, u godinama
\(x_5\) - Pulse - puls
\(x_6\) - Stress - indeks stresa

Sa sledećeg grafika možemo videti korelisanost prediktora. Osim toga možemo izračunati korelacionu matricu svih kolona.

bloodpress <- read.table("Multikolinearnost/bloodpress.txt", header=T)
attach(bloodpress)

pairs(bloodpress[,c(2:5)])

pairs(bloodpress[,c(2,6:8)])

round(cor(bloodpress[,c(2:8)]),3)

##           BP   Age Weight   BSA   Dur Pulse Stress
## BP     1.000 0.659  0.950 0.866 0.293 0.721  0.164
## Age    0.659 1.000  0.407 0.378 0.344 0.619  0.368
## Weight 0.950 0.407  1.000 0.875 0.201 0.659  0.034
## BSA    0.866 0.378  0.875 1.000 0.131 0.465  0.018
## Dur    0.293 0.344  0.201 0.131 1.000 0.402  0.312
## Pulse  0.721 0.619  0.659 0.465 0.402 1.000  0.506
## Stress 0.164 0.368  0.034 0.018 0.312 0.506  1.000

detach(bloodpress)

Vidimo da je BP delimično korelisan sa Weight i BSA, a nekorelisan sa Stress. Vidimo i da su Weight i BSA jako korelisani, a Stress i BSA nekorelisani.

Nekorelisani prediktori

Veštački primer

Vrednosti prediktora su namešteni tako da korelacija između njih bude nula. Napravićemo više linearnih modela, da vidimo kako nekorelisani prediktori utiču na model.

uncorrpreds <- read.table("Multikolinearnost/uncorrpreds.txt",header=T)
attach(uncorrpreds)

pairs(uncorrpreds)

cor(x1,x2)

## [1] 0

model.1 <- lm(y ~ x1)
summary(model.1)

## 
## Call:
## lm(formula = y ~ x1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.500 -2.500  0.500  1.938  3.750 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   52.750      3.346  15.764 4.13e-06 ***
## x1            -1.625      1.058  -1.536    0.176    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.993 on 6 degrees of freedom
## Multiple R-squared:  0.2821, Adjusted R-squared:  0.1625 
## F-statistic: 2.358 on 1 and 6 DF,  p-value: 0.1755

anova(model.1)

## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value Pr(>F)
## x1         1 21.125 21.1250  2.3581 0.1755
## Residuals  6 53.750  8.9583

Primetimo da je regresiona suma kvadrata \(SSR(x_1)=8\), da je ocenjeni koeficijent \(b_1=-1\), a njegova standardna greška \(se(b_1)=1.47\). Pravimo drugi model gde je prediktor \(x_2\)

model.2 <- lm(y ~ x2)
summary(model.2)

## 
## Call:
## lm(formula = y ~ x2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.500 -1.688  0.125  1.000  3.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  62.1250     4.7888  12.973 1.29e-05 ***
## x2           -2.3750     0.7873  -3.017   0.0235 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.227 on 6 degrees of freedom
## Multiple R-squared:  0.6027, Adjusted R-squared:  0.5364 
## F-statistic: 9.101 on 1 and 6 DF,  p-value: 0.02349

anova(model.2)

## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## x2         1 45.125  45.125  9.1008 0.02349 *
## Residuals  6 29.750   4.958                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dobijamo da je \(SSR(x_2)=24.5\), da je ocenjeni koeficijent \(b_2=-1.75\), a njegova standardna greška \(se(b_2)=1.35\). Pravimo i treći model koji sadrži oba prediktora.

model.12 <- lm(y ~ x1 + x2)
summary(model.12)

## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##      1      2      3      4      5      6      7      8 
##  0.125 -0.875  1.875 -1.125  1.375 -0.625  0.125 -0.875 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  67.0000     3.1494  21.274 4.25e-06 ***
## x1           -1.6250     0.4644  -3.499  0.01729 *  
## x2           -2.3750     0.4644  -5.115  0.00372 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.313 on 5 degrees of freedom
## Multiple R-squared:  0.8848, Adjusted R-squared:  0.8387 
## F-statistic:  19.2 on 2 and 5 DF,  p-value: 0.004504

anova(model.12)

## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## x1         1 21.125  21.125  12.246 0.017294 * 
## x2         1 45.125  45.125  26.159 0.003724 **
## Residuals  5  8.625   1.725                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ocenjeni koeficijenti su \(b_1=-1\), \(b_2=-1.75\), sa standardnim greškama \(se(b_1)=1.47\) i \(se(b_2)=1.35\). Sekvencijalna regresiona suma kvadrata je \(SSR(x_2| x_1)=24.5\). Da bismo dobili i drugu regresionu sumu kvadrata pravimo i model sa prediktorima \(x_2\) i \(x_1\) u tom poretku.

model.21 <- lm(y ~ x2 + x1)
summary(model.21)

## 
## Call:
## lm(formula = y ~ x2 + x1)
## 
## Residuals:
##      1      2      3      4      5      6      7      8 
##  0.125 -0.875  1.875 -1.125  1.375 -0.625  0.125 -0.875 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  67.0000     3.1494  21.274 4.25e-06 ***
## x2           -2.3750     0.4644  -5.115  0.00372 ** 
## x1           -1.6250     0.4644  -3.499  0.01729 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.313 on 5 degrees of freedom
## Multiple R-squared:  0.8848, Adjusted R-squared:  0.8387 
## F-statistic:  19.2 on 2 and 5 DF,  p-value: 0.004504

anova(model.21)

## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## x2         1 45.125  45.125  26.159 0.003724 **
## x1         1 21.125  21.125  12.246 0.017294 * 
## Residuals  5  8.625   1.725                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dobijamo \(SSR(x_1|x_2)=8\). Sada ćemo da uporedimo ova tri modela. Vidimo da

Ocenjene vrednosti koeficijenata su iste nezavisno da li su dobijene iz modela sa jednim ili oba prediktora. Ovo znači da uticaj koji jedan prediktor ima na model ne zavisi od drugog prediktora.
Standardne greške koeficijenata se manjaju ali jako malo u zavisnosti od toga da li je model sa jednim ili oba prediktora.
Sekvencijalne regresione sume kvadrata, u slučaju modela sa dva prediktora su iste kao obične regresione sume kvadrata modela sa jednim prediktorom! Važi \(SSR(x_1)=SSR(x_1|x_2)\) i \(SSR(x_2)=SSR(x_2|x_1)\). Ovo znači da doprinos prediktora \(x_1\) za smanjivanje varijabilnosti u zavisnoj promenljivoj ne zavisi od prediktora \(x_2\) i obrnuto.

library(car)

## Warning: package 'car' was built under R version 3.4.4

scatter3d(y ~ x1 + x2)

## Loading required namespace: rgl

detach(uncorrpreds)

Realni podaci

Vratimo se bazi bloodpress.txt. Hoćemo da ispitamo veza između zavisne promenljive \(y=BP\) i prediktora \(x_3=BSA\) i \(x_6=Stress\)

attach(bloodpress)

pairs(bloodpress[,c(2,5,8)])

round(cor(bloodpress[,c(2,5,8)]),3)

##           BP   BSA Stress
## BP     1.000 0.866  0.164
## BSA    0.866 1.000  0.018
## Stress 0.164 0.018  1.000

Sa grafika i iz korelacione matrice vidimo da su BP i BSA jako linearno zavisni, BP i Stress slabo zavisni, a prediktori BSA i Stress gotovo nekorelisani. Napravićemo više modela i uporediti ih, kao u prethodnom primeru. Prvi model je sa jednim prediktorom \(x_6=Stress\).

model.1 <- lm(BP ~ Stress)
summary(model.1)

## 
## Call:
## lm(formula = BP ~ Stress)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.6394 -3.3014  0.0722  2.2181  9.9287 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 112.71997    2.19345  51.389   <2e-16 ***
## Stress        0.02399    0.03404   0.705     0.49    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.502 on 18 degrees of freedom
## Multiple R-squared:  0.02686,    Adjusted R-squared:  -0.0272 
## F-statistic: 0.4969 on 1 and 18 DF,  p-value: 0.4899

anova(model.1)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stress     1  15.04  15.044  0.4969 0.4899
## Residuals 18 544.96  30.275

Ocenjeni koeficijent je \(b_6=0.024\) sa standardnom greškom \(se(b_6)=0.034\) a regresiona suma kvadrata modela je \(SSR(x_6)=15.04\). Drugi model je takođe sa jednim prediktorom \(x_3=BSA\)

model.2 <- lm(BP ~ BSA)
summary(model.2)

## 
## Call:
## lm(formula = BP ~ BSA)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.314 -1.963 -0.197  1.934  4.831 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   45.183      9.392   4.811  0.00014 ***
## BSA           34.443      4.690   7.343 8.11e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.79 on 18 degrees of freedom
## Multiple R-squared:  0.7497, Adjusted R-squared:  0.7358 
## F-statistic: 53.93 on 1 and 18 DF,  p-value: 8.114e-07

anova(model.2)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## BSA        1 419.86  419.86  53.927 8.114e-07 ***
## Residuals 18 140.14    7.79                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ocenjeni koeficijent je \(b_3=34.44\) sa standardnom greškom \(se(b_3)=4.69\) a regresiona suma kvadrata modela je \(SSR(x_3)=419.858\). Treći model, sa oba prediktora, u redosledu \(b_6=Stress, b_3=BSA\) je

model.12 <- lm(BP ~ Stress + BSA)
summary(model.12)

## 
## Call:
## lm(formula = BP ~ Stress + BSA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8992 -1.6483 -0.1643  1.7790  3.8524 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 44.24452    9.26104   4.777 0.000175 ***
## Stress       0.02166    0.01697   1.277 0.218924    
## BSA         34.33423    4.61110   7.446 9.56e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.743 on 17 degrees of freedom
## Multiple R-squared:  0.7716, Adjusted R-squared:  0.7448 
## F-statistic: 28.72 on 2 and 17 DF,  p-value: 3.534e-06

anova(model.12)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Stress     1  15.04   15.04  1.9998    0.1754    
## BSA        1 417.07  417.07 55.4430 9.561e-07 ***
## Residuals 17 127.88    7.52                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ocenjeni koeficijenti su \(b_6=0.0217\) je \(b_3=34.33\) sa standardnim greškama \(se(b_6)=0.017\) i \(se(b_3)=4.61\) a sekvencijalna regresiona suma kvadrata \(SSR(x_3|x_6)=417.07\). Pravimo još jedan model da bismo dobili drugu sekvencijalnu sumu \(SSR(x_6|x_3)=12.26\).

model.21 <- lm(BP ~ BSA + Stress)
summary(model.21)

## 
## Call:
## lm(formula = BP ~ BSA + Stress)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8992 -1.6483 -0.1643  1.7790  3.8524 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 44.24452    9.26104   4.777 0.000175 ***
## BSA         34.33423    4.61110   7.446 9.56e-07 ***
## Stress       0.02166    0.01697   1.277 0.218924    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.743 on 17 degrees of freedom
## Multiple R-squared:  0.7716, Adjusted R-squared:  0.7448 
## F-statistic: 28.72 on 2 and 17 DF,  p-value: 3.534e-06

anova(model.21)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## BSA        1 419.86  419.86 55.8132 9.149e-07 ***
## Stress     1  12.26   12.26  1.6296    0.2189    
## Residuals 17 127.88    7.52                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dakle, ako imamo prediktore koji su skoro nekorelisani

Dobijamo slične ocenjene koeficijente nagiba nezavisno od prediktora u modelu.
Obične i odgovarajuće sekvencijalne regresione sume kvadrata nisu iste, ali su slične: \(SSR(x_6|x_3)\approx SSR(x_6)\) i \(SSR(x_3|x_6)\approx SSR(x_3)\).

Znači, efekat koji jedan prediktor ima na zavisnu promenljivu je sličan nezavisno od ostalih prediktora u modelu - u smislu slične vrednosti koeficijenta, a i doprinos prediktora modelu je sličan nezavisno od ostalih prediktora - u smislu količine objašnjene disperzije.

scatter3d(BP ~ Stress + BSA)

detach(bloodpress)

Jako korelisani prediktori

U ovom delu ćemo objasniti zašto nije dobro imati jako korelisane prediktore. U istoj bazi bloodpress.txt posmatramo zavisnu promenljivu \(BP\) i prediktore \(x_2=Weight\) i \(x_3=BSA\).

attach(bloodpress)

pairs(bloodpress[,c(2,5,4)])

cor(bloodpress[,c(2,5,4)])

##               BP       BSA    Weight
## BP     1.0000000 0.8658789 0.9500677
## BSA    0.8658789 1.0000000 0.8753048
## Weight 0.9500677 0.8753048 1.0000000

Vidimo jaku korelisanost prediktora sa zavisnom promenljivom a i među prediktorima što je i očekivano. Opet pravimo tri modela. Prvi ima samo prediktor \(x_2=Weight\)

model.1 <- lm(BP ~ Weight)
summary(model.1)

## 
## Call:
## lm(formula = BP ~ Weight)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6933 -0.9318 -0.4935  0.7703  4.8656 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.20531    8.66333   0.255    0.802    
## Weight       1.20093    0.09297  12.917 1.53e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.74 on 18 degrees of freedom
## Multiple R-squared:  0.9026, Adjusted R-squared:  0.8972 
## F-statistic: 166.9 on 1 and 18 DF,  p-value: 1.528e-10

anova(model.1)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Weight     1 505.47  505.47  166.86 1.528e-10 ***
## Residuals 18  54.53    3.03                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dakle, ocenjeni koeficijent je \(b_2 = 1.2009\), sa standardnom greškom \(se(b_2) = 0.0930\), i regresionom sumom kvadrata \(SSR(x_2) = 505.472\). Sledeći model ima prediktor \(x_3=BSA\).

model.2 <- lm(BP ~ BSA)
summary(model.2)

## 
## Call:
## lm(formula = BP ~ BSA)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.314 -1.963 -0.197  1.934  4.831 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   45.183      9.392   4.811  0.00014 ***
## BSA           34.443      4.690   7.343 8.11e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.79 on 18 degrees of freedom
## Multiple R-squared:  0.7497, Adjusted R-squared:  0.7358 
## F-statistic: 53.93 on 1 and 18 DF,  p-value: 8.114e-07

anova(model.2)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## BSA        1 419.86  419.86  53.927 8.114e-07 ***
## Residuals 18 140.14    7.79                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ocenjeni koeficijent je \(b_3 = 34.44\), sa standardnom greškom \(se(b_4) = 4.69\), i regresionom sumom kvadrata \(SSR(x_3) = 419.858\). Treći model ima oba prediktora, \(x_2=Weight\) i \(x_3=BSA\) u tom poretku, a onda pravimo i još jedan model da bismo dobili drugu sekvencijalnu sumu kvadrata

model.12 <- lm(BP ~ Weight + BSA)
summary(model.12)

## 
## Call:
## lm(formula = BP ~ Weight + BSA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8932 -1.1961 -0.4061  1.0764  4.7524 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.6534     9.3925   0.602    0.555    
## Weight        1.0387     0.1927   5.392 4.87e-05 ***
## BSA           5.8313     6.0627   0.962    0.350    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.744 on 17 degrees of freedom
## Multiple R-squared:  0.9077, Adjusted R-squared:  0.8968 
## F-statistic: 83.54 on 2 and 17 DF,  p-value: 1.607e-09

anova(model.12)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq  F value    Pr(>F)    
## Weight     1 505.47  505.47 166.1648 3.341e-10 ***
## BSA        1   2.81    2.81   0.9251    0.3496    
## Residuals 17  51.71    3.04                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model.21 <- lm(BP ~ BSA + Weight)
summary(model.21)

## 
## Call:
## lm(formula = BP ~ BSA + Weight)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8932 -1.1961 -0.4061  1.0764  4.7524 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.6534     9.3925   0.602    0.555    
## BSA           5.8313     6.0627   0.962    0.350    
## Weight        1.0387     0.1927   5.392 4.87e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.744 on 17 degrees of freedom
## Multiple R-squared:  0.9077, Adjusted R-squared:  0.8968 
## F-statistic: 83.54 on 2 and 17 DF,  p-value: 1.607e-09

anova(model.21)

## Analysis of Variance Table
## 
## Response: BP
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## BSA        1 419.86  419.86 138.021 1.391e-09 ***
## Weight     1  88.43   88.43  29.069 4.871e-05 ***
## Residuals 17  51.71    3.04                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sada su koeficijenti \(b_2=1.039\), \(b_3=5.83\), sa standardnim greškama \(se(b_2)=0.193\), \(se(b_3)=6.06\) i sekvencijalnim regresionim sumama \(SSR(x_3|x_2)=2.814\) i \(SSR(x_2|x_3)=88.43\). Poredimo dobijene rezultate i zaključujemo

Ako su prediktori korelisani, vrednost ocenjenog koeficijenata se puno razlikuje u zavisnosti od ostalih prediktora uključenih u model. U našem primeru se veoma značajno promenila vrednost \(b_3\) kada se dodao prediktor \(x_2\). Dakle, prema modelu sa jednim prediktorom tvrdimo da se za jedinično povećanje \(x_3=BSA\) krvni pritisak (BP) poveća za 34.4. Prema modelu sa oba prediktora tvrdimo da se za jedinično povećanje \(BSA\), i konstantnu vrednost drugog prediktora, \(BP\) poveća za samo 5.83. Razlog za ovo je što, u prvom slučaju, kada se poveća \(BSA\) poveća se i \(Weight\) a oba prediktora utiču na povećanje \(BP\). U drugom slučaju povećavamo \(BSA\) ali \(Weight\) držimo konstantnim pa je manje povećanje u \(BP\). Osim ovog može se desiti i da prediktori koje nismo ni uključili u model, a korelisani su sa onima u modelu utiču na ocenjene koeficijente, tako da oni nemaju smisla. Na primer, farmaceutska kompanija ispituje prodaju njenih proizvoda u raznim gradovima u zavisnosti od veličine populacije i prosečnog prihoda. Očekivano je da što je veća populacija to bude veća prodaja, međutim dobili smo da je koeficijent uz veličinu populacije negativan, što bi značilo da što je veće mesto, to je manja prodaja. U stvari, desilo se da nismo uzeli u obzir konkurentsku firmu, koja je prisutnija u većim mestima. Utvrđeno je da je jako korelisano prisustvo konkurentske firme sa veličinom populacije. Dakle čak i prediktor koji nismo uključili u model, a korelisan je sa nekim prisutnim prediktorom utiče na model i zaključke koje izvodimo na osnovu modela. Zaključak je da, ako dobijamo koeficijente koji nemaju smisla, to može biti posledica korelisanih prediktora, prisutnih u modelu, ili ne!
Osim različitih ocena, menjaju se i njihove standardne greške. Primetimo da su mnogo veće kada imamo model sa dva prediktora. Ovo znači i šire intervale poverenja i manje precizne ocene parametara. Ovo se može objasniti na sledeći način. Na prvoj 3D slici, sa potpuno nekorelisanim prediktorima, podaci su razbacani na četiri strane i čak i da se malo promeni vrednost zavisne promenljive imamo jaku bazu za tu ravan. Na drugoj 3D slici, zbog nekorelisanosti tačke su raspoređene svuda po ravni. Čak i da se malo promene vrednosti zavisne promenljive malo menjaju, opet imamo jaku bazu i parametri provučene ravni se neće puno razlikovati. Ovo upravo znači male standardne greške koeficijenata. Sada pogledajmo treću 3D sliku

scatter3d(BP ~ Weight + BSA)

Zbog korelisanosti, čini se kao da su podaci skoncentrisani oko jedne prave u prostoru. To nije stabilna baza za našu ravan, za manje promene u zavisnoj promenljivoj mogu značajno da se promene koeficijenti ravni, odnosno ocenjeni koeficijenti modela. Ta nestabilnost ogleda se u velikim standardnim greškama

Sledeći efekat multikolinearnosti je u doprinosu smanjivanja sume kvadrata grešaka. U našem primeru dobili smo \(SSR(x_2)=505.472\) i \(SSR(x_2|x_3) = 88.43\). Radi se o tome da dva korelisana prediktora hvataju isti deo varijabilnosti u zavisnoj promenljivoj. Zato, ako već imamo \(x_2=Weight\) u modelu i onda dodamo korelisani prediktor \(x_3=BSA\), njime ne možemo puno popraviti situaciju. Isti je slučaj i da smo gledali drugu sekvencijalnu sumu.
Sledeći problem je da multikolinearnost dovodi do različitih rezultata testiranja hipoteza o značajnosti prediktora u zavisnosti od toga koji su drugi prediktori u modelu. U oba modela sa po jednim prediktorom p-vrednosti testova su jako male - oba prediktora pojedinačno su značajni. Ipak, u zajedničkom modelu, dobijamo da \(BSA\) nije značajan prediktor. Kada već imamo prediktor \(Weight\), \(BSA\) ne objašnjava više puno varijabilnosti.
Osim loših osobina , dobra vest je da multikolinearnost ne utiče na predviđanja. To možemo proveriti u Ru.

predict(model.1, interval="prediction",
        newdata=data.frame(Weight=92))

##       fit     lwr     upr
## 1 112.691 108.938 116.444

predict(model.2, interval="prediction",
        newdata=data.frame(BSA=2))

##        fit      lwr      upr
## 1 114.0689 108.0619 120.0758

predict(model.12, interval="prediction",
        newdata=data.frame(Weight=92, BSA=2))

##        fit      lwr      upr
## 1 112.8794 109.0801 116.6787

Objašnjenje za ovo je da, pošto se tačka (2,92) približno nalazi na pravoj koja opisuje vezu BSA i Weight, predviđena vrednost će biti dobra. Za male promene u zavisnoj promenljivoj, ravan može oscilirati oko prave zavisnosti prediktora, ali ta prava ostaje na ravni, pa i predviđanja za tačke oko te prave.

plot(BSA,Weight)

detach(bloodpress)

Detektovanje Multikolinearnosti korišćenjem VIF

Neki od načina za detektovanje multikolinearnosti su

Promene u vrednostima koeficijenata ili p-vrednostima testova u zavisnosti od parametara u modelu
Korelacije među parovima prediktora

Nije dovoljno ispitivati korelaciju među parovima jer se može desiti da je više prediktora u linearoj vezi, na primer \(X_3=2X_1+5X_2\). Zato nam za ispitivanje multikolinearnosti koristi VIF (Variance Inflation Factor). U prethodnom delu pomenuli smo da su u prisustvu multikolinearnosti standardne greške ocena veće. Za svaki prediktor definisaćemo \(VIF_k\) koji meri to povećanje u disperziji koeficijenta u prisustvu multikolinearnosti. Posmatrajmo model SLR \[y_i=\beta_0+\beta_kx_{ik}+\epsilon_i.\] Na prvom času dobili smo da je disperzija koeficijenta jednaka \[Var(b_k)_{min}=\frac{\sigma^2}{\sum_{i=1}^{n}(x_{ik}-\bar{x}_k)^2}.\] Kada imamo MLR model \[y_i=\beta_0+\beta_1x_{i1}+ \cdots + \beta_kx_{ik}+\cdots +\beta_{p-1}x_{i, p-1} +\epsilon_i\] gde postoje prediktori koji su korelisani sa \(x_k\) onda je disperzija \(b_k\) veća. Može se pokazati da važi sledeći izraz za disperziju ocene \[Var(b_k)=\frac{\sigma^2}{\sum_{i=1}^{n}(x_{ik}-\bar{x}_k)^2}\times \frac{1}{1-R_{k}^{2}},\] gde je \(R_{k}^{2}\) vrednost \(R^2\) za model kod kog je \(x_k\) zavisna promenljiva a prediktori su svi ostali prediktori iz početnog MLR modela. Vidimo da što je veći \(R_{k}^{2}\), odnosno što je jača linearna zavisnost među prediktorima, to je veća disperzija \(Var(b_k)\). Definišemo \(VIF_k\) kao količnik dve disperzije

\[\frac{Var(b_k)}{Var(b_k)_{min}}=\frac{\left( \frac{\sigma^2}{\sum(x_{ik}-\bar{x}_k)^2}\times \frac{1}{1-R_{k}^{2}} \right)}{\left( \frac{\sigma^2}{\sum(x_{ik}-\bar{x}_k)^2} \right)}=\frac{1}{1-R_{k}^{2}}\]

Vrednosti \(VIF_k\) blizu 1 govore da nema multikolinearnost, a velike vrednosti da ima. Vrednosti \(VIF\)-ova veće od 4 govore da bi trebalo dalje istražiti multikolinearnost, a one veće 10 govore da postoji jaka multikolinearnost.

Vraćamo se primeru od ranije, bloodpress.txt, gde smo imali jako korelisane prediktore \(Weight\) i \(BSA\). Prvo pravimo model sa svim prediktorima.

attach(bloodpress)

model.1 <- lm(BP ~ Age + Weight + BSA + Dur + Pulse + Stress)
summary(model.1)

## 
## Call:
## lm(formula = BP ~ Age + Weight + BSA + Dur + Pulse + Stress)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93213 -0.11314  0.03064  0.21834  0.48454 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -12.870476   2.556650  -5.034 0.000229 ***
## Age           0.703259   0.049606  14.177 2.76e-09 ***
## Weight        0.969920   0.063108  15.369 1.02e-09 ***
## BSA           3.776491   1.580151   2.390 0.032694 *  
## Dur           0.068383   0.048441   1.412 0.181534    
## Pulse        -0.084485   0.051609  -1.637 0.125594    
## Stress        0.005572   0.003412   1.633 0.126491    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4072 on 13 degrees of freedom
## Multiple R-squared:  0.9962, Adjusted R-squared:  0.9944 
## F-statistic: 560.6 on 6 and 13 DF,  p-value: 6.395e-15

# library(car)
vif(model.1)

##      Age   Weight      BSA      Dur    Pulse   Stress 
## 1.762807 8.417035 5.328751 1.237309 4.413575 1.834845

Očekivano, dobijamo velike vrednosti \(VIF\) za prediktore koje smo i ispitivali, \(VIF_{Weight}=8.42\), \(VIF_{BSA}=5.33\), \(VIF_{Pulse}=4.41\). Možemo proveriti i računanje vrednosti \(VIF\) na primer za \(Weight\).

model.2 <- lm(Weight ~ Age + BSA + Dur + Pulse + Stress)
summary(model.2)

## 
## Call:
## lm(formula = Weight ~ Age + BSA + Dur + Pulse + Stress)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7697 -1.0120  0.1960  0.6955  2.7035 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.674438   9.464742   2.079  0.05651 .  
## Age         -0.144643   0.206491  -0.700  0.49510    
## BSA         21.421654   3.464586   6.183 2.38e-05 ***
## Dur          0.008696   0.205134   0.042  0.96678    
## Pulse        0.557697   0.159853   3.489  0.00361 ** 
## Stress      -0.022997   0.013079  -1.758  0.10052    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.725 on 14 degrees of freedom
## Multiple R-squared:  0.8812, Adjusted R-squared:  0.8388 
## F-statistic: 20.77 on 5 and 14 DF,  p-value: 5.046e-06

1/(1-summary(model.2)$r.squared) # 8.417035

## [1] 8.417035

Kada naiđemo na velike vrednosti \(VIF\) pitanje je kako popraviti situaciju. Možemo izbaciti neke prediktore iz modela. Vraćamo se na korelacionu matricu

cor(bloodpress)

##                 Pt         BP        Age     Weight         BSA       Dur
## Pt      1.00000000 0.03113499 0.04269354 0.02485650 -0.03128800 0.1762455
## BP      0.03113499 1.00000000 0.65909298 0.95006765  0.86587887 0.2928336
## Age     0.04269354 0.65909298 1.00000000 0.40734926  0.37845460 0.3437921
## Weight  0.02485650 0.95006765 0.40734926 1.00000000  0.87530481 0.2006496
## BSA    -0.03128800 0.86587887 0.37845460 0.87530481  1.00000000 0.1305400
## Dur     0.17624551 0.29283363 0.34379206 0.20064959  0.13054001 1.0000000
## Pulse   0.11228508 0.72141316 0.61876426 0.65933987  0.46481881 0.4015144
## Stress  0.34315169 0.16390139 0.36822369 0.03435475  0.01844634 0.3116398
##            Pulse     Stress
## Pt     0.1122851 0.34315169
## BP     0.7214132 0.16390139
## Age    0.6187643 0.36822369
## Weight 0.6593399 0.03435475
## BSA    0.4648188 0.01844634
## Dur    0.4015144 0.31163982
## Pulse  1.0000000 0.50631008
## Stress 0.5063101 1.00000000

Možemo izabrati da izbacimo BSA ili Weight, a pošto je i \(Pulse\) dosta korelisan sa ostalim prediktorima izbacićemo i njega. Ostaje nam model sa četiri prediktora

model.3 <- lm(BP ~ Age + Weight + Dur + Stress)
summary(model.3)

## 
## Call:
## lm(formula = BP ~ Age + Weight + Dur + Stress)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11359 -0.29586  0.01515  0.27506  0.88674 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -15.869829   3.195296  -4.967 0.000169 ***
## Age           0.683741   0.061195  11.173 1.14e-08 ***
## Weight        1.034128   0.032672  31.652 3.76e-15 ***
## Dur           0.039889   0.064486   0.619 0.545485    
## Stress        0.002184   0.003794   0.576 0.573304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5505 on 15 degrees of freedom
## Multiple R-squared:  0.9919, Adjusted R-squared:  0.9897 
## F-statistic: 458.3 on 4 and 15 DF,  p-value: 1.764e-15

#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -15.869829   3.195296  -4.967 0.000169 ***
# Age           0.683741   0.061195  11.173 1.14e-08 ***
# Weight        1.034128   0.032672  31.652 3.76e-15 ***
# Dur           0.039889   0.064486   0.619 0.545485    
# Stress        0.002184   0.003794   0.576 0.573304    
# ---
# Residual standard error: 0.5505 on 15 degrees of freedom
# Multiple R-squared:  0.9919,  Adjusted R-squared:  0.9897 
# F-statistic: 458.3 on 4 and 15 DF,  p-value: 1.764e-15

vif(model.3)

##      Age   Weight      Dur   Stress 
## 1.468245 1.234653 1.200060 1.241117

#      Age   Weight      Dur   Stress 
# 1.468245 1.234653 1.200060 1.241117 

detach(bloodpress)

Vrednosti \(VIF\) u ovom modelu su sve blizu jedan, a adjusted \(R^2\) se minimalno smanjio u odnosu na početni model, iako smo izbacili dva prediktora.

Smanjivanje strukturne multikolinearnosti

Kao što je rečeno na početku, ovde se radi o multikolinearnosti kao posledici dodavanja novih prediktora, na primer \(x\) i \(x^2\).

Primer

Ispitujemo uticaj fizičke aktivnosti na imuni sistem čoveka. Za ispitivanje ovog fenomena mogu se izabrati različiti prediktori, jedan mogući izbor je merenje zavisnosi između

\(y\) - maksimalnog unosa kiseonika (ml/kg)
\(x\) - količine imunoglobina u krvi

U bazi exerimmun.txt nalaze se podaci o 30 osoba.

exerimmun <- read.table("Multikolinearnost/exerimmun.txt", header=T)
attach(exerimmun)

plot(oxygen, igg)

Sa grafika vidimo blagu zakrivljenost u zavisnosti pa ćemo pokušati da je objasnimo kvadratnom funkcijom

\[y_i=\beta_0+\beta_1x_i+\beta_{11}x_{i}^{2}+\varepsilon_i\]

oxygensq <- oxygen^2

model.1 <- lm(igg ~ oxygen + oxygensq)

plot(x=oxygen, y=igg,
     panel.last = lines(sort(oxygen), fitted(model.1)[order(oxygen)]))

summary(model.1)

## 
## Call:
## lm(formula = igg ~ oxygen + oxygensq)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -185.375  -82.129    1.047   66.007  227.377 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1464.4042   411.4012  -3.560  0.00140 ** 
## oxygen         88.3071    16.4735   5.361 1.16e-05 ***
## oxygensq       -0.5362     0.1582  -3.390  0.00217 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 106.4 on 27 degrees of freedom
## Multiple R-squared:  0.9377, Adjusted R-squared:  0.9331 
## F-statistic: 203.2 on 2 and 27 DF,  p-value: < 2.2e-16

#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -1464.4042   411.4012  -3.560  0.00140 ** 
# oxygen         88.3071    16.4735   5.361 1.16e-05 ***
# oxygensq       -0.5362     0.1582  -3.390  0.00217 ** 
# ---
# Residual standard error: 106.4 on 27 degrees of freedom
# Multiple R-squared:  0.9377,  Adjusted R-squared:  0.9331 
# F-statistic: 203.2 on 2 and 27 DF,  p-value: < 2.2e-16

Interpretacija koeficijenata u ovakvom modelu je sledeća. \(b_0\) je vrednost zavisne promenljive ako je prediktor jednak nuli - u ovom slučaju ova interpretacija nema smisla. Koeficijent \(b_1\) je koeficijent pravca tangente u \(x=0\), a \(b_2\) govori o konkavnosti funkcije. Ispitujemo multikolinearnost modela.

vif(model.1)

##   oxygen oxygensq 
## 99.94261 99.94261

#   oxygen oxygensq 
# 99.94261 99.94261 

plot(oxygen, oxygensq)

cor(oxygen, oxygensq)

## [1] 0.9949846

Vidimo jako velike vrednosti VIF za dva prediktora. Jasno je da su ova dva prediktora jako korelisana, sa grafika i po vrednosti korelacije. Za razliku od prošlog primera, ovde smo mi krivi za multikolinearnost.

U ovom slučaju problem možemo rešiti centriranjem prediktora.

oxcent <- oxygen-mean(oxygen)
oxcentsq <- oxcent^2

plot(oxcent, oxcentsq)

cor(oxcent, oxcentsq)

## [1] 0.2195179

Ovim smo značajno smanjili korelaciju u prediktorima. Pravimo novi model sa jednačinom

\[y_i=\beta_{0}^{*}+\beta_{1}^{*}(x_i-\bar{x})+\beta_{11}^{*}(x_i-\bar{x})^{2}+\epsilon_i.\]

model.2 <- lm(igg ~ oxcent + oxcentsq)
summary(model.2)

## 
## Call:
## lm(formula = igg ~ oxcent + oxcentsq)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -185.375  -82.129    1.047   66.007  227.377 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1632.1962    29.3486   55.61  < 2e-16 ***
## oxcent        33.9995     1.6890   20.13  < 2e-16 ***
## oxcentsq      -0.5362     0.1582   -3.39  0.00217 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 106.4 on 27 degrees of freedom
## Multiple R-squared:  0.9377, Adjusted R-squared:  0.9331 
## F-statistic: 203.2 on 2 and 27 DF,  p-value: < 2.2e-16

#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 1632.1962    29.3486   55.61  < 2e-16 ***
# oxcent        33.9995     1.6890   20.13  < 2e-16 ***
# oxcentsq      -0.5362     0.1582   -3.39  0.00217 ** 
# ---
# Residual standard error: 106.4 on 27 degrees of freedom
# Multiple R-squared:  0.9377,  Adjusted R-squared:  0.9331 
# F-statistic: 203.2 on 2 and 27 DF,  p-value: < 2.2e-16

vif(model.2)

##   oxcent oxcentsq 
## 1.050628 1.050628

Vidimo kako su se zaista smanjile vrednosti VIF. Pošto smo promenili prediktore, menja se i značenje koeficijenata modela. \(b_0\) sada predstavlja vrednost zavisne promenljive kada je prediktor jednak uzoračkoj sredini \(\bar{x}\); \(b_1\) predstavlja pravac tangente u tački \(\bar{x}\); značenje \(b_2\) je isto. Na grafiku vidimo ocenjenu kvadratnu zavisnost. Kako smo napravili centriranjem dobili da prediktori više nisu puno korelisani, očekujemo da i ocenjeni koeficijent uz \(oxcent\) bude sličan u kvadratnom i u linearnom modelu. U narednom kodu vidimo da oni zaista jesu slični.

plot(x=oxcent, y=igg,
     panel.last = lines(sort(oxcent), fitted(model.2)[order(oxcent)]))

model.3 <- lm(igg ~ oxcent)
summary(model.3)

## 
## Call:
## lm(formula = igg ~ oxcent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -228.16  -79.96  -11.78   83.75  211.93 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1557.633     22.782   68.37  < 2e-16 ***
## oxcent        32.743      1.932   16.95 2.97e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 124.8 on 28 degrees of freedom
## Multiple R-squared:  0.9112, Adjusted R-squared:  0.908 
## F-statistic: 287.2 on 1 and 28 DF,  p-value: 2.973e-16

plot(x=oxcent, y=igg,
     panel.last = lines(sort(oxcent), fitted(model.3)[order(oxcent)]))

Vratimo se na dva modela sa centriranim i necentriranim prediktorima. Oni su zapravo ekvivalentni, dovode do istog grafika, prediktori su samo linearno transformisani a time su se promenili i koeficijenti. Veza između koeficijenata je, prema uvedenim oznakama, sledeća \[\begin{align} b_{0}&=b_{0}^{*}-b_{1}^{*}\bar{x}+b_{11}^{*}\bar{x}^{2}\\ b_{1}&= b_{1}^{*}-2b_{11}^{*}\bar{x}\\ b_{11}&= b_{11}^{*}\\ \end{align}\]

Testiraćemo rečeno; na osnovu koeficijenata centriranog modela računamo koeficijente necentriranog i poredimo sa vrednostima dobijenim funkcijom lm.

coef(model.2)[1]-coef(model.2)[2]*mean(oxygen)+coef(model.2)[3]*mean(oxygen)^2

## (Intercept) 
##   -1464.404

coef(model.2)[2]-2*coef(model.2)[3]*mean(oxygen)

##  oxcent 
## 88.3071

coef(model.2)[3]

##   oxcentsq 
## -0.5362473

coef(model.1)

##   (Intercept)        oxygen      oxygensq 
## -1464.4042284    88.3070970    -0.5362473

Iako su ova dva modela ekvivalentna poželjno je koristiti centrirani zbog grešaka koje su manje u ovom slučaju. Sada ispitujemo reziduale modela

plot(x=fitted(model.2), y=residuals(model.2),
     panel.last = abline(h=0, lty=2))

qqnorm(residuals(model.2), main="", datax=TRUE)
qqline(residuals(model.2), datax=TRUE)

Reziduali su dobri, možemo koristiti model za predviđanje.

predict(model.2, interval="prediction",
        newdata=data.frame(oxcent=70-mean(oxygen), oxcentsq=(70-mean(oxygen))^2))

##        fit      lwr      upr
## 1 2089.481 1849.931 2329.031

detach(exerimmun)

Napomena: Kada se prave modeli polinomijalne regresije pravilo je da se prvo fituje model višeg stepena, pa onda da se izbacuju viši stepeni. Takođe, ako je na primer, koeficijent uz kvadrat značajan, koristićemo regresionu funkciju koja sadrži i niži stepen -linearni, njega ne izbacujemo.