sábado, 12 de mayo de 2012

The R Journal >> Current Issue

The R Journal >> Current Issue

 RSS Feed
ISSN: 2073-4859
The  Journal

Volume 3/2, December 2011

Download complete issueRefereed articles may be downloaded individually using the links below. [Bibliography of refereed articles]

Table of Contents

Editorial3

Contributed Research Articles

Creating and Deploying an Application with (R)Excel and R 
Thomas Baier, Erich Neuwirth and Michele De Meo
5
glm2: Fitting Generalized Linear Models with Convergence Problems 
Ian C. Marschner
12
Implementing the Compendium Concept with Sweave and DOCSTRIP 
Michael Lundholm
16
Watch Your Spelling! 
Kurt Hornik and Duncan Murdoch
22
Ckmeans.1d.dp: Optimal k-means Clustering in One Dimension by Dynamic Programming 
Haizhou Wang and Mingzhou Song
29
Nonparametric Goodness-of-Fit Tests for Discrete Null Distributions 
Taylor B. Arnold and John W. Emerson
34
Using the Google Visualisation API with R 
Markus Gesmann and Diego de Castillo
40
GrapheR: a Multiplatform GUI for Drawing Customizable Graphs in R 
Maxime Hervé
45
rainbow: An R Package for Visualizing Functional Time Series 
Han Lin Shang
54

Programmer's Niche

Portable C++ for R Packages 
Martyn Plummer
60

News and Notes

R's Participation in the Google Summer of Code 201164
Conference Report: useR! 201168
Forthcoming Events: useR! 201270
Changes in R72
Changes on CRAN84
News from the Bioconductor Project86
R Foundation News87
Read more...

R-project (R news & tutorials)

Read more...

Shareaholic for Safari - The best way to share, tweet, bookmark, save and e-mail webpages

Read more...

lunes, 21 de noviembre de 2011

RClimate Images, Links to Data & R Script Files

Esta es una página interesante sobre análisis de series temporales de clima con R. Encontrarán datos y los scripts con los que generaron cada gráfico.
Espero les sea útil, saludos!

RClimate Images, Links to Data & R Script Files


Climate Trends
RClimate Images, R Script and Data File Links

Return to: Climate Charts & Graphs Blog

Topic Click Image to Enlarge Links
GISS Land & Ocean Temperature Anomaly Trends

Shows monthly GISS Land & Sea Temperature Anomaly - C trend from 1880 until latest month available from NASA GISSTemp. Base period is 1951-1980.

Decadal means shown in blue steps and last monthly value highlighted in red. Decadal means show steep rise since mid 1970s.

Data

R script

GISS Land & Ocean Temperature Anomaly Map

Users can use NASA's interactive tool to make a global map of temperature anomalies on a 2x2 degree grid at this link.

Users can also download a text file of the data.

I have developed an RClimate script to map the NASA GISTemp anomaly data using readily available R mapping packages.

Data
R script

RSS Global Land & Ocean Temperature Anomaly Trends

Shows monthly satellite based RSS Land & Sea Temperature Anomaly - C from Nov., 1978 to last available month.

Annual means are shown in blue steps and last monthly value highlighted in red.

Data

R script

UAH Channel 5 Long Term and Current Month Trends

Uses R's figure inside figure capability to show both the long term anomaly trend as well as the current month daily trend.

R script includes data retrieval, conversion of wide format table to long format, calculation of monthly averages and production of 2 figures in one chart.

Data

R script

Consolidated Global Land and Ocean Temperature Anomaly Series Since 1880

Shows monthly temperature anomaly data for 5 major temperature L&O anomaly series: GISS, NOAA, Hadley, RSS and UAH. First 3 use station data, the last 2 use satellite based observations.

The CSV data file is updated monthly from source agency data files. All source agency data is downloaded so that any changes in source data will be reflected in the consolidated file.

RClimate.txt includes a series of custom functions() that I have written to facilitate LOTA download and analysis. Interested RUsers can source my RClimate functions with this link: http://processtrends.com/Files/RClimate.txt

Data

Consolidated Global Land and Ocean Temperature Anomaly Series for Selected Month Since 1880

This panel chart and table shows the monthly anomaly trends for the 5 major global temperature anomaly series and a table that shows how the selected month ranks over the entire record for each series.

Data
R Script

Year-to-data Consolidated Global Land and Ocean Temperature Anomaly Series Since 1880

This 5 panel chart and table shows the year-to-date anomaly trends for the 5 major global temperature anomaly series and a table that shows how the current year ranks over the entire record for each series.

Data
R Script
Common Baseline for 5 LOTA Series

The 5 global land-ocean temperature anomaly (LOTA) series use different baseline periods, making direct comparisons between the series more difficult than it would be if each series had the same baseline period.

This post shows how to convert the 5 major LOTA series to a common baseline. Links to on-line source data file and RClimate script are provided. Here is long term LOTA trends using a 133 month moving average and 1979-2008 baseline.

Data

R script

SST Anomaly (Monthly)

Shows monthly NOAA - NCDC monthly Sea Surface Temperature - C (SST) Anomaly trends, starting in 1880 to most current month.

Decadal means shown in blue steps and last monthly value highlighted in red.

Data

R script

Mean Sea Level Anomaly

Shows monthly global mean sea level trends of NOAA's Laboratory of Satellite Altimetry global mean sea level data. The chart shows the missions (TOPEX, Jason-1, Jason-2) by color code and highlights the latest monthly data point.

The seasonal signal has been removed and an inverted barometer has been applied. No glacial isostatic adjustment (GIA) has been made. The overall rate of increase is 2.9 mm/yr. GIA estimates are in the 0.2 - 0.5 mm/yr range.

This RClimate script requires user to download and save NetCDF file to own PC. Make sure to adjust source data file link in R script to properly locate the data file on your PC.

Data

R Script

Mean Sea Level Anomaly

Shows monthly global mean sea level trends of University of Colorado - Boulder's satellite based data. Source data has not been updated since Sept., 2009. I will resume updating when source data is updated.

Chart shows the monthly change in mean sea level in red, 60-days smoothing value in blue and overall trend in green. The overall rate of increase is 3.1 mm/yr, with a +- error band of 0.4 mm/yr.

Data
CO2 Trends - Mauna Loa

Shows monthly CO2 (ppmv). Data from 1958 to last available month.

Last monthly value highlighted in red.

Data

R script

Total Solar Irradiance (TSI)

Shows NASA Solar Radiation & Climate Experiment (SORCE) Mission Data from 2/25/03 to present. Data updated weekly.

Data

R script

Nino_34 SST Anomaly

Shows NOAA's weekly trend in Nino 3.4 SST anomalies from January, 1990 to most recent week. Data reflects week centered on Wednesdays.

Periods with negative anomalies (La Nina like conditions) are shown in blue and periods with positive anomalies (El Nino like conditions) are shown in red.

Most recent reading is highlighted in black.

Data

R script

Pacific Decadal Oscillation

Downloads and plots University of Washington, Joint Institute for the Study of Atmosphere and Ocean (JISAO) monthly Pacific Decadal Oscillation (PDO) from January, 1900 to latest available month at time script is run.

Negative PDO months are shown in light blue and positive months are shown in pink. The most recent reading is highlighted in red.

The title and legend include the latest data month.

Data

R script

Polar Amplification: 2000 to 2009

Shows NASA GISS temperature anomaly by latitude zone for 2000-2009 compared to baseline period (1951-1980).

Chart demonstrates significant Northern Hemisphere temperature anomaly increase with increasing latitudes. This means that Arctic regions are warming much more rapidly than the global mean. While the global mean temperature anomaly increased 0.6 C in 2000-2009, the Arctic region anomalies increased nearly 1.8 C.

Data

R script

Arctic Sea ice Extent: Comparison of 2010 and 2007 Trends

Retrieves JAXA daily Arctic Sea Ice Extent (SIE) csv file, plots 2007 and 2010 trends, highlights the maximum and current values for both series, calculates the decline rate for selected time periods and displays as table in plot.

R script

Data

Arctic Sea Ice Extent Daily Change Comparison: 2007 vs 2010
Arctic Sea Ice Extent by Month

Shows National Snow and Ice Data Center's (NSIDC) monthly Arctic Sea Ice extent trends in millions of square kilometers from 1979 to latest month.

Arctic sea ice extent includes areas with at least 15% sea ice concentrations.

Since sea ice extent varies seasonally, the chart shows the trend for each month. The latest month trend shown is shown in red. All months show declining sea ice extent, which is consistent with polar amplification, rising SSTA and rising global mean temperature anomaly trends.

2nd chart shows just latest month. Trend line added to display OLS trend line for monthly data.

Data

R script

NSIDC Annual Arctic Sea ice Extent Max, Min & Melt

Data

R script

Arctic Sea Ice Extent Forecast - Sept 2011

I used a quadratic regression model to fit 1980 - 2010 NSDIC September sea ice extent (SIE) data. It hen forecast the Sept., 2011 SIE with this model.

See Climate Charts & Graphs post.

Data
R script

Daily JAXA Arctic Sea Ice Extent for latest Month

Shows JAXA daily Arctic Sea Ice extent trends in millions of square kilometers from 2002 to latest month.

Arctic sea ice extent includes areas with at least 15% sea ice concentrations.

Since sea ice extent varies seasonally, the chart shows the trend for the current month for each year from 2002 to the present. The latest daily reading is shown in red for the current year and previous years.

R script

Data

Arctic Sea Ice Area Anomaly

The University of Illinois at Urbana Champaign (UIUC)'s The Cryosphere Today (link) tracks Arctic sea ice extent and area on a daily basis, They maintain a chart called of daily sea ice area anomalies extending back to 1979. They call it the Tail of the Tape. Since it is so long, it can be difficult to read (link).

I have developed a similar chart that displays in a single window so that readers can see the entire record in one glance. While the Cryoshpere Today gives more detail, I think my chart gives a better overview of the long term SIA anomaly trends.

Data
Temperature Sounding Profile

Here's an RClimate based plot of Washington DC area atmospheric temperature soundings that will be updated daily.

R Script
Arctic Oscillation

NOAA updates monthly Arctic Oscillation Index (link).

Recent AO data for past 120 days.

Data
R script

Data
Climate Oscillations and GISS Anomalies

This series constructs monthly file of GISS anomaly, PDO,AMO and ENSO data series and performs detailed analysis of the consolidated data file.

R script to develop consolidated file of monthly GISS, PDO, AMO, ENSO data series

R script

Data

Data

R script

Read more...

viernes, 11 de noviembre de 2011

Una charla sobre la estadística espacial... UAL, 10/11/2012.

Esta es la presentación de una charla sobre modelos espacio-temporales aplicados en ecología, que he dado el 10 de noviembre en la UAL en el Máster en Evaluación del Cambio Global. El objetivo fue hablar de la importancia de considerar las escalas espacial y temporal en los análisis de datos de ecología, haciendo énfasis en la consideración de estructuras de dependencia o autocorrelación. Dada la amplia variedad de técnicas estadísticas para este tipo de análisis, la presentación pretendió dar una breve guía donde resumir las principales técnicas disponibles según el tipo de datos que estemos analizando. Se trabajó con ejemplos ecológicos aplicados en el software R, utilizando varias librerías.
A pesar de haber sido una charla demasiado ambiciosa, espero que haya servido al menos para promover nuevas ideas sobre cómo abordar los datos espaciales y cómo hacerse nuevas preguntas desde un planteo formal y riguroso.
Aquí les dejo el material para cualquier curioso, si desean más ejemplos de aplicación en R, basta con mirar en el blog o enviarme un mail con preguntas (que espero poder responder).
Saludos!


presentacion
Read more...

martes, 22 de febrero de 2011

Códigos y ejemplos en R de aplicación geoestadísstica. Spatial prediction of soil moisture - spatial-analyst.net

Spatial prediction of soil moisture - spatial-analyst.net: "- Enviado mediante la barra Google"

Códigos y ejemplos de aplicación geoestadísstica.
Read more...

lunes, 10 de enero de 2011

esTcena | Santander Meteorology Group

esTcena | Santander Meteorology Group: "- Enviado mediante la barra Google"
Read more...

viernes, 10 de diciembre de 2010

Análisis geoestadístico con R. Paquete geoR. parte II.

Más funciones para el análisis de datos geoestadístico según el libro"Model-based Geostatistics" de Peter J. Diggle & Paulo Justiniano Ribeiro Jr, 2007.
Iré agregando material teórico y práctico sobre el tema.
Book data set: http://www.leg.ufpr.br/doku.php/pessoais:paulojus:mbgbook:datasets

##############################################################
#### Diggle Model Based Geoestadistics. Paquete geoR en R ####
##############################################################
library(geoR)
data(elevation)
###2.8 perspectiva general de modelos geoestadísiticos:
#análisis exploratorio no-espacial
plot(elevation)


with(elevation, hist(data, main = "", xlab = "elevation"))
with(elevation, plot(coords[, 1], data, xlab = "W-E", ylab = "elevation data", pch = 20, cex = 0.7))
lines(lowess(elevation$data ~ elevation$coords[, 1]))
with(elevation, plot(coords[, 2], data, xlab = "S-N", ylab = "elevation data", pch = 20, cex = 0.7))
lines(with(elevation, lowess(data ~ coords[, 2])))
#modelo lineal con covariables. The values "1st" and "2nd" passed to the argument trend are aliases to indicate first-degree and second-degree polynomials on the coordinates.
points(elevation, cex.max = 2.5)
points(elevation, trend = "1st", pt.div = 2, abs = T, cex.max = 2.5)
points(elevation, trend = "2nd", pt.div = 2, abs = T, cex.max = 2.5)
#calculo de los variogramas empíricos para los datos y los residuales
plot(variog(elevation, uvec = seq(0, 5, by = 0.5)), type = "b")
res1.v <- variog(elevation, trend = "1st", uvec = seq(0, 5, by = 0.5))
plot(res1.v, type = "b")
res2.v <- variog(elevation, trend = "2nd", uvec = seq(0, 5, by = 0.5))
lines(res2.v, type = "b", lty = 2)
#variograma residual y envolventes de simulación bajo permutación aleatoria de los residuales.
set.seed(231)
mc1 <- variog.mc.env(elevation, obj = res1.v)
plot(res1.v, env = mc1, xlab = "u")
mc2 <- variog.mc.env(elevation, obj = res2.v)
plot(res2.v, env = mc2, xlab = "u")


#estimaciones por ML (maximum likelihood) del modelo Gaussiano sin y con el término de tendencia.
ml0 <- likfit(elevation, ini = c(3000, 2), cov.model = "matern", kappa = 1.5); ml0
ml1 <- likfit(elevation, trend = "1st", ini = c(1300, 2), cov.model = "matern", kappa = 1.5); ml1
#interpolación espacial con kriging
locs <- pred_grid(c(0, 6.3), c(0, 6.3), by = 0.1)
KC <- krige.control(type = "sk", obj.mod = ml0)
sk <- krige.conv(elevation, krige = KC, loc = locs)
KCt <- krige.control(type = "sk", obj.mod = ml1, trend.d = "1st", trend.l = "1st")
skt <- krige.conv(elevation, krige = KCt, loc = locs)
#selección de funciones gráficas para producir mapas usando argumentos opcionales para asegurar que los pares usan la misma escala gris
pred.lim <- range(c(sk$pred, skt$pred))
sd.lim <- range(sqrt(c(sk$kr, skt$kr)))
image(sk, col = gray(seq(1, 0, l = 51)), zlim = pred.lim)
contour(sk, add = T, nlev = 6)
points(elevation, add = TRUE, cex.max = 2)
image(skt, col = gray(seq(1, 0, l = 51)), zlim = pred.lim)
contour(skt, add = T, nlev = 6)
points(elevation, add = TRUE, cex.max = 2)
image(sk, value = sqrt(sk$krige.var), col = gray(seq(1, 0, l = 51)), zlim = sd.lim)
contour(sk, value = sqrt(sk$krige.var), levels = seq(10, 27, by = 2), add = T)
points(elevation$coords, pch = "+")
image(skt, value = sqrt(skt$krige.var), col = gray(seq(1, 0, l = 51)), zlim = sd.lim)
contour(skt, value = sqrt(skt$krige.var), levels = seq(10, 27, by = 2), add = T)
points(elevation$coords, pch = "+")
#dejar fuera el efecto de la covariable: remoción de tendencia
plot(elevation, low = TRUE, trend = ~coords[, 2], qt.col = 1)

### 3.13 Modelos Gaussianos para datos geoestadísticos
#cálculo  y gráfico de la función de correlación estándar. Se elige la familia de correlaciones con el argumento cov.model.
x <- seq(0, 1, l = 101)
plot(x, cov.spatial(x, cov.model = "mat", kappa = 0.5, cov.pars = c(1, 0.25)), type = "l", xlabel = "u", ylabel = expression(rho(u)), ylim = c(0, 1))
lines(x, cov.spatial(x, cov.model = "mat", kappa = 1.5, cov.pars = c(1, 0.16)), lty = 2)
lines(x, cov.spatial(x, cov.model = "mat", kappa = 2.5, cov.pars = c(1, 0.13)), lty = 3)

#generar simulaciones de procesos Gaussianos 2D. La función grf() especifica el modelo y las localizaciones para lso valores simulados requeridos.
set.seed(159) #asegura que las simulaciones se gneran con el mismo número aleatorio de semilla, y por tanto las realizaciones simuladas tienen diferencias solo por los valores diferentes de los parámetros del modelo.
image(grf(100^2, grid = "reg", cov.pars = c(1, 0.25)), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
set.seed(159)
image(grf(100^2, grid = "reg", cov.pars = c(1, 0.13), cov.model = "mat", kappa = 2.5), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
#simulaciones de modelos anisotrópicos
set.seed(421)
image(grf(201^2, grid = "reg", cov.pars = c(1, 0.25), aniso.pars = c(pi/3, 4)), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
set.seed(421)
image(grf(201^2, grid = "reg", cov.pars = c(1, 0.25), aniso.pars = c(3 * pi/4, 2)), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
#In this case, we have used a stationary Gaussian model with mean mu = 850, nugget variance tau^2 = 100, signal variance sigma^2 = 3500 and MatLern correlation function with k = 2.5 and phi= 0.8.
sim1 <- grf(100, cov.pars = c(1, 0.15), cov.model = "matern", kappa = 1.5)
points(sim1)
data(elevation)
sim2 <- grf(grid = elevation$coords, cov.pars = c(3500, 0.8), nugget = 100)
sim2$data <- sim2$data + 850
points(sim2)
### 4.6 Modelos lineales generalizados (GLM) para datos geoestadísticos
 #simulación con el GLM de un modelo de Poisson
 #simulate a realisation of the Gaussian process at these locations with mu= 0.5, sigma^2 = 2 and MatLern correlation function with k = 1.5, phi = 0.2.
  set.seed(371)
  cp <- expand.grid(seq(0, 1, l = 10), seq(0, 1, l = 10))
  s <- grf(grid = cp, cov.pars = c(2, 0.2), cov.model = "mat", kappa = 1.5)
  image(s, col = gray(seq(1, 0.2, l = 21)))
  lambda <- exp(0.5 + s$data)
  y <- rpois(length(s$data), lambda = lambda)
  text(cp[, 1], cp[, 2], y, cex = 1.5)
 #extensión del modelo, incluyendo variación no-espacial extra de Poisson
  lambda <- exp(s$data + tau * rnorm(length(s$data)))
 #simulación con GLM del modelo Bernoulli
  set.seed(34)
  locs <- seq(0, 1, l = 401)
  s <- grf(grid = cbind(locs, 1), cov.pars = c(5, 0.1), cov.model = "matern", kappa = 1.5)
  p <- exp(s$data)/(1 + exp(s$data))
  ind <- seq(1, 401, by = 8)
  y <- rbinom(p[ind], size = 1, prob = p)
  plot(locs[ind], y, xlab = "locations", ylab = "data")
  lines(locs, p)

 #simulación con GLM del modelo Binomial. Con [Y (x)|S] ¡« Bin{n, p(x)} with n = 5 and p(x) = exp{mu + S(x)}/[1 + exp{mu + S(x)}],
  set.seed(23)
  s <- grf(60, cov.pars = c(5, 0.25))
  p <- exp(2 + s$data)/(1 + exp(2 + s$data))
  y <- rbinom(length(p), size = 5, prob = p)
  points(s)
  text(s$coords, label = y, pos = 3, offset = 0.3)


### Estimación de parámetros clásica: estimaciónd e tendencia, variogramas, estimación de la estructura de covarianza (OLS, WLS, etc), ML, estimación paramétrica por modelos lineales generalizados, MCML, etc.
 #calculo de variograma con diferentes opciones del tamaño bin
   data(elevation)
   plot(variog(elevation, option = "cloud"), xlab = "u", ylab = "V(u)")
   plot(variog(elevation, uvec = seq(0, 8, by = 0.5)), xlab = "u", ylab = "V(u)")
   plot(variog(elevation, trend = "2nd", max.dist = 6.5), xlab = "u", ylab = "V(u)")


    data(s100) #datos simulados con un modelo con media=cero, varianza=uno, nugget=0, y función de correlación exponencial con phi=.3.
    v1 <- variog(s100); plot(v1)
    v2 <- variog(s100, uvec = seq(0, 1, by = 0.1)); plot(v2)
    v3 <- variog(s100, max.dist = 1); plot(v3)
    round(v1$u, dig = 2) #The midpoints of the bins obtained from the commands listed above
    round(v2$u, dig = 2)
    round(v3$u, dig = 2)

    data(ca20)
    plot(ca20,trend= area+altitude, low=T) #sugiere tendencia cuadratica
    plot(variog(ca20, max.dist = 510)) #construimos distintos variogramas: observamos que incorporando sub-regoses al modelo este se comporta aproximadamente estacionario en los residuales
    plot(variog(ca20, trend = ~area, max.dist = 510))
    plot(variog(ca20, trend = ~area + altitude, max.dist = 510))
    t.all <- trend.spatial(ca20, trend = ~area + altitude, add = "2nd")
    plot(variog(ca20, trend = ~t.all, max.dist = 510))
set.seed(83) #variogramas de las simulaciones del proceso Gaussiano, cada gráfico indica el variograma verdadero y la curva suavizada, así como los variogramas empíricos de las tres simulaciones del proceso en 100 localizaciones. sim1 <- grf(100, cov.pars = c(1, 0.05), cov.model = "mat", kap = 1.5, nsim = 3) #Parameter values are k 1.5, sigma^2= 1, tau^2 = 0 in both panels. In the left-hand panel, phi= 0.05, in the right-hand panel phi= 0.2. plot(variog(sim1, max.dist = 1), type = "l", lty = 1:3, col = 1) lines.variomodel(seq(0, 1, l = 100), cov.model = "mat", kap = 1.5, cov.pars = c(1, 0.05), nug = 0) set.seed(83) sim2 <- grf(100, cov.pars = c(1, 0.2), cov.model = "mat", kap = 1.5, nsim = 3) plot(variog(sim2, max.dist = 1), type = "l", lty = 1:3, col = 1) lines.variomodel(seq(0, 1, l = 100), cov.model = "mat", kap = 1.5, cov.pars = c(1, 0.2), nug = 0)
#estimación de parámetros en modelos gaussianos con particular foco en los parámetros que definen la estructura de covarianza del modelo. s100.v <- variog(s100, max.dist = 1) #estimación mediante variograma plot(s100.v); lines.variomodel(seq(0, 1, l = 100), cov.pars = c(0.9, 0.2), cov.model = "mat", kap = 1.5, nug = 0.2)

wls <- variofit(s100.v, ini = c(0.9, 0.2), cov.model = "mat", kap = 1.5, nug = 0.2); wls #estimación mediante WLS ml <- likfit(s100, cov.model = "mat", kap = 1.5, ini = c(0.9, 0.2), nug = 0.2); ml #estimación por ML reml <- likfit(s100, cov.model = "mat", kap = 1.5, ini = c(0.9, 0.2), nug = 0.2, met = "reml"); reml #estimación por ML restringido #ajuste de 5 modelos alternaticos a los datos de suelo y sus covariables. data(ca20) m1 <- likfit(ca20, ini = c(100, 200), nug = 50) m2 <- likfit(ca20, trend = ~area, ini = c(60, 100), nug = 40) m3 <- likfit(ca20, trend = ~area + altitude, ini = c(60, 100), nug = 40) m4 <- likfit(ca20, trend = ~area + coords, ini = c(60, 100), nug = 40) m5 <- likfit(ca20, trend = ~area + altitude + coords, ini = c(60, 100), nug = 40) ### Predicción espacial #predicción por "conventional kriging" data(SIC) ml <- likfit(sic.all, ini = c(100, 40), nug = 10, lambda = 0.5, kappa = 1) #estimjación de los parámetros del modelo gr <- pred_grid(sic.borders, by = 7.5) #definimos una grilla para lso puntos de predicción KC <- krige.control(obj.model = ml) #pasa los parámetros del modelo mediante "ordinary krigning" OC <- output.control(n.pred = 1000, simul = TRUE, thres = 250) #incluye opciones para generar y guardar 1000 simulaciones de las distribuciones condicionales de S dado Y y definir un valor umbral de 250 que será usado para calcular la probabilidad de excedencia en cada localización de la predicción. set.seed(2419) #permite la reproducción de los resultados de la simulación pred <- krige.conv(sic.all, loc = gr, borders = sic.borders,krige = KC, out = OC) image(pred, col = gray(seq(1, 0.1, l = 21)), zlim = predlim, x.leg = c(0, 350), y.leg = c(-60, -30)) image(pred, loc = gr, val = sqrt(pred$krige.var), zlim = selim, col = gray(seq(1, 0.1, l = 21)), x.leg = c(0, 350), y.leg = c(-60, -30)) dim(pred$simulations) A200 <- apply(pred$simul, 2, function(y) sum(y > 200)/length(y)) hist(A200, main = "") image(pred, val = 1 - pred$prob, col = gray(seq(0.9, 0.1, l = 41)), x.leg = c(0, 350), y.leg = c(-60, -30))



#el mismo análisis para los datos de calcio en el suelo
   data(ca20)
   fit <- likfit(ca20, ini = c(100, 60), trend = ~area)
   gr <- pred_grid(ca20$borders, by = 10)
   gr0 <- polygrid(gr, borders = ca20$border, bound = T)
   ind.reg <- numeric(nrow(gr0))
   ind.reg[.geoR_inout(gr0, ca20$reg1)] <- 1
   ind.reg[.geoR_inout(gr0, ca20$reg2)] <- 2
   ind.reg[.geoR_inout(gr0, ca20$reg3)] <- 3
   ind.reg <- as.factor(ind.reg)
   KC <- krige.control(trend.d = ~area, trend.l = ~ind.reg, obj.model = fit)
   ca20pred <- krige.conv(ca20, loc = gr, krige = KC)
   par(mar = c(2.8, 2.5, 0.5, 0.5), mgp = c(1.8, 0.7, 0), mfrow = c(1, 2))
   image(ca20pred, loc = gr, col = gray(seq(1, 0, l = 21)), x.leg = c(4930, 5350), y.leg = c(4790, 4840))
   polygon(ca20$reg1)
   polygon(ca20$reg2)
   polygon(ca20$reg3)
   image(ca20pred, loc = gr, val = sqrt(ca20pred$krige.var), col = gray(seq(1, 0, l = 21)), x.leg = c(4930, 5350), y.leg = c(4790, 4840))
   polygon(ca20$reg1)
   polygon(ca20$reg2)
   polygon(ca20$reg3)
   coords <- cbind(c(0.2, 0.25, 0.6, 0.7), c(0.1, 0.8, 0.9, 0.3))
   KC <- krige.control(ty = "ok", cov.model = "mat", kap = 1.5, nug = 0.1, cov.pars = c(1, 0.1))
   krweights(coords, c(0.5, 0.5), KC)


Read more...

Libros para descargar (gratis) sobre Modelos Espacio-Temporales