Tento dokument spája viacero R Markdown súborov do jedného výstupu (v štýle publikovateľnom na RPubs). Kód je zachovaný podľa pôvodných riešení, iba mierne upravovaný tak, aby bol konzistentný v jednom spoločnom súbore (napr. odstránené duplicitné „setup“ časti a zjednotené načítanie knižníc).

Cvičenie 4 – Práca s databázou - import údajov, grafy, štatistiky

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"            "Humidity"              
## [4] "WindSpeed"              "GeneralDiffuseFlows"    "DiffuseFlows"          
## [7] "PowerConsumption_Zone1" "PowerConsumption_Zone2" "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

Cvičenie 6 – Econometrics in R - cvičenie 5 (Power consumption)

Pri práci budeme používať tieto knižnice:

Úvod do problému a hypotézy

Cieľom je vysvetliť spotrebu elektriny v zóne 1 (PowerConsumption_Zone1) pomocou meteorologických premenných a pravidelných časových efektov.

Budeme uvažovať tieto vysvetľujúce premenné:

Temperature (teplota)
Humidity (vlhkosť)
WindSpeed (rýchlosť vetra)
GeneralDiffuseFlows, DiffuseFlows (difúzne slnečné žiarenie)
časové efekty: hodina dňa a deň v týždni (spotreba má typicky silný denný a týždenný cyklus)

Pracovne očakávame, že časové premenné budú výrazne významné (zachytia denný/týždenný režim) a počasie prinesie dodatočnú vysvetľujúcu informáciu.

Príprava databázy, čistenie a úprava údajov

Dáta sú v súbore CSV. Predpokladáme, že je uložený v priečinku udaje pod názvom powerconsumption.csv.

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

## Prevod časovej premennej
udaje$Datetime <- as.POSIXct(udaje$Datetime, format = "%m/%d/%Y %H:%M")

## Rýchla kontrola štruktúry dát
str(udaje)

## 'data.frame':    52416 obs. of  9 variables:
##  $ Datetime              : POSIXct, format: "2017-01-01 00:00:00" "2017-01-01 00:10:00" ...
##  $ 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 ...

summary(udaje)

##     Datetime                    Temperature        Humidity       WindSpeed    
##  Min.   :2017-01-01 00:00:00   Min.   : 3.247   Min.   :11.34   Min.   :0.050  
##  1st Qu.:2017-04-01 23:57:30   1st Qu.:14.410   1st Qu.:58.31   1st Qu.:0.078  
##  Median :2017-07-01 23:55:00   Median :18.780   Median :69.86   Median :0.086  
##  Mean   :2017-07-01 23:55:00   Mean   :18.810   Mean   :68.26   Mean   :1.959  
##  3rd Qu.:2017-09-30 23:52:30   3rd Qu.:22.890   3rd Qu.:81.40   3rd Qu.:4.915  
##  Max.   :2017-12-30 23:50:00   Max.   :40.010   Max.   :94.80   Max.   :6.483  
##  GeneralDiffuseFlows  DiffuseFlows     PowerConsumption_Zone1
##  Min.   :   0.004    Min.   :  0.011   Min.   :13896         
##  1st Qu.:   0.062    1st Qu.:  0.122   1st Qu.:26311         
##  Median :   5.035    Median :  4.456   Median :32266         
##  Mean   : 182.697    Mean   : 75.028   Mean   :32345         
##  3rd Qu.: 319.600    3rd Qu.:101.000   3rd Qu.:37309         
##  Max.   :1163.000    Max.   :936.000   Max.   :52204         
##  PowerConsumption_Zone2 PowerConsumption_Zone3
##  Min.   : 8560          Min.   : 5935         
##  1st Qu.:16981          1st Qu.:13129         
##  Median :20823          Median :16415         
##  Mean   :21043          Mean   :17835         
##  3rd Qu.:24714          3rd Qu.:21624         
##  Max.   :37409          Max.   :47598

Výber premenných a tvorba časových premenných

Keďže ide o časový rad, pridáme: - hodinu dňa (Hour) - deň v týždni (DayOfWeek)

Aby boli výstupy prehľadnejšie, v tomto cvičení budeme pracovať s jedným mesiacom (jún 2017). Ak chceš pracovať s celým rokom, filter jednoducho zakomentuj.

udaje_model <- udaje[, c("Datetime",
                         "PowerConsumption_Zone1",
                         "Temperature", "Humidity", "WindSpeed",
                         "GeneralDiffuseFlows", "DiffuseFlows")]

udaje_model$Hour <- as.integer(format(udaje_model$Datetime, "%H"))
udaje_model$DayOfWeek <- factor(format(udaje_model$Datetime, "%a"),
                                levels = c("Mon","Tue","Wed","Thu","Fri","Sat","Sun"))

## Filter na jún 2017 (voliteľné)
udaje_model <- subset(
  udaje_model,
  Datetime >= as.POSIXct("2017-06-01 00:00:00") &
    Datetime <= as.POSIXct("2017-06-30 23:50:00")
)

## Kontrola chýbajúcich hodnôt
colSums(is.na(udaje_model))

##               Datetime PowerConsumption_Zone1            Temperature 
##                      0                      0                      0 
##               Humidity              WindSpeed    GeneralDiffuseFlows 
##                      0                      0                      0 
##           DiffuseFlows                   Hour              DayOfWeek 
##                      0                      0                      0

Imputácia chýbajúcich hodnôt (ak by sa vyskytli)

V tejto databáze bývajú chýbajúce hodnoty zriedkavé. Napriek tomu je dobré mať korektný postup. Ak sa niekde vyskytne NA, doplníme ho mediánom danej numerickej premennej.

if (anyNA(udaje_model)) {
  num_cols <- sapply(udaje_model, is.numeric)
  med <- sapply(udaje_model[, num_cols, drop = FALSE], median, na.rm = TRUE)

  for (cn in names(med)) {
    udaje_model[[cn]][is.na(udaje_model[[cn]])] <- med[cn]
  }
}

Kontrola tvaru údajov (boxploty)

Boxploty nám pomôžu rýchlo identifikovať extrémne hodnoty alebo neštandardné rozdelenia.

vars <- c("PowerConsumption_Zone1", "Temperature", "Humidity", "WindSpeed",
          "GeneralDiffuseFlows", "DiffuseFlows")

par(mfrow = c(2, 3))
par(mar = c(4, 4, 2, 1))

for (v in vars) {
  boxplot(udaje_model[[v]],
          main = v,
          ylab = "Hodnota",
          col = "lightblue")
}

par(mfrow = c(1, 1))

Lineárna regresia

Odhadneme lineárny model pomocou lm().

V prvom modeli použijeme všetky meteorologické premenné a pridáme časové efekty (hodina a deň v týždni). Časové efekty modelujeme ako faktorové premenné, aby model vedel zachytiť rozdielne priemery v jednotlivých hodinách/dňoch.

model1 <- lm(
  PowerConsumption_Zone1 ~ Temperature + Humidity + WindSpeed +
    GeneralDiffuseFlows + DiffuseFlows +
    factor(Hour) + DayOfWeek,
  data = udaje_model
)

summary(model1)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows + factor(Hour) + 
##     DayOfWeek, data = udaje_model)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13732.4  -1228.7     19.2   1470.4   7187.2 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          4.770e+04  6.895e+02  69.176  < 2e-16 ***
## Temperature         -1.875e+02  2.804e+01  -6.688 2.56e-11 ***
## Humidity            -5.302e+01  3.352e+00 -15.818  < 2e-16 ***
## WindSpeed           -2.844e+02  1.919e+01 -14.819  < 2e-16 ***
## GeneralDiffuseFlows  1.980e+00  3.552e-01   5.574 2.64e-08 ***
## DiffuseFlows        -3.643e+00  4.106e-01  -8.872  < 2e-16 ***
## factor(Hour)1       -1.324e+03  2.567e+02  -5.157 2.62e-07 ***
## factor(Hour)2       -8.817e+02  2.570e+02  -3.431 0.000607 ***
## factor(Hour)3       -6.175e+03  2.577e+02 -23.965  < 2e-16 ***
## factor(Hour)4       -1.199e+04  2.583e+02 -46.409  < 2e-16 ***
## factor(Hour)5       -1.655e+04  2.589e+02 -63.910  < 2e-16 ***
## factor(Hour)6       -1.698e+04  2.591e+02 -65.529  < 2e-16 ***
## factor(Hour)7       -1.590e+04  2.661e+02 -59.733  < 2e-16 ***
## factor(Hour)8       -1.381e+04  3.006e+02 -45.941  < 2e-16 ***
## factor(Hour)9       -1.144e+04  3.474e+02 -32.920  < 2e-16 ***
## factor(Hour)10      -9.447e+03  3.696e+02 -25.562  < 2e-16 ***
## factor(Hour)11      -7.655e+03  3.712e+02 -20.620  < 2e-16 ***
## factor(Hour)12      -6.561e+03  3.795e+02 -17.288  < 2e-16 ***
## factor(Hour)13      -5.888e+03  3.831e+02 -15.367  < 2e-16 ***
## factor(Hour)14      -5.350e+03  3.668e+02 -14.587  < 2e-16 ***
## factor(Hour)15      -5.289e+03  3.488e+02 -15.165  < 2e-16 ***
## factor(Hour)16      -5.390e+03  3.294e+02 -16.364  < 2e-16 ***
## factor(Hour)17      -4.018e+03  3.166e+02 -12.690  < 2e-16 ***
## factor(Hour)18      -1.297e+03  3.077e+02  -4.213 2.57e-05 ***
## factor(Hour)19       3.778e+03  2.749e+02  13.741  < 2e-16 ***
## factor(Hour)20       5.556e+03  2.634e+02  21.096  < 2e-16 ***
## factor(Hour)21       5.157e+03  2.595e+02  19.876  < 2e-16 ***
## factor(Hour)22       4.926e+03  2.576e+02  19.124  < 2e-16 ***
## factor(Hour)23       2.986e+03  2.568e+02  11.629  < 2e-16 ***
## DayOfWeekTue        -4.510e+02  1.446e+02  -3.119 0.001829 ** 
## DayOfWeekWed         1.452e+02  1.480e+02   0.981 0.326639    
## DayOfWeekThu         1.139e+02  1.419e+02   0.803 0.421980    
## DayOfWeekFri         3.682e+02  1.414e+02   2.603 0.009275 ** 
## DayOfWeekSat         1.329e+03  1.445e+02   9.197  < 2e-16 ***
## DayOfWeekSun        -7.597e+02  1.440e+02  -5.276 1.38e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2435 on 4285 degrees of freedom
## Multiple R-squared:  0.8902, Adjusted R-squared:  0.8893 
## F-statistic:  1021 on 34 and 4285 DF,  p-value: < 2.2e-16

Diagnostické grafy

Štandardné diagnostické grafy nám pomáhajú posúdiť: - nelinearity, - heteroskedasticitu, - odľahlé a vplyvné pozorovania, - približnú normalitu rezíduí.

par(mfrow = c(2, 2))
plot(model1)

Diagnostické grafy regresného modelu (model1)

par(mfrow = c(1, 1))

Testy: normalita rezíduí a odľahlé pozorovania

Pri veľkých vzorkách je bežné, že test normality vyjde ako „zamietnutý“ aj pri menších odchýlkach. Preto výsledok testu vždy interpretujeme spolu s Q-Q grafom.

res1 <- residuals(model1)

jarque.bera.test(res1)

## 
##  Jarque Bera Test
## 
## data:  res1
## X-squared = 1187.6, df = 2, p-value < 2.2e-16

outlierTest(model1)

##        rstudent unadjusted p-value Bonferroni p
## 25459 -5.685671         1.3897e-08   6.0033e-05
## 25460 -5.421674         6.2293e-08   2.6910e-04
## 25603 -5.119478         3.1982e-07   1.3816e-03
## 25461 -5.079365         3.9480e-07   1.7056e-03
## 25604 -4.874049         1.1329e-06   4.8941e-03
## 25462 -4.709636         2.5601e-06   1.1060e-02

Úprava modelu: log-transformácia žiarenia

Premenné žiarenia (GeneralDiffuseFlows, DiffuseFlows) bývajú výrazne asymetrické (mnoho nízkych hodnôt, občas vysoké). Použijeme transformáciu log(1 + x) (v R: log1p(x)), aby sme sa vyhli problému log(0).

model2 <- lm(
  PowerConsumption_Zone1 ~ Temperature + Humidity + WindSpeed +
    I(log1p(GeneralDiffuseFlows)) + I(log1p(DiffuseFlows)) +
    factor(Hour) + DayOfWeek,
  data = udaje_model
)

summary(model2)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + I(log1p(GeneralDiffuseFlows)) + I(log1p(DiffuseFlows)) + 
##     factor(Hour) + DayOfWeek, data = udaje_model)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13673.6  -1229.7     23.5   1471.0   7188.5 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    47442.32     679.80  69.789  < 2e-16 ***
## Temperature                     -175.44      27.57  -6.363 2.18e-10 ***
## Humidity                         -52.36       3.36 -15.585  < 2e-16 ***
## WindSpeed                       -282.95      19.20 -14.737  < 2e-16 ***
## I(log1p(GeneralDiffuseFlows))    488.64      96.11   5.084 3.85e-07 ***
## I(log1p(DiffuseFlows))          -771.21      72.87 -10.583  < 2e-16 ***
## factor(Hour)1                  -1318.24     256.72  -5.135 2.95e-07 ***
## factor(Hour)2                   -872.53     256.95  -3.396 0.000691 ***
## factor(Hour)3                  -6164.95     257.63 -23.929  < 2e-16 ***
## factor(Hour)4                 -11968.97     258.20 -46.354  < 2e-16 ***
## factor(Hour)5                 -16515.61     259.14 -63.732  < 2e-16 ***
## factor(Hour)6                 -16439.43     337.58 -48.699  < 2e-16 ***
## factor(Hour)7                 -14934.54     453.37 -32.942  < 2e-16 ***
## factor(Hour)8                 -12797.71     520.90 -24.568  < 2e-16 ***
## factor(Hour)9                 -10443.56     556.63 -18.762  < 2e-16 ***
## factor(Hour)10                 -8394.69     570.70 -14.710  < 2e-16 ***
## factor(Hour)11                 -6356.62     575.51 -11.045  < 2e-16 ***
## factor(Hour)12                 -5197.98     577.35  -9.003  < 2e-16 ***
## factor(Hour)13                 -4516.20     580.21  -7.784 8.77e-15 ***
## factor(Hour)14                 -4052.55     571.55  -7.090 1.56e-12 ***
## factor(Hour)15                 -4106.16     562.12  -7.305 3.29e-13 ***
## factor(Hour)16                 -4350.30     550.03  -7.909 3.27e-15 ***
## factor(Hour)17                 -3012.67     534.24  -5.639 1.82e-08 ***
## factor(Hour)18                  -293.25     516.47  -0.568 0.570197    
## factor(Hour)19                  4786.46     451.68  10.597  < 2e-16 ***
## factor(Hour)20                  5986.12     307.41  19.473  < 2e-16 ***
## factor(Hour)21                  5125.05     259.31  19.764  < 2e-16 ***
## factor(Hour)22                  4904.81     257.53  19.045  < 2e-16 ***
## factor(Hour)23                  2979.92     256.78  11.605  < 2e-16 ***
## DayOfWeekTue                    -430.79     144.72  -2.977 0.002930 ** 
## DayOfWeekWed                     196.60     148.22   1.326 0.184779    
## DayOfWeekThu                     125.31     141.51   0.886 0.375922    
## DayOfWeekFri                     396.80     141.02   2.814 0.004920 ** 
## DayOfWeekSat                    1331.95     144.37   9.226  < 2e-16 ***
## DayOfWeekSun                    -771.53     143.96  -5.359 8.79e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2435 on 4285 degrees of freedom
## Multiple R-squared:  0.8902, Adjusted R-squared:  0.8893 
## F-statistic:  1021 on 34 and 4285 DF,  p-value: < 2.2e-16

par(mfrow = c(2, 2))
plot(model2)

Diagnostické grafy regresného modelu (model2)

par(mfrow = c(1, 1))

res2 <- residuals(model2)

jarque.bera.test(res2)

## 
##  Jarque Bera Test
## 
## data:  res2
## X-squared = 1211.7, df = 2, p-value < 2.2e-16

outlierTest(model2)

##        rstudent unadjusted p-value Bonferroni p
## 25459 -5.661457         1.5991e-08   0.00006908
## 25460 -5.392039         7.3413e-08   0.00031715
## 25603 -5.105855         3.4360e-07   0.00148430
## 25461 -5.047628         4.6589e-07   0.00201260
## 25604 -4.861572         1.2063e-06   0.00521110
## 25462 -4.679678         2.9619e-06   0.01279600

Heteroskedasticita

Heteroskedasticita (nekonštantný rozptyl rezíduí) spôsobuje, že klasické štandardné chyby a t-testy môžu byť nespoľahlivé. Budeme ju skúmať: 1. vizuálne – pomocou grafov štvorcov rezíduí, 2. formálne – Breusch–Pagan testom (bptest).

Vizuálna kontrola (štvorce rezíduí)

Ako potenciálne problematické premenné v tejto databáze zmysluplne preveríme napríklad teplotu a vlhkosť, prípadne aj vyrovnané hodnoty (fitted).

tmp <- udaje_model
tmp$sr2 <- resid(model2)^2
tmp$fitted <- fitted(model2)

p1 <- ggplot(tmp, aes(x = Temperature, y = sr2)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "Temperature", y = "Squared residuals", title = "Squared residuals vs Temperature")

p2 <- ggplot(tmp, aes(x = Humidity, y = sr2)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "Humidity", y = "Squared residuals", title = "Squared residuals vs Humidity")

p1 + p2

Vizuálna kontrola heteroskedasticity (model2)

ggplot(tmp, aes(x = fitted, y = sr2)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "Fitted values", y = "Squared residuals", title = "Squared residuals vs fitted values")

Squared residuals vs fitted values (model2)

Breusch–Pagan test

bptest(model1)

## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 1066.7, df = 34, p-value < 2.2e-16

bptest(model2)

## 
##  studentized Breusch-Pagan test
## 
## data:  model2
## BP = 1064.2, df = 34, p-value < 2.2e-16

Robustné (heteroskedasticity-consistent) štandardné chyby

Aj keď heteroskedasticita nie je výrazná, v praxi sa často reportujú robustné štandardné chyby (White/HC). Nižšie uvádzame odhady pre model2 s HC1 korekciou.

coeftest(model2, vcov. = vcovHC(model2, type = "HC1"))

## 
## t test of coefficients:
## 
##                                  Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)                    47442.3217    598.7108  79.2408 < 2.2e-16 ***
## Temperature                     -175.4384     24.1693  -7.2587 4.617e-13 ***
## Humidity                         -52.3578      2.9974 -17.4679 < 2.2e-16 ***
## WindSpeed                       -282.9507     22.4675 -12.5938 < 2.2e-16 ***
## I(log1p(GeneralDiffuseFlows))    488.6411     82.5838   5.9169 3.536e-09 ***
## I(log1p(DiffuseFlows))          -771.2077     66.4586 -11.6043 < 2.2e-16 ***
## factor(Hour)1                  -1318.2374    326.1968  -4.0412 5.410e-05 ***
## factor(Hour)2                   -872.5329    363.4400  -2.4008  0.016403 *  
## factor(Hour)3                  -6164.9508    330.7945 -18.6368 < 2.2e-16 ***
## factor(Hour)4                 -11968.9716    273.4080 -43.7770 < 2.2e-16 ***
## factor(Hour)5                 -16515.6105    239.4738 -68.9663 < 2.2e-16 ***
## factor(Hour)6                 -16439.4311    285.7244 -57.5360 < 2.2e-16 ***
## factor(Hour)7                 -14934.5411    372.4581 -40.0972 < 2.2e-16 ***
## factor(Hour)8                 -12797.7136    429.3244 -29.8090 < 2.2e-16 ***
## factor(Hour)9                 -10443.5591    453.2561 -23.0412 < 2.2e-16 ***
## factor(Hour)10                 -8394.6861    464.2972 -18.0804 < 2.2e-16 ***
## factor(Hour)11                 -6356.6240    467.1706 -13.6066 < 2.2e-16 ***
## factor(Hour)12                 -5197.9775    470.7731 -11.0414 < 2.2e-16 ***
## factor(Hour)13                 -4516.2052    484.7641  -9.3163 < 2.2e-16 ***
## factor(Hour)14                 -4052.5488    479.3877  -8.4536 < 2.2e-16 ***
## factor(Hour)15                 -4106.1628    470.6623  -8.7242 < 2.2e-16 ***
## factor(Hour)16                 -4350.3018    458.6032  -9.4860 < 2.2e-16 ***
## factor(Hour)17                 -3012.6718    448.9005  -6.7112 2.182e-11 ***
## factor(Hour)18                  -293.2530    457.3322  -0.6412  0.521411    
## factor(Hour)19                  4786.4585    411.4654  11.6327 < 2.2e-16 ***
## factor(Hour)20                  5986.1194    284.1544  21.0664 < 2.2e-16 ***
## factor(Hour)21                  5125.0494    263.6862  19.4362 < 2.2e-16 ***
## factor(Hour)22                  4904.8127    262.5919  18.6785 < 2.2e-16 ***
## factor(Hour)23                  2979.9242    285.9998  10.4193 < 2.2e-16 ***
## DayOfWeekTue                    -430.7847    170.4356  -2.5276  0.011522 *  
## DayOfWeekWed                     196.6000    158.8200   1.2379  0.215829    
## DayOfWeekThu                     125.3096    147.1084   0.8518  0.394363    
## DayOfWeekFri                     396.7987    150.0967   2.6436  0.008232 ** 
## DayOfWeekSat                    1331.9508    141.7264   9.3980 < 2.2e-16 ***
## DayOfWeekSun                    -771.5285    166.1685  -4.6430 3.536e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Záver

Model spotreby elektriny dáva zmysel, ak popri meteorologických premenných zohľadníme aj časové efekty (hodina dňa, deň v týždni). Transformácia log(1 + x) pri premenných žiarenia pomáha stabilizovať rozdelenie a znižovať vplyv extrémnych hodnôt. Pri interpretácii normality rezíduí treba brať do úvahy veľkosť vzorky a opierať sa aj o diagnostické grafy. Ak by sa objavila heteroskedasticita, vhodným riešením je reportovať robustné štandardné chyby.

Cvičenie 7 – Špecifikácia modelu – spotreba elektriny

## Balíky pre regresiu, testy a diagnostiku
library(lmtest)
library(sandwich)
library(car)

## Pomocné balíky na prácu s dátumom a úpravy dát
library(dplyr)
library(lubridate)

V tomto zošite skúmame, ako je spotreba elektriny (zvolená zóna) ovplyvnená meteorologickými premennými (teplota, vlhkosť, rýchlosť vetra) a premennými súvisiacimi so slnečným žiarením (GeneralDiffuseFlows, DiffuseFlows).

Budeme postupovať nasledovne: 1. odhad základného lineárneho modelu, 2. otestovanie funkčnej formy (RESET), 3. diagnostika (reziduá, C+R grafy), 4. návrhy nelineárnych úprav (kvadrát, prípadne zlom cez dummy premennú).

1. Príprava údajov

## Načítanie dát
## (ak máš súbor v inom priečinku, uprav si cestu)
udaje <- read.csv("powerconsumption.csv", sep = ",", dec = ".", header = TRUE)

## Prevod dátumu/času
udaje$Datetime <- mdy_hm(udaje$Datetime)

## Pomocné časové premenné (môžu byť užitočné v rozšíreniach modelu)
udaje <- udaje %>%
  mutate(
    Hour = hour(Datetime),
    Weekday = wday(Datetime, label = TRUE, week_start = 1)
  )

## Rýchla kontrola štruktúry
str(udaje)

## 'data.frame':    52416 obs. of  11 variables:
##  $ Datetime              : POSIXct, format: "2017-01-01 00:00:00" "2017-01-01 00:10:00" ...
##  $ 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 ...
##  $ Hour                  : int  0 0 0 0 0 0 1 1 1 1 ...
##  $ Weekday               : Ord.factor w/ 7 levels "Mon"<"Tue"<"Wed"<..: 7 7 7 7 7 7 7 7 7 7 ...

summary(udaje)

##     Datetime                    Temperature        Humidity       WindSpeed    
##  Min.   :2017-01-01 00:00:00   Min.   : 3.247   Min.   :11.34   Min.   :0.050  
##  1st Qu.:2017-04-01 23:57:30   1st Qu.:14.410   1st Qu.:58.31   1st Qu.:0.078  
##  Median :2017-07-01 23:55:00   Median :18.780   Median :69.86   Median :0.086  
##  Mean   :2017-07-01 23:55:00   Mean   :18.810   Mean   :68.26   Mean   :1.959  
##  3rd Qu.:2017-09-30 23:52:30   3rd Qu.:22.890   3rd Qu.:81.40   3rd Qu.:4.915  
##  Max.   :2017-12-30 23:50:00   Max.   :40.010   Max.   :94.80   Max.   :6.483  
##                                                                                
##  GeneralDiffuseFlows  DiffuseFlows     PowerConsumption_Zone1
##  Min.   :   0.004    Min.   :  0.011   Min.   :13896         
##  1st Qu.:   0.062    1st Qu.:  0.122   1st Qu.:26311         
##  Median :   5.035    Median :  4.456   Median :32266         
##  Mean   : 182.697    Mean   : 75.028   Mean   :32345         
##  3rd Qu.: 319.600    3rd Qu.:101.000   3rd Qu.:37309         
##  Max.   :1163.000    Max.   :936.000   Max.   :52204         
##                                                              
##  PowerConsumption_Zone2 PowerConsumption_Zone3      Hour       Weekday   
##  Min.   : 8560          Min.   : 5935          Min.   : 0.00   Mon:7488  
##  1st Qu.:16981          1st Qu.:13129          1st Qu.: 5.75   Tue:7488  
##  Median :20823          Median :16415          Median :11.50   Wed:7488  
##  Mean   :21043          Mean   :17835          Mean   :11.50   Thu:7488  
##  3rd Qu.:24714          3rd Qu.:21624          3rd Qu.:17.25   Fri:7488  
##  Max.   :37409          Max.   :47598          Max.   :23.00   Sat:7488  
##                                                                Sun:7488

Poznámka: V týchto dátach typicky nebývajú chýbajúce hodnoty; ak by sa objavili, riešenie doplnenia (imputácie) by bolo vhodné doplniť až podľa zadania.

2. Základná lineárna regresia

Ako základ použijeme lineárny model, kde spotrebu (Zone1) vysvetľujeme počasím a difúznym žiarením.

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

summary(model0)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15613.5  -4933.0   -653.3   4113.8  17982.7 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.675e+04  2.138e+02 125.132   <2e-16 ***
## Temperature          5.349e+02  6.418e+00  83.331   <2e-16 ***
## Humidity            -5.651e+01  2.130e+00 -26.525   <2e-16 ***
## WindSpeed           -1.487e+02  1.358e+01 -10.951   <2e-16 ***
## GeneralDiffuseFlows -1.702e+00  1.464e-01 -11.628   <2e-16 ***
## DiffuseFlows        -8.731e-02  2.721e-01  -0.321    0.748    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6349 on 52410 degrees of freedom
## Multiple R-squared:  0.2073, Adjusted R-squared:  0.2072 
## F-statistic:  2741 on 5 and 52410 DF,  p-value: < 2.2e-16

Ak chceš rovnaký model pre inú zónu, stačí zmeniť ľavú stranu napr. na PowerConsumption_Zone2 alebo PowerConsumption_Zone3.

3. Test správnej funkčnej formy – Ramsey RESET

RESET test kontroluje, či lineárna špecifikácia nie je „príliš jednoduchá“ (napr. chýba nelinearita alebo dôležité premenné).

Hypotézy: - H0: model je funkčne správne špecifikovaný, - H1: model je nesprávne špecifikovaný (typicky nelinearita / vynechané premenné).

resettest(model0)

## 
##  RESET test
## 
## data:  model0
## RESET = 92.704, df1 = 2, df2 = 52408, p-value < 2.2e-16

Interpretácia: - ak p-hodnota < 0.05, model je pravdepodobne funkčne zle špecifikovaný, - ak p-hodnota ≥ 0.05, nevidíme silný dôkaz proti lineárnej špecifikácii.

4. Grafická diagnostika

4.1 Residuals vs Fitted

plot(model0, which = 1)

Ak vidíš zakrivenie alebo systematický vzor, často to signalizuje nelinearitu (napr. potrebu log/kvadrátu/interakcie).

4.2 Component + Residual (C+R) grafy

C+R grafy pomôžu identifikovať, ktorá konkrétna premenná má nelineárny vzťah k vysvetľovanej premennej.

crPlots(model0)

5. Nelineárna špecifikácia (kvadratické členy)

Pri spotrebe elektriny býva nelinearita často pri teplote (napr. kúrenie/chladenie). Skúsime teda pridať kvadratický člen teploty.

model_quad <- lm(
  PowerConsumption_Zone1 ~ Temperature + I(Temperature^2) +
    Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows,
  data = udaje
)

summary(model_quad)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + I(Temperature^2) + 
##     Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows, 
##     data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15556.9  -4897.7   -658.7   4118.7  18073.0 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         25687.2407   278.2829  92.306  < 2e-16 ***
## Temperature           683.1004    25.6667  26.614  < 2e-16 ***
## I(Temperature^2)       -4.0173     0.6735  -5.965 2.46e-09 ***
## Humidity              -59.3702     2.1831 -27.196  < 2e-16 ***
## WindSpeed            -138.6930    13.6779 -10.140  < 2e-16 ***
## GeneralDiffuseFlows    -1.5517     0.1485 -10.451  < 2e-16 ***
## DiffuseFlows           -0.3505     0.2755  -1.272    0.203    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6347 on 52409 degrees of freedom
## Multiple R-squared:  0.2078, Adjusted R-squared:  0.2077 
## F-statistic:  2291 on 6 and 52409 DF,  p-value: < 2.2e-16

## Porovnanie so základným modelom
anova(model0, model_quad)

## RESET po rozšírení modelu
resettest(model_quad)

## 
##  RESET test
## 
## data:  model_quad
## RESET = 262.5, df1 = 2, df2 = 52407, p-value < 2.2e-16

6. Zlom cez dummy premennú (deň vs noc) + zmena sklonu

Vzťahy v spotrebe sa často líšia cez deň a v noci. Ako jednoduchý indikátor dňa/noci použijeme GeneralDiffuseFlows: - DAY = 1, ak GeneralDiffuseFlows > 0 - DAY = 0, inak

Najprv otestujeme posun (zlom v intercept-e), potom zmenu sklonu pri teplote.

udaje$DAY <- ifelse(udaje$GeneralDiffuseFlows > 0, 1, 0)

## (A) zlom v autonómnom člene (posun úrovne spotreby cez deň)
model_day_shift <- lm(
  PowerConsumption_Zone1 ~ DAY + Temperature + Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows,
  data = udaje
)
summary(model_day_shift)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ DAY + Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15613.5  -4933.0   -653.3   4113.8  17982.7 
## 
## Coefficients: (1 not defined because of singularities)
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.675e+04  2.138e+02 125.132   <2e-16 ***
## DAY                         NA         NA      NA       NA    
## Temperature          5.349e+02  6.418e+00  83.331   <2e-16 ***
## Humidity            -5.651e+01  2.130e+00 -26.525   <2e-16 ***
## WindSpeed           -1.487e+02  1.358e+01 -10.951   <2e-16 ***
## GeneralDiffuseFlows -1.702e+00  1.464e-01 -11.628   <2e-16 ***
## DiffuseFlows        -8.731e-02  2.721e-01  -0.321    0.748    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6349 on 52410 degrees of freedom
## Multiple R-squared:  0.2073, Adjusted R-squared:  0.2072 
## F-statistic:  2741 on 5 and 52410 DF,  p-value: < 2.2e-16

## (B) zlom v sklone: teplota môže pôsobiť inak cez deň a inak v noci
model_day_slope <- lm(
  PowerConsumption_Zone1 ~ Temperature + I(DAY * Temperature) +
    Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows,
  data = udaje
)
summary(model_day_slope)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + I(DAY * Temperature) + 
##     Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows, 
##     data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15613.5  -4933.0   -653.3   4113.8  17982.7 
## 
## Coefficients: (1 not defined because of singularities)
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           2.675e+04  2.138e+02 125.132   <2e-16 ***
## Temperature           5.349e+02  6.418e+00  83.331   <2e-16 ***
## I(DAY * Temperature)         NA         NA      NA       NA    
## Humidity             -5.651e+01  2.130e+00 -26.525   <2e-16 ***
## WindSpeed            -1.487e+02  1.358e+01 -10.951   <2e-16 ***
## GeneralDiffuseFlows  -1.702e+00  1.464e-01 -11.628   <2e-16 ***
## DiffuseFlows         -8.731e-02  2.721e-01  -0.321    0.748    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6349 on 52410 degrees of freedom
## Multiple R-squared:  0.2073, Adjusted R-squared:  0.2072 
## F-statistic:  2741 on 5 and 52410 DF,  p-value: < 2.2e-16

## Porovnanie so základom
anova(model0, model_day_slope)

## RESET test pre model so zlomom v sklone
resettest(model_day_slope)

## 
##  RESET test
## 
## data:  model_day_slope
## RESET = 92.703, df1 = 2, df2 = 52407, p-value < 2.2e-16

Interpretácia interakcie I(DAY * Temperature): - koeficient pri Temperature je efekt teploty v noci (DAY = 0), - koeficient pri I(DAY * Temperature) je rozdiel v efekte teploty cez deň oproti noci.

7. Box–Cox transformácia (voliteľné)

Ak by model trpel nelinearitou v závislej premennej, niekedy pomôže transformácia spotreby. Toto je voliteľné, keďže môže zhoršiť interpretáciu.

library(MASS)
boxcox(model0)

Ak by vhodná transformácia vyšla blízko logaritmu, dá sa skúsiť:

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

## 
## Call:
## lm(formula = log(PowerConsumption_Zone1) ~ Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65223 -0.15051 -0.00796  0.13554  0.51379 
## 
## Coefficients:
##                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)          1.018e+01  6.669e-03 1525.693  < 2e-16 ***
## Temperature          1.661e-02  2.002e-04   82.955  < 2e-16 ***
## Humidity            -1.723e-03  6.646e-05  -25.927  < 2e-16 ***
## WindSpeed           -4.780e-03  4.236e-04  -11.284  < 2e-16 ***
## GeneralDiffuseFlows -2.224e-05  4.566e-06   -4.871 1.11e-06 ***
## DiffuseFlows         3.983e-05  8.487e-06    4.693 2.70e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1981 on 52410 degrees of freedom
## Multiple R-squared:  0.2201, Adjusted R-squared:   0.22 
## F-statistic:  2958 on 5 and 52410 DF,  p-value: < 2.2e-16

resettest(model_log)

## 
##  RESET test
## 
## data:  model_log
## RESET = 269.93, df1 = 2, df2 = 52408, p-value < 2.2e-16

8. Robustné štandardné chyby (voliteľné)

Pri dátach s časovou štruktúrou býva užitočné pozrieť aj robustné (HC) štandardné chyby:

coeftest(model0, vcov = vcovHC(model0, type = "HC1"))

## 
## t test of coefficients:
## 
##                        Estimate  Std. Error  t value Pr(>|t|)    
## (Intercept)          2.6750e+04  2.2087e+02 121.1130   <2e-16 ***
## Temperature          5.3486e+02  6.3420e+00  84.3361   <2e-16 ***
## Humidity            -5.6510e+01  2.1607e+00 -26.1539   <2e-16 ***
## WindSpeed           -1.4870e+02  1.3402e+01 -11.0953   <2e-16 ***
## GeneralDiffuseFlows -1.7020e+00  1.1080e-01 -15.3611   <2e-16 ***
## DiffuseFlows        -8.7307e-02  1.7969e-01  -0.4859   0.6271    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cvičenie 8 – Zhluková analýza (hierarchické klastrovanie) – spotreba elektriny

Úvod

Cieľom je ukázať hierarchickú zhlukovú analýzu na dátach o spotrebe elektriny. Databáza obsahuje merania v 10-minútovom kroku (teplota, vlhkosť, vietor, difúzne žiarenie a spotreba v troch zónach).
Keďže ide o časový rad s veľkým počtom riadkov, zhlukovanie na úrovni jednotlivých 10-minútových meraní by bolo neprehľadné. Preto budeme zhlukovať dni: každý deň opíšeme priemernými hodnotami premenných a následne zhlukujeme podobné dni.

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

Načítanie a príprava dát

## 1) Načítanie dát
data_raw <- read.csv("powerconsumption.csv", stringsAsFactors = FALSE)

## 2) Prevod dátumu a času
## Formát v súbore vyzerá napr. "1/1/2017 0:00"
data_raw$Datetime <- as.POSIXct(data_raw$Datetime, format = "%m/%d/%Y %H:%M")

## 3) Kontrola základných informácií
str(data_raw)

## 'data.frame':    52416 obs. of  9 variables:
##  $ Datetime              : POSIXct, format: "2017-01-01 00:00:00" "2017-01-01 00:10:00" ...
##  $ 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 ...

Pre ďalšie kroky vytvoríme denný dataset (priemery za deň). Spotrebu spojíme aj do jednej premennej TotalPower (súčet zón).

data_raw <- data_raw %>%
  mutate(
    Date = as.Date(Datetime),
    TotalPower = PowerConsumption_Zone1 + PowerConsumption_Zone2 + PowerConsumption_Zone3
  )

## Denné priemery
data_day <- data_raw %>%
  group_by(Date) %>%
  summarise(
    across(
      .cols = c(Temperature, Humidity, WindSpeed, GeneralDiffuseFlows, DiffuseFlows,
                PowerConsumption_Zone1, PowerConsumption_Zone2, PowerConsumption_Zone3, TotalPower),
      .fns = ~mean(.x, na.rm = TRUE)
    ),
    .groups = "drop"
  )

## Nastavenie riadkových mien pre prehľadnejšie výstupy
data_day_df <- as.data.frame(data_day)
data_day_df <- data_day_df %>% 
  filter(!is.na(Date))

## Rýchla ukážka (len prvých 10 dní)
kable(head(data_day, 10), caption = "Tab. 1: Ukážka denných priemerov (prvých 10 dní)") %>%
  kable_styling(full_width = FALSE)

Tab. 1: Ukážka denných priemerov (prvých 10 dní)
Date	Temperature	Humidity	WindSpeed	GeneralDiffuseFlows	DiffuseFlows	PowerConsumption_Zone1	PowerConsumption_Zone2	PowerConsumption_Zone3	TotalPower
2017-01-01	9.675299	68.51931	0.3151458	121.39077	25.99392	28465.23	17737.79	17868.80	64071.82
2017-01-02	12.476875	71.45632	0.0765625	120.40449	27.22741	28869.49	19557.73	17820.76	66247.98
2017-01-03	12.100000	74.98167	0.0767153	120.68601	28.57466	30562.45	20057.27	17620.80	68240.52
2017-01-04	10.509479	75.45979	0.0824167	122.95932	28.82722	30689.83	20102.08	17673.69	68465.60
2017-01-05	10.866444	71.04049	0.0838958	118.74986	29.74144	30802.91	20033.94	17664.18	68501.03
2017-01-06	11.685208	77.25778	0.1459514	63.70045	59.33356	30712.24	20139.77	17481.08	68333.09
2017-01-07	12.018819	76.21708	0.0829306	88.99876	48.84979	30612.70	19252.43	17985.94	67851.08
2017-01-08	14.724653	69.10924	0.0813056	71.08221	46.38776	29559.92	17221.12	17810.60	64591.64
2017-01-09	15.490139	68.33236	0.0741736	110.43367	55.82734	30540.25	19102.63	17417.79	67060.68
2017-01-10	14.019583	75.67840	0.0782917	112.49842	47.99520	31863.29	19405.27	17425.34	68693.90

Škálovanie premenných

Pri zhlukovaní používame vzdialenosti. Aby premenné s veľkými jednotkami (napr. spotreba) nedomninovali všetko ostatné, použijeme z-škálovanie:

\[ z = \frac{x-\mu}{\sigma} \]

data_complete <- data_day_df %>%
  dplyr::select(-Date) %>%   # výslovne dplyr::select
  na.omit()

data_scaled <- scale(data_complete)

Boxploty škálovaných premenných

num_vars <- as.data.frame(data_scaled)
num_plots <- ncol(num_vars)

par(mfrow = c(ceiling(sqrt(num_plots)), ceiling(num_plots / ceiling(sqrt(num_plots)))))
par(mar = c(4, 4, 2, 1))

for (col in names(num_vars)) {
  boxplot(
    num_vars[[col]],
    main = col,
    col = "lightblue",
    horizontal = TRUE
  )
}

mtext("Boxploty škálovaných denných premenných", outer = TRUE, cex = 1.2, font = 2)

Korelácie premenných

Ak by boli niektoré premenné extrémne korelované, jedna z dvojice by zbytočne „zdvojovala“ informáciu. Pozrieme sa na korelačnú maticu.

cor_mat <- cor(data_scaled, use = "pairwise.complete.obs")
cor_mat <- round(cor_mat, 2)

kable(cor_mat, caption = "Tab. 2: Korelačná matica (škálované denné dáta)") %>%
  kable_styling(full_width = FALSE)

Tab. 2: Korelačná matica (škálované denné dáta)
	Temperature	Humidity	WindSpeed	GeneralDiffuseFlows	DiffuseFlows	PowerConsumption_Zone1	PowerConsumption_Zone2	PowerConsumption_Zone3	TotalPower
Temperature	1.00	-0.31	0.56	0.64	0.08	0.73	0.42	0.61	0.73
Humidity	-0.31	1.00	-0.17	-0.39	-0.04	-0.21	-0.21	-0.26	-0.29
WindSpeed	0.56	-0.17	1.00	0.37	-0.08	0.47	0.30	0.41	0.49
GeneralDiffuseFlows	0.64	-0.39	0.37	1.00	0.35	0.54	0.11	0.58	0.55
DiffuseFlows	0.08	-0.04	-0.08	0.35	1.00	0.11	-0.25	0.22	0.08
PowerConsumption_Zone1	0.73	-0.21	0.47	0.54	0.11	1.00	0.39	0.71	0.87
PowerConsumption_Zone2	0.42	-0.21	0.30	0.11	-0.25	0.39	1.00	0.23	0.59
PowerConsumption_Zone3	0.61	-0.26	0.41	0.58	0.22	0.71	0.23	1.00	0.90
TotalPower	0.73	-0.29	0.49	0.55	0.08	0.87	0.59	0.90	1.00

Matica vzdialeností (ukážka)

Vzdialenosť dvoch dní počítame euklidovsky zo všetkých použitých premenných. Celá matica je veľká (máme veľa dní), preto ukážeme len prvých 12 × 12.

dist_obj <- dist(data_scaled, method = "euclidean")

## Ukážka prvých 12 dní
dist_mat_small <- as.matrix(dist_obj)[1:12, 1:12]
dist_mat_small <- round(dist_mat_small, 2)

kable(dist_mat_small, caption = "Tab. 3: Ukážka matice vzdialeností (prvých 12 dní)") %>%
  kable_styling(full_width = FALSE)

Tab. 3: Ukážka matice vzdialeností (prvých 12 dní)
1	2	3	4	5	6	7	8	9	10	11	12
0.00	0.95	1.47	1.46	1.38	1.92	1.52	1.34	1.69	1.94	1.54	2.30
0.95	0.00	0.77	0.92	0.86	1.41	1.05	1.29	1.19	1.37	1.25	1.76
1.47	0.77	0.00	0.32	0.43	1.06	0.72	1.58	1.19	0.84	1.17	1.29
1.46	0.92	0.32	0.00	0.40	1.08	0.79	1.75	1.41	0.99	1.30	1.39
1.38	0.86	0.43	0.40	0.00	1.15	0.84	1.61	1.21	0.99	1.05	1.37
1.92	1.41	1.06	1.08	1.15	0.00	0.54	1.59	1.27	0.94	1.21	0.87
1.52	1.05	0.72	0.79	0.84	0.54	0.00	1.25	1.02	0.69	0.94	0.92
1.34	1.29	1.58	1.75	1.61	1.59	1.25	0.00	1.00	1.49	1.03	1.74
1.69	1.19	1.19	1.41	1.21	1.27	1.02	1.00	0.00	0.91	0.51	1.07
1.94	1.37	0.84	0.99	0.99	0.94	0.69	1.49	0.91	0.00	0.92	0.62
1.54	1.25	1.17	1.30	1.05	1.21	0.94	1.03	0.51	0.92	0.00	1.07
2.30	1.76	1.29	1.39	1.37	0.87	0.92	1.74	1.07	0.62	1.07	0.00

Hierarchické zhlukovanie (Wardova metóda)

Wardova metóda je aglomeratívna: začína s jednočlennými zhlukmi a spája ich tak, aby bol nárast vnútornej variability čo najmenší. Prakticky často vytvára kompaktné a interpretovateľné zhluky.

hc <- hclust(dist_obj, method = "ward.D2")

plot(
  hc,
  labels = FALSE,
  main = "Hierarchické zhlukovanie dní (Ward.D2)",
  xlab = "Dni (bez popisov pre čitateľnosť)",
  sub = ""
)

k <- 3
rect.hclust(hc, k = k, border = "red")

Príslušnosť dní do klastrov

cluster_membership <- cutree(hc, k = k)

clusters_df <- data.frame(
  Date = rownames(data_complete),
  Cluster = factor(cluster_membership)
)

kable(head(clusters_df, 20), caption = "Tab. 4: Prvých 20 dní a ich klaster") %>%
  kable_styling(full_width = FALSE)

Tab. 4: Prvých 20 dní a ich klaster
Date	Cluster
1	1
2	1
3	1
4	1
5	1
6	1
7	1
8	1
9	1
10	1
11	1
12	1
13	1
14	1
15	1
16	1
17	1
18	1
19	1
20	1

Centroidy (priemery premenných v klastroch)

Tu už pracujeme s pôvodnými (neškálovanými) dennými priemermi, aby boli hodnoty interpretovateľné.

data_day_for_desc <- as.data.frame(data_complete)
data_day_for_desc$Cluster <- factor(cluster_membership)

centroids <- data_day_for_desc %>%
  group_by(Cluster) %>%
  summarise(across(where(is.numeric), ~mean(.x, na.rm = TRUE)))

kable(centroids, digits = 2, caption = "Tab. 5: Centroidy – priemerné hodnoty premenných v klastroch") %>%
  kable_styling(full_width = FALSE)

Tab. 5: Centroidy – priemerné hodnoty premenných v klastroch
Cluster	Temperature	Humidity	WindSpeed	GeneralDiffuseFlows	DiffuseFlows	PowerConsumption_Zone1	PowerConsumption_Zone2	PowerConsumption_Zone3	TotalPower
1	14.19	71.19	0.76	113.64	65.35	30272.39	20701.00	15304.36	66277.75
2	20.79	67.95	2.17	225.55	87.75	32957.74	19918.76	16876.61	69753.11
3	26.31	61.40	4.59	265.71	71.84	36322.49	24400.04	26459.26	87181.79

Záver

Zhluková analýza rozdelila dni do niekoľkých typov podľa kombinácie počasia (teplota, vlhkosť, vietor, difúzne žiarenie) a spotreby elektriny (zónovo aj spolu). Takéto zhluky sa dajú interpretovať napríklad ako „chladné dni s vyššou spotrebou“, „teplé dni s nižšou spotrebou“, prípadne dni so špecifickým priebehom žiarenia a odlišnou spotrebou. V praxi môže byť výsledok využiteľný na segmentáciu dní pre plánovanie výroby, nákupu energie alebo detekciu atypických období.

Cvičenie 9 – Ekonometria v R – cvičenie 9 (autokorelácia rezíduí)

1. Úvod a údaje

V tomto cvičení pracujeme s databázou powerconsumption.csv, ktorá obsahuje merania po 10 minútach. Cieľom je odhadnúť lineárny model spotreby elektriny a následne overiť predpoklad nezávislosti rezíduí (autokoreláciu).

Budeme modelovať spotrebu v zóne 1 (PowerConsumption_Zone1) v závislosti od: - teploty (Temperature), - vlhkosti (Humidity), - rýchlosti vetra (WindSpeed), - difúznych tokov (GeneralDiffuseFlows, DiffuseFlows).

1.1 Pracovná hypotéza

Predpokladáme, že meteorologické podmienky a intenzita difúznych tokov súvisia so spotrebou elektriny. Znamenko vplyvu (kladné/záporné) nefixujeme „nasilu“, pretože spotreba môže reagovať rozdielne v zime a v lete (kúrenie vs. chladenie). Dôležitá je pre nás najmä štatistická významnosť a autokorelácia rezíduí v časovom rade.

1.2 Načítanie a príprava dát

library(dplyr)
library(lubridate)
library(ggplot2)
library(lmtest)
library(sandwich)

## nacitanie dat (desatinna bodka, oddelovac ciarka)
## predpoklad: subor je v priecinku "udaje/" pod nazvom "powerconsumption.csv"
udaje <- read.csv("powerconsumption.csv", sep = ",", dec=".", header = TRUE)

## parsovanie casu
udaje <- udaje %>%
  mutate(
    Datetime = mdy_hm(Datetime)
  ) %>%
  arrange(Datetime)

## kontrola
str(udaje)

## 'data.frame':    52416 obs. of  9 variables:
##  $ Datetime              : POSIXct, format: "2017-01-01 00:00:00" "2017-01-01 00:10:00" ...
##  $ 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 ...

head(udaje)

Poznámka: dáta sú po 10 minútach, takže pri priamom OLS modeli bude autokorelácia takmer istá. Pre prehľadnosť si údaje zhrnieme na hodinové priemery (časový rad je stále dostatočne dlhý, ale čitateľnejší).

udaje_h <- udaje %>%
  mutate(Datetime_h = floor_date(Datetime, unit = "hour")) %>%
  group_by(Datetime_h) %>%
  summarise(
    Temperature = mean(Temperature),
    Humidity = mean(Humidity),
    WindSpeed = mean(WindSpeed),
    GeneralDiffuseFlows = mean(GeneralDiffuseFlows),
    DiffuseFlows = mean(DiffuseFlows),
    PC_Zone1 = mean(PowerConsumption_Zone1),
    PC_Zone2 = mean(PowerConsumption_Zone2),
    PC_Zone3 = mean(PowerConsumption_Zone3),
    .groups = "drop"
  )

summary(udaje_h)

##    Datetime_h                   Temperature        Humidity    
##  Min.   :2017-01-01 00:00:00   Min.   : 3.602   Min.   :12.71  
##  1st Qu.:2017-04-01 23:45:00   1st Qu.:14.404   1st Qu.:58.32  
##  Median :2017-07-01 23:30:00   Median :18.759   Median :69.82  
##  Mean   :2017-07-01 23:30:00   Mean   :18.810   Mean   :68.26  
##  3rd Qu.:2017-09-30 23:15:00   3rd Qu.:22.867   3rd Qu.:81.35  
##  Max.   :2017-12-30 23:00:00   Max.   :39.695   Max.   :94.75  
##    WindSpeed       GeneralDiffuseFlows  DiffuseFlows         PC_Zone1    
##  Min.   :0.05467   Min.   :  0.019     Min.   :  0.0400   Min.   :14329  
##  1st Qu.:0.07817   1st Qu.:  0.064     1st Qu.:  0.1242   1st Qu.:26293  
##  Median :0.08550   Median :  9.947     Median :  8.2413   Median :32342  
##  Mean   :1.95949   Mean   :182.697     Mean   : 75.0280   Mean   :32345  
##  3rd Qu.:4.91533   3rd Qu.:326.488     3rd Qu.:105.8833   3rd Qu.:37318  
##  Max.   :5.93367   Max.   :953.350     Max.   :861.0000   Max.   :51844  
##     PC_Zone2        PC_Zone3    
##  Min.   : 8686   Min.   : 6191  
##  1st Qu.:17017   1st Qu.:13148  
##  Median :20787   Median :16428  
##  Mean   :21043   Mean   :17835  
##  3rd Qu.:24678   3rd Qu.:21598  
##  Max.   :36255   Max.   :47224

2. Lineárna regresia v základnom tvare

Odhadujeme model (pre zónu 1):

\[ PC\_Zone1_t = \beta_0 + \beta_1 Temperature_t + \beta_2 Humidity_t + \beta_3 WindSpeed_t + \beta_4 GeneralDiffuseFlows_t + \beta_5 DiffuseFlows_t + e_t \]

model <- lm(PC_Zone1 ~ Temperature + Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows,
            data = udaje_h)

summary(model)

## 
## Call:
## lm(formula = PC_Zone1 ~ Temperature + Humidity + WindSpeed + 
##     GeneralDiffuseFlows + DiffuseFlows, data = udaje_h)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14359.7  -4907.0   -581.5   4124.6  17464.5 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         26818.9010   519.9802  51.577  < 2e-16 ***
## Temperature           535.9826    15.6051  34.347  < 2e-16 ***
## Humidity              -57.7531     5.1953 -11.116  < 2e-16 ***
## WindSpeed            -149.8401    33.0769  -4.530 5.98e-06 ***
## GeneralDiffuseFlows    -1.8396     0.3638  -5.057 4.34e-07 ***
## DiffuseFlows            0.2148     0.6912   0.311    0.756    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6280 on 8730 degrees of freedom
## Multiple R-squared:  0.2111, Adjusted R-squared:  0.2107 
## F-statistic: 467.2 on 5 and 8730 DF,  p-value: < 2.2e-16

Graf: empíria vs. fitted

udaje_h$fitted <- fitted(model)

ggplot(udaje_h, aes(x = Datetime_h, y = PC_Zone1)) +
  geom_line() +
  geom_line(aes(y = fitted), linewidth = 0.8) +
  labs(
    title = "Spotreba elektriny (Zone1): skutočné hodnoty vs. vyrovnané hodnoty",
    x = "Čas (hodina)",
    y = "Power consumption (Zone1)"
  ) +
  theme_minimal()

3. Autokorelácia rezíduí

V časových radoch často platí, že rezíduá nie sú nezávislé. Ak je chyba v čase \(t\) prepojená s chybou v čase \(t-1\), hovoríme o autokorelácii.

3.1 Rezíduá a ACF

res <- residuals(model)

acf(res, lag.max = 48, main = "ACF rezíduí (do 48 hodín)")

3.2 Durbin–Watsonov test

DW test je klasická kontrola autokorelácie 1. rádu.

dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.19372, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

3.3 Breusch–Godfrey test

BG test umožňuje testovať autokoreláciu aj pre viac lagov. Pri hodinových dátach dáva zmysel skúsiť napr. 24 hodín (1 deň).

bgtest(model, order = 24)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  model
## LM test = 8283.8, df = 24, p-value < 2.2e-16

4. Ako riešiť autokoreláciu (prakticky)

Ukážeme dve bežné možnosti: 1) dynamický model s oneskorenou závislou premennou, 2) robustné štandardné chyby (Newey–West).

4.1 Dynamizácia: lag závislej premennej (Koyckov typ uvažovania)

Pridáme \(PC\_Zone1_{t-1}\) do regresie. Tým zachytíme zotrvačnosť spotreby.

udaje_h <- udaje_h %>%
  arrange(Datetime_h) %>%
  mutate(PC_Zone1_lag1 = lag(PC_Zone1))

model_dyn <- lm(PC_Zone1 ~ Temperature + Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows + PC_Zone1_lag1,
                data = udaje_h)

summary(model_dyn)

## 
## Call:
## lm(formula = PC_Zone1 ~ Temperature + Humidity + WindSpeed + 
##     GeneralDiffuseFlows + DiffuseFlows + PC_Zone1_lag1, data = udaje_h)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10137.3  -1555.1   -320.4   1366.6  10766.9 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.927e+03  2.458e+02  11.907  < 2e-16 ***
## Temperature          4.413e+01  6.923e+00   6.373 1.94e-10 ***
## Humidity            -1.301e+01  2.172e+00  -5.991 2.17e-09 ***
## WindSpeed            5.274e+00  1.378e+01   0.383    0.702    
## GeneralDiffuseFlows  8.011e-01  1.518e-01   5.276 1.35e-07 ***
## DiffuseFlows         3.796e+00  2.880e-01  13.182  < 2e-16 ***
## PC_Zone1_lag1        8.976e-01  4.393e-03 204.335  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2612 on 8728 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.8636, Adjusted R-squared:  0.8635 
## F-statistic:  9210 on 6 and 8728 DF,  p-value: < 2.2e-16

dwtest(model_dyn)

## 
##  Durbin-Watson test
## 
## data:  model_dyn
## DW = 0.82121, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

bgtest(model_dyn, order = 24)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  model_dyn
## LM test = 7770.8, df = 24, p-value < 2.2e-16

Poznámka: pri dynamickom modeli sa často zmení významnosť „počasných“ premenných, lebo časť variability vysvetlí samotná zotrvačnosť spotreby.

4.2 Newey–West robustné štandardné chyby

Ak chceme ponechať pôvodný OLS odhad koeficientov, ale opraviť štandardné chyby pre autokoreláciu a heteroskedasticitu, použijeme Newey–West.

Pri hodinových dátach je bežná voľba napr. lag 24 (1 deň).

coeftest(model, vcov = NeweyWest(model, lag = 24, prewhite = FALSE))

## 
## t test of coefficients:
## 
##                        Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)         26818.90096   885.12630 30.2995 < 2.2e-16 ***
## Temperature           535.98259    25.81092 20.7657 < 2.2e-16 ***
## Humidity              -57.75314     9.35556 -6.1731 6.993e-10 ***
## WindSpeed            -149.84009    56.57687 -2.6484  0.008101 ** 
## GeneralDiffuseFlows    -1.83957     0.40275 -4.5675 5.003e-06 ***
## DiffuseFlows            0.21485     0.79960  0.2687  0.788169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5. Záver

Na dátach spotreby elektriny (časový rad) je autokorelácia rezíduí prirodzená a očakávaná. Formálne testy (DW, BG) a ACF graf zvyčajne potvrdia, že predpoklad nezávislosti rezíduí v základnom OLS modeli nie je splnený.

Riešenia: - dynamický model (zachytí zotrvačnosť), - alebo robustné odhady štandardných chýb (Newey–West), ak nám ide hlavne o korektnú inferenciu.

Cvičenie 10 – Cvičenie 10 – Multikolinearita v regresných modeloch (Power Consumption)

1. Čo je multikolinearita

Multikolinearita znamená, že niektoré vysvetľujúce premenné v regresii sú medzi sebou silne (takmer lineárne) závislé. Odhady koeficientov metódou najmenších štvorcov zostávajú neskreslené, ale:

rastú štandardné chyby koeficientov (koeficienty “nevychádzajú významné”),
koeficienty sú nestabilné (menia sa pri malej zmene dát),
interpretácia jednotlivých vplyvov je nespoľahlivá.

2. Dáta a východiskový model

Budeme pracovať s dátami o spotrebe elektriny a počasí. Databáza obsahuje 10-minútové merania a tri zóny spotreby.

Budeme odhadovať model:

\[ PowerConsumption\_Zone1_t = \beta_0 + \beta_1 Temperature_t + \beta_2 Humidity_t + \beta_3 WindSpeed_t + \beta_4 GeneralDiffuseFlows_t + \beta_5 DiffuseFlows_t + \beta_6 PowerConsumption\_Zone2_t + u_t \]

Premennú PowerConsumption_Zone2 pridávame zámerne – je silno korelovaná so Zone1, takže vie vytvoriť výraznú multikolinearitu.

Predpoklad: súbor powerconsumption.csv je v rovnakom priečinku ako tento .Rmd. Ak nie, uprav cestu v read.csv().

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

## Datetime si pre istotu konvertujeme na časový typ (nie je nutné pre regresiu, ale hodí sa pri kontrole)
udaje$Datetime <- as.POSIXct(udaje$Datetime, format = "%m/%d/%Y %H:%M")

str(udaje)

## 'data.frame':    52416 obs. of  9 variables:
##  $ Datetime              : POSIXct, format: "2017-01-01 00:00:00" "2017-01-01 00:10:00" ...
##  $ 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 ...

summary(udaje)

##     Datetime                    Temperature        Humidity       WindSpeed    
##  Min.   :2017-01-01 00:00:00   Min.   : 3.247   Min.   :11.34   Min.   :0.050  
##  1st Qu.:2017-04-01 23:57:30   1st Qu.:14.410   1st Qu.:58.31   1st Qu.:0.078  
##  Median :2017-07-01 23:55:00   Median :18.780   Median :69.86   Median :0.086  
##  Mean   :2017-07-01 23:55:00   Mean   :18.810   Mean   :68.26   Mean   :1.959  
##  3rd Qu.:2017-09-30 23:52:30   3rd Qu.:22.890   3rd Qu.:81.40   3rd Qu.:4.915  
##  Max.   :2017-12-30 23:50:00   Max.   :40.010   Max.   :94.80   Max.   :6.483  
##  GeneralDiffuseFlows  DiffuseFlows     PowerConsumption_Zone1
##  Min.   :   0.004    Min.   :  0.011   Min.   :13896         
##  1st Qu.:   0.062    1st Qu.:  0.122   1st Qu.:26311         
##  Median :   5.035    Median :  4.456   Median :32266         
##  Mean   : 182.697    Mean   : 75.028   Mean   :32345         
##  3rd Qu.: 319.600    3rd Qu.:101.000   3rd Qu.:37309         
##  Max.   :1163.000    Max.   :936.000   Max.   :52204         
##  PowerConsumption_Zone2 PowerConsumption_Zone3
##  Min.   : 8560          Min.   : 5935         
##  1st Qu.:16981          1st Qu.:13129         
##  Median :20823          Median :16415         
##  Mean   :21043          Mean   :17835         
##  3rd Qu.:24714          3rd Qu.:21624         
##  Max.   :37409          Max.   :47598

3. Odhad základného regresného modelu

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

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows + PowerConsumption_Zone2, 
##     data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17249.8  -2493.9     71.9   2444.0  15359.2 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             5.803e+03  1.455e+02  39.875  < 2e-16 ***
## Temperature             1.926e+02  4.013e+00  48.002  < 2e-16 ***
## Humidity                6.218e+00  1.295e+00   4.801 1.58e-06 ***
## WindSpeed              -5.599e+01  8.156e+00  -6.865 6.72e-12 ***
## GeneralDiffuseFlows    -3.408e-01  8.797e-02  -3.875 0.000107 ***
## DiffuseFlows            1.439e+00  1.634e-01   8.811  < 2e-16 ***
## PowerConsumption_Zone2  1.072e+00  3.514e-03 305.066  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3811 on 52409 degrees of freedom
## Multiple R-squared:  0.7144, Adjusted R-squared:  0.7144 
## F-statistic: 2.185e+04 on 6 and 52409 DF,  p-value: < 2.2e-16

4. Korelačná matica a scatterploty

Korelácie zachytávajú párové vzťahy. Pri multikolinearite nás zaujíma hlavne korelácia medzi vysvetľujúcimi premennými.

xvars <- udaje[, c("Temperature", "Humidity", "WindSpeed",
                   "GeneralDiffuseFlows", "DiffuseFlows",
                   "PowerConsumption_Zone2")]

round(cor(xvars), 3)

##                        Temperature Humidity WindSpeed GeneralDiffuseFlows
## Temperature                  1.000   -0.460     0.477               0.460
## Humidity                    -0.460    1.000    -0.136              -0.468
## WindSpeed                    0.477   -0.136     1.000               0.134
## GeneralDiffuseFlows          0.460   -0.468     0.134               1.000
## DiffuseFlows                 0.197   -0.257    -0.001               0.565
## PowerConsumption_Zone2       0.382   -0.295     0.146               0.157
##                        DiffuseFlows PowerConsumption_Zone2
## Temperature                   0.197                  0.382
## Humidity                     -0.257                 -0.295
## WindSpeed                    -0.001                  0.146
## GeneralDiffuseFlows           0.565                  0.157
## DiffuseFlows                  1.000                  0.045
## PowerConsumption_Zone2        0.045                  1.000

pairs(xvars, main = "Scatterplotová matica – vysvetľujúce premenné")

5. VIF (Variance Inflation Factor)

VIF pre premennú \(x_j\):

\[ VIF_j = \frac{1}{1 - R_j^2} \]

Pravidlá (orientačne): - VIF > 5: problém môže byť relevantný, - VIF > 10: silná multikolinearita.

library(car)
vif(model)

##            Temperature               Humidity              WindSpeed 
##               1.965281               1.464191               1.324643 
##    GeneralDiffuseFlows           DiffuseFlows PowerConsumption_Zone2 
##               1.952539               1.486298               1.205864

6. Condition Number (podmienenosť)

Condition number je indikátor, ktorý vie zachytiť aj “cyklickú” multikolinearitu medzi viacerými premennými naraz. Pravidlo (orientačne): - < 10 nízka, - 10–30 mierna, - 30–100 silná, - > 100 veľmi vážna.

Dôležité: Condition number je citlivý na mierku premenných (jednotky), preto pri porovnávaní často pomáha škálovanie.

X <- model.matrix(model)[, -1]   # bez interceptu
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

## [1] 10756.87

7. Riešenia multikolinearity

7.1 Vynechanie problematickej premennej

Najčastejší “čistý” zásah je vynechať regresor, ktorý spôsobuje multikolinearitu. Tu skúsime vynechať PowerConsumption_Zone2 a porovnať výsledok.

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

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15613.5  -4933.0   -653.3   4113.8  17982.7 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.675e+04  2.138e+02 125.132   <2e-16 ***
## Temperature          5.349e+02  6.418e+00  83.331   <2e-16 ***
## Humidity            -5.651e+01  2.130e+00 -26.525   <2e-16 ***
## WindSpeed           -1.487e+02  1.358e+01 -10.951   <2e-16 ***
## GeneralDiffuseFlows -1.702e+00  1.464e-01 -11.628   <2e-16 ***
## DiffuseFlows        -8.731e-02  2.721e-01  -0.321    0.748    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6349 on 52410 degrees of freedom
## Multiple R-squared:  0.2073, Adjusted R-squared:  0.2072 
## F-statistic:  2741 on 5 and 52410 DF,  p-value: < 2.2e-16

vif(model_noZone2)

##         Temperature            Humidity           WindSpeed GeneralDiffuseFlows 
##            1.811650            1.427287            1.322804            1.947516 
##        DiffuseFlows 
##            1.484903

7.2 Škálovanie (štandardizácia) vysvetľujúcich premenných

Škálovanie nemení “informatívnosť” dát, ale výrazne pomôže pri numerickej stabilite a pri Condition number. Koeficienty sa potom interpretujú v štandardizovaných jednotkách (menej intuitívne v pôvodných jednotkách).

udaje_scaled <- udaje

cols_to_scale <- c("Temperature", "Humidity", "WindSpeed",
                   "GeneralDiffuseFlows", "DiffuseFlows",
                   "PowerConsumption_Zone2")

udaje_scaled[cols_to_scale] <- scale(udaje_scaled[cols_to_scale], center = TRUE, scale = TRUE)

model_scaled <- lm(
  PowerConsumption_Zone1 ~ Temperature + Humidity + WindSpeed +
    GeneralDiffuseFlows + DiffuseFlows + PowerConsumption_Zone2,
  data = udaje_scaled
)

summary(model_scaled)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows + PowerConsumption_Zone2, 
##     data = udaje_scaled)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17249.8  -2493.9     71.9   2444.0  15359.2 
## 
## Coefficients:
##                        Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)            32344.97      16.65 1943.198  < 2e-16 ***
## Temperature             1120.13      23.33   48.002  < 2e-16 ***
## Humidity                  96.70      20.14    4.801 1.58e-06 ***
## WindSpeed               -131.52      19.16   -6.865 6.72e-12 ***
## GeneralDiffuseFlows      -90.12      23.26   -3.875 0.000107 ***
## DiffuseFlows             178.80      20.29    8.811  < 2e-16 ***
## PowerConsumption_Zone2  5576.18      18.28  305.066  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3811 on 52409 degrees of freedom
## Multiple R-squared:  0.7144, Adjusted R-squared:  0.7144 
## F-statistic: 2.185e+04 on 6 and 52409 DF,  p-value: < 2.2e-16

vif(model_scaled)

##            Temperature               Humidity              WindSpeed 
##               1.965281               1.464191               1.324643 
##    GeneralDiffuseFlows           DiffuseFlows PowerConsumption_Zone2 
##               1.952539               1.486298               1.205864

X <- model.matrix(model_scaled)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)
condition_number_scaled <- sqrt(max(eig$values) / min(eig$values))
condition_number_scaled

## [1] 2.755976

7.3 Jednoduchá zmena jednotiek (zachová interpretáciu)

Ak je problém hlavne v mierke (jednotkách), dá sa prepočítať len konkrétna premenná. Tu napríklad prepočítame spotrebu v Zone2 na “tisíce jednotiek” (len príklad – reálne jednotky závisia od datasetu).

udaje_unit <- udaje
udaje_unit$Zone2_1000 <- udaje_unit$PowerConsumption_Zone2 / 1000

model_unit <- lm(
  PowerConsumption_Zone1 ~ Temperature + Humidity + WindSpeed +
    GeneralDiffuseFlows + DiffuseFlows + Zone2_1000,
  data = udaje_unit
)

summary(model_unit)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + GeneralDiffuseFlows + DiffuseFlows + Zone2_1000, 
##     data = udaje_unit)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17249.8  -2493.9     71.9   2444.0  15359.2 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5803.08166  145.53221  39.875  < 2e-16 ***
## Temperature          192.61151    4.01255  48.002  < 2e-16 ***
## Humidity               6.21787    1.29518   4.801 1.58e-06 ***
## WindSpeed            -55.99210    8.15616  -6.865 6.72e-12 ***
## GeneralDiffuseFlows   -0.34084    0.08797  -3.875 0.000107 ***
## DiffuseFlows           1.43949    0.16338   8.811  < 2e-16 ***
## Zone2_1000          1072.04024    3.51412 305.066  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3811 on 52409 degrees of freedom
## Multiple R-squared:  0.7144, Adjusted R-squared:  0.7144 
## F-statistic: 2.185e+04 on 6 and 52409 DF,  p-value: < 2.2e-16

vif(model_unit)

##         Temperature            Humidity           WindSpeed GeneralDiffuseFlows 
##            1.965281            1.464191            1.324643            1.952539 
##        DiffuseFlows          Zone2_1000 
##            1.486298            1.205864

X <- model.matrix(model_unit)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)
condition_number_unit <- sqrt(max(eig$values) / min(eig$values))
condition_number_unit

## [1] 169.4135

7.4 PCA (voliteľné)

PCA vytvorí komponent (kombináciu) z korelovaných premenných a tým vie znížiť multikolinearitu. Nižšie spravíme PCA z dvojice GeneralDiffuseFlows a DiffuseFlows a použijeme prvú hlavnú komponentu.

X_pca <- scale(udaje[, c("GeneralDiffuseFlows", "DiffuseFlows")], center = TRUE, scale = TRUE)
pca_res <- prcomp(X_pca)

summary(pca_res)

## Importance of components:
##                           PC1    PC2
## Standard deviation     1.2509 0.6598
## Proportion of Variance 0.7824 0.2176
## Cumulative Proportion  0.7824 1.0000

udaje$PC1_flows <- pca_res$x[, 1]

model_pca <- lm(
  PowerConsumption_Zone1 ~ Temperature + Humidity + WindSpeed + PC1_flows + PowerConsumption_Zone2,
  data = udaje
)

summary(model_pca)

## 
## Call:
## lm(formula = PowerConsumption_Zone1 ~ Temperature + Humidity + 
##     WindSpeed + PC1_flows + PowerConsumption_Zone2, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17417.7  -2482.1     65.5   2439.3  15346.9 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             5.853e+03  1.436e+02  40.765  < 2e-16 ***
## Temperature             1.869e+02  3.932e+00  47.527  < 2e-16 ***
## Humidity                7.646e+00  1.280e+00   5.973 2.34e-09 ***
## WindSpeed              -5.643e+01  8.160e+00  -6.916 4.70e-12 ***
## PC1_flows               7.979e+01  1.511e+01   5.282 1.28e-07 ***
## PowerConsumption_Zone2  1.072e+00  3.515e-03 305.060  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3813 on 52410 degrees of freedom
## Multiple R-squared:  0.7141, Adjusted R-squared:  0.7141 
## F-statistic: 2.619e+04 on 5 and 52410 DF,  p-value: < 2.2e-16

vif(model_pca)

##            Temperature               Humidity              WindSpeed 
##               1.885141               1.428806               1.324566 
##              PC1_flows PowerConsumption_Zone2 
##               1.287453               1.205603

7.5 Ridge Regression (voliteľné)

Ridge regresia stabilizuje odhady pomocou penalizácie. Koeficienty sú potom skreslené (biased), ale často s nižšou variabilitou.

library(MASS)

ridge_df <- udaje[, c("PowerConsumption_Zone1",
                      "Temperature", "Humidity", "WindSpeed",
                      "GeneralDiffuseFlows", "DiffuseFlows",
                      "PowerConsumption_Zone2")]

## Ridge sa správa lepšie na škálovaných premenných
ridge_df[, -1] <- scale(ridge_df[, -1], center = TRUE, scale = TRUE)

ridge_fit <- lm.ridge(
  PowerConsumption_Zone1 ~ .,
  data = ridge_df,
  lambda = seq(0, 50, 5)
)

ridge_fit

##             Temperature Humidity WindSpeed GeneralDiffuseFlows DiffuseFlows
##  0 32344.97    1120.128 96.69522 -131.5177           -90.11784     178.8008
##  5 32344.97    1120.127 96.52650 -131.4471           -90.06920     178.7398
## 10 32344.97    1120.127 96.35786 -131.3765           -90.02056     178.6787
## 15 32344.97    1120.126 96.18928 -131.3060           -89.97193     178.6177
## 20 32344.97    1120.126 96.02077 -131.2354           -89.92331     178.5568
## 25 32344.97    1120.126 95.85233 -131.1649           -89.87469     178.4958
## 30 32344.97    1120.125 95.68396 -131.0944           -89.82608     178.4349
## 35 32344.97    1120.125 95.51565 -131.0239           -89.77747     178.3740
## 40 32344.97    1120.124 95.34742 -130.9534           -89.72887     178.3132
## 45 32344.97    1120.124 95.17925 -130.8829           -89.68028     178.2523
## 50 32344.97    1120.123 95.01115 -130.8124           -89.63169     178.1915
##    PowerConsumption_Zone2
##  0               5576.181
##  5               5575.584
## 10               5574.987
## 15               5574.391
## 20               5573.795
## 25               5573.198
## 30               5572.602
## 35               5572.007
## 40               5571.411
## 45               5570.815
## 50               5570.220

plot(ridge_fit)

8. Zhrnutie

Multikolinearita nespôsobuje skreslenie koeficientov, ale zvyšuje ich štandardné chyby a robí koeficienty nestabilnými.
VIF zachytáva problém po premenných, condition number vie zachytiť aj širšie prepojenia medzi viacerými premennými.
Praktické riešenia: vynechanie premennej, škálovanie, zmena jednotiek, PCA alebo ridge.

Ekonometria master dokument

Adam Michalec

17. 12. 2025

Cvičenie 4 – Práca s databázou - import údajov, grafy, štatistiky

Import údajov z .csv

Grafy

Scatter plot

Boxplot

Základné štatistiky

Testovanie hypotéz

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

ANOVA – rozdiely medzi mesiacmi

Lineárna regresia

Cvičenie 6 – Econometrics in R - cvičenie 5 (Power consumption)

Úvod do problému a hypotézy

Príprava databázy, čistenie a úprava údajov

Výber premenných a tvorba časových premenných

Imputácia chýbajúcich hodnôt (ak by sa vyskytli)

Kontrola tvaru údajov (boxploty)

Lineárna regresia

Diagnostické grafy

Testy: normalita rezíduí a odľahlé pozorovania

Úprava modelu: log-transformácia žiarenia

Heteroskedasticita

Vizuálna kontrola (štvorce rezíduí)

Breusch–Pagan test

Robustné (heteroskedasticity-consistent) štandardné chyby

Záver

Cvičenie 7 – Špecifikácia modelu – spotreba elektriny

1. Príprava údajov

2. Základná lineárna regresia

3. Test správnej funkčnej formy – Ramsey RESET

4. Grafická diagnostika

4.1 Residuals vs Fitted

4.2 Component + Residual (C+R) grafy

5. Nelineárna špecifikácia (kvadratické členy)

6. Zlom cez dummy premennú (deň vs noc) + zmena sklonu

7. Box–Cox transformácia (voliteľné)

8. Robustné štandardné chyby (voliteľné)

Cvičenie 8 – Zhluková analýza (hierarchické klastrovanie) – spotreba elektriny

Úvod

Načítanie a príprava dát

Škálovanie premenných

Boxploty škálovaných premenných

Korelácie premenných

Matica vzdialeností (ukážka)

Hierarchické zhlukovanie (Wardova metóda)

Príslušnosť dní do klastrov

Centroidy (priemery premenných v klastroch)

Záver

Cvičenie 9 – Ekonometria v R – cvičenie 9 (autokorelácia rezíduí)

1. Úvod a údaje

1.1 Pracovná hypotéza

1.2 Načítanie a príprava dát

2. Lineárna regresia v základnom tvare

Graf: empíria vs. fitted

3. Autokorelácia rezíduí

3.1 Rezíduá a ACF

3.2 Durbin–Watsonov test

3.3 Breusch–Godfrey test

4. Ako riešiť autokoreláciu (prakticky)

4.1 Dynamizácia: lag závislej premennej (Koyckov typ uvažovania)

4.2 Newey–West robustné štandardné chyby

5. Záver

Cvičenie 10 – Cvičenie 10 – Multikolinearita v regresných modeloch (Power Consumption)

1. Čo je multikolinearita

2. Dáta a východiskový model

3. Odhad základného regresného modelu

4. Korelačná matica a scatterploty

5. VIF (Variance Inflation Factor)

6. Condition Number (podmienenosť)

7. Riešenia multikolinearity

7.1 Vynechanie problematickej premennej

7.2 Škálovanie (štandardizácia) vysvetľujúcich premenných

7.3 Jednoduchá zmena jednotiek (zachová interpretáciu)

7.4 PCA (voliteľné)

7.5 Ridge Regression (voliteľné)

8. Zhrnutie