Lineární regrese a korelace
Opakování: Načtení a výběr dat
Data ke stažení: pr1.zip
Instalace balíčků v R
Balíčky (packages) zavádíme příkazem
library(knihovna)Pokud obdržíme chybu
Error in library(knihovna) : there is no package called ‘knihovna’je nutné balíček jednorázově nainstalovat příkazem (vyžaduje internet)
install.packages("knihovna")Načtení tidy data
Nutno nastavit pracovní složku, aby fungovaly relativní cesty k souborům: Session - Set Working Directory - Choose Directory… (Ctrl+Shift+H)
library(readxl)
data <- read_excel("pr1.xlsx") # nepovinné parametry defaultně: sheet = 1, col_names = TRUEKontrola tabulky
Vykreslení
plot(data$x, data$y)
barplot(table(data$class))
Faktorové proměnné
Je velice vhodné specifikovat, které proměnné jsou faktorové. Obzvlášť, pokud se jedná o číselné a nikoliv textové hodnoty. Je rozdíl, zda bude model chápat hodnoty 1, 2 a 3 jako stejně vzdálená čísla na numerické ose, a nebo jen jako označení tří různých skupin.
str(data) # ověření typů proměnných## tibble [100 × 3] (S3: tbl_df/tbl/data.frame)
## $ x : num [1:100] 0.827 0.716 0.0962 0.8931 0.8173 ...
## $ y : num [1:100] 1.583 1.49 0.443 2.052 1.853 ...
## $ class: chr [1:100] "blue" "blue" "blue" "blue" ...data$class <- factor(data$class)
str(data) # ověření typů proměnných## tibble [100 × 3] (S3: tbl_df/tbl/data.frame)
## $ x : num [1:100] 0.827 0.716 0.0962 0.8931 0.8173 ...
## $ y : num [1:100] 1.583 1.49 0.443 2.052 1.853 ...
## $ class: Factor w/ 2 levels "blue","red": 1 1 1 1 1 1 1 1 1 1 ...Řadovou proměnnou je možné vytvořit dalšími parametry.
data$class <- factor(data$class, ordered = TRUE, levels = c("blue", "red"))
str(data) # ověření typů proměnných## tibble [100 × 3] (S3: tbl_df/tbl/data.frame)
## $ x : num [1:100] 0.827 0.716 0.0962 0.8931 0.8173 ...
## $ y : num [1:100] 1.583 1.49 0.443 2.052 1.853 ...
## $ class: Ord.factor w/ 2 levels "blue"<"red": 1 1 1 1 1 1 1 1 1 1 ...Výběr řádků na základě podmínek
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionvyber <- data %>% filter(class == "red")
vyber## # A tibble: 50 × 3
## x y class
## <dbl> <dbl> <ord>
## 1 0.429 2.03 red
## 2 0.855 3.02 red
## 3 0.589 2.14 red
## 4 0.760 1.89 red
## 5 0.442 1.76 red
## 6 0.0160 0.755 red
## 7 0.531 1.63 red
## 8 0.962 2.61 red
## 9 0.258 1.37 red
## 10 0.817 2.21 red
## # ℹ 40 more rowsLineární regrese
Zpracováno na motivy přednášky Regression Models v rámci specializace Data Science na Coursera (prof. Brian Caffo z Johns Hopkins Bloomberg School of Public Health).
Pracujme s daty v interní proměnné mtcars (více info získáte příkazem ?mtcars).
plot(mpg ~ wt, mtcars) # závislost mpg (míle na galon) na hmotnosti auta (wt), patrný je trend
fit <- lm(mpg ~ wt, data = mtcars) # lineární model mpg v závislosti na wt
abline(fit) # přidání přímky podle lineárního modelu fit
fit # výpis základních parametrů modelu##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Coefficients:
## (Intercept) wt
## 37.285 -5.344Význam parametrů
- Intercept (počátek, výchozí hodnota, průsečík s vertikální osou): hodnota mpg při wt = 0, neboli vertikální posun přímky při nulové hodnotě vstupního parametru (tzv. prediktoru či nezávislé proměnné). Tato hodnota je spíše hypotetická, protože automobil o nulové hmotnosti neexistuje, ale je důležitá pro samotný model a jednoduše můžeme dopočítat konkrétní hodnoty pro reálné rozsahy hmotností.
- wt: sklon regresní přímky (tzv. slope). S každým zvýšením wt o 1 klesne mpg o -5.344.
Skvělou vlastností lineárních modelů je jejich jednoduchá interpretovatelnost. Vidíme, jak určitá změna jednoho parametru změní hodnotu jiného parametru.
Důležité je však uvědomit si, že vazba mezi proměnnými neznamená kauzalitu. Nelze tak na základě této analýzy jednoduše říct, že jedna věc způsobuje druhou. Vztah může být klidně zcela opačný, obě proměnné se mohou ovlivňovat navzájem, a nebo být ovlivněny zcela jinou proměnnou v pozadí, kterou jsme do modelu nezahrnuli. O těchto problémech se zmíníme dále a určitým řešením může být použití vícenásobné regrese (zahrnutí více proměnných do modelu), ale ani v takovém případě nemůžeme mluvit o analýze kauzality (odhalování příčin a následků).
Koeficienty modelu jsou hledány tak, aby byly minimalizovány odchylky skutečných hodnot od proloženého trendu (model). Nejčastěji se tak děje tzv. metodou nejmenších čtverců, která se snaží minimalizovat součet jednotlivých odchylek umocněných na druhou.
Proložení dat křivkou je vždy jen odhadem skutečného trendu základního souboru, můžeme proto zobrazit podrobnější statistiku.
coef(fit) # jen samotné koeficienty (bodový odhad)## (Intercept) wt
## 37.285126 -5.344472confint(fit) # konfidenční intervaly odhadu koeficientů## 2.5 % 97.5 %
## (Intercept) 33.450500 41.119753
## wt -6.486308 -4.202635Konfidenční pásmo pro regresní přímku též automaticky vypočítá a vykreslí nástroj v knihovně ggplot2.
library(ggplot2)
ggplot(mtcars, aes(x = wt, y = mpg)) + geom_point() + stat_smooth(method = "lm")## `geom_smooth()` using formula = 'y ~ x'
summary(fit) # na konci řádku vidíme pro každý koeficient p-hodnotu Pr(>|t|) pro H0: koeficient = 0##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10Vidíme tedy, že oba koeficienty mají velice nízkou p-hodnotu, takže např. na hladině významnosti \(\alpha = 0.05\) můžeme zamítnout H0, že Intercept (výchozí hodnota mpg při nulové wt) i sklon přímky (wt) jsou nulové.
Důležitý je parametr R-squared (\(r^2\) nebo \(R^2\)). Nabývá hodnot mezi 0 a 1 a odpovídá poměru rozptylu popsaného modelem ku celkovému rozptylu modelované veličiny (regresní rozptyl / celkový rozptyl dat). Čím vyšší číslo, tím větší část variability dat model popisuje. Platí totiž celkový rozptyl dat = rozptyl odchylek (viz další kapitola) + regresní rozptyl, tedy regresní rozptyl = celkový rozptyl - rozptyl odchylek.
Čím menší rozptyl mají zbytkové odchylky od regresního trendu, tím větší část rozptylu dat model popisuje.
Poznámka: pokud model zahrnuje člen Intercept, R-squared je rovna druhé mocnině korelace mezi vstupními proměnnými (která se často označuje r), viz kapitola o korelaci dole v tomto dukumentu.
Adjusted R-squared je pokusem provést korekci jevu, kdy ve vícenásobné regresi (více vstupních proměnných - regresorů) přidávání dalších (i nerelevantních) proměnných způsobuje nárůst základního koeficientu R-squared. Na rozdíl od základního R-squared narůstá tato upravená verze pouze v případě, když je přínos nové proměnné nad rámec náhody. Tato hodnota nevychází vyšší než základní R-squared a může nabývat i záporných hodnot. Nemá proto jednoduchou interpretaci jako základní R-squared a používá se spíše jako měřítko pro vyhodnocení, zda přidávání dalších proměnných do modelu má skutečně smysl. Pro takové aplikace však existují i jiná kritéria, jako např. BIC či AIC.
Predikce
Chceme pomocí modelu predikovat očekávané hodnoty např. pro wt = 3, 4 a 5. Tím dostáváme hodnoty ležící přímo na regresní přímce.
noveWT <- c(3, 4, 5)
predict(fit, newdata = data.frame(wt = noveWT))## 1 2 3
## 21.25171 15.90724 10.56277Funkce predict() nabízí výpočet intervalů spolehlivosti jak pro určení možného výskytu regresní přímky (střed), tak skutečných běžných dat (tzv. predikční interval - zohledňuje variabilitu dat).
Dosadíme celý rozsah hodnot wt a z intervalů spolehlivosti se nám stanou pásy spolehlivosti.
rozsahWT <- seq(min(mtcars$wt), max(mtcars$wt), length.out = 100)
noveWT <- data.frame(wt = rozsahWT)
# a) intervals for the regression line at point wt
p1 <- predict(fit, noveWT, interval = ("confidence"))
# b) the prediction of what a mpg would be at point wt (prediction interval)
p2 <- predict(fit, noveWT, interval = ("prediction"))
plot(mtcars$wt, mtcars$mpg, pch=21, col="black", bg="lightblue", cex=1)
abline(fit, lwd = 2)
lines(rozsahWT, p1[,2]); lines(rozsahWT, p1[,3])
lines(rozsahWT, p2[,2], col = "lightblue"); lines(rozsahWT, p2[,3], col = "lightblue")
Knihovna ggplot2 přímo neumožňuje zobrazit u statistiky predikční pásmo, proto ho musíme nejdříve samostatně vyčíslit.
mtcars_pred <- cbind(mtcars, predict(fit, data.frame(wt = mtcars$wt), interval = "prediction"))
library(ggplot2)
ggplot(mtcars_pred, aes(x = wt, y = mpg)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "blue", alpha = 0.2) +
geom_point() +
stat_smooth(method = "lm")## `geom_smooth()` using formula = 'y ~ x'
Residuals (odchylky)
Základním předpokladem modelu je, že odchylky od trendu se řídí normálním rozdělením, jedná se tedy o představu „trendu“ a „šumu s normálním rozdělením.“
Než začneme model skutečně používat či interpretovat, měli bychom si tento předpoklad ověřit.
par(mfrow = c(1, 2))
plot(mtcars$wt, resid(fit)) # odchylky v závislosti na vstupní proměnné náhodně rozptýleny okolo středu
plot(predict(fit), resid(fit)) # závislost predikovaných hodnot (trend) na zbytkových hodnotách, viz korelační graf dále - tvarem by se mělo jednat o "mlhovinu"
shapiro.test(resid(fit)) # ani test na normalitu nemůže vyloučit normální rozdělení##
## Shapiro-Wilk normality test
##
## data: resid(fit)
## W = 0.94508, p-value = 0.1044Jakmile by se v grafu zbytkových odchylek ukazovaly nějaké vzory, bude v datech pravděpodobně nelineární souvislost, kterou náš model nedocháže podchytit.
Dalším problémem může být tzv. heteroskedasticita, kdy je velikost rozptylu závislá na velikosti vstupní proměnné. U takových dat nemají smysl p-hodnoty a další věci, které předpokládají rovnoměrný rozptyl (to se týká ale například i t-testů). V takovém případě jsou také pravidla modelu porušena a je dobré přemýšlet, zda netransformovat proměnné do jiných veličin, které by byly pro model vhodnější (např. zkusit data logaritmovat).
par(mfrow = c(1, 2))
x <- runif(100, 0, 6); y <- x + rnorm(100, mean = 0, sd = .001*x);
fit <- lm(y ~ x)
plot(x, y); abline(fit) # vypadá jako perfektní model fit
plot(x, resid(fit)) # vidíme velice jasně heteroskedasticity
abline(h = 0)
Outliers (body ležící výrazně mimo)
Některé body vadí více, jiné méně.
- Levý dolní roh: naprosto „hodný“ bod, leží „nevybočuje z davu“ a neovlivňuje parametry modelu.
- Levý horní roh: typický outlier. Vzhledem k tomu, že jeho souřadnice X se nachází v oblasti, kde je i velké množství „běžných“ bodů, tato většina ho „převálcuje“ a jeden takový outlier příliš parametry modelu nezmění.
- Pravý horní roh: sice je to také na první pohled outlier pro svou netypickou hodnotu X (a rozhodně by stálo za to, zamyslet se, proč tato hodnota leží zcela jinde než zbytek ostatních), ale naštěstí leží na regresní přímce, takže koeficienty modelu neovlivní.
- Pravý dolní roh: tento bod silně vybočuje, a proto má velký vliv, jeho zahrnutí do modelu by koeficienty silně ovlivnilo.
Pro analýzu vlastností outliers můžeme v R používat celou řadu diagnostických nástrojů. Například funkce hatvalues ukazuje pro každou položku dat její odchylku od modelu, konkrétně 1 – poměr residual modelu s a bez outlieru (čím blíže 1, tím větší vliv). Můžeme tedy hledat body, které leží daleko od přímky. Funkce dfbetas ukazuje, jak hodně by se změnily koeficienty modelu, když vyloučíme ten který bod (tímto nalezneme ojedinělé body se silným efektem). Můžeme také použít přímo funkci plot, které dáme jako parametr náš model, a ona zobrazí hned několik diagnostických grafů.
Příklad: normálně rozdělená data a jeden outlier se souřadnicemi [10,10] (schválně ho umístíme hned na první pozici v datech)
n <- 100; x <- c(10, rnorm(n)); y <- c(c(10, rnorm(n)))
plot(x, y, frame = FALSE, cex = 2, pch = 21, bg = "lightblue",
col = "black")
fit <- lm(y ~ x); abline(fit, lwd=2)
fit2 <- lm(y[-1] ~ x[-1]); abline(fit2) # fit všech dat kromě 1. bodu
round(hatvalues(fit)[1 : 10], 3) # prvních 10 hodnot - 1. bod (náš outlier [10,10] trochu vybočuje)## 1 2 3 4 5 6 7 8 9 10
## 0.484 0.020 0.012 0.010 0.010 0.021 0.011 0.017 0.010 0.028round(dfbetas(fit)[1 : 10, 2], 3) # prvních 10 hodnot, 2. koeficient (sklon) je nejvíce ovlivněn 1. bodem## 1 2 3 4 5 6 7 8 9 10
## 5.958 0.012 -0.035 -0.019 0.019 -0.198 -0.010 -0.068 0.022 0.012par(mfrow = c(2, 2))
plot(fit) # na všech grafech vidíme, že bod s indexem 1 je zcela mimo
Problémy
Problémy nastávají, pokud do lineárního modelu nezahrneme všechny proměnné, které mají na výsledek vliv (což u jednoduché lineární regrese ani nemůžeme, potřebovali bychom vícenásobný regresní model).
Simpson’s paradox
Jasný trend v oddělených skupinách zmizí (nebo se zcela převrátí), když se skupiny sloučí. Pokud do modelu vložíme skutečné kauzální vztahy, Simpsonův paradox zmizí.
Příklady simulací – dva způsoby léčby
Ukázky, kdy může vícenásobná regrese poskytnout výrazně lepší náhled na souvislosti v datech než např. t-testy nebo ANOVA.
Př. 1 X nabývá u obou skupin podobného rozsahu hodnot. X souvisí s Y, ale je rozdílný intercept (výchozí hodnota neboli posun) podle skupiny, mezi skupinami je tedy konstantní posun v Y. Sloučit obě skupiny dohromady by byla chyba. Dostali bychom sice krásnou regresi, ale s velkým nevysvětlitelným rozptylem. Přitom skupina krásně vysvětlí, jestli jsme nahoře nebo dole.
Př. 2 Rozsah X mají obě skupiny zcela jiný. X souvisí s Y, intercept a slope nezávisí na skupině. Kdybychom vzali jen průměry, jedna skupina bude působit výrazně lépe než druhá. Není to ale celé jen tím, že jsme jednu skupinu měřili pro jiný rozsah X než tu druhou? Není to celé vlastně jen jedna skupina? Není efekt skupiny nulový? A nebo možná naopak – každá skupina leží celkově někde jinde, ale z koeficientů lineární regrese to nepoznáme.
Př. 3 Z hlediska průměrů vychází červená výše než modrá. Když však extrapolujeme modrá data pro větší X, je jasné, že modrá leží na vyšší přímce! Rozsah X má každá skupina jiný (modrá má nižší X, červená vyšší). Bylo by chybou skupiny spojit do jedné, pak by vyšel nesmyslný sklon.
Př. 4 X nesouvisí se skupinou. Průměr červené a modré je prakticky stejný (na rozdíl od příkladu 1). Silný posun v Y ve vztahu ke skupině. Při zafixovaném X je dost důkazů pro porovnání obou skupin.
Př. 5 Průměry skupin jsou stejné. Není tu konstantní posun skupin (není tu tzv. group effect, který bychom nalezli t-testem). Skupina však zcela převrací směr v závislosti na X. Intercept i slope závisejí na skupině. Rozsah X mají obě skupiny stejný. Zafixujeme-li X, je zde dost informace o efektu skupiny (pokud se zrovna netrefíme do X = 0), akorát pro různá X pokaždé jiný…
Př. 6 Tady nejsou diskrétní skupiny, ale spojitá veličina (barva), která data někam posouvá. Ale vždy při konkrétní zafixované hodnotě této veličiny je jiný intercept a zdá se, že sklon je stejný. Kdybychom tuto veličinu do modelu nezahrnuli, sklon sice najdeme, ale vzhledem k rozptylu dat na výšku nám přijde příliš málo významný (bude velice malé \(R^2\)).
Korelace
Korelační koeficient je sklon regresní přímky, když jsou obě proměnné normalizovány ve stylu z-skóre. Korelační koeficient je tedy úzce svázán s lineární regresí.
Korelace však neznamená kauzalitu! To, že nalezneme souvislost mezi dvěma proměnnými, vůbec nemusí znamenat, že jedna proměnná ovlivňuje druhou. Souvislost může být zapříčiněna společným faktorem (nebo faktory) „v pozadí“.
Nulová korelace neznamená, že mezi proměnnými není vazba! Korelační koeficienty jsou definovány jen pro určité situace. Pokud splníme jejich podmínky, funguje dobře. Pokud mezi proměnnými není žádná vazba, korelace je nulová. Opačné tvrzení ale automaticky nefunguje. Vazba může být komplikovanější (např. nelineární vztah, jiné rozdělení náhodných veličin apod.) a korelační koeficient ji neodhalí. Je proto vždy velice důležité před samotným výpočtem korelačního koeficientu zobrazit korelační graf, viz dále.
Kladné znaménko korelace říká, že pokud se jedna hodnota zvýší, druhá hodnota se také trochu zvýší (přestože obě obsahují stále i silnou náhodnou složku, tak tento trend je dlouhodobě z mnoha párů patrný).
Korelace může mít záporné znaménko (opačný směr vztahu – když jde jedna proměnná nahoru, druhá opačně dolů).
Korelace, kterou jsme poznali v akustice, slouží k detekci podobností mezi signály, pokud je vzájemně různě posunujeme. Obdobně korelační koeficient ve statistice detekuje souvislosti mezi daty. Vzhledem k tomu, že mají náhodný (stochastický) charakter, nemá smysl se bavit o časové souslednosti, zkoumáme vždy jen příslušné páry hodnot (bez posunů).
Podle typu proměnných rozlišujeme různé vzorce na výpočet korelačního koeficientu.
- Pearsonův korelační koeficient: pro metrická data, normálně rozdělená.
- Spearmanův korelační koeficient: pro data bez normálního rozdělení, vhodný i pro ordinální proměnné (vyhodnocuje se pouze pořadí, nikoliv velikost).
- Kendallovo tau: opět vyhodnocuje shodu mezi pořadím, jiný přístup.
- Koeficient fí: pro binární data (veličiny nabývají jen dvou hodnot).
První tři typy lze v R vypočítat funkcí cor(x, y), s nepovinným parametrem method = … a možnými hodnotami “pearson” (defaultní, není nutno zadávat), “kendall” a “spearman”.
Poslední koeficient fí lze spočítat přímo metodou “pearson” ve funkci cor(x, y).
Korelační graf (scatter plot)
Jedná se o bodový X-Y graf, kde vynášíme společně vždy příslušné dvojice hodnot x a y.
Pakliže spolu obě veličiny vůbec nejsou vázány, oba páry budou kreslit „mlhovinu“ bodů dle jejich pravděpodobnostního rozdělení - v případě normálního rozdělení očekáváme elipsu, kde ve středu budou hodnoty nejčastěji, k okrajům méně často. Při rovnoměrném rozdělení obdržíme rovnoměrně rozdělený obdélník. V případě nějaké vazby mezi veličinami se budou tyto obrazce stáčet (stále poloha bodů vůči sobě obsahuje náhodnou složku, ale pokud jde jeden bod jedním směrem, druhý díky vazbě také půjde buď stejným či opačným směrem). Pokud by byla korelace 100 %, body budou ležet na šikmé přímce. Pomocí tohoto grafu je možné objevit i nelineární vztahy.
# Příklad nekorelovaných proměnných, vykreslení scatter plot a výpočet korelace
n <- 10000
x1 <- rnorm(n)
y1 <- rnorm(n)
cor(x1, y1)## [1] 0.001317395# Lineární vazba
x2 <- rnorm(n)
y2 <- 1.5*x2 + rnorm(n)
cor(x2, y2)## [1] 0.8338571# Nelineární vazba - korelační koeficient nic nepozná
x3 <- rnorm(n)
y3 <- 1.5*x3^2 + rnorm(n)
cor(x3, y3)## [1] -0.02323575par(mfrow = c(1, 3))
plot(x1, y1, pch = 1, cex=0.3, asp = 1)
plot(x2, y2, pch = 1, cex=0.3, asp = 1)
plot(x3, y3, pch = 1, cex=0.3, asp = 1)
Různé vazby a jejich Pearsonův korelační koeficient
Všechny tyto případy ukazují, že korelační koeficient nemůže nahradit sílu vizuálního vyhodnocení dat.
První řádek ukazuje použití korelačního koeficientu tak, jak „byl myšlen“, tedy bez porušení podmínek.
Druhý řádek ukazuje jeden speciální případ stoprocentní vazby proměnných (jedna proměnná je stoprocentně závislá na druhé proměnné bez další náhodné složky), jen s různým koeficientem (např \(y = 3x\)). Vzhledem k normování uvnitř korelace vychází koeficient korelace stále 1, jen v prostředním případu vychází NaN (not-a-number), prože se jedná o vztah ve stylu \(y = 0x\).
Třetí řádek znázorňuje několik nelineárních vazeb, které korelační koeficient neodhalí.
by
DenisBoigelot,
original uploader was
Imagecreator
- Own work, original
uploader was
Imagecreator.
Licensed under
CC0
via Wikimedia Commons.
Anscombeho kvartet (Anscombe’s quartet)
Čtveřice zcela rozdílných dat, která přesto mají společné souhrnné ukazatele: průměr, rozptyl, korelaci a regresní přímku. Dostupná např. v interní proměnné anscombe.
anscombe## x1 x2 x3 x4 y1 y2 y3 y4
## 1 10 10 10 8 8.04 9.14 7.46 6.58
## 2 8 8 8 8 6.95 8.14 6.77 5.76
## 3 13 13 13 8 7.58 8.74 12.74 7.71
## 4 9 9 9 8 8.81 8.77 7.11 8.84
## 5 11 11 11 8 8.33 9.26 7.81 8.47
## 6 14 14 14 8 9.96 8.10 8.84 7.04
## 7 6 6 6 8 7.24 6.13 6.08 5.25
## 8 4 4 4 19 4.26 3.10 5.39 12.50
## 9 12 12 12 8 10.84 9.13 8.15 5.56
## 10 7 7 7 8 4.82 7.26 6.42 7.91
## 11 5 5 5 8 5.68 4.74 5.73 6.89c(mean(anscombe$x1), mean(anscombe$x2), mean(anscombe$x3), mean(anscombe$x4))## [1] 9 9 9 9c(mean(anscombe$y1), mean(anscombe$y2), mean(anscombe$y3), mean(anscombe$y4))## [1] 7.500909 7.500909 7.500000 7.500909c(var(anscombe$x1), var(anscombe$x2), var(anscombe$x3), var(anscombe$x4))## [1] 11 11 11 11c(var(anscombe$y1), var(anscombe$y2), var(anscombe$y3), var(anscombe$y4))## [1] 4.127269 4.127629 4.122620 4.123249c(cor(anscombe$x1, anscombe$y1), cor(anscombe$x2, anscombe$y2), cor(anscombe$x3, anscombe$y3), cor(anscombe$x4, anscombe$y4))## [1] 0.8164205 0.8162365 0.8162867 0.8165214example(anscombe)
© 10. 5. 2015 Tomáš Bořil, borilt@gmail.com