GeoFit: Maximum-Likelihood-Based Fitting of Gaussian and non-Gaussian...

In GeoModels: Procedures for Gaussian and Non Gaussian Geostatistical (Large) Data Analysis

Maximum-Likelihood-Based Fitting of Gaussian and non-Gaussian random fields.

Description

Maximum weighted (stochastic) 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.

Usage

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="bernoulli", type="Pairwise", upper=NULL, 
 varest=FALSE, weighted=FALSE,
 X=NULL, spobj=NULL, spdata=NULL)

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 spatio-temporal realisation) or a (d \times d \times t \times n)-array (a single spatio-temporal realisation on regular grid). See Details for the accepted data layouts.
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 random field 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; see GeoCovmatrix for the list of implemented correlation models.
distance	String; the name of the spatial distance. Default is "Eucl" (Euclidean distance). See Details for the accepted options.
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: "angle" and "ratio", i.e. the anisotropy angle and the anisotropy ratio, respectively.
est.aniso	A bivariate logical vector providing which anisotropy parameters must be estimated.
grid	Logical; if FALSE (default) the data are interpreted as spatial or spatio-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 Details for the accepted options.
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 computation. See Details for more information.
neighb	Numeric; an optional positive integer indicating the order of neighborhood in the composite likelihood computation. See 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 (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/linear algebra backend used in likelihood computations. Default is "cholesky". Another possible choice is "svd" (when available).
model	String; the type of random field (and associated density) used in the likelihood objects. Default is "Gaussian"; see Details for the accepted options.
n	Numeric; number of trials in a binomial random field; number of successes in a negative binomial random field.
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). Default is "Nelder-Mead". Other possible choices are "nlm", "BFGS", "SANN", "L-BFGS-B", "nlminb", "bobyqa". For "L-BFGS-B", "nlminb" and "bobyqa" bounds can be passed via lower and upper. In the one-dimensional case, optimize is used.
radius	Numeric; the radius of the sphere in the case of lon-lat coordinates. Default is 1.
score	Logical; if TRUE the score function is computed. Default is FALSE.
sensitivity	Logical; if TRUE the sensitivity matrix is computed. Default is FALSE.
sparse	Logical; if TRUE then maximum likelihood / composite likelihood may exploit sparse-matrix algorithms (e.g., spam). Typically used with compactly supported covariance models. Default is FALSE.
start	An optional named list with initial values for parameters to be estimated. Default is NULL (see Details).
thin_method	String; thinning scheme in stochastic weighted pairwise likelihood (used when p_neighb < 1). Default is "bernoulli" (independent Bernoulli thinning). An alternative thinning scheme is "FixedBudget".
type	String; the type of likelihood objects. If "Pairwise" (default) the composite likelihood is formed by pairwise components (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 GeoVarest. Default is FALSE.
weighted	Logical; if TRUE the likelihood objects are weighted; see Details for the accepted options. Default is FALSE.
X	Numeric; matrix of spatio(temporal) covariates in the linear mean specification.
spobj	An object of class sp or spacetime.
spdata	Character; the name of the data component in the sp or spacetime object.

Details

GeoFit provides weighted and stochastic weighted composite-likelihood estimation based on pairs for Gaussian and non-Gaussian random fields, including nearest-neighbor and stochastic nearest-neighbor pairwise likelihoods; see Caamaño-Carrillo et al. (2024) and Bevilacqua et al. (2026). It also provides independence composite-likelihood estimation. 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:

1. Eucl, Euclidean distance (default);

2. Chor, chordal distance;

3. Geod, geodesic distance.

The likelihood parameter represents the composite-likelihood configuration:

1. Conditional, composite likelihood formed by conditionals;

2. Marginal, composite likelihood formed by marginals (default);

3. Full, standard likelihood.

It must be coupled with type:

1. Pairwise, composite likelihood based on pairs;

2. Independence, composite likelihood based on independence;

3. 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 "bernoulli").

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 random field is bivariate, otherwise FALSE.
clic	The composite information criterion after a GeoVarest call; 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/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 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/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.

Author(s)

Moreno Bevilacqua, [e-mail],https://sites.google.com/view/moreno-bevilacqua/home, Víctor Morales Oñate, [e-mail], https://sites.google.com/site/moralesonatevictor/, Christian Caamaño-Carrillo, [e-mail],https://www.researchgate.net/profile/Christian-Caamano

References

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. and Clarke, J. (2017). Estimating covariance functions of multivariate skew-Gaussian random fields on the sphere. Spatial Statistics, 22, 388–402. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.spasta.2017.05.003")}

Alegría, A., Bevilacqua, M. and 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., Caamaño-Carrillo, C. and Gaetan, C. (2020). On modelling positive continuous data with spatio-temporal dependence. Environmetrics, 31(7), e2628. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/env.2628")}

Bevilacqua, M., Caamaño-Carrillo, C., Arellano-Valle, R. B. and Morales-Oñate, V. (2021). Non-Gaussian geostatistical modeling using (skew) t processes. Scandinavian Journal of Statistics, 48(1), 212–245. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1111/sjos.12447")}

Blasi, F., Caamaño-Carrillo, C., Bevilacqua, M. and Furrer, R. (2022). A selective view of climatological data and likelihood estimation. Spatial Statistics, 50, 100596. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.spasta.2022.100596")}

Bevilacqua, M., Caamaño-Carrillo, C., Arellano-Valle, R. B. and Gómez, C. (2022). A class of random fields with two-piece marginal distributions for modeling point-referenced data with spatial outliers. TEST, 31(3), 644–674. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s11749-021-00797-5")}

Morales-Navarrete, D., Bevilacqua, M., Caamaño-Carrillo, C. and Castro, L. M. (2024). Modelling point referenced spatial count data: A Poisson process approach. Journal of the American Statistical Association, 119(545), 664–677. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.2022.2140053")}

Bevilacqua, M., Alvarado, E. and Caamaño-Carrillo, C. (2024). A flexible Clayton-like spatial copula with application to bounded support data. Journal of Multivariate Analysis, 201, 105277. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.jmva.2023.105277")}

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 and Methods, 107, 268–280. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.2011.646928")}

Bevilacqua, M. and Gaetan, C. (2015). Comparing composite likelihood methods based on pairs for spatial Gaussian random fields. Statistics and Computing, 25(5), 877–892. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s11222-014-9471-3")}

Caamaño-Carrillo, C., Bevilacqua, M., López, C. and Morales-Oñate, V. (2024). 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, 191, 107887. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.csda.2023.107887")}

Bevilacqua, M., Cuevas-Pacheco, F. and Caamaño-Carrillo, C. (2026). Fast stochastic nearest neighbor pairwise composite likelihood for massive spatial datasets. arXiv preprint, arXiv:2607.06142. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2607.06142")}

See Also

GeoCovmatrix for covariance matrix construction, GeoSim for simulation, GeoKrig for prediction, GeoVariogram and GeoWLS for variogram-based tools, GeoVarest for bootstrap variance estimation.

Examples

library(GeoModels)

###############################################################
############ Examples of spatial Gaussian random fields ################
###############################################################


# 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 field:
data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param,X=X)$data



################################################################
###
### Example 0. Maximum independence composite likelihood fitting of
### a Gaussian random field (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 field 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. Maximum stochastic marginal pairwise likelihood 
### fitting of a Gaussian random field using Nelder-Mead
### 
###############################################################

#N=100000 # number of location sites
#x <- runif(N, 0, 1)
#y <- runif(N, 0, 1)
#coords <- cbind(x,y)
#X=cbind(rep(1,N),runif(N))
#data <- GeoSimapprox(coordx=coords,method="TB",L=20000,parallel=TRUE,
#	  corrmodel=corrmodel, param=param,X=X)$data
#fixed<-list(nugget=nugget)
#start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill,smooth=smooth)

#fit2 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,
# neighb=3,likelihood="Marginal",optimizer="Nelder-Mead",
# p_neighb=0.2, thin_method="bernoulli",
# type="Pairwise", start=start,fixed=fixed,X=X)
#print(fit2)

################################################################
###
### Example 3. Standard Maximum likelihood fitting of
### a Gaussian random field 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 fields #############
###############################################################


################################################################
###
### Example 4. Maximum pairwise likelihood fitting of a Weibull random field 
### 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 field:
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 %*% unlist(fit$param)[1:2])
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 5. 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 6. Maximum pairwise likelihood fitting of 
### a Bernoulli random field 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)) # matrix 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 field:
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 %*% unlist(fit$param))
fixed <- list(nugget=0,power2=4,mean=Xb)
start <- list(scale=0.2)
lower <- list(scale=0)
upper <- list(scale=2)
# Maximum Marginal pairwise likelihood:
fit1 <- GeoFit(data=data, coordx=coords, corrmodel=corrmodel, n=1,
 likelihood="Marginal", type="Pairwise", neighb=3,
 model="Binomial", start=start, fixed=fixed,
 lower=list(scale=0.001), upper=list(scale=1))
 
print(fit1)


## estimating jointly marginal and dependence 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="Marginal", type="Pairwise", neighb=3,
 model="Binomial", start=start, fixed=fixed)
 
print(fit2)


###############################################################
######### Examples of Gaussian spatio-temporal random fields ###########
###############################################################
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 spatio-temporal Gaussian random field:
data <- GeoSim(coordx=coords,coordt=time,corrmodel="Matern_Matern",
 param=param)$data

################################################################
###
### Example 7. 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 field ###########
###############################################################

################################################################
### Example 8. 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 field:
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)

GeoModels documentation built on July 29, 2026, 5:06 p.m.