Cieľom tejto práce je odhadnúť lineárny regresný model na vlastných dátach - basketbalové dáta euroligy a následne overiť, či sú splnené základné predpoklady klasického lineárneho modelu. Osobitne sa zameriavam na:
Model GLS odhadujeme bez ohľadu na výsledky diagnostických testov, v súlade so zadaním, aby sme demonštrovali samotnú formuláciu generalized least squares a porovnali výsledky s OLS odhadom.
library(lmtest)
library(sandwich)
library(nlme)Použitý dataset obsahuje výsledky zápasov Euroleague v časovom období rokov 2003 až 2019. Keďže dáta majú časový rozmer, pre účely diagnostiky autokorelácie sme ich usporiadali podľa dátumu a následne agregovali na dennú úroveň.
df <- read.csv("C:/Users/marti/Downloads/euroleague_dset_csv_03_19.csv", stringsAsFactors = FALSE)df$DATE <- as.Date(df$DATE, format = "%d.%m.%Y")
# ponecháme len zmysluplné pozorovania
df <- subset(df, !is.na(DATE) & P2T > 0 & HTT > 0)
# zoradenie podľa dátumu
df <- df[order(df$DATE), ]Na účely modelovania vytvoríme denný dataset. Budeme pracovať s:
mean_P2T = priemerný počet bodov v druhom polčase za
deň,mean_HTT = priemerný počet bodov v prvom polčase za
deň,games = počet zápasov v daný deň.daily_means <- aggregate(cbind(P2T, HTT) ~ DATE, data = df, FUN = mean)
daily_counts <- aggregate(P2T ~ DATE, data = df, FUN = length)
daily <- merge(daily_means, daily_counts, by = "DATE")
names(daily)[names(daily) == "P2T.x"] <- "mean_P2T"
names(daily)[names(daily) == "HTT"] <- "mean_HTT"
names(daily)[names(daily) == "P2T.y"] <- "games"
daily <- daily[order(daily$DATE), ]
daily$t <- 1:nrow(daily)
head(daily)## DATE mean_P2T mean_HTT games t
## 1 2003-11-03 62.00000 58.00000 1 1
## 2 2003-11-05 82.75000 75.25000 4 2
## 3 2003-11-06 79.71429 78.71429 7 3
## 4 2003-11-12 71.80000 69.20000 5 4
## 5 2003-11-13 75.14286 71.57143 7 5
## 6 2003-11-19 76.50000 78.16667 6 6summary(daily)## DATE mean_P2T mean_HTT games
## Min. :2003-11-03 Min. : 40.00 Min. : 55.50 Min. : 1.000
## 1st Qu.:2008-02-27 1st Qu.: 72.75 1st Qu.: 72.73 1st Qu.: 3.000
## Median :2013-01-17 Median : 76.75 Median : 76.31 Median : 4.000
## Mean :2012-05-14 Mean : 76.75 Mean : 76.54 Mean : 4.181
## 3rd Qu.:2016-04-18 3rd Qu.: 80.33 3rd Qu.: 80.62 3rd Qu.: 5.000
## Max. :2019-05-19 Max. :110.00 Max. :103.50 Max. :10.000
## t
## Min. : 1.0
## 1st Qu.:222.8
## Median :444.5
## Mean :444.5
## 3rd Qu.:666.2
## Max. :888.0Z agregovaných dát vzniklo 888 denných pozorovaní.
Budeme odhadovať nasledujúci lineárny regresný model:
\[ mean\_P2T_t = \beta_0 + \beta_1 mean\_HTT_t + \beta_2 games_t + u_t \]
kde:
m_ols <- lm(mean_P2T ~ mean_HTT + games, data = daily)
summary(m_ols)##
## Call:
## lm(formula = mean_P2T ~ mean_HTT + games, data = daily)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.635 -3.988 0.053 3.585 26.922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.8261 2.7772 21.182 < 2e-16 ***
## mean_HTT 0.2474 0.0353 7.009 4.75e-12 ***
## games -0.2414 0.1255 -1.924 0.0546 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.482 on 885 degrees of freedom
## Multiple R-squared: 0.05726, Adjusted R-squared: 0.05513
## F-statistic: 26.87 on 2 and 885 DF, p-value: 4.668e-12Odhadli sme lineárny regresný model, v ktorom sme priemerný počet
bodov v druhom polčase za deň (mean_P2T) vysvetľovali
priemerným počtom bodov v prvom polčase za deň (mean_HTT) a
počtom zápasov v daný deň (games).
Koeficient pri premennej mean_HTT je kladný a
štatisticky významný (Estimate = 0.2474, p < 0.001), čo naznačuje, že
dni s vyšším priemerným počtom bodov v prvom polčase majú tendenciu mať
aj vyšší priemerný počet bodov v druhom polčase.
Koeficient pri premennej games je záporný, avšak na 5 %
hladine významnosti nie je štatisticky významný (p = 0.0546), takže jeho
vplyv nemožno jednoznačne potvrdiť.
Celková vysvetľovacia schopnosť modelu je pomerne nízka
(R^2 = 0.0573), čo znamená, že model vysvetľuje len menšiu
časť variability premennej mean_P2T.
Breusch–Paganov test používame na overenie, či je rozptyl rezíduí konštantný.
bp_res <- bptest(m_ols)
bp_res##
## studentized Breusch-Pagan test
##
## data: m_ols
## BP = 75.991, df = 2, p-value < 2.2e-16Breusch–Paganov test poskytol hodnotu testovacej štatistiky BP = 75.991 a p-hodnotu menšiu než 2.2e-16. Na hladine významnosti 5 % preto zamietame nulovú hypotézu homoskedasticity. Usudzujeme teda, že v modeli je prítomná heteroskedasticita, čiže rozptyl rezíduí nie je konštantný.
Breusch–Godfreyho test používame na overenie prítomnosti autokorelácie rezíduí.
bg_res_1 <- bgtest(m_ols, order = 1, order.by = ~ t, data = daily)
bg_res_2 <- bgtest(m_ols, order = 2, order.by = ~ t, data = daily)
bg_res_1##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: m_ols
## LM test = 0.83934, df = 1, p-value = 0.3596bg_res_2##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: m_ols
## LM test = 2.1247, df = 2, p-value = 0.3456Breusch–Godfreyho test pre autokoreláciu do rádu 1 poskytol p-hodnotu 0.3596, a preto na hladine významnosti 5 % nulovú hypotézu neprítomnosti autokorelácie nezamietame.
Rovnako aj Breusch–Godfreyho test pre autokoreláciu do rádu 2 poskytol p-hodnotu 0.3456, takže ani v tomto prípade nulovú hypotézu neprítomnosti autokorelácie nezamietame.
Na základe týchto výsledkov sme v modeli nezistili štatisticky významnú autokoreláciu rezíduí.
Keďže Breusch–Paganov test potvrdil heteroskedasticitu, upravíme inferenciu pomocou Whiteových robustných štandardných chýb.
white_hc1 <- coeftest(m_ols, vcov = vcovHC(m_ols, type = "HC1"))
white_hc1##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.826083 3.247896 18.1121 < 2.2e-16 ***
## mean_HTT 0.247409 0.043021 5.7509 1.222e-08 ***
## games -0.241431 0.152611 -1.5820 0.114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Po aplikovaní Whiteových robustných štandardných chýb sa odhadnuté koeficienty regresie nezmenili, zmenili sa však ich štandardné chyby, t-štatistiky a p-hodnoty.
Premenná mean_HTT ostala aj po robustnej korekcii
štatisticky významná (p < 0.001), čo znamená, že jej vplyv na
mean_P2T je robustný aj pri zohľadnení
heteroskedasticity.
Naopak, premenná games po robustnej korekcii štatisticky
významná nie je (p = 0.114). To naznačuje, že klasické OLS štandardné
chyby pri tejto premennej mohli byť skreslené a pôvodná hraničná
významnosť nebola dostatočne spoľahlivá.
Pre úplnosť odhadneme aj HAC štandardné chyby, ktoré sú robustné voči heteroskedasticite aj autokorelácii.
hac_nw <- coeftest(m_ols, vcov = NeweyWest(m_ols))
hac_nw##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.826083 3.178576 18.5071 < 2.2e-16 ***
## mean_HTT 0.247409 0.042456 5.8274 7.881e-09 ***
## games -0.241431 0.136280 -1.7716 0.07681 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Aj po použití HAC (Newey–West) štandardných chýb ostali odhady koeficientov nezmenené, upravila sa však ich inferencia tak, aby bola robustná voči heteroskedasticite aj autokorelácii.
Premenná mean_HTT ostala štatisticky významná (p <
0.001), zatiaľ čo premenná games nie je štatisticky
významná ani po tejto korekcii (p = 0.0768).
Výsledky HAC štandardných chýb sú v súlade s výsledkami BP a BG testov: heteroskedasticita sa v modeli prejavila, ale autokorelácia sa výrazne nepotvrdila.
V tejto verzii modelujeme heteroskedasticitu pomocou variance
function varPower.
gls_het <- gls(
mean_P2T ~ mean_HTT + games,
data = daily,
weights = varPower(form = ~ mean_HTT)
)
summary(gls_het)## Generalized least squares fit by REML
## Model: mean_P2T ~ mean_HTT + games
## Data: daily
## AIC BIC logLik
## 5844.956 5868.884 -2917.478
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~mean_HTT
## Parameter estimates:
## power
## 0.7698487
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 59.06272 2.6622650 22.185139 0.000
## mean_HTT 0.24403 0.0345724 7.058535 0.000
## games -0.23624 0.1234803 -1.913220 0.056
##
## Correlation:
## (Intr) mn_HTT
## mean_HTT -0.978
## games -0.163 -0.031
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.794474462 -0.599397265 0.006042298 0.557809026 3.910901155
##
## Residual standard error: 0.229116
## Degrees of freedom: 888 total; 885 residualNásledne sme model odhadli metódou GLS so špecifikáciou
heteroskedasticity pomocou váh
varPower(form = ~ mean_HTT).
Odhadnutý parameter variance function (power = 0.7698)
naznačuje, že rozptyl chýb sa mení v závislosti od úrovne premennej
mean_HTT. Týmto spôsobom GLS priamo modeluje
heteroskedasticitu, a nemení teda iba štandardné chyby, ale samotný
spôsob odhadu modelu.
Aj v GLS modeli ostáva premenná mean_HTT štatisticky
významná, zatiaľ čo premenná games zostáva na 5 % hladine
významnosti nevýznamná (p = 0.056).
Spomedzi GLS modelov dosiahol tento model najnižšiu hodnotu AIC, čo naznačuje, že zo skúmaných GLS špecifikácií najlepšie vystihuje štruktúru dát.
V tejto verzii modelujeme autokoreláciu chýb pomocou štruktúry AR(1).
gls_ar1 <- gls(
mean_P2T ~ mean_HTT + games,
data = daily,
correlation = corAR1(form = ~ t)
)
summary(gls_ar1)## Generalized least squares fit by REML
## Model: mean_P2T ~ mean_HTT + games
## Data: daily
## AIC BIC logLik
## 5853.874 5877.802 -2921.937
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.03310598
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 59.17182 2.7903105 21.206176 0.0000
## mean_HTT 0.24275 0.0354691 6.843850 0.0000
## games -0.23875 0.1268518 -1.882082 0.0602
##
## Correlation:
## (Intr) mn_HTT
## mean_HTT -0.978
## games -0.218 0.029
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.804003504 -0.615512703 0.005713452 0.554130691 4.170486623
##
## Residual standard error: 6.482498
## Degrees of freedom: 888 total; 885 residualModel GLS s autokorelačnou štruktúrou AR(1) bol odhadnutý pomocou
corAR1(form = ~ t).
Odhadnutý parameter autokorelácie bol veľmi nízky
(Phi = 0.0331), čo naznačuje len minimálnu sériovú
závislosť chýb. Tento výsledok je v súlade s Breusch–Godfreyho testami,
ktoré autokoreláciu nepotvrdili.
Aj v tomto modeli ostáva premenná mean_HTT štatisticky
významná, zatiaľ čo premenná games zostáva nevýznamná (p =
0.0602).
V porovnaní s GLS modelom pre heteroskedasticitu má tento model vyššiu hodnotu AIC, a preto sa javí ako menej vhodný.
Napokon odhadneme GLS model, ktorý súčasne pripúšťa heteroskedasticitu aj autokoreláciu.
gls_both <- gls(
mean_P2T ~ mean_HTT + games,
data = daily,
correlation = corAR1(form = ~ t),
weights = varPower(form = ~ mean_HTT)
)
summary(gls_both)## Generalized least squares fit by REML
## Model: mean_P2T ~ mean_HTT + games
## Data: daily
## AIC BIC logLik
## 5846.056 5874.77 -2917.028
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.03224993
## Variance function:
## Structure: Power of variance covariate
## Formula: ~mean_HTT
## Parameter estimates:
## power
## 0.7667278
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 59.36847 2.6740129 22.202014 0.0000
## mean_HTT 0.23987 0.0347386 6.905074 0.0000
## games -0.23367 0.1248323 -1.871866 0.0616
##
## Correlation:
## (Intr) mn_HTT
## mean_HTT -0.977
## games -0.162 -0.034
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.792432822 -0.596159894 0.002494362 0.559237441 3.914109784
##
## Residual standard error: 0.232252
## Degrees of freedom: 888 total; 885 residualNapokon sme odhadli GLS model, ktorý súčasne zohľadňuje
heteroskedasticitu aj autokoreláciu, teda
corAR1(form = ~ t) a
weights = varPower(form = ~ mean_HTT).
Odhadnutý parameter autokorelácie ostal veľmi nízky
(Phi = 0.0322), zatiaľ čo parameter variance function
(power = 0.7667) potvrdzuje prítomnosť heteroskedasticity.
To znamená, že dominantným problémom v modeli je skôr heteroskedasticita
než autokorelácia.
Aj v tomto modeli je mean_HTT štatisticky významná a
games zostáva nevýznamná (p = 0.0616).
Hoci tento model umožňuje najvšeobecnejšiu špecifikáciu chýb, jeho AIC je mierne vyššie než pri GLS modeli, ktorý modeluje iba heteroskedasticitu. Preto sa ako najvhodnejší javí práve GLS model s heteroskedasticitou.
coef(m_ols)## (Intercept) mean_HTT games
## 58.8260831 0.2474087 -0.2414311coef(gls_het)## (Intercept) mean_HTT games
## 59.0627181 0.2440304 -0.2362450coef(gls_ar1)## (Intercept) mean_HTT games
## 59.1718157 0.2427451 -0.2387455coef(gls_both)## (Intercept) mean_HTT games
## 59.3684718 0.2398724 -0.2336693Pri OLS a OLS s robustnými štandardnými chybami ostávajú odhadnuté koeficienty rovnaké. Pri GLS sa však môžu meniť aj samotné koeficienty, pretože GLS zohľadňuje štruktúru chýb už pri odhade modelu.
V našom prípade sa koeficienty medzi OLS a GLS menia iba mierne, čo naznačuje, že problém modelu spočíva predovšetkým v inferencii a v heteroskedasticite, nie v zásadne nesprávnej špecifikácii lineárneho vzťahu.
Na základe výsledkov diagnostických testov sme zistili, že v modeli je prítomná heteroskedasticita, zatiaľ čo autokorelácia rezíduí sa nepotvrdila ani pri teste rádu 1, ani pri teste rádu 2.
Tomu zodpovedajú aj výsledky robustných štandardných chýb a GLS
odhadov. Whiteove robustné štandardné chyby aj HAC štandardné chyby
potvrdili, že premenná mean_HTT ostáva štatisticky
významná, zatiaľ čo premenná games nie je na 5 % hladine
významnosti robustne významná.
GLS modely sme odhadli bez ohľadu na výsledky testov, v súlade so zadaním. Spomedzi nich sa ako najvhodnejší ukázal GLS model so špecifikovanou heteroskedasticitou, pretože najlepšie zachytil štruktúru chýb a dosiahol najnižšiu hodnotu AIC. Naopak, AR(1) zložka sa ukázala ako veľmi slabá, čo je v súlade s výsledkami Breusch–Godfreyho testu.
Z výsledkov teda vyplýva, že hlavným porušením klasických predpokladov v tomto modeli je heteroskedasticita, nie autokorelácia.