Geo S Manual

Geos-manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 85

Package ‘GeoModels’

October 11, 2018

Version 1.0.3-4

Date 2018-01-15

Title A Package for Geostatistical Gaussian and non Gaussian Data

Analysis

Author Moreno Bevilacqua [aut, cre],

Víctor Morales-Onate [aut]

Maintainer Moreno Bevilacqua <moreno.bevilacqua@uv.cl>

Depends R (>= 2.12.0)

Description This package provides a set of procedures for a) simulation and estimation of some spa-

tial and spatio-temporal random ﬁelds using standard likelihood and a likelihood approxima-

tion method called composite likelihood and b) prediction using best linear unbiased prediction.

Spatio (temporal) bivariate data estimation involves estimation of both regression and covari-

ance parameters.

Gaussian and some non Gaussian Random ﬁelds can be analyzed using the GeoModels pack-

age. Among them gamma, Weibull, log-Gaussian, binomial, negative

binomial, skewGaussian, StudendT and circular random ﬁelds can be analyzed.

Imports methods

Suggests spam, scatterplot3d, ﬁelds, mapproj, gpuR, pbivnorm, plot3D,

shape, numDeriv, dfoptim, hypergeo, sphereplot,

License GPL (>= 2)

URL https://vmoprojs.github.io/GeoModels-page/

Repository GitHub

Encoding UTF-8

BugReports https://github.com/vmoprojs/GeoModels/issues

Rtopics documented:

anomalies .......................................... 2

CheckBiv .......................................... 3

CheckDistance ....................................... 3

CheckSph.......................................... 4

CheckST .......................................... 4

CkCorrModel........................................ 5

CkInput ........................................... 5

CkLikelihood ........................................ 7

2anomalies

CkModel .......................................... 7

CkType ........................................... 8

CkVarType ......................................... 8

CompLik .......................................... 9

CorrelationPar........................................ 11

CorrParam.......................................... 11

DeviceInfo ......................................... 12

GeoCovariogram ...................................... 13

GeoCovmatrix........................................ 18

GeoFit............................................ 29

GeoKrig........................................... 43

GeoResiduals ........................................ 53

GeoSim ........................................... 54

GeoTests .......................................... 61

GeoVarestbootstrap..................................... 64

GeoVariogram........................................ 65

GeoWLS .......................................... 70

Lik.............................................. 73

MatDecomp......................................... 75

MatSqrt,MatInv,MatLogDet................................ 76

NuisParam.......................................... 77

Prscores........................................... 78

StartParam.......................................... 79

winds ............................................ 81

winds.coords ........................................ 81

WlsStart........................................... 82

Index 84

anomalies Annual precipitation anomalies in U.S.

Description

A (7252x3)-matrix containing lon/lat and yearly total precipitation anomalies registered at 7.352 lo-

cation sites in USA. For more details see http://www.image.ucar.edu/Data/precip_tapering/.

Usage

data(anomalies)

Format

A numerical matrix of dimension 7252x3.

Source

Kaufman, C.G., Schervish, M.J., Nychka, D.W. (2008) Covariance tapering for likelihood-based

estimation in large spatial data sets. Journal of the American Statistical Association, Theory &

Methods,103, 1545–1555.

CheckBiv 3

CheckBiv Checking Bivariate covariance models

Description

The procedure control if the correlation model is bivariate.

Usage

CheckBiv(numbermodel)

Arguments

numbermodel numeric; the number associated to a given correlation model.

Value

Returns TRUE or FALSE depending if the correlation model is bivariate or not.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

Examples

library(GeoModels)

CheckBiv(CkCorrModel("Bi_matern_sep"))

CheckDistance Checking Distance

Description

The procedure controls the type of distance.

Usage

CheckDistance(distance)

Arguments

distance String; the type of distance, for the description see GeoCovmatrix. Default is

Eucl. Other possible values are Geod and Chor that is euclidean, geodesic and

chordal distance.

Value

Returns 0,1,2 for euclidean,geodesic, chordal distances respectively. Otherwise returns NULL.

4CheckST

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CheckSph Checking if a covariance is valid only on sphere

Description

Subroutine called by InitParam. The procedure controls if a covariance model is valid only on the

sphere.

Usage

CheckSph(numbermodel)

Arguments

numbermodel Numeric; the code number for the covariance model.

Value

Returns TRUE or FALSE

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

CheckST Checking SpaceTime covariance models

Description

The procedure control if the correlation model is spacetime.

Usage

CheckST(numbermodel)

Arguments

numbermodel numeric; the number associated to a given correlation model.

Value

Returns TRUE or FALSE depending if the correlation model is spatiotemporal or not.

CkCorrModel 5

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

Examples

library(GeoModels)

CheckST(CkCorrModel("gneiting"))

CkCorrModel Checking Correlation Model

Description

The procedure controls if the correlation model inserted is correct.

Usage

CkCorrModel(corrmodel)

Arguments

corrmodel String; the name of a correlation model, for the description see GeoCovmatrix.

Value

Return a number associated to a given correlation model if the model is considered in the package.

Otherwise return NULL.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

CkInput Checking Input

Description

Subroutine called by the ﬁtting procedures. The procedure controls the the validity of the input

inserted by the users.

Usage

CkInput(coordx, coordy, coordt, coordx_dyn,corrmodel, data, distance,

fcall, fixed, grid,likelihood, maxdist, maxtime,

model, n, optimizer, param, radius,

start, taper, tapsep, type, varest, vartype,

weighted,X)

6CkInput

Arguments

coordx A numeric (d×2)-matrix (where dis the number of points) assigning 2-dimensions

of coordinates or a numeric vector assigning 1-dimension of coordinates.

coordy A numeric vector assigning 1-dimension of coordinates; coordy is interpreted

only if coordx is a numeric vector otherwise it will be ignored.

coordt A numeric vector assigning 1-dimension of temporal coordinates.

corrmodel String; the name of a correlation model, for the description see GeoFit.

coordx_dyn A list of mnumeric (dt×2)-matrices containing dynamical (in time) spatial

coordinates. Optional argument, the default is NULL

data A numeric vector or a (n×d)-matrix or (d×d×n)-matrix of observations.

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See the Section Details.

fcall String; Fitting to call the ﬁtting procedure and simulation to call the simula-

tion.

fixed A named list giving the values of the parameters that will be considered as

known values. The listed parameters for a given correlation function will be

not estimated, i.e. if list(nugget=0) the nugget effect is ignored.

grid Logical; if FALSE (the default) the data are interpreted as a vector or a (n×d)-

matrix, instead if TRUE then (d×d×n)-matrix is considered.

likelihood String; the conﬁguration of the composite likelihood. Marginal is the default.

maxdist Numeric; an optional positive value indicating the maximum spatial distance

considered in the composite-likelihood computation.

maxtime Numeric; an optional positive value indicating the maximum temporal lag sepa-

ration in the composite-likelihood.

radius Numeric; the radius of the sphere in the case of lon-lat coordinates. The default

is 6371, the radius of the earth.

model String; the density associated to the likelihood objects. Gaussian is the default.

nNumeric; the number of trials in a binomial random ﬁelds. Default is 1.

optimizer String; the optimization algorithm (see optim for details). ’Nelder-Mead’ is the

default.

param A numeric vector of parameters, needed only in simulation. See GeoSim.

start A named list with the initial values of the parameters that are used by the nu-

merical routines in maximization procedure. NULL is the default.

taper String; the name of the tapered correlation function.

tapsep Numeric; an optional value indicating the separabe parameter in the space time

quasi taper (see Details).

type String; the type of the likelihood objects. If Pairwise (the default) then the

marginal composite likelihood is formed by pairwise marginal likelihoods.

varest Logical; if TRUE the estimate’ variances and standard errors are returned. FALSE

is the default.

vartype String; the type of estimation method for computing the estimate variances, see

GeoFit.

weighted Logical; if TRUE the likelihood objects are weighted. If FALSE (the default) the

composite likelihood is not weighted.

XNumeric; Matrix of space-time covariates in the linear mean speciﬁcation.

CkLikelihood 7

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CkLikelihood Checking Composite-likelihood Type

Description

Subroutine called by InitParam. The procedure controls the type of the composite-likelihood in-

serted by the users.

Usage

CkLikelihood(likelihood)

Arguments

likelihood String; the conﬁguration of the composite likelihood. Marginal is the default.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CkModel Checking Random Field type

Description

Subroutine called by InitParam. The procedure controls the type of random ﬁeld inserted by the

users.

Usage

CkModel(model)

Arguments

model String; the density associated to the likelihood objects. Gaussian is the default.

8CkVarType

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CkType Checking Likelihood Objects

Description

Subroutine called by InitParam. The procedure controls the type of likelihood objects inserted by

the users.

Usage

CkType(type)

Arguments

type String; the type of the likelihood objects. If Pairwise (the default) then the

marginal composite likelihood is formed by pairwise marginal likelihoods.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CkVarType Checking Variance Estimates Type

Description

Subroutine called by InitParam. The procedure controls the method used to compute the estimates’

variances.

Usage

CkVarType(type)

Arguments

type String; the method used to compute the estimates’ variances. If SubSamp the

estimates’ variances are computed by the sub-sampling method, see GeoFit.

CompLik 9

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CompLik Optimizes the Composite log-likelihood

Description

Subroutine called by GeoFit. The procedure estimates the model parameters by maximisation of

the composite log-likelihood.

Usage

CompLik(bivariate, coordx, coordy ,coordt, coordx_dyn, corrmodel, data, distance,

flagcorr, flagnuis, fixed, GPU, grid,likelihood,local,lower,

model, n, namescorr, namesnuis, namesparam,

numparam, numparamcorr, optimizer, onlyvar, param,

spacetime, type,upper, varest, vartype,

weigthed, winconst, winstp,winconst_t, winstp_t, ns, X,sensitivity)

Arguments

bivariate Logical; if TRUE then the data come froma a bivariate random ﬁeld. Otherwise

from a univariate random ﬁeld.

coordx A numeric (d×2)-matrix (where dis the number of points) assigning 2-dimensions

of coordinates or a numeric vector assigning 1-dimension of coordinates.

coordy A numeric vector assigning 1-dimension of coordinates; coordy is interpreted

only if coordx is a numeric vector otherwise it will be ignored.

coordt A numeric vector assigning 1-dimension of temporal coordinates. Optional ar-

gument, the default is NULL then a spatial random ﬁeld is expected.

coordx_dyn A list of mnumeric (dt×2)-matrices containing dynamical (in time) spatial

coordinates. Optional argument, the default is NULL

corrmodel Numeric; the id of the correlation model.

data A numeric vector or a (n×d)-matrix or (d×d×n)-matrix of observations.

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See the Section Details.

flagcorr A numeric vector of binary values denoting which paramerters of the correlation

function will be estimated.

flagnuis A numeric vector of binary values denoting which nuisance paramerters will be

estimated.

fixed A numeric vector of parameters that will be considered as known values.

GPU Numeric; if NULL (the default) no GPU computation is performed.

10 CompLik

grid Logical; if FALSE (the default) the data are interpreted as a vector or a (n×d)-

matrix, instead if TRUE then (d×d×n)-matrix is considered.

likelihood String; the conﬁguration of the compositelikelihood, see GeoFit.

local Numeric; number of local work-items of the GPU

lower A numeric vector with the lower bounds of the parameters’ ranges.

model Numeric; the id value of the density associated to the likelihood objects.

nNumeric; number of trials in a binomial random ﬁelds.

namescorr String; the names of the correlation parameters.

namesnuis String; the names of the nuisance parameters.

namesparam String; the names of the parameters to be maximised.

numparam Numeric; the number of parameters to be maximised.

numparamcorr Numeric; the number of correlation parameters.

optimizer String; the optimization algorithm (see optim for details). ’Nelder-Mead’ is the

default.

onlyvar Logical; if TRUE (and varest is TRUE) only the variance covariance matrix is

computed without optimizing. FALSE is the default.

param A numeric vector of parameters’ values.

spacetime Logical; if TRUE the random ﬁeld is spatial-temporal otherwise is a spatial ﬁeld.

type String; the type of the likelihood objects. If Pairwise (the default) then the

marginal composite likelihood is formed by pairwise marginal likelihoods.

upper A numeric vector with the upper bounds of the parameters’ ranges.

varest Logical; if TRUE the estimate’ variances and standard errors are returned. FALSE

is the default.

vartype String; the type of estimation method for computing the estimate variances, see

GeoFit.

weigthed Logical; if TRUE then decreasing weigths coming from a compactly supported

correlation function with compact support maxdist (maxtime)are used.

winconst Numeric; a positive value for computing the spatial sub-window in the sub-

sampling procedure.

winstp Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

spatial sub-sampling procedure.

winconst_t Numeric; a positive value for computing the temporal sub-window in the sub-

sampling procedure.

winstp_t Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

temporal sub-sampling procedure.

ns Numeric; Number of (dynamical) temporal instants.

XNumeric; Matrix of space-time covariates in the linear mean speciﬁcation.

sensitivity Logical; if TRUE then the sensitivy matrix is computed

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CorrelationPar 11

CorrelationPar Lists the Parameters of a Correlation Model

Description

Subroutine called by InitParam and other procedures. The procedure returns a list with the param-

eters of a given correlation model.

Usage

CorrelationPar(corrmodel)

Arguments

corrmodel Integer; an integer associated to a given correlation model.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

CorrParam Lists the Parameters of a Correlation Model

Description

Subroutine called by InitParam and other procedures. The procedure returns a list with the param-

eters of a given correlation model.

Usage

CorrParam(corrmodel)

Arguments

corrmodel String; the name of a correlation model.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoCovmatrix

12 DeviceInfo

Examples

require(GeoModels)

################################################################

###

### Example 1. Parameters of the Matern model

###

###############################################################

CorrParam("Matern")

################################################################

###

### Example 2. Parameters of the Generalized Wendland model

###

###############################################################

CorrParam("GenWend")

################################################################

###

### Example 3. Parameters of the Generalized Wendland model

###

###############################################################

CorrParam("GenCauchy")

################################################################

###

### Example 4. Parameters of the space time Gneiting model

###

###############################################################

CorrParam("Gneiting")

################################################################

###

### Example 5. Parameters of the bi-Matern separable model

###

###############################################################

CorrParam("Bi_Matern_sep")

DeviceInfo Prints Device Information

Description

Prints the device details available in your computer. Device name, Max compute units, whether it

supports double precision, among others.

GeoCovariogram 13

Usage

DeviceInfo()

Details

The user can take this information into account so that the local parameter is set up in GeoFit

when GPU computation is chosen.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

Examples

library(GeoModels)

DeviceInfo()

GeoCovariogram Computes ﬁtted covariance and/or variogram

Description

The procedure computes and plots covariance or variogram estimated ﬁtting a Gaussian, and non

Gaussian spatio (temporal) bivariate random ﬁelds. Allows to add the empirical estimates in order

to compare them with the ﬁtted model.

Usage

GeoCovariogram(fitted, distance="Eucl",answer.cov=FALSE,

answer.vario=FALSE, answer.range=FALSE, fix.lags=NULL,

fix.lagt=NULL, show.cov=FALSE, show.vario=FALSE,

show.range=FALSE, add.cov=FALSE, add.vario=FALSE,

pract.range=95, vario, ...)

Arguments

fitted A ﬁtted object obtained from the GeoFit or GeoWLS procedures.

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See GeoFit.

answer.cov Logical; if TRUE a vector with the estimated covariance function is returned; if

FALSE (the default) the covariance is not returned.

answer.vario Logical; if TRUE a vector with the estimated variogram is returned; if FALSE (the

default) the variogram is not returned.

answer.range Logical; if TRUE the estimated pratical range is returned; if FALSE (the default)

the pratical range is not returned.

fix.lags Integer; a positive value denoting the spatial lag to consider for the plot of the

temporal proﬁle.

14 GeoCovariogram

fix.lagt Integer; a positive value denoting the temporal lag to consider for the plot of the

spatial proﬁle.

show.cov Logical; if TRUE the estimated covariance function is plotted; if FALSE (the de-

fault) the covariance function is not plotted.

show.vario Logical; if TRUE the estimated variogram is plotted; if FALSE (the default) the

variogram is not plotted.

show.range Logical; if TRUE the estimated pratical range is added on the plot; if FALSE (the

default) the pratical range is not added.

add.cov Logical; if TRUE the vector of the estimated covariance function is added on the

current plot; if FALSE (the default) the covariance is not added.

add.vario Logical; if TRUE the vector with the estimated variogram is added on the current

plot; if FALSE (the default) the correlation is not added.

pract.range Numeric; the percent of the sill to be reached.

vario AVariogram object obtained from the GeoVariogram procedure.

... other optional parameters which are passed to plot functions.

Value

The returned object is eventually a list with:

covariance The vector of the estimated covariance function;

variogram The vector of the estimated variogram function;

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

References

Cressie, N. A. C. (1993) Statistics for Spatial Data. New York: Wiley.

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

See Also

GeoFit.

Examples

library(GeoModels)

library(scatterplot3d)

################################################################

###

### Example 1. Plot of fitted covariance and fitted

### and empirical variograms from a Gaussian RF

### with Matern correlation.

###

###############################################################

set.seed(21)

# Set the coordinates of the points:

x <- runif(300, 0, 1)

GeoCovariogram 15

y <- runif(300, 0, 1)

coords=cbind(x,y)

# Set the model's parameters:

corrmodel <- "Matern"

model <- "Gaussian"

mean <- 0

sill <- 1

nugget <- 0

scale <- 0.2/3

smooth=0.5

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

# Simulation of the Gaussian random field:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, model=model,param=param)$data

start=list(mean=0,scale=scale,sill=sill)

fixed=list(nugget=nugget,smooth=smooth)

# Maximum composite-likelihood fitting of the Gaussian random field:

fit <- GeoFit(data=data,coordx=coords, corrmodel=corrmodel,model=model,

likelihood="Marginal",type='Pairwise',start=start,

fixed=fixed,maxdist=0.05)

# Empirical estimation of the variogram:

vario <- GeoVariogram(data=data,coordx=coords,maxdist=0.5)

# Plot of covariance and variogram functions:

GeoCovariogram(fit, show.cov=TRUE,show.vario=TRUE, vario=vario,pch=20)

################################################################

###

### Example 2. Plot of fitted covariance and fitted

### and empirical variograms from a Binomial

### RF with exponential correlation.

###

###############################################################

set.seed(2111)

model="Binomial";n=20

# Set the coordinates of the points:

x <- runif(500, 0, 1)

y <- runif(500, 0, 1)

coords=cbind(x,y)

# Set the model's parameters:

corrmodel <- "exponential"

mean <- 0

sill <- 1

nugget <- 0

scale <- 0.2/3

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

# Simulation of the Gaussian RF:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, model=model,param=param,n=n)$data

start=list(mean=0,scale=scale,sill=sill)

fixed=list(nugget=nugget)

16 GeoCovariogram

# Maximum composite-likelihood fitting of the BinomGaussian random field:

fit <- GeoFit(data,coordx=coords, corrmodel=corrmodel,model=model,

likelihood="Marginal",type='Pairwise',start=start,n=n,

fixed=fixed,maxdist=0.03)

# Empirical estimation of the variogram:

vario <- GeoVariogram(data,coordx=coords,maxdist=0.5)

# Plot of covariance and variogram functions:

GeoCovariogram(fit, show.cov=TRUE,show.vario=TRUE, vario=vario,pch=20)

################################################################

###

### Example 3. Plot of fitted covariance and fitted

### and empirical variograms from a RF

### RF with Wend0 correlation.

###

###############################################################

set.seed(211)

model="Gamma";shape=4

# Set the coordinates of the points:

x <- runif(700, 0, 1)

y <- runif(700, 0, 1)

coords=cbind(x,y)

# Set the model's parameters:

corrmodel <- "Wend0"

mean <- 0

sill <- 1

nugget <- 0

scale <- 0.3

power2=4

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

# Simulation of the Gaussian RF:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, model=model,param=param)$data

start=list(mean=0,scale=scale,shape=shape)

fixed=list(nugget=nugget,sill=sill,power2=power2)

# Maximum composite-likelihood fitting of the BinomGaussian random field:

fit <- GeoFit(data,coordx=coords, corrmodel=corrmodel,model=model,

likelihood="Marginal",type='Pairwise',start=start,

fixed=fixed,maxdist=0.03)

# Empirical estimation of the variogram:

vario <- GeoVariogram(data,coordx=coords,maxdist=0.5)

# Plot of covariance and variogram functions:

GeoCovariogram(fit, show.cov=TRUE,show.vario=TRUE, vario=vario,pch=20)

################################################################

###

### Example 4. Plot of fitted and empirical variograms

### from a space time Gaussian random fields

GeoCovariogram 17

### with double exponential correlation.

###

###############################################################

set.seed(92)

# Define the spatial-coordinates of the points:

x <- runif(50, 0, 1)

y <- runif(50, 0, 1)

coords=cbind(x,y)

# Define the temporal sequence:

time <- seq(0, 15, 1)

# Simulation of the spatio-temporal Gaussian random field:

data <- GeoSim(coordx=coords, coordt=time, corrmodel="Exp_Exp",param=list(mean=mean,

nugget=nugget,scale_s=0.5/3,scale_t=2/2,sill=sill))$data

fixed=list(nugget=0, mean=0)

start=list(scale_s=0.2, scale_t=0.5, sill=1)

# Maximum composite-likelihood fitting of the space-time Gaussian random field:

fit <- GeoFit(data, coordx=coords, coordt=time, corrmodel="Exp_Exp", maxtime=2,

maxdist=0.1, likelihood="Marginal", type="Pairwise",

fixed=fixed, start=start)

# Empirical estimation of spatio-temporal covariance:

vario <- GeoVariogram(data,coordx=coords, coordt=time, maxtime=5,maxdist=0.5)

# Plot of the fitted space-time variogram

GeoCovariogram(fit,vario=vario,show.vario=TRUE)

# Plot of covariance, variogram and spatio and temporal profiles:

GeoCovariogram(fit,vario=vario,fix.lagt=1,fix.lags=1,show.vario=TRUE,pch=20)

################################################################

###

### Example 5. Plot of parametric and empirical variograms

### estimated from a Bivariate Gaussian random fields with

### Matern correlation.

###

###############################################################

# Simulation of a bivariate spatial Gaussian random field:

set.seed(892)

# Define the spatial-coordinates of the points:

x <- runif(200, -1, 1)

y <- runif(200, -1, 1)

coords=cbind(x,y)

# Simulation of a bivariate Gaussian Random field

# with matern (cross) covariance function

scale_1 = 0.25/3

scale_2 = 0.2/3

scale_12 = 0.15/3

sill_1=1

sill_2=1

18 GeoCovmatrix

smooth=0.5

pcol=0.3

param=list(mean_1=0,mean_2=0,scale_1=scale_1,scale_2=scale_2,scale_12=scale_12,

sill_1=sill_1,sill_2=sill_2,nugget_1=0,nugget_2=0,

smooth_1=smooth,smooth_12=smooth,smooth_2=smooth,pcol=pcol)

data <- GeoSim(coordx=coords, corrmodel="Bi_Matern", param=param)$data

# Empirical bivariate variogram estimation:

biv_vario=GeoVariogram(data,coordx=coords, bivariate=TRUE,maxdist=c(1,1,1))

# selecting fixed and estimating parameters

fixed=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,

smooth_1=smooth,smooth_12=smooth,smooth_2=smooth)

start=list(sill_1=var(data[1,]),sill_2=var(data[2,]),

scale_1=scale_1,scale_2=scale_2,scale_12=scale_12,

pcol=cor(data[1,],data[2,]))

# Maximum likelihood fitting of the bivariate random field:

fit<- GeoFit(data, coordx=coords, corrmodel="Bi_Matern",likelihood="Marginal",

type="Pairwise",start=start,fixed=fixed,maxdist=c(0.1,0.1,0.1))

GeoCovariogram(fit, vario=biv_vario,show.vario=TRUE,pch=20)

GeoCovmatrix Spatial and Spatio-temporal (tapered) Covariance Matrix

Description

The function computes the covariance matrix associated to a spatial or spatio-temporal or a bivariate

spatial Gaussian or non Gaussian randomm ﬁeld with given covariance model and a set of spatial

location sites and temporal instants.

Usage

GeoCovmatrix(coordx, coordy=NULL, coordt=NULL, coordx_dyn=NULL, corrmodel,

distance="Eucl", grid=FALSE, maxdist=NULL, maxtime=NULL,

model="Gaussian", n=1, param, radius=6371, sparse=FALSE,

taper=NULL, tapsep=NULL, type="Standard",X=NULL)

Arguments

coordx A numeric (d×2)-matrix (where dis the number of spatial sites) giving 2-

dimensions of spatial coordinates or a numeric d-dimensional vector giving 1-

dimension of spatial coordinates. Coordinates on a sphere for a ﬁxed radius

radius are passed in lon/lat format expressed in decimal degrees.

coordy A numeric vector giving 1-dimension of spatial coordinates; coordy is inter-

preted only if coordx is a numeric vector or grid=TRUE otherwise it will be

ignored. Optional argument, the default is NULL then coordx is expected to be

numeric a (d×2)-matrix.

GeoCovmatrix 19

coordt A numeric vector giving 1-dimension of temporal coordinates. At the moment

implemented only for the Gaussian case. Optional argument, the default is NULL

then a spatial random ﬁeld is expected.

coordx_dyn A list of Tnumeric (dt×2)-matrices containing dynamical (in time) coordinates.

Optional argument, the default is NULL

corrmodel String; the name of a correlation model, for the description see the Section De-

tails.

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See GeoFit.

grid Logical; if FALSE (the default) the data are interpreted as spatial or spatial-

temporal realisations on a set of non-equispaced spatial sites (irregular grid).

See GeoFit.

maxdist Numeric; an optional positive value indicating the marginal spatial compact sup-

port in the case of tapered covariance matrix. See GeoFit.

maxtime Numeric; an optional positive value indicating the marginal temporal compact

support in the case of spacetime tapered covariance matrix. See GeoFit.

nNumeric; the number of trials in a binomial random ﬁelds. Default is 1.

model String; the type of RF. See GeoFit.

param A list of parameter values required for the covariance model.

radius Numeric; a value indicating the radius of the sphere when using covariance

models valid using the great circle distance. Default value is the radius of the

earth in Km (i.e. 6371)

sparse Logical; if TRUE the function return an object of class spam. This option should

be used when a parametric compactly supporte covariance is used. Default is

FALSE.

taper String; the name of the taper correlation function if type is Tapering, see the

Section Details.

tapsep Numeric; an optional value indicating the separabe parameter in the space-time

non separable taper or the colocated correlation parameter in a bivariate spatial

taper (see Details).

type String; the type of covariance matrix Standard (the default) or Tapering for

tapered covariance matrix.

XNumeric; Matrix of space-time covariates.

Details

In the spatial case, the covariance matrix of the random vector

[Z(s1), . . . , Z(sn,)]T

with a speciﬁc spatial covariance model is computed. Here nis the number of the spatial location

sites.

In the space-time case, the covariance matrix of the random vector

[Z(s1, t1), Z(s2, t1), . . . , Z(sn, t1), . . . , Z(sn, tm)]T

with a speciﬁc space time covariance model is computed. Here mis the number of temporal in-

stants.

20 GeoCovmatrix

In the bivariate case, the covariance matrix of the random vector

[Z1(s1), Z2(s1), . . . , Z1(sn), Z2(sn)]T

with a speciﬁc spatial bivariate covariance model is computed.

The location site sican be a point in the d-dimensional euclidean space with d= 2 or a point (given

in lon/lat degree format) on a sphere of arbitrary radius.

Here there is the list of all the implemented space and space-time and bivariate correlation models.

The argument param is a list including all the parameters of a given correlation model speciﬁed

by the argument corrmodel. For each correlation model one can check the associated correlation

parameters using CorrParam. In what follows κ > 0,β > 0,α, αs, αt∈(0,2], and γ∈[0,1]. The

associated parameters in the argument param are smooth,power2,power,power_s,power_t and

sep respectively. Moreover let 1(A)=1when Ais true and 0otherwise.

• Spatial correlation models:

1. Cauchy deﬁned as:

R(h) = 1 + h2−β/2

It is a special case of the Gencauchy model.

2. Exp deﬁned as:

R(h) = e−h

This model is a special case of the Matern and the Stable model.

3. Gauss deﬁned as:

R(h) = e−h2

This model is a special case of the stable model.

4. GenCauchy (generalised Cauchy) deﬁned as:

R(h) = (1 + hα)−β/α

If his the geodesic distance then α∈(0,1].

5. Matern deﬁned as:

R(h) = 21−κΓ(κ)−1hκKκ(h)

If his the geodesic distance then κ∈(0,0.5]

6. Stable deﬁned as:

R(h) = e−hα

If his the geodesic distance then α∈(0,1].

7. W ave deﬁned as:

R(h) = sin(h)/h

This model is valid only for dimensions less than or equal to 3.

8. W end0deﬁned as:

R(h) = (1 −h)µ1(h∈[0,1])

where µ≥0.5(d+ 1). If his the geodesic distance then µ≥2.

9. W end1deﬁned as:

R(h) = (1 −h)µ+1(1 + (µ+ 1)h)1(h∈[0,1])

where µ≥0.5(d+ 1) + 1. If his the geodesic distance then µ≥4.

GeoCovmatrix 21

10. W end2deﬁned as:

R(h) = (1 −h)µ+2(1 + (µ+ 2)h+ (1/3)((µ+ 1)2−1)h2)1(h∈[0,1])

where µ≥0.5(d+ 1) + 2. If his the geodesic distance then µ≥6.

11. GenW end (Generalized Wendland) deﬁned as:

R(h) = Z1

[(1 −x)µ−1(x2−h2)κ−11(h∈[0,1])]dx/B(2κ+ 1, µ)

where µ≥0.5(d+ 1) + κ. The cases κ= 0,1,2correspond to the W end0,W end1and

W end2respectively.

12. Multiquadric deﬁned as:

R(h) = (1 −α0.5)2β/(1 + (α0.5)2−αcos(h))β, h ∈[0, π]

This model is valid on the unit sphere and his the geodesic distance.

13. Sinpower deﬁned as:

R(h)=1−(sin(h/2))α, h ∈[0, π]

This model is valid on the unit sphere and his the geodesic distance.

14. Smoke deﬁned as:

R(h) = K∗1F2(1/α, 1/α + 0.5,2/α + 0.5 + κ), h ∈[0, π]

where K= (Γ(a)Γ(i))/Γ(i)Γ(o)). This model is valid on the unit sphere and his the

geodesic distance.

• Spatio-temporal correlation models.

–Non-separable models:

1. Gneiting deﬁned as:

R(h, u) = e−hαs/((1+uαt)0.5γαs)/(1 + uαt)

2. Gneiting_GC

R(h, u) = e−uαt/((1+hαs)0.5γαt)/(1 + hαs)

where hcan be both the euclidean and the geodesic distance

3. Iacocesare

R(h, u) = (1 + hαs+uα

t)−β

4. P orcu

R(h, u) = (0.5(1 + hαs)γ+ 0.5(1 + uαt)γ)−γ−1

5. P orcu1

R(h, u)=(e−hαs(1+uαt)0.5γαs)/((1 + uαt)1.5)

6. Stein

R(h, u)=(hψ(u)Kψ(u)(h))/(2ψ(u)Γ(ψ(u) + 1))

where ψ(u) = ν+u0.5αt

7. W enx_space,x= 0,1,2deﬁned as:

R(h, u) = φ(u)3.5+2xW enx(h/φ(u), µs), x = 0,1,2

where φ(u) = (1 + u0.5αt)−γ,0< γ ≤αt/2,µs≥0.5(d+ 5) + x.

22 GeoCovmatrix

8. W enx_time,x= 0,1,2deﬁned as:

R(h, u) = φ(h)3.5+2xW enx(u/φ(h); µt), x = 0,1,2

where φ(h) = (1 + h0.5αs)−γ,0< γ ≤αs/2,µt≥0.5(d+ 5) + x.

9. Multiquadric_st deﬁned as:

R(h, u) = ((1 −0.5αs)2/(1 + (0.5αs)2−αsψ(u)cos(h)))as, h ∈[0, π]

where ψ(u) = (1 + (u/at)αt)−1. This model is valid on the unit sphere and his the

geodesic distance.

10. Sinpower_st deﬁned as:

R(h, u) = (eαscos(h)ψ(u)/as(1 + αscos(h)ψ(u)/as))/k

where ψ(u) = (1+(u/at)αt)−1and k= (1+αs/as)exp(αs/as), h ∈[0, π]This

model is valid on the unit sphere and his the geodesic distance.

–Separable models.

Space-time separable correlation models are easly obtained as the product of a spatial and

a temporal correlation model, that is

R(h, u) = R(h)R(u)

Several combinations are possible:

1. Exp_Exp deﬁned as:

R(h, u) = Exp(h)Exp(u)

2. Matern_Matern deﬁned as:

R(h, u) = Matern(h;κs)Matern(u;κt)

3. Stable_Stable deﬁned as:

R(h, u) = Stable(h;αs)Stable(u;αt)

4. W endx_W endy deﬁned as

R(h, u) = W endx(h;µs)W endy(u;µt), x, y = 0,1,2

Note that some models are nested. (The Exp_Exp with M atern_Matern for instance.)

• Spatial bivariate correlation models (see below):

1. Bi_Matern (Bivariate full Matern model)

2. Bi_Matern_contr (Bivariate Matern model with contrainsts)

3. Bi_Matern_sep (Bivariate separable Matern model )

4. Bi_LMC (Bivariate linear model of coregionalization)

5. Bi_LMC_contr (Bivariate linear model of coregionalization with constraints )

6. Bi_W endx (Bivariate full Wendland model)

7. Bi_W endx_contr (Bivariate Wendland model with contrainsts)

8. Bi_W endx_sep (Bivariate separable Wendland model)

• Spatial taper.

For spatial covariance tapering the taper functions are:

GeoCovmatrix 23

1. Bohman deﬁned as:

T(h) = (1 −h)(sin(2πh)/(2πh)) + (1 −cos(2πh))/(2π2h)1[0,1](h)

2. W endlandx, x = 0,1,2deﬁned as:

T(h) = W endx(h;x+ 2), x = 0,1,2

• Spatio-temporal tapers.

For spacetime covariance tapering the taper functions are:

1. W endlandx_W endlandy (Separable tapers) x, y = 0,1,2deﬁned as:

T(h, u) = W endx(h;x+ 2)W endy(h;y+ 2), x, y = 0,1,2.

2. W endlandx_time (Non separable temporal taper) x= 0,1,2deﬁned as: W enx_time,

x= 0,1,2assuming αt= 2,µs= 3.5 + xand γ∈[0,1] to be ﬁxed using tapsep.

3. W endlandx_space (Non separable spatial taper) x= 0,1,2deﬁned as: W enx_space,

x= 0,1,2assuming αs= 2,µt= 3.5 + xand γ∈[0,1] to be ﬁxed using tapsep.

• Spatial bivariate taper (see below).

1. Bi_W endlandx, x = 0,1,2

Remarks:

In what follows we assume σ2, σ2

1, σ2

2, τ2, τ 2

1, τ2

2, a, as, at, a11, a22, a12, κ11, κ22, κ12, f11, f12, f21, f22

positive.

The associated parameters in param are sill,sill_1,sill_2,nugget,nugget_1,nugget_2,scale,scale_s,scale_t,

scale_1,scale_2,scale_12,smooth_1,smooth_2,smooth_12,a_1,a_12,a_21,a_2 respectively.

Let R(h)be a spatial correlation model given in standard notation. Then the covariance model

applied with arbitrary variance, nugget and scale equals to:

C(h)=(σ2+τ21(h= 0))R(h/a, , ...), h ≥0

Similarly if R(h, u)is a spatio-temporal correlation model given in standard notation, then the

covariance model is:

C(h, u)=(σ2+τ21(h= 0, u = 0))R(h/as, u/at, ...)h≥0, u ≥0

Here ‘...’ stands for additional parameters.

Let R(h)be a spatial taper given in standard notation. Then the taper function applied with an

arbitrary compact support (ds) equals to:

T(h) = R(h/ds)

Then the tapered covariance function is given by:

Ctap(h) = T(h)C(h)

Similarly if R(h, u)is a spatio-temporal taper given in standard notation, then the taper function

applied with arbitrary compact supports (ds, dt)Tequals to:

T(h, u) = R(h/ds, u/dt)

Then the tapered covariance function is given by:

Ctap(h, u) = T(h, u)C(h, u)

Compact supports dsand dtcan be set by the user with maxdist and maxtime.

The bivariate models implemented are the following :

24 GeoCovmatrix

1. Bi_Matern deﬁned as:

Cij (h) = ρij (σiσj+τ2

i1(i=j, h = 0))M atern(h/aij , κij )i, j = 1,2. h ≥0

where ρ=ρ12 =ρ21 is the correlation colocated parameter and ρii = 1. The model

Bi_Matern_sep (separable matern) is a special case when a=a11 =a12 =a22 and

κ=κ11 =κ12 =κ22. The model Bi_Matern_contr (constrained matern) is a special case

when a12 = 0.5(a11 +a22)and κ12 = 0.5(κ11 +κ22)

2. Bi_W endx (x= 0,1,2) deﬁned as:

Cij (h) = ρij (σiσj+τ2

i1(i=j, h = 0))W endx(h/aij , νij + 1) i, j = 1,2. h ≥0

where ρ=ρ12 =ρ21 is the correlation colocated parameter and ρii = 1. The model

Bi_W endx_sep (separable wendland) is a special case when a=a11 =a12 =a22 and

µ=µ11 =µ12 =µ22. The model Bi_W endx_contr (constrained matern) is a special case

when a12 = 0.5(a11 +a22)and µ12 = 0.5(µ11 +µ22)

3. Bi_LMC deﬁned as:

Cij (h) =

k=1

(fikfjk +τ2

i1(i=j, h = 0))R(h/ak)

where R(h)is a correlation model. The model Bi_LM C_contr is a special case when f=

f12 =f21. Bivariate LMC models, in the current version of the package, is obtained with

R(h)equal to the exponential correlation model.

The bivariate spatial tapers implemented are the following :

1. Bi_W endlandx,x= 0,1,2deﬁned as:

Tij (h) = rij W endx(h/dij , x), i, j = 1,2x= 0,1,2h≥0

with rii = 1 and r12 =r21 to be ﬁxed using tapsep.

If Tij (h)is a bivariate taper, Then the tapered bivariate covariance function is given by:

Ctap

ij (h) = Tij (h)Cij (h)

Compact supports d11, d12, d22 can be set by the user with maxdist.

Value

Returns an object of class CovMat. An object of class CovMat is a list containing at most the

following components:

bivariate Logical:TRUE if the Gaussian random ﬁeld is bivariaete otherwise FALSE;

coordx Ad-dimensional vector of spatial coordinates;

coordy Ad-dimensional vector of spatial coordinates;

coordt At-dimensional vector of temporal coordinates;

coordx_dyn A list of tmatrices of spatial coordinates;

covmatrix The covariance matrix if type isStandard. An object of class spam if type is

Tapering or Standard and sparse is TRUE.

corrmodel String: the correlation model;

distance String: the type of spatial distance;

GeoCovmatrix 25

grid Logical:TRUE if the spatial data are in a regular grid, otherwise FALSE;

nozero In the case of tapered matrix the percentage of non zero values in the covariance

matrix. Otherwise is NULL.

maxdist Numeric: the marginal spatial compact support if type is Tapering;

maxtime Numeric: the marginal temporal compact support if type is Tapering;

nThe number of trial for Binomial RFs

namescorr String: The names of the correlation parameters;

numcoord Numeric: the number of spatial coordinates;

numtime Numeric: the number the temporal coordinates;

model The type of RF, see GeoFit.

param Numeric: The covariance parameters;

tapmod String: the taper model if type is Tapering. Otherwise is NULL.

spacetime TRUE if spatio-temporal and FALSE if spatial covariance model;

sparse Logical: is the returned object of class spam? ;

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

References

Daley J. D., Porcu E., Bevilacqua M. (2015) Classes of compactly supported covariance functions

for multivariate random ﬁelds. Stochastic Environmental Research and Risk Assessment. 29 (4),

1249–1263.

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

Gneiting, T. (2013), Strictly and Non-Strictly Positive Deﬁnite Functions on Spheres Bernoulli, 19,

1327-1349.

Gneiting, T. (2002). Nonseparable, stationary covariance functions for space-time data. Journal of

the American Statistical Association, 97, 590–600.

Gneiting T, Kleiber W., Schlather M. 2010. Matern cross-covariance functions for multivariate

random ﬁelds. Journal of the American Statistical Association, 105, 1167–1177.

Porcu, E.,Bevilacqua, M. and Genton M. (2015) Spatio-Temporal Covariance and Cross-Covariance

Functions of the Great Circle Distance on a Sphere. Journal of the American Statistical Association.

DOI: 10.1080/01621459.2015.1072541

See Also

GeoKrig,GeoSim,GeoFit

Examples

library(GeoModels)

library(spam)

################################################################

###

### Example 1. Spatial covariance matrix associated to

### a Matern correlation model

26 GeoCovmatrix

###

###############################################################

# Define the spatial-coordinates of the points:

x <- runif(500, 0, 1)

y <- runif(500, 0, 1)

coords <- cbind(x,y)

# Correlation Parameters for Matern model

CorrParam("Matern")

# Matern Parameters

param=list(smooth=0.5,sill=1,scale=0.2,nugget=0)

matrix1 <- GeoCovmatrix(coordx=coords, corrmodel="Matern", param=param)

dim(matrix1$covmatrix)

################################################################

###

### Example 2. Spatial tapered Covariance matrix associated to

### a Matern correlation model

###

###############################################################

matrix2 <- GeoCovmatrix(coordx=coords, corrmodel="Matern", param=param,

maxdist=0.2,taper="Wendland1",type="Tapering")

# Tapered covariance matrix

as.matrix(matrix2$covmatrix)[1:15,1:15]

# Percentage of no zero values in the tapered matrix

matrix2$nozero

################################################################

###

### Example 3. Spatial covariance matrix associated to

### a Generalized Wendland correlation model

###

###############################################################

# Gen Wendland Parameters

param=list(sill=1,scale=0.2,nugget=0,smooth=0,power2=4)

matrix3 <- GeoCovmatrix(coordx=coords, corrmodel="GenWend", param=param,

,sparse=TRUE,maxdist=0.1)

# Percentage of no zero values in the tapered matrix

matrix3$nozero

################################################################

###

### Example 4. Spatial covariance matrix associated to

### a Generalized Cauchy correlation model

GeoCovmatrix 27

###

###############################################################

# Gen Cauchy Parameters

param=list(sill=1,scale=0.2,nugget=0,power1=1,power2=1)

# Correlation Parameters for Gen Cauchy model

CorrParam("GenCauchy")

matrix4 <- GeoCovmatrix(coordx=coords, corrmodel="GenCauchy", param=param)

matrix4$covmatrix[1:4,1:4]

################################################################

###

### Example 5. Covariance matrix associated to

### a space-time double exponential correlation model

###

###############################################################

# Define the temporal-coordinates:

times <- seq(1, 4, 1)

# Define covariance parameters

param=list(scale_s=0.3,scale_t=0.5,sill=1)

# Correlation Parameters for double exp model

CorrParam("Exp_Exp")

# Simulation of a spatial Gaussian random field:

matrix5 <- GeoCovmatrix(coordx=coords, coordt=times, corrmodel="Exp_Exp",

param=param)

dim(matrix5$covmatrix)

################################################################

###

### Example 6. Covariance matrix associated to

### a skew gaussian RF with Exp correlation model

###

###############################################################

param=list(sill=1,scale=0.3/3,nugget=0,skew=4)

# Simulation of a spatial Gaussian random field:

matrix6 <- GeoCovmatrix(coordx=coords, corrmodel="Exp", param=param,

model="SkewGaussian")

# covariance matrix

matrix6$covmatrix[1:10,1:10]

################################################################

###

### Example 7. Covariance matrix associated to

### a Weibull RF with Genwend correlation model

28 GeoCovmatrix

###

###############################################################

param=list(sill=1,scale=0.3,nugget=0,shape=4,mean=0,smooth=1,power2=5)

# Simulation of a spatial Gaussian random field:

matrix7 <- GeoCovmatrix(coordx=coords, corrmodel="GenWend", param=param,

sparse=TRUE,model="Weibull")

# covariance matrix

matrix7$nozero

################################################################

###

### Example 8. Covariance matrix associated to

### a binomial gaussian RF with Wendland correlation model

###

###############################################################

param=list(sill=1,scale=0.2,nugget=0,power2=4)

# Simulation of a spatial Gaussian random field:

matrix8 <- GeoCovmatrix(coordx=coords, corrmodel="Wend0", param=param,n=5,

model="Binomial")

# covariance matrix

matrix8$covmatrix[1:10,1:10]

################################################################

###

### Example 9. Covariance matrix associated to

### a bivariate Matern exponential correlation model

###

###############################################################

set.seed(8)

# Define the spatial-coordinates of the points:

x <- runif(10, -1, 1)

y <- runif(10, -1, 1)

coords <- cbind(x,y)

# Parameters

param=list(scale=0.3,mean_1=0,mean_2=0,sill_1=1,sill_2=1,

nugget_1=0,nugget_2=0,smooth=0.5,pcol=-0.25)

# Covariance matrix

matrix9 <- GeoCovmatrix(coordx=coords, corrmodel="Bi_matern_sep", param=param)$covmatrix

head(matrix8)

GeoFit 29

GeoFit Max-Likelihood-Based Fitting of Gaussian and non Gaussian RFs.

Description

Maximum weighted composite-likelihood ﬁtting for Gaussian and some Non-Gaussian univariate

spatial, spatio-temporal and bivariate spatial RFs The function returns the model parameters’ esti-

mates and the estimates’ variances. Moreover the function allows to ﬁx any of the parameters and

setting upper/lower bound in the optimization.

Usage

GeoFit(data, coordx, coordy=NULL, coordt=NULL,coordx_dyn=NULL, corrmodel, distance="Eucl",

fixed=NULL,GPU=NULL, grid=FALSE, likelihood='Marginal', local=c(1,1),lower=NULL,

maxdist=NULL, maxtime=NULL, method="cholesky", model='Gaussian',

n=1, onlyvar=FALSE ,optimizer='Nelder-Mead',radius=6371, sensitivity=FALSE,

sparse=FALSE, start=NULL, taper=NULL, tapsep=NULL,type='Pairwise',

upper=NULL, varest=FALSE, vartype='SubSamp', weighted=FALSE, winconst=NULL,

winstp=NULL,winconst_t=NULL, winstp_t=NULL,X=NULL)

Arguments

data Ad-dimensional vector (a single spatial realisation) or a (d×d)-matrix (a single

spatial realisation on regular grid) or a (t×d)-matrix (a single spatial-temporal

realisation) or an (d×d×t×n)-array (a single spatial-temporal realisation on

regular grid). For the description see the Section Details.

coordx A numeric (d×2)-matrix (where dis the number of spatial sites) assigning 2-

dimensions of spatial coordinates or a numeric d-dimensional vector assigning

1-dimension of spatial coordinates. Coordinates on a sphere for a ﬁxed radius

radius are passed in lon/lat format expressed in decimal degrees.

coordy A numeric vector assigning 1-dimension of spatial coordinates; coordy is in-

terpreted only if coordx is a numeric vector or grid=TRUE otherwise it will be

ignored. Optional argument, the default is NULL then coordx is expected to be

numeric a (d×2)-matrix.

coordt A numeric vector assigning 1-dimension of temporal coordinates. Optional ar-

gument, the default is NULL then a spatial RF is expected.

coordx_dyn A list of mnumeric (dt×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 De-

tails.

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See the Section Details.

fixed An optional named list giving the values of the parameters that will be consid-

ered as known values. The listed parameters for a given correlation function will

be not estimated.

GPU Numeric; if NULL (the default) no OpenCL computation is performed. The user

can choose the device to be used. Use DeviceInfo() function to see available

devices, only double precision devices are allowed

30 GeoFit

grid Logical; if FALSE (the default) the data are interpreted as spatial or spatial-

temporal realisations on a set of non-equispaced spatial sites (irregular grid).

likelihood String; the conﬁguration of the composite likelihood. Marginal is the default,

see the Section Details.

local Numeric; number of local work-items of the OpenCL setup

lower An optional named list giving the values for the lower bound of the space pa-

rameter when the optimizer is L-BFGS-B or optimize. T he names of the list

must be the same of the names in the start list.

maxdist Numeric; an optional positive value indicating the maximum spatial distance

considered in the composite or tapered likelihood computation. See the Section

Details for more information.

maxtime Numeric; an optional positive value indicating the maximum temporal separa-

tion considered in the composite or tapered likelihood computation (see De-

tails).

method String; the type of matrix decomposition used in the simulation. Default is

cholesky. The other possible choices is svd.

model String; the type of RF and therefore the densities associated to the likelihood

objects. Gaussian is the default, see the Section Details.

nNumeric; number of trials in a binomial RF; number of successes in a negative

binomial RF

optimizer String; the optimization algorithm (see optim for details). Nelder-Mead is the

default. Other possible choices are nlm and L-BFGS-B. In this last case upper

and lower bounds can be passed by the user. In the case of one-dimensional

optimization, the function optimize is used.

onlyvar Logical; if TRUE (and varest is TRUE) only the variance covariance matrix is

computed without optimizing. FALSE is the default.

radius Numeric; the radius of the sphere in the case of lon-lat coordinates. The default

is 6371, the radius of the earth.

sensitivity Logical; if TRUE then the sensitivy matrix is computed

sparse Logical; if TRUE then maximum likelihood is computed using sparse matrices

algorithms (spam packake).It should be used with compactly supported covari-

ance models.FALSE is the default.

start An optional named list with the initial values of the parameters that are used

by the numerical routines in maximization procedure. NULL is the default (see

Details).

taper String; the name of the type of covariance matrix. It can be Standard (the

default value) or Tapering for taperd covariance matrix.

tapsep Numeric; an optional value indicating the separabe parameter in the space time

adaptive taper (see Details).

type String; the type of the likelihood objects. If Pairwise (the default) then the

marginal composite likelihood is formed by pairwise marginal likelihoods (see

Details).

upper An optional named list giving the values for the upper bound of the space pa-

rameter when the optimizer is or L-BFGS-B or optimize. The names of the list

must be the same of the names in the start list.

varest Logical; if TRUE the estimates’ variances and standard errors are returned. FALSE

is the default.

GeoFit 31

vartype String; (SubSamp the default) the type of method used for computing the esti-

mates’ variances, see the Section Details.

weighted Logical; if TRUE the likelihood objects are weighted, see the Section Details. If

FALSE (the default) the composite likelihood is not weighted.

winconst Numeric; a positive value for computing the spatial sub-window in the sub-

sampling procedure. See Details for more information.

winstp Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

spatial sub-sampling procedure. The case 1correspond to no overlapping. See

Details for more information.

winconst_t Numeric; a positive value for computing the temporal sub-window in the sub-

sampling procedure. See Details for more information.

winstp_t Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

temporal sub-sampling procedure. The case 1correspond to no overlapping.

See Details for more information.

XNumeric; Matrix of spatio(temporal)covariates in the linear mean speciﬁcation.

Details

Note, that the standard likelihood may be seen as particular case of the composite likelihood. In

this respect GeoFit provides maximum likelihood ﬁtting. Only composite likelihood estimation

based on pairs are considered. Speciﬁcally marginal pairwise likelihood is considered for each type

of random ﬁeld (Gaussian and not Gaussian). In the Gaussian case other types of estimation are

available: conditional and difference pairwise likelihood, covariance tapering and cross-validation

method.

The optimization method is speciﬁed using optimizer. The default method is Nelder-mead and

other available methods are nlm,BFGS and L-BFGS-B. In this last case upper and lower bounds

constraints in the optimization can be speciﬁed using lower and upper.

Depending on the dimension of data and on the name of the correlation model, the observations

are assumed as a realization of a spatial, spatio-temporal or bivariate RF. Speciﬁcally, with data,

coordx,coordy,coordt parameters:

• If data is a numeric d-dimensional vector, coordx and coordy are two numeric d-dimensional

vectors (or coordx is (d×2)-matrix and coordy=NULL), then the data are interpreted as a single

spatial realisation observed on dspatial sites;

• If data is a numeric (t×d)-matrix, coordx and coordy are two numeric d-dimensional vectors

(or coordx is (d×2)-matrix and coordy=NULL), coordt is a numeric t-dimensional vector,

then the data are interpreted as a single spatial-temporal realisation of a RF observed on d

spatial sites and for ttimes.

• If data is a numeric (2×d)-matrix, coordx and coordy are two numeric d-dimensional vectors

(or coordx is (d×2)-matrix and coordy=NULL), then the data are interpreted as a single spatial

realisation of a bivariate RF observed on dspatial sites.

• If data is a list, coordxdyn is a list and coordt is a numeric t-dimensional vector, then the

data are interpreted as a single spatial-temporal realisation of a RF observed on dynamical

spatial sites (different locations sites for each temporal instants) and for ttimes.

It is possible consider data on a regular grid using option grid=TRUE. Speciﬁcally:

• If data is a numeric (d×d)-matrix, coordx and coordy are two numeric d-dimensional

vectors, grid=TRUE, then the data are interpreted as a single spatial RF realisation observed

on dequispaced spatial sites (named regular grid).

32 GeoFit

• If data is a numeric (d×d×t)-array, coordx and coordy are two numeric d-dimensional

vectors, coordt is a numeric t-dimensional vector, grid=TRUE, then the data are interpreted

as a single spatial-temporal realisation of a RF observed on dequispaced spatial sites and for

ttimes.

• If data is a numeric (d×d×2)-array, coordx and coordy are two numeric d-dimensional

vectors, grid=TRUE, then the data are interpreted as a single realisation of a bivariate RF

observed on dequispaced spatial sites.

When grid=FALSE is is possible to specify a matrix of covariates using X. Speciﬁcally:

• In the spatial case Xmust be a (d×k) covariates matrix associated to data a numeric d-

dimensional vector;

• In the spatiotemporal case Xmust be a (N×k) covariates matrix associated to data a numeric

(t×d)-matrix, where N=t×d;

• In the spatiotemporal case Xmust be a (N×k) covariates matrix associated to data a numeric

(t×d)-matrix, where N= 2 ×d;

The corrmodel parameter allows to select a speciﬁc correlation function for the RF. (See GeoCovmatrix

The distance parameter allows to consider differents kinds of spatial distances. The settings alter-

natives are:

1. Eucl, the euclidean distance (default value);

2. Chor, the chordal distance;

3. Geod, the geodesic distance;

The likelihood parameter represents the composite-likelihood conﬁgurations. The settings alter-

natives are:

1. Conditional, the composite-likelihood is formed by conditionals likelihoods;

2. Marginal, the composite-likelihood is formed by marginals likelihoods;

3. Full, the composite-likelihood turns out to be the standard likelihood;

The model parameter indicates the type of RF considered. The available options are:

RF with marginal symmetric distribution:

•Gaussian, for a Gaussian RF

•Logistic, for a Logistic RF (see Bevilacqua, M. Caamano C., Gaetan , C. 2018)

•StudentT, for a StudentT RF (see Bevilacqua, M. Caamano C., 2018)

RF with positive right skewed marginal distribution:

•Gamma for a Gamma RF (see Bevilacqua, M. Caamano C., 2018)

•Gamma for a Weibull RF (see Bevilacqua, M. Caamano C., 2018)

•LogGaussian for a LogGaussian RF

•LogLogistics for a LogLogistic RF (see Bevilacqua, M. Caamano C., Gaetan , C. 2018)

RF with with possibly asymmetric marginal distribution:

•SkewGaussian for a skew Gaussian RF (see Hao Zhang, H. and El-Shaarawi A. (2010))

•SinhAsinh for Sinh-arcsinh RF (see Bevilacqua, M. Caamano C., 2018)

GeoFit 33

RF with for directional data

•Wrapped for a wrapped Gaussian RF (see Alegria A., Bevilacqua, M., Porcu, E. (2016))

Rf with marginal counts data

•Binomial for a Binomial RF (see Bevilacqua, M. Caamano C., Gaetan C. (2018))

•BinomialNeg for a negative Binomial RF(see Bevilacqua, M. Caamano C., Gaetan C. (2018))

For a given model the associated parameters are given by nuisance and covariance parameters. In

order to obtain the nuisance parameter use NuisParam. In order to obtain the covariance parameter

associated to a given covariance model use CorrParam.

All the nuisance and covariance parameters must be speciﬁed by the user using the start and the

fixed parameter. Speciﬁcally:

The start parameter allows to specify (as starting values for the optimization) the parameters to be

estimated. The fixed parameter allows to ﬁx some of the parameters.

Regression parameters in the linear specﬁcation must be speciﬁed as mean,mean1,..meank (see

NuisParam). In this case a matrix of covariates with suitable dimension can be speciﬁed using the

parameter X. In the case of a single mean then Xshould not be speciﬁed and it is interpreted as a

vector of ones.

The taper parameter, optional in case that type=Tapering, indicates the type of taper correlation

model. (See GeoCovmatrix)

The type parameter represents the type of likelihood used in the composite-likelihood deﬁnition.

The possible alternatives are listed in the following scheme.

If a Gaussian or (any) non Gaussian RF is considered then the possible combination is marginal

pairwise likelihoods (likelihood=Marginal) and type="Pairwise").

If a Gaussian RF is considered (model=Gaussian) then:

• If the composite is formed by marginal likelihoods (likelihood=Marginal):

–Pairwise, the composite-likelihood is deﬁned by the pairwise likelihoods;

–Difference, the composite-likelihood is deﬁned by likelihoods which are obtained as

difference of the pairwise likelihoods.

• If the composite is formed by conditional likelihoods (likelihood=Conditional)

–Pairwise, the composite-likelihood is deﬁned by the pairwise conditional likelihoods.

• If the composite is formed by a full likelihood (likelihood=Full):

–Standard, the objective function is the classical multivariate likelihood;

–Restricted, the objective function is the restricted version of the full likelihood (e.g.

Harville 1977, see References);

–Tapering, the objective function is the tapered 2 (unbiased estimating equation) version

of the full likelihood (e.g. Kaufman et al. 2008, see References);

–Tapering1, the objective function is the tapered 1 (biased estimating equation) version

of the full likelihood (e.g. Kaufman et al. 2008, see References);

–CV, the objective function is the cross validation estimation method ;

The maxdist parameter set the maximum spatial distance below which pairs of sites with inferior

distances are considered in the composite-likelihood. This can be lower of the maximum spatial

distance. Note that this corresponds to use a weighted composite-likelihood with binary weights.

Pairs with distance less than maxdist have weight 1 and are included in the likelihood computation,

instead those with greater distance have weight 0 and then excluded. The default is NULL, in this

case the effective maximum spatial distance between sites is considered.

34 GeoFit

The same arguments of maxdist are valid for maxtime but here the weigthed composite-likelihood

regards the case of spatial-temporal ﬁeld. At the moment is implemented only for Gaussian RFs.

The default is NULL, in this case the effective maximum temporal lag between pairs of observations

is considered.

In the case of tapering likelihood maxdist and maxtime describes the spatial and temporal compact

support of the taper (see GeoCovmatrix). If they are not speciﬁed then the maximum spatial and

temporal distances are considered. In the case of space time adaptive taper the tapsep parameter

allows to specify the spatio temporal compact support (see GeoCovmatrix).

The varest parameter speciﬁes if the standard error estimation of the estimated parameters must

be computed. For Gaussian RF and standard (restricted) likelihood estimation, standard errors are

computed as square root of the diagonal elements of the Fisher Information matrix (asymptotic

covariance matrix of the estimates under increasing domain). For Gaussian RF and tapered and

composite likelihood estimation, standard errors estimate are computed as square root of the diag-

onal elements of the Godambe Information matrix. (asymptotic covariance matrix of the estimates

under increasing domain (see Shaby, B. and D. Ruppert (2012) for tapering and Bevilacqua et.

al. (2012) , Bevilacqua and Gaetan (2013) for weighted composite likelihood)). The vartype pa-

rameter speciﬁes the method used to compute the estimates’ variances in the composite likelihood

case. In particular for estimating the variability matrix Jin the Godambe expression matrix. This

parameter is considered if varest=TRUE. The options are:

•SubSamp (the default), indicates the Sub-Sampling method;

The weighted parameter speciﬁes if the likelihoods forming the composite-likelihood must be

weighted. If TRUE the weights are selected by opportune procedures that improve the efﬁcient

of the maximum composite-likelihood estimator (not implemented yet). If FALSE the efﬁcient im-

provement procedure is not used.

For computing the standard errors by the sub-sampling procedure, winconst and winstp param-

eters represent respectively a positive constant used to determine the sub-window size and the the

step with which the sub-window moves.

In the spatial case (subset of R2), the domain is seen as a rectangle B×H, therefore the size of the

sub-window side bis given by b=winconst ×p(B)(similar is of h). For a complete description

see Lee and Lahiri (2002). By default winconst is set B/(4 ×p(B)). This way we ensure to

have at least eight spatial blocks. The winstp parameter is used to determine the sub-window step.

The latter is given by the proportion of the sub-window size, so that when winstp=1 there is not

overlapping between contiguous sub-windows. In the spatial case by default winstp=0.5. The

sub-window is moved by successive steps in order to cover the entire spatial domain. Observations,

that fall in disjoint or overlapping windows are considered indipendent samples.

In the spatio-temporal case the subsampling is meant only in time as described by Li et al. (2007).

Thus, winconst_t represents the lenght of the temporal sub-window. By default the size of the

sub-window is computed following the rule established in Li et al. (2007). By default winstp is the

time step.

Value

Returns an object of class GeoFit. An object of class GeoFit is a list containing at most the

following components:

bivariate Logical:TRUE if the Gaussian RF is bivariate, otherwise FALSE;

clic The composite information criterion, if the full likelihood is considered then it

coincides with the Akaike information criterion;

coordx Ad-dimensional vector of spatial coordinates;

GeoFit 35

coordy Ad-dimensional vector of spatial coordinates;

coordt At-dimensional vector of temporal coordinates;

coordx_dyn A list of dynamical (in time) spatial coordinates;

convergence A string that denotes if convergence is reached;

corrmodel The correlation model;

data The vector or matrix or array of data;

distance The type of spatial distance;

fixed The vector of ﬁxed parameters;

iterations The number of iteration used by the numerical routine;

likelihood The conﬁguration of the composite likelihood;

logCompLik The value of the log composite-likelihood at the maximum;

maxdist The maximum spatial distance used in the weigthed composite likelihood. If no

spatial distance is speciﬁed then it is NULL;

maxtime The maximum temporal distance used in the weigthed composite likelihood. If

no spatial distance is speciﬁed then it is NULL;

message Extra message passed from the numerical routines;

model The density associated to the likelihood objects;

nThe number of trials in a binominal RF;the number of successes in a negative

Binomial RFs;

nozero In the case of tapered likelihood the percentage of non zero values in the covari-

ance matrix. Otherwise is NULL.

numcoord The number of spatial coordinates;

numtime The number of the temporal realisations of the RF;

param The vector of parameters’ estimates;

param The radius of the sphere in the case of great circle distance;

stderr The vector of standard errors;

sensmat The sensitivity matrix;

varcov The matrix of the variance-covariance of the estimates;

varimat The variability matrix;

vartype The method used to compute the variance of the estimates;

type The type of the likelihood objects.

winconst The constant used to compute the window size in the spatial sub-sampling;

winstp The step used for moving the window in the spatial sub-sampling;

winconst_t The constant used to compute the window size in the spatio-temporal sub-sampling;

winstp_ The step used for moving the window in the spatio-temporal sub-sampling;

XThe matrix of covariates;

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

36 GeoFit

References

Maximum Restricted Likelihood Estimator:

Harville, D. A. (1977) Maximum Likelihood Approaches to Variance Component Estimation and

to Related Problems. Journal of the American Statistical Association,72, 320–338.

Tapered likelihood:

Kaufman, C. G., Schervish, M. J. and Nychka, D. W. (2008) Covariance Tapering for Likelihood-

Based Estimation in Large Spatial Dataset. Journal of the American Statistical Association,103,

1545–1555.

Shaby, B. and D. Ruppert (2012). Tapered covariance: Bayesian estimation and asymptotics. J.

Comp. Graph. Stat.,21-2, 433–452.

Composite-likelihood:

Varin, C., Reid, N. and Firth, D. (2011). An Overview of Composite Likelihood Methods. Statistica

Sinica,21, 5–42.

Varin, C. and Vidoni, P. (2005) A Note on Composite Likelihood Inference and Model Selection.

Biometrika,92, 519–528.

Weighted Composite-likelihood for binary RFs:

Patrick, J. H. and Subhash, R. L. (1998) A Composite Likelihood Approach to Binary Spatial Data.

Journal of the American Statistical Association, Theory & Methods,93, 1099–1111.

Weighted Composite-likelihood for wrapped directional RFs:

Alegria A., Bevilacqua, M., Porcu, E. (2016) Likelihood-based inference for multivariate space-

time wrapped-Gaussian ﬁelds. Journal of Statistical Computation and Simulation.86(13), 2583–

2597.

Weighted Composite-likelihood for Gaussian RFs:

Bevilacqua, M. Gaetan, C., Mateu, J. and Porcu, E. (2012) Estimating space and space-time co-

variance functions for large data sets: a weighted composite likelihood approach. Journal of the

American Statistical Association, Theory & Methods,107, 268–280.

Bevilacqua, M. Gaetan, C. (2013) On composite likelihood inference based on pairs for spatial

Gaussian RFs Techical Report, Department of Statistics, de Valparaiso University.

Sub-sampling estimation:

Carlstein, E. (1986) The Use of Subseries Values for Estimating the Variance. The Annals of Statis-

tics,14, 1171–1179.

Heagerty, P. J. and Lumley T. (2000) Window Subsampling of Estimating Functions with Applica-

tion to Regression Models. Journal of the American Statistical Association, Theory & Methods,95,

197–211.

Lee, Y. D. and Lahiri S. N. (2002) Variogram Fitting by Spatial Subsampling. Journal of the Royal

Statistical Society. Series B,64, 837–854.

Li, B., Genton, M. G. and Sherman, M. (2007). A nonparametric assessment of properties of space-

time covariance functions. Journal of the American Statistical Association,102, 736–744

Examples

library(GeoModels)

library(fields)

###############################################################

############ Examples of spatial Gaussian RFs ################

GeoFit 37

###############################################################

################################################################

###

### Example 1, 2, 3: Estimation of a spatial Gaussian RF with

### exponential correlation using and pairwise likelihood

### maximum likelihood and tapering likelihood

###############################################################

# Define the spatial-coordinates of the points:

set.seed(3)

N=400 # number of location sites

x <- runif(N, 0, 1)

set.seed(6)

y <- runif(N, 0, 1)

coords <- cbind(x,y)

# Define spatial matrix covariates

X=cbind(rep(1,N),runif(N))

# Set the covariance model's parameters:

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)

# Simulation of the spatial Gaussian RF:

data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param,X=X)$data

fixed<-list(nugget=nugget)

start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill)

################################################################

###

### Example 1. Maximum pairwise likelihood fitting of

### Gaussian RFs with exponential correlation.

###

###############################################################

fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,

maxdist=0.05,likelihood="Marginal",type="Pairwise",

start=start,fixed=fixed,X=X)

print(fit1)

################################################################

###

### Example 2. Standard Maximum likelihood fitting of

### Gaussian RFs with exponential correlation.

###

###############################################################

fit2 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,

likelihood="Full",type="Standard",

start=start,fixed=fixed,X=X)

print(fit2)

38 GeoFit

################################################################

###

### Example 3. Tapered Maximum likelihood fitting of

### Gaussian RFs with exponential correlation.

###

###############################################################

fit3 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,

likelihood="Full",type="Tapering",taper="Wendland1",

maxdist=0.1,start=start,fixed=fixed,X=X)

print(fit3)

###############################################################

############ Examples of spatial non-Gaussian RFs #############

###############################################################

################################################################

###

### Example 4. Maximum pairwise likelihood fitting of spatial

### Gamma and Weibull RFs with Wendland correlation

###

###############################################################

set.seed(524)

# Define the spatial-coordinates of the points:

N=500

x <- runif(N, 0, 1)

y <- runif(N, 0, 1)

coords <- cbind(x,y)

X=cbind(rep(1,N),runif(N))

mean=1; mean1=2 # regression parameters

nugget=0

shape=2

scale=0.2

model="Weibull"

corrmodel="Wend0"

param=list(mean=mean,mean1=mean1,sill=1-nugget,scale=scale,shape=shape,nugget=nugget,power2=4)

# Simulation of a non stationary weibull RF:

data <- GeoSim(coordx=coords, corrmodel=corrmodel,model=model,X=X,

param=param)$data

fixed<-list(nugget=nugget,power2=4,sill=1-nugget)

start<-list(mean=mean,mean1=mean1,scale=scale,shape=shape)

# Maximum pairwise composite-likelihood fitting of the RF:

fit <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,

maxdist=.05,likelihood="Marginal",type="Pairwise",X=X,

start=start,fixed=fixed)

print(fit$param)

model="Gamma"

start<-list(mean=mean,mean1=mean1,scale=scale)

fixed<-list(nugget=nugget,power2=4,sill=1-nugget,shape=6)

GeoFit 39

# Maximum pairwise composite-likelihood fitting of the RF:

fit <- GeoFit(data=data,coordx=coords,corrmodel="Wend0", model=model,

maxdist=.05,likelihood="Marginal",type="Pairwise",X=X,

start=start,fixed=fixed)

print(fit$param)

################################################################

###

### Example 5. Maximum pairwise likelihood fitting of

### StudendT spatial RFs with Wendland correlation

###

###############################################################

set.seed(15274)

# Define the spatial-coordinates of the points:

N=300

x <- runif(N, 0, 1)

y <- runif(N, 0, 1)

coords <- cbind(x,y)

X=cbind(rep(1,N),runif(N))

mean=1; mean1=2 # regression parameters

nugget=0

sill=0.5

scale=0.2

df=4 ## degrees of freedom

model="StudentT"

corrmodel="Wend0"

param=list(mean=mean,mean1=mean1,sill=sill,scale=scale,df=1/df,nugget=nugget,power2=4)

# Simulation of a studentT RF:

data <- GeoSim(coordx=coords, corrmodel=corrmodel,model=model,X=X,

param=param)$data

## estimation assuming df unknown

fixed<-list(nugget=nugget,power2=4)

start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill,df=1/df)

# Maximum pairwise composite-likelihood fitting of the RF:

fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,

maxdist=.05,likelihood="Marginal",type="Pairwise",X=X,

start=start,fixed=fixed)

print(fit1$param)

## df must be rounded and fixed

df=round(1/(as.numeric(fit1$param['df'])))

fixed<-list(nugget=nugget,power2=4,df=1/df)

start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill)

# Maximum pairwise composite-likelihood fitting of the RF:

fit <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,

maxdist=.05,likelihood="Marginal",type="Pairwise",X=X,

start=start,fixed=fixed)

print(fit$param)

################################################################

###

### Example 6. Maximum pairwise likelihood fitting of

### SinhAsinh-Gaussian spatial RFs with Wendland correlation

###

###############################################################

40 GeoFit

set.seed(261)

model="SinhAsinh"

# Define the spatial-coordinates of the points:

x <- runif(500, 0, 1)

y <- runif(500, 0, 1)

coords <- cbind(x,y)

corrmodel="Wend0"

mean=0;nugget=0

sill=1

skew=-0.5

tail=1.5

power2=4

c_supp=0.2

# model parameters

param=list(power2=power2,skew=skew,tail=tail,

mean=mean,sill=sill,scale=c_supp,nugget=nugget)

data <- GeoSim(coordx=coords, corrmodel=corrmodel,model=model, param=param)$data

plot(density(data))

fixed=list(power2=power2,nugget=nugget)

start=list(scale=c_supp,skew=skew,tail=tail,mean=mean,sill=sill)

# Maximum pairwise likelihood:

fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,

maxdist=0.05,likelihood="Marginal",type="Pairwise",

start=start,fixed=fixed)

print(fit1$param)

# Maximum likelihood:

fit2 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,

likelihood="Full",type="Standard",

start=start,fixed=fixed)

print(fit2$param)

################################################################

###

### Example 7. Maximum pairwise likelihood fitting of

### Binomial and negative Binomial RFs

### with exponential correlation.

###

###############################################################

set.seed(422)

N=350

x <- runif(N, 0, 1)

y <- runif(N, 0, 1)

coords <- cbind(x,y)

mean=0.1; mean1=0.8; mean2=-0.5 # regression parameters

X=cbind(rep(1,N),runif(N),runif(N)) # marix covariates

corrmodel <- "Wend0"

param=list(mean=mean,mean1=mean1,mean2=mean2,sill=1,nugget=0,scale=0.2,power2=4)

# Simulation of the spatial Binomial-Gaussian RF:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, model="Binomial", n=10,X=X,

param=param)$data

fixed <- list(nugget=nugget,power2=4,sill=1)

GeoFit 41

start <- list(scale=0.2,mean=mean,mean1=mean1,mean2=mean2)

# Maximum pairwise likelihood:

fit1 <- GeoFit(data=data, coordx=coords, corrmodel=corrmodel,n=10, X=X,

maxdist=0.05,model="Binomial", fixed=fixed,

start=start)

print(fit1)

set.seed(220)

# Simulation of the spatial Negative Binomial-Gaussian RF:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, model="BinomialNeg", n=5,X=X,

param=param)$data

# Maximum pairwise likelihood:

fit2 <- GeoFit(data=data, coordx=coords, corrmodel=corrmodel, n=5,X=X,

maxdist=0.05,model="BinomialNeg", fixed=fixed,

start=start)

print(fit2)

###############################################################

######### Examples of spatio-temporal RFs ###########

###############################################################

set.seed(52)

# Define the temporal sequence:

time <- seq(1, 10, 1)

# Define the spatial-coordinates of the points:

x <- runif(20, 0, 1)

set.seed(42)

y <- runif(20, 0, 1)

coords=cbind(x,y)

# Set the covariance model's parameters:

corrmodel="Exp_Exp"

scale_s=0.2/3

scale_t=1

sill=1

nugget=0

mean=0

param<-list(mean=0,scale_s=scale_s,scale_t=scale_t,

sill=sill,nugget=nugget)

# Simulation of the spatial-temporal Gaussian RF:

data <- GeoSim(coordx=coords,coordt=time,corrmodel=corrmodel,

param=param)$data

################################################################

###

### Example 8. Maximum pairwise likelihood fitting of a

### space time Gaussian RF with double-exponential correlation.

###

###############################################################

# Fixed parameters

fixed<-list(nugget=nugget)

# Starting value for the estimated parameters

start<-list(mean=mean,scale_s=scale_s,scale_t=scale_t,sill=sill)

42 GeoFit

# Maximum composite-likelihood fitting of the RF:

fit <- GeoFit(data=data,coordx=coords,coordt=time,

corrmodel="Exp_Exp",maxtime=1,maxdist=0.3,

likelihood="Marginal",type="Pairwise",

start=start,fixed=fixed)

print(fit)

################################################################

###

### Example 9. Maximum standard likelihood fitting of a

### space time Gaussian RF observed on dynamical spatial coordinates

### with double-exponential correlation.

###############################################################

maxN=50

coordx_dyn=list()

set.seed(31)

for(k in 1:length(time))

{

NN=sample(1:maxN,size=1)

x <- runif(NN, 0, 1)

y <- runif(NN, 0, 1)

coordx_dyn[[k]]=cbind(x,y)

}

data <- GeoSim(coordx_dyn=coordx_dyn, coordt=time, corrmodel="Exp_Exp",

param=param)$data

fit <- GeoFit(data=data,coordx_dyn=coordx_dyn,coordt=time,

corrmodel="Exp_Exp",

likelihood="Full",type="Standard",

start=start,fixed=fixed)

print(fit)

###############################################################

######### Examples of spatial bivariate RFs ###########

###############################################################

################################################################

###

### Example 10. Maximum, and pairwise likelihood fitting of a

### bivariate Gaussian RF with separable Bivariate matern

### (cross) correlation model.

###

###############################################################

# Define the spatial-coordinates of the points:

set.seed(5)

x <- runif(250, 0, 1)

y <- runif(250, 0, 1)

coords=cbind(x,y)

# parameters

param=list(mean_1=0,mean_2=0,scale=0.1,smooth=0.5,sill_1=1,sill_2=1,

GeoKrig 43

nugget_1=0,nugget_2=0,pcol=0.2)

# Simulation of a spatial Gaussian RF:

data <- GeoSim(coordx=coords, corrmodel="Bi_Matern_sep",

param=param)$data

# selecting fixed and estimated parameters

fixed=list(nugget_1=0,nugget_2=0,smooth=0.5)

start=list(mean_1=0,mean_2=0,sill_1=var(data[1,]),sill_2=var(data[2,]),

scale=0.1,pcol=cor(data[1,],data[2,]))

# Maximum pairwise likelihood

fitcl<- GeoFit(data=data, coordx=coords, corrmodel="Bi_Matern_sep",

likelihood="Marginal",type="Pairwise",start=start,

fixed=fixed,maxdist=c(0.05,0.05,0.05))

print(fitcl)

# Maximum likelihood :

fitml<- GeoFit(data=data, coordx=coords, likelihood="Full",

corrmodel="Bi_Matern_sep", type="Standard",

start=start, fixed=fixed)

print(fitml)

GeoKrig Spatial (bivariate) and spatio temporal optimal linear prediction for

Gaussian and non Gaussian RFs.

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(data, coordx, coordy=NULL, coordt=NULL, coordx_dyn=NULL, corrmodel,distance="Eucl",

grid=FALSE, loc, maxdist=NULL, maxtime=NULL, method="cholesky",

model="Gaussian", n=1,nloc=NULL,mse=FALSE, lin_opt=TRUE,

param, radius=6371, sparse=FALSE,taper=NULL,tapsep=NULL,

time=NULL, type="Standard",type_mse=NULL,

type_krig="Simple",weigthed=TRUE,which=1, X=NULL,Xloc=NULL)

Arguments

data Ad-dimensional vector (a single spatial realisation) or a (d×d)-matrix (a single

spatial realisation on regular grid) or a (t×d)-matrix (a single spatial-temporal

realisation) or an (d×d×t×n)-array (a single spatial-temporal realisation on

regular grid) giving the data used for prediction.

coordx A numeric (d×2)-matrix (where dis the number of spatial sites) giving 2-

dimensions of spatial coordinates or a numeric d-dimensional vector giving 1-

dimension of spatial coordinates used for prediction. qndd-dimensional vector

giving 1-dimension of spatial coordinates. Coordinates on a sphere for a ﬁxed

radius radius are passed in lon/lat format expressed in decimal degrees.

44 GeoKrig

coordy A numeric vector giving 1-dimension of spatial coordinates used for prediction;

coordy is interpreted only if coordx is a numeric vector or grid=TRUE other-

wise it will be ignored. Optional argument, the default is NULL then coordx is

expected to be numeric a (d×2)-matrix.

coordt A numeric vector giving 1-dimension of temporal coordinates used for predic-

tion; the default is NULL then a spatial random ﬁeld is expected.

coordx_dyn A list of mnumeric (dt×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 De-

tails.

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).

lin_opt Logical;If TRUE (default) then optimal (pairwise) linear kriging is computed.

Otherwise optimal (pairwise) kriging is computed in the mean square sense.

loc A numeric (n×2)-matrix (where nis the number of spatial sites) giving 2-

dimensions of spatial coordinates to be predicted.

maxdist Numeric; an optional positive value indicating the maximum spatial compact

support in the case of covariance tapering kriging.

maxtime Numeric; an optional positive value indicating the maximum temporal compact

support in the case of covasriance tapering kriging.

method String; the type of matrix decomposition used in the simulation. Default is

cholesky. The other possible choices is svd.

nNumeric; the number of trials in a binomial random ﬁelds. Default is 1.

nloc Numeric; the number of trials of the locations sites to be predicted in a binomial

random ﬁelds type II. 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.

radius Numeric: the radius of the sphere if coordinates are passed in lon/lat format;

sparse Logical; if TRUE kriging is computed with sparse matrices algorithms using spam

package. Default is FALSE. It should be used with compactly supported covari-

ances.

taper String; the name of the taper correlation function, see the Section Details.

tapsep Numeric; an optional value indicating the separabe parameter in the space time

quasi taper (see Details).

time A numeric (m×1) vector (where mis the number of temporal instants) giving the

temporal instants to be predicted; the default is NULL then only spatial prediction

is performed.

type String; if Standard then standard kriging is performed;if Tapering then krig-

ing with covariance tapering is performed;if Pairwise then pairwise kriging is

performed

GeoKrig 45

type_mse String; if Theoretical then theoretical MSE pairwise kriging is computed. If

SubSamp then an estimation based on subsampling is computed.

type_krig String; the type of kriging. If Simple (the default) then simple kriging is per-

formed. (See the Section Details).

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 (tapered) cokriging it indicates which variable

to predict. It can be 1 or 2

XNumeric; Matrix of spatio(temporal)covariates in the linear mean speciﬁcation.

Xloc Numeric; Matrix of spatio(temporal)covariates in the linear mean speciﬁcation

associated to predicted locations.

Details

Best linear unbiased predictor and associated mean square error is computed for Gaussian and some

non Gaussian cases. Speciﬁcally, 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 ﬁxed

parameters and given the type of RF (model) the function computes simple kriging, for the speciﬁed

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 func-

tion. The list param speciﬁes 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 speciﬁed by the parameter mean in the list param (See examples). In the

Gaussian case, it is possible to perform kriging based on covariance tapering for simple kriging

(Furrer et. al, 2008). In this case, space or space-time tapered function and spatial or spatio- tem-

poral compact support must be speciﬁed. For the choice of a space or space-time tapered function

see GeoCovmatrix. When performing kriging with covariance tapering, sparse matrix algorithms

are exploited using the package spam.

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 Ad-dimensional vector of spatial coordinates used for prediction;

coordy Ad-dimensional vector of spatial coordinates used for prediction;

coordt At-dimensional vector of temporal coordinates used for prediction;

corrmodel String: the correlation model;

covmatrix The covariance matrix if type is Standard. An object of class spam if type is

Tapering

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, other-

wise FALSE;

loc A (n×2)-matrix of spatial locations to be predicted.

nThe number of trial for Binomial RFs

46 GeoKrig

nozero In the case of tapered simple kriging the percentage of non zero values in the

covariance matrix. Otherwise is NULL.

numcoord Numeric:he number dof spatial coordinates used for prediction;

numloc Numeric: the number nof spatial coordinates to be predicted;

numtime Numeric: the number dof the temporal instants used for prediction;

numt Numeric: the number mof the temporal instants to be predicted;

model The type of RF, see GeoFit.

param Numeric: The covariance parameters;

pred A (m×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;

tapmod String: the taper model if type is Tapering. Otherwise is NULL.

time Am-dimensional vector of temporal coordinates to be predicted;

type String: the type of kriging (Standard or Tapering).

type_krig String: the type of kriging.

mse A (m×n)-matrix of spatio or spatio temporal mean square error kriging predic-

tion;

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

References

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

Furrer R., Genton, M.G. and Nychka D. (2006). Covariance Tapering for Interpolation of Large

Spatial Datasets. Journal of Computational and Graphical Statistics, 15-3, 502–523.

See Also

GeoCovmatrix

Examples

library(GeoModels)

library(fields)

################################################################

############### Examples of Spatial kriging ###################

################################################################

# Define the spatial-coordinates of the points:

set.seed(79)

x <- runif(200, 0, 1)

y <- runif(200, 0, 1)

coords<-cbind(x,y)

# Set the exponential cov parameters:

corrmodel_1 <- "exponential"

GeoKrig 47

mean<-0

sill<-1

nugget<-0

scale<-0.3/3

param<-list(mean=mean,sill=sill,nugget=nugget,scale=scale)

# Set the wendland parameters (two compatible correlations):

corrmodel_2 <- "Wend0"

mean<-0

sill<-1

nugget<-0

power2=3

c_supp<-0.3

param_wen<-list(mean=mean,sill=sill,nugget=nugget,scale=c_supp,power2=power2)

# Simulation of the spatial Gaussian random field:

data <- GeoSim(coordx=coords, corrmodel=corrmodel_1,

param=param)$data

# locations to predict

xx<-seq(0,1,0.03)

loc_to_pred<-as.matrix(expand.grid(xx,xx))

################################################################

###

### Example 1. Spatial simple kriging of n sites of a

### Gaussian random fields with exponential correlation.

###

###############################################################

pr<-GeoKrig(loc=loc_to_pred,coordx=coords,corrmodel=corrmodel_1,

param= param, data=data,mse=TRUE)

################################################################

###

### Example 2. Spatial tapered simple kriging of n sites of a

### Gaussian random fields with exponential correlation.

###

###############################################################

pr_tap=GeoKrig(loc=loc_to_pred,coordx=coords,corrmodel=corrmodel_1,data=data,

param= param,type="Tapering",maxdist=0.2,taper="Wendland1",mse=TRUE)

################################################################

###

### Example 3. Spatial simple kriging of n sites of a

### Gaussian random fields using a compatible Wendland model

###

###############################################################

param<-list(mean=mean,sill=sill,nugget=nugget,scale=scale,power2=power2)

pr_wen=GeoKrig(loc=loc_to_pred,coordx=coords,corrmodel=corrmodel_2,data=data,

param=param_wen,sparse=TRUE,mse=TRUE)

colour <- rainbow(100)

par(mfrow=c(3,2))

zlim=c(-2.6,2.6)

48 GeoKrig

# simple kriging map prediction

image.plot(xx, xx, matrix(pr$pred,ncol=length(xx)),col=colour,

zlim=zlim,xlab="",ylab="",

main="Simple Kriging with exponential model ")

# simple kriging map prediction variance

image.plot(xx, xx, matrix(pr$mse,ncol=length(xx)),col=colour,

xlab="",ylab="",main="Std error")

# simple tapered kriging map prediction

image.plot(xx, xx, matrix(pr_tap$pred,ncol=length(xx)),col=colour,

zlim=zlim,xlab="",ylab="",main="Simple Tapered Kriging")

# simple taperd kriging map prediction variance

image.plot(xx, xx, matrix(pr_tap$mse,ncol=length(xx)),col=colour,

xlab="",ylab="",main="Std error")

# simple tapered kriging map prediction

image.plot(xx, xx, matrix(pr_wen$pred,ncol=length(xx)),col=colour,

zlim=zlim,xlab="",ylab="",main="Simple Kriging with Wendland model")

# simple kriging map prediction variance

image.plot(xx, xx, matrix(pr_wen$mse,ncol=length(xx)),col=colour,

xlab="",ylab="",main="Std error")

################################################################

###

### Example 4. Spatial simple kriging of a binomial

### random field

###

###############################################################

set.seed(312)

model="Binomial";n=6

# Define the spatial-coordinates of the points:

x <- runif(800)

y <- runif(800)

coords=cbind(x,y)

#### mean and covariance parameters ###

mean<-0

sill<-1

nugget<-0

scale<-0.2

#################

param<-list(mean=mean,sill=sill,nugget=nugget,scale=scale,power2=4)

# Simulation of the Binomial Gaussian random field:

data <- GeoSim(coordx=coords, corrmodel="Wend0",model=model,n=n,

sparse=TRUE,param=param)$data

par(mfrow=c(1,2))

#### map of simulated data

quilt.plot(x, y, data,nlevel=n+1,col=rainbow(n+1),zlim=c(0,n), main="Data")

GeoKrig 49

## estimation with pairwise likelihood

fixed=list(nugget=nugget,power2=4,sill=1)

start=list(mean=0,scale=scale)

# Maximum pairwise likelihood fitting :

fit <- GeoFit(data, coordx=coords, corrmodel="Wend0",model=model,n=n,

likelihood='Marginal', type='Pairwise',maxdist=0.03,

start=start,fixed=fixed)

# locations to predict

xx<-seq(0,1,0.03)

loc_to_pred<-as.matrix(expand.grid(xx,xx))

## simple kriging

pr<-GeoKrig(data=data, coordx=coords,loc=loc_to_pred,corrmodel="Wend0",model=model,n=n,

sparse=TRUE,param= as.list(c(fit$param,fixed)))

#standard binomial kriging

map_binom=matrix(round(pr$pred),ncol=length(xx))

image.plot(xx, xx, map_binom,nlevel=n+1,col=rainbow(n+1),zlim=c(0,n),

xlab="",ylab="",main="Simple Kriging ")

################################################################

###

### Example 5. Spatial simple kriging of a Weibull

### random field

###

###############################################################

set.seed(312)

model="Weibull";

# Define the spatial-coordinates of the points:

x <- runif(800)

y <- runif(800)

coords=cbind(x,y)

#### mean and covariance parameters ###

mean<-0

sill<-1

nugget<-0

scale<-0.2

shape=0.8

#################

param<-list(mean=mean,sill=sill,shape=shape,nugget=nugget,scale=scale,power2=4)

# Simulation of the Weibull random field:

data <- GeoSim(coordx=coords, corrmodel="Wend0",model=model,

sparse=TRUE,param=param)$data

par(mfrow=c(1,2))

#### map of simulated data

quilt.plot(x, y, data, main="Data")

50 GeoKrig

## estimation with pairwise likelihood

fixed=list(nugget=nugget,power2=4,sill=1)

start=list(mean=0,scale=scale,shape=shape)

# Maximum pairwise likelihood fitting :

fit <- GeoFit(data, coordx=coords, corrmodel="Wend0",model=model,

likelihood='Marginal', type='Pairwise',maxdist=0.03,

start=start,fixed=fixed)

# locations to predict

xx<-seq(0,1,0.04)

loc_to_pred<-as.matrix(expand.grid(xx,xx))

#optimal linear kriging

pr<-GeoKrig(data=data, coordx=coords,loc=loc_to_pred,corrmodel="Wend0",model=model,

sparse=TRUE,param= as.list(c(fit$param,fit$fixed)))

map_weibull=matrix(pr$pred,ncol=length(xx))

image.plot(xx, xx, map_weibull,

xlab="",ylab="",main="Simple Kriging ")

################################################################

###

### Example 5. Spatial simple kriging of a t

### random field

###

###############################################################

model="StudentT"

df=6

corrmodel <- "Wend0"

nsel=800

coords=cbind(runif(nsel),runif(nsel))

mean <- 0

sill <- 1.9

nugget <- 0

power2=4

scale <- 0.2

# Starting value for the estimated parameters

set.seed(3132)

param=list(nugget=nugget,mean=mean, scale=scale,sill=sill,df=1/df,power2=power2)

data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param,

model=model,sparse=TRUE)$data

fixed<-list(nugget=nugget,power2=4,df=1/df)

start<-list(mean=mean, scale=scale,sill=sill)

# Maximum pairwise composite-likelihood fitting of the RF:

fit <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,

maxdist=0.02,likelihood="Marginal",type="Pairwise",

start=start,fixed=fixed, model = model)

GeoKrig 51

# locations to predict

xx<-seq(0,1,0.04)

loc_to_pred<-as.matrix(expand.grid(xx,xx))

#optimal linear kriging

pr<-GeoKrig(data=data, coordx=coords,loc=loc_to_pred,corrmodel="Wend0",model=model,

sparse=TRUE,param= as.list(c(fit$param,fit$fixed)))

par(mfrow=c(1,2))

#### map of simulated data

quilt.plot(coords[,1], coords[,2], data, main="Data")

map_t=matrix(pr$pred,ncol=length(xx))

image.plot(xx, xx, map_t,

xlab="",ylab="",main="Simple Kriging ")

################################################################

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

################################################################

# Define the spatial-coordinates of the points:

x <- runif(80, 0, 1)

y <- runif(80, 0, 1)

coords<-cbind(x,y)

times<-1:8

# Define model correlation and associated parameters

corrmodel<-"exp_exp"

param<-list(nugget=0,mean=0,scale_s=0.2/3,scale_t=1,sill=1)

# Simulation of the space time Gaussian random field:

set.seed(31)

data<-GeoSim(coordx=coords,coordt=times,corrmodel=corrmodel,

param=param)$data

# Maximum pairwise likelihood fitting of the space time random field:

start <- list(scale_s=0.2/3,scale_t=1,sill=1,mean=0)

fixed <- list(nugget=0)

fit <- GeoFit(data, coordx=coords, coordt=times,

corrmodel=corrmodel, likelihood='Marginal',

type='Pairwise',start=start,fixed=fixed,

maxdist=0.1,maxtime=1)

################################################################

###

### Example 6. Spatio temporal simple kriging of n locations

### sites and m temporal instants for a Gaussian random fields

### with estimated double exponential correlation.

###

###############################################################

set.seed(3)

param<-as.list(c(fit$param,fit$fixed))

# locations to predict

xx<-seq(0,1,0.04)

loc_to_pred<-as.matrix(expand.grid(xx,xx))

# Define the times to predict

52 GeoKrig

times_to_pred<-1:2

pr<-GeoKrig(loc=loc_to_pred,time=times_to_pred,coordx=coords,coordt=times,

corrmodel=corrmodel, param=param,data=data,mse=TRUE)

par(mfrow=c(2,2))

zlim=c(-2.5,2.5)

colour <- rainbow(100)

for(i in 1:2) {

image.plot(xx, xx, matrix(pr$pred[i,],ncol=length(xx)),col=colour,

zlim=zlim, main = paste(" Kriging Time=" , i),ylab="")

image.plot(xx, xx, matrix(pr$mse[i,],ncol=length(xx)),col=colour,

main = paste("Std err Time=" , i),ylab="")

}

################################################################

########### Examples of spatial bivariate cokriging ############

################################################################

###

### Example 7. Bivariate simple cokriging of n locations

### for a Gaussian random fields with separable Matern correlation

###

###############################################################

# Define the spatial-coordinates of the points:

x <- runif(80, 0, 1)

y <- runif(80, 0, 1)

coords<-cbind(x,y)

# Simulation of a spatial bivariate Gaussian random field

# with Matern separable covariance model

set.seed(12)

param=list(scale=0.3/3,mean_1=0,mean_2=0,sill_1=1,sill_2=1,

nugget_1=0,nugget_2=0,pcol=0.7,smooth=0.5)

data <- GeoSim(coordx=coords, corrmodel="Bi_matern_sep", param=param)$data

fixed=list(nugget_1=0,nugget_2=0,smooth=0.5,mean_1=0,mean_2=0)

start=list(sill_1=var(data[,1]),sill_2=var(data[,2]),scale=0.3/3,

pcol=cor(data[1,],data[2,]))

# Maximum Composite likelihood fitting of the random field:

fitcl<- GeoFit(data, coordx=coords, corrmodel="Bi_matern_sep",

likelihood="Marginal",type="Pairwise",maxdist=0.1,

start=start,fixed=fixed)

# locations to predict

xx<-seq(0,1,0.04)

loc_to_pred<-as.matrix(expand.grid(xx,xx))

pr1<-GeoKrig(loc=loc_to_pred,coordx=coords,corrmodel="Bi_matern_sep",

param= as.list(c(fitcl$param,fitcl$fixed)), data=data,which=1,mse=TRUE)

GeoResiduals 53

pr2<-GeoKrig(loc=loc_to_pred,coordx=coords,corrmodel="Bi_matern_sep",

param= as.list(c(fitcl$param,fitcl$fixed)), data=data,which=2,mse=TRUE)

par(mfrow=c(2,2))

# simple kriging map prediction of the first variable

image.plot(xx, xx, matrix(pr1$pred,ncol=length(xx)),col=colour,

xlab="",ylab="",main="First Simple coKriging")

# simple kriging map prediction variance of the first variable

image.plot(xx, xx, matrix(pr1$mse,ncol=length(xx)),col=colour,

xlab="",ylab="",main="Std error")

# simple kriging map prediction of the second variable

image.plot(xx, xx, matrix(pr2$pred,ncol=length(xx)),col=colour,

xlab="",ylab="",main="Second Simple coKriging")

# simple kriging map prediction variance of the second variable

image.plot(xx, xx, matrix(pr2$mse,ncol=length(xx)),col=colour,

xlab="",ylab="",main="Std error")

GeoResiduals Computes ﬁtted covariance and/or variogram

Description

The procedure return a GeoFit object associated to the estimated residuals

Usage

GeoResiduals(fit)

Arguments

fit A ﬁtted object obtained from the GeoFit.

Value

A GeoFit object with the estimated residuals

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit.

54 GeoSim

Examples

library(GeoModels)

set.seed(211)

model="Weibull";shape=4

# Set the coordinates of the points:

x <- runif(700, 0, 1)

y <- runif(700, 0, 1)

coords=cbind(x,y)

# Set the model's parameters:

corrmodel <- "Wend0"

mean <- 5

sill <- 1

nugget <- 0

scale <- 0.3

power2=4

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

# Simulation of the Gaussian RF:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, model=model,param=param)$data

start=list(mean=mean,scale=scale,shape=shape)

fixed=list(nugget=nugget,sill=sill,power2=power2)

# Maximum composite-likelihood fitting of the BinomGaussian random field:

fit <- GeoFit(data,coordx=coords, corrmodel=corrmodel,model=model,

likelihood="Marginal",type='Pairwise',start=start,

fixed=fixed,maxdist=0.1)

res=GeoResiduals(fit)

# Empirical estimation of the variogram for the residuals:

vario <- GeoVariogram(res$data,coordx=coords,maxdist=0.5)

# Plot of covariance and variogram functions:

GeoCovariogram(res, show.vario=TRUE, vario=vario,pch=20)

GeoSim Simulation of Gaussian and non Gaussian Random Fields.

Description

Simulation of Gaussian and some non Gaussian spatial, spatio-temporal and spatial bivariate ran-

dom ﬁelds. The function return a realization of a Random Field for a given covariance model and

covariance parameters. Simulation is based on Cholesky decomposition.

Usage

GeoSim(coordx, coordy=NULL, coordt=NULL, coordx_dyn=NULL, corrmodel, distance="Eucl",

GPU=NULL, grid=FALSE, local=c(1,1),method="cholesky", model='Gaussian', n=1, param,

radius=6371, sparse=FALSE, X=NULL)

GeoSim 55

Arguments

coordx A numeric (d×2)-matrix (where dis the number of spatial sites) giving 2-

dimensions of spatial coordinates or a numeric d-dimensional vector giving 1-

dimension of spatial coordinates. Coordinates on a sphere for a ﬁxed radius

radius are passed in lon/lat format expressed in decimal degrees.

coordy A numeric vector giving 1-dimension of spatial coordinates; coordy is inter-

preted only if coordx is a numeric vector or grid=TRUE otherwise it will be

ignored. Optional argument, the default is NULL then coordx is expected to be

numeric a (d×2)-matrix.

coordt A numeric vector giving 1-dimension of temporal coordinates. Optional argu-

ment, the default is NULL then a spatial RF is expected.

coordx_dyn A list of mnumeric (dt×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 De-

tails.

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See the Section Details of GeoFit.

GPU Numeric; if NULL (the default) no GPU computation is performed.

grid Logical; if FALSE (the default) the data are interpreted as spatial or spatial-

temporal realisations on a set of non-equispaced spatial sites (irregular grid).

local Numeric; number of local work-items of the GPU

method String; the type of matrix decomposition used in the simulation. Default is

cholesky. The other possible choices is svd.

model String; the type of RF and therefore the densities associated to the likelihood

objects. Gaussian is the default, see the Section Details.

nNumeric; the number of trials for binomial RFs. The number of successes in the

negative Binomial RFs. Default is 1.

param A list of parameter values required in the simulation procedure of RFs, see Ex-

amples.

radius Numeric; a value indicating the radius of the sphere when using the great circle

distance. Default value is the radius of the earth in Km (i.e. 6371)

sparse Logical; if TRUE then cholesky decomposition is performed using sparse ma-

trices algorithms (spam packake). It should be used with compactly supported

covariance models.FALSE is the default.

XNumeric; Matrix of space-time covariates.

Value

Returns an object of class GeoSim. An object of class GeoSim is a list containing at most the

following components:

bivariate Logical:TRUE if the Gaussian RF is bivariate, otherwise FALSE;

coordx Ad-dimensional vector of spatial coordinates;

coordy Ad-dimensional vector of spatial coordinates;

coordt At-dimensional vector of temporal coordinates;

coordx_dyn A list of dynamical (in time) spatial coordinates;

56 GeoSim

corrmodel The correlation model; see GeoCovmatrix.

data The vector or matrix or array of data, see GeoFit;

distance The type of spatial distance;

model The type of RF, see GeoFit.

nThe number of trial for Binomial RFs;the number of successes in a negative

Binomial RFs;

numcoord The number of spatial coordinates;

numtime The number the temporal realisations of the RF;

param The vector of parameters’ estimates;

radius The radius of the sphere if coordinates are passed in lon/lat format;

randseed The seed used for the random simulation;

spacetime TRUE if spatio-temporal and FALSE if spatial RF;

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

Examples

library(GeoModels)

library(mapproj)

library(fields)

################################################################

###

### Example 1. Simulation of a spatial Gaussian RF on a regular grid

###

###############################################################

# Define the spatial-coordinates of the points:

x <- seq(0,1,0.045)

y <- seq(0,1,0.045)

set.seed(261)

# Simulation of a spatial Gaussian RF with Matern correlation function

data1 <- GeoSim(x,y,grid=TRUE, corrmodel="Matern", param=list(smooth=0.5,

mean=0,sill=1,scale=0.4/3,nugget=0))$data

# Simulation of a spatial Gaussian RF with Generalized Wendland correlation function

# using sparse alghorithm matrices

set.seed(261)

data2 <- GeoSim(x,y,grid=TRUE, corrmodel="GenWend", param=list(smooth=0,

power2=4,mean=0,sill=1,scale=0.4,nugget=0))$data

par(mfrow=c(1,2))

image.plot(x,y,data1,col=terrain.colors(100),main="Matern",xlab="",ylab="")

image.plot(x,y,data2,col=terrain.colors(100),main="Wendland",xlab="",ylab="")

################################################################

###

### Example 2. Simulation of a spatial binomial RF based on

### the latent Gaussian RF with exponential correlation

GeoSim 57

### on a regular grid

###

################################################################

# Define the spatial-coordinates of the points:

x <- seq(0, 1, 0.022)

y <- seq(0, 1, 0.022)

coords <- cbind(x,y)

set.seed(251)

n=5

# Simulation of a spatial Binomial RF:

sim <- GeoSim(x,y,grid=TRUE, corrmodel="Wend0",

model="Binomial",n=n,sparse=TRUE,

param=list(nugget=0,mean=0,scale=.2,sill=1,power2=4))

image.plot(x,y,sim$data,nlevel=n+1,col=terrain.colors(n+1),zlim=c(0,n))

################################################################

###

### Example 3. Simulation of a spatial Weibull RF

### with exponential correlation

###

###############################################################

# Define the spatial-coordinates of the points:

x <- seq(0,1,0.032)

y <- seq(0,1,0.032)

set.seed(261)

# Simulation of a spatial Gaussian RF with Matern correlation function

data1 <- GeoSim(x,y,grid=TRUE, corrmodel="Exponential",model="Weibull",

param=list(shape=1.2,mean=0,sill=1,scale=0.3/3,nugget=0))$data

image.plot(x,y,data1,col=terrain.colors(200),main="Weibull RF",xlab="",ylab="")

################################################################

###

### Example 4. Simulation of a spatial t RF

### with exponential correlation

###

###############################################################

# Define the spatial-coordinates of the points:

x <- seq(0,1,0.03)

y <- seq(0,1,0.03)

set.seed(268)

# Simulation of a spatial Gaussian RF with Matern correlation function

data1 <- GeoSim(x,y,grid=TRUE, corrmodel="GenWend",model="StudentT", sparse=TRUE,

param=list(df=1/4,mean=0,sill=1,scale=0.3,nugget=0,smooth=1,power2=5))$data

image.plot(x,y,data1,col=terrain.colors(100),main="Student-t RF",xlab="",ylab="")

################################################################

###

### Example 5. Simulation of a sinhasinh RF

### with Wend0 correlation.

58 GeoSim

###

###############################################################

# Define the spatial-coordinates of the points:

x <- runif(800, 0, 2)

y <- runif(800, 0, 2)

coords <- cbind(x,y)

set.seed(261)

corrmodel="Wend0"

# Simulation of a spatial Gaussian RF:

param=list(power2=4,skew=0,tail=1,

mean=0,sill=1,scale=0.2,nugget=0) ## gaussian case

data0 <- GeoSim(coordx=coords, corrmodel=corrmodel,

model="SinhAsinh", param=param,sparse=TRUE)$data

plot(density(data0),xlim=c(-7,7))

param=list(power2=4,skew=0,tail=0.7,

mean=0,sill=1,scale=0.2,nugget=0) ## heavy tails

data1 <- GeoSim(coordx=coords, corrmodel=corrmodel,

model="SinhAsinh", param=param,sparse=TRUE)$data

lines(density(data1),lty=2)

param=list(power2=4,skew=0.5,tail=1,

mean=0,sill=1,scale=0.2,nugget=0) ## asymmetry

data2 <- GeoSim(coordx=coords, corrmodel=corrmodel,

model="SinhAsinh", param=param,sparse=TRUE)$data

lines(density(data2),lty=3)

################################################################

###

### Example 6. Simulation of a bivariate Gaussian RF

### with separable bivariate exponential correlation model

### on a regular grid.

###

###############################################################

# Define the spatial-coordinates of the points:

x <- seq(-1,1,0.08)

y <- seq(-1,1,0.08)

# Simulation of a bivariate spatial Gaussian RF:

# with a separable Bivariate Matern

set.seed(12)

param=list(mean_1=0,mean_2=0,scale=0.12,smooth=0.5,

sill_1=1,sill_2=1,nugget_1=0,nugget_2=0,pcol=0.5)

data <- GeoSim(x,y,grid=TRUE,corrmodel="Bi_matern_sep",

param=param)$data

par(mfrow=c(1,2))

image.plot(x,y,data[,,1],col=terrain.colors(100),main="1",xlab="",ylab="")

image.plot(x,y,data[,,2],col=terrain.colors(100),main="2",xlab="",ylab="")

GeoSim 59

################################################################

###

### Example 7. Simulation of a spatio temporal Gaussian RF.

### observed on dynamical location sites with double exponential correlation

###

###############################################################

# Define the dynamical spatial-coordinates of the points:

coordt=1:5

coordx_dyn=list()

maxN=40

set.seed(8)

for(k in 1:length(coordt))

{

NN=sample(1:maxN,size=1)

x <- runif(NN, 0, 1)

y <- runif(NN, 0, 1)

coordx_dyn[[k]]=cbind(x,y)

}

coordx_dyn

param<-list(nugget=0,mean=0,scale_s=0.2/3,scale_t=2/3,sill=1)

data <- GeoSim(coordx_dyn=coordx_dyn, coordt=coordt, corrmodel="Exp_Exp",

param=param)$data

## spatial realization at first temporal instants

data[[1]]

## spatial realization at third temporal instants

data[[3]]

################################################################

###

### Example 8. Simulation of a Gaussian RF

### with a Wend0 correlation on the planet earth

###

###############################################################

require(plot3D)

require(sphereplot)

distance="Geod";radius=6371

NN=4500 ## total point on the sphere on lon/lat format

set.seed(80)

coords=pointsphere(NN,c(-180,180),c(-90,90),c(radius,radius))

## Set the wendland parameters

corrmodel <- "Wend0"

param<-list(mean=0,sill=1,nugget=0,scale=radius*0.3,power2=3)

# Simulation of a spatial Gaussian RF on the sphere

#set.seed(2)

data <- GeoSim(coordx=coords,corrmodel=corrmodel,sparse=TRUE,

distance=distance,radius=radius,param=param)$data

#converting in 3d cartesian coordinates

b=sph2car(coords[,1], coords[,2], radius = radius, deg = TRUE)

60 GeoSim

x0 = b[,1]; y0 = b[,2]; z0 = b[,3]

# plotting

scatter3D(x0,y0,z0,colvar=data,clim=c(min(data),max(data)),pch=20,cex=0.8,

col=rainbow(200))

################################################################

###

### Example 9. Simulation of a Gaussian RF

### with Wend0 model on USA

###

###############################################################

distance="Geod";radius=6378.88

NN=40

x=seq(-125,-64,length.out=NN)

y=seq(27,50, length.out =NN)

nrow(expand.grid(x,y))

## Set the wendland parameters

corrmodel <- "Wend0"

param<-list(mean=0,sill=1,nugget=0,scale=radius*0.3,power2=3)

# Simulation of a spatial Gaussian RF on the sphere

#set.seed(2)

data <- GeoSim(x,y,grid=TRUE,corrmodel=corrmodel,sparse=TRUE,

distance=distance,radius=radius,param=param)$data

image.plot(x,y,data,col=terrain.colors(100),xlab="",ylab="")

map("usa", add = TRUE)

################################################################

###

### Example 10. Simulation of a Wrapped RF

### with exponential correlation

### on a regular grid

###

###############################################################

# Define the spatial-coordinates of the points:

x <- runif(200,0, 1)

y <- runif(200,0, 1)

coords <- cbind(x,y)

set.seed(251)

# Simulation of a spatial wrapped RF:

sim <- GeoSim(coordx=coords, corrmodel="Exp",

model="Wrapped",

param=list(nugget=0,mean=0,scale=.1,sill=1))$data

long <- 0.08;

x1 <- coords[,1] + long*cos(sim)

y1 <- coords[,2] + long*sin(sim)

eps <- 0.1

plot(0,xlim=c(0-eps,1+eps),ylim=c(0-eps,1+eps));

GeoTests 61

require(shape)

Arrows(coords[,1], coords[,2], x1, y1, arr.length = 0.2, code = 2, arr.type = "triangle",col=1)

GeoTests Statistical Hypothesis Tests for Nested Models

Description

The function performs statistical hypothesis tests for nested models based on composite or standard

likelihood versions of: Wald-type and Wilks-type (likelihood ratio) statistics.

Usage

GeoTests(object1, object2, ..., statistic)

Arguments

object1 An object of class GeoFit.

object2 An object of class GeoFit that is a nested model within object1.

... Further successively nested objects.

statistic String; the name of the statistic used within the hypothesis test (see Details).

Details

The implemented hypothesis tests for nested models are based on the following statistics:

1. Wald-type (Wald);

2. Likelihood ratio or Wilks-type (Wilks under standard likelihood); For composite likelihood

available variants of the basic version are:

• Rotnitzky and Jewell adjustment (WilksRJ);

• Satterhwaite adjustment (WilksS);

• Chandler and Bate adjustment (WilksCB);

• Pace, Salvan and Sartori adjustment (WilksPSS);

More speciﬁcally, consider an p-dimensional random vector Ywith probability density function

f(y;θ), where θ∈Θis a q-dimensional vector of parameters. Suppose that θ= (ψ, τ)can

be partitioned in a q0-dimensional subvector ψand q00-dimensional subvector τ. Assume also to be

interested in testing the speciﬁc values of the vector ψ. Then, one can use some statistical hypothesis

tests for testing the null hypothesis H0:ψ=ψ0against the alternative H1:ψ6=ψ0. Composite

likelihood versions of ’Wald’ and ’score’ statistics have the usual asymptotic chi-square distribution

with q0degree of freedom. The Wald-type statistic is

W= ( ˆ

ψ−ψ0)T(Gψψ)−1(ˆ

θ)( ˆ

ψ−ψ0),

where Gψψ is the q0×q0submatrix of the Godambe or Fisher information pertaining to ψand ˆ

θis

the maximum likelihood estimator from the full model. This statistic can be called from the routine

GeoTests assigning at the argument statistic the value: Wald.

62 GeoTests

Alternatively to the Wald-type statistic one can use the composite version of the Wilks-type or

likelihood ratio statistic, given by

W= 2[C`(ˆ

θ;y)−C`{ψ0,ˆτ(ψ0); y}].

In the composite likelihood case, the asymptotic distribution of the composite likelihood ratio statis-

tic is given by

W˙∼X

λiχ2,

for i= 1, . . . , q0, where χ2

iare q0iid copies of a chi-square one random variable and λ1, . . . , λq0

are the eigenvalues of the matrix (Hψψ )−1Gψψ. There exist several adjustments to the composite

likelihood ratio statistic in order to get an approximated χ2

q0. For example, Rotnitzky and Jewell

(1990) proposed the adjustment W0=W/¯

λwhere ¯

λis the average of the eigenvalues λi. This

statistic can be called within the routine by the value: WilksRJ. A better solution is proposed by

Satterhwaite (1946) deﬁning W00 =νW/(q0¯

λ), where ν= (Piλ)2/Piλ2

ifor i= 1 ...,q0, is the

effective number of the degree of freedom. Note that in this case the distribution of the likelihood

ratio statistic is a chi-square random variable with νdegree of freedom. This statistic can be called

from the routine assigning the value: WilksS. For the adjustments suggested by Chandler and Bate

(2007) and Pace, Salvan and Sartori (2011) we refere to the articles (see References), these versions

can be called from the routine assigning respectively the values: WilksCB and WilksPSS.

Value

An object of class c("data.frame"). The object contain a table with the results of the tested

models. The rows represent the responses for each model and the columns the following results:

Num.Par The number of the model’s parameters.

Diff.Par The difference between the number of parameters of the model in the previous

row and those in the actual row.

Df The effective number of degree of freedom of the chi-square distribution.

Chisq The observed value of the statistic.

Pr(>chisq) The p-value of the quantile Chisq computed using a chi-squared distribution

with Df degrees of freedom.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

References

Chandler, R. E., and Bate, S. (2007). Inference for Clustered Data Using the Independence log-

likelihood. Biometrika,94, 167–183.

Pace, L., Salvan, A. and Sartori, N. (2011). Adjusting Composite Likelihood Ratio Statistics. Sta-

tistica Sinica,21, 129–148.

Rotnitzky, A. and Jewell, N. P. (1990). Hypothesis Testing of Regression Parameters in Semipara-

metric Generalized Linear Models for Cluster Correlated Data. Biometrika,77, 485–497.

Satterthwaite, F. E. (1946). An Approximate Distribution of Estimates of Variance Components.

Biometrics Bulletin,2, 110–114.

Varin, C., Reid, N. and Firth, D. (2011). An Overview of Composite Likelihood Methods. Statistica

Sinica,21, 5–42.

GeoTests 63

See Also

GeoFit.

Examples

library(GeoModels)

################################################################

###

### Example 1. Parametric test of Gaussianity

### using SinhAsinh random fields using standard likelihood

###

###############################################################

set.seed(314)

model="SinhAsinh"

# Define the spatial-coordinates of the points:

NN=500

x <- runif(NN, 0, 1)

y <- runif(NN, 0, 1)

coords <- cbind(x,y)

# Parameters

mean=0; nugget=0; sill=1

### skew and tail parameters

skew=0;tail=1 ## H0 i.e Gaussianity

# underlying model correlation

corrmodel="Wend0"

power2=4;c_supp=0.2

# simulation

param=list(power2=power2,skew=skew,tail=tail,

mean=mean,sill=sill,scale=c_supp,nugget=nugget)

data <- GeoSim(coordx=coords, corrmodel=corrmodel,model=model, param=param)$data

##### H1

fixed=list(power2=power2,nugget=nugget,mean=mean)

start=list(scale=c_supp,skew=skew,tail=tail,sill=sill)

# Maximum pairwise composite-likelihood fitting of the RF:

fitH1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,

likelihood="Full",type="Standard",varest=TRUE,

start=start,fixed=fixed)

fitH1$param

##### H0: Gaussianity (i.e tail1=1, skew=0 fixed)

fixed=list(power2=power2,nugget=nugget,mean=mean,tail=1,skew=0)

start=list(scale=c_supp,sill=sill)

# Maximum pairwise composite-likelihood fitting of the RF:

fitH0 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,

likelihood="Full",type="Standard",varest=TRUE,

start=start,fixed=fixed)

# Wald statistic test

GeoTests(fitH1, fitH0 ,statistic='Wald')

64 GeoVarestbootstrap

# likelihood ratio statistic test

GeoTests(fitH1, fitH0 , statistic='Wilks')

################################################################

###

### Example 2. Composite likelihood-based hypothesis testing.

### Testing significance of a covariate parameter

###

###############################################################

set.seed(31)

# Define the spatial-coordinates of the points:

N=1500

x <- runif(N, 0, 1)

y <- runif(N, 0, 1)

X=cbind(rep(1,N),runif(N))

coords=cbind(x,y)

# Set the model's parameters:

corrmodel <- "Exp"

mean <- 1; mean1=0.3 ## H0: mean1=0

sill <- 1

nugget <- 0

scale <- 0.15/3

# Simulation of the spatial Gaussian random field:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, param=list(mean=mean,mean1=mean1,

sill=sill,nugget=nugget,scale=scale),X=X)$data

# Pairwise-likelihood fitting of the random field, full model:

start=list(mean=mean,mean1=mean1,sill=sill,scale=scale)

fixed=list(nugget=nugget)

fitH1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, maxdist=0.01,

varest=TRUE,likelihood="Marginal",type="Pairwise",

fixed=fixed,start=start,X=X)

# Pairwise-likelihood fitting of the random field, with a nested model:

start=list(mean=mean,sill=sill,scale=scale)

fixed=list(nugget=nugget,mean1=0) # setting H0 restriction

fitH0 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, maxdist=0.01,

varest=TRUE,likelihood="Marginal",type="Pairwise",

fixed=fixed,start=start,X=X)

# Hypothesis testing results:

# composite Wald statistic test:

GeoTests(fitH1, fitH0 ,statistic='Wald')

# composite likelihood ratio statistic test with PSS adjustment:

GeoTests(fitH1, fitH0, statistic='WilksPSS')

GeoVarestbootstrap Update a GeoFit object using parametric bootstrap for std error esti-

mation

GeoVariogram 65

Description

The procedure update a GeoFit object estimating stderr estimation using parametric bootstrap.

Usage

GeoVarestbootstrap(fit,K=100,sparse=FALSE,GPU=NULL,local=c(1,1))

Arguments

fit A ﬁtted object obtained from the GeoFit.

KThe number of simulations in the parametric bootstrap.

sparse Logical; if TRUE then cholesky decomposition is performed using sparse matri-

ces algorithms (spam packake).

GPU Numeric; if NULL (the default) no OpenCL computation is performed. The user

can choose the device to be used. Use DeviceInfo() function to see available

devices, only double precision devices are allowed

local Numeric; number of local work-items of the OpenCL setup

Value

Returns an object of class GeoFit.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit.

GeoVariogram Empirical Variogram(variants) estimation

Description

The function returns an empirical estimate of the variogram for spatio (temporal) and bivariate

random ﬁelds.

Usage

GeoVariogram(data, coordx, coordy=NULL, coordt=NULL, coordx_dyn=NULL, cloud=FALSE,

distance='Eucl', grid=FALSE, maxdist=NULL,

maxtime=NULL, numbins=NULL, radius=6371,

type='variogram',bivariate=FALSE)

66 GeoVariogram

Arguments

data Ad-dimensional vector (a single spatial realisation) or a (n×d)-matrix (niid

spatial realisations) or a (d×d)-matrix (a single spatial realisation on regular

grid) or an (d×d×n)-array (niid spatial realisations on regular grid) or a

(t×d)-matrix (a single spatial-temporal realisation) or an (t×d×n)-array

(niid spatial-temporal realisations) or or an (d×d×t×n)-array (a single

spatial-temporal realisation on regular grid) or an (d×d×t×n)-array (niid

spatial-temporal realisations on regular grid). See GeoFit for details.

coordx A numeric (d×2)-matrix (where dis the number of spatial sites) assigning 2-

dimensions of spatial coordinates or a numeric d-dimensional vector assigning

1-dimension of spatial coordinates. Coordinates on a sphere for a ﬁxed radius

radius are passed in lon/lat format expressed in decimal degrees.

coordy A numeric vector assigning 1-dimension of spatial coordinates; coordy is in-

terpreted only if coordx is a numeric vector or grid=TRUE otherwise it will be

ignored. Optional argument, the default is NULL then coordx is expected to be

numeric a (d×2)-matrix.

coordt A numeric vector assigning 1-dimension of temporal coordinates. Optional ar-

gument, the default is NULL then a spatial random ﬁeld is expected.

coordx_dyn A list of mnumeric (dt×2)-matrices containing dynamical (in time) spatial

coordinates. Optional argument, the default is NULL

cloud Logical; if TRUE the variogram cloud is computed, otherwise if FALSE (the de-

fault) the empirical (binned) variogram is returned.

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 are interpreted as spatial or spatial-

temporal realisations on a set of non-equispaced spatial sites.

maxdist A numeric value denoting the spatial maximum distance, see the Section De-

tails.

maxtime A numeric value denoting the temporal maximum distance, see the Section De-

tails.

numbins A numeric value denoting the numbers of bins, see the Section Details.

radius Numeric; a value indicating the radius of the sphere when using the great circle

distance. Default value is the radius of the earth in Km (i.e. 6371)

type A String denoting the type of variogram. The option available is : variogram.

bivariate Logical; if FALSE (the default) the data are interpreted as univariate spatial or

spatial-temporal realisations. Otherwise they are intrepreted as a a realization

from a bivariate ﬁeld.

Details

We brieﬂy report the deﬁnitions of semi-variogram used in this function. In the case of a spatial

Gaussian random ﬁeld the sample variogram estimator is deﬁned by

ˆγ(h)=0.5X

xi,xj∈N(h)

(Z(xi)−Z(xj))2/|N(h)|

where N(h)is the set of all the sample pairs whose distances fall into a tolerance region with size

h(equispaced intervalls are considered).

GeoVariogram 67

In the case of a spatio-temporal Gaussian random ﬁeld the sample variogram estimator is deﬁned

ˆγ(h, u)=0.5X

(xi,l),(xj,k)∈N(h,u)

(Z(xi, l)−Z(xj, k))2/|N(h, u)|

where N(h, u)is the set of all the sample pairs whose spatial distances fall into a tolerance region

with size hand |k−l|=u. Note, that Z(xi, l)is the observation at site xiand time l.

The numbins parameter indicates the number of adjacent intervals to consider in order to grouped

distances with which to compute the (weighted) lest squares.

The maxdist parameter indicates the maximum spatial distance below which the shorter distances

will be considered in the calculation of the (weigthed) least squares.

The maxtime parameter indicates the maximum temporal distance below which the shorter dis-

tances will be considered in the calculation of the (weigthed) least squares.

Value

Returns an object of class Variogram. An object of class Variogram is a list containing at most the

following components:

bins Adjacent intervals of grouped spatial distances if cloud=FALSE. Otherwise if

cloud=TRUE all the spatial pairwise distances;

bint Adjacent intervals of grouped temporal distances if cloud=FALSE. Otherwise if

cloud=TRUE all the temporal pairwise distances;

cloud If the variogram cloud is returned (TRUE) or the empirical variogram (FALSE);

centers The centers of the spatial bins;

distance The type of spatial distance;

lenbins The number of pairs in each spatial bin;

lenbinst The number of pairs in each spatial-temporal bin;

lenbint The number of pairs in each temporal bin;

maxdist The maximum spatial distance used for the calculation of the variogram. If no

spatial distance is speciﬁed then it is NULL;

maxtime The maximum temporal distance used for the calculation of the variogram. If

no temporal distance is speciﬁed then it is NULL;

variograms The empirical spatial variogram;

variogramst The empirical spatial-temporal variogram;

variogramt The empirical temporal variogram;

type The type of estimated variogram

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

References

Cressie, N. A. C. (1993) Statistics for Spatial Data. New York: Wiley.

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

68 GeoVariogram

See Also

GeoFit

Examples

library(GeoModels)

################################################################

###

### Example 1. Empirical estimation of the semi-variogram from a

### spatial Gaussian random field with exponential correlation.

###

###############################################################

set.seed(514)

# Set the coordinates of the sites:

x <- runif(200, 0, 1)

y <- runif(200, 0, 1)

coords <- cbind(x,y)

# Set the model's parameters:

corrmodel <- "exponential"

mean <- 0

sill <- 1

nugget <- 0

scale <- 0.3/3

# Simulation of the spatial Gaussian random field:

data <- GeoSim(coordx=coords, corrmodel=corrmodel, param=list(mean=mean,

sill=sill, nugget=nugget, scale=scale))$data

# Empirical spatial semi-variogram estimation:

fit <- GeoVariogram(coordx=coords,data=data,maxdist=0.6)

# Results:

plot(fit$centers, fit$variograms, xlab='h', ylab=expression(gamma(h)),

ylim=c(0, max(fit$variograms)), pch=20,

main="Semi-variogram")

################################################################

###

### Example 2. Empirical estimation of the variogram from a

### spatio-temporal Gaussian random fields with Gneiting

### correlation function.

###

###############################################################

set.seed(331)

# Define the temporal sequence:

# Set the coordinates of the sites:

x <- runif(400, 0, 1)

y <- runif(400, 0, 1)

coords <- cbind(x,y)

times <- seq(1,5,1)

# Simulation of a spatio-temporal Gaussian random field:

data <- GeoSim(coordx=coords, coordt=times, corrmodel="gneiting",

GeoVariogram 69

param=list(mean=0,scale_s=0.1,scale_t=0.1,sill=1,

nugget=0,power_s=1,power_t=1,sep=0.5))$data

# Empirical spatio-temporal semi-variogram estimation:

fit <- GeoVariogram(data=data, coordx=coords, coordt=times, maxtime=5,maxdist=0.5)

# Results: Marginal spatial empirical semi-variogram

par(mfrow=c(2,2), mai=c(.5,.5,.3,.3), mgp=c(1.4,.5, 0))

plot(fit$centers, fit$variograms, xlab='h', ylab=expression(gamma(h)),

ylim=c(0, max(fit$variograms)), xlim=c(0, max(fit$centers)),

pch=20,main="Marginal spatial semi-variogram",cex.axis=.8)

# Results: Marginal temporal empirical semi-variogram

plot(fit$bint, fit$variogramt, xlab='t', ylab=expression(gamma(t)),

ylim=c(0, max(fit$variogramt)),xlim=c(0,max(fit$bint)),

pch=20,main="Marginal temporal semi-variogram",cex.axis=.8)

# Building space-time semi-variogram

st.vario <- matrix(fit$variogramst,length(fit$centers),length(fit$bint))

st.vario <- cbind(c(0,fit$variograms), rbind(fit$variogramt,st.vario))

# Results: 3d Spatio-temporal semi-variogram

require(scatterplot3d)

st.grid <- expand.grid(c(0,fit$centers),c(0,fit$bint))

scatterplot3d(st.grid[,1], st.grid[,2], c(st.vario),

highlight.3d=TRUE, xlab="h",ylab="t",

zlab=expression(gamma(h,t)), pch=20,

main="Space-time semi-variogram",cex.axis=.7,

mar=c(2,2,2,2), mgp=c(0,0,0),

cex.lab=.7)

# A smoothed version

par(mai=c(.2,.2,.2,.2),mgp=c(1,.3, 0))

persp(c(0,fit$centers), c(0,fit$bint), st.vario,

xlab="h", ylab="u", zlab=expression(gamma(h,u)),

ltheta=90, shade=0.75, ticktype="detailed", phi=30,

theta=30,main="Space-time semi-variogram",cex.axis=.8,

cex.lab=.8)

################################################################

###

### Example 3. Empirical estimation of the (cross) semivariograms

### from a bivariate Gaussian random fields with Matern

### correlation function.

###

###############################################################

# Simulation of a bivariate spatial Gaussian random field:

set.seed(29)

# Define the spatial-coordinates of the points:

x <- runif(200, 0, 1)

set.seed(7)

y <- runif(200, 0, 1)

coords=cbind(x,y)

# Simulation of a bivariate Gaussian Random field

# with matern (cross) covariance function

param=list(mean_1=0,mean_2=0,scale_1=0.15/3,scale_2=0.2/3,scale_12=0.15/3,

70 GeoWLS

sill_1=1,sill_2=1,nugget_1=0,nugget_2=0,

smooth_1=0.5,smooth_12=0.5,smooth_2=0.5,pcol=-0.45)

data <- GeoSim(coordx=coords, corrmodel="Bi_matern", param=param)$data

# Empirical semi-(cross)variogram estimation:

biv_vario=GeoVariogram(data,coordx=coords, bivariate=TRUE,maxdist=c(0.5,0.5,0.5))

# Variograms plots

par(mfrow=c(2,2))

plot(biv_vario$centers,biv_vario$variograms[1,],pch=20,xlab="h",ylim=c(0,1.2),

ylab="",main=expression(gamma[11](h)))

plot(biv_vario$centers,biv_vario$variogramst,pch=20,xlab="h",

ylab="",main=expression(gamma[12](h)))

plot(biv_vario$centers,biv_vario$variogramst,pch=20,xlab="h",ylab="",

main=expression(gamma[21](h)))

plot(biv_vario$centers,biv_vario$variograms[2,],pch=20,xlab="h",,ylim=c(0,1.2),

ylab="",main=expression(gamma[22](h)))

GeoWLS WLS of Random Fields

Description

the function returns the parameters’ estimates and the estimates’ variances of a random ﬁeld ob-

tained by the weigthed least squares estimator.

Usage

GeoWLS(data, coordx, coordy=NULL, coordt=NULL, coordx_dyn=NULL, corrmodel,

distance="Eucl", fixed=NULL, grid=FALSE, maxdist=NULL,

maxtime=NULL, model='Gaussian', optimizer='Nelder-Mead',

numbins=NULL, radius=6371, start=NULL, weighted=FALSE)

Arguments

data Ad-dimensional vector (a single spatial realisation) or a (d×d)-matrix (a single

spatial realisation on regular grid) or an (d×d×n)-array (niid spatial realisa-

tions on regular grid) or a (t×d)-matrix (a single spatial-temporal realisation)

or an (d×d×t×n)-array (a single spatial-temporal realisation on regular grid).

See GeoFit for details.

coordx A numeric (d×2)-matrix (where dis the number of spatial sites) giving 2-

dimensions of spatial coordinates or a numeric d-dimensional vector giving 1-

dimension of spatial coordinates.

coordy A numeric vector giving 1-dimension of spatial coordinates; coordy is inter-

preted only if coordx is a numeric vector or grid=TRUE otherwise it will be

ignored. Optional argument, the default is NULL then coordx is expected to be

numeric a (d×2)-matrix.

coordt A numeric vector giving 1-dimension of temporal coordinates. Optional argu-

ment, the default is NULL then a spatial random ﬁeld is expected.

coordx_dyn A list of mnumeric (dt×2)-matrices containing dynamical (in time) spatial

coordinates. Optional argument, the default is NULL

GeoWLS 71

corrmodel String; the name of a correlation model, for the description (see GeoFit).

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See the Section Details of GeoFit.

fixed A named list giving the values of the parameters that will be considered as

known values. The listed parameters for a given correlation function will be

not estimated, i.e. if list(nugget=0) the nugget effect is ignored.

grid Logical; if FALSE (the default) the data are interpreted as a vector or a (n×d)-

matrix, instead if TRUE then (d×d×n)-matrix is considered.

maxdist A numeric value denoting the maximum distance, see Details and GeoFit.

maxtime Numeric; an optional positive value indicating the maximum temporal lag con-

sidered in the composite-likelihood computation (see GeoFit.

model String; the type of random ﬁeld. Gaussian is the default, see GeoFit for the

different types.

optimizer String; the optimization algorithm (see optim for details). ’Nelder-Mead’ is the

default.

numbins A numeric value denoting the numbers of bins, see the Section Details

radius Numeric; a value indicating the radius of the sphere when using the great circle

distance. Default value is the radius of the earth in Km (i.e. 6371)

start A named list with the initial values of the parameters that are used by the nu-

merical routines in maximization procedure. NULL is the default (see GeoFit).

weighted Logical; if TRUE then the weighted least square estimator is considered. If FALSE

(the default) then the classic least square is used.

Details

The numbins parameter indicates the number of adjacent intervals to consider in order to grouped

distances with which to compute the (weighted) lest squares.

The maxdist parameter indicates the maximum distance below which the shorter distances will be

considered in the calculation of the (weigthed) least squares.

Value

Returns an object of class WLS. An object of class WLS is a list containing at most the following

components:

bins Adjacent intervals of grouped distances;

bint Adjacent intervals of grouped temporal separations

centers The centers of the bins;

coordx The vector or matrix of spatial coordinates;

coordy The vector of spatial coordinates;

coordt The vector of temporal coordinates;

convergence A string that denotes if convergence is reached;

corrmodel The correlation model;

data The vector or matrix of data;

distance The type of spatial distance;

fixed The vector of ﬁxed parameters;

72 GeoWLS

iterations The number of iteration used by the numerical routine;

maxdist The maximum spatial distance used for the calculation of the variogram used in

least square estimation. If no spatial distance is speciﬁed then it is NULL;

maxtime The maximum temporal distance used for the calculation of the variogram used

in least square estimation. If no temporal distance is speciﬁed then it is NULL;

message Extra message passed from the numerical routines;

model The type of random ﬁelds;

numcoord The number of spatial coordinates;

numtime The number the temporal realisations of the random ﬁeld;

param The vector of parameters’ estimates;

variograms The empirical spatial variogram;

variogramt The empirical temporal variogram;

variogramst The empirical spatial-temporal variogram;

weighted A logical value indicating if its the weighted method;

wls The value of the least squares at the minimum.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

References

Cressie, N. A. C. (1993) Statistics for Spatial Data. New York: Wiley.

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

See Also

GeoFit,optim

Examples

library(GeoModels)

# Set the coordinates of the sites:

set.seed(211)

x <- runif(200, 0, 1)

set.seed(98)

y <- runif(200, 0, 1)

coords <- cbind(x,y)

################################################################

###

### Example 1. Least square fitting of a Gaussian random field

### with exponential correlation.

###

###############################################################

# Set the model's parameters:

Lik 73

corrmodel <- "Exponential"

mean <- 0

sill <- 1

nugget <- 0

scale <- 0.15/3

param <- list(mean=0,sill=sill, nugget=nugget, scale=scale)

# Simulation of the Gaussian random field:

set.seed(2)

data <- GeoSim(coordx=coords, corrmodel=corrmodel, param=param)$data

fixed=list(nugget=0,mean=mean)

start=list(scale=scale,sill=sill)

# Least square fitting of the random field:

fit <- GeoWLS(data=data,coordx=coords, corrmodel=corrmodel,

fixed=fixed,start=start,maxdist=0.5)

# Results:

print(fit)

################################################################

###

### Example 3. Least square fitting of a spatio-temporal

### Gaussian random field with double exponential correlation.

###

###############################################################

# Define the temporal sequence:

time <- seq(1, 10, 1)

mean <- 0

sill <- 1

scale_s <- 0.15/3

scale_t <- 2/3

param <- list(mean=0,scale_s=scale,scale_t=scale_t,sill=sill,nugget=nugget)

# Simulation of the Gaussian random field:

set.seed(35)

data <- GeoSim(coordx=coords,coordt=time, corrmodel="exp_exp",

param=param)$data

fixed<-list(nugget=nugget,mean=0)

start<-list(scale_s=scale_s,scale_t=scale_t,sill=1)

# Weighted least square estimation:

fit <- GeoWLS(data=data, coordx=coords,coordt=time, corrmodel="exp_exp",

,maxdist=0.5,maxtime=3,fixed=fixed,start=start)

# Results

print(fit)

Lik Optimizes the Log Likelihood

74 Lik

Description

Subroutine called by GeoFit. The procedure estimates the model parameters by maximization of

the log-likelihood.

Usage

Lik(bivariate,coordx,coordy,coordt,coordx_dyn,corrmodel,data,fixed,flagcor,flagnuis,

grid,lower,mdecomp,model,namescorr,namesnuis,namesparam,numcoord,

numpairs,numparamcor,numtime,optimizer,onlyvar,param,radius,setup,

spacetime,sparse,varest,taper,type,upper,ns,X)

Arguments

bivariate Logical; if TRUE then the data come froma a bivariate random ﬁeld. Otherwise

from a univariate random ﬁeld.

coordx A numeric (d×2)-matrix (where dis the number of spatial sites) assigning 2-

dimensions of spatial coordinates or a numeric d-dimensional vector assigning

1-dimension of spatial coordinates.

coordy A numeric vector assigning 1-dimension of spatial coordinates; coordy is in-

terpreted only if coordx is a numeric vector or grid=TRUE otherwise it will be

ignored. Optional argument, the default is NULL then coordx is expected to be

numeric a (d×2)-matrix.

coordt A numeric vector assigning 1-dimension of temporal coordinates. Optional ar-

gument, the default is NULL then a spatial random ﬁeld is expected.

coordx_dyn A list of mnumeric (dt×2)-matrices containing dynamical (in time) spatial

coordinates. Optional argument, the default is NULL

corrmodel Numeric; the id of the correlation model.

data A numeric vector or a (n×d)-matrix or (d×d×n)-matrix of observations.

flagcor A numeric vector of ﬂags denoting which correlation parameters have to be

estimated.

flagnuis A numeric verctor of ﬂags denoting which nuisance parameters have to esti-

mated.

fixed A numeric vector of parameters that will be considered as known values.

grid Logical; if FALSE (the default) the data are interpreted as a vector or a (n×d)-

matrix, instead if TRUE then (d×d×n)-matrix is considered.

lower A numeric vector with the lower bounds of the parameters’ ranges.

model Numeric; the id value of the density associated to the likelihood objects.

namescorr String; the names of the correlation parameters.

namesnuis String; the names of the nuisance parameters.

namesparam String; the names of the parameters to be maximised.

numcoord Numeric; the number of coordinates.

numpairs Numeric; the number of pairs.

numparamcor Numeric; the number of the correlation parameters.

numtime Numeric; the number of temporal observations.

mdecomp String; the type of matrix decomposition used in the simulation. Default is

cholesky. The other possible choices is svd (Singular values decomposition).

MatDecomp 75

optimizer String; the optimization algorithm (see optim for details). ’Nelder-Mead’ is the

default.

onlyvar Logical; if TRUE (and varest is TRUE) only the variance covariance matrix is

computed without optimizing. FALSE is the default.

param A numeric vector of parameters.

sparse Logical; if TRUE then maximum likelihood is computed using sparse matrices

algorithms.FALSE is the default.

radius Numeric; the radius of the sphere when considering data on a sphere.

ns Numeric: vector of number of location sites for each temporal instants

setup A List of useful components for the estimation based on the maximum tapered

likelihood.

spacetime Logical; if the random ﬁeld is spatial (FALSE) or spatio-temporal (TRUE).

varest Logical; if TRUE the estimate’ variances and standard errors are returned. FALSE

is the default.

taper String; the name of the taper correlation function.

type String; the type of the likelihood objects. If Pairwise (the default) then the

marginal composite likelihood is formed by pairwise marginal likelihoods.

upper A numeric vector with the upper bounds of the parameters’ ranges.

XNumeric; Matrix of spatio(temporal)covariates in the linear mean speciﬁcation.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

MatDecomp Matrix decomposition

Description

Matrix decomposition.

Usage

MatDecomp(mtx, method)

Arguments

mtx numeric; a square positive or semipositive deﬁnite matrix.

method string; the type of matrix decomposition. Two possible choices: cholesky and

svd.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

76 MatSqrt, MatInv, MatLogDet

MatSqrt, MatInv, MatLogDet

Square root, inverse and log determinant of a (semi)positive deﬁnite

matrix, given a matrix decomposition.

Description

Square root, inverse and log determinant of a (semi)positive deﬁnite matrix, given a matrix decom-

position.

Usage

MatSqrt(mat.decomp,method)

MatInv(mat.decomp,method)

MatLogDet(mat.decomp,method)

Arguments

mat.decomp numeric; a matrix decomposition.

method string; the type of matrix decomposition. Two possible choices: cholesky and

svd.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

MatDecomp

Examples

library(GeoModels)

################################################################

###

### Example 1. Inverse of Covariance matrix associated to

### a Matern correlation model

###

###############################################################

# Define the spatial-coordinates of the points:

x <- runif(15, 0, 1)

y <- runif(15, 0, 1)

coords <- cbind(x,y)

# Matern Parameters

param=list(smooth=0.5,sill=1,scale=0.2,nugget=0)

a=matrix <- GeoCovmatrix(coordx=coords, corrmodel="Matern", param=param)

## decomposition with cholesky method

b=MatDecomp(a$covmat,method="cholesky")

## inverse of covariance matrix

inverse=MatInv(b,method="cholesky")

NuisParam 77

NuisParam Lists the Nuisance Parameters of a Random Field

Description

Subroutine called by InitParam and other procedures. The procedure returns a list with the nuisance

parameters of a given random ﬁeld model.

Usage

NuisParam(model, bivariate,num_betas)

Arguments

model String; the name of a random ﬁeld.

bivariate Logical; if FALSE (the default) the correlation model is univariate spatial or

spatial-temporal. Otherwise is bivariate.

num_betas Numerical; the nunber of mean parameters in the linear speciﬁcation.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

Examples

library(GeoModels)

NuisParam("Gaussian", FALSE,1)

## note that in the bivariate case sill and nugget are considered as correlation parameteres

NuisParam("Gaussian", TRUE,1)

NuisParam("BinomialGaussian", FALSE,1)

NuisParam("Chisq", FALSE,2)

NuisParam("SkewGaussian", FALSE,3)

NuisParam("SinhAsinhGaussian",FALSE,1)

78 Prscores

Prscores Computation of three predictive scores: RMSE, LSCORE, CRPS for

spatial, spatiotemporal and bivariate Gaussian RF.

Description

The function computes RMSE, LSCORE, CRPS predictive scores.

Usage

Prscores(data, method="cholesky", matrix)

Arguments

data Ad-dimensional vector (a single spatial realisation) or a a(t×d)-matrix (a single

spatial-temporal realisation). or a a(2×d)-matrix (a single bivariate realisation).

method String; the type of matrix decomposition used in the computation of the predic-

tive scores. Default is cholesky. The other possible choices is svd.

matrix An object of class matrix. See the Section Details.

Details

For a given covariance matrix object (GeoCovmatrix) and a given spatial, spatiotemporal or bivari-

are realization from a Gaussian random ﬁeld, the function computes three predictive scores.

Value

Returns a list containing the following informations:

RMSE Root-mean-square error predictive score

LSCORE Logarithmic predictive score

CRPS Continuous ranked probability predictive score

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

References

Zhang H. and Wang Y. (2010). Kriging and cross-validation for massive spatial data. Environ-

metrics, 21, 290–304.

Gneiting T. and Raftery A. Strictly Proper Scoring Rules, Prediction, and Estimation. Journal of

the American Statistical Association, 102

See Also

GeoCovmatrix

StartParam 79

Examples

library(GeoModels)

library(fields)

################################################################

######### Examples of predictive score computation ############

################################################################

# Define the spatial-coordinates of the points:

x <- runif(500, 0, 2)

y <- runif(500, 0, 2)

coords=cbind(x,y)

matrix1 <- GeoCovmatrix(coordx=coords, corrmodel="Matern", param=list(smooth=0.5,

sill=1,scale=0.2,nugget=0))

data <- GeoSim(coordx=coords, corrmodel="Matern", param=list(mean=0,smooth=0.5,

sill=1,scale=0.2,nugget=0))$data

Pr_scores <- Prscores(data,matrix=matrix1)

Pr_scores

StartParam Initializes the Parameters for Estimation Procedures

Description

Subroutine called by the ﬁtting procedures. The procedure initializes the parameters for the ﬁtting

procedure.

Usage

StartParam(coordx, coordy, coordt,coordx_dyn, corrmodel, data, distance, fcall,

fixed, grid,likelihood, maxdist, maxtime, model, n, param,

parscale, paramrange, radius, start, taper, tapsep,

type,typereal,varest, vartype, weighted, winconst,

winstp, winconst_t, winstp_t,X)

Arguments

coordx A numeric (d×2)-matrix (where dis the number of points) assigning 2-dimensions

of coordinates or a numeric vector assigning 1-dimension of coordinates.

coordy A numeric vector assigning 1-dimension of coordinates; coordy is interpreted

only if coordx is a numeric vector otherwise it will be ignored.

coordt A numeric vector assigning 1-dimension of temporal coordinates.

coordx_dyn A list of mnumeric (dt×2)-matrices containing dynamical (in time) spatial

coordinates. Optional argument, the default is NULL

corrmodel String; the name of a correlation model.

data A numeric vector or a (n×d)-matrix or (d×d×n)-matrix of observations.

80 StartParam

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See the Section Details.

fcall String; "ﬁtting" to call the ﬁtting procedure and "simulation" to call the simula-

tion procedure.

fixed A named list giving the values of the parameters that will be considered as

known values.

grid Logical; if FALSE (the default) the data are interpreted as a vector or a (n×d)-

matrix, instead if TRUE then (d×d×n)-matrix is considered.

likelihood String; the conﬁguration of the composite likelihood.

maxdist Numeric; an optional positive value indicating the maximum spatial distance

considered in the composite-likelihood computation.

maxtime Numeric; an optional positive value indicating the maximum temporal lag con-

sidered in the composite-likelihood computation.

radius Numeric; the radius of the sphere in the case of lon-lat coordinates. The default

is 6371, the radius of the earth.

model String; the density associated to the likelihood objects. Gaussian is the default.

nNumeric; number of trials for binomial random ﬁelds.

param A numeric vector of parameter values required in the simulation procedure of

random ﬁelds.

parscale A numeric vector of scaling factor to improve the maximizing procedure, see

optim.

paramrange A numeric vector of parameters ranges, see optim.

start A named list with the initial values of the parameters that are used by the nu-

merical routines in maximization procedure.

taper String; the name of the type of covariance matrix. It can be Standard (the

default value) or Tapering for taperd covariance matrix.

tapsep Numeric; an optional value indicating the separabe parameter in the space time

adaptive taper (see Details).

type String; the type of likelihood objects. Temporary value set to be "WLeast-

Square" (weigthed least-square) in order to compute the starting values.

typereal String; the real type of likelihood objects. See GeoFit.

varest Logical; if TRUE the estimates’ variances and standard errors are returned. FALSE

is the default.

vartype String; the type of estimation method for computing the estimate variances, see

the Section Details.

weighted Logical; if TRUE the likelihood objects are weighted, see GeoFit.

winconst Numeric; a positive value for computing the spatial sub-window in the sub-

sampling procedure.

winstp Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

spatial sub-sampling procedure.

winconst_t Numeric; a positive value for computing the temporal sub-window in the sub-

sampling procedure.

winstp_t Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

temporal sub-sampling procedure.

XNumeric; Matrix of space-time covariates.

winds 81

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit

winds Irish Daily Wind Speeds

Description

A matrix containing daily wind speeds, in kilometers per hour, from 1961 to 1978 at 12 sites in

Ireland.

Usage

data(irishwinds)

Format

A (6574 ×11)-matrix containing wind speed observations.

Source

Haslett, J. and Raftery, A. E. (1989), Space-time modelling with long-memory dependence: assess-

ing Ireland’s wind-power resource (with discussion), Applied Statistics, 38, 1–50.

winds.coords Weather Stations of the Irish Daily Wind Speeds

Description

A data frame containing information about the weather stations where the data are recorded in

Ireland.

Usage

data(irishwinds)

Format

A data frame containing site - the name of the city (character), abbr - the abbrevation (character),

elev - the elevation (numeric), lat - latitude (numeric) and lon - longitude.

Source

Haslett, J. and Raftery, A. E. (1989), Space-time modelling with long-memory dependence: assess-

ing Ireland’s wind-power resource (with discussion), Applied Statistics, 38, 1–50.

82 WlsStart

WlsStart Computes Starting Values based on Weighted Least Squares

Description

Subroutine called by GeoFit. The function returns opportune starting values for the composite-

likelihood ﬁtting procedure based on weigthed least squares.

Usage

WlsStart(coordx, coordy, coordt, coordx_dyn, corrmodel, data, distance, fcall,

fixed, grid,likelihood, maxdist, maxtime, model, n, param,

parscale, paramrange, radius,start, taper, tapsep, type, varest,

vartype, weighted, winconst,winconst_t, winstp_t, winstp,X)

Arguments

coordx A numeric (d×2)-matrix (where dis the number of points) assigning 2-dimensions

of coordinates or a numeric vector assigning 1-dimension of coordinates.

coordy A numeric vector assigning 1-dimension of coordinates; coordy is interpreted

only if coordx is a numeric vector otherwise it will be ignored.

coordt A numeric vector assigning 1-dimension of temporal coordinates.

coordx_dyn A list of mnumeric (dt×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.

data A numeric vector or a (n×d)-matrix or (d×d×n)-matrix of observations.

distance String; the name of the spatial distance. The default is Eucl, the euclidean

distance. See the Section Details.

fcall String; "ﬁtting" to call the ﬁtting procedure and "simulation" to call the simula-

tion procedure.

fixed A named list giving the values of the parameters that will be considered as

known values.

grid Logical; if FALSE (the default) the data are interpreted as a vector or a (n×d)-

matrix, instead if TRUE then (d×d×n)-matrix is considered.

likelihood String; the conﬁguration of the composite likelihood.

maxdist Numeric; an optional positive value indicating the maximum spatial distance

considered in the composite-likelihood computation.

maxtime Numeric; an optional positive value indicating the maximum temporal separa-

tion considered in the composite-likelihood computation.

model String; the name of the model. Here the default is NULL.

nNumeric; number of trials in a binomial random ﬁeld.

param A numeric vector of parameter values required in the simulation procedure of

random ﬁelds.

parscale A numeric vector with scaling values for improving the maximisation routine.

paramrange A numeric vector with the range of the parameter space.

WlsStart 83

radius Numeric; a value indicating the radius of the sphere when using the great circle

distance. Default value is the radius of the earth in Km (i.e. 6371)

start A numeric vector with starting values.

taper String; the name of the type of covariance matrix. It can be Standard (the

default value) or Tapering for taperd covariance matrix.

tapsep Numeric; an optional value indicating the separabe parameter in the space time

quasi taper (see Details).

type String; the type of estimation method.

varest Logical; if TRUE the estimates’ variances and standard errors are returned. FALSE

is the default.

vartype String; the type of estimation method for computing the estimate variances, see

the Section Details.

weighted Logical; if TRUE the likelihood objects are weighted, see GeoFit.

winconst Numeric; a positive value for computing the spatial sub-window in the sub-

sampling procedure.

winstp Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

spatial sub-sampling procedure.

winconst_t Numeric; a positive value for computing the temporal sub-window in the sub-

sampling procedure.

winstp_t Numeric; a value in (0,1] for deﬁning the the proportion of overlapping in the

temporal sub-sampling procedure.

XNumeric; Matrix of spatio(temporal)covariates in the linear mean speciﬁcation.

Author(s)

Moreno Bevilacqua, <moreno.bevilacqua@uv.cl>,https://sites.google.com/a/uv.cl/moreno-bevilacqua/

home, Víctor Morales Oñate, <victor.morales@uv.cl>

See Also

GeoFit.

Index

∗Topic Composite

CheckBiv,3

CheckDistance,3

CheckSph,4

CheckST,4

CkCorrModel,5

CkInput,5

CkLikelihood,7

CkModel,7

CkType,8

CkVarType,8

CompLik,9

CorrelationPar,11

CorrParam,11

GeoCovariogram,13

GeoFit,29

GeoKrig,43

GeoResiduals,53

GeoVarestbootstrap,64

Lik,73

MatDecomp,75

MatSqrt, MatInv, MatLogDet,76

NuisParam,77

StartParam,79

∗Topic Devices

DeviceInfo,12

∗Topic LeastSquare

GeoWLS,70

WlsStart,82

∗Topic Predictive scores

Prscores,78

∗Topic Simulation

GeoCovmatrix,18

GeoSim,54

∗Topic Variogram

GeoVariogram,65

∗Topic datasets

anomalies,2

winds,81

winds.coords,81

∗Topic spatial

GeoTests,61

anomalies,2

CheckBiv,3

CheckDistance,3

CheckSph,4

CheckST,4

CkCorrModel,5

CkInput,5

CkLikelihood,7

CkModel,7

CkType,8

CkVarType,8

CompLik,9

CorrelationPar,11

CorrParam,11,20,33,45

DeviceInfo,12

GeoCovariogram,13

GeoCovmatrix,3,5,11,18,32–34,45,46,56,

GeoFit,4,6–11,13,14,19,25,29,44,46,53,

55,56,63,65,66,68,70–72,75,77,

80,81,83

GeoKrig,25,43

GeoResiduals,53

GeoSim,6,25,54

GeoTests,61

GeoVarestbootstrap,64

GeoVariogram,14,65

GeoWLS,13,70

Lik,73

MatDecomp,75,76

MatInv (MatSqrt, MatInv, MatLogDet),76

MatLogDet (MatSqrt, MatInv, MatLogDet),

MatSqrt (MatSqrt, MatInv, MatLogDet),76

MatSqrt, MatInv, MatLogDet,76

NuisParam,33,77

optim,6,10,30,71,72,75,80

print.GeoFit (GeoFit),29

print.GeoSim (GeoSim),54

INDEX 85

print.GeoWLS (GeoWLS),70

Prscores,78

StartParam,79

winds,81

winds.coords,81

WlsStart,82

Geo S Manual

Geos-manual

Navigation menu

Versions of this User Manual:

Views

Navigation