### Geostatistika ### ##################### library(geoR) # 1. nacteni dat a jejich pocatecni pruzkum: # namerene vysky v 52 mistech ctvercoveho okna o strane 6.7 jednotek # (1 jednotka = 50 stop = 15,24 metru) data(elevation) ?elevation # vysky v rozmezi 690 az 960 jednotek (1 jednotka = 10 stop = 3,05 metru): points(elevation,cex.min=1,cex.max=4) summary(elevation) # stredni pozorovana vyska je 827,1 plot(elevation,lowess=TRUE) # a) rozdeleni vysek do ctyr skupin (kvantily) # b) prubeh dat v zavislosti na menicich se souradnicich # (moznost odhalit prostorovy trend) # c) histogram dat - mirna asymetrie, gaussovsky model muze dobre aproximovat data; # obecne se jedna o histogram zkonstruovany z korelovanych pozorovani => mene stabilni # nez u nahodneho vyberu # Data povazujeme za realizaci nahodneho pole Z(x) v danych bodech x_1,...,x_n. # Z(x) znaci skutecnou vysku v bode x. # 2. pruzkumna analyza: # Z obrazku je videt prevaha vyssich hodnot v jizni casti oblasti => mozny prostorovy trend: # Z(x) = mu(x) + e(x). # Rovnez pozorujeme vyssi hodnoty koncentrovane na zapadni a vychodni strane, zatimco nizsi uprostred. ?points.geodata # Rezidua po nafitovani linearniho trendu: mu(x) = beta_0+beta_1x_1+beta_2x_2. points(elevation, trend="1st", abs=T, col="gray") # Rezidua po nafitovani kvadratickeho trendu: # mu(x) = beta_0 + beta_1x_1 + beta_2x_2 + beta_3x_1^2 + beta_4x_2^4 + beta_5x_1x_2. points(elevation, trend="2nd", abs=T, col="gray") # Kladna (cerne) a zaporna (bile) rezidua se vyskytuji u sebe - indikuje prostorove autokorelace. # Porovnani kvality fitovaneho trendu neni jednoduche: xmat1 <- unclass(trend.spatial(trend = "1st", geodata = elevation)) data1 <- as.vector(lm(elevation$data ~ xmat1 + 0)$residuals) summary(data1) xmat2 <- unclass(trend.spatial(trend = "2nd", geodata = elevation)) data2 <- as.vector(lm(elevation$data ~ xmat2 + 0)$residuals) summary(data2) par(mfrow=c(1,2), mar=c(1,2,1,1)) boxplot(data1, ylim=c(-70,120), main="Linearni trend") boxplot(data2, ylim=c(-70,120), main="Kvadraticky trend") dev.off() # Vykreslime rezidua: plot(elevation,lowess=TRUE,trend="1st") # Vyrazna zavislost rezidui na souradnici x. plot(elevation,lowess=TRUE,trend="2nd") # Uz neni jasny prostorovy trend v hodnotach rezidui. # Pro posouzeni prostorovych autokorelaci odhadneme variogram. # Vypocet (semi)variogramu je zalozen na hodnotach (Z(x_a) - Z(x_b))^2 / 2, zde je vykreslime proti # hodnotam vzdalenosti |x_a - x_b|: plot(variog(elevation, option="cloud")) # Oblak neni moc prehledny, ale da se nekdy vyuzit k identifikaci odlehlych pozorovani v tom smyslu, # ze pozorovani v danem bode bude hodnotou "od vsech ostatnich daleko" bez ohledu na prostorovou # (euklidovskou) vzdalenost mist pozorovani. # Variogram pro puvodni data: plot(variog(elevation, uvec=seq(0,5,by=0.5)), type="b") # Rezidua z modelu s linearnim trendem: res1.v <- variog(elevation, trend="1st", uvec=seq(0,5,by=0.5)) plot(res1.v, type="b") # Rezidua z modelu s kvadratickym trendem: res2.v <- variog(elevation, trend="2nd", uvec=seq(0,5,by=0.5)) lines(res2.v, type="b", lty=2) # Pro puvodni data variogram roste => pokud bychom modelovali stacionarnim nahodnym polem, # pak rozsah autokorelaci presahuje velikost studovane oblasti. # Proto je lepsi zahrnout prostorove se menici stredni hodnotu do modelu. # Pokud by empiricky variogram vykazoval pouze male prostorove korelace, muzeme se ptat, # jestli lze data modelovat jako Z(x_i) = mu(x_i) + e_i, kde e_i jsou nekorelovane. # Spocteme Monte Carlo obalky tak, ze provedeme 99 nahodnych permutaci rezidui napric polohami: mc1 <- variog.mc.env(elevation, obj=res1.v) plot(res1.v, env=mc1) mc2 <- variog.mc.env(elevation, obj=res2.v) plot(res2.v, env=mc2) # Grafy potvrzuji pritomnost prostorovych autokorelaci - rostouci tvar empirickeho variogramu # vypada jako vyznamny. Statistickou vyznamnost bychom vsak museli testovat na zaklade vetsiho # mnozstvi permutaci. # 3. parametricky model: # Predpokladejme, ze Z(x) je stacionarni gaussovske nahodne pole se stredni hodnotou beta # a Maternovou kovariancni funkci s parametrem kappa=1.5. # Model pro kovariancni funkci pole Z: Z(x) = Y*beta + S(x) + e, kde # Z(x) jsou pozorovane hodnoty, Y kovariaty, Y*beta je trend, S(x) je stacionarni gaussovske nahodne # pole, e je chybovy clen (napr. chyba mereni nezavisla na hodnotach pole, opakovana mereni ve stejnem # miste mohou dat ruzne hodnoty). # beta: vektor parametru popisujici trend, v nejjednodussim pripade primo (konstantni) str. hodnota # sigma^2: rozptyl hodnot pole S(x) (partial sill) # tau^2: rozptyl chyboveho clenu (nugget variance) # phi: skalovaci parametr kovariancni funkce S(x) # Odhady dalsich kovariancnich parametru (tau^2, sigma^2, phi) metodou maximalni verohodnosti jsou: ml0 <- likfit(elevation, ini=c(3000,2), cov.model="matern", kappa=1.5) ml0 # (pocatecni hodnoty pro numerickou maximalizaci byly tau^2=0, sigma^2=3000, phi=2) # Predchozi analyza ukazala, ze vyhodnejsi bude zahrnout prostorovy trend do modelu ml1 <- likfit(elevation, trend="1st", ini=c(1300,2), cov.model="matern", kappa=1.5) ml1 # Protoze trend vysvetlil cast prostorove variace, je odhad sigma^2 daleko mensi # nez u stacionarniho modelu. Zkusme jit jeste dale: ml2 <- likfit(elevation, trend="2nd", ini=c(1300,2), cov.model="matern", kappa=1.5) ml2 ################################ # # Ukol: provedte podobnou analyzu pro data s obsahem vapniku ve vzorcich pudy data(ca20) # Zkuste uvazovat trend obsahujici kovariaty area (druh oblasti - typ pudy) # a altitude (nadmorska vyska), napr. trend=~altitude nebo trend=~area+altitude apod. # Zahrnuli byste trend do modelu? Jaky? # Existuji prukazne prostorove autokorelace? plot(ca20) ca20 # Rozdeleni zkoumane oblasti na tri casti podle typu terenu: plot(ca20$borders, type="l") points(ca20$reg1, type="l") points(ca20$reg2, type="l") points(ca20$reg3, type="l")