GeoFit2: Max-Likelihood-Based Fitting of Gaussian and non Gaussian RFs.

Description

Maximum weighted composite-likelihood fitting for Gaussian and some Non-Gaussian univariate spatial, spatio-temporal and bivariate spatial RFs. A first preliminary estimation is performed using independence composite-likelihood for the marginal parameters of the model. The estimates are then used as starting values in the second final estimation step. The function allows fixing any of the parameters and setting upper/lower bounds in the optimization.

Usage

GeoFit2(data, coordx, coordy=NULL, coordz=NULL, coordt=NULL,
        coordx_dyn=NULL, copula=NULL, corrmodel, distance="Eucl",
        fixed=NULL, anisopars=NULL, est.aniso=c(FALSE,FALSE),
        grid=FALSE, likelihood="Marginal", lower=NULL, maxdist=Inf,
        neighb=NULL, p_neighb=1, maxtime=Inf, memdist=TRUE,
        method="cholesky", model="Gaussian", n=1, onlyvar=FALSE,
        optimizer="Nelder-Mead", radius=1, score=FALSE,
        sensitivity=FALSE, sparse=FALSE, start=NULL,
        thin_method="iid", type="Pairwise", upper=NULL,
        varest=FALSE, weighted=FALSE, X=NULL, nosym=FALSE,
        spobj=NULL, spdata=NULL)

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 AIC.
coordx: A \(d\)-dimensional vector of spatial coordinates.
coordy: A \(d\)-dimensional vector of spatial coordinates.
coordt: A \(t\)-dimensional vector of temporal coordinates.
coordx_dyn: A list of dynamical (in time) spatial coordinates.
conf.int: Confidence intervals for standard maximum likelihood estimation.
convergence: A string that denotes if convergence is reached.
copula: The type of copula.
corrmodel: The correlation model.
data: The vector or matrix or array (or list) of data.
distance: The type of spatial distance.
fixed: A list of fixed parameters.
iterations: The number of iterations used by the numerical routine.
likelihood: The configuration of the composite likelihood.
logCompLik: The value of the log composite-likelihood at the maximum.
maxdist: The maximum spatial distance used in the weighted composite likelihood (or NULL).
maxtime: The order of temporal neighborhood in the composite likelihood computation.
message: Extra message passed from the numerical routines.
model: The density associated to the likelihood objects.
missp: TRUE if a misspecified Gaussian model is used in the composite likelihood.
n: The number of trials in a binomial RF; the number of successes in a negative binomial RF.
neighb: The order of spatial neighborhood in the composite likelihood computation.
ns: The number of (different) location sites in the bivariate case.
numcoord: The number of spatial coordinates.
numtime: The number of temporal realisations of the RF.
param: A list of parameter estimates.
radius: The radius of the sphere in the case of great-circle distance.
stderr: Standard errors for standard maximum likelihood estimation.
sensmat: The sensitivity matrix.
varcov: The variance-covariance matrix of the estimates.
type: The type of the likelihood objects.
X: The matrix of covariates.

Arguments

data: A \(d\)-dimensional vector (a single spatial realisation) or a (\(d \times d\))-matrix (a single spatial realisation on regular grid) or a (\(t \times d\))-matrix (a single spatial-temporal realisation) or an (\(d \times d \times t \times n\))-array (a single spatial-temporal realisation on regular grid). For the description see the Section Details.
coordx: A numeric (\(d \times 2\))-matrix or (\(d \times 3\))-matrix. Coordinates on a sphere for a fixed radius radius are passed in lon/lat format expressed in decimal degrees.
coordy: A numeric vector giving 1-dimension of spatial coordinates; optional argument, default is NULL.
coordz: A numeric vector giving 1-dimension of spatial coordinates; optional argument, default is NULL.
coordt: A numeric vector assigning 1-dimension of temporal coordinates. Optional argument, default is NULL: if NULL then a spatial RF is expected.
coordx_dyn: A list of \(m\) numeric (\(d_t \times 2\))-matrices containing dynamical (in time) spatial coordinates. Optional argument, default is NULL.
copula: String; the type of copula. It can be "Clayton" or "Gaussian".
corrmodel: String; the name of a correlation model, for the description see the Section Details.
distance: String; the name of the spatial distance. Default is "Eucl" (euclidean distance). See the Section Details.
fixed: An optional named list giving the values of the parameters that will be considered as known values. The listed parameters for a given correlation function will not be estimated.
anisopars: A list of two elements: "angle" and "ratio", i.e. the anisotropy angle and the anisotropy ratio, respectively.
est.aniso: A bivariate logical vector providing which anisotropic parameters must be estimated.
grid: Logical; if FALSE (default) the data are interpreted as spatial or spatial-temporal realisations on a set of non-equispaced spatial sites (irregular grid).
likelihood: String; the configuration of the composite likelihood. "Marginal" is the default; see the Section Details.
lower: An optional named list giving lower bounds for parameters when the optimizer is L-BFGS-B, nlminb, bobyqa or optimize. Names must match those in start.
maxdist: Numeric; an optional positive value indicating the maximum spatial distance considered in the composite likelihood computation. See the Section Details for more information.
neighb: Numeric; an optional positive integer indicating the order of neighborhood in the composite likelihood computation. See the Section Details for more information.
p_neighb: Numeric; a value in \((0,1]\) specifying the expected fraction of nearest-neighbor pairs retained through stochastic thinning. If 1 (default), no thinning is applied and all nearest-neighbor pairs are used. If <1, pairs are randomly retained using independent Bernoulli sampling.
maxtime: Numeric; an optional positive integer indicating the order of temporal neighborhood in the composite likelihood computation.
memdist: Logical; if TRUE then all distances useful in the composite likelihood estimation are computed before the optimization. FALSE is deprecated.
method: String; the type of matrix decomposition used in the likelihood computation. Default is "cholesky". Another possible choice is "svd" (when available).
model: String; the type of RF and therefore the densities associated to the likelihood objects. "Gaussian" is the default; see the Section Details.
n: Numeric; number of trials in a binomial RF; number of successes in a negative binomial RF.
onlyvar: Logical; if TRUE (and varest=TRUE) only the variance-covariance matrix is computed without optimizing. Default is FALSE.
optimizer: String; the optimization algorithm (see optim for details). "Nelder-Mead" is the default. Other possible choices are "nlm", "BFGS", "SANN", "L-BFGS-B", "nlminb", "bobyqa". In these last three cases upper and lower bounds can be passed by the user. In the one-dimensional case, optimize is used.
radius: Numeric; the radius of the sphere in the case of lon-lat coordinates. Default value is 1.
score: Logical; should score function be computed? Default is FALSE.
sensitivity: Logical; if TRUE then the sensitivity matrix is computed.
sparse: Logical; if TRUE then maximum likelihood is computed using sparse matrix algorithms (e.g., spam). It should be used with compactly supported covariance models. Default is FALSE.
start: An optional named list with initial values for parameters used by the numerical routines in the maximization procedure. Default is NULL (see Details).
thin_method: String; thinning scheme used when p_neighb<1. Default is "iid" (independent Bernoulli thinning).
type: String; the type of the likelihood objects. If "Pairwise" (default) then the marginal composite likelihood is formed by pairwise marginal likelihoods (see Details).
upper: An optional named list giving upper bounds for parameters when the optimizer is L-BFGS-B, nlminb, bobyqa or optimize. Names must match those in start.
varest: Logical; if TRUE the estimates' variances and standard errors are returned. For composite likelihood estimation it is deprecated. Use sensitivity=TRUE and update the object using GeoVarestbootstrap. Default is FALSE.
weighted: Logical; if TRUE the likelihood objects are weighted; see the Section Details. Default is FALSE.
X: Numeric; matrix of spatio(temporal) covariates in the linear mean specification.
nosym: Logical; if TRUE symmetric weights are not considered. This allows faster but less efficient CL estimation.
spobj: An object of class sp or spacetime.
spdata: Character; the name of data in the sp or spacetime object.

Author

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

Details

The function GeoFit2 is similar to GeoFit. However, GeoFit2 performs a preliminary estimation using maximum independence composite likelihood for the marginal parameters of the model and then uses the obtained estimates as starting values in the final weighted composite likelihood estimation (that includes both marginal and dependence parameters). This provides robust starting values for the marginal parameters in the optimization algorithm.

Stochastic thinning of nearest-neighbor pairs can be enabled via p_neighb<1. The argument thin_method controls the thinning scheme (default "iid").

Examples

Run this code



library(GeoModels)


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

################################################################
###
### Example 1 : Maximum pairwise conditional likelihood fitting 
###  of a Gaussian RF with Matern correlation
###
###############################################################
model="Gaussian"
# 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 <- "Matern"
mean <- 0.2
mean1 <- -0.5
sill <- 1
nugget <- 0
scale <- 0.2/3
smooth=0.5
param<-list(mean=mean,mean1=mean1,sill=sill,nugget=nugget,scale=scale,smooth=smooth)

# Simulation of the spatial Gaussian RF:
data <- GeoSim(coordx=coords,model=model,corrmodel=corrmodel, param=param,X=X)$data

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

################################################################
###
### Maximum pairwise likelihood fitting of
### Gaussian RFs with exponential correlation.
### 
###############################################################
fit1 <- GeoFit2(data=data,coordx=coords,corrmodel=corrmodel, 
                    optimizer="BFGS",neighb=3,likelihood="Conditional",
                    type="Pairwise", start=start,fixed=fixed,X=X)
print(fit1)




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


################################################################
###
### Example 2. Maximum pairwise likelihood fitting of 
### a LogGaussian  RF with Generalized 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
sill=0.5
scale=0.2
smooth=0

model="LogGaussian"
corrmodel="GenWend"
param=list(mean=mean,mean1=mean1,sill=sill,scale=scale,
                    nugget=nugget,power2=4,smooth=smooth)
# Simulation of a  non stationary LogGaussian RF:
data <- GeoSim(coordx=coords, corrmodel=corrmodel,model=model,X=X,
           param=param)$data

fixed<-list(nugget=nugget,power2=4,smooth=smooth)
start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill)
I=Inf
lower<-list(mean=-I,mean1=-I,scale=0,sill=0)
upper<-list(mean= I,mean1= I,scale=I,sill=I)

# Maximum pairwise composite-likelihood fitting of the RF:
fit <- GeoFit2(data=data,coordx=coords,corrmodel=corrmodel, model=model,
                    neighb=3,likelihood="Conditional",type="Pairwise",X=X,
                    optimizer="nlminb",lower=lower,upper=upper,
                    start=start,fixed=fixed)
print(unlist(fit$param))


################################################################
###
### Example 3. Maximum pairwise likelihood fitting of
### SinhAsinh   RFs with Wendland0 correlation
###
###############################################################
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 <- GeoFit2(data=data,coordx=coords,corrmodel=corrmodel, model=model,
                    neighb=3,likelihood="Marginal",type="Pairwise",
                    start=start,fixed=fixed)
print(unlist(fit1$param))

Run the code above in your browser using DataLab