--- title: "Planning a uSDT study" description: > How many subjects, and how many trials in each task? The power of the three hypotheses, the cheapest design that reaches it, the smallest effect a design detects and a check by simulation, with the contextual cuing data of Vadillo, Malejka and Shanks (2025) as a pilot study. output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Planning a uSDT study} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r} #| label: setup #| include: false knitr::opts_chunk$set( collapse = FALSE, comment = "", fig.width = 8, fig.height = 5.2, fig.align = "center", out.width = "100%", dpi = 150 ) # The printed output rules the tables with box characters in a UTF-8 locale # and falls back to ASCII when there is none, so the vignette asks for UTF-8 # to read the same wherever it is built. It only asks when the session is not # already there, and the last chunk puts the width and the locale back as they # were. old_options <- options(width = 80) old_ctype <- Sys.getlocale("LC_CTYPE") if (!isTRUE(l10n_info()[["UTF-8"]])) { for (loc in c("en_US.UTF-8", "C.UTF-8", "English_United States.utf8")) { suppressWarnings(Sys.setlocale("LC_CTYPE", loc)) if (isTRUE(l10n_info()[["UTF-8"]])) break } # Windows accepts locale names it does not actually have, so trust # l10n_info() to verify the locale, and undo the attempts when none # of the names produced a UTF-8 session. if (!isTRUE(l10n_info()[["UTF-8"]])) { invisible(suppressWarnings(Sys.setlocale("LC_CTYPE", old_ctype))) } } ``` ```{r} #| label: library library(uSDT) ``` ## Three numbers to choose A study of unconscious processing has three sizes to choose: how many subjects, how many trials of the awareness task and how many of the influence task. Trials affect measurement error in each subject's sensitivity ($d'$). Subjects affect the precision of population estimates. The hypotheses depend on them in different ways: - **H1**, the difference between mean sensitivities, gains from both subjects and trials, like any comparison of means. - **H2**, the latent correlation between task sensitivities, needs enough trials to separate between-subject variation from measurement error. - **H3**, the latent intercept, estimates influence sensitivity at zero awareness sensitivity. Its precision also depends on the slope. The workflow is to specify a population, compute power, choose a feasible design, vary assumptions and check the chosen design by simulation. ## The population Planning starts from the values the population is assumed to have. A pilot study is a natural source. The data of Vadillo, Malejka and Shanks (2025), which the package includes, serve as one here: ```{r} #| label: pilot d <- usdt_data_tasks( awareness = vadillo_awareness, influence = vadillo_cuing, subject_col = "subj", condition_col = "condition", condition_levels = c(signal = "old", noise = "new"), response_col = list(awareness = "judged.old", influence = "rt"), response_levels = list(awareness = c(signal = 1, noise = 0), influence = c(signal = "faster", noise = "slower")), dichotomize = list(awareness = FALSE, influence = TRUE) ) pilot <- hsdt(d) ``` [`usdt_population()`](../reference/usdt_population.html) takes the estimates of a fit with `from`, and any value given in the call replaces the estimate. For illustration, we target an intercept of 0.1. Choose the smallest effect of interest from substantive arguments. Neither this example nor the pilot determines it. Giving the intercept uses the latent regression parameterization, which keeps the awareness mean, both standard deviations and the correlation, and moves the influence mean so that the intercept is 0.1: ```{r} #| label: population pop <- usdt_population(from = pilot, intercept = 0.1) pop ``` The arguments have the names of [`usdt_priors()`](../reference/usdt_priors.html) and take numbers: `dprime`, `sd_dprime`, `cor_dprime`, `criterion`, `sd_criterion` and `intercept`. A population can also be written from scratch, in either parameterization, as `usdt_population(dprime = c(awareness = 0.2, influence = 0.3), sd_dprime = 0.3, cor_dprime = 0.4)`. The plot shows the distribution of the true $d'$ across subjects, with the latent regression line and its intercept. Given the trials of each task, it adds, as a dashed ellipse, the distribution of the $d'$ those trials measure: ```{r} #| label: population-plot #| fig.width: 6.5 #| fig.height: 6.5 #| out.width: "75%" plot(pop, trials = c(awareness = 64, influence = 384), dichotomize = "influence") ``` The subjects of this population differ very little in their awareness sensitivity (a standard deviation of 0.12), much less than 64 trials can tell apart. The measured awareness $d'$ spreads almost three times as widely as the true one: the reliability of the awareness task is 0.13. Most of what a study would see in the awareness task is measurement error, and that is what limits H2 and H3. ## The power of a design [`usdt_power()`](../reference/usdt_power.html) gives the power of each hypothesis for a number of subjects and of trials per task. The hypotheses are written as in [`usdt_equivalence()`](../reference/usdt_equivalence.html), so they can be two-sided, one-sided, equivalence or minimal-effect tests. Here we use the pilot's design, testing the correlation two-sided and the difference and intercept one-sided: ```{r} #| label: power usdt_power(pop, subjects = 104, trials = c(awareness = 64, influence = 384), hypothesis = c("diff < 0", "rho = 0", "intercept > 0"), dichotomize = "influence", combine = "all") ``` The `all` row is the probability of rejecting all three nulls in the same study. It is not the product of their separate powers. If a population value lies in H0, the reported rejection rate is the type I error rate. The pilot's design would detect an intercept of 0.1 in one study out of three. The `Estimable` column estimates how often the sensitivity covariance stays off the boundary in the approximation. Non-estimable draws count as failures to reject H2 and H3. Several numbers of subjects and of trials give a grid of designs, and `plot()` draws power against the subjects: ```{r} #| label: power-curves curves <- usdt_power(pop, subjects = seq(50, 800, 25), trials = list(awareness = c(64, 256), influence = 384), hypothesis = c("rho = 0", "intercept > 0"), dichotomize = "influence") plot(curves) ``` Four times as many awareness trials do more for H2 and H3 than many more subjects: they raise the reliability of the awareness task from 0.13 to 0.37. ### How the power is computed Power uses a sampling-distribution approximation by default. Measurement error in each subject's $d'$ has the variance that [`usdt_reliability()`](../reference/usdt_reliability.html) reports, which falls with the number of trials. The measured $d'$ of $N$ subjects then have a known covariance, and their sample mean and covariance a known distribution. `usdt_power()` draws sample means and covariances 20,000 times, subtracts measurement variance, applies approximate tests and counts rejections. Invalid sensitivity covariance estimates count as boundary outcomes for H2 and H3. This approximates the hierarchical analysis. It does not fit a model to each draw. With `method = "simulation"`, `usdt_power()` instead simulates whole studies, splits the influence responses at each subject's median as [`meyen_split()`](../reference/meyen_split.html) does, fits each one with `hsdt()` and tests it. This slower method checks the chosen design below. ## Choosing subjects and trials [`usdt_plan()`](../reference/usdt_plan.html) chooses the design. A range lets it choose a value and a single number fixes it. With subjects and trials free, many designs reach the same power, trading subjects for trials, and what each costs decides among them. Here the cost is the time of the study in seconds: 600 per subject for instructions and breaks, 5 per trial of the recognition test and 2 per trial of visual search. ```{r} #| label: plan plan <- usdt_plan(pop, hypothesis = "intercept > 0", power = 0.8, subjects = c(50, 1000), trials = list(awareness = c(32, 512), influence = c(192, 768)), cost = c(subject = 600, awareness = 5, influence = 2), dichotomize = "influence") plan ``` The cheapest design is a small study with long tasks, and both trial counts end at the limit of their ranges: with an awareness task this unreliable, more trials are cheaper than more subjects. A recognition test of 500 trials, however, may not be feasible, since testing the displays so many times could change what is remembered. With at most 128 awareness trials: ```{r} #| label: plan-limited plan <- usdt_plan(pop, hypothesis = "intercept > 0", power = 0.8, subjects = c(50, 2000), trials = list(awareness = c(32, 128), influence = c(192, 768)), cost = c(subject = 600, awareness = 5, influence = 2), dichotomize = "influence") plan plot(plan) ``` The plot shows the cost of the cheapest design for each number of awareness trials. The grey line is the closed form that the search starts from. The green points are the designs the search checked with the draws of `usdt_power()`. The circle is the chosen one. With an unreliable awareness task the closed form is optimistic, so the search does not trust it and checks a lattice of designs and the neighbours of the best. ### Budget and precision targets With fixed trials, these calls choose the number of subjects. `budget` maximizes power within a cost limit. `precision` targets an expected confidence interval width instead of power. ```{r} #| label: plan-targets usdt_plan(pop, hypothesis = "intercept > 0", budget = 30 * 3600, subjects = c(50, 2000), trials = c(awareness = 128, influence = 384), cost = c(subject = 600, awareness = 5, influence = 2), dichotomize = "influence") usdt_plan(pop, precision = c(intercept = 0.2), subjects = c(50, 2000), trials = c(awareness = 128, influence = 384), dichotomize = "influence") ``` ## Other assumptions The plan rests on the population. Several values in an argument of `usdt_population()` give one scenario each, and [`update()`](../reference/usdt_population.html) changes them. With the awareness task at 128 trials and the influence task at 384, a larger intercept or a weaker correlation than assumed needs another number of subjects: ```{r} #| label: scenarios scenarios <- update(pop, intercept = c(0.1, 0.15), cor_dprime = c(0.3, 0.49)) usdt_plan(scenarios, hypothesis = "intercept > 0", power = 0.8, subjects = c(50, 3000), trials = c(awareness = 128, influence = 384), dichotomize = "influence") ``` An intercept of 0.15 needs about half the subjects that 0.1 needs. A weaker correlation needs fewer subjects, not more: with a flatter slope, the intercept depends less on how far the regression is extrapolated to an awareness sensitivity of zero. ## The smallest effect the design detects [`usdt_mdes()`](../reference/usdt_mdes.html) finds, for each hypothesis, the smallest effect a design detects with the given power. It moves the tested quantity away from H0 with the rest of the population fixed, as its parameterization says. With the chosen design: ```{r} #| label: mdes mdes <- usdt_mdes(pop, subjects = plan$design$subjects, trials = c(awareness = plan$design$trials_awareness, influence = plan$design$trials_influence), hypothesis = c("diff < 0", "rho = 0", "intercept > 0"), dichotomize = "influence") mdes plot(mdes) ``` The power of H2 does not keep rising with the correlation: near 1, the estimated correlation often reaches the boundary, where the test cannot be run. For an equivalence test, `usdt_mdes()` gives instead the narrowest equivalence region the design can establish. ## Checking the design by simulation A design chosen with the approximation should be checked by simulating it. `usdt_power()` takes the plan and simulates its designs, analysing each study with `hsdt()`. The following is optional and is not run when building the article: ```{r} #| label: check #| eval: false check <- usdt_power(plan, method = "simulation", nsim = 1000, ncores = 2, seed = 2026) check summary(check) ``` The report gives the simulated power with its Monte Carlo error and 95% interval, and how many studies failed to converge. `summary()` reports the precision of each quantity across the simulated studies: the mean estimate and its bias, the empirical standard error (the standard deviation of the estimates) against the mean standard error the fits report, and the width and coverage of the intervals. Each study has its own random number stream, so the results do not depend on the number of cores. ## How close is the approximation? The figure compares, for eight designs of two populations, the power of the three point nulls from the closed form, from the approximation of `usdt_power()` and from 1000 simulated studies analysed with `hsdt()`. Population A is a plausible case of unconscious processing, with reliabilities from 0.4 to 0.85. The Vadillo population has an awareness reliability of 0.13 with 64 trials. ```{r} #| label: validation #| echo: false #| fig.height: 6.2 validation <- utils::read.csv(system.file("extdata", "power-validation.csv", package = "uSDT")) validation$design <- sprintf("%s subjects, %s + %s trials", validation$subjects, validation$trials_awareness, validation$trials_influence) validation$design <- factor(validation$design, levels = rev(unique(validation$design))) validation$population <- factor(validation$population, levels = c("A", "Vadillo")) validation$hypothesis <- factor( validation$hypothesis, levels = c("diff = 0", "rho = 0", "intercept = 0"), labels = c("H1: difference = 0", "H2: correlation = 0", "H3: intercept = 0")) long <- rbind( data.frame(validation[c("population", "design", "hypothesis")], engine = "Closed form", power = validation$closed), data.frame(validation[c("population", "design", "hypothesis")], engine = "Approximation of usdt_power()", power = validation$approximate), data.frame(validation[c("population", "design", "hypothesis")], engine = "Simulation with hsdt()", power = validation$simulation)) long$engine <- factor(long$engine, levels = c("Simulation with hsdt()", "Approximation of usdt_power()", "Closed form")) ggplot2::ggplot(long, ggplot2::aes(x = power, y = design)) + ggplot2::geom_errorbar( data = validation, ggplot2::aes(xmin = simulation_low, xmax = simulation_high, y = design), inherit.aes = FALSE, orientation = "y", width = 0.45, colour = "#292D30", linewidth = 0.5) + ggplot2::geom_point(ggplot2::aes(shape = engine, colour = engine, fill = engine), size = 2.6, stroke = 0.9) + ggplot2::scale_shape_manual(values = c(21, 24, 23), name = NULL) + ggplot2::scale_colour_manual(values = c("#292D30", "#116B60", "#B6543A"), name = NULL) + ggplot2::scale_fill_manual(values = c("#292D30", "#3FA88E", "white"), name = NULL) + ggplot2::scale_x_continuous(limits = c(0, 1), breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) + ggplot2::facet_grid(population ~ hypothesis, scales = "free_y", space = "free_y") + ggplot2::labs( x = "Power", y = NULL, caption = paste0( "1000 simulated studies per design. The bars are 95% Wilson intervals ", "of the simulated power. Studies whose fit ends on the boundary count ", "as failures to reject the correlation and the intercept.")) + ggplot2::theme_classic(base_size = 12) + ggplot2::theme( legend.position = "top", panel.spacing.x = grid::unit(14, "pt"), strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_text(colour = "#252A2D"), panel.grid.major.x = ggplot2::element_line(colour = "#E4E5E1", linewidth = 0.4), axis.text = ggplot2::element_text(colour = "#3E4347"), plot.caption = ggplot2::element_text(colour = "#687076", size = 9.5, hjust = 0)) ``` For these designs, H1 agrees across methods. For H2 and H3, the closed form diverges when awareness reliability is low because it ignores boundary outcomes. The sampling-distribution approximation follows the simulation more closely, but remains somewhat optimistic for the Vadillo population. Check chosen designs by simulation, particularly with low reliability. The approximation assumes equal signal/noise variances. To check sensitivity to unequal variances, set `sd_signal` in `usdt_population()` and use `method = "simulation"`. Those studies are still analysed with the **equal-variance** frequentist model. ## Summary | Question | Function | |---|---| | What population do I assume? | `usdt_population()`, from values or from a pilot fit | | What power does my design have? | `usdt_power()` | | How many subjects, trials, or both? | `usdt_plan()`, with `cost`, `budget` or `precision` | | What if my assumptions are wrong? | several values in `usdt_population()` or `update()` | | What is the smallest effect I can detect? | `usdt_mdes()` | | Does the chosen design hold up? | `usdt_power(plan, method = "simulation")` and `summary()` | ## References Vadillo, M. A., Malejka, S., & Shanks, D. R. (2025). Mapping the reliability multiverse of contextual cuing. *Journal of Experimental Psychology: Learning, Memory, and Cognition*, 51(6), 910-927. ```{r} #| label: cleanup #| include: false options(old_options) invisible(suppressWarnings(Sys.setlocale("LC_CTYPE", old_ctype))) ```