Špecifikácia modelu
V. Gazda
library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())Pokračujeme v predchádzajúcich úlohách, kde s pomocou regresie skúmame, nakoľko je očakávaná dĺžka života ovplyvnená priemernou dobou dochádzky do školy, úrovňou HDP na obyvateľa, respektíve priemernou úrovňou BMI v danej krajine. Výsledok tejto regresie je uvedený nižšie
#######################################################################
# PRIPRAVA UDAJOV
#######################################################################
udaje <- read.csv("udaje/Life_Expectancy_Data.csv",dec=".",sep=",",header = TRUE)
# select just the record from 2015
udaje.2015 <- udaje[udaje$Year==2015,c("Life.expectancy","BMI","GDP","Schooling")]
# data imputation
# Compute column medians
#column_medians <- sapply(udaje.2015, median, na.rm = TRUE)
# Impute missing values with column medians
# Compute column medians
column_medians <- sapply(udaje.2015, median, na.rm = TRUE)
# Impute missing values with column medians
udaje_imputed <- udaje.2015
for (col in names(udaje.2015)) {
udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}
udaje.2015 <- udaje_imputed
udaje <- udaje.2015
################################################################################
# ZAKLADNA REGRESIA
################################################################################
attach(udaje)
model <- lm(Life.expectancy ~ +1 + BMI + GDP + Schooling, data = udaje)
summary(model)
Call:
lm(formula = Life.expectancy ~ +1 + BMI + GDP + Schooling, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-17.6471 -2.4024 0.4563 3.2355 12.8591
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.419e+01 1.839e+00 24.033 <2e-16 ***
BMI 4.948e-02 2.144e-02 2.308 0.0222 *
GDP 6.972e-05 3.845e-05 1.813 0.0715 .
Schooling 1.921e+00 1.645e-01 11.676 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.033 on 179 degrees of freedom
Multiple R-squared: 0.6224, Adjusted R-squared: 0.6161
F-statistic: 98.36 on 3 and 179 DF, p-value: < 2.2e-16Testovať, či je model v správnej funkčnej forme (t. j. či je lineárna špecifikácia vhodná, alebo či by ste mali transformovať premenné, napríklad pomocou logaritmov alebo mocninami), možno vykonať viacerými spôsobmi.
1. Test RESET (test chyby špecifikácie Ramseyho regresnej rovnice - Ramsey Reset Test)
Ide o najštandardnejší formálny test nesprávnej špecifikácie funkčnej formy ale dá sa použiť aj pre prípady, ak sme nešpecifikovali všetky vysvetľujúce premenné.
Myšlienka: Nech pôvodný model má tvar \[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + u_t\] Ak je váš model správne špecifikovaný, potom pridaním mocnín vyrovnaných hodnôt (napr. \(\hat y_t^2\), \(\hat{y}_t^3\)) by sa pôvodný model nemal podstatne zlepšiť, teda budeme testovať pôvodný model uvedený vyššie a model
\[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + \gamma_2\hat y_t^2 + \gamma_3\hat{y}_t^3 + u_t\]
Budeme testovať hypotézu
\(H_0:\) model je správne špecifikovaný (\(\gamma_2 = \gamma_3 = 0\))
oproti
\(H_1:\) model je nesprávne špecifikovaný (\(\gamma_2 \ne 0 \quad \text{alebo} \quad \gamma_3 \ne 0\))
# Suppose your model is:
model <- lm(Life.expectancy ~ +1 + BMI + GDP + Schooling, data = udaje)
# RESET test from 'lmtest' package:
library(lmtest)
resettest(model)
RESET test
data: model
RESET = 6.9899, df1 = 2, df2 = 177, p-value = 0.001197Interpretácia:
Ak p-hodnota < 0.05, prijímame alternatívnu hypotézu, a teda model je zrejme zle špecifikovaný, t.j. chýbajú mu niektoré dodatočné vysvetľujúce premenné, alebo je potrebné urobiť nelineárnu transformáciu používaných premenných.
2. Grafická analýza
Graf Residuals vs. Fitted
Grafická analýza vzťahu medzi vyrovnanými hodnotami náhodnej premennej a rezíduami:
plot(model, which = 1)
Ak rezíduá vykazujú nenáhodný vzor (napr. zakrivenie), model nemusí byť lineárny v premenných a môže pomôcť nejaká ich funkčná transformácia. Mali by sme ale zvážiť aj začlenenie ďalšej premennej
Grafy C+R
Táto analýza nám môže pomôcť pri hľadaní odpovede na otázku, ktorú premennú by sme mali transfomovať pomocou nejakej známej funkcie. Vychádzajme z pôvodnej rovnice
\[y_t = \beta_0 + \beta_1 x_{t1} + \dots +\beta_k x_{tk} + u_t\] Túto rovnicu najprv odhadneme a potom vykresľujeme grafy, kde výraz component+residual (C+R) plot vykresľuje na zvislej osi \(\hat{\beta}_ix_{ti}+e_t\) a na vodorovnej osi vykresľuje hodnoty \(x_{ti}\)
Tieto grafy pomáhajú identifikovať nelineárne vzťahy pre každý regresor:
car::crPlots(model)
Ak sa krivka výrazne odchyľuje od priamky, zvážime transformáciu posudzovanej premennej (logaritmus, druhá mocnina atď.). Povedzme si, že v našm prípade je najväčší odklon od linearity v prípade premenne Schooling, ale možno aj BMI.
3. Nelineárna špecifikácia
Častokrát môžeme aj zložitejšie nelineárne vzťahy modelovať s pomocou ich aproximácie s polynómom, teda v prípade kvadratických členov
\[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + \dots + \gamma_i x_{ik}^2 + \dots + \gamma_j x_{jk}^2 + \dots + u_t\] Príklad na túto modifikáciu uvidíme nižšie.
3. Porovnanie základného a modifikovaného modelu
Predpokladajme, že sa pri nelineárnych úpravách pôvodnej rovnice dostaneme k zavedeniu kvadrátu premennej Schooling a BMI nakoľko sme motivovaní práve Component+Residual Plotom uvedenýn s vysvetlením pomocou týchto premenných, čo je uvedené vyššie. Tu sa vyrovnaná krivka asi najviac líši od priamky.
Ak má transformovaný model vyšší upravený koefcient determinácie \(R^2_{adj}\) a pri RESET test prijmeme alternatívnu hypotézu, odporúčame si výsledky potvrdiť s pomocou Anova testu oboch modelov a prípadne opakovaného Reset Testu uplatneneného na nelineárne transformovaný model
model <- lm(Life.expectancy ~ +1 + BMI + GDP + Schooling)
model_schooling_kvadr <- lm(Life.expectancy ~ +1 + BMI + GDP + Schooling + I(Schooling^2) + I(BMI^2) )
summary(model_schooling_kvadr)
Call:
lm(formula = Life.expectancy ~ +1 + BMI + GDP + Schooling + I(Schooling^2) +
I(BMI^2))
Residuals:
Min 1Q Median 3Q Max
-16.293 -2.546 0.429 3.310 12.758
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.926e+01 5.217e+00 9.443 < 2e-16 ***
BMI -2.498e-01 8.656e-02 -2.886 0.004388 **
GDP 6.582e-05 3.813e-05 1.726 0.086008 .
Schooling 1.901e+00 8.008e-01 2.374 0.018689 *
I(Schooling^2) -5.065e-03 3.187e-02 -0.159 0.873892
I(BMI^2) 3.957e-03 1.110e-03 3.563 0.000471 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.889 on 177 degrees of freedom
Multiple R-squared: 0.6477, Adjusted R-squared: 0.6378
F-statistic: 65.09 on 5 and 177 DF, p-value: < 2.2e-16anova(model, model_schooling_kvadr)Analysis of Variance Table
Model 1: Life.expectancy ~ +1 + BMI + GDP + Schooling
Model 2: Life.expectancy ~ +1 + BMI + GDP + Schooling + I(Schooling^2) +
I(BMI^2)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 179 4534.9
2 177 4231.4 2 303.6 6.3498 0.002171 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1resettest(model_schooling_kvadr)
RESET test
data: model_schooling_kvadr
RESET = 14.715, df1 = 2, df2 = 175, p-value = 1.239e-06Takto modifikovaný modeluvádza kvadrát premennej BMI ako štatisticky významný a tiež dochádza k ďalšiemu rastu kupraveného koeficientu determinácie \(R^2_{adj}\). Ak vylúčime kvadrát premennej Schooling, ktorý je štatisticky nevýznamný, potom dotaneme nasledovný výsledok
model_schooling_kvadr <- lm(Life.expectancy ~ +1 + BMI + GDP + Schooling + I(BMI^2) )
summary(model_schooling_kvadr)
Call:
lm(formula = Life.expectancy ~ +1 + BMI + GDP + Schooling + I(BMI^2))
Residuals:
Min 1Q Median 3Q Max
-16.2631 -2.5357 0.4467 3.3402 12.7876
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.000e+01 2.412e+00 20.730 < 2e-16 ***
BMI -2.496e-01 8.631e-02 -2.892 0.004309 **
GDP 6.463e-05 3.728e-05 1.734 0.084694 .
Schooling 1.776e+00 1.645e-01 10.800 < 2e-16 ***
I(BMI^2) 3.951e-03 1.107e-03 3.570 0.000459 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.876 on 178 degrees of freedom
Multiple R-squared: 0.6477, Adjusted R-squared: 0.6397
F-statistic: 81.8 on 4 and 178 DF, p-value: < 2.2e-16pričom upravený koeficient determinácie sa nezmenil. Pokračujeme teda v ďalšom využívaní tohto modelu.
5. Použitie rozšíreného RESET testu (NEPOTREBNÉ)
model_rozsireny <- lm(Life.expectancy ~ BMI + GDP + Schooling + I(BMI^2) + I(GDP^2) + I(Schooling^2) )
summary(model_rozsireny)
Call:
lm(formula = Life.expectancy ~ BMI + GDP + Schooling + I(BMI^2) +
I(GDP^2) + I(Schooling^2))
Residuals:
Min 1Q Median 3Q Max
-16.2301 -2.4487 0.5409 3.4535 12.8283
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.974e+01 5.248e+00 9.476 < 2e-16 ***
BMI -2.549e-01 8.682e-02 -2.936 0.003764 **
GDP 1.557e-04 1.101e-04 1.414 0.159259
Schooling 1.836e+00 8.047e-01 2.282 0.023696 *
I(BMI^2) 4.006e-03 1.113e-03 3.600 0.000413 ***
I(GDP^2) -1.794e-09 2.062e-09 -0.870 0.385620
I(Schooling^2) -4.089e-03 3.191e-02 -0.128 0.898172
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.893 on 176 degrees of freedom
Multiple R-squared: 0.6492, Adjusted R-squared: 0.6373
F-statistic: 54.29 on 6 and 176 DF, p-value: < 2.2e-16Tu testujeme, či koeficienty pri kvadratických členoch sú štatisticky významné, ak áno, potom musíme urobiť túto kvadratickú modifikáciu. Anova testom môžeme potom sledovať, či v porovnaní s pôvodným modelom došlo k zlepšeniu.
anova(model,model_rozsireny)Analysis of Variance Table
Model 1: Life.expectancy ~ +1 + BMI + GDP + Schooling
Model 2: Life.expectancy ~ BMI + GDP + Schooling + I(BMI^2) + I(GDP^2) +
I(Schooling^2)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 179 4534.9
2 176 4213.2 3 321.71 4.4796 0.004665 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1resettest(model_rozsireny)
RESET test
data: model_rozsireny
RESET = 14.649, df1 = 2, df2 = 174, p-value = 1.319e-06Pri reset teste tu vidíme, že problém neprávnej špecifikácie sme ešte neodstránili.
6. Transformácia pomocou dummy premennej a lineárnej lomenej funkcie
Predpokladajme, že máme dummy premennú \(D_t\), ktorá obsahuje len nuly a jedničky. Takáto premenná sa dá využiť modelovanie zlomov, a to napr.
- zlom v autonómnom člene \(\beta_0\) a to nasledovnou špecifikáciou \[y_t = \beta_0 + \beta_D D+ \beta_1 x_{t1} + \dots +\beta_k x_{tk} + u_t\] čo interpretujeme ako posun regresnej priamky (regresnej nadroviny) o \(\beta_D\) jednotiek pozdĺž zvislej osi a to len v pozorovaniach, ak je splnená podmienka \(D_t = 1\)
- zlom v sklone regresnej priamky (nadroviny) a to len v pozorovaniach, ak je splnená podmienka \(D_t = 1\), čo dosiahneme nasledovnou špecifikáciou \[y_t = \beta_0 + \beta_1 x_{t1} + \dots + \beta_{i}x_{ti} + \beta_{Di}D_tx_{ti}+ \dots + \beta_k x_{tk} + u_t\] kde teda sklon priamky pozdĺž premenne \(x_{ti}\) je \(\beta_i\) ale len v prípade \(D_t=0\), inak je ten sklon rovný \(\beta_i+\beta_{D_i}\).
Ak si pozrieme C+R graf v prípade vysvetľujúcej veličiny BMI, všimneme si, že vyrovnaná krivka nadobudne približne pre hodnoty \(BMI>27\) iný sklon. Zaveďme preto premen DUM, pre ktorú bude platiť \(DUM_t = 1\) ak \(BMI_t > 27\) a pokúsme sa analyzovať vplyv BMI separátne pre obidve úrovne.
Kvantifikujme teraz model so zlomom v autonómnom člene, ak teda \(BMI_t<27\) a teda \(DUM_t = 0\) a \(BMI_t\geq 27\) a teda \(DUM_t = 1\):
udaje$DUM <- ifelse(udaje$BMI < 27, 0, 1)
modelD_auto <- lm(Life.expectancy ~ +1 + DUM + BMI + GDP + Schooling,data=udaje )
summary(modelD_auto)
Call:
lm(formula = Life.expectancy ~ +1 + DUM + BMI + GDP + Schooling,
data = udaje)
Residuals:
Min 1Q Median 3Q Max
-18.8082 -2.6521 0.6349 3.1410 12.9335
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.451e+01 1.848e+00 24.086 <2e-16 ***
DUM 2.016e+00 1.445e+00 1.395 0.1647
BMI 1.294e-02 3.381e-02 0.383 0.7023
GDP 6.844e-05 3.836e-05 1.784 0.0761 .
Schooling 1.910e+00 1.643e-01 11.627 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.02 on 178 degrees of freedom
Multiple R-squared: 0.6265, Adjusted R-squared: 0.6181
F-statistic: 74.65 on 4 and 178 DF, p-value: < 2.2e-16V tomto prípade sa nepotvrdilo, že by existoval signifikantný zlom v posune regresnej nadroviny.
modelD_sklon <- lm(Life.expectancy ~ +1 + BMI + I(DUM*BMI) + GDP + Schooling,data=udaje )
summary(modelD_sklon)
Call:
lm(formula = Life.expectancy ~ +1 + BMI + I(DUM * BMI) + GDP +
Schooling, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-19.6911 -2.0156 0.3699 2.7548 11.3676
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.995e+01 2.134e+00 23.405 < 2e-16 ***
BMI -2.449e-01 6.633e-02 -3.692 0.000296 ***
I(DUM * BMI) 2.407e-01 5.164e-02 4.661 6.14e-06 ***
GDP 6.667e-05 3.641e-05 1.831 0.068762 .
Schooling 1.753e+00 1.599e-01 10.966 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.765 on 178 degrees of freedom
Multiple R-squared: 0.6635, Adjusted R-squared: 0.6559
F-statistic: 87.75 on 4 and 178 DF, p-value: < 2.2e-16Tentokrát upravený koeficient determinácie narástol, čo považujeme za veľmi dobrú vlastnosť modelu. V prípade krajín s nízkymi hodnotami BMI indexu môžeme konštatovať, že nárast BMI o jednu jednotku je u nich sprevádzané s poklesom očakávanej dĺžky živoat o približne 0.2449 roka a v prípade ostatných krajín nárast BMI o jednu jednotku vyvolá nárast očakávanej dĺžky života približne o -0.2449+0.2407 = 0.0042 roka. Pokúsme sa teraz otestovať významnosť zavedenia zlomu v sklone nadroviny Porovnanie pôvodného modelu a modelu s premenlivým sklonom nadroviny vykonáme s pomocou anova testu:
anova(model, modelD_sklon)Analysis of Variance Table
Model 1: Life.expectancy ~ +1 + BMI + GDP + Schooling
Model 2: Life.expectancy ~ +1 + BMI + I(DUM * BMI) + GDP + Schooling
Res.Df RSS Df Sum of Sq F Pr(>F)
1 179 4534.9
2 178 4041.6 1 493.37 21.729 6.137e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1resettest(modelD_sklon)
RESET test
data: modelD_sklon
RESET = 5.7851, df1 = 2, df2 = 176, p-value = 0.003687Na základe anova analýzy môžeme teda konštatovať, že model zlepšil svoje štatistické vlastnosti, ale Ramsey Reset test naďalej považuje chybu špecifikácie modelu za významnú.
7. Box-Coxov transformačný test (nepotrebné)
Ide o systematickejší test na zistenie, ako by sa závislá premenná mala transformovať pomocou nejakej funkcie (napr. \(\ln(Y)\), \(\sqrt(Y)\), \(1/Y\)). Konkrétna transformácia sa robí podľa vzorca
\[ Y_i^{(\lambda)} = \begin{cases} \dfrac{Y_i^{\lambda} - 1}{\lambda}, & \text{if } \lambda \neq 0, \\[10pt] \ln(Y_i), & \text{if } \lambda = 0. \end{cases} \]
library(MASS)
boxcox(model)
Interpretácia:
Maximum funkcie vierohodnosti definuje najlepšiu hodnotu \(\lambda\) zvislou bodkovanou čiarou, pričom ak:
- \(\lambda \approx 1\), potom nie je potrebná žiadna transformácia;
- \(\lambda \approx 0\), potom je potrebná logaritmická transformácia premennej \(\log(Y)\);
- \(\lambda \approx 0.5\),potom je potrebná odmocninová transformácia premennej \(\sqrt(Y)\);
- \(\lambda \approx -1\),potom je potrebná reciproká transformácia premennej \(1/Y\).
Konkrétna transformácia sa potom vykoná podľa vzorca:
\[ Y_i^{(\lambda)} = \begin{cases} \dfrac{Y_i^{\lambda} - 1}{\lambda}, & \text{if } \lambda \neq 0, \\[10pt] \ln(Y_i), & \text{if } \lambda = 0. \end{cases} \]
U nás je hodnota \(\lambda\) niekde na úrovni 1.8, teda v regresii upravíme vysvetľovanú premennú Life.expectancy a opäť spustíme regresiu.
model_lambda <- lm(I((Life.expectancy^1.8-1)/1.8) ~ +1 + BMI + GDP + Schooling, data = udaje)
summary(model_lambda)
Call:
lm(formula = I((Life.expectancy^1.8 - 1)/1.8) ~ +1 + BMI + GDP +
Schooling, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-492.26 -78.37 9.22 92.60 407.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.997e+02 5.393e+01 7.411 4.79e-12 ***
BMI 1.369e+00 6.288e-01 2.177 0.0308 *
GDP 2.304e-03 1.128e-03 2.043 0.0425 *
Schooling 5.798e+01 4.825e+00 12.016 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 147.6 on 179 degrees of freedom
Multiple R-squared: 0.6347, Adjusted R-squared: 0.6286
F-statistic: 103.7 on 3 and 179 DF, p-value: < 2.2e-16resettest(model_lambda)
RESET test
data: model_lambda
RESET = 6.9812, df1 = 2, df2 = 177, p-value = 0.001207# tu anova porovnanie modelu nemozeme pouzit, lebo vysvetľovaná premenná sa v oboch modeloch líšiTáto úprava nám síce nepatrne zvýžila \(R^2_{adj}\) ale problém chybnej špecifikácie sa neodstránil (pozri p-value reset testu). Okrem toho sme stratili príjemnú interpretáciu modelu a preto od Cox-Boxovej transformácie upustíme.