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

## ----setup--------------------------------------------------------------------
library(trialSizing)

## ----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)
}

## ----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)")

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

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

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