GeoKrig: Spatial (bivariate) and spatio temporal optimal linear prediction for Gaussian and non Gaussian random fields.

Description

For a given set of spatial location sites (and temporal instants), the function computes optimal linear prediction and associated mean square error for the Gaussian and non Gaussian case.

Usage

GeoKrig(estobj=NULL,data, coordx, coordy=NULL, coordz=NULL, coordt=NULL, 
coordx_dyn=NULL, corrmodel,distance="Eucl",
    grid=FALSE, loc, 
    method="cholesky", model="Gaussian", n=1,nloc=NULL,mse=FALSE, 
     param, anisopars=NULL,radius=1, sparse=FALSE,
     time=NULL,  type_krig="Simple",weigthed=TRUE,which=1,
     copula=NULL, X=NULL,Xloc=NULL,Mloc=NULL,spobj=NULL,spdata=NULL,varcov=NULL)

Value

Returns an object of class Kg. An object of class Kg is a list containing at most the following components:

bivariate: TRUE if spatial bivariate cokriging is performed, otherwise FALSE;
coordx: A \(d\)-dimensional vector of spatial coordinates used for prediction;
coordy: A \(d\)-dimensional vector of spatial coordinates used for prediction;
coordz: A \(d\)-dimensional vector of spatial coordinates used for prediction;
coordt: A \(t\)-dimensional vector of temporal coordinates used for prediction;
corrmodel: String: the correlation model;
covmatrix: The covariance matrix.
data: The vector or matrix or array of data used for prediction
distance: String: the type of spatial distance;
grid: TRUE if the spatial data used for prediction are observed in a regular grid, otherwise FALSE;
loc: A (\(n \times 2\))-matrix of spatial locations to be predicted.
n: The number of trial for Binomial RFs
nozero: In the case of sparse simple kriging the percentage of non zero values in the covariance matrix. Otherwise is NULL.
numcoord: Numeric:he number \(d\) of spatial coordinates used for prediction;
numloc: Numeric: the number \(n\) of spatial coordinates to be predicted;
numtime: Numeric: the number \(d\) of the temporal instants used for prediction;
numt: Numeric: the number \(m\) of the temporal instants to be predicted;
model: The type of RF, see GeoFit.
param: Numeric: The covariance parameters;
pred: A (\(m \times n\))-matrix of spatio or spatio temporal kriging prediction;
radius: Numeric: the radius of the sphere if coordinates are pssed in lon/lat format;
spacetime: TRUE if spatio-temporal kriging and FALSE if spatial kriging;
time: A \(m\)-dimensional vector of temporal coordinates to be predicted;
type: String: the type of kriging (Standard).
type_krig: String: the type of kriging (simple or universal)
mse: A (\(m \times n\))-matrix of spatio or spatio temporal mean square error kriging prediction;

Arguments

estobj: An object of class Geofit that includes information about data, model and estimates.
data: A \(d\)-dimensional vector (a single spatial realisation) or a (\(d \times d\))-matrix (a single spatial realisation on regular grid) or a (\(t \times d\))-matrix (a single spatial-temporal realisation) or an (\(d \times d \times t \times n \))-array (a single spatial-temporal realisation on regular grid) giving the data used for prediction.
coordx: A numeric (\(d \times 2\))-matrix or (\(d \times 3\))-matrix Coordinates on a sphere for a fixed radius radius are passed in lon/lat format expressed in decimal degrees.
coordy: A numeric vector giving 1-dimension of spatial coordinates; Optional argument, the default is NULL.
coordz: A numeric vector giving 1-dimension of spatial coordinates; Optional argument, the default is NULL.
coordt: A numeric vector giving 1-dimension of temporal coordinates used for prediction; the default is NULL then a spatial random field is expected.
coordx_dyn: A list of \(m\) numeric (\(d_t \times 2\))-matrices containing dynamical (in time) spatial coordinates. Optional argument, the default is NULL
corrmodel: String; the name of a correlation model, for the description see the Section Details.
distance: String; the name of the spatial distance. The default is Eucl, the euclidean distance. See the Section Details of GeoFit.
grid: Logical; if FALSE (the default) the data used for prediction are interpreted as spatial or spatial-temporal realisations on a set of non-equispaced spatial sites (irregular grid).
loc: A numeric (\(n \times 2\))-matrix (where n is the number of spatial sites) giving 2-dimensions of spatial coordinates to be predicted.
method: String; the type of matrix decomposition used in the simulation. Default is cholesky. The other possible choices is svd.
n: Numeric; the number of trials in a binomial random fields. Default is \(1\).
nloc: Numeric; the number of trials of the locations sites to be predicted in a binomial random fields. Default is \(1\).
mse: Logical; if TRUE (the default) MSE of the kriging predictor is computed
model: String; the type of RF and therefore the densities associated to the likelihood objects. Gaussian is the default, see the Section Details.
param: A list of parameter values required for the correlation model.See the Section Details.
anisopars: A list of two elements: "angle" and "ratio" i.e. the anisotropy angle and the anisotropy ratio, respectively.
radius: Numeric: the radius of the sphere if coordinates are passed in lon/lat format;Default value is 1.
sparse: Logical; if TRUE kriging is computed with sparse matrices algorithms using spam package. Default is FALSE. It should be used with compactly supported covariances.
time: A numeric (\(m \times 1\)) vector (where m is the number of temporal instants) giving the temporal instants to be predicted; the default is NULL then only spatial prediction is performed.
type_krig: String; the type of kriging. If Simple (the default) then simple kriging is performed. Otherwise Universal.
weigthed: Logical; if TRUE then decreasing weigths coming from a compactly supported correlation function with compact support maxdist (maxtime)are used in the pairwise kriging.
which: Numeric; In the case of bivariate cokriging it indicates which variable to predict. It can be 1 or 2
copula: String; the type of copula. It can be "Clayton" or "Gaussian"
X: Numeric; Matrix of spatio(temporal)covariates in the linear mean specification.
Xloc: Numeric; Matrix of spatio(temporal)covariates in the linear mean specification associated to predicted locations.
Mloc: Numeric; Vector of spatio(temporal) estimated means associated to predicted locations.
spobj: An object of class sp or spacetime
spdata: Character:The name of data in the sp or spacetime object
varcov: Inverse of the covariance matrix estimates.

Author

Moreno Bevilacqua, moreno.bevilacqua89@gmail.com,https://sites.google.com/view/moreno-bevilacqua/home, Víctor Morales Oñate, victor.morales@uv.cl, https://sites.google.com/site/moralesonatevictor/, Christian", Caamaño-Carrillo, chcaaman@ubiobio.cl,https://www.researchgate.net/profile/Christian-Caamano

Details

Best linear unbiased predictor and associated mean square error is computed for Gaussian and some non Gaussian cases. Specifically, for a spatial or spatio-temporal or spatial bivariate dataset, given a set of spatial locations and temporal istants and a correlation model corrmodel with some fixed parameters and given the type of RF (model) the function computes simple kriging, for the specified spatial locations loc and temporal instants time, providing also the respective mean square error. For the choice of the spatial or spatio temporal correlation model see details in GeoCovmatrix function. The list param specifies mean and covariance parameters, see CorrParam and GeoCovmatrix for details. The type_krig parameter indicates the type of kriging. In the case of simple kriging, the known mean can be specified by the parameter mean in the list param (See examples).

References

Gaetan, C. and Guyon, X. (2010) Spatial Statistics and Modelling. Spring Verlang, New York.

Examples

Run this code


library(GeoModels)
################################################################
########### Examples of spatial kriging ############
################################################################

################################################################
###
### Example 1. Spatial  kriging of a
### Gaussian random fields with Gen wendland correlation.
###
################################################################

model="Gaussian"
set.seed(79)
x = runif(300, 0, 1)
y = runif(300, 0, 1)
coords=cbind(x,y)
# Set the exponential cov parameters:
corrmodel = "GenWend"
mean=0; sill=5; nugget=0
scale=0.2;smooth=0;power2=4

param=list(mean=mean,sill=sill,nugget=nugget,scale=scale,smooth=smooth,power2=power2)

# Simulation of the spatial Gaussian random field:
data = GeoSim(coordx=coords, corrmodel=corrmodel,model=model,
              param=param)$data

## estimation with pairwise likelihood
fixed=list(nugget=nugget,smooth=0,power2=power2)
start=list(mean=0,scale=scale,sill=1)
I=Inf
lower=list(mean=-I,scale=0,sill=0)
upper=list(mean= I,scale=I,sill=I)
# Maximum pairwise likelihood fitting :
fit = GeoFit(data, coordx=coords, corrmodel=corrmodel,model=model,
             likelihood='Marginal', type='Pairwise',neighb=3,
             optimizer="nlminb", lower=lower,upper=upper,
             start=start,fixed=fixed)

# locations to predict
xx=seq(0,1,0.03)
loc_to_pred=as.matrix(expand.grid(xx,xx))

## first option
#param=append(fit$param,fit$fixed)
#pr=GeoKrig(loc=loc_to_pred,coordx=coords,corrmodel=corrmodel,
#       model=model,param=param,data=data,mse=TRUE)

## second option using object GeoFit
pr=GeoKrig(fit,loc=loc_to_pred,mse=TRUE)


colour = rainbow(100)

opar=par(no.readonly = TRUE)
par(mfrow=c(1,3))
quilt.plot(coords,data,col=colour)
# simple kriging map prediction
image.plot(xx, xx, matrix(pr$pred,ncol=length(xx)),col=colour,
           xlab="",ylab="",
           main=" Kriging ")

# simple kriging MSE map prediction variance
image.plot(xx, xx, matrix(pr$mse,ncol=length(xx)),col=colour,
           xlab="",ylab="",main="Std error")
par(opar)

################################################################
###
### Example 2. Spatial  kriging of a Skew
### Gaussian random fields with Matern correlation.
###
################################################################
model="SkewGaussian"
set.seed(79)
x = runif(300, 0, 1)
y = runif(300, 0, 1)
coords=cbind(x,y)
# Set the exponential cov parameters:
corrmodel = "Matern"
mean=0
sill=2
nugget=0
scale=0.1
smooth=0.5
skew=3
param=list(mean=mean,sill=sill,nugget=nugget,scale=scale,smooth=smooth,skew=skew)

# Simulation of the spatial skew Gaussian random field:
data = GeoSim(coordx=coords, corrmodel=corrmodel,model=model,
              param=param)$data

fixed=list(nugget=nugget,smooth=smooth)
start=list(mean=0,scale=scale,sill=1,skew=skew)
I=Inf
lower=list(mean=-I,scale=0,sill=0,skew=-I)
upper=list(mean= I,scale=I,sill=I,skew=I)
# Maximum pairwise likelihood fitting :
fit = GeoFit2(data, coordx=coords, corrmodel=corrmodel,model=model,
              likelihood='Marginal', type='Pairwise',neighb=3,
              optimizer="nlminb", lower=lower,upper=upper,
              start=start,fixed=fixed)

# locations to predict
xx=seq(0,1,0.03)
loc_to_pred=as.matrix(expand.grid(xx,xx))
## optimal linear kriging
pr=GeoKrig(fit,loc=loc_to_pred,mse=TRUE)

colour = rainbow(100)

opar=par(no.readonly = TRUE)
par(mfrow=c(1,3))
quilt.plot(coords,data,col=colour)
# simple kriging map prediction
image.plot(xx, xx, matrix(pr$pred,ncol=length(xx)),col=colour,
           xlab="",ylab="",
           main=" Kriging ")

# simple kriging MSE map prediction variance
image.plot(xx, xx, matrix(pr$mse,ncol=length(xx)),col=colour,
           xlab="",ylab="",main="Std error")
par(opar)

################################################################
###
### Example 3. Spatial  kriging  of a 
###   Gamma random field with mean spatial regression
###
###############################################################
set.seed(312)
model="Gamma"
corrmodel = "GenWend"  
# Define the spatial-coordinates of the points:
NN=300
coords=cbind(runif(NN),runif(NN))
## matrix covariates
a0=rep(1,NN)
a1=runif(NN,0,1)
X=cbind(a0,a1)
##Set model parameters
shape=2
## regression parameters
mean = 1;mean1= -0.2
# correlation parameters
nugget = 0;power2=4
scale = 0.3;smooth=0    

## simulation
param=list(shape=shape,nugget=nugget,mean=mean,mean1=mean1, 
           scale=scale,power2=power2,smooth=smooth)
data = GeoSim(coordx=coords,corrmodel=corrmodel, param=param,
              model=model,X=X)$data

#####starting and fixed parameters
fixed=list(nugget=nugget,power2=power2,smooth=smooth)
start=list(mean=mean,mean1=mean1, scale=scale,shape=shape)

## estimation with pairwise likelihood
fit2 = GeoFit(data=data,coordx=coords,corrmodel=corrmodel,X=X,
              neighb=3,likelihood="Conditional",type="Pairwise",
              start=start,fixed=fixed, model = model)

# locations to predict with associated covariates
xx=seq(0,1,0.03)
loc_to_pred=as.matrix(expand.grid(xx,xx))
NP=nrow(loc_to_pred)
a0=rep(1,NP)
a1=runif(NP,0,1)
Xloc=cbind(a0,a1)

#optimal linear  kriging 
pr=GeoKrig(fit2,loc=loc_to_pred,Xloc=Xloc,sparse=TRUE,mse=TRUE)

## map 
opar=par(no.readonly = TRUE)
par(mfrow=c(1,3))
quilt.plot(coords,data,main="Data")
map=matrix(pr$pred,ncol=length(xx))
mapmse=matrix(pr$mse,ncol=length(xx))
image.plot(xx, xx, map,
           xlab="",ylab="",main="Kriging ")

image.plot(xx, xx, mapmse,
           xlab="",ylab="",main="MSE")
par(opar)


################################################################
########### Examples of spatio temporal kriging ############
################################################################

################################################################
###
### Example 4. Spatio temporal simple kriging of n locations
### sites and m temporal instants for a Gaussian random fields
### with estimated double Wendland correlation.
###
###############################################################
model="Gaussian"
# Define the spatial-coordinates of the points:
x = runif(300, 0, 1)
y = runif(300, 0, 1)
coords=cbind(x,y)
times=1:4

# Define model correlation modek and associated parameters
corrmodel="Wend0_Wend0"
param=list(nugget=0,mean=0,power2_s=4,power2_t=4,
           scale_s=0.2,scale_t=2,sill=1)

# Simulation of the space time Gaussian random field:
set.seed(31)
data=GeoSim(coordx=coords,coordt=times,corrmodel=corrmodel,sparse=TRUE,
            param=param)$data

# Maximum pairwise likelihood fitting of the space time random field:
start = list(scale_s=0.15,scale_t=2,sill=1,mean=0)
fixed = list(nugget=0,power2_s=4,power2_t=4)

fit = GeoFit(data, coordx=coords, coordt=times, model=model, corrmodel=corrmodel, 
             likelihood='Conditional', type='Pairwise',start=start,fixed=fixed,
             neighb=3,maxtime=1)

# locations to predict
xx=seq(0,1,0.04)
loc_to_pred=as.matrix(expand.grid(xx,xx))
#  Define the times to predict
times_to_pred=2

pr=GeoKrig(fit,loc=loc_to_pred,time=times_to_pred,sparse=TRUE,mse=TRUE)

opar=par(no.readonly = TRUE)
par(mfrow=c(1,3))
zlim=c(-2.5,2.5)
colour = rainbow(100)


quilt.plot(coords,data[2,] ,col=colour,main = paste(" data  at Time 2"))   
image.plot(xx, xx, matrix(pr$pred,ncol=length(xx)),col=colour,
           main = paste(" Kriging at Time 2"),ylab="")
image.plot(xx, xx, matrix(pr$mse,ncol=length(xx)),col=colour,
           main = paste("Std err Time at time 2"),ylab="")


par(opar)


################################################################
###
### Example r. Spatial bivariate  simple cokriging of n locations
### sites  for a bivariate Gaussian random fields
### with estimated Matern correlation.
###
###############################################################
#set.seed(6)
#NN=1500 # number of spatial locations 
#x = runif(NN, 0, 1); 
#y = runif(NN, 0, 1) 
#coords=cbind(x,y) 

## setting parameters
#mean_1 = 2; mean_2= -1
#nugget_1 =0;nugget_2=0
#sill_1 =0.5; sill_2 =1; 

### correlation parameters
#CorrParam("Bi_Matern")
#scale_1=0.2/3; scale_2=0.15/3; scale_12=0.5*(scale_2+scale_1) 
#smooth_1=smooth_2=smooth_12=0.5
#pcol = -0.4
#param= list(nugget_1=nugget_1,nugget_2=nugget_2,
#          sill_1=sill_1,sill_2=sill_2,
#          mean_1=mean_1,mean_2=mean_2,
#          smooth_1=smooth_1, smooth_2=smooth_2,smooth_12=smooth_12,
#          scale_1=scale_1, scale_2=scale_2,scale_12=scale_12,
#          pcol=pcol)

## simulation
#data = GeoSim(coordx=coords, corrmodel="Bi_Matern",model=model,param=param)$data

#fixed=list(mean_1=mean_1,mean_2=mean_2, nugget_1=nugget_1,nugget_2=nugget_2, 
#       smooth_1=smooth_1, smooth_2=smooth_2,smooth_12=smooth_12)

#start=list( sill_1=sill_1,sill_2=sill_2,
#        scale_1=scale_1,scale_2=scale_2,scale_12=scale_12, pcol=pcol)

## estimation with maximum likelihood 
#fit = GeoFit(data=data,coordx=coords, corrmodel="Bi_Matern",
  #likelihood="Marginal",type="Pairwise",optimizer="BFGS",neighb=5,
  #start=start,fixed=fixed)

###### co-kriging for the fist component ##############
#xx=seq(0,1,0.022)
#loc_to_pred=as.matrix(expand.grid(xx,xx))
#pr1 = GeoKrig(fit,which=1,mse=TRUE,loc=loc_to_pred)

#opar=par(no.readonly = TRUE)
#par(mfrow=c(1,2))
#zlim=c(-2.5,2.5)
#colour = rainbow(100)
#quilt.plot(coords,data[1,] ,col=colour,main = paste(" Fist component"))   
#quilt.plot(loc_to_pred,pr1$pred,col=colour,
#           main = paste(" Kriging first component"),ylab="")
#par(opar)

Run the code above in your browser using DataLab