FitComposite {CompRandFld}  R Documentation 
Maximum weighted compositelikelihood fitting of spatiotemporal Gaussian, binary and spatial maxstable random fields. For the spatiotemporal Gaussian random field, (restricted) maximum likelihood and tapered likelihood fitting is also avalable. The function returns the model parameters' estimates and the estimates' variances and allows to fix any of the parameters.
FitComposite(data, coordx, coordy=NULL, coordt=NULL, corrmodel, distance='Eucl', fixed=NULL, grid=FALSE, likelihood='Marginal', margins='Gev', maxdist=NULL, maxtime=NULL, model='Gaussian', optimizer='NelderMead', replicates=1, start=NULL, taper=NULL, tapsep=NULL, threshold=NULL,type='Pairwise', varest=FALSE, vartype='SubSamp', weighted=FALSE, winconst, winstp)
data 
A ddimensional vector (a single spatial realisation) or a (n x d)matrix (n iid spatial realisations) or a (d x d)matrix (a single spatial realisation on regular grid) or an (d x d x n)array (n iid spatial realisations on regular grid) or a (t x d)matrix (a single spatialtemporal realisation) or an (t x d x n)array (n iid spatialtemporal realisations) or or an (d x d x t)array (a single spatialtemporal realisation on regular grid) or an (d x d x t x n)array (n iid spatialtemporal realisations on regular grid). For the description see the Section Details. 
coordx 
A numeric (d x 2)matrix (where

coordy 
A numeric vector assigning 1dimension of
spatial coordinates; 
coordt 
A numeric vector assigning 1dimension of
temporal coordinates. At the moment implemented only for the
Gaussian case. Optional argument, the default is 
corrmodel 
String; the name of a correlation model, for the description see the Section Details. 
distance 
String; the name of the spatial distance. The default is 
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 be not estimated, i.e. if

grid 
Logical; if 
likelihood 
String; the configuration of the composite
likelihood. 
margins 
String; the type of the marginal distribution of the
maxstable field. 
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 separation considered in the composite or tapered likelihood computation (see Details). 
model 
String; the type of random field and therefore the densities associated to the likelihood
objects. 
optimizer 
String; the optimization algorithm
(see 
replicates 
Numeric; a positive integer denoting the number of independent and identically distributed (iid) replications of a spatial or spatialtemporal random field. Optional argument, the default value is 1 then a single realisation is considered. 
start 
An optional named list with the initial values of the
parameters that are used by the numerical routines in maximization
procedure. 
taper 
String; the name of the type of covariance matrix.
It can be 
tapsep 
Numeric; an optional value indicating the separabe parameter in the space time quasi taper (see Details). 
threshold 
Numeric; a value indicating a threshold for the
binary random field. Optional in the case that 
type 
String; the type of the likelihood objects. If 
varest 
Logical; if 
vartype 
String; ( 
weighted 
Logical; if 
winconst 
Numeric; a positive value for computing the subwindow
size where observations are sampled in the subsampling procedure (if 
winstp 
Numeric; a value in (0,1] for computing the subwindow step (in the subsampling procedure). This value denote the proportion of the subwindow size. Optional argument, the default is 0.5. See Details for more information. 
Note, that the standard likelihood may be seen as particular case of the
composite likelihood. In this respect FitComposite
provides maximum (restricted) likelihood
fitting. Only composite likelihood estimation based on pairs are considered. Specifically marginal pairwise,
conditional pairwise and difference pairwise. Covariance tapering is considered only for Gaussian random fields.
With data
, coordx
, coordy
, coordt
, grid
and replicates
parameters:
If data
is a numeric ddimensional vector, coordx
and coordy
are two
numeric ddimensional vectors (or coordx
is (d x 2)matrix and coordy=NULL
),
coordt=NULL
, grid=FALSE
and replicates=1
, then the data are interpreted as a single spatial
realisation observed on d spatial sites;
If data
is a numeric (n x d)matrix, coordx
and coordy
are two
numeric ddimensional vectors (or coordx
is (d x 2)matrix and coordy=NULL
),
coordt=NULL
, grid=FALSE
and replicates=n
, then the data are interpreted as n
iid replications of a spatial
random field observed on d
spatial sites.
If data
is a numeric (d x d)matrix, coordx
and coordy
are two
numeric ddimensional vectors, coordt=NULL
, grid=TRUE
and replicates=1
, then the data are interpreted as
a single spatial random field realisation observed on d
equispaced spatial sites (named regular grid).
If data
is a numeric (d x d x n)array, coordx
and coordy
are two
numeric ddimensional vectors, coordt=NULL
, grid=TRUE
and replicates=n
, then the data are interpreted as
n
iid realisations of a spatial random field observed on d
equispaced spatial sites.
If data
is a numeric (t x d)matrix, coordx
and coordy
are two
numeric ddimensional vectors (or coordx
is (d x 2)matrix and coordy=NULL
),
coordt
is a numeric tdimensional vector, grid=FALSE
and replicates=1
, then the data
are interpreted as a single spatialtemporal realisation of a random field observed on d
spatial sites and for t
times.
If data
is a numeric (t x d x n)array, coordx
and coordy
are two
numeric ddimensional vectors (or coordx
is (d x 2)matrix and coordy=NULL
),
coordt
is a numeric tdimensional vector, grid=FALSE
and replicates=n
, then the data
are interpreted as n
iid realisations of a spatialtemporal random field observed on d
spatial sites and for t
times.
If data
is a numeric (d x d x t)array, coordx
and coordy
are two
numeric ddimensional vectors, coordt
is a numeric tdimensional vector, grid=TRUE
and
replicates=1
, then the data are interpreted as a single spatialtemporal realisation of a random field observed on
d
equispaced spatial sites and for t
times.
If data
is a numeric (d x d x t x n)array, coordx
and coordy
are two
numeric ddimensional vectors, coordt
is a numeric tdimensional vector, grid=TRUE
and
replicates=n
, then the data are interpreted as n
iid realisation of a spatialtemporal random field observed on
d
equispaced spatial sites and for t
times.
The corrmodel
parameter allows to select a specific correlation
function for the random field. (See Covmatrix
).
The distance
parameter allows to consider differents kinds of spatial distances.
The settings alternatives are:
Eucl
, the euclidean distance (default value);
Chor
, the chordal distance;
Geod
, the geodesic distance;
The likelihood
parameter represents the compositelikelihood
configurations. The settings alternatives are:
Conditional
, the compositelikelihood is formed by
conditionals likelihoods;
Marginal
, the compositelikelihood is formed by
marginals likelihoods;
Full
, the compositelikelihood turns out to be the standard likelihood;
The margins
parameter concerns only maxstable fields and indicates how the margins are
considered. The options are Gev
or Frechet
, where in the
former case the marginals are supposed generalized
extreme value distributed and in the latter case unit Frechet
distributed.
The maxdist
parameter set the maximum
spatial distance below which pairs of sites with inferior distances
are considered in the compositelikelihood. This can be
inferior of the effective maximum spatial distance. Note that
this corresponds to use a weighted compositelikelihood 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.
The same arguments of maxdist
are valid for maxtime
but
here the weigthed compositelikelihood regards the case of
spatialtemporal field. At the moment is
implemented only for Gaussian random fields. 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 model (see Covmatrix
). If they are not specified then the maximum spatial and temporal distances are considered.
In the case of space time quasi taper the tapsep
parameter allows to specify the spatio temporal compact support (see Covmatrix
).
The model
paramter indicates the type of random field
considered, for instance model=Gaussian
denotes a Gaussian random field.
Accordingly, this also determines the analytical expression of the finite dimensional distribution associated with the random field.
The available options are:
Gaussian
, for a Gaussian random field (see
i.e. Wackernagel, H. 1998);
BinaryGauss
, for a Binary random field (see Heagerty
and Lele 1998)
BrowResn
, for a BrownResnick maxstable random field (see Kabluchko, Z. et al. 2009);
ExtGauss
, for an Extremal Gaussian maxstable random
field (known also as Schlather model) (see Schlather, M. 2002);
ExtT
, for an Extremal t
maxstable random field (see Davison, A. C. et al. 2012);
Note, that only for the Gaussian
case the estimation procedure is implemented for spatial and spatialtemporal
random fields.
The start
parameter allows to specify starting values.
If start
is omitted the routine is computing the
starting values using the weighted moment estimator.
The taper
parameter, optional in case that
type=Tapering
, indicates the type of taper
correlation model. (See Covmatrix
)
The threshold
parameter indicates the value (common for all
the spatial sites) above which the values of the underlying Gaussian latent process
are considered sucesses events (values below are instead
failures). See e.g. Heagerty and Lele (1998) for more details.
The type
parameter represents the type of likelihood used in the
compositelikelihood definition. The possible alternatives are listed
in the following scheme.
If a Gaussian random field is considered
(model=Gaussian
):
If the composite is formed by marginal likelihoods (likelihood=Marginal
):
Pairwise
, the compositelikelihood is defined by
the pairwise likelihoods;
Difference
, the compositelikelihood is defined by
likelihoods
which are obtained as difference of the pairwise likelihoods.
If the composite is formed by conditional likelihoods
(likelihood=Conditional
)
Pairwise
, the compositelikelihood is defined 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
version of the full likelihood (e.g. Kaufman et al. 2008, see References).
The varest
parameter specifies if the standard error estimation of the estimated parameters must be computed.
For Gaussian random field and standard (restricted) likelihood estimation, standard errrors 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 random field and tapered and composite likelihood estimation, standard errors estimate are computed as square root of the diagonal 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
parameter specifies the method used to compute the estimates' variances in the composite likelihood case. In particular
for estimating the variability matrix J
in the Godambe expression matrix.
This parameter is considered if varest=TRUE
. The options are:
SubSamp
(the default), indicates the SubSampling method;
Sampling
, indicates that the variability matrix
is estimated by the sample contropart (available only for n iid
replications of the random field, i.e. replicates=n
);
The weighted
parameter specifies if the likelihoods forming the
compositelikelihood must be weighted. If TRUE
the weights are
selected by opportune procedures that improve the efficient of the
maximum compositelikelihood estimator (not implemented yet). If
FALSE
the efficient improvement procedure is not used.
For computing the standard errors by the subsampling procedure,
winconst
and winstp
parameters represent respectively a positive constant used to
determine the subwindow size and the the step with which the
subwindow moves.
In the spatial case (subset of R^2), the domain is seen as
a rectangle BxH, therefore the size of the
subwindow side b is given by b=winconst x sqrt(B) (similar is of h).
For a complete description see Lee and Lahiri (2002).
By default winconst
is set B / (2
x sqrt(B)).
The winstp
parameter is used to determine the subwindow step. The latter is given by the
proportion of the subwindow size, so that when winstp=1
there
is not overlapping between contiguous subwindows.
In the spatial case by default winstp=0.5
.
The subwindow 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 spatiotemporal case the subsampling is meant only in time as
described by Li et al. (2007). Thus, winconst
represents
the lenght of the temporal subwindow. By default the size of the
subwindow is computed following the rule established in Li et al. (2007).
By default winstp
is the time step.
Observe that in the spatiotemporal case, the returned values by
srange
and trange
, represent respectively the minimum and maximum
of the marginal spatial distances and those of the temporal
separations. Thus, the minimum being not the overall (i.e. considering
the spatiotemporal coordinates) is not zero, as one could be expect and
the latter can be easily added by the user.
Returns an object of class FitComposite
.
An object of class FitComposite
is a list containing
at most the following components:
clic 
The composite information criterion, if the full likelihood is considered then it coincides with the Akaike information criterion; 
coordx 
A ddimensional vector of spatial coordinates; 
coordy 
A ddimensional vector of spatial coordinates; 
coordt 
A tdimensional vector of temporal 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 fixed parameters; 
iterations 
The number of iteration used by the numerical routine; 
likelihood 
The configuration of the composite likelihood; 
logCompLik 
The value of the log compositelikelihood at the maximum; 
message 
Extra message passed from the numerical routines; 
model 
The density associated to the likelihood objects; 
nozero 
In the case of tapered likelihood the percentage of non zero values in the covariance matrix. Otherwise is NULL. 
numcoord 
The number of spatial coordinates; 
numrep 
The number of the iid replicatations of the random field; 
numtime 
The number the temporal realisations of the random field; 
param 
The vector of parameters' estimates; 
srange 
The minimum and maximum spatial distance (see Details). The maximum
is 
stderr 
The vector of standard errors; 
sensmat 
The sensitivity matrix; 
varcov 
The matrix of the variancecovariance of the estimates; 
varimat 
The variability matrix; 
vartype 
The method used to compute the variance of the estimates; 
trange 
The minimum and maximum temporal separation (see Details). The maximum
is 
threshold 
The threshold used in the binary random field. 
type 
The type of the likelihood objects. 
winconst 
The constant use to compute the window size in the subsampling procedure; 
winstp 
The step used for moving the window in the subsampling procedure 
Simone Padoan, simone.padoan@unibocconi.it, http://faculty.unibocconi.it/simonepadoan; Moreno Bevilacqua, moreno.bevilacqua@uv.cl, https://sites.google.com/a/uv.cl/morenobevilacqua/home.
Padoan, S. A. and Bevilacqua, M. (2015). Analysis of Random Fields Using CompRandFld. Journal of Statistical Software, 63(9), 1–27.
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 LikelihoodBased 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., 212, 433–452.
Compositelikelihood:
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 Compositelikelihood for binary random fields:
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 Compositelikelihood for maxstable random fields:
Davison, A. C. and Gholamrezaee, M. M. (2012) Geostatistics of extremes. Proceedings of the Royal Society of London, series A, 468, 581–608.
Padoan, S. A. (2008). Computational Methods for Complex Problems in Extreme Value Theory. PhD. Thesis, Department of Statistics, University of Padua.
Padoan, S. A. Ribatet, M. and Sisson, S. A. (2010) LikelihoodBased Inference for MaxStable Processes. Journal of the American Statistical Association, Theory & Methods, 105, 263–277.
Weighted Compositelikelihood for Gaussian random fields:
Bevilacqua, M. Gaetan, C., Mateu, J. and Porcu, E. (2012) Estimating space and spacetime 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. (2014) Comparing composite likelihood methods based on pairs for spatial Gaussian random fields Statistics and Computing.DOI:10.1111/2041210X.12167
Spatial Extremes:
Davison, A. C., Padoan, S. A., and Ribatet, M. (2012) Statistical Modelling of Spatial Extremes, with discussion. Statistical Science, 27, 161–186.
de Haan, L., and Pereira, T. T. (2006) Spatial Extremes: Models for the Stationary Case. The Annals of Statistics, 34, 146–168.
Kabluchko, Z. (2010) Extremes of Independent Gaussian Processes. Extremes, 14, 285–310.
Kabluchko, Z., Schlather, M., and de Haan, L. (2009) Stationary maxstable fields associated to negative definite functions. The Annals of Probability, 37, 2042–2065.
Schlather, M. (2002) Models for Stationary MaxStable Random Fields. Extremes, 5, 33–44.
Smith, R. L. (1990) MaxStable Processes and Spatial Extremes. Unpublished manuscript, University of North California.
Subsampling estimation:
Carlstein, E. (1986) The Use of Subseries Values for Estimating the Variance. The Annals of Statistics, 14, 1171–1179.
Heagerty, P. J. and Lumley T. (2000) Window Subsampling of Estimating Functions with Application 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 spacetime covariance functions. Journal of the American Statistical Association, 102, 736–744
Covmatrix
, WLeastSquare
, optim
library(CompRandFld) library(RandomFields) library(spam) set.seed(3132) ############################################################### ############ Examples of spatial random fields ################ ############################################################### # Define the spatialcoordinates of the points: x < runif(100, 0, 10) y < runif(100, 0, 10) # Set the covariance model's parameters: corrmodel < "exponential" mean < 0 sill < 1 nugget < 0 scale < 1.5 param<list(mean=mean,sill=sill,nugget=nugget,scale=scale) coords<cbind(x,y) # Simulation of the spatial Gaussian random field: data < RFsim(coordx=coords, corrmodel=corrmodel, param=param)$data # Fixed parameters fixed<list(mean=mean,nugget=nugget) # Starting value for the estimated parameters start<list(scale=scale,sill=sill) ################################################################ ### ### Example 1. Maximum likelihood fitting of ### Gaussian random fields with exponential correlation. ### One spatial replication. ### Likelihood setting: composite with ### marginal pairwise likelihood objects. ### ############################################################### # Maximum compositelikelihood fitting of the random field: fit < FitComposite(data, coordx=coords, corrmodel=corrmodel, maxdist=2, likelihood="Marginal",type="Pairwise",varest=TRUE, start=start,fixed=fixed) # Results: print(fit) ################################################################ ### ### Example 2. Maximum likelihood fitting of ### Gaussian random fields with exponential correlation. ### One spatial replication. ### Likelihood setting: standard full likelihood. ### ############################################################### # Maximum compositelikelihood fitting of the random field: fit < FitComposite(data, coordx=coords, corrmodel=corrmodel,likelihood="Full", type="Standard",varest=TRUE,start=start,fixed=fixed) # Results: print(fit) ################################################################ ### ### Example 3. Maximum likelihood fitting of ### Gaussian random fields with exponetial correlation. ### One spatial replication. ### Likelihood setting: tapered full likelihood. ### ############################################################### # Maximum tapered likelihood fitting of the random field: fit < FitComposite(data, coordx=coords, corrmodel=corrmodel,likelihood="Full", type="Tapering",taper="Wendland1",maxdist=1.5, varest=TRUE,start=start,fixed=fixed) # Results: print(fit) ################################################################ ### ### Example 4. Maximum compositelikelihood fitting of ### maxstable random fields. Extremal Gaussian model with ### exponential correlation. n iid spatial replications. ### Likelihood setting: composite with marginal pairwise ### likelihood objects. ### ############################################################### # Simulation of a maxstable random field in the specified points: data < RFsim(x, y, corrmodel=corrmodel, model="ExtGauss", replicates=30, param=list(mean=mean,sill=sill,nugget=nugget,scale=scale))$data # Maximum compositelikelihood fitting of the random field: fit < FitComposite(data, x, y, corrmodel=corrmodel, model="ExtGauss", replicates=30, varest=TRUE, vartype="Sampling", margins="Frechet",start=list(sill=sill,scale=scale)) # Results: print(fit) ################################################################ ### ### Example 5. Maximum likelihood fitting of ### BinaryGaussian random fields with exponential correlation. ### One spatial replication. ### Likelihood setting: composite with marginal pairwise ### likelihood objects. ### ############################################################### #set.seed(3128) #x < runif(200, 0, 10) #y < runif(200, 0, 10) # Simulation of the spatial BinaryGaussian random field: #data < RFsim(coordx=coords, corrmodel=corrmodel, model="BinaryGauss", # threshold=0, param=list(mean=mean,sill=.8, # nugget=nugget,scale=scale))$data # Maximum compositelikelihood fitting of the random field: #fit < FitComposite(data, coordx=coords, corrmodel=corrmodel, threshold=0, # model="BinaryGauss", fixed=list(nugget=nugget, # mean=0),start=list(scale=.1,sill=.1)) # Results: #print(fit) ############################################################### ######### Examples of spatiotemporal random fields ########### ############################################################### # Define the temporal sequence: #time < seq(1, 80, 1) # Define the spatialcoordinates of the points: #x < runif(10, 0, 10) #y < runif(10, 0, 10) #coords=cbind(x,y) # Set the covariance model's parameters: #corrmodel="exp_exp" #scale_s=0.5 #scale_t=1 #sill=1 #nugget=0 #mean=0 #param<list(mean=0,scale_s=1,scale_t=1,sill=sill,nugget=nugget) # Simulation of the spatialtemporal Gaussian random field: #data < RFsim(coordx=coords,coordt=time,corrmodel=corrmodel, # param=param)$data # Fixed parameters #fixed<list(mean=mean,nugget=nugget) # Starting value for the estimated parameters #start<list(scale_s=scale_s,scale_t=scale_t,sill=sill) ################################################################ ### ### Example 6. Maximum likelihood fitting of ### Gaussian random field with doubleexponential correlation. ### One spatiotemporal replication. ### Likelihood setting: composite with conditional pairwise ### likelihood objects. ### ############################################################### # Maximum compositelikelihood fitting of the random field: #fit < FitComposite(data=data,coordx=coords,coordt=time,corrmodel="exp_exp", # maxtime=2,maxdist=1,likelihood="Marginal",type="Pairwise", # start=start,fixed=fixed) # Results: #print(fit)