## ----knitr-opts, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  fig.width = 7.2, fig.height = 4.6, out.width = "100%",
  dpi = 110, fig.align = "center"
)

## ----setup--------------------------------------------------------------------
library(BorderEffect)

## ----montecarlo---------------------------------------------------------------
fl <- field_layout(nx = 12, ny = 8, ntrt = 3, nblk = 2,
                   arrangement = "triangular", width = 550, height = 210,
                   seed = 3000)
W  <- border_weights(fl)
X  <- model.matrix(~ trt + blk, fl)

iter <- 300  # the paper uses 3000; reduced here for a fast vignette
sim_kappa <- function(mean_fun, sd_fun) {
  replicate(iter, {
    y <- rnorm(nrow(fl), mean_fun(fl), sd_fun(fl))
    kappa_hat(y, W, X)
  })
}
k0 <- sim_kappa(function(d) 2.21, function(d) 0.06)                       # no border
k1 <- sim_kappa(function(d) ifelse(d$outer1 == 1, 2.21, 1.74),
                function(d) ifelse(d$outer1 == 1, 0.06, 0.035))          # 1 border

c_no <- "#0072B2"; c_edge <- "#D55E00"   # Okabe-Ito blue / vermillion (CVD-safe)
d0 <- density(k0); d1 <- density(k1)
plot(d0, main = "kappa: no border vs single outer edge", xlab = "kappa",
     xlim = range(k0, k1, 0), ylim = c(0, max(d0$y, d1$y) * 1.05),
     lwd = 2, col = c_no)
lines(d1, lwd = 2, col = c_edge)
abline(v = 0, lty = 3, col = "grey55")
legend("topright", c("no border", "single edge"), col = c(c_no, c_edge),
       lwd = 2, bty = "n")