######################### Pacote "ramps" ################################ # Version: 0.6-10 # # Date: 2011-08-16 # # Author: Brian J. Smith , # # Jun Yan , e # # Mary Kathryn Cowles # # # # Description: Bayesian geostatistical modeling of Gaussian processes # # using a reparameterized and marginalized posterior sampling (RAMPS) # # algorithm designed to lower autocorrelation in MCMC samples. # # Package performance is tuned for large spatial datasets. # ######################## ################################ ### Spatial Correlation Structur (estruturas de correlação espacial): ### # corRCauchy: Cauchy Spatial Correlation Structure: # # the correlation between two observations a distance d # # apart is 1/(1 + (d/r)^2). Letting r denote the range. # # corRExp: Exponential Spatial Correlation Structur # # corRExp2: Non-Separable Exponential Spatio-Temporal Correlation # # Structure # # corRExpwr: Powered Exponential Spatial Correlation Structure # # corRExpwr2: Non-Separable Powered Exponential Spatio-Temporal # # Correlation Structure # # corRExpwr2Dt: Non-Separable Temporally Integrated Powered Exponential # # Spatial Correlation Structure # # corRGaus: Gaussian Spatial Correlation Structure # # corRGneit: Gneiting Spatial Correlation Structure # # corRLin: Linear Spatial Correlation Structure # # corRMatern: Matern Spatial Correlation Structure # # corRSpher: Spherical Spatial Correlation Structure # # corRWave: Sine Wave Spatial Correlation Structure # ######################## ################################ # ######## corRCauchy: Cauchy Spatial Correlation Structure ############ # # corRCauchy(value = numeric(0), form = ~ 1, # # metric = c("euclidean", "maximum", "manhattan", "haversine"), # # radius = 3956) # # # # Argumentos: value: O padrão é numeric(0). # # form: fórmula unilateral form ~ S1 +...+ Sp, # # especificando # # covariáveis espaciais S1 através Sp. # # O padrão é 1 (corresponde a usar a ordem das # # observações nos dados como covariável) # # metric: especifica a distância métrica a ser utilizada # # euclidiana: é a padrão # # maximum: para a diferença máxima # # manhattan: para a soma das diferenças absolutas e # # Haversine: para a distância do grande-círculo # # (milhas) entre longitude/latitude. # # radius: utilizado na fórmula de Haversine. # # O padrão é o raio da Terra de 3956 milhas. # ######################## ################################ # Exemplo 1: # sp1 <- corRCauchy(form = ~ x + y + z) # sp1 require(ramps) spatDat <- data.frame(x = (0:4)/4, y = (0:4)/4) spatDat cs1Cauchy <- corRCauchy(1, form = ~ x + y) # cs1Cauchy cs1Cauchy <- Initialize(cs1Cauchy, spatDat) corMatrix(cs1Cauchy) cs2Cauchy <- corRCauchy(1, form = ~ x + y, metric = "man") cs2Cauchy <- Initialize(cs2Cauchy, spatDat) corMatrix(cs2Cauchy) # Exemplo 2: sp1 <- corRWave(form = ~ x + y + z) spatDat <- data.frame(x = (0:4)/4, y = (0:4)/4) spatDat cs1Wave <- corRWave(1, form = ~ x + y) cs1Wave <- Initialize(cs1Wave, spatDat) corMatrix(cs1Wave) cs2Wave <- corRWave(1, form = ~ x + y, metric = "man") cs2Wave <- Initialize(cs2Wave, spatDat) corMatrix(cs2Wave) ### DIC - Deviance Information Criterion ### # # DIC(object) # # Arguments : object: object returned by georamps # ### expand.chain - Expand MCMC Samples for georamps Model Fits ### # # expand.chain(object, n) # # Arguments: object: object returned by georamps. # n: additional number of times to iterate the MCMC sampler. ### georamps - Bayesian Geostatistical Model Fitting with RAMP ### # # General function for ???tting Bayesian geostatistical models using the # reparameterized and marginalized posterior sampling (RAMPS) algorithm # of Yan et al. (2007). # # georamps(fixed, random, correlation, data, subset, weights, # variance = list(fixed = ~ 1, random = ~ 1, spatial = ~ 1), # aggregate = list(grid = NULL, blockid = ""), kmat = NULL, # control = ramps.control(...), contrasts = NULL, ...) ## Load the included uranium datasets for use in this example data(NURE) ## Geostatistical analysis of areal measurements (Análise geoestatística para medições de área) NURE.ctrl1 <- ramps.control( iter = 25, beta = param(0, "flat"), sigma2.e = param(1, "invgamma", shape = 2.0, scale = 0.1, tuning = 0.75), phi = param(10, "uniform", min = 0, max = 35, tuning = 0.50), sigma2.z = param(1, "invgamma", shape = 2.0, scale = 0.1) ) NURE.fit1 <- georamps(log(ppm) ~ 1, correlation = corRExp(form = ~ lon + lat, metric = "haversine"), weights = area, data = NURE, subset = (measurement == 1), aggregate = list(grid = NURE.grid, blockid = "id"), control = NURE.ctrl1 ) print(NURE.fit1) summary(NURE.fit1) ## Analysis of point-source measurements (Análise do ponto ) NURE.ctrl2 <- ramps.control( iter = 25, beta = param(0, "flat"), sigma2.e = param(1, "invgamma", shape = 2.0, scale = 0.1, tuning = 0.75), phi = param(10, "uniform", min = 0, max = 35, tuning = 0.5), sigma2.z = param(1, "invgamma", shape = 2.0, scale = 0.1) ) NURE.fit2 <- georamps(log(ppm) ~ 1, correlation = corRExp(form = ~ lon + lat, metric = "haversine"), data = NURE, subset = (measurement == 2), control = NURE.ctrl2 ) print(NURE.fit2) summary(NURE.fit2) ## Joint analysis of areal and point-source measurements with ## prediction only at grid sites (Análise conjunta) NURE.ctrl <- ramps.control( iter = 25, beta = param(rep(0, 2), "flat"), sigma2.e = param(rep(1, 2), "invgamma", shape = 2.0, scale = 0.1, tuning = 0.75), phi = param(10, "uniform", min = 0, max = 35, tuning = 0.5), sigma2.z = param(1, "invgamma", shape = 2.0, scale = 0.1),z.monitor = NURE.grid ) NURE.fit <- georamps(log(ppm) ~ factor(measurement) - 1, correlation = corRExp(form = ~ lon + lat, metric = "haversine"), variance = list(fixed = ~ measurement), weights = area * (measurement == 1) + (measurement == 2), data = NURE, aggregate = list(grid = NURE.grid, blockid = "id"), control = NURE.ctrl ) print(NURE.fit) summary(NURE.fit) ## Discard initial 5 MCMC samples as a burn-in sequence ## (Descarte inicial de 5 amostras MCMC como uma seqüência burn-in) fit <- window(NURE.fit, iter = 6:25) print(fit) summary(fit) ## Deviance Information Criterion DIC(fit) ## Prediction at unmeasured sites (Previsão em locais não medidos) ct <- map("state", "connecticut", plot = FALSE) lon <- seq(min(ct$x, na.rm = TRUE), max(ct$x, na.rm = TRUE), length = 20) lat <- seq(min(ct$y, na.rm = TRUE), max(ct$y, na.rm = TRUE), length = 15) grid <- expand.grid(lon, lat) newsites <- data.frame(lon = grid[,1], lat = grid[,2], measurement = 1) pred <- predict(fit, newsites) plot(pred, func = function(x) exp(mean(x)), database = "state", regions = "connecticut", resolution = c(200, 150), bw = 5, main = "Posterior Mean", legend.args = list(text = "ppm", side = 3, line = 1)) plot(pred, func = function(x) exp(sd(x)), database = "state", regions = "connecticut", resolution = c(200, 150), bw = 5, main = "Posterior Standard Deviation", legend.args = list(text = "ppm", side = 3, line = 1)) ### os dados ### data(NURE) ## Map areal and point-source measurements ppm1 <- NURE$ppm[NURE$measurement == 1] level <- (max(ppm1) - ppm1) / diff(range(ppm1)) map("county", "connecticut", fill = TRUE, col = gray(level)) title("Connecticut Uranium Measurements") points(NURE$lon, NURE$lat) ## Map grid sites map("county", "connecticut") title("Regular Grid of Coordinates") points(NURE.grid$lon, NURE.grid$lat) ### param Initialization of georamps Model Parameter ### param(rep(0, 2), "flat") ## Random generation of initial values for an inverse-gamma prior param(rep(NA, 2), "invgamma", shape = 2.0, scale = 0.1) ## Independent normal priors param(rep(0, 2), "normal", mean = c(0, 0), variance = c(100, 100)) ## Correlated normal priors npv <- rbind(c(100, 25), c(25, 100)) param(rep(0, 2), "normal", mean = c(0, 0), variance = npv) ## Uniform prior and MCMC tuning parameter specification param(10, "uniform", min = 0, max = 100, tuning = 0.5) data(simJSS) ## Map areal and point-source measurements y <- simIowa$y[simIowa$areal == 1] level <- (max(y) - y) / diff(range(y)) map("county", "iowa", fill = TRUE, col = gray(level)) title("Simulated Iowa Measurements") points(simIowa$lon, simIowa$lat) ## Map grid sites map("county", "iowa") title("Regular Grid of Coordinates") points(simGrid$lon, simGrid$lat)