Simulating Covariance Data to Plan a Study

The displayed results and figures in this guide were computed from the executable R Markdown source shipped with the package. Routine package checks use these precomputed results to keep the CRAN check time modest. To rerun every code chunk, copy system.file("vignette-source", "simulation-design.Rmd", package = "terradish") to a writable directory and render that copy with rmarkdown::render().

Goal. Simulate a response under stated assumptions and plan how to compare sampling designs.

The data-generation guide explains the landscapes, SLiM histories, sampling, and matched markers used in these examples. Open it in R with vignette("data-generation", package = "terradish").

How many ponds would you sample?

Suppose field effort allows visits to 12, 18, or all 24 candidate ponds, and before sampling you want to know whether those designs can recover a specified forest association. This guide sets up that question on the Ohio pond landscape from vignette("getting-started", package = "terradish").

The response is drawn directly from a covariance model. These are new statistical draws, not the SLiM SNPs or microsatellites. They give a fast check of recovery when the analysis model is correct, whereas SLiM generates ancestry and mutations through a biological process that need not obey the fitted likelihood.

The guide runs one draw and one fit, and larger design comparisons are copyable recipes that routine builds skip. One successful fit cannot establish statistical power, interval coverage, or field sample size.

Keep the landscape and specify a new response

library(terradish)
library(terra)
data(ohio_ponds)
ponds <- ohio_ponds$ponds

# Reuse the forest context, full graph, and candidate pond locations.
# saveStack retains the raster needed when designs subset focal ponds.
x <- scale_covariates(c(unwrap(ohio_ponds$forest500),
                        unwrap(ohio_ponds$covariates)[["slope"]]))
surface <- conductance_surface(x, ponds[c("x", "y")],
  directions = 8, saveStack = TRUE)
c(candidate_ponds = length(surface$demes), resolution_m = res(x)[1])
#> candidate_ponds    resolution_m
#>              24             100

# A declared one-SD association, not an estimate from the SLiM data.
theta_true <- c(forest500 = 0.6)
tau_true <- 2
variance_true <- 1

How to read the setup. Each design selects ponds from the same 24 candidates while retaining routes through the full 100 m landscape. The forest coefficient 0.6 implies a conductance ratio of exp(0.6), about 1.82, per standard deviation of forest context. Fixing context at 500 m isolates the site-design question from scale estimation.

The simulator uses

\[\Sigma=\tau E(\theta)+vI,\qquad \nu S\sim\mathrm{Wishart}(\nu,\Sigma).\]

Here \(E\) is graph covariance, \(\tau\) weights it, and \(v\) is additional diagonal variance. The simulation call’s sigma takes \(v\) itself, but the fitted measurement parameter named sigma is \(\log(v)\). Neither is the Gaussian raster-smoothing scale.

Generating with nu = 100 sets how noisy each simulated response is, not how many SNPs would supply that information. The simulator requires nu of at least the number of sampled ponds (24 here).

Generate once and inspect recovery

sim <- simulate_covariance_response(theta_true, ~ forest500, surface,
  tau = tau_true, sigma = variance_true, nu = 100, seed = 1)
#> Warning in assemble_model_matrix(formula, x): Removed unused spatial
#> covariates: slope
#> Warning in assemble_model_matrix(formula, x): Removed unused spatial
#> covariates: slope
#> Warning in assemble_model_matrix(formula, x): Removed unused spatial
#> covariates: slope
control <- NewtonRaphsonControl(maxit = 200, ctol = 1e-8, ftol = 1e-8)
fit <- terradish(sim$covariance ~ forest500, surface,
  measurement_model = wishart_covariance, nu = 100,
  optimizer = "bfgs", control = control)
#> Warning in assemble_model_matrix(formula, x): Removed unused spatial
#> covariates: slope
stopifnot(fit$convergence$code == 0,
          fit$fit$subproblem$convergence == 0)
rbind(generating = theta_true, estimated = coef(fit))
#>            forest500
#> generating 0.6000000
#> estimated  0.5924285
confint(fit)
#>               2.5 %    97.5 %
#> forest500 0.3257009 0.8591562

How to read the output. The estimated forest coefficient is 0.59, against the generating 0.6, and its 95% interval (0.33 to 0.86) contains the true value. Only the design-sweep recipe below can estimate how often such intervals do. The harmless warnings note that this formula ignores slope, which stays in the raster stack for the multivariate recipes.

plot(fit, type = "fit")

One simulated covariance: observed versus fitted values.

How to read the plot. The upper cluster holds the 24 pond variances and the lower cluster the covariances between ponds. The blue line is a least-squares fit, not the 1:1 line. Most visual agreement comes from separating variances from covariances, so judge the fit within the lower cluster. plot(fit, type = "marginal", data = surface) shows the fitted forest-conductance curve for comparison with exp(0.6 * x).

Compare site designs without changing the question

Before a design sweep, choose a recovery target, plausible effect sizes, a range of information assumptions, and a computation budget. The recipe below compares spatially spread and randomly selected ponds. sample_sizes counts ponds, not animals per pond or loci, so the helper does not simulate the benefit of collecting 15 versus 30 animals at a pond.

# Start with nsim = 5 to check runtime and failures, then increase replication.
# All candidates use the same response within each simulated replicate.
power <- covariance_response_power(
  theta = theta_true, formula = ~ forest500, data = surface,
  sample_sizes = c(12, 18, 24), strategies = c("spacefill", "random"),
  fit_models = list(forest = list(formula = ~ forest500),
                    forest_slope = list(formula = ~ forest500 + slope)),
  tau = tau_true, sigma = variance_true,
  nu = 100, nu_fit = 100, nsim = 200, n_designs = 2, seed = 905,
  optimizer = "bfgs", cores = 1, control = control)
power$summary
power$parameter_summary

How to read the eventual results. These tables appear only when you run the recipe. Read fitting success before recovery, because failures are part of the assessment, not rows to discard. Slope has no generating effect in this benchmark, so its selection rate shows whether model comparison favors unneeded structure.

Quantity	Meaning	Decision it supports
`fit_rate`	Fraction of attempted fits that converge	Is the proposed design numerically usable?
Parameter `power`	Fraction of attempts detecting the generating sign at the stated threshold	Is the specified association detected often enough?
`bias`, `rmse`	Average signed error and root mean squared error	Are estimates accurate as well as detectable?
`coverage`	Fraction of available intervals containing the generating value	Is uncertainty calibrated among the intervals being summarized?
`conductance_power`	Fraction recovering the spatial surface above the specified correlation threshold	Is the fitted surface pattern recovered?
`selected_AIC_rate`	Frequency a candidate ranks first in this candidate set	Does model comparison favor the generating association?

Check denominators alongside these summaries. Failed or boundary fits count as non-detections in parameter power, while bias and coverage depend on available estimates and intervals. High nominal power with poor coverage is not a sound design. Surface correlation measures relative spatial pattern, not absolute conductance or movement.

Separate data information from assumed precision

nu sets the variability of the generated data, and nu_fit sets the information assumed when fitting. Raising only nu_fit narrows reported intervals without improving the data. Compare the two explicitly:

design <- list(theta = theta_true, formula = ~ forest500, data = surface,
  sample_sizes = 18, strategies = "spacefill",
  tau = tau_true, sigma = variance_true, nu = 100,
  nsim = 200, seed = 905, optimizer = "bfgs", cores = 1, control = control)
matched <- do.call(covariance_response_power, design)
overstated <- do.call(covariance_response_power,
  c(design, list(nu_fit = 1000)))
rbind(matched = matched$parameter_summary,
      overstated = overstated$parameter_summary)

How to read this comparison. Common seeds pair the simulated draws, so compare recovery, coverage, and fitting success together. A higher detection rate under overstated precision may simply reflect intervals that are too narrow. Neither setting establishes the effective nu of the bundled SNP panel.

What can you take to the field?

Report performance across plausible weak and strong associations instead of choosing the scenario that makes the smallest sample sufficient. For a proportion estimated from \(B\) independent replicates, approximate Monte Carlo standard error is \(\sqrt{p(1-p)/B}\), about 0.028 with 200 replicates and \(p=0.8\). Reused site designs add another layer to consider when summarizing uncertainty across a sweep.

Save settings, seeds, all fit outcomes, package versions, and the complete assessment object. Keep graph resolution, buffer, and forest scaling fixed when comparing site counts. If the goal includes choosing a scale or spline complexity, simulate and evaluate that selection procedure too.

Summary of key functions

Function	Purpose
`conductance_surface()`	Retain the landscape and candidate focal sites
`simulate_covariance_response()`	Draw one response under declared statistical assumptions
`terradish()`	Fit and inspect one draw before a larger assessment
`covariance_response_power()`	Compare recovery across specified sampling designs