######################################## Zpracování dat fyzikálních měření - 6. 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("ukazka1.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 ukazka1.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 # pole napište jméno funkce) nebo zkuste najít odpoveď na internetu nebo napište autorovi seriálu na mail nozicka@fykos.cz ########## Instalace potřebných knihoven ########## install.packages("investr") # Instalovat stačí spustit jednou. install.packages("splines") # Instalovat stačí spustit jednou install.packages("sandwich") # Instalovat stačí spustit jednou library(investr) # Importování knihovny musíte provést pokaždé, když znovu zapínáte RStudio. library(splines) # Importování knihovny musíte provést pokaždé, když znovu zapínáte RStudio. library(sandwich) # Importování knihovny musíte provést pokaždé, když znovu zapínáte RStudio. ########## Nelineární regrese #################### # Do proměnné mereni je potřeba načíst data ze souboru ukazka1.csv #### Graf s naměřenými daty plot(mereni$x, mereni$y) ### Graf s naměřenými daty, legendou a popiskami os plot(mereni$x, mereni$y, xlab = "hodnoty x", ylab = "hodnoty y", main = "První graf") legend("topleft", legend = c("naměřené hodnoty"), pch = c(1)) # Tyto funkce byly podrobně popsány ve vzorovém skriptu k minulé sérii (zde uvádíme jen pro připomenutí). #### Odhadování parametrů v nelineárním regresním modelu # O odhadování parametrů se stará funkce nls() (nls = non-linear solver). model <- nls(y ~ log(a + b * x), data = mereni, start = list(a = 1, b = 1)) # První argument je tvar prokládané funkce. # - vždy se požaduje tvar: ~ # - mohou se zde vyskytovat jen nezávisle proměnná, regresní koeficienty (libovolně pojmenované, zde např. a, b) a vestavěné # matematické funkce (jako *, /, ^2, sqrt(), log(), sin(), cos(), exp(), atd.) # - v případě lineární regrese se vždy implicitně uvažoval absolutní člen, nyní už tomu tak není # 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é). # Třetí argument start = list(...) specifikuje počáteční body pro numerickou metodu hledání odhadů regreních koeficientů (gradient descent) # - proměnná start vždy požaduje jen a pouze objekt typu list, který se vytvoří příkazem list(...) # - uvnitř příkazu list(...) specifikujeme počáteční body pro všechny regresní koeficienty # - podle tohoto příkazu algoritmus mimo jiné pozné, co má chápat jako regresní koeficient a co nikoliv, proto je důležité # uvést opravdu všechny regresní koeficienty ### Grafické znázornění prokládané funkce (tento příkaz potřebuje knihovnu investr) plotFit(model) # Stejný příkaz jsme používali v minulém vzorovém skriptu, zde uvádíme jen pro zopakování. #### Pro shrnutí odhadů v našem nelineárním regresním modelu použijeme funkci summary() summary(model) # Nyní něco k interpretaci toho, co dostaneme na výstupu: # - Formula: popisuje, jakou funkci jsme naměřenými daty prokládali. # - Parameters: udávají odhady regresních koeficientů. Každý řádek odpovídá jednomu regresnímu koeficientu, # jména jsou stejná, jaká jsme uvedli v příkazu nls(). # - Ve sloupci Estimate je hodnota odhadu příslušného koeficentu. # - Ve sloupci Std. Error je nejistota měření tohoto koeficientu # - Sloupce t value a Pr(>|t|) udávají testovou statistiku (t value) a p-hodnotu (Pr(>|t|)) testu, že příslušný # regresní koeficient je roven 0 (pro jinou hypotézu a alternativu je potřeba spočítat testovou statistiku ručně). # - Řádek Signif. codes se stále vztahuje k testům o nenulovosti regresních koeficientů -> zvýrazňuje ty testy s malou p-hodnotou, # pokud je na konci řádku s regresním koeficientem symbol *** znamená to, že p-hodnota příslušného testu je velice malá (menší než 0.001), # pokud je na konci řádku s regresním koeficientem symbol ** znamená to, že p-hodnota příslušného testu je menší než 0.01 atd. # - 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ů. # - Number of iterations to converge udává počet kroků numerického algoritmu, který musel proběhnout. # - Achieved convergence tolerance udává přesnost, se kterou numerický algoritmus počítal. ########## Kreslení grafických výstupů #################### # Hezčí graf plotFit(model, xlab = "hodnoty x", ylab = "hodnoty y", main = "Ukázka nelineární regrese") # Všechny tyto funkce byly podrobně popsány v minulém dílu seriálu, jen stručně zopakujeme: # - xlab, ylab udávají popisky os # - main udává nadpis grafu # - pro připomenutí: pokud chceme psát matematické znaky použijeme příkaz expression(...) # např. expression(lambda) - napíše řecké písmeno lambda (podobně alpha, beta, gamma, atd.) # expression(x^2), expression(sqrt(x)) - pro matematické výrazy # pro spojení textu a matematických výrazů použijeme příkaz paste(..., ..., ...) # např. expression(paste("hodnoty ", lambda, " [", m^2, "]")) #### Intervalový odhad pro prokládanou funkci # Intervalový odhad pro prokládanou funkci vykreslený do stejného grafu plotFit(model, xlab = "hodnoty x", ylab = "hodnoty y", main = "Ukázka nelineární regrese", interval = "confidence", level = 0.95) # Intervalový odhad se do grafu přidá argumentem interval = "confidence" # Argument level je hladina spolehlivosti tohoto intervalu (pokud není zadáno, bere se hodnota 0.95) # To stejné za použití jiných barev a stylů čar plotFit(model, xlab = "hodnoty x", ylab = "hodnoty y", main = "Ukázka nelineární regrese", interval = "confidence", level = 0.95, col.fit = "red", col.conf = "blue", lty.fit = 1, lty.conf = 3, pch = 4) # Argument col.fit (resp. col.conf) udává barvu proložené funkce (resp. intervalu spolehlivosti). # Argument lty.fit (resp. lty.conf) udává typ čáry proložené funkce (resp. intervalového odhadu). # Argument pch udává, jakými symboly se mají znázornit měřená data. #### Přidávání legendy legend("topleft", legend = c(expression(paste(y, " = ", log(a + bx))), "naměřená data", "proložená funkce", "95 % interval spol."), lty = c(NA, NA, 1, 3), pch = c(NA, 4, NA, NA), col = c(NA, "black", "red", "blue")) # První argument je poloha legendy v grafu # Argument legend = c(..., ..., ...) udává popisky legendy (může se uvnitř používat příkaz expression() pro matematické symboly) # Argument lty = c(..., ..., ...) udává typy čar vykreslených na jednotlivých řádcích legendy (NA = žádná čára na příslušném řádku) # Argument pch = c(..., ..., ...) udává typy symbolů vykreslených na jednotlivých řádcích legendy (NA = žádný symbol na příslušném řádku) # Argument col = c(..., ..., ...) udává barvy jednotlivých symbolů a čar na příslušných řádncích legendy ### 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 #################### # Jak jsme psali v textu seriálu, regresní diagnostika se v případě nelineární regrese dělá naprosto stejně, jako v případě regrese # lineární. Je založena na residuích modelu a jejich vlastnostech. Residua se získají pomocí příkazu residuals(...) residua <- residuals(model) # Základním nástrojem je graf naměřených dat a proložené funkce. Pokročilejší metody využívají různé grafy residuí. ### Graf residuí oproti nezávisle proměnné plot(mereni$x, residua) # Měli bychom vidět náhodný shluk bodů kolem osy x. Pokud ne, je něco špatně. V tomto případě je všechno v pořádku. ### Graf residuí oproti hodnotám prokládané funkce prolozene <- fitted(model) plot(prolozene, residua) # Měli bychom vidět náhodný shluk bodů kolem osy x. Pokud ne, je něco špatně. V tomto případě je všechno v pořádku. ### Graf residuí oproti posunutým residuím posunuta <- c(0, residua) neposunuta <- c(residua, 0) plot(neposunuta, posunuta) abline(a = 0, b = 0, lty = 2) abline(a = 0, b = 100000, lty = 2) # Měli bychom vidět náhodný shluk bodů kolem počátku souřadnic bez tendence shlukování v jednotlivých kvadrantech. # V tomto případě je vše v pořádku. # Ke konstrukci tohoto grafu jsme opět použili drobný trik (první a poslední bod je nepřesný kvůli nulové hodnotě), # celkový výsledek tímto ale nebude nijak zásadně ovlivněn. ########## Whiteův (sandvičový) odhad kovarianční matice #################### #### Načtení nových dat pro tento úkol (načtěte si data ze souboru ukazka3.csv) mereni <- read.csv("ukazka3.csv", sep = ",") ### Grafické znázornění dat plot(mereni$x, mereni$y) ### Proložení paraboly klasickým způsobem model <- lm(y ~ x + I(x^2), data = mereni) plotFit(model) ### Ověření homoskedasticity dat # Vidíme už z tohoto grafu, že data jsou pravděpodobně heteroskedasticitní. # Graf residuí oproti hodnotám nezávisle proměnné plot(mereni$x, residuals(model)) # Nyní už jednoznačně vidíme heteroskedasticitu. ### Whiteův (sendvičový) odhad kovarianční matice # Kovarianční matice odhadnutá klasickým způsobem cov <- vcov(model) cov # Odmocniny z diagonálních prvků této matice jsou nejistoty měření regresních koeficientů, # můžeme se o tom jednoduše přesvědčit. # Diagonální prvky diag <- diag(cov) # Odmocniny z diagonálních prvků sqrt(diag) # Srovnejme s std. error, který poskytne klasický regresní model. summary(model) # Musíme si na tomto místě uvědomit, že pokud je porušena homoskedasticita, jsou tyto nejistoty získané # klasickou cestou špatně a nemůžeme je používat! # Musíme použít jinou kovarianční matici získanou pomocí sendvičového odhadu covTRUE <- sandwich(model) covTRUE # Odmocniny z diagonálních prvků této matice jsou správné nejistoty měření regresních koeficientů. diagTRUE <- diag(covTRUE) # Porovnejte orpavdové nejistoty měření a nejistoty měření získané špatnou metodou. # Špatné sqrt(diag) # Správné sqrt(diagTRUE) # Vidíme, že moc velký rozdíl v nejistotách měření nepozorujeme, i když heteroskedasticita byla opravdu # výrazná. Toto ospravedlňuje zanedbání heteroskedasticity, pokud je jen malá (obdrželi bychom prakticky # stejné nejistoty měření regresních koeficientů oběma metodami). # Všechny funkce v tomto odstavci fungují stejně pro lineární i nelineární regresní modely (model můžeme # získat příkazem "model <- lm(...)" nebo "model <- nls(...)", Whiteův odhad získáme stejným způsobem) ### Další příklad - data z minulé seriálové úlohy # V minulé seriálové úloze jsme odhalili (alespoň to byl náš úkol) mírnou heteroskedasticitu # v datech z příkladu c). V diskuzi stačilo zmínit, že ji můžeme zanedbat a výsledky se téměř nezmění, # nyní si to můžeme ověřit. # Načteme data linreg2.csv (dostupné u zadání seriálové úlohy v 5. sérii) mereni <- read.csv("linreg2.csv", sep = ",") # Základní grafické znázornění dat plot(mereni$x, mereni$y) # Použitý model model <- lm(y ~ I(log(x)), data = mereni) summary(model) plotFit(model) # Nejistoty měření při zanedbání heteroskedasticity sqrt(diag(vcov(model))) # Nejistoty měření pomocí Whiteova odhadu kovarianční matice sqrt(diag(sandwich(model))) # Prakticky žádný rozdíl u koeficientu odpovídajícího log(x). # Jen malý rozdíl u Interceptu (tento koeficient bude heteroskedasticitou vždy ovlivněn nejvíce a stejně # se příslušná nejistota změnila jen mírně). ########## Regresní spliny ###################### # Do promenne mereni je potreba nacist data ze souboru ukazka2.csv mereni <- read.csv("ukazka2.csv", sep = ",") # Grafické znázornění dat plot(y ~ x, data = mereni) # Napadá Vás nějaká "klasická" funkce, která by se tímto dala proložit? # Odpověď je ne, musíme použít regresní spliny. #### Proložení dat regresními spliny model <- lm(y ~ bs(x, knots = c(1, 1.8, 2.2, 3.4, 3.6, 4.5, 5), degree = 3), data = mereni) # K proložení dat regresními spliny využíváme lineární regresi (příkaz lm()), nikoliv nelineární regresi! # Zavoláme funkci lm() tak, jak jsme zvyklí, akorát první argument (prokládaná funkce) musí být jiný než obvykle. # - Regresní spliny budeme prokládat, pokud vzorec prokládané funkce bude mít tvar y ~ bs(...). Příkaz bs() zavolá # regresní spliny (bs = basis spline) # - První argument příkazu bs() je nezávisle proměnná v našem modelu # - Druhým argumentem (knots) je seznam uzlů (navíc se automaticky vytvoří uzel v minimální a maximální hodnotě nezávisle proměnné) # - Třetím argumentem (degree) příkazu bs() je stupeň polynomů, které chceme prokládat. # Grafické znázornění proložené funkce plotFit(model) # Stupeň polynomu většinou volíme jako 3, vyšší stupeň už nepřinese výrazně lepší výsledky. Můžete si sami vyzkoušet vliv tohoto parametru. # Uzly volíme jako body, ve kterých vidíme změnu chování měřených hodnot. Volbu uzlů zkoušíme metodou pokus-omyl, dokud nebude graf # proložené funkce hezky vypadat (někdy je nutné dát 2 uzly hned vedle sebe - zejména v bodech, kde se výrazně mění chování měřených hodnot). # Musíme si dát pozor, abychom nezvolili příliš velký počet uzlů! Pro dobré fungování popsaných metod je ideální mít alespoň 5 krát # více měření než uzlů (čím méně uzlů, tím lépe, ale nesmí to být na úkor správnosti prokládané funkce). #### Odhady regresních koeficientů summary(model) # Je nutné si uvědomit, že tyto odhady nemají žádnou fyzikální interpretaci. # Jediné, čeho docílíme je hezky vyadající graf! #### Graf se všemi náležitostmi plotFit(model, xlab = "hodnoty x", ylab = "hodnoty y", main = "Ukázka regresních splinů", interval = "confidence") legend("topright", legend = c("regresní spliny stupně 3", "naměřené hodnoty", "proložená funkce", "95 % interval spol."), lty = c(NA, NA, 1, 2), pch = c(NA, 1, NA, NA)) ########## Chybové úsečky ########################## # Na tomto místě jen stručně ukažme, jak se nakreslí graf s chybovými úsečkami. #### Generování náhodných dat (pro použití v grafu) x <- c(1:10) y <- 2 * x + rnorm(10, mean = 0, sd = 0.5) stdev <- runif(10, min = 0.5, max = 2) #### Kreslení grafu # předpokládejme, že v proměnné x máme hodnoty nezávisle proměnné, v proměnné y průměry měření závisle proměnné odpovídající hodnotám # nezávisle proměnné a v proměnné stdev máme hodnoty výběrových směrodatných odchylek průměrů měření nezávisle proměnné odpovídajících # hodnotám závisle proměnné. # V tomto příkladě jsou tyto hodnoty náhodně vygenerované, v praxi je musíte spočítat z naměřených dat! (mohou vám v tom pomoci # vestavěné funkce mean() a sd(), které jsme popsali v 2. a 3. dílu seriálu) # Odhad lineárního regresního modelu (prokládáme přímku) model <- lm(y ~ x, data = data.frame(x, y)) # V případě nelineární regrese bychom zavolali funkci nls() a dále bychom pokračovali úplně stejně. # Vykreslení grafu proložené funkce bez naměřených hodnot. plotFit(model, pch = NA, xlab = "hodnoty x", ylab = "hodnoty y", main = "Ukázka grafu s chybovými úsečkami") # Nechceme do grafu kreslit přímo naměřené hodnoty (to zařídíme pomocí argumentu pch, který musí mít hodnotu NA), protože pozdějí # do grafu vyneseme chybové úsečky. # Přidáme body odpovídající průměrům hodnot nezávisle proměnné. points(x, y) ### Přidáme chybové úsečky. arrows(x, y - stdev, x, y + stdev, angle = 90, code = 3, length = 0.05) # Chybové úsečky se přidávají pomocí příkazu arrows (tento příkaz původně sloužil na kreslení šipek, ale vhodnou volbou argumentů pomocí # něj vykreslíme chybové úsečky). # - první 4 argumenty popisují konce chybových úseček (po řadě x-ová a y-ová souřadnice začátku úsečky a konce úsečky). Tento argument bere # buď jednotlivá čísle nebo celý vektor čísel, jako v této ukázce (potom vykreslí více chybvých úseček - jedna chybová úsečka bude odpovídat # prvkům na stejném místě ve 4 zadaných vektorech) # - argument angle ovlivňuje jak budou vypadat konce chybových úseček, přesněji řečeno úhel mezi "tělem úsečky" a "koncovou čárou" # (pro hodnoty menší než 90 to budou spíše šipky - můžete si vyzkoušet) # - argument code ovlivňuje styl chybové úsečky (hodnota 3 udává stejné zakončení na obou stranách, vždy používejte tento styl). # - Argument length udává dělku zakončení chybových úseček (těch malých vodorovných čar na koncích úseček), tento argument nemá žádný vliv # na velikost samotné úsečky. # Pokud bychom chtěli přidat i vodorovné chybové úsečky, postupovali bychom analogicky stdevX <- runif(10, min = 0.15, max = 0.4) # Náhodně vygenerované délky chybových úseček, v praxi se tyto hodnoty musí vypočítat z měřených dat! arrows(x - stdevX, y, x + stdevX, y, angle = 90, code = 3, length = 0.05) # Musíme si jen uvědomit, jaké budou nyní souřadnice začátků a konců chybových úseček.