Práca s databázou - import údajov, grafy, štatistiky

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE
)

Import údajov z otv. databáz

1. Mendeley Data Tuto sa dostaneme z Mendeley Data, kde si údaje viete voľne stiahnuť. Údaje sa vzťahujú k už publikovaným článkom vo vydavateľstve Elsevier. Výber sa dá urobiť jednoducho zadaním kľúčových slov.

2. Kaggle Data Tuto sa dostaneme z Kaggle Datasets, kde si údaje viete voľne stiahnuť. Údaje sa vzťahujú k projektom podporovaným Kaggle. Výber sa dá urobiť jednoducho zadaním kľúčových slov.

3. Databázy knižníc R – .library(datasets). alebo .library(wooldridge).ale aj iné – stačí si dať príkaz data().

library(datasets)

# dostupné datasety v 'datasets' balíku
ds <- as.data.frame(utils::data(package = "datasets")$results)[, c("Item","Title")]
knitr::kable(head(ds, 20), col.names = c("Dataset", "Title"))

Dataset	Title
AirPassengers	Monthly Airline Passenger Numbers 1949-1960
BJsales	Sales Data with Leading Indicator
BJsales.lead (BJsales)	Sales Data with Leading Indicator
BOD	Biochemical Oxygen Demand
CO2	Carbon Dioxide Uptake in Grass Plants
ChickWeight	Weight versus age of chicks on different diets
DNase	Elisa assay of DNase
EuStockMarkets	Daily Closing Prices of Major European Stock Indices, 1991-1998
Formaldehyde	Determination of Formaldehyde
HairEyeColor	Hair and Eye Color of Statistics Students
Harman23.cor	Harman Example 2.3
Harman74.cor	Harman Example 7.4
Indometh	Pharmacokinetics of Indomethacin
InsectSprays	Effectiveness of Insect Sprays
JohnsonJohnson	Quarterly Earnings per Johnson & Johnson Share
LakeHuron	Level of Lake Huron 1875-1972
LifeCycleSavings	Intercountry Life-Cycle Savings Data
Loblolly	Growth of Loblolly Pine Trees
Nile	Flow of the River Nile
Orange	Growth of Orange Trees

# príklad vstavaného datasetu
head(CO2)

Môžeme použiť aj databázu určenú pre ekonometriu – balík wooldridge:

# install.packages("wooldridge")
library(wooldridge)
ds <- as.data.frame(utils::data(package = "wooldridge")$results)[, c("Item","Title")]
knitr::kable(head(ds, 20), col.names = c("Dataset", "Title")) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE,
    position = "center"
  )

Error in knitr::kable(head(ds, 20), col.names = c("Dataset", "Title")) %>%  : 
  could not find function "%>%"

Import údajov z .csv (vlastná databáza dopravných nehôd)

V tejto práci použijem databázu o dopravných nehách, uloženú v súbore premavka.csv.csv. Dataset obsahuje informácie o dátume a čase nehody, type ovládania dopravy (semafór, značky), poveternostných podmienkach, osvetlení (deň / noc), type zrážky, stave vozovky, počte zúčastnených vozidiel a počte zranených osôb.

Súbor mám uložený v tom istom priečinku ako tento RMarkdown. Importujem ho pomocou funkcie read.csv(), pričom oddeľovač položiek je čiarka a desatinným oddeľovačom je bodka.

# Import vlastnej CSV databázy dopravných nehôd
udaje <- read.csv("premavka.csv.csv",
                  header = TRUE,
                  sep = ",",
                  dec = ".")

# náhľad na prvé riadky
head(udaje)

# názvy premenných
colnames(udaje)

 [1] "crash_date"                    "traffic_control_device"       
 [3] "weather_condition"             "lighting_condition"           
 [5] "first_crash_type"              "trafficway_type"              
 [7] "alignment"                     "roadway_surface_cond"         
 [9] "road_defect"                   "crash_type"                   
[11] "intersection_related_i"        "damage"                       
[13] "prim_contributory_cause"       "num_units"                    
[15] "most_severe_injury"            "injuries_total"               
[17] "injuries_fatal"                "injuries_incapacitating"      
[19] "injuries_non_incapacitating"   "injuries_reported_not_evident"
[21] "injuries_no_indication"        "crash_hour"                   
[23] "crash_day_of_week"             "crash_month"                  

Pre ďalšiu analýzu budem pracovať hlavne s premennými:

• weather_condition – poveternostné podmienky (CLEAR, RAIN, SNOW, FOG/SMOKE/HAZE, atď.)
• lighting_condition – osvetlenie (DAYLIGHT, DARKNESS, LIGHTED ROAD, atď.)
• num_units – počet vozidiel / jednotiek zúčastnených na nehode
• injuries_total – celkový počet zranených pri nehode
• crash_hour – hodina nehody
• crash_day_of_week – deň v týždni (1–7)
• crash_month – mesiac (1–12)

Grafy

ggplot2 – knižnica pre grafy

Najprv si vyberiem podmnožinu nehôd, napríklad iba nehody, ktoré sa stali cez deň (DAYLIGHT), aby som vedela lepšie vizualizovať vzťahy medzi premennými.

library(dplyr)
library(ggplot2)

udaje.daylight <- udaje %>%
  filter(lighting_condition == "DAYLIGHT") %>%
  select(num_units, injuries_total, crash_hour, crash_month, weather_condition)

head(udaje.daylight)

Knižnica ggplot2 je v súčasnosti najčastejšie používaná grafická knižnica v R. Predpripravené príklady grafov sú uvedené napríklad v R Graph Gallery.

Basic Scatter plot: Počet zranení podľa hodiny dňa (nehody cez deň)

ggplot(udaje.daylight, aes(x = crash_hour, y = injuries_total)) +
  geom_point(alpha = 0.6) +
  theme_minimal() +
  labs(
    title = "Počet zranení podľa hodiny dňa (nehody cez deň)",
    x = "Hodina nehody",
    y = "Celkový počet zranených"
  )

Tento graf ukazuje, ako sa počet zranených pri nehodách mení v priebehu dňa. Väčšina nehôd má nulový alebo nízky počet zranených bez ohľadu na hodinu, ale skupiny bodov naznačujú, že k závažnejším nehodám (viac zranených) môže dochádzať v popoludňajších a večerných hodinách, keď je na cestách viac áut a premávka je hustejšia.

Boxplot: Počet zranení podľa počasia

# Bar plot with grouping
library(ggplot2)

library(ggplot2)


ggplot(udaje, aes(x = factor(weather_condition), y = injuries_total)) +
  geom_boxplot(fill = "lightblue", color = "darkblue") +
  labs(
    title = "Počet zranení podľa poveternostných podmienok",
    x = "Počasie",
    y = "Celkový počet zranených"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Boxplot porovnáva rozdelenie počtu zranených pri rôznych typoch počasia. Vidno, že pri jasnom počasí (CLEAR) je medián počtu zranených nízky, ale vzhľadom na veľký počet nehôd sa občas vyskytnú aj prípady s vyšším počtom zranených. Pri daždi (RAIN), snehu (SNOW) či námraze sa rozdelenie mierne posúva, čo môže naznačovať vyššie riziko zranení pri zhoršených podmienkach, aj keď väčšina nehôd má stále 0 alebo 1 zranenú osobu.

Základné štatistiky

knitr – jednoduchá tabuľka

V tejto časti spočítam základné štatistiky pre premennú injuries_total podľa poveternostných podmienok (weather_condition). Zaujíma ma, ako sa líši počet zranených pri rôznych typoch počasia.

library(knitr)

injury.stats <- udaje %>%
  group_by(weather_condition) %>%
  summarise(
    n      = n(),
    mean   = mean(injuries_total, na.rm = TRUE),
    sd     = sd(injuries_total, na.rm = TRUE),
    min    = min(injuries_total, na.rm = TRUE),
    q25    = quantile(injuries_total, 0.25, na.rm = TRUE),
    median = median(injuries_total, na.rm = TRUE),
    q75    = quantile(injuries_total, 0.75, na.rm = TRUE),
    max    = max(injuries_total, na.rm = TRUE),
    .groups = "drop"
  )

kable(injury.stats, digits = 2,
      caption = "Základné štatistiky počtu zranených podľa poveternostných podmienok")

Základné štatistiky počtu zranených podľa poveternostných podmienok

weather_condition	n	mean	sd	min	q25	median	q75	max
BLOWING SAND, SOIL, DIRT	1	0.00	NA	0	0	0	0	0
BLOWING SNOW	127	0.31	0.75	0	0	0	0	4
CLEAR	164700	0.39	0.81	0	0	0	1	21
CLOUDY/OVERCAST	7533	0.37	0.75	0	0	0	1	11
FOG/SMOKE/HAZE	360	0.44	0.83	0	0	0	1	5
FREEZING RAIN/DRIZZLE	510	0.46	0.86	0	0	0	1	6
OTHER	627	0.48	0.82	0	0	0	1	5
RAIN	21703	0.41	0.81	0	0	0	1	15
SEVERE CROSS WIND GATE	32	0.25	0.76	0	0	0	0	4
SLEET/HAIL	308	0.46	0.85	0	0	0	1	6
SNOW	6871	0.31	0.70	0	0	0	0	7
UNKNOWN	6534	0.17	0.49	0	0	0	0	6

Tabuľka ukazuje, ako sa priemerný počet zranených a rozptyl zranení líšia medzi jednotlivými typmi počasia. Napríklad pri počasí RAIN a FREEZING RAIN/DRIZZLE sú priemerné hodnoty injuries_total mierne vyššie ako pri CLEAR, čo je v súlade s očakávaním, že zhoršené podmienky zvyšujú riziko zranení. Zároveň vidno, že pri kategórii UNKNOWN sú priemery veľmi nízke, čo môže súvisieť aj s nekvalitným alebo neúplným záznamom údajov.

library(dplyr)
library(knitr)
library(kableExtra)

# Základné štatistiky počtu zranených podľa poveternostných podmienok
injury.stats <- udaje %>%
  group_by(weather_condition) %>%
  summarise(
    n      = n(),
    mean   = mean(injuries_total, na.rm = TRUE),
    sd     = sd(injuries_total, na.rm = TRUE),
    min    = min(injuries_total, na.rm = TRUE),
    q25    = quantile(injuries_total, 0.25, na.rm = TRUE),
    median = median(injuries_total, na.rm = TRUE),
    q75    = quantile(injuries_total, 0.75, na.rm = TRUE),
    max    = max(injuries_total, na.rm = TRUE),
    .groups = "drop"
  )

# krajšia tabuľka s kableExtra
injury.stats %>%
  kable(
    digits  = 2,
    caption = "Základné štatistiky počtu zranených podľa poveternostných podmienok"
  ) %>%
  kable_styling(
    full_width        = FALSE,
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  column_spec(1, bold = TRUE) %>%                             # prvý stĺpec (počasie) tučne
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%        # hlavička tučná so sivým pozadím
  add_header_above(c(" " = 1, "Štatistiky počtu zranených" = 8))

Základné štatistiky počtu zranených podľa poveternostných podmienok

	Štatistiky počtu zranených
weather_condition	n	mean	sd	min	q25	median	q75	max
BLOWING SAND, SOIL, DIRT	1	0.00	NA	0	0	0	0	0
BLOWING SNOW	127	0.31	0.75	0	0	0	0	4
CLEAR	164700	0.39	0.81	0	0	0	1	21
CLOUDY/OVERCAST	7533	0.37	0.75	0	0	0	1	11
FOG/SMOKE/HAZE	360	0.44	0.83	0	0	0	1	5
FREEZING RAIN/DRIZZLE	510	0.46	0.86	0	0	0	1	6
OTHER	627	0.48	0.82	0	0	0	1	5
RAIN	21703	0.41	0.81	0	0	0	1	15
SEVERE CROSS WIND GATE	32	0.25	0.76	0	0	0	0	4
SLEET/HAIL	308	0.46	0.85	0	0	0	1	6
SNOW	6871	0.31	0.70	0	0	0	0	7
UNKNOWN	6534	0.17	0.49	0	0	0	0	6

NA

## Korelačná matica a heatmapa

install.packages("reshape2")

Error in install.packages : Updating loaded packages

library(dplyr)
library(ggplot2)
library(reshape2)

# výber numerických premenných, ktoré dávajú zmysel
num_data <- udaje %>%
  select(
    num_units,
    injuries_total,
    injuries_fatal,
    injuries_incapacitating,
    injuries_non_incapacitating,
    injuries_reported_not_evident,
    injuries_no_indication,
    crash_hour,
    crash_day_of_week,
    crash_month
  )

# korelačná matica (Pearsonova korelácia)
cor_mat <- cor(num_data, use = "pairwise.complete.obs")

# pre ggplot potrebujeme "dlhý" tvar
cor_df <- melt(cor_mat, varnames = c("Var1", "Var2"), value.name = "correlation")

ggplot(cor_df, aes(x = Var1, y = Var2, fill = correlation)) +
  geom_tile() +
  scale_fill_gradient2(
    limits = c(-1, 1),
    na.value = "white",
    name = "Korelácia"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  labs(
    title = "Heatmapa korelačnej matice numerických premenných"
  )

Heatmapa ukazuje, ako spolu súvisia jednotlivé číselné premenné v databáze nehôd. Vidno, že počet zranených rastie najmä vtedy, keď je do nehody zapojených viac vozidiel a keď sú zranenia vážnejšie. Časové premenné ako hodina, deň alebo mesiac majú len slabý vzťah k počtu zranených, takže čas nehody nezohráva až takú dôležitú úlohu ako samotný rozsah a závažnosť nehody.

Testovanie hypotéz

t-test: Porovnanie počtu zranených – Deň vs. Noc

Porovnám priemerný počet zranených pri nehodách, ktoré sa stali:

• cez deň (DAYLIGHT)
• v tme na osvetlenej ceste (DARKNESS, LIGHTED ROAD)

Použijem dvojvýberový t-test.

t.test.result <- t.test(
  udaje$injuries_total[udaje$lighting_condition == "DAYLIGHT"],
  udaje$injuries_total[udaje$lighting_condition == "DARKNESS, LIGHTED ROAD"]
)

print(t.test.result)

    Welch Two Sample t-test

data:  udaje$injuries_total[udaje$lighting_condition == "DAYLIGHT"] and udaje$injuries_total[udaje$lighting_condition == "DARKNESS, LIGHTED ROAD"]
t = -22.186, df = 85327, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.10675978 -0.08942783
sample estimates:
mean of x mean of y 
0.3593793 0.4574731 

Výsledok t-testu ukazuje, že priemerný počet zranených pri nehodách cez deň (DAYLIGHT) a v tme na osvetlenej ceste (DARKNESS, LIGHTED ROAD) sa štatisticky významne líši (p-hodnota je prakticky nulová). Priemerný počet zranených je vyšší pri nehodách v tme na osvetlenej ceste, čo naznačuje, že aj pri umelom osvetlení sú nočné podmienky rizikovejšie ako denné.

ANOVA: Vplyv poveternostných podmienok na počet zranených

Teraz ma zaujíma, či sa priemerný počet zranených líši medzi rôznymi typmi počasia (weather_condition). Na to použijem jednofaktorovú ANOVA.

anova.result <- aov(injuries_total ~ weather_condition, data = udaje)
summary(anova.result)

                      Df Sum Sq Mean Sq F value Pr(>F)    
weather_condition     11    388   35.24   55.25 <2e-16 ***
Residuals         209294 133474    0.64                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANOVA testuje, či je priemerný počet zranených rovnaký pre všetky typy poveternostných podmienok. Výsledky ukazujú veľmi nízku p-hodnotu (p < 0.001), čo znamená, že počasie má štatisticky významný vplyv na priemerný počet zranených. Nehovorí nám to však, ktoré konkrétne kategórie sa líšia medzi sebou – iba to, že aspoň jedna skupina (napríklad RAIN alebo SNOW) sa od ostatných významne odlišuje.

Lineárna regresia: faktory ovplyvňujúce počet zranených

V tejto časti vytvorím lineárny regresný model, ktorý sa pokúsi vysvetliť počet zranených (injuries_total) pomocou:

• num_units – počet zúčastnených vozidiel,
• crash_hour – hodina nehody,
• crash_day_of_week – deň v týždni,
• crash_month – mesiac nehody.

model <- lm(injuries_total ~ num_units + crash_hour + crash_day_of_week + crash_month,
            data = udaje)

summary(model)

Call:
lm(formula = injuries_total ~ num_units + crash_hour + crash_day_of_week + 
    crash_month, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9330 -0.3719 -0.3518  0.3171 20.9325 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.2478411  0.0110628 -22.403  < 2e-16 ***
num_units          0.3233581  0.0043563  74.227  < 2e-16 ***
crash_hour        -0.0024418  0.0003085  -7.916 2.46e-15 ***
crash_day_of_week -0.0059035  0.0008787  -6.718 1.84e-11 ***
crash_month        0.0030261  0.0005033   6.013 1.82e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7892 on 209301 degrees of freedom
Multiple R-squared:  0.02626,   Adjusted R-squared:  0.02624 
F-statistic:  1411 on 4 and 209301 DF,  p-value: < 2.2e-16

Regresný model ukazuje, ako spolu súvisí počet zranených s počtom zúčastnených jednotiek a časom nehody. Premenná num_units má kladný a štatisticky významny koeficient, čo znamená, že s rastúcim počtom vozidiel v nehode rastie aj očakávaný počet zranených. Časové premenné (hodina, deň v týždni, mesiac) síce vychádzajú štatisticky významne, ale ich koeficienty sú malé, takže zmena času má v priemere oveľa menší vplyv ako samotný rozsah nehody (počet vozidiel).

Pekná tabuľka koeficientov (broom + kableExtra)

# install.packages(c("broom", "kableExtra", "dplyr", "stringr"))
library(broom)
library(stringr)

coef.tbl <- tidy(model, conf.int = TRUE) %>%
  mutate(
    term = dplyr::recode(term,
      "(Intercept)"        = "Intercept",
      "num_units"          = "Number of units",
      "crash_hour"         = "Crash hour",
      "crash_day_of_week"  = "Day of week",
      "crash_month"        = "Month"
    ),
    stars = dplyr::case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.1   ~ "·",
      TRUE            ~ ""
    )
  ) %>%
  transmute(
    Term = term,
    Estimate = estimate,
    `Std. Error` = std.error,
    `t value` = statistic,
    `p value` = p.value,
    `95% CI` = str_c("[", round(conf.low, 3), ", ", round(conf.high, 3), "]"),
    Sig = stars
  )

coef.tbl %>%
  kable(
    digits = 3,
    caption = "OLS Regression Coefficients (injuries_total ~ num_units + crash_hour + crash_day_of_week + crash_month)"
  ) %>%
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover", "condensed")) %>%
  column_spec(1, bold = TRUE) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%
  footnote(
    general = "Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.",
    threeparttable = TRUE
  )

OLS Regression Coefficients (injuries_total ~ num_units + crash_hour + crash_day_of_week + crash_month)

Term	Estimate	Std. Error	t value	p value	95% CI	Sig
Intercept	-0.248	0.011	-22.403	0	[-0.27, -0.226]	***
Number of units	0.323	0.004	74.227	0	[0.315, 0.332]	***
Crash hour	-0.002	0.000	-7.916	0	[-0.003, -0.002]	***
Day of week	-0.006	0.001	-6.718	0	[-0.008, -0.004]	***
Month	0.003	0.001	6.013	0	[0.002, 0.004]	***
Note:
Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.

Štatistiky prispôsobenia modelu

library(broom)
library(dplyr)
library(knitr)
library(kableExtra)

fit.tbl <- glance(model) %>%
  transmute(
    `R-squared`      = r.squared,
    `Adj. R-squared` = adj.r.squared,
    `F-statistic`    = statistic,
    `F p-value`      = p.value,
    `AIC`            = AIC,
    `BIC`            = BIC,
    `Num. obs.`      = nobs
  )

fit.tbl %>%
  kable(digits = 3, caption = "Model Fit Statistics") %>%
  kable_styling(full_width = FALSE,
                bootstrap_options = c("condensed"))

Model Fit Statistics

R-squared	Adj. R-squared	F-statistic	F p-value	AIC	BIC	Num. obs.
0.026	0.026	1410.997	0	494869.2	494930.7	209306

Štatistiky prispôsobenia modelu ukazujú, že hodnota R-squared je síce pomerne nízka (okolo 0.026), čo znamená, že model vysvetľuje iba malú časť variability v počte zranených. Napriek tomu sú všetky testované premenné štatisticky významne, takže model je vhodný ako jednoduchý ilustračný nástroj na pochopenie trendov – no na presné predikcie by bolo treba zahrnúť aj ďalšie faktory (napr. typ križovatky, rýchlosť, typ vozidla a pod.).