######################################## Zpracování dat fyzikálních měření - 5. díl ################################### ########## Načítání vstupních dat ########## # Načtení vstupních dat ze souboru formátu csv a úprava, aby byly všechna data chápána jako čísla (typ numeric). # Místo ... doplňte cestu k souboru. Pro načítání vstupních dat vždy používejte tuto dvojici příkazů. mereni <- read.csv("vzor.csv", sep = ",") mereni <- data.frame(lapply(mereni, as.numeric)) mereni # Pro účely simulace použitých metod se v tomto skriptu budou používat vzorová data (soubor vzor.csv). # Celý vzorový skript předpokládá, že máte načtený tento datový soubor, pro vyřešení seriálových úloh # stačí načíst jiná data a mírně upravit příslušné příkazy. # V případ jakýchkoliv nejasností zkuste použít Help (v pravém spodním okně záložka Help a do příslušného textového # pole napište jméno funkce) nebo zkuste najít odpoveď na internetu nebo napište autorovi seriálu na mail nozicka@fykos.cz # Na většinu metod popsaných v seriálu exitují jednoduché funkce, takže # ani není třeba pamatovat si složité vzorce nebo do nich dosazovat. # Důležité je umět použít správný postup a umět správně interpretovat výsledky. # Na toto se zaměříme ve vzorovém skriptu. ########## Instalace potřebných knihoven ########## install.packages("investr") # Instalovat stačí spustit jednou. library(investr) # Importování knihovny musíte provést pokaždé, když znovu zapínáte RStudio. ########## Lineární regrese ########## #### Graf s naměřenými daty plot(mereni$x, mereni$y0) # první argument jsou hodnoty na ose x, druhý argument jsou hodnoty na ose y # další volitelné argumenty funkce plot() jsou: # - xlab = "...", ylab = "..." - popisky os # - main = "..." - nadpis grafu plot(mereni$x, mereni$y0, xlab = "hodnoty x", ylab = "hodnoty y", main = "První graf") # legenda se přidá pomocí příkazu legend(...) legend("topleft", legend = c("naměřené hodnoty"), pch = c(1)) # funkce legend má 2 povinné argumenty # - první argument udává polohu v grafu (možnosti "topleft", "bottomright" atd.). Alternativně je možné legendu # umístit pomocí argumentů x = ..., y = ... (kde x,y udávají polohu horního levého rohu legendy). # - argument legend = ... udává, jaké popisky se mají v legendě objevit. Pokud je popisků více, musíme zadat jako # pole, o což se stará příkaz c("popisek1", "popisek2", "popisek3"). Každý popisek bude na jednom řádku. # - argument pch = ... udává typ zobrazeného symbolu (pch = point character). Můžete vyzkoušet různé typy symbolů # pomocí zadávání různých čísel (stejný argument bere i příkaz plot(), pomocí něhož kreslí různé symboly). # Pokud má být symbolů více, musíme je zadat jako pole pomocí příkazu c(..., ...). #### Odhadování parametrů v lineárním regresním modelu # O odhadování parametrů se stará funkce lm() (lm = linear model). model1 <- lm(y0 ~ x, data = mereni) # První argument je tvar prokládané funkce. # - vždy se požaduje tvar: ~ f1() + f2() + ... + fk() # - vžde se implicitně uvažuje, že je zahrnutý absolutní člen beta_0 (značení jako v textu seriálu), pokud bychom # nechtěli mít absolutní člen zahrnutý musíme připsat ještě -1 # Druhý argument data = ... specifikuje, jaká data se mají pro odhad použít (příslušná proměnná musí obsahovat sloupce # pojmenované stejně jako použité závisle a nezávisle proměnné). # Grafické znázornění prokládané funkce (tento příkaz potřebuje knihovnu investr) plotFit(model1) # Tomuto příkazu se budeme ještě dále věnovat. #### Pro shrnutí odhadů v našem lineárním regresním modelu použijeme funkci summary() summary(model1) # Nyní něco k interpretaci toho, co dostaneme na výstupu: # - Call: popisuje, jakou funkci jsme zavolali. # - Residuals: popisuje vlastnosti residuí našeho modelu. 1Q udává 1. kvartil (pozor na rozdíl kvantil vs. kvartil), # což je taková hodnota, že 1/4 residuí je menší než toto číslo a 3/4 residuí jsou větší než toto číslo. # Podobně Median a 3Q (3. kvartil) udávají hodnoty v polovině resp. ve 3/4 residuí. # - Coefficients: udávají odhady regresních koeficientů. Každý řádek odpovídá jednomu regresnímu koeficientu, # přičemž (Intercept) značí absolutní člen (beta_0 ve značení jako v seriálu). # - Ve sloupci Estimate je hodnota odhadu příslušného koeficentu. # - Ve sloupci Std. Error je nejistota určení tohoto koeficientu (spočítaná pomocí vzorců uvedených v seriálu) # - Sloupce t value a Pr(>|t|) se vztahují ke statistickému testu o nulovosti příslušného regresního koefecientu, # o těchto testech si povíme v příštím dílu seriálu. Zatím můžeme tyto dva sloupce ignorovat. # - Řádek Signif. codes se stále vztahuje k testům o nenulovosti regresních koeficientů -> prozatím ignorujme # - Residual standard error udává odhad parametru sigma (podle značení v seriálu). degrees of freedom udává počet # měření minus počet regresních parametrů. # - Multiple R-squared udává koeficient determinace R^2 modelu, Adjusted R-squared udává upravedný koeficient # determinace (přesně jak bylo popsáno v textu seriálu). # - Poslední řádek F-statistics: budeme ignorovat. #### Konstrukce intervalových odhadů pro regresní koeficienty je jednoduchá, použijeme zabudovanou funkci. # Výstup je spočítán pomocí vzorců uvedených v seriálu. # Příklad pro obecnou hladinu (1 - alpha) a oba regresní koeficienty najednou alpha <- 0.05 confint(model1, level = 1 - alpha) #### Složitější regresní modely # Pokud chceme použít regresní model, ve kterém budeme prokládat složitější funkci, bude se syntaxe trochu lišit. # Představme si například, že chceme prokládat naměřenými daty hyperbolu, tedy funkci beta0 + beta1 * (1 / x). model1 <- lm(y1 ~ I(1/x), data = mereni) # Abychom od sebe oddělili jednotlivé funkce (vzpomeňte si, že v seriálu se psalo, že vždy prokládáme funkci typu # beta0 + beta1 * f1(x) + beta2 * f2(x) + ... betak * fk(x), potřebujeme tedy jasně říci, co je f1, f2, ..., fk) # používáme symbol I(...). Jedině pro funkci f(x) = x není nutné psát I(x), ale stačí prostě x. # Každá jednotlivá funkce f1(x), ... , fk(x), tedy musí být zadána uvnitř I(...) # Pokud bychom chtěli prokládat např. kubickou funkci, museli bychom psát y ~ I(x) + I(x^2) + I(x^3) nebo alternativně # by stačilo y ~ x + I(x^2) + I(x^3) # Shrnutí odhadů regresních koeficientů (interpretace stejná jako u předchozího modelu) summary(model1) #### Kreslení grafů # Nyní si podrobněji rozebereme, jak se kreslí pěkné grafy. # Základem je funkce plotFit plotFit(model1) # K tomuto základnímu příkazu se může přidat mnoho dalších argumentů a vznikne hezčí graf. # Mnoho těchto dodatečných argumentů je naprosto stejných jako u příkazu plot() plotFit(model1, xlab = "hodnoty x", ylab = "hodnoty y", main = "Graf proložené funkce") legend("topright", pch = c(1,NA), lty = c(NA,1), legend = c("naměřené hodnoty", "proložená funkce")) # Argument lty (lty = line type) udává jaký typ čáry se má zobrazit v legendě (1 = plná čára). # Když se používají argumenty lty a/nebo pch musí vždy mít stejný počet hodnot jako je řádků legendy (tj. kolik # jednotlivých stringů (string je to mezi uvozovkami) je v proměnné legend). Pokud nechceme aby se někam nakreslila # čára nebo symbol (jako třeba v tomto případě, kdy chceme na první řádce pouze symbol a na druhé pouze čáru), # zadáme místo čísla NA (NA = not availible = není dostupný) # Dále můžeme využívat různé barvy proložené funkce, o což se stará argument col.fit = "...". # Jednotlivé hodnoty tohoto argumentu jsou prostě jména barev, např. "red", "blue", ... plotFit(model1, xlab = "hodnoty x", ylab = "hodnoty y", main = "Graf proložené funkce", col.fit = "red", pch = 4) legend("topright", pch = c(NA,4,NA), lty = c(NA,NA,1), legend = c(expression(paste(y, " = ", a/x + b)),"naměřené hodnoty", "proložená funkce"), col = c(NA, "black", "red")) # Pokud chceme do legendy psát matematické výrazy, použijeme funkci expression(...) místo stringu "...". Pokud chceme spojit matematické # výrazy a stringy musíme použít funkci paste(..., ...) (paste = slepit). Použití vidíte na tomto příkladu, správná syntaxe je # expression(paste(..., ..., ...)) #### Intervalový odhad pro prokládanou funkci # Pro výpočet intervalového odhadu prokládané funkce složí funkce predict. Argument funkce predict musí být proměnná # typu data.frame s hodnotou (nebo hodnotami) nezávisle proměnné, ve které chceme spočítat intervalový odhad pro # prokládanou funkci. Druhý argument je interval, který má hodnotu "confidence" (confidence interval = interval spolehlivost) # Představme si, že chceme interval spolehlivosti pro prokládanou funkci v bodě x = 5 na hladině spolehlivosti # (1 - alpha) = 0.95 (argument level = ...) predFit(model1, newdata = data.frame(x = 5), interval = "confidence", level = 0.95) # fit udává hodnotu proložené funkce, lwr (lwr = lower = spodní) a upr (upr = upper = horní) udávají spodní a horní # hranici intervalu spolehlivosti. # to samé pro více hodnot najednou predFit(model1, newdata = data.frame(x = c(5, 6, 7, 8)), interval = "confidence", level = 0.95) # U této funkce je velice důležité, aby se shodovali jména proměnných v regresním modelu a jména nových dat. Pokud # se neshodují funkce zahlásí error. Velice často se stává, že na toto člověk zapomene, proto si na to dávejte pozor. # To samé znázorněno graficky pomocí funkce plotFit (nemělo by nás překvapit, že mnoho naměřených hodnot padne # mimo interval spolehlivosti, je to tím, že díky velkému počtu měření máme prokládanou funkci určenou poměrně přesně) plotFit(model1, interval = "confidence", xlab = "hodnoty x", ylab = "hodnoty y", main = "Graf proložené funkce", col.fit = "red", col.conf = "blue") # Abychom dostali takovýto graf museli jsme přidat argument interval = "confidence". Dále jsme upravili barvy # pomocí příkazů col.fit = ... (barva proložené funkce) a col.conf = ... (barva intervalu spolehlivosti) # Případně bychom mohli ještě upravovat pomocí příkazů pch = ... (typ symbolu zobrazujících jednotlivá měření) # lty.fit = ... (typ čáry proložené funkce, defaultní je hodnota 1), lty.conf = ... (typ čáry intervalu spolehlivosti, # defaultní je 2) # Přidáme legendu pomocí známého příkazu legend() legend("topright", legend = c(expression(paste(y, " = ", a/x + b)), "naměřené hodnoty", "proložená funkce", "interval spol. (95 %)"), pch = c(NA, 1, NA, NA), lty = c(NA, NA, 1, 2), col = c(NA, "black", "red", "blue")) ### TAKTO NĚJAK BY MĚL IDEÁLNĚ VYPADAT SPRÁVNÝ GRAF. Co hlavně nesmí chybět je: # - popisky os (včetně jednotek, ve kterých jsou vyneseny jednotlivé osy - v těchto příkladech bylo vynecháno...) # - dostatečně informativní legenda (včetně tvaru prokládané funkce) # - různé čáry a symboly musí být dostatečně odlišeny (barva, typ čáry) a popsány v legendě ### Na některých operačních systémech a verzích softwaru R a RStudio bývá někdy problém s češtinou a s rámečkem legendy při ### exportu obrázků. Tento problém se (většinou) dá vyřešit tím, že grafy kreslíme rovnou do souborů, což se provede následujícími příkazy. # Pro otevření kreslení grafů do souborů je příkaz pdf("...", width = ..., height = ..., encoding = "..."). # Místo prvního argumentu "..." patří jméno souboru (musí končit .pdf). Argumenty width a height upravují šířku a výšku výsledného obrázku # (dobrá volba je width = 7, height = 5). Argument encoding upravuje kódování textu v obrázcích (v závislosti na vašem operačním systému # je potřeba vyzkoušet možnosti "ISOLatin" a "ISOLatin2"). Soubor se vytvoří (nebo přepíše, pokud už existoval) v současném pracovním adresáři. # Současný pracovní adresář můžete zobrazit pomocí příkazu getwd() (getwd = get working directory) a změnit pomocí příkazu setwd(), který bere # jediný argument, kterým je cesta v adresářové struktuře vašeho počítače. getwd() # setwd("...") # Místo ... je případně potřeba doplnit cestu do požadovaného adresáře. pdf("...", width = 7, height = 5, encoding = "ISOLatin2") # Od příkazu pdf() se všechny grafické výstupy kreslí do příslušného souboru, nikoliv do grafického okna v RStudiu, jak jste zvyklí # (ani je tam tedy neuvidíte). # ZDE PATŘÍ PŘÍKAZY NA KRESLENÍ GRAFŮ A LEGENDY... dev.off() # Příkaz dev.off() přepíná zpět do režimu kreslení do grafického okna v RStudiu. # Pokud by byl s tímto nějaký problém, neváhejte napsat na nozicka@fykos.cz a nějak to společně vyřešíme. ########## Regresní diagnostika ########## # Nyní se budeme věnovat regresní diagnostice tak, jak byla popsána v textu seriálu #### Správnost prokládané funkce - grafické nástroje # Budeme uvažovat dva regresní modely pro naše data, jeden správný a jeden špatný, a ukážeme si, jak z residuí # můžeme poznat, že pracujeme se špatným modelem. dobryModel <- lm(y2 ~ x + I(x^2), data = mereni) spatnyModel <- lm(y2 ~ x, data = mereni) # Z odhadů regresních koeficientů nic nepoznáme, ale ve špatném modelu (tam kde prokládáme špatnou funkci) jsou # všechny odhady špatně. Jen z těchto informací ale nedokážeme poznat, který model je ten špatný. summary(dobryModel) summary(spatnyModel) # Zkusíme se podívat na graf naměřených dat a proložené funkce plotFit(dobryModel) plotFit(spatnyModel) # Zde je jasně vidět rozdíl. U špatného modelu je regresní funkce na krajích grafu výrazně pod měřenými hodnotami # a uprostřed grafu zase výrazně nad měřenými hodnotami (což značí špatnou volbu prokládané funkce). # U dobrého modelu prokládaná funkce krásně kopíruje měřená data. # Totéž bychom poznali z grafu residuí. Residua modelu jsou součástí výstupu funkce lm(), dostaneme se k nim pomocí # příkazu $residuals plot(mereni$x, dobryModel$residuals) plot(mereni$x, spatnyModel$residuals) # Na ose x jsou měření nezávisle proměnné a na ose y jsou residua modelu. U správného modelu můžeme vidět náhodně # rozeseté body kolem osy x. U spatneho modelu vidíme, že residua jsou výrazně kladná na krajích grafu a výrazně # záporná uprostřed grafu (což značí špatnou volbu prokládané funkce). # Totéž s použitím trochu záludnějších dat dobryModel <- lm(y3 ~ x + I(x^2), data = mereni) spatnyModel <- lm(y3 ~ x, data = mereni) # grafy proložené závislosti plotFit(dobryModel) plotFit(spatnyModel) # Zde toho moc není vidět. # Zkusíme grafy residuí plot(mereni$x, dobryModel$residuals) plot(mereni$x, spatnyModel$residuals) # Na grafech residuí už vidíme výrazný rozdíl mezi správným a špatným modelem. #### Statistiký test správnosti prokládané funkce (Lack of fit test) # Tento test můžeme použít pouze v případě, že pro každou hodnotu nezávisle proměnné máme naměřeno více (alespoň 5) # hodnot závisle proměnné. Naše data tedy musí vypadat nějak takto: plot(y4 ~ x2, data = mereni) # Opět budeme porovnávat 2 modely. dobryModel <- lm(y4 ~ I(x2) + I(x2^2), data = mereni) spatnyModel <- lm(y4 ~ x2, data = mereni) # Mohli bychom použít stejné grafické metody jako výše, ale v tomto speciálním případě můžeme provést i statistický # test. Podoba tohoto testu byla popsána v seriálu, proto ho zde už nebudeme opakovat. # Pro tento test bohužel neexistuje vhodná vestavěná funkce, zde ho provedeme pomocí více příkazů. plnyModel <- lm(y4 ~ as.factor(x2), data = mereni) anova(dobryModel, plnyModel) anova(spatnyModel, plnyModel) # Jen stručné vysvětlení: plnyModel je model, kde pro každou hodnotu nezávisle proměnné uvažujeme vlastní hodnotu závisle proměnné. # Tento model nám pomůže při určování různých součtů čtverců, které se vyskytují ve vzorce na testovou statistiku tohoto testu. # Příkaz anova (anova = analysis of variance) následně spočte všechny příslušné součty čtverců a testovou statistiku. # Ve sloupci Res.df můžeme vidět členy (n - p) a (N - n), ve sloupci RSS (RSS = residual sums of squares) můžeme vidět součty # čtverců vyskytujících se v testové statistice. Ve sloupci F je hodnota testové statistiky a v posledním sloupci je p-hodnota # testu (podle ní se orientujte). # Je vidět, že pro dobryModel testovou hypotézu nezamítáme a pro spatnyModel zamítáme. # Ke stejnému výsledku bychom dospěli i pomocí grafických metod plotFit(dobryModel) plotFit(spatnyModel) #### Homoscedasticita - grafické nástroje # Opět budeme porovnávat modely a ukážeme si, na co se zaměřit. dobryModel <- lm(y5 ~ x + I(x^2), data = mereni) OKModel <- lm(y6 ~ x + I(x^2), data = mereni) spatnyModel <- lm(y7 ~ x + I(x^2), data = mereni) # Porovnáme grafy proložených funkcí plotFit(dobryModel) plotFit(OKModel) plotFit(spatnyModel) # Na těchto grafech jsou vidět nějaké náznaky toho, že u některých modelů by nemusel být splněn předpoklad # homoscedasticity, ale lépe to bude vidět z grafu residuí. plot(mereni$x, dobryModel$residuals) # Tento model je bezproblémový. Residua jsou stejnoměrně rozprostřeny kolem osy x, předpoklad homoscedasticity # je splněn. plot(mereni$x, OKModel$residuals) # Zde už vidíme, že resiua odpovídající nízkým hodnotám nezávisle proměnné jsou menší než residua odpovídající # velkým hodnotám nezávisle proměnné. Nicméně tato závislost není nijak velká a takovýto model bude tedy poskytovat # poměrně přesné závěry. plot(mereni$x, spatnyModel$residuals) # Zde už vidíme, že residua silně závisí na hodnotě nezávisle proměnné. Takovýto model nebude poskytovat přesné # závěry, neboť předpoklad homoscedasticity je výrazně porušen. # Grafy residuí mohou na první pohled vypadat dost podobně (zejména poslední 2), ale musíme sledovat měřítko osy y! # Podobně bychom mohli zkonstruovat grafy hodnot proložené funkce oproti residuím. Měli bychom tam vidět podobnou # závislost. Hodnoty proložené funkce jsou uloženy ve výstupu z funkce lm() v proměnné $fitted.values plot(dobryModel$fitted.values, dobryModel$residuals) plot(OKModel$fitted.values, OKModel$residuals) plot(spatnyModel$fitted.values, spatnyModel$residuals) # Z těchto grafů dojdeme k podobným závěrům. Opět musíme sledovat měřítko osy y! # Je nutné poznamenat, že s porušeným předpokladem homoscedasticity toho moc nenaděláme (jediné řešení je přeměřit # data jinou metodou), můžeme jen v diskuzi upozornit na možné nepřesnosti, které mohou vzniknout. #### Grafický test nezávislých měření # Nezávislá měření se v praxi ověřují velice těžko, je nutné na toto dostatečně myslet při měření dat a zajistit, # aby byl tento předpoklad splněn. Jen jedna ukázka grafu residuí vůči posunutým residuím. dobryModel <- lm(y2 ~ x + I(x^2), data = mereni) spatnyModel <- lm(y2 ~ x, data = mereni) # residua modelů dobraRes <- dobryModel$residuals spatnaRes <- spatnyModel$residuals # Posunutá residua vytvoříme tak, že jednou přidáme nulu na začátek vektoru residuí a jednou na konec. Tímto na # odpovídajících si místech budou residua o jedno posunutá (první resp. poslední prvek vektoru je potom nula, # ale tuto malou nepřesnost můžeme ignorovat) dobraResPos <- c(dobraRes, 0) dobraResNepos <- c(0, dobraRes) spatnaResPos <- c(spatnaRes, 0) spatnaResNepos <- c(0, spatnaRes) # Grafy residuí vůči posunutým residuím plot(dobraResNepos, dobraResPos) abline(a = 0, b = 0, lty = 2) abline(a = 0, b = 1000000, lty = 2) plot(spatnaResNepos, spatnaResPos) abline(a = 0, b = 0, lty = 2) abline(a = 0, b = 1000000, lty = 2) # U dobrého modelu je vidět náhodný shluk bodů kolem počátku souřadnic. U špatného modelu je vidět koncentrace # bodů především v 1. a 3. kvadrantu. # Lepší je si do grafů vynést také osy odpovídající nulovým hodnotám x a y, což se provede pomocí funkce abline(a = ..., b = ...), # tato funkce do grafu nakreslí přímku o rovnici y(x) = a + b*x. Dodatečným argumentem lty můžeme upravit typ kreslené přímky. # Jak jste si určitě všimli, trochu podvádíme. Vertikální čára není úplně vertikální, ale je to přímka o rovnici y(x) = 100000 * x. # Na výsledném grafu ji ale běžný smrtelník nedokáže odlišit od vertikální čáry :D.