Validation on simulated fields

The idea

On real uniformity trials the optimal plot size is unknown – that is the whole reason to estimate it. To validate the methods we need the opposite: data whose spatial structure is known by construction, so an estimate can be compared against a truth. This article simulates such fields and checks what the package recovers.

The fields are drawn from a Gaussian process with a known exponential variogram. Two quantities are then known exactly:

• the range of the variogram, a, the distance over which basic units stay correlated; and
• the lag-1 autocorrelation between adjacent units, \(\rho = \exp(-3/a)\), which follows from the effective-range convention the package uses (\(\gamma\) reaches 95% of the sill at a).

The generator

A field is one draw from a multivariate normal whose covariance decays exponentially with distance. Building it needs nothing beyond base R:

gen_field <- function(nr, nc, range, mu = 100, psill = 350, nugget = 50) {
  xy <- expand.grid(r = seq_len(nr), c = seq_len(nc))
  D  <- as.matrix(dist(xy))
  S  <- psill * exp(-3 * D / range)   # structured covariance, 95% decay at `range`
  diag(S) <- psill + nugget           # nugget: variance with no spatial structure
  z  <- mu + t(chol(S)) %*% rnorm(nr * nc)
  matrix(as.numeric(z), nrow = nr, ncol = nc)
}

A larger range means smoother fields: neighbours resemble one another over a longer distance. The three fields below, on a 16 x 16 grid, share one colour scale, so the increasing spatial structure is visible directly.

set.seed(99)
fields <- lapply(c(2, 4, 6), function(a) gen_field(16, 16, range = a))
names(fields) <- paste("range =", c(2, 4, 6))

plot(check_trial(fields), title = "Simulated fields (kriged surface)")
#> Checking 3 trial(s).

Recovering the autocorrelation

The lag-1 autocorrelation is the quantity the Paranaíba method rests on, and the one a short trial estimates most stably. For each range we draw several fields, estimate \(\rho\) with calc_paranaiba() (the average of its two directional walks), and compare it with the truth \(\exp(-3/a)\).

set.seed(2026)
ranges <- c(2, 4, 6)
rho <- do.call(rbind, lapply(ranges, function(a) {
  do.call(rbind, lapply(1:10, function(i) {
    s <- suppressMessages(calc_paranaiba(gen_field(16, 16, a)))$summary
    data.frame(range = a, rho_true = exp(-3 / a),
               rho_est = mean(c(s$rho_row, s$rho_col)))
  }))
}))

agg <- aggregate(cbind(rho_true, rho_est) ~ range, rho, mean)
agg$bias <- agg$rho_est - agg$rho_true
round(agg, 3)
#>   range rho_true rho_est   bias
#> 1     2    0.223   0.193 -0.030
#> 2     4    0.472   0.380 -0.092
#> 3     6    0.607   0.498 -0.108

The estimate tracks the truth across the whole range, with a small downward bias – the familiar attenuation of a lag-1 autocorrelation estimated from a short series. The relationship is tight across individual draws, not just on average:

library(ggplot2)

ggplot(rho, aes(rho_true, rho_est, colour = factor(range))) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, colour = "grey50") +
  geom_point(size = 2.4, alpha = 0.8) +
  labs(title = "Estimated vs true lag-1 autocorrelation",
       subtitle = sprintf("correlation across all draws: %.2f",
                          cor(rho$rho_true, rho$rho_est)),
       x = expression("true  " * rho == exp(-3/a)),
       y = expression("estimated  " * rho),
       colour = "range") +
  coord_equal() +
  theme_minimal(base_size = 12)

Points sit just below the identity line: the estimator is slightly conservative but unbiased in ordering, which is what matters when \(\rho\) feeds a plot-size formula.

Response of the CV-curve methods

The plateau methods do not estimate \(\rho\); they read the optimal plot size off the CV-versus-size curve built by calc_cv_shapes(). A field with more spatial structure keeps rewarding larger plots for longer, so the optimum should grow with the range. It does:

set.seed(7)
cv <- do.call(rbind, lapply(ranges, function(a) {
  do.call(rbind, lapply(1:4, function(i) {
    tab <- suppressMessages(calc_cv_shapes(gen_field(16, 16, a)))
    lrp <- fit_lrp(tab, x = "x", cv = "cv", step = 0.25)
    qrp <- fit_qrp(tab, x = "x", cv = "cv", step = 0.25)
    data.frame(range = a,
               LRP = unname(lrp$parameters["Breakpoint"]),
               QRP = unname(qrp$parameters["Breakpoint"]))
  }))
}))
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).
#> Using x = 'x', cv = 'cv' (single series).

aggregate(cbind(LRP, QRP) ~ range, cv, mean)
#>   range     LRP     QRP
#> 1     2 10.6875 16.7500
#> 2     4 18.3125 27.1250
#> 3     6 20.3125 33.0625

Both optima increase monotonically with the generating range: the methods respond to a controlled change in spatial structure in the expected direction.

What this does and does not show

The simulation validates two things cleanly. First, the autocorrelation machinery underneath check_trial() and calc_paranaiba() recovers a known \(\rho\). Second, the CV-curve optima move in the right direction as spatial structure is dialled up.

It also makes an honest limitation visible. There is no single “true optimal plot size” to recover: the methods target different definitions and answer different questions. The plateau methods (LRP, QRP) grow with the range, while the Paranaíba optimum moves the other way – it is largest at \(\rho = 0\) and shrinks as dependence strengthens, by construction of its formula (see vignette("paranaiba")). That divergence is a documented property, not an error, and it is exactly what vignette("compare") is for.

A second caveat worth stating: the point estimate of the variogram range from a single small trial is noisy, so it is not used as a validation target here. The nugget-to-sill ratio and the autocorrelation are the more stable readings, which is why check_trial() leads with them.

Extending this to your own scenarios

The generator is the hook. Vary its arguments to build the cases you care about:

• raise nugget toward psill to simulate a field with weak spatial dependence, where every method should return a small optimum and check_trial() should report a high nugget-to-sill ratio;
• change the grid dimensions to study how few basic units the methods tolerate;
• replace the exponential covariance with a spherical or Gaussian one to check robustness to the variogram shape.

Wrapping the loops above in a function that returns estimate-minus-truth is the whole of a validation harness. See vignette("check_trial") for the diagnostics and vignette("compare") for reading several methods against each other.