--- title: "Precision and ROPE design analysis" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Precision and ROPE design analysis} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = FALSE, comment = "") # Console colour carries no meaning on a rendered page. pkgdown turns it on for # its own build, and the escape sequences then reach the reader as literal text, # so colour is switched off here for a plain vignette render and a site build # alike. The fixed width keeps printed output inside the documentation column. options(cli.num_colors = 1, cli.hyperlink = FALSE, crayon.enabled = FALSE, width = 80) ``` ```{r setup} library(pilotr) ``` When sample sizes are large, or when the question of interest is whether an effect is *practically* meaningful, power against a point null is often the wrong target. pilotr implements a precision-based design analysis against a region of practical equivalence (ROPE). This is a fast, frequentist analogue of the Bayesian approach that compares a highest-density interval with a ROPE. > As in the power vignette, the fixed-*N* analysis below uses a small `n_sims` so that the > vignette builds quickly, and the sample-size sweep, which needs far more model fits, was > precomputed with exactly the code shown and shipped with the package. For real planning we > recommend `n_sims >= 200`, the value the sweep uses. ## The idea Over Monte Carlo replicates, and for each focal fixed effect, pilotr records the 95% confidence interval and whether it falls determinately outside the ROPE (the effect is practically meaningful) or entirely inside it (practical equivalence to zero), together with the expected interval width. The interval is a Wald approximation, the estimate plus or minus 1.96 standard errors, chosen for speed and for comparability across replicates. Sweeping sample size then locates the minimum *N* at which a focal effect reaches a determinate decision with a target probability. ## A worked design We reuse the crossed priming design, which represents a small priming effect on log reaction time. ```{r} spec_c <- build_spec(list( name = "priming", seed = 1, design_kind = "within", include_items = TRUE, n_subject = 24, n_item = 18, factor_name = "condition", lev1 = "related", lev2 = "unrelated", intercept = 6, effect = 0.05, subj_int_sd = 0.12, subj_slope_sd = 0.04, subj_corr = 0.2, item_int_sd = 0.08, item_slope_sd = 0.02, item_corr = -0.1, family = "shifted_lognormal", resp_name = "RT", sigma = 0.3, shift = 200)) ``` So that the analysis can run from the specification alone, pilotr auto-derives the analysis model, including the contrast columns, the response transform and the mixed-model formula. ```{r} model_formula(spec_c) ``` The companion `model_data()` adds the other two derivations to a simulated data set, the numeric `effect` contrast column the formula's terms refer to and the `.y` response, here the log-transformed RT because the family is `shifted_lognormal`. ```{r} head(model_data(spec_c, simulate_design(spec_c))) ``` ## Precision at a fixed *N* We declare the focal effect (its coefficient name and true value) and a ROPE half-width. Here the effect lives on the log scale, and we treat anything smaller than 0.02 as practically equivalent to zero. ```{r} pr <- precision_design( spec_c, focal = c(effect = 0.05), rope = 0.02, n_sims = 25) pr ``` The columns are interpreted as follows. `p_meaningful` is the probability the 95% CI lands entirely outside the ROPE, a determinate 'meaningful' decision. `p_equivalent` is the probability it lands entirely inside, a determinate 'equivalent' decision. `mean_ci_width` is the expected precision. ## Sweeping sample size The same ROPE has to be carried into the sweep. Leaving `rope` at its default would compare the interval against a region as wide as the effect itself, and the decision probability would then fall away as *N* grew, which is backwards. ```{r, eval = FALSE} prc <- precision_curve(spec_c, focal = c(effect = 0.05), subject_ns = c(15, 30, 60, 100, 140, 180, 220, 260), rope = 0.02, n_sims = 200) prc[, c("n_subject", "p_meaningful", "p_equivalent", "mean_ci_width", "n_converged")] ``` ```{r, include = FALSE} # The result of running exactly the chunk above, precomputed and shipped with # the package so that the vignette builds within CRAN's check-time budget. All # 200 replicates converged at every sample size. prc <- read.csv("precision-curve-cache.csv") ``` ```{r, echo = FALSE} prc[, c("n_subject", "p_meaningful", "p_equivalent", "mean_ci_width", "n_converged")] ``` As *N* grows the interval tightens and `p_meaningful` rises, so the design reaches a determinate decision more reliably. That is a more informative criterion than power against a point null, because it asks whether the study can distinguish the effect from a negligible one, where power asks only whether it can be distinguished from zero. Each estimate rests on 200 replicates, all of which converged, so its Monte Carlo standard error is at most about 0.035, which the figure below shows as a band. Scanning the table for the smallest simulated *N* that meets a target of '`p_meaningful` ≥ 0.90' puts this design at 220 subjects, but that answer is a property of where the grid points happen to fall. `solve_curve()` fits the decision probability against the sample size and inverts the fit, so all eight points inform the answer, where scanning the table uses only the two either side of the target. ```{r} solved <- solve_curve(prc, target = 0.9) unlist(solved[c("value", "lo", "hi", "dispersion")]) ``` The sweep read by eye therefore asks for some thirty-five subjects more than the fitted crossing does, and 220 sits outside the interval altogether. The fit is a binomial probit regression weighted by the replicates behind each point, and the interval is the delta-method interval that `MASS::dose.p()` computes for a fitted `glm`. The interval matters more here than the point it surrounds: with 200 replicates per point this design needs somewhere between about 165 and 207 subjects to reach a determinate ROPE decision nine times in ten, and a preregistration quoting a single figure would be claiming more than the simulation supports. The reported `dispersion` above 1 is part of that width. It says the eight points scatter about the fitted curve rather more than their replicate counts alone would explain, so the interval has been widened to match, which is the heterogeneity factor of classical probit analysis doing its work. Nothing is extrapolated either: a curve that never reaches the target within the sizes swept is refused, and the refusal reports the range the sweep did cover. `p_equivalent` stays at 0 throughout, and correctly so. The true effect of 0.05 lies outside the ROPE by construction, so no interval should ever land entirely inside it. The column earns its place in designs where practical equivalence is the hypothesis of interest. ```{r, fig.width = 6.5, fig.height = 3.2, dev.args = list(bg = "transparent")} library(ggplot2) # Each probability is a proportion over the converged replicates, so it carries # a binomial Monte Carlo standard error, and the band is its 95% interval. The # two panels are on different scales, hence the free y axis. prc$se <- sqrt(prc$p_meaningful * (1 - prc$p_meaningful) / prc$n_converged) panels <- c("P(CI outside the ROPE)", "Mean 95% CI width") long <- rbind( data.frame(n_subject = prc$n_subject, panel = panels[1], y = prc$p_meaningful, lo = pmax(0, prc$p_meaningful - 1.96 * prc$se), hi = pmin(1, prc$p_meaningful + 1.96 * prc$se)), data.frame(n_subject = prc$n_subject, panel = panels[2], y = prc$mean_ci_width, lo = NA, hi = NA)) long$panel <- factor(long$panel, levels = panels) # 0.90 is the target decision probability. 0.06 is the width at which a CI # centred on the true effect just clears the ROPE, that is 2 * (0.05 - 0.02). refs <- data.frame(panel = factor(panels, levels = panels), ref = c(0.90, 0.06)) # The solved sample size and its interval, drawn only on the decision panel, # since the width panel is on a different scale and answers a different question. solve_band <- data.frame(panel = factor(panels[1], levels = panels), lo = solved$lo, hi = solved$hi, at = solved$value) ggplot(long, aes(n_subject, y)) + geom_hline(data = refs, aes(yintercept = ref), linetype = 2, colour = "grey60") + geom_rect(data = solve_band, inherit.aes = FALSE, aes(xmin = lo, xmax = hi, ymin = -Inf, ymax = Inf), fill = "grey60", alpha = .15) + geom_vline(data = solve_band, aes(xintercept = at), linetype = 2, colour = "grey60") + geom_ribbon(aes(ymin = lo, ymax = hi), alpha = .15, fill = "#2C6FB0", na.rm = TRUE) + geom_line(colour = "#2C6FB0", linewidth = 0.8) + geom_point(colour = "#2C6FB0", size = 2.2) + facet_wrap(~ panel, scales = "free_y") + labs(x = expression(italic(N) ~ "subjects"), y = NULL) + theme_minimal(base_size = 12) + # theme_minimal still paints a white plot.background over the transparent # canvas, so both surfaces have to be cleared for the page colour to reach the # figure. theme(plot.background = element_rect(fill = NA, colour = NA), panel.background = element_rect(fill = NA, colour = NA), panel.grid = element_line(colour = "grey80"), strip.background = element_rect(fill = NA, colour = NA)) ``` The two panels answer the same question from either side. On the left, the horizontal dashed line marks the 0.90 target, and the vertical band the sample size solved against it. On the right the dashed line marks a width of 0.06, which is where a confidence interval centred on the true effect just clears a ROPE of 0.02, since the interval has to keep its lower limit above the ROPE and so must be narrower than 2 × (0.05 − 0.02). Width falls below that line well before the decision probability reaches 0.90, because an interval centred exactly on the true effect is the best case and sampling variation moves the centre about. ## See also The [power vignette](power-analysis.html) covers simulation-based power and Type S/M errors. The [getting-started vignette](getting-started.html) covers the core simulate-and-inspect loop.