Modelamiento

La clase de modelos estables de correlación están datos por:

\[ \rho(||\boldsymbol{h}||,\alpha,\beta) = exp\left\{ -\left(\frac{||\boldsymbol{h}||}{\alpha}\right)^{\beta} \right\} \]

donde \(||\boldsymbol{h}||\geq0\),\(\alpha>0\),\(0<\beta\leq2\).

En GeoModels (Bevilacqua, Morales-Oñate, and Caamaño-Carrillo (2018)), este modelo de correlación corresponde a:

Exponential, exponential, exp o exponential son las especificaciones en el parámetro corrmodel. Esto es \(\beta = 1\)
Gauss, gauss, Gaussian and gaussian son las especificaciones en el parámetro corrmodel. Esto es \(\beta = 2\)
Stable o stable son las especificaciones en el parámetro corrmodel. En este caso ambos parámetros son libres.

CorFunStable <- function(lag,R_power,scale)
{
  rho=exp(-R_pow(lag/scale,R_power))
  return (rho)
}

Ejemplo:

scale = 1.2/3
R_power = 1
curve(CorFunStable(x,R_power,scale), ylab=expression(paste(rho,"(",h,")")),0,2,xlab = "h")
abline(v=scale*3)

GeoModels

Usando la librería GeoModels

library(GeoModels)

###############################################################
############ Ejemplos de RFs Gaussianos ################
###############################################################

################################################################
###
### Estimación de un RF Gaussiano con 
### correlación exponencial usando verosimilitud por parejas
###############################################################

# Se define las coordenadas espaciales de los puntos:
set.seed(3)
N <- 400  # número de locaciones espaciales
x <- runif(N, 0, 1)
set.seed(6)
y <- runif(N, 0, 1)
coords <- cbind(x,y)

plot(coords)

# Se define una matriz de covariables
X <- cbind(rep(1,N),runif(N))

# Se fijan los parámetros del modelo de covarianza:
corrmodel <- "Exp"
mean <- 0.2
mean1 <- -0.5
sill <- 1
nugget <- 0
scale <- 0.2/3
param <- list(mean=mean,mean1=mean1,sill=sill,nugget=nugget,scale=scale)

# Simulación del RF espacial gausiano:
data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param,X=X)$data

# Parámetros fijos
fixed <- list(nugget=nugget)
# Valores iniciales para los parámetros a estimarse
start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill)

# Ajuste:
fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, 
                    maxdist=0.05,likelihood="Marginal",type="Pairwise",
                    start=start,fixed=fixed,X=X)
print(fit1)

## 
## ##################################################################
## Maximum  Composite-Likelihood Fitting of Gaussian Random Fields
## 
## Setting: Marginal Composite-Likelihood 
## 
## Model: Gaussian 
## 
## Type of the likelihood objects: Pairwise 
## 
## Covariance model: Exp 
## 
## Optimizer: Nelder-Mead 
## 
## Number of spatial coordinates: 400 
## Number of dependent temporal realisations: 1 
## Type of the random field: univariate 
## Number of estimated parameters: 4 
## 
## Type of convergence: Successful 
## Maximum log-Composite-Likelihood value: -1481.59
## 
## Estimated parameters:
##     mean     mean1     scale      sill  
##  0.16799  -0.41887   0.08109   1.10017  
## 
## ##################################################################

nn <- 50
vls <- seq(0,1,length.out = 50)
loc_to_pred <- as.matrix(expand.grid(vls,vls))
Xp <- cbind(rep(1,nn),runif(nn))

param_est <- as.list(c(fixed,fit1$param))
pr_geo <- GeoKrig(data=data,loc=loc_to_pred,coordx=coords,corrmodel=corrmodel,param=param_est,mse=TRUE,X=Xp,Xloc = loc_to_pred)

La matriz de locaciones, predicciones y error cuadrático medio puede obtenerse como:

predictions <- cbind(loc_to_pred,pr_geo$pred,pr_geo$mse)
head(predictions)

##            Var1 Var2                    
## [1,] 0.00000000    0 0.4310890 0.7245282
## [2,] 0.02040816    0 0.4435910 0.6234961
## [3,] 0.04081633    0 0.4500887 0.5383983
## [4,] 0.06122449    0 0.4742722 0.4503595
## [5,] 0.08163265    0 0.5265860 0.2832902
## [6,] 0.10204082    0 0.8254361 0.4093110

colour = rainbow(100)
image.plot(vls, vls, matrix(predictions[,3],ncol=length(vls)),col=colour,
           xlab="",ylab="",main="Simple Kriging con modelo Exponential")

image.plot(vls, vls, matrix(predictions[,4],ncol=length(vls)),col=colour,
           xlab="",ylab="",main="Std error")

Un ejemplo completo

Objetivo

Explicar cómo analizar datos esféricos (esto es, datos espaciales en la superficie de la esfera) con el paquete GeoModels

Denotamos con \(\mathbb{S}^2\) la esfera de \(\text{R}^3\) con radio arbitrario \(R>0\), definida como \[\mathbb{S}^2=\{ \boldsymbol{x}\in \text{R}^3: ||\boldsymbol{x}||=R \}.\]

Empezamos por cargar las librerías requeridas:

rm(list=ls())
library(GeoModels)
library(mapproj)
library(globe)
library(fields)
library(sphereplot)

## Loading required package: rgl

radius=1
set.seed(1891)

Modelos gaussianos definidos en la esfera

Consideramos un campo aleatorio gaussiano con media cero y varianza uno, \(Y=\{Y(\boldsymbol{x}), \boldsymbol{x} \in \mathbb{S}^2 \}\), definido en una esfera de radio \(R>0\). Un punto en \(\mathbb{S}^2\) puede ser parametrizado con coordenadas esféricas usando radianes como \((R,lon^*,lat^*)\), donde \(lon^* \in [-\pi,\pi]\) es la longitud y \(lat^* \in [-\pi/2,\pi/2]\) es la latitud. De forma alternativa, uno puede usar grados decimales \((R,lon,lat)\) donde \(lon \in [-180,180]\) y \(lat \in [-90,90]\) (el paquete GeoModels usa la segunda parametrización).

Se usan distancias como:

Euclideana: previo un cómputo de proyeccción
Ortodómica (The great circle distance): corresponde al camino más corto entre dos puntos en la superficie de la tierra y suele ser la forma natural de trabajo con datos en la esfera.
Distancia de Chordal (CH): es el segmento debajo del arco que une dos puntos para un par ubicado en la esfera.

Primero fijamos las coordenadas esféricas en grados decimales en una esfera unitaria con la función pointsphere del paquete sphereplot.

NN=1500 ## number of location sites on the earth 
coords=pointsphere(NN,c(-180,180),c(-90,90),c(radius,radius))[,1:2]

Como la esfera de radio \(1\) es una representación reescalada satisfactoria del planeta Tierra, podemos visualizar algunos puntos usando el paquete globe

globeearth(eye=place("newyorkcity"))
globepoints(loc=coords,pch=20,cex=0.4)

globeearth(eye=place("everest"))
globepoints(loc=coords,pch=20,cex=0.4)

Ahora usamos un modelo de covarianza Smoke definido como:

\[ \mathcal{F}_{1/\alpha,1/\alpha+0.5,\nu}(d_{GC} )= \frac{B(1/\alpha+0.5,\nu+1/\alpha)}{B(1/\alpha+0.5,\nu)} {}_{2}F_{1}(1/\alpha, 1/\alpha+0.5, 2/\alpha+0.5+\nu;\cos (d_{GC} )) \]

corrmodel = "Smoke" ## modelo de correlación y parámetros
scale = 0.3
smooth = 0.5
sill = 1
nugget = 0
mean = 0

Aquí el parámetro scale corresponde a \(\alpha\), sill es el parámetro de varianza \(\sigma^2\) y smooth es el parámetro de suavizamiento \(\nu\).

Estamos listos para simular los datos en la esfera de radio \(1\) usando la función GeoSim

param <- list(mean =mean, sill =sill, nugget =nugget, smooth =smooth,
scale = scale)
data <- GeoSim(coordx=coords, corrmodel=corrmodel, param=param,
             distance="Geod",sparse=TRUE,radius=radius)$data

Estimación

Ahora usamos la función GeoFit para hacer la estimación. En primer lugar, con verosimilitud completa.

fixed <-list(nugget = nugget )
start <-list(mean =mean , scale =scale , sill =sill , smooth = smooth )
fit_geo_ml <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, 
              likelihood="Full",type="Standard",
              start=start,fixed=fixed,distance="Geod",radius=1)

El cálculo de las distancias en el gran círculo que se usan en la matriz de covarianza puede ser obtenido usando la función rdist.earth del paquete fields.

geod_ds <- rdist.earth(coords,miles =F,R=1)
max_geod <- max(geod_ds)

También podemos hacer la estimación usando verosimilitud compuesta:

fit_geo <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, 
             likelihood="Marginal",type="Pairwise",maxdist=max_geod/40,
             start=start,fixed=fixed,distance="Geod",radius=radius)

Podemos ver los resultados de ambas estimaciones:

fit_geo_ml$param

##      mean     scale      sill    smooth 
## 0.1069231 0.2180537 0.7426665 0.4704625

fit_geo$param

##       mean      scale       sill     smooth 
## 0.08363007 0.26731324 0.77793523 0.43087207

El cálculo de la distancia de Chordal que se usa en la matriz de covarianza es:

chor_ds <- 2*sin(geod_ds/2)
max_chor <- max(chor_ds)

Se puede usar este tipo de distancia en GeoModels, lo hacemos con estimación de máxima verosimilitud en parejas para un modelo de correlación Matérn:

distance <- "Chor"
fit_chor <- GeoFit(data=data,coordx=coords,corrmodel="Matern", 
                 likelihood="Marginal",type="Pairwise",maxdist=max_chor/40,
                    start=start,fixed=fixed,distance=distance,radius=1)
print(fit_chor)

## 
## ##################################################################
## Maximum  Composite-Likelihood Fitting of Gaussian Random Fields
## 
## Setting: Marginal Composite-Likelihood 
## 
## Model: Gaussian 
## 
## Type of the likelihood objects: Pairwise 
## 
## Covariance model: Matern 
## 
## Optimizer: Nelder-Mead 
## 
## Number of spatial coordinates: 1500 
## Number of dependent temporal realisations: 1 
## Type of the random field: univariate 
## Number of estimated parameters: 4 
## 
## Type of convergence: Successful 
## Maximum log-Composite-Likelihood value: -1375.52
## 
## Estimated parameters:
##    mean    scale     sill   smooth  
## 0.09924  0.21283  0.79911  0.42123  
## 
## ##################################################################

También podemos considerar el modelo Matérn usando una proyección de los puntos en la esfera. Usamos la función mapprojec del paquete mapproj para obtener las coordenadas proyectadas usando una proyección sinusoidal:

prj <- mapproject(coords[,1], coords[,2], projection="sinusoidal")
coords_prj <- cbind(prj$x,prj$y)

Podemos visualizar los puntos proyectados:

sinusoidal.proj <- map(database= "world", ylim=c(-90,90), xlim=c(-180,180),
 col="grey80", fill=TRUE, plot=FALSE, projection="sinusoidal")
map(sinusoidal.proj)
points(coords_prj,pch=20,cex=0.4)

eucl_ds <- dist (coords_prj)
max_eucl <- max(eucl_ds)

distance <- "Eucl"
fit_eucl <- GeoFit(data =data ,coordx = coords_prj,corrmodel ="Matern",
maxdist = max_eucl/40, likelihood ="Marginal",type ="Pairwise",
start =start , fixed =fixed , distance ="Eucl",radius =1)

Veamos los resultados de estas estimaciones:

fit_chor$param

##       mean      scale       sill     smooth 
## 0.09924241 0.21282982 0.79911224 0.42122767

fit_eucl$param

##       mean      scale       sill     smooth 
## 0.08159125 0.13138696 0.75755533 0.49424735

Un gráfico de caja y bigotes puede usarse para mostrar el comportamiento de los diferentes tipos de distancias

boxplot (c(geod_ds),c( chor_ds),c(eucl_ds),
names = c(" Greatcircle "," Chordal "," Euclidean "))

Ahora podemos calcular los residuos:

res <- GeoResiduals(fit_geo_ml)

Con este cálculo, revisamos el supuesto de normalidad:

qqnorm(res$data)
abline(0,1)

También podemos comparar el variograma empírico y el calculado usando GeoVariogram y GeoCovariogram

vario <- GeoVariogram(data=res$data,coordx=coords,maxdist=max_geod/2,distance="Geod",radius = 1)
GeoCovariogram(res,vario=vario,show.vario=TRUE,pch=20)

Predicción

Primero necesitamos especificar los puntos donde se desea predecir. Supongamos que deseamos predecir en campo aleatorio gausseano alrededor del polo norte:

loc_to_pred <- as.matrix(expand.grid(seq(-180,180,2),seq(60,90,1)))

Obtenemos las predicciones como sigue:

param_est <- as.list(c(fixed,fit_geo_ml$param))
pr_geo <- GeoKrig(data=data,loc=loc_to_pred,coordx=coords,corrmodel=corrmodel,
            radius=radius,distance="Geod",param=param_est,mse=TRUE)

La matriz de locaciones, predicciones y error cuadrático medio puede obtenerse como:

predictions <- cbind(loc_to_pred,pr_geo$pred,pr_geo$mse)
head(predictions)

##      Var1 Var2                    
## [1,] -180   60 0.2254704 0.3290589
## [2,] -178   60 0.2180641 0.3483527
## [3,] -176   60 0.2137519 0.3569404
## [4,] -174   60 0.2186007 0.3531381
## [5,] -172   60 0.2319532 0.3356831
## [6,] -170   60 0.2494899 0.3061777

Finalmente, visualizamos los resultados:

globeearth(eye=place("newyorkcity"))
globepoints(loc=loc_to_pred,pch=20,col =  heat.colors(length(pr_geo$pred),alpha=0.1)[rank(pr_geo$pred)])

globeearth(eye=place("newyorkcity"))
globepoints(loc=loc_to_pred,pch=20,col =  cm.colors(length(pr_geo$mse),alpha=0.1)[rank(pr_geo$mse)])

GeoModels: una introducción

Víctor Morales Oñate
Sitio personal
ResearchGate
GitHub
LinkedIn

30 de octubre de 2020

Modelamiento

GeoModels

Un ejemplo completo

Objetivo

Modelos gaussianos definidos en la esfera

Estimación

Predicción

Algunas reflexiones

Referencias

GeoModels: una introducción

Víctor Morales Oñate Sitio personal ResearchGate GitHub LinkedIn

30 de octubre de 2020

Modelamiento

GeoModels

Un ejemplo completo

Objetivo

Modelos gaussianos definidos en la esfera

Estimación

Predicción

Algunas reflexiones

Referencias

Víctor Morales Oñate
Sitio personal
ResearchGate
GitHub
LinkedIn