Metody v klimatologii III. Využití R v klimatologii Úvod do R https://crQn.r-proiecT.orq/ • Interaktivní „výpočetní prostředí" pro statistické výpočty a grafiku • Volně šiřitelný pod SNU GPL • Nezávislý na platformě • Řádkový interpret, objektově orientovaný • Umožňuje jednoduché skriptování • Snadno rozšiřitelný ( > 2000 rozšiřujících knihoven na WWW) • Nepřeberné množství výukových materiálů • Počáteční investice do ovládání jazyka R • Nutnost často upgradovat R i některé knihovny Úvod do R Úvod do R Grafické uživatelské rozhraní - R studio, Rand «JdUmOrnJPTHingH- . '■■ IToriMlPIM.il ■ .dMWdl« a' war* i cite History 3 E] Sown on Sure T 3( i» _*Sourct ■ |_ lata- y un- import oaiiMt> /cwriui 1 HBrary^ggploti) 2 sourceC'plots/f oriMXPlot. R") 1 ^ 5 siMMry( sumniaryfdianiaridsSprice) win. lit qu. tsedian Mean 3rd Qu. Max. 326 930 2401 393 3 5324 1S620 > avesije <- round(niean(di amends Scarat), 4) > clarity <- levels [diamonds JelaHty) > p *- qplctCcarat, price, + data-diamonds, color-clarity, + xlah-"carat", /lali-",Priee''l + main-"Oismgnd Pricing") •J > foriwt. piotcp. sire-M) ■"i Carat R v klimatologii 1) Speciální knihovny (libraries) dpIR, treeclim - dendroklimatologie ncdf4 - analýza polí ve formátu NetCĎF climatol 2) www stránky https://www.benjaminbell.co.uk/2018/01/getting-climate-data.html https://rclimate.wordpress.com/ imate Using R andData to Understand dibits Change Základy práce v R S využitím dodaných skriptů obsahujících vybrané příkazy jazyka R a datových souborů si projděte základy práce v programovém prostředí R a zopakujte si základní úlohy popisné statistiky • Nainstalujte program R: https://cran.r-pro iect.org/ • Vytvořte si pracovní adresář, např. C:/data/AAF(5cviceni ■ Ďo tohoto adresáře si z ISu zkopírujte všechny soubory, které najdete ve složce R_cviceni • Spusťte program R a otevřete si úvodní ukázkový skript: File -> Open scrípt... scr_l_zaklad_2020. R • Na 4. řádku upravte cestu k vašemu pracovnímu adresáři : setwd("C:/data/AAF(5cviceni") 3 Základy práce v R S RGui 02-bit] File Edit Packages Windows Help [^lÍHlP[nl W í R Console rc.ean of y. rc.ean of y Q..65317204 0.03322581 > f statistické testy: paiovy t-test (porovnáni dvou nezávislých souboru) > if t. test [nao_index[, ľ], nao index|, 8], paired — T] I > > £ statistické testy: ±-test (lisi se významne variabilita hodnot indexu NAC v? :> var . test [~ao_lr.dex .", 2" , r.ao index1.31) ssrc.ple estirc.ate, ratio of variant 1. 273' boxplot [nao it R C:\dataVR_data\MFGcviceni\scr_1_zaklad_2G2G.R- R Editor t Uvod dc S; verze 2020/10 data: nao_indeJ nUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU 1.2737_ njir ^ nastaveni pracovního adresáre alternativě hyp< ■•'"■^ ' "W»™m"*""j*"^'IWW" ' 95 percent conf: getwd[) # výpise pracovni adresár C.9539535 1.70_ # DATA, DATOVÉ STRUKTURY # vytvor! prorr.énnou a zápise do ni hodnoty x. <- c (1,2,3,4,5,6,7,8) sršíš rc.eaľi [xj length(x) # nápověda : ?ir.ean () is.vector[x) # vektor je základni datovou strukturou v R, konstanta je vektor delky 1 Základy práce v R • Klikněte na začátek prvního řádku. Postupně tiskněte CTRL-R. Program bude vykonávat jednotlivé příkazy kódu • Sledujte výstupy příkazů na konzole a později také v grafickém okně • Příkazy můžete provádět opakovaně v závislosti, kde umístíte kurzor či jakou část kódu vyberete • U jednotlivých příkazů můžete měnit nastavení parametrů či dodávat parametry další s pomocí helpu: ? nazev_f unkce() • Projdete-li všechny příkazy, můžete vyzkoušet script scr_2_graf .R, který vytvoří přiložený graf. Část jeho kódu využijete ve cvičení 1800 1850 1900 1950 2000 Roky - chladné • Explorační (průzkumová) analýza s cílem otestovat vlastnosti časové řady (homogenizace, autokorelace, přítomnost odlehlých hodnot, normalita, homoskedascita, ...) ■ Analýza trendové složky ■ Analýza cyklického chování • Modelování časových řad (např. ARIMA modely), • Analýza extrémních hodnot (viz část IV) • Hledání změn v chování časové řady (change point detection) - nejen jako součást procesu homogenizace, ale též jako indikátor změny režimu chování (structural changes) • Superposed Epoch Analysis,... Analýza trendu Lineárni trend - MNC Smernice: 0.26 A /\ —i-1— 1990 2000 ttSrok Neměli bychom si zkontrolovat, zda je MNČ vůbec použitelná? Vypočtený trend nevystihuje dobře kolísání v celé délce časové řady Histogram of ts.tt 0 10 20 30 40 5 Neparametrický odhad lineárního trendu Může existovat řada důvodů, proč je MNC k odhadu lineárního trendu nevhodná 12 ř-6 š 2A 1935 IMn 1595 2M0 20nS Year /"uiw Model/'tuhí sen Mooei https:/ /www. researchgate.net/figure/Comparis on-of-linear-and-Theil-Sen-regression-slopes-for-a-pixel-in-the-Sahel_fig3_266705601 • Hledáme metody, která je využitelná z hlediska vlastností analyzovaných dat • Nezávisí na rozdělení hodnot • Není ovlivněna přítomností odlehlých hodnot, Neparametrický odhad lineárního trendu • Nejprve testujeme, zda řada obsahuje významný trend Mann- Kendallův test • Následně odhadujeme parametry přímky - Theill-Senová metoda odhadu parametrů • M-K test je testem na přítomnost monotónního trendu v časové řadě • Lze ho využít, pokud nejsou splněny podmínky pro použití MNC • Neparametrický odhad - nevyžaduje splnění předpokladu normality • Není citlivý k odlehlým hodnotám • Lze ho aplikovat i na časovou řadu s chybějícími hodnotami • © Vyžaduje, aby hodnoty členů časové řady byly nekorelované 6 Neparametrický odhad lineárního trendu Mann-Kendallův test analyzuje rozdíly ve znaménkách mezi dřívějšími a pozdějšími datovými body. Pokud je v řadě přítomen trend, hodnoty znamének budou mít tendenci neustále růst resp. neustále klesat. The MK test tests whether to reject the null hypothesis (ifo) and accept the alternative hypothesis (ff0), where fro: No monotonie trend Ha : Monotonie trend is present The initial assumption of the MK test is that the Hq i& true and that the data must be convincing beyond a reasonable doubt before Hq is rejected and Ha is accepted The MK test is conducted as follows (from Gilbert 1987, pp. 209-213) : 1. List the data in the order in which they were collected over tnne.iii, ij... .. xn, which denote the measurements obtained at times 1, 2... . tn, respectively 2. Determine the sign of alliifji — l)/2 possible differences Xj — x^._ where j > k. These differences are X% — Xt ,X3 — XI f . . . , X„ — X1,X5 — X2,X& — X$, . . . , X„ — Xn-57Xn — Xn_l 3. Let SgIl(aTj — x^) be an indicator function that takes on the values 1. 0= or -1 according to the sign of xj — Xj,., that is: sgn(itj — jiit) - 1 if Xj — Xf. > 0 = 0 ifxj — Xfc = 0= or if the sign of x j — x^. determined due to non-detects 4. Compute -1 if ar3- - xk < 0 which is the number of positive differences minus the number of negative differences. If is apositive number, observations obtained later in time tend to be larger than observations made earlier. If is a negative number, then observations made later in time tend to be smaller than observations made earlier. Neparametrický odhad lineárního trendu Mann-Kendallův test 5. If n < 10, follow the procedure described in Gilbert (1987, page 209, Section 16.4.1) by looking up S in a table of probabilities (Gilbert 1987, Table A18, page 272). If this probability is less than a (the probability of concluding a trend exists when there is none), then reject the null hypothesis and conclude the trend exists. If n cannot be found in the table of probabilities (which can happen if there are tied data values), the next value farther from zero in the table is used. For example, if S = 12 and there is no value for S = 12 in the table, it is handled the same as S = 13. If n > 10, continue with steps 6 through 10 to determine whether a trend exists. This follows the procedure described in Gilbert (1987, page 211, Section 16.4.2). 6. Compute the variance of S as follows: VAH(S) = [„(„ - 1)(2„+ 5) - £ *„(', " 1)(2*p + 5)] where g is the number of tied groups and tp is the number of observations in the p th group. 7. Compute the MK test statistic, 5 as follows: ZMK = —1 if-? > 0 MK -jVARiS) = 0 if S = 0 A positive (negative) value of Z,\jk indicates that the data tend to increase (decrease) with time. Neparametrický odhad lineárního trendu Mann-Kendallův test 8. Suppose we want to te&t the null hypothesis J/o: No monotonie trend versus the alternative hypothesis Ha: Upward monotonie trend at the Type I error rate a , where 0 < ťt < D.5_ (Kote that a is the tolerable probability that the MK te&t will falsely reject the null hypothesis.) Then ffo is rejected and Ha is accepted if Z\jk > Zi_a, where Z^_a is the 100(1 — a)131 percentile of the standard normal distribution. These percentiles are provided in many statistics book (for example Gilbert 19S7. Table Al_ page 254) and statistical software packages. 9. To test Hq above versus Ha: Downward monotonie trend at the Type I error rate a, H$ is rejected and Ha is accepted if Zjjjt < — Z^_a. 10. To test the H$ above versus H a : Upward oi downward monotonie trend atthe Type I errorrate a, Ha is rejected and Ha is accepted if \Zmk | ^ Z\_Qj% where the vertical bars denote ab&olute value. Počítačové programy běžně poskytují p-hodnoty příslušející vypočtenému testovacímu kritériu, pode které lze výsledek testu jednoznačně interpretovat Mann-Kendallův test v R # vyžaduj e knihovnu if(!require(Kendall)){install.packages("Kendall")} library(Kendall) HKD <- HannKendall(tt$Tx30! summary(HKD) Score = 652 , Var(Score) = 24489.33 denominator = 1738.724 tau = 0.375, 2-sided pvalue =3.1829e-05 V naší v časové radě se vyskytuje rostoucí statisticky významný monotónní trend. Jeho parametry zjistíme a trend vykreslíme Theil-Senovou metodou Metoda odhadu parametru lineárního trendu (Theil-Sen estimate) Opět se jedná o neparametrickou metodu Line Segments Connecting Pairs of Points 0 5 10 15 X 1) Spojíme úsečkami všechny body časové řady 2) Pro každou úsečku vypočteme její směrnici 3) Směrnici časové řady vypočteme jako medián všech směrnic jednotlivých úseček 4) Stejným způsobem vypočteme i absolutní člen jako druhý parametr přímky Různé metody odhadu trendu Běžný způsob odhadu lineárního trendu (LRM) může být ovlivněn výskytem odlehlých hodnot. ■ V takovém případě je vhodnější odhad lineárního trendu neparametrickými metodami (např. M-Kči QR) Adaptivní trend lze sestavit pomocí lokální regrese či pomocí splínových funkcí Krátkodobé tendence v časové řadě lze vystihnout jejím :-shlazením klouzavými průměry, (Saussovým filtrem či řadou dalších metod. CASOVA RADA ADAPTIVNÍ TREND leploly_brnolYEAR LINERARNI TREND 1350 1900 1950 teploty_brnc*YEAR SHLAZOVANI 1350 1900 1950 1eplo1y_bmoíYEAR V programu R byl sestaven skript pro aplikaci různých metod analýzy časové řady (scr_3_casove_rady.R), vstupní data obsahuje soubor brno_t.csv 9 Príklady analýzy - cyklická složka (wavelet analysis) Odhad spektrální hustoty - předpokládá, že zastoupení Vínková (wavelet) transformace cyklu v čase se nemění Príklady analýzy - cyklická složka (wavelet analysis) „Materská" vlnka provádí po časové řadě dva druhy pohybu: • posun • školování (roztažení a smrštění) https ://www.youtube .com/watch ?v=QX1-xGVFqmw Řada šířky letokruhů dubu, ve které se v některých obdobích objevuje statisticky významný cyklus délky 16 roků Příklady analýzy - cyklická složka (wavelet analysis) Temperature _ Precipitation _SPEI TSCO« leůO 1700 1800 1900 JOOöTsOOh 36cö 1700 1800 1900 ZOOoTsOOtE 1«» 1700 1800 19g0 ÍO0O Standardized cross-wavelet power spectrum 9.8e-04 3-te-02 l.Oe+00 3.2e+01 l.Oe+03 „Cross-wavelet" spektrum mezi řadami teploty, srážek a SPEI a řadami vysvětlujících proměnných (sluneční aktivita a NAO index). Šipky indikují fázový posun mezi dvěma řadami (ukazuje-li šipka vpřed, jsou obě řady ve fázi) 11 Príklady analýzy - „regime shift detection' Kolář et al. (2022) Effects of social and climatic factors on building activity in the Czech lands between 1450 and 1950 Calendar Year (AD) (A) Variability of standardized felling date indices in the period 1450-1950 completed with the model of regime shifts indicating individual breakpoints using Sequential T-test Analysis of Regime Shifts (bold line) and Pettitt test (arrows). The most serious wars, epidemics and famines are indicated. (B) Reconstruction of annual temperature (Dobrovolný et al. 2010) and precipitation (Dobrovolný et al. 2015) anomalies (w.r.t. 1961-1990) smoothed by Gaussian low pass filter and Pearson's correlation coefficients between these climate parameters and standardized number of felling dates for the common period 1501-1950. (C) Number of recorded local windstorms with highlighted widespread windstorms (black dots) based on Brázdil et al. (2004). (Ď) Evolution of population and forest cover in the Czech lands (CSO). Note: WAS - War of the Austrian Succession, 7YW- Seven Years' War._ Příklad analýzy prostorové informace Jaký je rozdíl mezi tlakovým polem v měsících, kdy je ve střední Evropě výrazné sucho (obrázek vlevo), v porovnání s normálními podmínkami (průměrné tlakové pole daného měsíce v období 1961-1990 - obrázek uprostřed). Vpravo je výsledek v podobě rozdílu (normál - sucho). Vyznačeny jsou gridové body, ve kterých je rozdíl statisticky významný May, mean SLP 1961-1990 mean SLP in 31 dry Mays since 1850 SLP anomaly (map2 - mapl) SLP at extremely dry Mays defined in the 1850-2010 period and SLP differences in extremely dry months compared to MSLP of the 1961-1990 reference period; black points mark statistically significant (a=0.05) SLP decrease/incerese w.r.t the 1961-1990 reference period V programu R byl sestaven skript pro vizualizaci sestavených polí (scr_4_slp_mapy.R), vstupní data obsahuje soubor au_sucho.csv