Nelineárne špecifikácie – COVID dáta

1 Úvod

V tomto dokumente odhadneme nelineárne regresné špecifikácie nad COVID dátami, podľa zadania z cvičenia 7.

Budeme porovnávať nasledujúce modely:

• Lineárny model
  
• Kvadratický model
  
• Logaritmická špecifikácia
  
• Exponenciálny model
  
• Polynomiálny model 3. stupňa

Cieľom je ukázať, ako nelinearita môže zlepšiť fit regresného modelu.

2 Balíčky

library(dplyr)
library(ggplot2)
library(lubridate)
library(broom)
library(tidyr)
library(magrittr)
library(kableExtra)
library(scales)

3 Dáta

covid <- read.csv("Covid.csv", stringsAsFactors=FALSE, sep = ";")
covid$date <- as.Date(covid$date)

glimpse(covid)

## Rows: 4,372
## Columns: 5
## $ country    <chr> "Czechia", "Czechia", "Czechia", "Czechia", "Czechia", "Cze…
## $ date       <date> 2020-01-04, 2020-01-05, 2020-01-06, 2020-01-07, 2020-01-08…
## $ cases      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ deaths     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ population <int> 10673216, 10673216, 10673216, 10673216, 10673216, 10673216,…

Po načítaní datasetu Covid.csv do prostredia R sa zobrazí jeho základná štruktúra pomocou funkcie glimpse(). Vidíme, že dataset obsahuje celkovo 4 372 riadkov a 5 premenných: country, date, cases, deaths a population.

Premenná country je typu chr, čo znamená, že obsahuje textové názvy krajín — v tomto prípade napríklad „Czechia“, pričom sa rovnaká krajina objavuje vo viacerých riadkoch zodpovedajúcich jednotlivým dňom. Premenná date je správne načítaná a teda dátumové dáta sú typu date.

Premenné cases a deaths sú typu integer, teda numerického charakteru. Z ukážky vidíme, že v úvodných dňoch pandémie majú všetky hodnoty rovnakú hodnotu 0, čo je konzistentné s tým, že v počiatočnom období sa v krajine ešte nevyskytovali potvrdené prípady ani úmrtia.

Posledná premenná population je taktiež numerická a vo všetkých riadkoch obsahuje identický počet obyvateľov danej krajiny — konkrétne približne 10,67 milióna, čo zodpovedá populácii Českej republiky.

Dataset je teda načítaný správne, jeho štruktúra je konzistentná, a na základe týchto informácií je možné pokračovať v ďalšej analýze a modelovaní nelineárnych regresných špecifikácií podľa zadania.

4 Lineárny model

m_lin <- lm(deaths ~ cases, data=covid)
summary(m_lin)

## 
## Call:
## lm(formula = deaths ~ cases, data = covid)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -37478  -2783  -1848  10053  27372 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.771e+03  2.683e+02   10.33   <2e-16 ***
## cases       1.678e-02  1.104e-04  151.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13170 on 4370 degrees of freedom
## Multiple R-squared:  0.8408, Adjusted R-squared:  0.8408 
## F-statistic: 2.308e+04 on 1 and 4370 DF,  p-value: < 2.2e-16

Interpretácia lineárneho modelu

Po odhadnutí lineárnej regresie typu
\[ \text{deaths} = \beta_0 + \beta_1 \cdot \text{cases} \]
vidíme výsledky z funkcie summary(m_lin).

Model odhaduje, že pri nulovom počte prípadov je očakávaný počet úmrtí približne
β₀ = 2771 (intercept), čo však v kontexte pandémie nemá priamu interpretačnú hodnotu, keďže nulové prípady sa vyskytovali iba na začiatku pandémie. Kľúčovým parametrom je premenná cases, ktorej koeficient je odhadnutý ako
β₁ = 0.01678.
To znamená, že každý ďalší potvrdený prípad je spojený s nárastom približne o 0.0167 úmrtia, čo zodpovedá približne 1 úmrtiu na 60 prípadov. Tento vzťah je veľmi silno štatisticky signifikantný, keďže p-hodnota je < 2e-16.

Reziduálna štandardná chyba modelu je 13170, čo naznačuje veľký rozptyl okolo regresnej priamky. Koeficient determinácie R² = 0.8408 ukazuje, že lineárna regresia dokáže vysvetliť približne 84 % variability úmrtí, čo je veľmi vysoká hodnota. Aj F-test modelu je extrémne signifikantný (p < 2e-16), čo znamená, že celkový model má výraznú vysvetľovaciu silu.

Celkovo teda lineárny model poskytuje dobrý základ na analýzu, avšak na základe veľkého rozptylu reziduí a výsledkov ďalších testov (napr. na heteroskedasticitu) je vhodné uvažovať aj o nelineárnych a transformovaných špecifikáciách modelu.

5 Kvadratický model

covid <- covid %>% mutate(cases2 = cases^2)
m_quad <- lm(deaths ~ cases + cases2, data=covid)
summary(m_quad)

## 
## Call:
## lm(formula = deaths ~ cases + cases2, data = covid)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -37332  -3227  -2239  10075  27667 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.224e+03  3.191e+02  10.102  < 2e-16 ***
## cases       1.594e-02  3.397e-04  46.910  < 2e-16 ***
## cases2      1.552e-10  5.927e-11   2.618  0.00888 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13160 on 4369 degrees of freedom
## Multiple R-squared:  0.8411, Adjusted R-squared:  0.841 
## F-statistic: 1.156e+04 on 2 and 4369 DF,  p-value: < 2.2e-16

Interpretácia kvadratického modelu

Dataset bol následne rozšírený o novú premennú cases2, ktorá predstavuje kvadratický člen v počte prípadov, teda \(\text{cases}^2\). Následne bol odhadnutý kvadratický model v tvare:

\[ \text{deaths} = \beta_0 + \beta_1 \cdot \text{cases} + \beta_2 \cdot \text{cases}^2. \]

Výstup modelu naznačuje, že všetky koeficienty sú štatisticky signifikantné, pričom člen cases má veľmi malú p-hodnotu (p < 2e-16), čo potvrdzuje jeho silný vzťah s úmrtiami. Kvadratický člen cases2 je taktiež signifikantný (p ≈ 0.0088), čo znamená, že medzi prípadmi a úmrtiami existuje aj nelineárny vzťah — konkrétne mierne zakrivenie, ktoré lineárny model nedokázal zachytiť.

Oproti jednoduchému lineárnemu modelu dochádza k miernemu zlepšeniu kvality modelu. Hodnota R-squared vzrástla z približne 0.8408 na 0.8411, čo naznačuje, že kvadratický model vysvetľuje o niečo väčší podiel variability v počte úmrtí. Aj keď je zlepšenie malé, model potvrdzuje, že nelineárna zložka má štatistický význam.

Reziduálny štandardný omyl tiež mierne klesol, čo podporuje lepšie prispôsobenie modelu.
Celkovo teda kvadratický model poskytuje mierne lepší fit ako lineárny, a to najmä vďaka tomu, že dokáže zachytiť jemné nelineárne vzory v dátach.

6 Logaritmický model

covid_log <- covid %>% filter(cases > 0)
m_log <- lm(deaths ~ log(cases), data=covid_log)
summary(m_log)

## 
## Call:
## lm(formula = deaths ~ log(cases), data = covid_log)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -26587 -18713  -9262   6178  74368 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -74368.3     1892.7  -39.29   <2e-16 ***
## log(cases)    8199.6      143.1   57.32   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24720 on 4127 degrees of freedom
## Multiple R-squared:  0.4432, Adjusted R-squared:  0.4431 
## F-statistic:  3285 on 1 and 4127 DF,  p-value: < 2.2e-16

Interpretácia logaritmický modelu

Logaritmický model skúma vzťah medzi premennými deaths a log(cases), teda medzi počtom úmrtí a logaritmicky transformovaným počtom prípadov. Takáto transformácia sa často používa v situáciách, kde pôvodné dáta vykazujú výraznú nelinearitu alebo explozívny rast, pretože logaritmus „stláča“ vysoké hodnoty a znižuje ich vplyv na regresný model.

Výstup funkcie summary(m_log) ukazuje, že model má tvar:

deaths = β₀ + β₁ · log(cases)

Pričom oba odhadnuté koeficienty sú štatisticky vysoko signifikantné (p-hodnoty < 2e-16), čo znamená, že logaritmus počtu prípadov má jednoznačne významný vplyv na počet úmrtí. Odhad parametra β₁ je približne 8199.6, čo znamená, že zvýšenie log(cases) o 1 jednotku (čo zodpovedá približne e-násobnému zvýšeniu počtu prípadov) je spojené s nárastom úmrtí o približne 8200.

Hodnota Multiple R-squared ≈ 0.4432 ukazuje, že model vysvetľuje približne 44 % variability v počte úmrtí. To je menej ako v lineárnom alebo kvadratickom modeli, čo naznačuje, že jednoduchý logaritmický tvar zachytáva iba obmedzenú časť štruktúry v dátach.

Residual standard error (≈ 24720) je výrazne vyšší než v predošlých modeloch, čo taktiež potvrdzuje slabšiu presnosť tohto modelu. Napriek tomu je logaritmická špecifikácia užitočná, pretože vizuálne odhalí zmenu šikmosti dát a poskytne odlišnú interpretáciu dynamiky medzi prípadmi a úmrtiami.

Celkovo teda logaritmický model potvrdzuje silnú štatistickú väzbu medzi premennými, ale poskytuje menej presné predikcie ako kvadratické či polynomiálne modely.

7 Exponenciálny model

covid_pos <- covid %>% filter(deaths > 0)
m_exp <- lm(log(deaths) ~ cases, data=covid_pos)
summary(m_exp)

## 
## Call:
## lm(formula = log(deaths) ~ cases, data = covid_pos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.6997 -1.0259  0.6541  1.2674  1.9086 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.699e+00  3.763e-02  204.56   <2e-16 ***
## cases       8.684e-07  1.494e-08   58.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.726 on 4063 degrees of freedom
## Multiple R-squared:  0.4541, Adjusted R-squared:  0.454 
## F-statistic:  3380 on 1 and 4063 DF,  p-value: < 2.2e-16

Interpretácia exponenciálny modelu

V tejto časti odhadujeme exponenciálny model, v ktorom je závislá premenná transformovaná pomocou prirodzeného logaritmu. Špecifikácia má tvar
log(deaths) = β₀ + β₁ * cases + ε,
čo znamená, že predpokladáme exponenciálny rast úmrtí v závislosti od počtu prípadov. Na rozdiel od predchádzajúceho logaritmického modelu sa tu loguje premenná deaths, nie cases.

Z výstupu vidíme, že model sa odhadol len na pozorovaniach, kde je počet úmrtí väčší ako nula, čo je nevyhnutná podmienka pre aplikáciu logaritmu.

Samotný výpis ukazuje:

• Intercept má hodnotu 7.6999, pričom jeho p-hodnota < 2e-16, čo znamená, že je vysoko štatisticky významný.
  
• Koeficient pri premennej cases má hodnotu 0.8684 × 10⁻⁶, taktiež s p-hodnotou < 2e-16, teda aj tento parameter je silne signifikantný.
  
• Interpretácia koeficientu je nasledovná: pri náraste cases o jednotku sa očakáva približne exp(β₁) - 1 ≈ β₁ percentuálna zmena v deaths, čo znamená veľmi malé, ale štatisticky významné zvýšenie logaritmu úmrtí.

Multiple R-squared = 0.454, čo je zreteľné zlepšenie oproti lineárnemu modelu bez transformácií. To naznačuje, že exponenciálna špecifikácia zachytáva časť variability dát lepšie ako jednoduchá lineárna regresia.

Model má tiež veľmi vysokú hodnotu F-statistiky (3380) s p-hodnotou < 2.2e-16, čo potvrdzuje, že regresný model je ako celok štatisticky významný.

Celkovo teda môžeme povedať, že exponenciálny model poskytuje lepšie prispôsobenie dátam ako základný lineárny model, pričom koeficienty sú stabilné, interpretovateľné a vysoko signifikantné.

8 Polynomiálny model (3. stupeň)

covid <- covid %>%
  mutate(
    cases2 = cases^2,
    cases3 = cases^3
  )
m_poly <- lm(deaths ~ cases + cases2 + cases3, data=covid)
summary(m_poly)

## 
## Call:
## lm(formula = deaths ~ cases + cases2 + cases3, data = covid)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23801  -5822    792   6834  41090 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.182e+02  3.079e+02  -2.658   0.0079 ** 
## cases        3.502e-02  6.383e-04  54.860   <2e-16 ***
## cases2      -9.805e-09  2.982e-10 -32.884   <2e-16 ***
## cases3       1.166e-15  3.437e-17  33.939   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11710 on 4368 degrees of freedom
## Multiple R-squared:  0.8742, Adjusted R-squared:  0.8742 
## F-statistic: 1.012e+04 on 3 and 4368 DF,  p-value: < 2.2e-16

Interpretácia polynomiálneho modelu

Polynomiálny model tretieho stupňa rozširuje základnú lineárnu špecifikáciu tým, že k pôvodnej premennej cases pridáva aj jej druhú a tretiu mocninu (cases2 a cases3). Tým umožňuje zachytiť zložitejšie, zakrivené vzťahy medzi počtom prípadov a úmrtí, ktoré jednoduchý lineárny model nedokáže vystihnúť.

Po odhade modelu deaths ~ cases + cases2 + cases3 vidíme, že všetky tri regresné koeficienty sú štatisticky vysoko signifikantné (p-hodnoty < 2e-16), čo znamená, že každá z týchto premenných významne prispieva k vysvetleniu variability v počte úmrtí.

Koeficient pri cases je kladný, pri cases2 záporný a pri cases3 opäť kladný, čo poukazuje na nelineárne zakrivenie: počet úmrtí najskôr rastie s prípadmi takmer lineárne, následne sa tempo rastu spomaľuje (vplyvom záporného členy druhého stupňa) a pri veľmi vysokých hodnotách môže opäť rásť rýchlejšie. Ide teda o typickú „vlnovitú“ polynomiálnu krivku tretieho stupňa.

Hodnota R-squared = 0.8742 je najvyššia spomedzi všetkých doteraz odhadnutých modelov, čo znamená, že tento model vysvetľuje približne 87,4 % variability úmrtí. Rovnako vysoká je aj upravená hodnota Adjusted R-squared, čo potvrdzuje, že pridané členy skutočne zlepšili kvalitu modelu a nejde iba o náhodné prispôsobenie.

Model má veľmi nízku hodnotu reziduálneho štandardného omylu (okolo 11710) vzhľadom na rozsah dát a globálny F-test je vysoko signifikantný (p-value < 2.2e-16), čo potvrdzuje celkovú štatistickú spoľahlivosť špecifikácie.

Celkovo možno povedať, že polynomiálny model tretieho stupňa poskytuje najpresnejší a najlepšie prispôsobený opis vzťahu medzi počtom prípadov COVID-19 a počtom úmrtí spomedzi všetkých testovaných modelov, a to vďaka schopnosti zachytiť nelineárne vzorce v dátach.

9 Porovnanie modelov (AIC & R²)

comparison <- tibble(
  model = c("Lineárny", "Kvadratický", "Logaritmický", "Exponenciálny", "Polynomiálny (3)"),
  AIC = c(AIC(m_lin), AIC(m_quad), AIC(m_log), AIC(m_exp), AIC(m_poly)),
  R2 = c(summary(m_lin)$r.squared,
         summary(m_quad)$r.squared,
         summary(m_log)$r.squared,
         summary(m_exp)$r.squared,
         summary(m_poly)$r.squared)
)

comparison %>%
  kable(caption="Porovnanie modelov podľa AIC a R²") %>%
  kable_classic(full_width=FALSE)

Porovnanie modelov podľa AIC a R²

model	AIC	R2
Lineárny	95354.66	0.8408266
Kvadratický	95349.80	0.8410759
Logaritmický	95254.66	0.4432058
Exponenciálny	15977.07	0.4541172
Polynomiálny (3)	94328.53	0.8742402

Interpretácia porovnania modelov

Po porovnaní jednotlivých regresných špecifikácií pomocou ukazovateľov AIC a R² vidíme výrazné rozdiely v tom, ako dobre jednotlivé modely opisujú vzťah medzi premennými cases a deaths.

Z tabuľky je zrejmé, že exponenciálny model (log(deaths) ~ cases) dosahuje najnižšiu hodnotu AIC (15977.07), čo znamená, že zo všetkých uvažovaných modelov poskytuje najlepšie vyváženie medzi kvalitou fitu a penalizáciou za počet parametrov. Súčasne však jeho R² = 0.454 je pomerne nízke, čo naznačuje, že model síce opisuje štruktúru dát efektívne, no miera vysvetlenej variability je obmedzená.

Naopak, polynomiálny model 3. stupňa (deaths ~ cases + cases² + cases³) má najvyššie R² = 0.874, čo znamená, že zo všetkých modelov zachytáva najväčšiu časť variability úmrtí. Jeho AIC (94328.53) je síce o niečo vyššie než pri kvadratickom či lineárnom modeli, no stále veľmi nízke v porovnaní s logaritmickým modelom. Tento model teda údajom vyhovuje najlepšie z hľadiska presnosti fitu, no je zároveň najzložitejší.

Lineárny model aj kvadratický model dosahujú veľmi podobné AIC a R², čo naznačuje, že v týchto dátach je nelinearita prítomná, no skôr v silnejšej (vyššej) forme, než akú umožňuje kvadratická aproximácia.

Logaritmický model má výrazne horšie AIC aj R², čo znamená, že log-transformácia cases nevystihuje povahu rastu deaths rovnako dobre ako ostatné alternatívy.

Celkovo teda môžeme povedať, že: - z hľadiska najnižšieho AIC je najlepším kandidátom exponenciálny model,
- z hľadiska najvyššieho R² je najvhodnejší polynomiálny model 3. stupňa,
- jednoduché špecifikácie (lineárny, kvadratický) poskytujú relatívne dobrý fit, ale nedokážu zachytiť silnejšie zakrivenie vzťahu,
- logaritmická transformácia nie je pre tieto dáta vhodná.

Výsledky teda potvrdzujú, že medzi počtom potvrdených prípadov a počtom úmrtí existuje jasne nelineárny vzťah, pričom najlepšie ho vystihujú vyššie polynomiálne špecifikácie alebo exponenciálne preformulovanie modelu.

10 Vizualizácia fitov

ggplot(covid, aes(cases, deaths)) +
  geom_point(alpha=0.2) +
  geom_smooth(method="lm", se=FALSE, color="blue") +
  geom_smooth(method="lm", formula=y~x+I(x^2), se=FALSE, color="red") +
  geom_smooth(method="lm", formula=y~poly(x,3,raw=TRUE), se=FALSE, color="green") +
  labs(title="Porovnanie lineárneho a nelineárnych fitov",
       x="Prípady (cases)", y="Úmrtia (deaths)",
       caption="Modrá = lineárny | Červená = kvadratický | Zelená = polynomiálny model")

Interpretácia vyzualizácie fitov

Vizualizácia porovnania lineárneho a nelineárnych modelov nám jasne ukazuje rozdiely v tom, ako dobre jednotlivé špecifikácie vystihujú vzťah medzi premennými cases a deaths. Graf obsahuje všetky reálne pozorovania (vyznačené sivými bodmi), ktoré tvoria charakteristickú, mierne zakrivenú trajektóriu úmrtí v závislosti od počtu prípadov.

Modrá krivka predstavuje lineárny model, ktorý síce vo všeobecnosti sleduje rastúci trend, ale zreteľne nedokáže zachytiť zakrivenie ani zmeny dynamiky v rôznych oblastiach rozsahu dát. V nižších hodnotách cases je príliš strmý a vo vyšších oblastiach zasa zaostáva za skutočným priebehom — model je teda systematicky nesprávne špecifikovaný.

Červená krivka (kvadratický model) už zakrivenie zachytáva oveľa lepšie. Vidíme, že tvar funkcie je bližšie k reálnym bodom a model flexibilnejšie reaguje na meniace sa tempo rastu úmrtí. Aj tak však v niektorých intervaloch buď mierne nadstreľuje, alebo podstreľuje skutočné hodnoty.

Zelená krivka (polynomiálny model 3. stupňa) je najbližšie k pozorovaným dátam. Jej zakrivenie je dynamickejšie, lepšie sleduje zrýchľovanie aj spomaľovanie rastu a najpresnejšie kopíruje reálny priebeh. Tento model má najnižší AIC spomedzi testovaných špecifikácií, čo potvrdzuje, že zvažované nelineárne riešenia poskytujú citeľne lepší fit než jednoduchá lineárna špecifikácia.

Celkovo teda graf výrazne ilustruje, že vzťah medzi prípadmi a úmrtiami nie je lineárny, ale má zreteľný nelineárny charakter, ktorý najlepšie zachytáva polynomiálna špecifikácia vyššieho stupňa.

11 Záver

• Kvadratický aj polynomiálny model lepšie zachytávajú zakrivenie vzťahu.
  
• Logaritmické a exponenciálne modely menia interpretáciu (percentuálny vs. exponenciálny rast).
  
• Porovnanie pomocou AIC a R² ukazuje, ktorý model má lepší fit.