GeoVarest: Score-based variance estimation for 'GeoFit' objects

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

Score-based variance estimation for GeoFit objects

Description

The function updates a fitted GeoFit object by estimating the variability matrix of the composite likelihood score through parametric simulation. The fitted model is used to generate K independent datasets. For each simulated dataset, the composite likelihood score is evaluated at the original estimate \hat\theta, without refitting the model. The empirical variance of these simulated scores provides an estimate of the variability matrix J. Together with the sensitivity matrix H, computed by GeoFit when sensitivity = TRUE, this yields the Godambe sandwich covariance matrix

G^{-1} = H^{-1} J H^{-1}.

The updated object contains standard errors, Wald confidence intervals, p-values, the estimated matrices J, H^{-1} and G^{-1}, and composite likelihood information criteria based on the penalty \mathrm{tr}(H^{-1}J).

Usage

GeoVarest(fit, K = 100, sparse = FALSE,
 method = c("cholesky", "TB", "CE"),
 alpha = 0.95, L = 10000,
 parallel = TRUE, ncores = NULL, progress = TRUE, seed = NULL)

Arguments

fit	A fitted object obtained from GeoFit. The object must contain the sensitivity matrix, hence GeoFit should be called with sensitivity = TRUE.
K	The number of simulations used in the parametric score bootstrap.
sparse	Logical; if TRUE, then Cholesky decomposition is performed using sparse matrix algorithms.
method	String; the method of simulation. The default is "cholesky". For large data sets the options "TB" and "CE" call approximate simulation methods; see GeoSimapprox.
alpha	Numeric; the level of the confidence interval.
L	Numeric; the number of lines in the turning bands method.
parallel	Logical; if TRUE, the score evaluation step is parallelized.
ncores	Numeric; number of cores involved in the parallel score evaluation. If NULL, the number of cores is chosen automatically as 4.
progress	Logical; if TRUE, progress information is shown.
seed	Optional integer seed for reproducibility of the simulated samples.

Details

Let cl(\theta) denote the composite log-likelihood and let U(\theta) = \nabla cl(\theta) be the corresponding composite score. The function simulates K data sets from the fitted model and evaluates the composite score at the fitted parameter value \hat\theta. The variability matrix is estimated as

\hat J = Var\{U_1(\hat\theta), \ldots, U_K(\hat\theta)\}.

If H is the sensitivity matrix stored in fit$sensmat, the inverse Godambe matrix is estimated by

\widehat{G^{-1}} = H^{-1} \hat J H^{-1}.

Standard errors are obtained from the square root of the diagonal of \widehat{G^{-1}}.

For composite likelihoods, the penalty used in the information criterion is

tr(H^{-1}\hat J),

and the composite likelihood information criterion is computed as

-2 cl(\hat\theta) + 2 tr(H^{-1}\hat J).

Differently from GeoVarestbootstrap, this function does not refit the model for each simulated data set. It estimates the variability matrix of the score and then computes the sandwich/Godambe covariance matrix.

Value

Returns an updated object of class GeoFit. The following components are added or updated:

stderr	Estimated standard errors obtained from the inverse Godambe matrix.
varcov	Estimated inverse Godambe matrix \widehat{G^{-1}}.
godambe	Estimated Godambe matrix.
Jmat	Estimated variability matrix of the composite score.
Hinv	Inverse, or generalized inverse, of the sensitivity matrix.
claic	Composite likelihood AIC-type criterion.
clic	Same value as claic.
clbic	Composite likelihood BIC-type criterion.
clic_penalty	Penalty term tr(H^{-1}\hat J).
conf.int	Wald-type confidence intervals based on the estimated standard errors.
pvalues	Wald-type p-values.
scores	Matrix of successful bootstrap score evaluations.
score_logCompLik	Composite log-likelihood values corresponding to the successful score evaluations.
score_failures	Data frame with failed score evaluations, if any.
estimates	Matrix containing the successful scores and the associated composite log-likelihood values.

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

See Also

GeoFit for the fitted objects used as input

Examples

library(GeoModels)

################################################################
###
### Example 1. Test on the parameter
### of a regression model using conditional composite likelihood
###
###############################################################
set.seed(342)
model="Gaussian" 
# Define the spatial-coordinates of the points:
NN=3500
x = runif(NN, 0, 1)
y = runif(NN, 0, 1)
coords = cbind(x,y)
# Parameters
mean=1; mean1=-1.25; # regression parameters
 sill=1 # variance

# matrix covariates
X=cbind(rep(1,nrow(coords)),runif(nrow(coords)))

# model correlation 
corrmodel="Matern"
smooth=0.5;scale=0.1; nugget=0;

# simulation
param=list(smooth=smooth,mean=mean,mean1=mean1,
 sill=sill,scale=scale,nugget=nugget)
data = GeoSim(coordx=coords, corrmodel=corrmodel,
 model=model, param=param,X=X)$data

I=Inf

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

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 = GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,
 likelihood="Conditional",type="Pairwise",sensitivity=TRUE,
 lower=lower,upper=upper,neighb=3,
 optimizer="nlminb",X=X,
 start=start,fixed=fixed)

unlist(fit$param)


#fit_update=GeoVarest(fit,K=100,parallel=TRUE)
#fit_update$stderr
#fit_update$conf.int
#fit_update$pvalues

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