12 Závislost dvou veličin

12.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 \]

Takový výsledek, produkt, není v praxi dobře interpretovatelný.

12.2 Korelace

Standardizovaný sdílený rozptyl dvou veličiny je korelace. Vyjadřuje, zda se jedna ze zkoumaných veličin 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{12.1}\]

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

\[ \rho = r_s = 1 - \dfrac{6\sum_{i=1}^{n}d_i^2}{n(n^2-1)} \tag{12.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 kvadrát 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.

12.3 Regresní model

Narozdíl od korelace, v regresní analýze sledujeme příčinnou závislost. 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.

12.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 původní zadání v podobě rovnice, dále identifikované koeficienty 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.

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.

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

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')

Srovnejme modely vzájemně

Kód

AIC(md1, md2)

    df      AIC
md1  3 412.8013
md2  4 413.0561

Cvičení

Prozkoumejte závislost odtoku na velikosti povodí

Povodí	velikost \(x_i\)	pořadí	odtok \(y_i\)	pořadí	rozdíl \(d_i\)	\(d_i^2\)
1
2
3
4
5
6
7
8
9
10