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

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

Import údajov z otvorených databáz

Pri práci v R je veľmi pohodlné využívať už hotové verejne dostupné databázy:

1. Mendeley Data – vedecké datasety.
2. Kaggle Datasets – rozsiahle databázy vrátane energetiky.
3. Databázy v R – balíky datasets, wooldridge atď.

library(datasets)
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

head(CO2)

library(wooldridge)
library(kableExtra)
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"), full_width = FALSE)

Dataset	Title
admnrev	admnrev
affairs	affairs
airfare	airfare
alcohol	alcohol
apple	apple
approval	approval
athlet1	athlet1
athlet2	athlet2
attend	attend
audit	audit
barium	barium
beauty	beauty
benefits	benefits
beveridge	beveridge
big9salary	big9salary
bwght	bwght
bwght2	bwght2
campus	campus
card	card
catholic	catholic

Import údajov z .csv

Pracujeme s databázou spotreby elektriny:

• Datetime
• meteorologické premenné (Temperature, Humidity, …)
• spotreba v troch zónach (PowerConsumption_Zone1, _Zone2, _Zone3)

udaje <- read.csv("powerconsumption.csv", header = TRUE)
head(udaje)

str(udaje)

'data.frame':   52416 obs. of  9 variables:
 $ Datetime              : chr  "1/1/2017 0:00" "1/1/2017 0:10" "1/1/2017 0:20" "1/1/2017 0:30" ...
 $ Temperature           : num  6.56 6.41 6.31 6.12 5.92 ...
 $ Humidity              : num  73.8 74.5 74.5 75 75.7 76.9 77.7 78.2 78.1 77.3 ...
 $ WindSpeed             : num  0.083 0.083 0.08 0.083 0.081 0.081 0.08 0.085 0.081 0.082 ...
 $ GeneralDiffuseFlows   : num  0.051 0.07 0.062 0.091 0.048 0.059 0.048 0.055 0.066 0.062 ...
 $ DiffuseFlows          : num  0.119 0.085 0.1 0.096 0.085 0.108 0.096 0.093 0.141 0.111 ...
 $ PowerConsumption_Zone1: num  34056 29815 29128 28229 27336 ...
 $ PowerConsumption_Zone2: num  16129 19375 19007 18361 17872 ...
 $ PowerConsumption_Zone3: num  20241 20131 19668 18899 18442 ...

colnames(udaje)

[1] "Datetime"               "Temperature"           
[3] "Humidity"               "WindSpeed"             
[5] "GeneralDiffuseFlows"    "DiffuseFlows"          
[7] "PowerConsumption_Zone1" "PowerConsumption_Zone2"
[9] "PowerConsumption_Zone3"

udaje$Datetime <- as.POSIXct(udaje$Datetime, format = "%m/%d/%Y %H:%M", tz = "UTC")
udaje$Date  <- as.Date(udaje$Datetime)
udaje$Month <- format(udaje$Datetime, "%m")
udaje$Hour  <- format(udaje$Datetime, "%H")

Grafy

library(dplyr)
library(ggplot2)

set.seed(123)

udaje.sample <- udaje %>%
  dplyr::slice_sample(n = 2000) %>%      # alebo dplyr::sample_n(2000)
  dplyr::select(
    Datetime, Temperature, Humidity, WindSpeed,
    GeneralDiffuseFlows, DiffuseFlows,
    PowerConsumption_Zone1, PowerConsumption_Zone2, PowerConsumption_Zone3
  )

head(udaje.sample)

Scatter plot

library(ggplot2)
ggplot(udaje.sample, aes(x = Temperature, y = PowerConsumption_Zone1)) +
  geom_point(alpha = 0.5) +
  theme_minimal() +
  labs(title = "Vzťah medzi teplotou a spotrebou elektriny (zóna 1)",
       x = "Teplota (°C)", y = "Spotreba v zóne 1")

Boxplot

ggplot(udaje, aes(x = factor(Hour), y = PowerConsumption_Zone1)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Spotreba elektriny podľa hodiny dňa (zóna 1)",
       x = "Hodina", y = "Spotreba zóna 1")

Základné štatistiky

library(knitr)
power.stats <- udaje %>%
  group_by(Hour) %>%
  summarise(
    n = n(),
    mean = mean(PowerConsumption_Zone1, na.rm = TRUE),
    sd = sd(PowerConsumption_Zone1, na.rm = TRUE),
    min = min(PowerConsumption_Zone1, na.rm = TRUE),
    q25 = quantile(PowerConsumption_Zone1, 0.25, na.rm = TRUE),
    median = median(PowerConsumption_Zone1, na.rm = TRUE),
    q75 = quantile(PowerConsumption_Zone1, 0.75, na.rm = TRUE),
    max = max(PowerConsumption_Zone1, na.rm = TRUE),
    .groups = "drop"
  )

kable(power.stats, digits = 2, caption = "Základné štatistiky spotreby v zóne 1 podľa hodín")

Základné štatistiky spotreby v zóne 1 podľa hodín

Hour	n	mean	sd	min	q25	median	q75	max
00	2184	30217.24	4530.07	21192.91	26855.24	28884.31	33237.01	45558.68
01	2184	27565.53	4586.72	18161.01	24289.48	25875.39	30084.44	44503.31
02	2184	26292.69	5011.50	16727.09	22923.48	24264.54	28502.54	44859.34
03	2184	25338.12	3996.48	16107.34	22202.69	24056.84	27835.98	40409.01
04	2184	24633.58	2927.81	15663.80	22113.73	24248.85	26931.80	32220.40
05	2184	23437.58	2002.03	15846.08	22009.89	23120.96	24596.41	28787.21
06	2184	23203.60	1881.38	16398.99	21948.19	23227.87	24534.74	28113.09
07	2184	24038.35	2687.71	13895.70	22351.71	23851.31	26010.25	31759.91
08	2184	26583.26	3250.58	14327.09	24424.51	26277.86	29017.56	35902.51
09	2184	29725.36	3406.73	17024.81	27509.05	29402.79	32060.71	39591.21
10	2184	32566.46	3434.23	20974.18	30229.18	32254.74	34719.90	42474.41
11	2184	34487.63	3434.61	24212.66	32118.70	34221.18	36432.63	43919.20
12	2184	35052.52	3406.06	16814.98	32713.09	34939.58	36945.30	44213.27
13	2184	35112.27	3282.76	26412.92	32899.90	34950.30	36931.86	44098.20
14	2184	34643.60	3294.29	25544.62	32460.48	34260.35	36373.79	43598.67
15	2184	33912.00	3283.99	25292.31	31723.12	33553.93	35755.77	42960.27
16	2184	33385.20	3085.33	25334.08	31356.22	32963.95	35299.37	41425.97
17	2184	34924.98	3243.09	25259.68	32464.34	34859.97	37375.43	45040.18
18	2184	38846.13	4012.83	26140.11	36199.90	39066.91	41420.39	48118.94
19	2184	42795.92	3510.00	27388.61	40123.71	43101.96	45280.80	51820.82
20	2184	43822.59	3543.75	35023.57	41912.54	44032.63	46155.79	52204.40
21	2184	42216.48	3790.12	33593.92	40015.87	42331.89	44686.96	50778.78
22	2184	39068.64	4235.92	29247.44	35823.61	38548.91	42420.07	48701.66
23	2184	34409.58	4639.96	24645.45	30900.47	33361.46	37770.49	48254.30

library(kableExtra)
power.stats %>%
  kable(digits = 2, caption = "Základné štatistiky spotreby v zóne 1 podľa hodín") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped","hover","condensed")) %>%
  column_spec(1, bold = TRUE) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%
  add_header_above(c(" " = 1, "Spotreba" = 8))

Základné štatistiky spotreby v zóne 1 podľa hodín

	Spotreba
Hour	n	mean	sd	min	q25	median	q75	max
00	2184	30217.24	4530.07	21192.91	26855.24	28884.31	33237.01	45558.68
01	2184	27565.53	4586.72	18161.01	24289.48	25875.39	30084.44	44503.31
02	2184	26292.69	5011.50	16727.09	22923.48	24264.54	28502.54	44859.34
03	2184	25338.12	3996.48	16107.34	22202.69	24056.84	27835.98	40409.01
04	2184	24633.58	2927.81	15663.80	22113.73	24248.85	26931.80	32220.40
05	2184	23437.58	2002.03	15846.08	22009.89	23120.96	24596.41	28787.21
06	2184	23203.60	1881.38	16398.99	21948.19	23227.87	24534.74	28113.09
07	2184	24038.35	2687.71	13895.70	22351.71	23851.31	26010.25	31759.91
08	2184	26583.26	3250.58	14327.09	24424.51	26277.86	29017.56	35902.51
09	2184	29725.36	3406.73	17024.81	27509.05	29402.79	32060.71	39591.21
10	2184	32566.46	3434.23	20974.18	30229.18	32254.74	34719.90	42474.41
11	2184	34487.63	3434.61	24212.66	32118.70	34221.18	36432.63	43919.20
12	2184	35052.52	3406.06	16814.98	32713.09	34939.58	36945.30	44213.27
13	2184	35112.27	3282.76	26412.92	32899.90	34950.30	36931.86	44098.20
14	2184	34643.60	3294.29	25544.62	32460.48	34260.35	36373.79	43598.67
15	2184	33912.00	3283.99	25292.31	31723.12	33553.93	35755.77	42960.27
16	2184	33385.20	3085.33	25334.08	31356.22	32963.95	35299.37	41425.97
17	2184	34924.98	3243.09	25259.68	32464.34	34859.97	37375.43	45040.18
18	2184	38846.13	4012.83	26140.11	36199.90	39066.91	41420.39	48118.94
19	2184	42795.92	3510.00	27388.61	40123.71	43101.96	45280.80	51820.82
20	2184	43822.59	3543.75	35023.57	41912.54	44032.63	46155.79	52204.40
21	2184	42216.48	3790.12	33593.92	40015.87	42331.89	44686.96	50778.78
22	2184	39068.64	4235.92	29247.44	35823.61	38548.91	42420.07	48701.66
23	2184	34409.58	4639.96	24645.45	30900.47	33361.46	37770.49	48254.30

Testovanie hypotéz

t-test – január vs. júl

t.test(udaje$PowerConsumption_Zone1[udaje$Month == "01"],
       udaje$PowerConsumption_Zone1[udaje$Month == "07"])

    Welch Two Sample t-test

data:  udaje$PowerConsumption_Zone1[udaje$Month == "01"] and udaje$PowerConsumption_Zone1[udaje$Month == "07"]
t = -31.524, df = 8895.4, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5097.477 -4500.643
sample estimates:
mean of x mean of y 
 31032.49  35831.55 

ANOVA – rozdiely medzi mesiacmi

anova.result <- aov(PowerConsumption_Zone1 ~ factor(Month), data = udaje)
summary(anova.result)

                 Df    Sum Sq   Mean Sq F value Pr(>F)    
factor(Month)    11 2.802e+11 2.547e+10   559.7 <2e-16 ***
Residuals     52404 2.385e+12 4.551e+07                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Lineárna regresia

model <- lm(PowerConsumption_Zone1 ~ Temperature + Humidity + WindSpeed + DiffuseFlows,
            data = udaje)
summary(model)

Call:
lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
    WindSpeed + DiffuseFlows, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-15476.7  -4922.1   -806.6   3977.6  17950.6 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  26513.6453   213.0782 124.431  < 2e-16 ***
Temperature    513.0160     6.1452  83.482  < 2e-16 ***
Humidity       -49.9528     2.0571 -24.283  < 2e-16 ***
WindSpeed     -142.7044    13.5866 -10.503  < 2e-16 ***
DiffuseFlows    -1.7212     0.2333  -7.379 1.62e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6357 on 52411 degrees of freedom
Multiple R-squared:  0.2052,    Adjusted R-squared:  0.2052 
F-statistic:  3383 on 4 and 52411 DF,  p-value: < 2.2e-16

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

# Your model (already fitted)
# model <- lm(ESG.INDEX ~ RETURN.ON.ASSETS + FIRM.SIZE + DEBT.TO.ASSET, data = udaje.2013)

coef.tbl <- tidy(model, conf.int = TRUE) %>%
  mutate(
    term = dplyr::recode(
      term,
      "(Intercept)"  = "Intercept",
      "Temperature"  = "Teplota",
      "Humidity"     = "Vlhkosť",
      "WindSpeed"    = "Vietor",
      "DiffuseFlows" = "Difúzne žiarenie"
    ),
    stars = 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 (ESG.INDEX ~ RETURN.ON.ASSETS + FIRM.SIZE + DEBT.TO.ASSET)"
  ) %>%
  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 (ESG.INDEX ~ RETURN.ON.ASSETS + FIRM.SIZE + DEBT.TO.ASSET)

Term	Estimate	Std. Error	t value	p value	95% CI	Sig
Intercept	26513.645	213.078	124.431	0	[26096.01, 26931.281]	***
Teplota	513.016	6.145	83.482	0	[500.971, 525.061]	***
Vlhkosť	-49.953	2.057	-24.283	0	[-53.985, -45.921]	***
Vietor	-142.704	13.587	-10.503	0	[-169.334, -116.074]	***
Difúzne žiarenie	-1.721	0.233	-7.379	0	[-2.178, -1.264]	***
Note:
Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.

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 = "Štatistiky prispôsobenia regresného modelu"
  ) %>%
  kableExtra::kable_styling(
    full_width = FALSE,
    bootstrap_options = c("condensed")
  )

Štatistiky prispôsobenia regresného modelu

R-squared	Adj. R-squared	F-statistic	F p-value	AIC	BIC	Num. obs.
0.205	0.205	3383.323	0	1066806	1066859	52416

NA