--- title: "Simulating Covariance Data to Plan a Study" author: "Bill Peterman" output: rmarkdown::html_vignette: number_sections: true toc: true self_contained: true mathjax: null pandoc_args: ["--mathml"] vignette: > %\VignetteIndexEntry{Simulating Covariance Data to Plan a Study} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- 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](data-generation.html) 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")`](getting-started.html). 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 ``` r 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 ``` r 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. ``` r plot(fit, type = "fit") ``` ![One simulated covariance: observed versus fitted values.](precomputed/simulation-design-recovery-plot-1.png) **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. ``` r # 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: ``` r 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 | ## See also - [`vignette("getting-started", package = "terradish")`](getting-started.html): the SLiM deme example and matched marker panels. - [`vignette("wishart-covariance", package = "terradish")`](wishart-covariance.html): response geometry and effective information. - [`vignette("model-comparison", package = "terradish")`](model-comparison.html): selection and prediction are different assessments. - [`vignette("ibe-ibr-workflow", package = "terradish")`](ibe-ibr-workflow.html): a simulated response with a known environmental component. - [`vignette("large-landscapes", package = "terradish")`](large-landscapes.html): solver settings when a design sweep requires many fits.