engine, taking an engine object from the new constructors inla_engine() and
jax_engine() (or the strings "inla", the default, and "jax" for their default settings).
Each constructor holds all settings of its engine, so hb_area() itself only describes the
model:
inla_engine(strategy, link, scale_model, prior_prec, prior_phi, prior_rho, prior_prec_time, prior_rho_time, prior_prec_subarea, draws, seed, ...), where ... goes to INLA::inla().
control.* lists given there are merged with the package defaults (before, they clashed with
them in hb_area() and replaced them in hb_unit()/hb_twofold()). Priors and draws left at NULL
take the default of the model function, so the defaults of each function are unchanged;
settings a function does not use are ignored with a warning.jax_engine(warmup, samples, chains, device, seed, prior_beta, robust, keep_draws, progress)
fits the same area-level models by full MCMC (NumPyro NUTS on JAX: CPU, NVIDIA GPU or the
Apple silicon GPU via device = "metal") through
'reticulate' (Suggests); the engine formerly shipped in the separate 'fastsaegpu' package now
lives in fastsae. It supports all six families, bym2/bym/besag/generic1 (alias
"leroux"), rw1/ar1/iid and every st_interaction, plus prior_beta = "horseshoe" and
robust = TRUE (Student-t area effects). Results have the same structure as the INLA engine,
with the posterior draws in draws, the largest split R-hat in rhat_max (convergence is
rhat_max < 1.05), and DIC/WAIC in goodness. Posterior self-benchmarking (also by
benchmark_group) adjusts the MCMC draws exactly as for INLA, and benchmark_sae() uses the
draws for the MSE matrix.engine_spec; engine holds the engine name.hb_unit() and hb_twofold() take the same engine = inla_engine(...) (INLA only). In
hb_twofold(), prior_prec_area is now inla_engine(prior_prec = ), n_samples is
inla_engine(draws = ), and the posterior draws for the area aggregates honour seed. Their
results carry engine and engine_spec too.jax_engine(device = "metal") now samples on the Apple GPU. jax-metal cannot compile two
operations of NumPyro's NUTS sampler (the popcnt bit count, and lax.cond statements whose
outputs are constants), so fits used to fall back to the CPU. fastsae now replaces both with
equivalent operations that Metal can compile, and steps the chains from Python; the draws match
the CPU sampler. On an Apple M4, Metal is faster than the CPU for large models (about 2x at
10,000 domains) and slower for small ones, where a fixed cost of about two minutes per fit
dominates (see the "hb_area() engines" benchmark article). device = "auto" keeps using the CPU
on macOS.setup_numpyro_env() and check_numpyro_available() to install and check the Python side.
Nothing is installed or downloaded unless setup_numpyro_env() is called. The check looks only at
an initialised Python, RETICULATE_PYTHON / RETICULATE_PYTHON_ENV and the r-numpyro-gpu
environment, so it answers in about a second instead of searching the whole system for Python
installations (which took over 20 s on Windows).method = "laplace" for non-spatial binomial and Poisson
models) has been removed; such models are now fitted by INLA. 'lme4' moved to Suggests.method (method = "laplace" becomes engine = inla_engine(strategy = "laplace")),
and passing strategy, link, scale_model, prior_*, benchmark_draws, benchmark_seed,
n_samples or INLA::inla() arguments directly to hb_area(), hb_unit() or hb_twofold().
They still work, with a warning, and are forwarded to inla_engine().print()/summary() show the engine, and the self-benchmarking line now names the method used
(posterior draws or Wang-Fuller-Qu) instead of always reporting Wang-Fuller-Qu.eblup_stfh() no longer requires a complete panel. Cells whose response is NA -- missing
periods of a sampled domain, or domains that are never sampled -- stay in the data and are
predicted. The model is defined on the full domain x time grid (the SAR effect on all domains
of W, the AR(1) effect on all periods); the variance components are estimated by REML from
the sampled cells (the marginal likelihood with the missing cells integrated out), and every
cell gets the best linear unbiased predictor given the sampled data (Saei and Chambers, 2005;
Benedetti et al., 2024). A missing period borrows from the other periods of its domain and
from the neighbours; a domain never sampled is predicted from its neighbours (spatial kriging).df_eblup gains a sampled column; print() reports the number of
predicted cells. The auxiliary variables must be known in every cell, and W must include
the unsampled domains.benchmarks/sim_stfh_unsampled.R checks the
MSE estimators by Monte Carlo.mys_panel and mys_proxmat (domains 21 and 25 are unsampled).mse_type, which includes "none" (no MSE; mse and
rse are NA): eblup_fh() ("analytical", "none"), eblup_sfh() ("analytical",
"pbmse", "npbmse", "none"), eblup_stfh() and eblup_twofold() ("analytical",
"bootstrap", "none") and eblup_bhf() ("bootstrap", "none"). The chosen value is
returned in mse_type.
eblup_sfh(mse_method = ) and compute_mse = of eblup_stfh() and eblup_bhf() are
deprecated and still work with a warning (compute_mse = FALSE is mse_type = "none").eblup_stfh() gives the analytical MSE and eblup_bhf()
the bootstrap MSE (both gave none before). Use mse_type = "none" for the previous output.eblup_fh(), eblup_sfh(),
eblup_stfh() and eblup_twofold() only printed a message and eblup_bhf() said nothing.
INLA fits of hb_area(), hb_unit() and hb_twofold() now read convergence from the INLA
mode search (ok, mode$mode.status; it was always TRUE), and hb_area(engine = "jax")
warns when the largest split R-hat is 1.05 or more.vardir must be finite and positive in every sampled area (a missing value gave an error
from the compiled code);NA), and at least two
sampled areas in eblup_twofold() and eblup_bhf();W of eblup_sfh() and eblup_stfh() must be non-negative; a zero row (a domain without
neighbours) or rows that do not sum to 1 give a warning, since the SAR model with
rho in (-1, 1) assumes a row-standardised W, as 'sae' requires;hb_area(),
hb_twofold() and in the Xpop means of hb_unit(): INLA treats it as 0, so the estimate
of that area was wrong without any message. hb_twofold() replaced a missing, zero or
negative vardir of a sampled subarea by 1, and hb_area() stopped with an INLA error on a
missing one; both now stop with the message of the EBLUP functions. hb_area() also rejects
linearly dependent covariates;eblup_bhf() names the domain (it gave an internal error on
the number of rows);domain values give a warning in eblup_fh(), eblup_sfh() and hb_area()
without time (in hb_area() such rows share one area effect).seed of eblup_sfh(), eblup_stfh(), eblup_bhf() and eblup_twofold() no longer
changes the random number stream of the session; the bootstrap MSEs are unchanged.eblup_twofold() uses B = 100 and seed = -1 (no seed) by default, as the other EBLUP
functions (it used B = 200 and seed = NULL; NULL still means no seed).hb_area(), hb_unit() and hb_twofold() are fitted in INLA's "classic"
mode (Rue, Martino and Chopin, 2009); the other families keep INLA's default "compact" mode
(van Niekerk et al., 2023), which leaves the linear predictors out of the latent field and is
designed for data-rich models. The choice rests on comparisons with the exact posterior
(numerical integration over the variances) and with long Stan runs with the same model and
priors (hb_area(): 12 configurations, 18 paired comparisons):
hb_area(), when the
area variance is small relative to vardir; it was accurate on mys), 8-15% (hb_unit(),
cornsoybean) and 10% (hb_twofold()), and hb_unit() gave different county means on
identical calls (up to 3 hectares and 13% in the SDs). In the classic mode the posterior means
were within 0.07 posterior SD and the SDs within 0.91-1.03 times the reference, at the same
speed.hb_unit() binomial with 4000
units: 72 s against 2.6 s).random_effect_var) is
biased downwards while the area estimates are accurate; this is documented.
inla_engine(inla.mode = ) selects a mode.spatial = "besag", "bym", "bym2", "leroux") of hb_area(), hb_unit() and
hb_twofold() use only the pattern W > 0. A dense W such as distance-decay weights becomes
a (nearly) complete graph and an almost unstructured spatial effect; it now gives a warning.
mys_proxmat is such a matrix (all 861 pairs linked): its documentation said it was based on
contiguous boundaries, and the examples and vignettes now fit CAR models on a 4-nearest-
neighbour graph of the sim_area centroids. mys_proxmat remains suited to the SAR models.hb_unit(family = "binomial" / "poisson") warns that the domain means are approximate. The
non-sampled part is predicted at the mean covariates, inv_link(xbar' beta + u_d), which is
biased for the mean of inv_link(x' beta + u_d) over the units: in a simulation -8%
(binomial) and -16% (Poisson) for the sampled domains, with 90% and 78% coverage of the 95%
intervals. Unbiased estimation needs unit-level population covariates (Molina and Rao, 2010);
see ?hb_unit. The Gaussian model is not affected.hb_twofold(weight = ) warns about missing or negative weights (they get weight 0 in
df_area, as before, but silently).compare_sae() between a Hierarchical Bayes and an EBLUP model no longer calls the model with
the smaller posterior variance "more efficient": it reports which one has the lower reported
uncertainty and notes that a posterior variance and an MSE are not the same quantity.inla_engine(seed = 1L) is the new
default (it was 0L, INLA's random seed), and the seed is also set on the R side, since INLA
draws the hyperparameter configurations with R's RNG; the session's random number stream is
restored. This affects the area aggregates of hb_twofold() (their SD changed by up to 9%
between identical calls), posterior self-benchmarking in hb_area() and benchmark_sae().eblup_twofold(),
hb_twofold()); the unit-level vignette covers hb_unit().hb_area(st_interaction = "type2" / "type3" / "type4") now fit with INLA. They always failed
before (the interaction reused the index of a main effect, which INLA rejects), and with
spatial = "none" the type III/IV interaction was silently dropped. The interaction is now a
separate term on copied indices: one temporal process per domain (type II), or an ICAR field
per period (type III) or evolving with temporal (type IV). Types III and IV need a
structured spatial effect and stop with a clear message otherwise.hb_area(spatial = "leroux", st_interaction = "separable") fits instead of crashing inside
INLA: the grouped Leroux field uses INLA's besagproper2, the same model as generic1.inla_engine(link = ) is used. It was accepted but ignored; it now sets the link of the
likelihood ("probit", "cloglog" and "loglog" for binomial and beta), checked against the
family. The posterior draws for self-benchmarking and benchmark_sae() use the matching inverse
link, and the link is returned in link.hb_area(spatial = "bym"): random_effect_var is the variance of the spatial component (it was
the iid one). Interaction terms no longer feed random_effect_var, rho or rho_time, and with
st_interaction = "separable" the AR(1) group correlation is returned in rho_time.hb_area() objects have class c("fastsae_hb_area", "fastsae_hb", "fastsae"), like
hb_unit() and hb_twofold().inla_engine(prior_prec = ) explains the different defaults: the PC prior
pc.prec(1, 0.01) of hb_area() and hb_twofold() assumes an area standard deviation of
order 1 (proportions, rates, log scales), while hb_unit() keeps the vague
loggamma(0.01, 0.01) because unit-level responses are often on a raw scale.goodness of eblup_fh(), eblup_sfh(), eblup_stfh() and eblup_twofold() follows one
convention. loglikelihood is the full Gaussian log-likelihood at the estimates, and AIC/BIC
are computed from it for ML and REML fits, as in 'sae' and Marhuenda, Morales and Pardo (2014):
the restricted likelihood depends on the covariates, so it cannot compare models with different
covariates and is not reported. goodness always has these three elements.
eblup_twofold(): loglikelihood, AIC and BIC of REML fits were based on the restricted
log-likelihood and now use the full one; BIC uses the number of sampled subareas (it used
the number of areas).eblup_stfh(): the element loglike is renamed loglikelihood, as in the other models.eblup_twofold(): a negative second-order MSE estimate (possible with the ML bias correction)
is set to NA with a warning instead of giving a NaN RSE.eblup_stfh() reports convergence = FALSE when the fit ends at invalid variance components
(it reported TRUE with no estimates). The documentation explains why the model is fitted by
REML only (as in Marhuenda, Molina and Morales 2013 and sae::eblupSTFH()) and that tvalue
is a z-statistic.eblup_stfh(): the analytical $g_3$ was computed with the wrong matrix orientation
($(K_a \Psi V \Psi K_b){ii}$ instead of $\psi_i^2 (K_a V K_b){ii}$), which is wrong whenever
vardir is heterogeneous. $g_3$ now uses the exact derivative of the EBLUP,
$\partial\tilde\theta/\partial\theta_a = \Psi P V_a P y$, and the REML information at the final estimates.
eblup_stfh() bootstrap MSE: the variance components are now re-estimated in every replicate
(before, they were held fixed, so the bootstrap missed the $g_3$ contribution), and the SAR effects
are generated with the lower Cholesky factor (before, $\mathrm{Var}(u_1) = RR' \ne \Omega_1$).
eblup_stfh() now sorts the panel internally into the domain-major order the model requires and returns
results in the caller's row order. Before, unsorted input was silently mis-specified -- including the
time-major mys_panel used in the examples, tests and vignette (where sae::eblupSTFH(), which assumes
domain-major data, gave the same wrong fit).
eblup_stfh(): when the classical Fisher scoring fails, or converges to a negative variance that would
only be truncated to 0, the fit switches to constrained REML (variances >= 0, |rho| <= 0.999; projected
scoring with step halving), which returns the actual constrained maximiser. Convergence in a simulation
went from 63% to 99% of the replicates; fits on which the classical iteration converges to an interior
point are unchanged (identical to sae::eblupSTFH()). A correlation whose variance is 0 is unidentified
and gets an NA standard error.
eblup_stfh() bootstrap: non-converged replicates are replaced by new draws until B valid replicates
are obtained (as in the SFH bootstrap).
eblup_twofold(): the analytical MSE estimator is now $g_1 + g_2 + 2 g_3$ (REML) and
$g_1 + g_2 + 2 g_3 - b_{ML}'\nabla g_1$ (ML); before it was $g_1 + g_2 + g_3$, which ignores the downward bias
of the plug-in $g_1$. The REML Fisher information now includes the cross-area terms.
eblup_sfh() with unsampled domains: the SAR process is defined on all domains of W, so the sampled
domains are now fitted with the marginal covariance $\sigma^2 \Sigma_{full}[s,s]$ (before: the SAR built
from the sub-matrix W[s,s]), consistent with the kriging of the unsampled domains and the bootstrap.
eblup_sfh(): an upper-boundary $\hat\rho = 0.999$ is no longer mapped to $1$ (singular for a
row-standardised W); bootstrap replicates are no longer rejected when estimates exceed 1e6
(scale-dependent check); npbmse no longer fails when $\hat\sigma^2_u = 0$.
hb_unit() with popsize_var: the non-sampled part is now predicted at the non-sampled covariate mean
$\bar x_{d,r} = (N_d \bar X_d - n_d \bar x_{d,s})/(N_d - n_d)$ (before: at $\bar X_d$), and its posterior
SD includes $\sigma_e^2/(N_d - n_d)$.
hb_unit() (Gaussian, non-spatial): INLA's posterior mode search now starts at the REML estimates
of the nested error model. Under the vague log-gamma prior it could stall near $\sigma^2_u = 0$
depending on the number of threads (e.g. $\hat\sigma^2_u \approx 1$ instead of $\approx 68$ on
cornsoybean). A control.mode passed through inla_engine() takes precedence.
benchmark_sae(method = "optimal"): the adjustment is now proportional to $w_i,\mathrm{MSE}_i$, the
minimiser of $\sum_i (\theta_i' - \hat\theta_i)^2/\mathrm{MSE}_i$ subject to the benchmark (before:
$\mathrm{MSE}_i / w_i$). Results are unchanged for equal weights.
eblup_twofold() with one subarea per area (only s2v + s2u identified): Fisher scoring uses the
pseudo-inverse of the singular information.
eblup_bhf(): the bootstrap MSE now uses the same population means (popnmean_xpop) as the estimate.
eblup_fh(transform = "log"): documentation now describes the back-transformation accurately (conditional-mean
correction $\exp(\tilde\mu + g_1/2)$ with delta-method MSE; not the full Slud-Maiti second-order correction),
and self-benchmarking warns that the constraint holds on the log scale only.
eblup_fh(transform = "log"): unsampled domains now always use the synthetic correction
$\exp(x'\hat\beta + \hat\sigma^2_u/2)$, even when their vardir is supplied.
HB models: reported variance components are posterior means of the variance, $E[1/\tau \mid y]$, instead of $1/E[\tau \mid y]$.
eblup_stfh() analytical MSE now includes the second-derivative term $g_4$ (V is nonlinear in
$\rho_1, \rho_2$), as the SFH MSE does: mse = g1 + g2 + 2 g3 - g4, with g4 returned. Simulation:
relative bias -1.7% (was +2.4% without $g_4$, -6.1% without $g_3$).
eblup_sfh(): when the classical Fisher scoring fails or ends with $\hat\sigma^2 < 0$, the (RE)ML
likelihood is maximised under the constraints (as for eblup_stfh()), also inside the bootstrap. An
estimate at $\rho = -0.999$ is reported as $-1$ only when $I + W$ is non-singular (not for bipartite
graphs such as a rook grid, where it made the bootstrap MSE explode).
eblup_sfh() unsampled (kriged) domains: the analytical MSE now includes $g_3$ (exact derivative of
the kriging predictor) and the second-order bias of the plug-in $g_1$ (and the ML bias term):
g1 + g2 + g3 - tr(I^{-1} Hess g1)/2. Simulation (REML): relative bias -3.8% (was -7.5%).
eblup_sfh() / eblup_stfh() analytical MSE: a correlation estimated on the boundary
($|\hat\rho| \ge 0.999$, or reported as $-1$) is treated as fixed -- only the information block of the
other parameters is inverted -- because the Prasad-Rao expansion is not valid there and could explode.
The SE of such a $\rho$ is NA. For SFH this departs from sae::mseSFH in that case only (interior
estimates still agree exactly).
eblup_sfh(method = "ML") analytical MSE: $g_3$ and the ML bias term $b_{ML} = \frac12 I^{-1} h$ now use
the ML information instead of the REML one (as sae::mseSFH does). With $\hat\sigma^2 \approx 0$ and
$\hat\rho \approx 1$ the REML choice made the bias term explode (MSE up to $10^6$ in simulations); with few
domains the ML information also reduced the median MSE bias (e.g. +3.6% vs +6.5%, $m = 36$). The ML
MSE therefore no longer equals sae::mseSFH; REML is unchanged.
eblup_stfh() is much faster: REML scores and Fisher information are computed from the block
structure of $V$ (domain blocks plus a rank-$D$ spatial part) without forming $M \times M$ matrices
($O(T D^3 + M(D+p)^2)$ instead of $O(M^3)$ per iteration; identical results). Fitting the
$40 \times 11$ mys_panel went from about 7.7 s to 0.4 s, and the parametric bootstrap (which
re-estimates the model in every replicate) from about 2.4 s to 0.4 s per replicate.
hb_area(self_benchmark = TRUE) now benchmarks every posterior draw by default
(benchmark_method = "posterior"; difference / logit / ratio adjustment by family), so the constraint
is exact for every likelihood and SD / credible intervals refer to the benchmarked estimator. Benchmark
targets are computed on the parameter scale (proportions for binomial, rates for counts with exposure);
before, the WFQ targets used counts while the estimates were proportions / rates. benchmark_method = "wfq" keeps the previous augmented-model approach.
eblup_fh(transform = "log", self_benchmark = TRUE): ratio benchmarking of the back-transformed
estimates on the original scale (exact), instead of a WFQ constraint that only held on the log scale.
benchmark_sae(): the group MSE used by the hierarchical "optimal" method now includes the covariances
between area prediction errors (full MSE matrix for FH, SFH, STFH and INLA area-level HB fits), and the
output has a new mse_benchmarked column, the first-order MSE of the benchmarked estimates
($mse_i - 2 k_i (Mw)_i + adj_i^2$; within about 1% of the empirical MSE in a simulation).
diagnose(): the Brown et al. (2001) goodness-of-fit statistic $W$ now uses an upper-tail
$\chi^2$ decision (as emdi::compare.fh); the former two-sided interval flagged every correct
FH fit in a simulation, because $W$ is conservative. The bias regression of the direct on the
model estimates is now descriptive (brown_test$bias_regression: alpha, beta, standard
errors, $R^2$) instead of an F test of $\alpha = 0, \beta = 1$, which rejected 93% of correct
FH fits (the slope exceeds 1 for shrinkage estimators); it no longer affects status.
simulation_metrics$rrmse is the relative RMSE across domains (mean_rrmse duplicated
mean_abs_rb). DIC/WAIC/marginal likelihood are now also read for hb_unit() and
hb_twofold() fits, and CPO/PIT are only used when they are per domain.
spatial_test(): the FH vs spatial FH comparison uses ML fits, so the likelihood ratio and
AIC/BIC are computed from ML log-likelihoods at ML estimates (before: at REML estimates). The
recommendation now rests on a single test at level alpha (decision_test: the likelihood
ratio test when vardir is given, otherwise the LM-error test); requiring any of three tests
recommended a spatial model for 9% of non-spatial data sets at alpha = 0.05. The Moran's I,
LM-error and LM-lag statistics are unchanged and agree with spdep.
diagnose(): standardised residuals are $(y - \hat\theta)/\sqrt{vardir - MSE}$ (the variance of
$y - \hat\theta$) instead of $/\sqrt{vardir}$, which shrank them (sd 0.76 under a correct FH
model) and hid outliers in the residual and Q-Q plots.
Reliability tiers (add_reliability_flags(), export_sae(), map_sae(), compare_sae(),
autoplot(type = "rse")) default to the BPS-Statistics Indonesia criteria: RSE <= 25% reliable,
25% < RSE <= 50% use with caution, RSE > 50% unreliable (was 20%/30%, documented as the BPS
standard). map_sae(type = "reliability") now uses its thresholds argument (they were ignored).
benchmark_sae(type = "total") in a hierarchical benchmark without group targets: the group
aggregates are weighted totals (they were weighted means, on the wrong scale for stage 1).
benchmark(), a plain alias of benchmark_sae(), is removed; use benchmark_sae().df_ebp/ebp result tables (from the removed ebp_area()) is removed
from the methods and utilities; results are in df_hb (HB) or df_eblup (EBLUP).diagnose() returns brown_test$bias_regression instead of brown_test$wald_test and
simulation_metrics$rrmse instead of mean_rrmse (see above).map_sae() works with any name of the sf geometry column (it assumed "geometry"), labels
HB fits correctly (they were labelled "EBP-compare_sae() and autoplot() plots label panel and sub-area rows uniquely (repeated domains
made some plots fail or overlap); autoplot() of an unnamed list of models names them.print()/summary() show the estimation method of eblup_fh() and eblup_twofold() (REML or
ML; was "eblup") and describe eblup_stfh(model = "S") correctly.export_sae() lists both variance components of two-fold models in the model summary.benchmark_sae(): the MSE covariance between areas is applied as a matrix-free operator
(FH: $O(n p^2)$, no $n \times n$ matrices). FH with $m = 2000$: 9.4 s / 490 MB -> 0.002 s / 1.7 MB;
$m = 20000$ now runs in about 0.01 s.eblup_bhf(): the REML/ML fit and the bootstrap MSE now run in C++ (src/eblup_bhf_core.cpp).
The likelihood is profiled to $\gamma = \sigma_u^2/\sigma_e^2 \ge 0$ and evaluated from within-domain
cross-products and domain means ($O(Dp^2)$ per evaluation, independent of the number of units); the
maximum is found from a log-grid scan plus Brent refinement, so $\sigma_u^2 = 0$ is handled exactly.
Estimates agree with lme4/sae up to their optimiser tolerance (deviance never worse), and the
bootstrap uses the same random draws as before, so MSEs agree with the previous lme4-based
implementation to about 1e-5 for a given seed. Bootstrap: cornsoybean $B = 200$ 0.64 s -> 0.005 s;
$D = 300$, $n = 6000$, $B = 100$ 0.96 s -> 0.035 s. n_threads now runs replicates in parallel
(results do not depend on it); for small models R's sequential RNG dominates the run time.eblup_bhf() output: fit$lme (the lmerMod object) is replaced by fit$residuals (unit-level
conditional residuals, used by residuals()); estcoef standard errors are $(X'\hat V^{-1}X)^{-1}$
as in lme4. Rank-deficient designs now give an explicit error.eblup_sfh() (fit, analytical MSE, kriging and the bootstrap refits) no longer forms any
$D \times D$ inverse or $D^3$ product. With $A = (I - \rho W')(I - \rho W)$ and
$M = A + \sigma^2 E'\Psi^{-1}E$ (both sparse), $V^{-1}$, $GV^{-1}$, $g_1 = \sigma^2 \mathrm{diag}(M^{-1})$,
the scores/information, $g_3$, $g_4$, the ML bias and the kriging terms follow from Cholesky
solves with $A$ and $M$ (src/sar_engine.h; Eigen sparse $LDL'$ with AMD ordering, dense
factorisations when $W$ itself is dense, e.g. distance weights). Results are unchanged (relative
differences below 2e-10). Analytical: $D = 400$ 2.2 s -> 0.10 s, $D = 1600$ 80 s -> 2.5 s,
$D = 2500$ about 6 s; bootstrap ($D = 400$, $B = 20$): 16 s -> 1.3 s. RcppEigen is added to
LinkingTo.eblup_stfh(): the final step (Fisher information at the estimates, SEs, $g_1$-$g_4$ of the analytical
MSE) no longer forms $M \times M$ matrices ($M = DT$); it uses the same structured representation of $V$
as the iterations ($A^{-1}$ block diagonal, $P = A^{-1} - U W U'$ with $U = [A^{-1}Z, V^{-1}X]$ whose
$A^{-1}Z$ part has one non-zero per row). The scoring step also exploits the diagonal blocks of $U'Z$ and
$U'EU$ and the sparsity of $C = 2\rho W'W - W - W'$. Results are unchanged (relative differences below
3e-11). $T = 5$: $D = 250$ 28 s -> 1.6 s, $D = 500$ about 4 min -> 10 s, $D = 1000$ 75 s.create_weights(): distances for knn, distance and idw are computed in blocks of rows (at most
about 1M distances at a time) instead of a full $D \times D$ distance matrix, and the weight matrix is
filled once from the neighbour pairs (no intermediate 0/1 matrix and no normalisation copies; idw
styles are applied per block). The result is the same matrix as before (identical, or within 1e-16
for idw). Peak memory at $D = 10000$: queen 2.3 GB -> 0.8 GB, knn 3.1 GB -> 0.8 GB, distance
3.8 GB -> 0.8 GB, i.e. essentially the returned matrix itself (764 MB).create_weights() with a coordinate matrix whose columns are named X/Y (e.g. from
sf::st_coordinates()) failed with "undefined columns selected"; coordinate column names are now
matched case-insensitively.eblup_twofold):
hb_twofold):
INLA.df_hb) and weighted area aggregates (df_area) with full Monte Carlo uncertainty propagation.hb_area):
INLA).hb_unit):
INLA.benchmark_sae):
diagnose(): Brown et al. (2001) goodness-of-fit and bias diagnostics, residual analysis, spatial autocorrelation tests (Moran's I), and official statistical reliability classification (BPS standards).compare_sae(): Tabular model comparison across information criteria and error metrics.map_sae(): Spatial thematic choropleth maps integrated with sf.export_sae(): Export estimation results with full statistical metadata to Excel (openxlsx/writexl) and CSV.sim_area_data(), sim_spatial_weights(), and sim_series_data().RcppArmadillo and OpenMP) with strict validation checks.