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 = FALSE, ncores = 6, progress = TRUE, seed = NULL,
 min_success_rate = 0.8)

Arguments

fit	A fitted object obtained from GeoFit. The object must contain the sensitivity matrix, hence GeoFit should be called with sensitivity = TRUE. Full/Standard likelihood fits are not accepted; use GeoFit(..., likelihood = "Full", type = "Standard", varest = TRUE) for Hessian-based standard errors.
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; default FALSE. If TRUE, the score evaluation step is parallelized.
ncores	Positive integer or NULL; default 6. With parallel=TRUE, an explicit integer requests that many workers, capped by detected cores and the number of bootstrap jobs. Set ncores=NULL for automatic selection, capped at six workers and normally leaving one detected core free.
progress	Logical; if TRUE, progress information is shown.
seed	Optional integer seed for reproducibility of the simulated samples.
min_success_rate	Minimum fraction of successful score-bootstrap replications required to return variance estimates. The default is 0.8.

Details

For spatio-temporal fits, the original numeric coordt values stored in fit are reused. Irregularly spaced times are supported with method = "cholesky"; approximate simulation methods retain the restrictions documented in GeoSimapprox.

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.

For stochastic nearest-neighbor fits, the realized retained-pair graph stored in the original fit is reused unchanged for every simulated data set. Thus the score bootstrap estimates variability conditional on the selected pair graph, in agreement with the current implementation described in the stochastic-NN methodology.

Parallel multisession workers are started with the package-library paths of the calling R session, including the library containing the loaded GeoModels installation.

For method = "TB" without a copula, parallel workers are persistent and replications are streamed one at a time (simulate, evaluate, release). This keeps peak memory bounded when both L and K are large while still reusing the selected pair graph and static fitting context within each worker. For fitted copula models, method="TB" is simulated through the turning-bands backend of GeoSimCopula; the copula path is currently not streamed replicate-by-replicate. On platforms where future reports forked multicore execution as safe, GeoModels uses it automatically to reduce worker startup and serialization overhead; otherwise it falls back to multisession with the current GeoModels library path propagated explicitly.

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.

Three-dimensional coordinates

For a purely spatial univariate fit with explicit irregular three-dimensional Euclidean coordinates, method = "TB" is supported when no anisopars were used. Method "CE" remains restricted to regular two-dimensional grids. Use method = "cholesky" for bivariate, spatio-temporal, anisotropic, or other unsupported three-dimensional cases.

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 Sept. 23, 2026, 5:07 p.m.