GeoR

Más información en:

• GeoR R-project
• Introducción a la geoestadística con el paquete GeoR - Ruben F Casal

El paquete geoR proporciona herramientas para el análisis de datos geoestadísticos en R (otra alternativa es el paquete gstat, por ejemplo…).

library(geoR)

Datos de prueba:

En la primera parte de esta sección consideraremos un proceso espacial sin tendencia, para ello usamos los datos s100 del paquete, que son especialmente simulados para la documentación del paquete.

data(s100) # Cargar datos estacionarios
class(s100)

[1] "geodata"

summary(s100)

Number of data points: 100 

Coordinates summary
        Coord.X    Coord.Y
min 0.005638006 0.01091027
max 0.983920544 0.99124979

Distance summary
        min         max 
0.007640962 1.278175109 

Data summary
      Min.    1st Qu.     Median       Mean    3rd Qu. 
-1.1676955  0.2729882  1.1045936  0.9307179  1.6101707 
      Max. 
 2.8678969 

Other elements in the geodata object
[1] "cov.model" "nugget"    "cov.pars"  "kappa"     "lambda"   

plot(s100)

También se pueden graficar los datos, con el color y tamaño en los puntos respecto al valor de cada dato.

points(s100, col="gray")

Análisis estructural

La primera etapa en el desarrollo de un análisis geoestadístico es la determinación de la dependencia espacial entre los datos medidos de una variable. Esta fase es también conocida como análisis estructural.

Variograma experimental

Estimador clásico y robusto
Si hay valores atípicos (o la distribción de los datos es asimétrica) puede ser preferible utilizar el estimador robusto. Se puede calcular este estimador estableciendo estimator.type = “modulus”

# Estimador clasico
ve = variog(s100, max.dist=0.6) # Variograma de los datos.

variog: computing omnidirectional variogram

names(ve)

 [1] "u"                "v"                "n"               
 [4] "sd"               "bins.lim"         "ind.bin"         
 [7] "var.mark"         "beta.ols"         "output.type"     
[10] "max.dist"         "estimator.type"   "n.data"          
[13] "lambda"           "trend"            "pairs.min"       
[16] "nugget.tolerance" "direction"        "tolerance"       
[19] "uvec"             "call"            

# Estimador robusto
#ve = variog(s100, max.dist=0.6, bin.cloud=T,)
ve1 = variog(s100, max.dist=0.6, estimator.type = "modulus")

variog: computing omnidirectional variogram

par(mfrow=c(1,2)) 
plot(ve, main="Estimador clásico") 
plot(ve1, main="Estimador robusto")

par(mfrow=c(1,1))

Variograma direccional

En el caso de anisotropía, también se pueden obtener variogramas direccionales con la función variog mediante los argumentos direction y tolerance. Por ejemplo, para calcular un variograma en la dirección de 60 grados (con la tolerancia angular por defecto de 22.5 grados):

vario.60 <- variog(s100, max.dist = 0.6, direction = pi/3) 

variog: computing variogram for direction = 60 degrees (1.047 radians)
        tolerance angle = 22.5 degrees (0.393 radians)

plot(vario.60)

Variograma en 4 direcciones

Para estudiar si hay anisotropía, se pueden cálcular de forma rápida variogramas direccionales con la función variog4. Por defecto calcula cuatro variogramas direccionales, correspondientes a los ángulos 0, 45, 90 y 135 grados:

var4 <- variog4(s100, max.dist=0.6)

variog: computing variogram for direction = 0 degrees (0 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 45 degrees (0.785 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 90 degrees (1.571 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 135 degrees (2.356 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing omnidirectional variogram

plot(var4)

Ajuste visual

“A ojo”: representando diferentes modelos sobre un variograma empírico.

# Variograma experimental
ve = variog(s100, max.dist=0.6, estimator.type = "modulus")

# Ajuste visual
va = eyefit(ve1)
va

Modelo exponencial - Ajuste teórico

Cuando se utilizan las funciones variofit y likfit para la estimación de parámetros, el efecto pepita (nugget) puede ser estimado o establecido a un valor fijo. Lo mismo ocurre con los parámetros de suavidad, anisotropía y transformación de los datos. También se dispone de opciones para incluir una tendencia. Las tendencias pueden ser polinomios en función de las coordenadas y/o funciones lineales de otras covariables.

Mínimos cuadrados ordinarios

spar = 0.82-0.16 # silla parcial
v_teorico_OLS = variofit(ve1, ini=c(spar,0.4), 
                     cov.model="exp", fix.nugget = F,
                     weights = "equal",nug = 32, max.dist = 5)

variofit: covariance model used is exponential 
variofit: weights used: equal 
variofit: minimisation function used: optim 

unreasonable initial value for sigmasq + nugget (too high)unreasonable initial value for nugget (too high)

Mínimos cuadrados ponderados

v_teorico_WLS = variofit(ve1, ini=c(spar,0.4), 
                     cov.model="exp", fix.nugget=F,
                     weights="cressie",nug=32, max.dist = 5)

variofit: covariance model used is exponential 
variofit: weights used: cressie 
variofit: minimisation function used: optim 

unreasonable initial value for sigmasq + nugget (too high)unreasonable initial value for nugget (too high)

Maxima verosimilitud

v_ML = likfit(s100, ini=c(spar,0.4)) 

---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.

Máxima verosimilitud restringida

v_REML <- likfit(s100, ini = c(1, 0.5), lik.method = "RML")

---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.

Resultados

v_teorico_OLS # OLS

variofit: model parameters estimated by OLS (ordinary least squares):
covariance model is: exponential
parameter estimates:
    tausq   sigmasq       phi 
   0.1747 2400.5764 1971.1998 
Practical Range with cor=0.05 for asymptotic range: 5905.187

variofit: minimised sum of squares = 0.0487

v_teorico_WLS # WLS

variofit: model parameters estimated by WLS (weighted least squares):
covariance model is: exponential
parameter estimates:
    tausq   sigmasq       phi 
   0.1895 6927.1545 5821.4660 
Practical Range with cor=0.05 for asymptotic range: 17439.55

variofit: minimised weighted sum of squares = 33.7187

v_ML # ML

likfit: estimated model parameters:
    beta    tausq  sigmasq      phi 
"0.7766" "0.0000" "0.7517" "0.1827" 
Practical Range with cor=0.05 for asymptotic range: 0.5473847

likfit: maximised log-likelihood = -83.57

v_REML # REML

likfit: estimated model parameters:
    beta    tausq  sigmasq      phi 
"0.7478" "0.0000" "0.8473" "0.2102" 
Practical Range with cor=0.05 for asymptotic range: 0.6296295

likfit: maximised log-likelihood = -81.53

Estimacion de los parametros con diferentes métodos,

plot(ve1, main = "Estimador empírico y modelo exponencial ")
lines(v_ML, max.dist = 0.6) # ML

lines(v_teorico_OLS, lty = 2, max.dist = 0.6) # OLS
lines(v_teorico_WLS, lty = 2, lwd = 2, max.dist = 0.6) # WLS

lines(v_REML, lty = 2, lwd = 2, max.dist = 0.6) # REML
legend(0.3, 0.3, legend = c("ML", "OLS", "WLS", "REML"),
        lty = c(1, 2, 2, 1), lwd = c(1,1,2,2))

Validación cruzada

Para verificar si un modelo de variograma describe adecuadamente la dependencia espacial de los datos (p.e. comparar modelos), se emplea normalmente la técnica de validación cruzada, función xvalid en geoR. Por defecto la validación se realiza sobre los datos eliminando cada observación (y utilizando las restantes para predecir), aunque se puede utilizar un conjunto diferente de posiciones (o de datos) mediante el argumento location.xvalid (y data.xvalid).

# Validacion cruzada

There were 27 warnings (use warnings() to see them)

s100.cv.ols <- xvalid(s100, model=v_teorico_OLS)

xvalid: number of data locations       = 100
xvalid: number of validation locations = 100
xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 
xvalid: end of cross-validation

s100.cv.ml <- xvalid(s100, model=v_ML)

xvalid: number of data locations       = 100
xvalid: number of validation locations = 100
xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 
xvalid: end of cross-validation

# Errores
RMSE_ML = round(sqrt(mean(s100.cv.ml$error^2)),3)
RMSE_OLS = round(sqrt(mean(s100.cv.ols$error^2)),3)
RMSE_ML # Este

[1] 0.479

RMSE_OLS 

[1] 0.529

Se pueden obtener envolventes (envelopes, i.e. valores máximos y mínimos aproximados) del variograma empírico mediante simulación: Bajo un modelo de variograma, para ilustrar la variabilidad del variograma empírico.

test_mc = variog.mc.env(s100, obj.variog=ve1)

variog.env: generating 99 simulations by permutating data values
variog.env: computing the empirical variogram for the 99 simulations
variog.env: computing the envelops

plot(ve1, envelope=test_mc)

Kriging

El paquete geoR dispone de opciones para los métodos kriging tradicionales.

Para obtener una rejilla discreta de predicción puede ser de utilidad la función expand.grid:

# Defining a prediction grid

There were 40 warnings (use warnings() to see them)

loci <- expand.grid(seq(0,1,l=21), seq(0,1,l=21))
plot(loci)

El comando para realizar kriging ordinario con variograma ml sería:

# predicting by ordinary kriging

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls
9: In readChar(file, size, TRUE) : truncating string with embedded nuls

kc <- krige.conv(s100, loc=loci,
                krige=krige.control(obj.model = v_ML))

krige.conv: model with constant mean
krige.conv: Kriging performed using global neighbourhood 

image(kc)

contour(kc)