library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
library(dplyr)
library(car)
library(tseries)Econometrics in R - cvičenie 5
Základná diagnostika modelu
Pokračujeme s databázou [Mean consumption expenditure by type of household (hbs_exp_t134)] (Eurostat, Household Budget Survey), v pripravenej podobe súboru HBS_regresia_udaje.csv.
Príprava databázy
Súbor obsahuje priemerné ročné výdavky domácností (Expenditure, PPS na domácnosť) pre dva typy domácností (jeden dospelý / dvaja dospelí) v európskych krajinách a v rokoch 1988 – 2020.
| Premenná | Význam |
|---|---|
Country |
krajina |
Year |
rok zisťovania |
Household |
typ domácnosti (One adult, Two adults) |
TwoAdults |
umelá (dummy) premenná: 1 = dvaja dospelí, 0 = jeden dospelý |
Expenditure |
priemerné ročné výdavky domácnosti |
Year_c |
centrovaný rok: Year - 2015 |
udaje <- read.csv("~/Desktop/ekonometria/HBS_regresia_udaje.csv", dec = ".", sep = ",", header = TRUE)
str(udaje)'data.frame': 348 obs. of 6 variables:
$ Country : chr "Austria" "Austria" "Austria" "Austria" ...
$ Year : int 1994 1994 1999 1999 2005 2005 2010 2010 2015 2015 ...
$ Household : chr "One adult" "Two adults" "One adult" "Two adults" ...
$ TwoAdults : int 0 1 0 1 0 1 0 1 0 1 ...
$ Expenditure: num 14226 22893 17200 26306 19530 ...
$ Year_c : int -21 -21 -16 -16 -10 -10 -5 -5 0 0 ...colSums(is.na(udaje)) # chýbajúce hodnoty Country Year Household TwoAdults Expenditure Year_c
0 0 0 0 0 0 table(udaje$Year) # počet pozorovaní podľa rokov
1988 1994 1999 2005 2010 2015 2020
18 30 54 58 60 68 60 length(unique(udaje$Country)) # počet krajín[1] 34Databáza neobsahuje chýbajúce hodnoty, preto nie je potrebná imputácia. Má 348 pozorovaní za 34 krajín.
Uloženie všetkých objektov workspace, čistenie workspace
- Aby sme nemuseli v každom projekte znova importovať a čistiť databázu, uložíme celý obsah pamäte (všetky objekty) do súboru:
save.image("pamat.Rdata")- Mazanie vybraných objektov:
aaa <- 6
bbb <- "text"
rm(list = c("aaa", "bbb"))
ls()[1] "udaje"- Vyčistenie celej pamäte:
rm(list = ls())
ls()character(0)- Opätovné načítanie uloženého workspace:
load("pamat.Rdata")load() používame na začiatku dokumentu a save.image() na konci, aby sme prenášali obsah pamäte z jedného cvičenia do druhého. Nepotrebné objekty priebežne mažeme príkazom rm().
Lineárna regresia
Odhadujeme model, ktorý vysvetľuje výdavky domácnosti typom domácnosti a časovým trendom:
\[ Expenditure_i = \beta_0 + \beta_1 TwoAdults_i + \beta_2 Year\_c_i + u_i \]
Premenná TwoAdults je umelá premenná, Year_c je rok centrovaný na 2015. Absolútny člen \(\beta_0\) má preto zmysluplnú interpretáciu: očakávané výdavky domácnosti s jedným dospelým v roku 2015.
model <- lm(Expenditure ~ 1 + TwoAdults + Year_c, data = udaje)Odhadnuté regresné koeficienty:
tidy(model) |>
mutate(
Significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
p.value < 0.10 ~ ".",
TRUE ~ ""
)
) |>
kable(
digits = 4,
col.names = c("Variable", "Estimate", "Std. Error", "t value", "p-value", "Significance"),
caption = "Regression coefficients"
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Variable | Estimate | Std. Error | t value | p-value | Significance |
|---|---|---|---|---|---|
| (Intercept) | 15225.9578 | 715.8501 | 21.2698 | 0 | *** |
| TwoAdults | 8302.3621 | 883.7075 | 9.3949 | 0 | *** |
| Year_c | 218.3794 | 47.4754 | 4.5998 | 0 | *** |
Interpretácia výsledkov regresného modelu
Funkcia broom::tidy(model) poskytuje odhady regresných koeficientov a výsledky testovania ich štatistickej významnosti.
| Stĺpec | Vysvetlenie |
|---|---|
| Variable | Premenná, ku ktorej patrí koeficient. (Intercept) je \(\hat\beta_0\). |
| Estimate | Bodový odhad \(\hat\beta_j\) – zmena očakávanej hodnoty vysvetľovanej premennej pri jednotkovej zmene vysvetľujúcej premennej, ceteris paribus. |
| Std. Error | Štandardná chyba odhadu \(SE(\hat\beta_j)\); čím menšia, tým presnejší odhad. |
| t value | Hodnota t-štatistiky na test nulovosti koeficientu. |
| p-value | Najmenšia hladina významnosti, pri ktorej by sme zamietli \(H_0\). Spravidla očakávame \(p<0.05\). |
| Significance | Hviezdičkové označenie významnosti. |
Pre každý koeficient testujeme \(H_0:\beta_j=0\) proti \(H_1:\beta_j\neq 0\) štatistikou
\[ t_j=\frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \]
ktorá má pri platnosti klasických predpokladov (vrátane normality chýb) Studentovo t-rozdelenie s \(n-k\) stupňami voľnosti (\(k\) = počet koeficientov vrátane absolútneho člena). Ak \(p<\alpha\), zamietame \(H_0\).
| Symbol | p-hodnota | Interpretácia |
|---|---|---|
*** |
\(p < 0.001\) | významný na hladine 0,1 % |
** |
\(p < 0.01\) | významný na hladine 1 % |
* |
\(p < 0.05\) | významný na hladine 5 % |
. |
\(p < 0.10\) | významný na hladine 10 % |
| bez symbolu | \(p \geq 0.10\) | nevýznamný |
Interpretácia koeficientov nášho modelu
Absolútny člen (Intercept): \(\hat\beta_0 \approx 15\,226\). Očakávané ročné výdavky domácnosti s jedným dospelým v roku 2015 (kedy je Year_c = 0) sú približne 15 226 PPS.
Koeficient pri TwoAdults: \(\hat\beta_1 \approx 8\,302\). Domácnosť s dvoma dospelými mala v rovnakom roku v priemere o približne 8 300 PPS vyššie ročné výdavky ako domácnosť s jedným dospelým (ceteris paribus). Koeficient je štatisticky významný na hladine 0,1 %.
Koeficient pri Year_c: \(\hat\beta_2 \approx 218\). S každým ďalším rokom rastú očakávané výdavky v priemere o približne 218 PPS, za predpokladu rovnakého typu domácnosti. Za 10 rokov je to zhruba 2 180 PPS. Koeficient je štatisticky významný na hladine 0,1 %.
Pozor: Údaje sú panelové – tie isté krajiny sa v dátach opakujú v rôznych rokoch. Pozorovania teda nie sú úplne nezávislé a štandardné chyby môžu byť podhodnotené. Úroveň výdavkov sa medzi krajinami (napr. Luxembursko vs. Rumunsko) líši oveľa viac než medzi rokmi, čo model s dvoma vysvetľujúcimi premennými nezachytáva.
Kvalita vyrovnania modelu
Kvalitu vyrovnania posúdime pomocou \(R^2\) a F-testu.
glance(model) |>
select(-deviance) |>
kable(
digits = 4,
col.names = c("R.squared", "Adj.R.squared", "Sigma", "F Stat.", "p.value",
"df", "logLik", "AIC", "BIC", "df.residual", "N"),
caption = "Model summary"
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| R.squared | Adj.R.squared | Sigma | F Stat. | p.value | df | logLik | AIC | BIC | df.residual | N |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.2408 | 0.2364 | 8242.675 | 54.7115 | 0 | 2 | -3630.228 | 7268.456 | 7283.865 | 345 | 348 |
Interpretácia súhrnných charakteristík
Funkcia broom::glance(model) charakterizuje model ako celok.
| Stĺpec | Vysvetlenie |
|---|---|
| R.squared | \(R^2=1-RSS/TSS\) – podiel variability vysvetľovanej premennej vysvetlený modelom. |
| Adj.R.squared | \(\bar R^2=1-(1-R^2)\frac{n-1}{n-k}\) – penalizuje počet premenných; môže byť záporný. |
| Sigma | \(\hat\sigma=\sqrt{RSS/(n-k)}\) – odhad smerodajnej odchýlky náhodnej zložky, v jednotkách vysvetľovanej premennej. |
| F Stat. | Test spoločnej významnosti: \(H_0:\beta_1=\beta_2=0\), \(F=\frac{R^2/(k-1)}{(1-R^2)/(n-k)}\) s \(k-1\) a \(n-k\) stupňami voľnosti. |
| p.value | P-hodnota celkového F-testu (nie jednotlivých koeficientov). |
| df | Počet vysvetľujúcich premenných bez absolútneho člena (\(k-1=2\)). |
| logLik | Logaritmus maximálnej vierohodnosti; použiteľný pri porovnávaní modelov s rovnakou vysvetľovanou premennou. |
| AIC | \(-2\ln L+2k\); nižšia hodnota = lepší model. |
| BIC | \(-2\ln L+k\ln n\); prísnejšie penalizuje zložitosť. |
| df.residual | \(n-k\). |
| N | Počet použitých pozorovaní. |
| deviance | \(RSS=\sum e_i^2\) (v tabuľke vynechané). |
Interpretácia pre náš model:
- \(R^2 \approx 0.24\) (upravené \(\bar R^2 \approx 0.24\)): model vysvetľuje približne 24 % variability výdavkov domácností. Zvyšných 76 % pripadá najmä na rozdiely medzi krajinami, ktoré v modeli nie sú. Upravený a obyčajný koeficient determinácie sú takmer rovnaké, lebo máme dosť pozorovaní a len dve vysvetľujúce premenné.
- \(\hat\sigma \approx 8\,243\): typická chyba modelu je približne 8 200 PPS, čo je pri priemerných výdavkoch okolo 17 800 PPS veľa.
- \(F \approx 54.7\) s \(p<0.001\): zamietame \(H_0\), že obe smernice sú nulové. Aspoň jedna vysvetľujúca premenná je významná (v našom prípade obe).
- \(df = 2\), \(df_{residual} = 348 - 3 = 345\), \(N = 348\).
Vysoké (resp. nízke) \(R^2\) samo osebe nehovorí nič o správnej špecifikácii ani o kauzalite.
Reziduálny súčet štvorcov (deviance) a reziduálnu štandardnú chybu si môžeme overiť:
nobs(model)[1] 348deviance(model)[1] 23439882950sigma(model)[1] 8242.675sqrt(deviance(model) / df.residual(model))[1] 8242.675Posledné dva výsledky sú rovnaké (až na zaokrúhlenie).
Záverečné zhrnutie: charakteristiky rozlišujeme na (1) vysvetľovaciu schopnosť (R.squared, Adj.R.squared), (2) významnosť a presnosť (Sigma, F Stat., p.value, df, df.residual), (3) porovnávanie modelov a ostatné (logLik, AIC, BIC, N, deviance).
Grafická analýza kvality modelu
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))Červená krivka v grafoch je LOWESS/LOESS vyhladenie (lokálne vážená metóda najmenších štvorcov).
Residuals vs. fitted
- Centrovanie okolo nuly: reziduály by mali kolísať okolo 0.
- Zakrivenie červenej čiary naznačuje nelinearitu.
- Rozptyl: kužeľovitý tvar by znamenal heteroskedasticitu. Pri našich dátach je to pravdepodobné, pretože v bohatších krajinách sú výdavky vyššie aj rozptýlenejšie (pozrite, či sa rezíduá s rastúcimi vyrovnanými hodnotami rozširujú).
- Odľahlé body majú pri ňom číslo riadku. Číslo nájdete cez
udaje[číslo, ].
Q-Q plot
Os X: teoretické kvantily normálneho rozdelenia; os Y: štandardizované rezíduá. Pri normalite ležia body pozdĺž 45° čiary. Rezíduá nášho modelu majú kladnú šikmosť (približne 0,62) a špicatosť približne 4,0 – pravý koniec grafu sa preto odchyľuje nahor od priamky (ťažší pravý chvost).
Scale-location plot
Os Y: druhá odmocnina absolútnych štandardizovaných rezíduí. Rovná červená čiara a rovnomerné rozloženie bodov = konštantný rozptyl (homoskedasticita). Rastúca čiara = heteroskedasticita.
Residuals vs. leverage
Páka meria, ako ďaleko je pozorovanie od stredu hodnôt vysvetľujúcich premenných. Kontúry Cookovej vzdialenosti (0,5 a 1) ukazujú vplyv na koeficienty. V našom modeli je najväčšia Cookova vzdialenosť približne 0,04, takže žiadne pozorovanie neprekračuje hranice 0,5 ani 1, a žiadne neovplyvňuje koeficienty neprimerane.
# najvplyvnejšie pozorovania
cd <- cooks.distance(model)
udaje[order(cd, decreasing = TRUE)[1:3], ] Country Year Household TwoAdults Expenditure Year_c
212 Luxembourg 2020 Two adults 1 53217 5
210 Luxembourg 2015 Two adults 1 55940 0
206 Luxembourg 2005 Two adults 1 53537 -10max(cd)[1] 0.04433336Testy normality rezíduí
Rezíduá sú lineárnou kombináciou náhodných zložiek:
\[ \mathbf e = \mathbf y - \hat{\mathbf y} = (\mathbf I - \mathbf X(\mathbf X^T\mathbf X)^{-1}\mathbf X^T)\,\mathbf u \]
Preto normalitu nepozorovateľných chýb \(\mathbf u\) testujeme sprostredkovane cez rezíduá \(\mathbf e\).
Jarqueov–Berov test
\(H_0\): rezíduá majú normálne rozdelenie; \(H_1\): nie.
\[ JB = \frac{n}{6}\left(S^2+\frac{(K-3)^2}{4}\right)\ \overset{a}{\sim}\ \chi^2(2) \]
kde \(S\) je šikmosť a \(K\) špicatosť rezíduí (pri normalite \(S=0\), \(K=3\)). Ak \(p<0.05\), normalitu zamietame; ak \(p\ge0.05\), nezamietame, no tým sme normalitu nedokázali.
res <- residuals(model)
jb_test <- tseries::jarque.bera.test(res)
tidy(jb_test) |>
select(statistic, p.value, method) |>
kable(digits = 4,
col.names = c("JB statistic", "p-value", "Test"),
caption = "Jarque-Bera test of residual normality") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| JB statistic | p-value | Test |
|---|---|---|
| 37.8142 | 0 | Jarque Bera Test |
JB štatistika vychádza približne 37,8 a p-hodnota je výrazne pod 0,001. Normalitu rezíduí zamietame.
Poznámka: Existujú aj ďalšie testy normality (Shapirov–Wilkov, Andersonov–Darlingov, Kolmogorovov–Smirnovov). JB test volíme pre jednoduchú interpretáciu cez šikmosť a špicatosť. Je asymptotický, preto je vhodné ho doplniť Q-Q grafom. Pri 348 pozorovaniach je jeho použitie primerané.
Test odľahlých pozorovaní – outlierTest()
Funkcia car::outlierTest() používa externe študentizované rezíduá
\[ t_i = \frac{e_i}{\hat{\sigma}_{(i)}\sqrt{1-h_{ii}}} \]
kde \(h_{ii}\) je páka a \(\hat\sigma_{(i)}\) odhad smerodajnej odchýlky po vynechaní \(i\)-teho pozorovania. \(H_0\): pozorovanie nie je odľahlé; \(H_1\): je odľahlé. Pri \(H_0\) má \(t_i\) Studentovo rozdelenie s \(n-k-1\) stupňami voľnosti. Pre viacnásobné testovanie sa používa Bonferroniho korekcia \(p_{i,Bonf}=\min(1,n\cdot p_i)\). Odľahlé pozorovanie identifikujeme, ak \(p_{i,Bonf}<0.05\).
outlier_test <- car::outlierTest(model)
outlier_table <- data.frame(
Observation = names(outlier_test$rstudent),
rstudent = unname(outlier_test$rstudent),
unadjusted.p = unname(outlier_test$p),
bonf.p = unname(outlier_test$bonf.p)
)
outlier_table |>
kable(digits = 4,
col.names = c("Observation", "Studentized residual",
"Unadjusted p-value", "Bonferroni p-value"),
caption = "Outlier test with Bonferroni correction") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Observation | Studentized residual | Unadjusted p-value | Bonferroni p-value |
|---|---|---|---|
| 210 | 4.0335 | 1e-04 | 0.0236 |
| 206 | 4.0016 | 1e-04 | 0.0268 |
# o ktoré pozorovanie ide
udaje[as.numeric(outlier_table$Observation), ] Country Year Household TwoAdults Expenditure Year_c
210 Luxembourg 2015 Two adults 1 55940 0
206 Luxembourg 2005 Two adults 1 53537 -10Poznámka: Odľahlé pozorovanie nemusí byť chybné a nemáme ho automaticky odstraňovať. Funkcia zobrazí len najextrémnejšie pozorovanie. Rozhoduje stĺpec
Bonferroni p-value(< 0,05). Vplyv na koeficienty posudzujeme Cookovou vzdialenosťou.
V našom modeli je najextrémnejším pozorovaním Luxembursko, 2015, dvaja dospelí (výdavky 55 940 PPS). Študentizované rezíduum je približne 4,03 a Bonferroniho p-hodnota približne 0,024, čo je pod 0,05. Pozorovanie je teda štatisticky významne odľahlé. Nejde však o chybu v dátach: Luxembursko je dlhodobo krajina s výrazne najvyššími výdavkami (podobne aj Luxembursko 2005 a 2020 s dvomi dospelými). Toto je systematický rozdiel medzi krajinami, ktorý model nezachytáva. Ako sme videli v Cookovej vzdialenosti, na odhad koeficientov to však výrazný vplyv nemá.
Úprava modelu
Rezíduá sú zošikmené doprava a rozptyl rastie s úrovňou výdavkov, preto skúsime logaritmickú transformáciu vysvetľovanej premennej. Výdavky sú kladné a pravostranne zošikmené, takže logaritmus je tu prirodzený:
\[ \ln(Expenditure_i) = \beta_0 + \beta_1 TwoAdults_i + \beta_2 Year\_c_i + u_i \]
Logaritmus mení interpretáciu koeficientov – pri malých hodnotách \(\hat\beta_j\) sú to približne percentuálne zmeny:
Year_c: za 1 rok sa výdavky zmenia o približne \(100\cdot\hat\beta_2\) %,TwoAdults: domácnosť s dvoma dospelými má výdavky vyššie o \(100\cdot(e^{\hat\beta_1}-1)\) % (pri väčších hodnotách \(\hat\beta_1\) sa presný vzorec líši od aproximácie).
model2 <- lm(log(Expenditure) ~ TwoAdults + Year_c, data = udaje)
tidy(model2) |>
kable(digits = 4,
col.names = c("Variable", "Estimate", "Std. Error", "t value", "Pr(>|t|)"),
caption = "Regression results – log(Expenditure)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Variable | Estimate | Std. Error | t value | Pr(>|t|) |
|---|---|---|---|---|
| (Intercept) | 9.4974 | 0.0439 | 216.1796 | 0 |
| TwoAdults | 0.4760 | 0.0542 | 8.7771 | 0 |
| Year_c | 0.0133 | 0.0029 | 4.5777 | 0 |
exp(coef(model2)["TwoAdults"]) - 1 # presný percentuálny rozdielTwoAdults
0.6096588 Interpretácia: \(\hat\beta_1 \approx 0.476\), teda domácnosť s dvoma dospelými má výdavky vyššie o \(e^{0.476}-1 \approx 61\) %. \(\hat\beta_2 \approx 0.0133\), takže výdavky rastú o približne 1,3 % ročne (za 10 rokov o približne 14 %). Obe premenné sú štatisticky významné na hladine 0,1 %. \(R^2 \approx 0.22\).
Grafická diagnostika upraveného modelu
par(mfrow = c(2, 2))
plot(model2)
par(mfrow = c(1, 1))Test normality rezíduí upraveného modelu
res2 <- residuals(model2)
jb_test2 <- tseries::jarque.bera.test(res2)
tidy(jb_test2) |>
select(statistic, p.value, method) |>
kable(digits = 4,
col.names = c("JB statistic", "p-value", "Test"),
caption = "Jarque-Bera test of residual normality – model 2") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| JB statistic | p-value | Test |
|---|---|---|
| 27.9249 | 0 | Jarque Bera Test |
JB štatistika klesla zo zhruba 37,8 na približne 27,9, p-hodnota je však stále pod 0,001. Logaritmická transformácia teda priblížila rozdelenie rezíduí k normálnemu, ale normalitu stále zamietame. Rezíduá sú teraz mierne zošikmené doľava (šikmosť približne −0,68), špicatosť je blízka 3 (približne 3,2).
Test odľahlých hodnôt upraveného modelu
outlier_test2 <- car::outlierTest(model2)
outlier_table2 <- data.frame(
Observation = names(outlier_test2$rstudent),
rstudent = unname(outlier_test2$rstudent),
unadjusted.p = unname(outlier_test2$p),
bonf.p = unname(outlier_test2$bonf.p)
)
outlier_table2 |>
kable(digits = 4,
col.names = c("Observation", "Studentized residual",
"Unadjusted p-value", "Bonferroni p-value"),
caption = "Outlier test – model 2") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Observation | Studentized residual | Unadjusted p-value | Bonferroni p-value |
|---|---|---|---|
| 277 | -3.4951 | 5e-04 | 0.1865 |
udaje[as.numeric(outlier_table2$Observation), ] Country Year Household TwoAdults Expenditure Year_c
277 Romania 1999 One adult 0 1903 -16Najextrémnejším pozorovaním je teraz Rumunsko, 1999, jeden dospelý (výdavky 1 903 PPS). Študentizované rezíduum je približne −3,50 a Bonferroniho p-hodnota približne 0,19. Pozorovanie už nie je významne odľahlé (p > 0,05). Logaritmická transformácia tak odstránila problém s Luxemburskom.
Malé úlohy
- Porovnajte oba modely pomocou AIC. Pozor: modely majú rôznu vysvetľovanú premennú (\(Expenditure\) vs. \(\ln Expenditure\)), preto sa ich AIC nedajú porovnať priamo. Priame hodnoty (7 268 vs. 518) by viedli k nezmyselnému záveru. Porovnanie je korektné až po prepočte AIC modelu s logaritmom na pôvodnú škálu (Jacobiho korekcia): \(AIC_{adj} = AIC(model2) + 2\sum \ln y_i\).
aic1 <- AIC(model)
aic2_adj <- AIC(model2) + 2 * sum(log(udaje$Expenditure))
c(model_levels = aic1, model_log_adj = aic2_adj) model_levels model_log_adj
7268.456 7225.814 Upravené AIC modelu s logaritmom vychádza približne 7 226, zatiaľ čo AIC modelu v úrovniach je približne 7 268. Podľa AIC je lepší logaritmický model (rozdiel je približne 43, teda veľmi výrazný).
- Urobte vlastnú modifikáciu modelu a výsledky porovnajte. Napríklad pridajte interakciu medzi typom domácnosti a časom (rozdiel medzi domácnosťami sa môže v čase meniť):
model3 <- lm(log(Expenditure) ~ TwoAdults * Year_c, data = udaje)
tidy(model3)# A tibble: 4 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 9.50 0.0489 194. 0
2 TwoAdults 0.465 0.0692 6.72 7.36e-11
3 Year_c 0.0141 0.00413 3.41 7.38e- 4
4 TwoAdults:Year_c -0.00143 0.00584 -0.245 8.06e- 1AIC(model2, model3) # rovnaká vysvetľovaná premenná – AIC je priamo porovnateľné df AIC
model2 4 518.2566
model3 5 520.1958Interakčný člen vychádza štatisticky nevýznamný (p-hodnota približne 0,35 v modeli v úrovniach) a AIC sa nezlepší, takže interakciu nezaraďujeme. Ďalšie nápady: pridať umelé premenné za krajiny (fixné efekty) alebo skúsiť iný rok / podmnožinu krajín.
Pozor: AIC/BIC porovnávame iba pri rovnakej vysvetľovanej premennej a rovnakých pozorovaniach (alebo po Jacobiho korekcii, ako vyššie).
Záver
V tejto úlohe sme kvantifikovali viacnásobnú lineárnu regresiu výdavkov domácností na type domácnosti a čase a demonštrovali základnú diagnostiku modelu (R², F-test, grafická diagnostika, test normality, test odľahlých pozorovaní).
Model v úrovniach vysvetľuje približne 24 % variability výdavkov. Domácnosť s dvoma dospelými má výdavky vyššie o približne 8 300 PPS a výdavky rastú o približne 218 PPS ročne. Diagnostika však ukázala zamietnutie normality (JB test) aj štatisticky významne odľahlé pozorovanie (Luxembursko). Logaritmická transformácia vysvetľovanej premennej tieto problémy zmiernila: odľahlé pozorovanie prestalo byť významné a po Jacobiho korekcii je podľa AIC lepší. Normalita rezíduí však zamietnutá zostala.
Spoločným obmedzením oboch modelov je panelová štruktúra dát (opakované pozorovania tých istých krajín), ktorú model s dvoma vysvetľujúcimi premennými nezachytáva. Preto by bolo vhodné v ďalšom kroku zahrnúť efekty krajín.
save.image("pamat.Rdata")
rm(list = ls())