---
title: "Validation on simulated fields"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Validation on simulated fields}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5,
                      fig.height = 4.0, dpi = 120)
```

```{r setup}
library(trialSizing)
```

## 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:

```{r generator}
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.

```{r fields, fig.height = 3.2}
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)")
```

## 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)$.

```{r rho-recovery}
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)
```

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:

```{r rho-plot, fig.height = 4.2}
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:

```{r cv-response}
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"]))
  }))
}))

aggregate(cbind(LRP, QRP) ~ range, cv, mean)
```

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.