11 Závislost dvou veličin

11.1 Kovariance

Míru těsnosti vztahu dvou proměnných lze vajádřit kovariancí

\[ \mathrm{cov}_{xy} = \dfrac{1}{n-1}\sum_{i=1}^{n}x_iy_i \]

Výsledek v podobě součinu není v praxi dobře interpretovatelný, proto se k porovnání těstnosti vztahu používá výpočtu korelace jako standardizované míry v rozmezí \(\langle-1;1\rangle\).

11.2 Korelace

Standardizovaný sdílený rozptyl dvou veličiny je korelace. Vyjadřuje, zda se jedna ze zkoumaných veličin v průměru mění, pokud se zároveň mění i druhá. Korelace neadresuje příčinnost, nelze tedy jednoznačně určit závislou a nezávislou proměnnou.

Míru korelace dvou veličin posuzujeme korelačním koeficientem. V případe normality u obou veličin lze použít Paersonův korelační koeficient,

\[ r_{xy} = \dfrac{\sum(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum(x_i - \bar{x})^2}\sqrt{\sum(y_i - \bar{y})^2}} \tag{11.1}\]

přičemž platí \(r_{xy}\in\langle-1;1\rangle\). Alternativně můžeme použít neparametrický výpočet Spearmanův, založený na diferencích pořadí pozorovaných hodnot, které definujeme jako \(d_i = x_{ri} - y_{ri}\)

\[ \rho_{xy} = r_s = 1 - \dfrac{6\sum_{i=1}^{n}d_i^2}{n(n^2-1)} \tag{11.2}\]

Cvičení

Napište vlastní funkci, které pro dané vektory x a y spočítá korelační Spearmanův koeficient.

Nejprve je nutné identifikovat pořadí,
dále spočítat druhou mocninu vzdáleností,
a po dosazení počtu měření \(n\) dosadit do rovnice.

Vyzkoušejte obě funkce na datech

Řešení:

Kód

sp_k <- function(x, y) {
  rank_x <- rank(x) 
  rank_y <- rank(y) 
  d <- rank_x - rank_y 
  d_sq <- sum(d^2) 
  n <- length(x) 
  1 - ((6 * d_sq) / (n*(n^2 - 1))) 
}

Kód

sp_k(_, _)
cor(_, _, method = "Spearman")

Pokuste se odhadnout korelační koeficient následujících veličin

s tolerancí 0.05. A) , B) , C) , D) . Pro desetinný rozvoj použijte tečku.

11.3 Regresní model

V regresní analýze posuzujeme příčinnou závislost (narozdíl od korelace). Identifikujeme nezávislou proměnnou a závislou proměnnou. Míru závislosti posuzujeme pomocí lineárního modelu.

\[ y = \beta_0 + \beta_1x \] kde \(y\) je závislá proměnná, \(x\) je nezávislá proměnná, \(\beta_0\) je počátek (intercept) a \(\beta_1\) je sklon přímky.

Lineární model regresního typu definujeme v R pomocí funkce lm(), které na první pozici dosazujeme výraz ve formátu rovnice formula, např. lm(y ~ x). Model pak vyhodnocujeme funkcí summary() a mezi sebou porovnáváme pomocí kritérií zohledňující vysvětlenou variabilitu/věrohodnost modelu a jeho komplexnost, např. Akaikeho - AIC() nebo Bayesovské informační kritérium BIC().

\[ \text{AIC} = 2k - 2\ln(\hat{L}) \] kde \(k\) je počet volných parametrů modelu a \(\hat{L}\) je věrohodnostní funkce modelu.

11.3.1 První regresní model

Uvažujme následující data, kde předpokládáme příčinnou závislost \(y\) na \(x\).

Kód

dfr <- data.frame(
  x <- c( 3.93,  3.83,  7.11,  6.98,  7.87,  9.11, 10.04,  9.99, 10.11, 10.93, 
         11.01, 11.93, 11.94, 11.88, 11.82, 12.98, 13.04, 13.04, 13.11, 13.99, 
         13.97, 13.99, 13.97, 15.08, 14.93, 14.86, 15.89, 16.15, 17.08, 17.18, 
         16.96, 17.95, 18.14, 18.09, 18.05, 19.04, 19.08, 19.07, 20.11, 20, 
         20,    20.14, 19.92, 22.13, 23.09, 23.9 , 23.95, 24.06, 23.95, 24.94
  ),
  y <- c( 1.93, 10.01,  3.73, 22.09, 15.86, 10.09, 18.09, 25.91, 34.09, 16.89, 
         28.04, 13.95, 19.98, 23.98, 27.78, 26.19, 33.85, 33.96, 45.88, 25.96, 
         35.89, 60.  , 80.03, 20.06, 25.92, 54.03, 31.96, 40.12, 32.03, 39.95, 
         50.2 , 41.91, 55.92, 75.92, 83.89, 35.94, 45.94, 67.85, 32.05, 48.11, 
         52.11, 55.95, 64.02, 65.9 , 53.94, 69.92, 92.1 , 92.78, 99.03, 85.09)
)

Sestrojíme lineární model funkcí lm() a uložíme do proměnné md1.

Kód

md1 <- lm(formula = y ~ x, data = dfr)

Struktura objektu md1 je komplexní, můžeme vybírat jednotlivé proměnné k dalším účelům, například pro tvorbu grafů.

Kód

str(md1, 1)

1: Vypíšeme strukturu objektu do první hierarchické úrovně

List of 12
 $ coefficients : Named num [1:2] -15.72 3.78
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
 $ residuals    : Named num [1:50] 2.79 11.25 -7.43 11.42 1.82 ...
  ..- attr(*, "names")= chr [1:50] "1" "2" "3" "4" ...
 $ effects      : Named num [1:50] -300.7778 140.5186 -9.3636 9.458 0.0675 ...
  ..- attr(*, "names")= chr [1:50] "(Intercept)" "x" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:50] -0.862 -1.24 11.163 10.672 14.037 ...
  ..- attr(*, "names")= chr [1:50] "1" "2" "3" "4" ...
 $ assign       : int [1:2] 0 1
 $ qr           :List of 5
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 48
 $ xlevels      : Named list()
 $ call         : language lm(formula = y ~ x, data = dfr)
 $ terms        :Classes 'terms', 'formula'  language y ~ x
  .. ..- attr(*, "variables")= language list(y, x)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..- attr(*, "term.labels")= chr "x"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(y, x)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
 $ model        :'data.frame':  50 obs. of  2 variables:
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ x
  .. .. ..- attr(*, "variables")= language list(y, x)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. ..- attr(*, "term.labels")= chr "x"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(y, x)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
 - attr(*, "class")= chr "lm"

Kód

with(data = dfr,
     expr = plot(x,
                 y,
                 xlim = c(0, 30),
                 ylim = c(md1$coefficients[1],
                          max(md1$fitted.values) + 2)))
abline(md1$coefficients, col = "orangered")
abline(v = 0, lty = "dashed")
abline(h = md1$coefficients[1], lty = "dashed")
text(x = 5,
     y = md1$coefficients[1] + 4, 
     labels = expression(paste(beta == round(md1$coefficients[1], 2))))

1: Určíme základní proměnné \(x\) a \(y\) v grafu a nastavíme limity okna.
2: Z výstupu objektu lineárního modelu vytáhneme hodnoty koeficientů \(\hat{\beta_0}\), \(\hat{\beta_1}\) a použijeme jako koeficienty přímky.
3: Dokreslíme průnik v počátku.

Kód

summary(md1)


Call:
lm(formula = y ~ x, data = dfr)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.274  -9.410  -1.639  10.063  42.925 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -15.7241     6.3224  -2.487   0.0164 *  
x             3.7816     0.3884   9.736 6.03e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.43 on 48 degrees of freedom
Multiple R-squared:  0.6639,    Adjusted R-squared:  0.6568 
F-statistic: 94.79 on 1 and 48 DF,  p-value: 6.032e-13

Souhrn obsahuje několik boduů k interpretaci kvality modelu.

Call: původní zadání závislosti v podobě rovnice.
Residuals: popisné statistiky residuí modelu.
3.Coefficients: identifikované koeficienty modelu
- Intercept odhad \(y\), když \(x=0\).
- x S každým zvýšením \(x\) se \(y\) zvýší o \(3.7816\); statisticky významný odhad ***
Residual standard error: Průměrná velikost residuí.
Multiple/Adujsted R-squared: Podíl vysvětlené variability v \(y\). Adjusted zohledňuje počet prediktorů a penalizuje složitější modely.
F-statistic: Testuje statistickou významnost modelu.

významnost závislosti indikuje přítomnost jedné nebo více * u vysvětlující proměnné; dále lze vyčíst podíl vysvětlené variabiltiy Adjusted R-squared: 0.6417.

\(R^2 = 1-\dfrac{\mathrm{MSE}}{\mathrm{MST}}\)

kde \(\mathrm{MSE}\) je střední kvadratická chyba (vzdálenost mezi skutečnými hodnotami a hodnotami predikovanými modelem) a \(\mathrm{MST}\) střední kvadratická odchylka (jak moc se jednotlivé hodnoty liší od průměru).

\[ \mathrm{MSE} = \dfrac{1}{n}\sum\limits_{i=1}^{n}(y_i-\hat{y_i})^2 \]

\[ \mathrm{MST} = \dfrac{1}{n}\sum\limits_{i=1}^{n}(y_i-\bar{y})^2 \]

Kód

par(mfrow = c(2, 2))
plot(md1, 1:4)

R dále poskytuje grafické nástroje k posouzení vhodnosti zvoleného modelu. Zde nás budou zajímat zejména residua modelu. Dle prvního grafu by mohl model se závislostí na polynomu být lepší volbou.

Residuals vs. fitted posuzuje homoskedasticitu (náhodný rozptyl residuí).
Normal Q-Q vynesení kvantilů residuí oproti kvantilům normálního rozdělení.
Zobrazení druhé odmocniny absolutních residuí vůči odhadnutým hodnotám k posouzení homoskedasticity.
Residuals vs. Levarage identifikuje hodnoty s výrazným pákovým efektem (schopností ovlivnit výslednou shodu modelu).

Pokusme se sestavit alternativní model s kvadratickou vysvětlující proměnnou.

Kód

md2 <- lm(y ~ poly(x, 2))
summary(md2)


Call:
lm(formula = y ~ poly(x, 2))

Residuals:
    Min      1Q  Median      3Q     Max 
-28.086  -9.010  -3.459   4.693  45.048 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   42.536      2.027  20.985  < 2e-16 ***
poly(x, 2)1  140.519     14.333   9.804 6.05e-13 ***
poly(x, 2)2   18.519     14.333   1.292    0.203    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.33 on 47 degrees of freedom
Multiple R-squared:  0.6754,    Adjusted R-squared:  0.6616 
F-statistic: 48.89 on 2 and 47 DF,  p-value: 3.29e-12

Call: ukazuje zadání modelu s přidaným polynomem až do 2. stupně.
Residuals:
Mají normální rozdělení Mají vyšší střední hodnotu než md1, což indikuje lepší shodu modelu Mají nižší rozptyl než md1, což indikuje lepší shodu modelu
.
3.Coefficients: identifikované koeficienty modelu
- V \(x=0\) má \(y\) hodnotu .
- poly(x, 1) S každým zvýšením \(x\) se \(y\) zvýší o \(3.7816\); statisticky významný odhad ***
Průměrná velikost residuí je: .
Model vysvětluje variability.

Kód

par(mfrow = c(2, 2))
plot(md2, 1:4)

Kód

with(data = dfr, 
     expr = plot(x, y, xlim = c(0, 30)))
pred <- predict(md2)
ix <- sort(dfr$x, index.return=T)$ix
lines(x[ix], pred[ix], col = 'orangered')
abline(v = 0, lty = "dashed")
abline(h = md2$coefficients[1], lty = "dashed")
text(x = 5,
     y = md2$coefficients[2] + 4, 
     labels = expression(paste(beta == round(md2$coefficients[2], 2))))

Srovnejme modely vzájemně

Kód

AIC(md1, md2)

    df      AIC
md1  3 412.8013
md2  4 413.0561