| GeoFit | R Documentation |
Maximum weighted composite-likelihood fitting for Gaussian and some Non-Gaussian univariate spatial, spatio-temporal and bivariate spatial random fields. The function allows fixing any of the parameters and setting upper/lower bounds in the optimization. Different optimization methods can be used.
GeoFit(data, coordx, coordy=NULL, coordz=NULL, coordt=NULL,
coordx_dyn=NULL, copula=NULL, corrmodel=NULL, 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, spobj=NULL, spdata=NULL)
data |
A |
coordx |
A numeric ( |
coordy |
A numeric vector giving 1-dimension of spatial coordinates; optional argument, default is |
coordz |
A numeric vector giving 1-dimension of spatial coordinates; optional argument, default is |
coordt |
A numeric vector assigning 1-dimension of temporal coordinates. Optional argument, default is |
coordx_dyn |
A list of |
copula |
String; the type of copula. It can be |
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 |
fixed |
An optional named list giving the values of the parameters that will be considered as known values. The listed parameters for a given model/correlation function will not be estimated. |
anisopars |
A list of two elements: |
est.aniso |
A bivariate logical vector providing which anisotropy parameters must be estimated. |
grid |
Logical; if |
likelihood |
String; the configuration of the composite likelihood. |
lower |
An optional named list giving lower bounds for parameters when the optimizer is |
maxdist |
Numeric; an optional positive value indicating the maximum spatial distance considered in the composite 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 |
maxtime |
Numeric; an optional positive integer indicating the order of temporal neighborhood in the composite likelihood computation. |
memdist |
Logical; if |
method |
String; the type of matrix decomposition/linear algebra backend used in likelihood computations.
Default is |
model |
String; the type of random field (and associated density) used in the likelihood objects. Default is |
n |
Numeric; number of trials in a binomial random field; number of successes in a negative binomial random field. |
onlyvar |
Logical; if |
optimizer |
String; the optimization algorithm (see |
radius |
Numeric; the radius of the sphere in the case of lon-lat coordinates. Default is |
score |
Logical; if |
sensitivity |
Logical; if |
sparse |
Logical; if |
start |
An optional named list with initial values for parameters to be estimated. Default is |
thin_method |
String; thinning scheme used when |
type |
String; the type of likelihood objects. If |
upper |
An optional named list giving upper bounds for parameters when the optimizer is |
varest |
Logical; if |
weighted |
Logical; if |
X |
Numeric; matrix of spatio(temporal) covariates in the linear mean specification. |
spobj |
An object of class |
spdata |
Character; the name of the data component in the |
GeoFit provides weighted composite likelihood estimation based on pairs and independence composite likelihood estimation
for Gaussian and non-Gaussian random fields. Specifically, marginal and conditional pairwise likelihoods are available
for each type of random field.
The optimization method is specified using optimizer. The default method is Nelder-Mead; other available methods are
nlm, BFGS, SANN, L-BFGS-B, bobyqa, and nlminb. In the last three cases,
bounds can be specified using lower and upper.
Depending on the dimension of data and on the name of the correlation model,
the observations are assumed to be a realization of a spatial, spatio-temporal or bivariate random field.
Specifically, with data, coordx, coordy, coordt:
If data is a numeric d-dimensional vector and coordx, coordy are two numeric d-dimensional vectors
(or coordx is a (d \times 2)-matrix and coordy=NULL), then the data are interpreted as a single spatial realisation observed on d spatial sites;
If data is a numeric (t \times d)-matrix and coordt is a numeric t-dimensional vector, then the data are interpreted as a single spatio-temporal realisation observed on d sites and t times;
If data is a numeric (2 \times d)-matrix, then the data are interpreted as a single bivariate spatial realisation observed on d spatial sites;
If data is a list, coordx_dyn is a list and coordt is a numeric t-dimensional vector, then the data are interpreted as a spatio-temporal realisation observed on dynamical spatial sites (different locations for each time) and for t times.
It is also possible to specify a matrix of covariates using X. Specifically:
In the spatial case, X must be a (d \times k) matrix associated to data (a d-vector);
In the spatio-temporal case, X must be a (N \times k) matrix associated to data (a t \times d-matrix), where N=t\times d;
In the bivariate case, X must be a (N \times k) matrix associated to data (a 2 \times d-matrix), where N=2\times d.
The distance parameter allows different kinds of spatial distances:
Eucl, euclidean distance (default);
Chor, chordal distance;
Geod, geodesic distance.
The likelihood parameter represents the composite-likelihood configuration:
Conditional, composite likelihood formed by conditionals;
Marginal, composite likelihood formed by marginals (default);
Full, standard likelihood.
It must be coupled with type:
Pairwise, composite likelihood based on pairs;
Independence, composite likelihood based on independence;
Standard, standard likelihood.
Stochastic thinning of nearest-neighbor pairs can be enabled via p_neighb<1. The argument thin_method
controls the thinning scheme (default "iid").
Returns an object of class GeoFit.
An object of class GeoFit is a list containing at most the following components:
bivariate |
Logical: |
clic |
The composite information criterion; if the full likelihood is considered then it coincides with AIC. |
coordx |
A |
coordy |
A |
coordt |
A |
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/matrix/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 |
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 |
|
n |
The number of trials/successes for (negative) binomial models. |
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. |
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 likelihood objects. |
X |
The matrix of covariates. |
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
General 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.
Non Gaussian random fields:
Alegrıa A., Caro S., Bevilacqua M., Porcu E., Clarke J. (2017) Estimating covariance functions of multivariate skew-Gaussian random fields on the sphere. Spatial Statistics 22 388-402
Alegria A., Bevilacqua, M., Porcu, E. (2016) Likelihood-based inference for multivariate space-time wrapped-Gaussian fields. Journal of Statistical Computation and Simulation. 86(13), 2583–2597.
Bevilacqua M., Caamano C., Gaetan C. (2020) On modeling positive continuous data with spatio-temporal dependence. Environmetrics 31(7)
Bevilacqua M., Caamaño C., Arellano Valle R.B., Morales-Oñate V. (2020) Non-Gaussian Geostatistical Modeling using (skew) t Processes. Scandinavian Journal of Statistics.
Blasi F., Caamaño C., Bevilacqua M., Furrer R. (2022) A selective view of climatological data and likelihood estimation Spatial Statistics 10.1016/j.spasta.2022.100596
Bevilacqua M., Caamaño C., Arellano-Valle R. B., Camilo Gomez C. (2022) A class of random fields with two-piece marginal distributions for modeling point-referenced data with spatial outliers. Test 10.1007/s11749-021-00797-5
Morales-Navarrete D., Bevilacqua M., Caamaño C., Castro L.M. (2022) Modelling Point Referenced Spatial Count Data: A Poisson Process Approach TJournal of the American Statistical Association To appear
Caamaño C., Bevilacqua M., López C., Morales-Oñate V. (2023) Nearest neighbours weighted composite likelihood based on pairs for (non-)Gaussian massive spatial data with an application to Tukey-hh random fields estimation Computational Statistics and Data Analysis To appear
Bevilacqua M., Alvarado E., Caamaño C., (2023) A flexible Clayton-like spatial copula with application to bounded support data Journal of Multivariate Analysis To appear
Weighted Composite-likelihood for (non-)Gaussian random fields:
Bevilacqua, M. Gaetan, C., Mateu, J. and Porcu, E. (2012) Estimating space and space-time covariance 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. (2015) Comparing composite likelihood methods based on pairs for spatial Gaussian random fields. Statistics and Computing, 25(5), 877-892.
Caamaño C., Bevilacqua M., López C., Morales-Oñate V. (2023) Nearest neighbours weighted composite likelihood based on pairs for (non-)Gaussian massive spatial data with an application to Tukey-hh random fields estimation Computational Statistics and Data Analysis To appear
library(GeoModels)
###############################################################
############ Examples of spatial Gaussian random fieldss ################
###############################################################
# Define the spatial-coordinates of the points:
set.seed(3)
N=300 # number of location sites
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)
coords <- cbind(x,y)
# Define spatial matrix covariates and regression parameters
X=cbind(rep(1,N),runif(N))
mean <- 0.2
mean1 <- -0.5
# Set the covariance model's parameters:
corrmodel <- "Matern"
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 random fields:
data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param,X=X)$data
################################################################
###
### Example 0. Maximum independence composite likelihood fitting of
### a Gaussian random fields (no dependence parameters)
###
###############################################################
# setting starting parameters to be estimated
start<-list(mean=mean,mean1=mean1,sill=sill)
fit1 <- GeoFit(data=data,coordx=coords,likelihood="Marginal",
type="Independence", start=start,X=X)
print(fit1)
################################################################
###
### Example 1. Maximum conditional pairwise likelihood fitting of
### a Gaussian random fields using Nelder-Mead
###
###############################################################
# setting fixed and starting parameters to be estimated
fixed<-list(nugget=nugget,smooth=smooth)
start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill)
fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,
neighb=3,likelihood="Conditional",optimizer="Nelder-Mead",
type="Pairwise", start=start,fixed=fixed,X=X)
print(fit1)
################################################################
###
### Example 2. Standard Maximum likelihood fitting of
### a Gaussian random fields using nlminb
###
###############################################################
# Define the spatial-coordinates of the points:
set.seed(3)
N=250 # number of location sites
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)
coords <- cbind(x,y)
param<-list(mean=mean,sill=sill,nugget=nugget,scale=scale,smooth=smooth)
data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param)$data
# setting fixed and parameters to be estimated
fixed<-list(nugget=nugget,smooth=smooth)
start<-list(mean=mean,scale=scale,sill=sill)
I=Inf
lower<-list(mean=-I,scale=0,sill=0)
upper<-list(mean=I,scale=I,sill=I)
fit2 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,
optimizer="nlminb",upper=upper,lower=lower,
likelihood="Full",type="Standard",
start=start,fixed=fixed)
print(fit2)
###############################################################
############ Examples of spatial non-Gaussian random fieldss #############
###############################################################
################################################################
###
### Example 3. Maximum pairwise likelihood fitting of a Weibull random fields
### with Generalized Wendland correlation with Nelder-Mead
###
###############################################################
set.seed(524)
# 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
shape=2
scale=0.2
smooth=0
model="Weibull"
corrmodel="GenWend"
param=list(mean=mean,mean1=mean1,scale=scale,
shape=shape,nugget=nugget,power2=4,smooth=smooth)
# Simulation of a non stationary weibull random fields:
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,shape=shape)
# Maximum independence likelihood:
fit <- GeoFit(data=data, coordx=coords, X=X,
likelihood="Marginal",type="Independence", corrmodel=corrmodel,
,model=model, start=start, fixed=fixed)
print(unlist(fit$param))
## estimating dependence parameter fixing vector mean parameter
Xb=as.numeric(X %*% c(mean,mean1))
fixed<-list(nugget=nugget,power2=4,smooth=smooth,mean=Xb)
start<-list(scale=scale,shape=shape)
# Maximum conditional composite-likelihood fitting of the random fields:
fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,
neighb=3,likelihood="Conditional",type="Pairwise",
optimizer="Nelder-Mead",
start=start,fixed=fixed)
print(unlist(fit1$param))
### joint estimation of the dependence parameter and mean parameters
fixed<-list(nugget=nugget,power2=4,smooth=smooth)
start<-list(mean=mean,mean1=mean1,scale=scale,shape=shape)
fit2 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,
neighb=3,likelihood="Conditional",type="Pairwise",X=X,
optimizer="Nelder-Mead",
start=start,fixed=fixed)
print(unlist(fit2$param))
################################################################
###
### Example 4. Maximum pairwise likelihood fitting of
### a Skew-Gaussian spatial random fields with Wendland correlation
###
###############################################################
set.seed(261)
model="SkewGaussian"
# 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=-4.5
power2=4
c_supp=0.2
# model parameters
param=list(power2=power2,skew=skew,
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,mean=mean,sill=sill)
lower=list(scale=0,skew=-I,mean=-I,sill=0)
upper=list(scale=I,skew=I,mean=I,sill=I)
# Maximum marginal pairwise likelihood:
fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,
neighb=3,likelihood="Marginal",type="Pairwise",
optimizer="bobyqa",lower=lower,upper=upper,
start=start,fixed=fixed)
print(unlist(fit1$param))
################################################################
###
### Example 5. Maximum pairwise likelihood fitting of
### a Bernoulli random fields with exponential correlation
###
###############################################################
set.seed(422)
N=250
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,nugget=0,scale=0.2,power2=4)
# Simulation of the spatial Binomial-Gaussian random fields:
data <- GeoSim(coordx=coords, corrmodel=corrmodel, model="Binomial", n=1,X=X,
param=param)$data
## estimating the marginal parameters using independence cl
fixed <- list(power2=4,scale=0.2,nugget=0)
start <- list(mean=mean,mean1=mean1,mean2=mean2)
# Maximum independence likelihood:
fit <- GeoFit(data=data, coordx=coords,n=1, X=X,
likelihood="Marginal",type="Independence",corrmodel=corrmodel,
,model="Binomial", start=start, fixed=fixed)
print(fit)
## estimating dependence parameter fixing vector mean parameter
Xb=as.numeric(X %*% c(mean,mean1,mean2))
fixed <- list(nugget=0,power2=4,mean=Xb)
start <- list(scale=0.2)
# Maximum conditional pairwise likelihood:
fit1 <- GeoFit(data=data, coordx=coords, corrmodel=corrmodel,n=1,
likelihood="Conditional",type="Pairwise", neighb=3
,model="Binomial", start=start, fixed=fixed)
print(fit1)
## estimating jointly marginal and dependnce parameters
fixed <- list(nugget=0,power2=4)
start <- list(mean=mean,mean1=mean1,mean2=mean2,scale=0.2)
# Maximum conditional pairwise likelihood:
fit2 <- GeoFit(data=data, coordx=coords, corrmodel=corrmodel,n=1, X=X,
likelihood="Conditional",type="Pairwise", neighb=3
,model="Binomial", start=start, fixed=fixed)
print(fit2)
###############################################################
######### Examples of Gaussian spatio-temporal random fieldss ###########
###############################################################
set.seed(52)
# Define the temporal sequence:
time <- seq(1, 9, 1)
# Define the spatial-coordinates of the points:
x <- runif(20, 0, 1)
y <- runif(20, 0, 1)
coords=cbind(x,y)
# Set the covariance model's parameters:
scale_s=0.2/3;scale_t=1
smooth_s=0.5;smooth_t=0.5
sill=1
nugget=0
mean=0
param<-list(mean=0,scale_s=scale_s,scale_t=scale_t,
smooth_t=smooth_t, smooth_s=smooth_s ,sill=sill,nugget=nugget)
# Simulation of the spatial-temporal Gaussian random fields:
data <- GeoSim(coordx=coords,coordt=time,corrmodel="Matern_Matern",
param=param)$data
################################################################
###
### Example 6. Maximum pairwise likelihood fitting of a
### space time Gaussian random fields with double-exponential correlation
###
###############################################################
# Fixed parameters
fixed<-list(nugget=nugget,smooth_s=smooth_s,smooth_t=smooth_t)
# Starting value for the estimated parameters
start<-list(mean=mean,scale_s=scale_s,scale_t=scale_t,sill=sill)
# Maximum composite-likelihood fitting of the random fields:
fit <- GeoFit(data=data,coordx=coords,coordt=time,
corrmodel="Matern_Matern",maxtime=1,neighb=3,
likelihood="Marginal",type="Pairwise",
start=start,fixed=fixed)
print(fit)
###############################################################
######### Examples of a bivariate Gaussian random fields ###########
###############################################################
################################################################
### Example 7. Maximum pairwise likelihood fitting of a
### bivariate Gaussian random fields with separable Bivariate matern
### (cross) correlation model
###############################################################
# Define the spatial-coordinates of the points:
set.seed(89)
x <- runif(300, 0, 1)
y <- runif(300, 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,
nugget_1=0,nugget_2=0,pcol=0.2)
# Simulation of a spatial bivariate Gaussian random fields:
data <- GeoSim(coordx=coords, corrmodel="Bi_Matern_sep",
param=param)$data
# selecting fixed and estimated parameters
fixed=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,smooth=0.5)
start=list(sill_1=var(data[1,]),sill_2=var(data[2,]),
scale=0.1,pcol=cor(data[1,],data[2,]))
# Maximum marginal pairwise likelihood
fitcl<- GeoFit(data=data, coordx=coords, corrmodel="Bi_Matern_sep",
likelihood="Marginal",type="Pairwise",
start=start,fixed=fixed,
neighb=3)
print(fitcl)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.