---
title: "Matching, dimensions and designs"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Matching, dimensions and designs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = FALSE,
  comment = "",
  R.options = list(
    cli.num_colors = 1,
    cli.hyperlink = FALSE,
    crayon.enabled = FALSE,
    width = 80
  )
)
# 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.
```

Choosing stimuli is a constrained selection problem. One lexical property is
manipulated, several others must not vary with it, and the words that satisfy
every constraint at once are a small and awkward subset of the language. This
vignette covers how `lexsync` states those constraints and how it resolves them:
the dimensions it computes, the filters that narrow a lexicon to a candidate
pool, the four matching methods and the situations each one suits, the tolerance
windows that bound an acceptable match, the continuous alternative to a matched
dichotomy, and the report that tells you what control the selection actually
achieved.

Everything below runs on the English example lexicon bundled with the package, a
3000-word slice that needs no download.

```{r setup}
library(lexsync)
schema <- yaml::read_yaml(
  system.file("extdata", "schema.yaml", package = "lexsync")
)
lex <- load_lexicon(
  system.file("extdata", "en_example.csv", package = "lexsync"),
  schema, language = "english"
)
head(
  lex[, c("word", "frequency", "length", "n_syllables", "n_density", "old20")],
  4
)
```

## The lexical dimensions

A dimension is any numeric column the matcher can equate. `lexsync` distinguishes
those read straight from the lexicon file from those computed from the
orthographic form at load time, and the schema records the distinction along with
each dimension's unit. Units matter when reading a report: a raw difference of 1
means one letter on `length` but a factor of ten in corpus occurrence on
`frequency`, which is why the comparisons are standardised.

```{r dimensions}
dims <- schema$dimensions
knitr::kable(data.frame(
  dimension = names(dims),
  type = vapply(dims, function(d) d$type, character(1)),
  unit = vapply(dims, function(d) d$unit, character(1)),
  row.names = NULL
), caption = "The dimensions declared in schema.yaml")
```

`length` is the character count of the orthographic form. `frequency` is the Zipf
value read from the lexicon's `freq_zipf` column, a log-scaled measure on which
roughly 1 is a word occurring once per hundred million tokens and 7 a word
occurring ten thousand times per million. The value is the base-10 logarithm of
the frequency per billion words (van Heuven et al., 2014). It is
preferred to raw counts because it is linear in the log space in which frequency
effects are approximately linear, and because it is comparable across corpora of
different sizes. `n_density` is Coltheart's N, the count of same-length words
differing by one letter substitution (Coltheart et al., 1977). `old20` is the
mean Levenshtein distance to the twenty nearest words in the reference set, which
unlike Coltheart's N is defined for long words, where substitution neighbours run
out (Yarkoni et al., 2008). A larger OLD20 means a more isolated word, so it runs
in the opposite direction to `n_density`.

`n_syllables` also arrives at load time, counting maximal runs of vowels, an
orthographic estimate with no pronunciation model behind it. Only
`bigram_freq` is computed on request, by `add_bigram_frequency()`: it is the mean
corpus probability of a word's adjacent letter bigrams, standing proxy for
phonotactic probability.

```{r syllables}
count_syllables(c("cat", "table", "rhythm", "being", "chocolate"))
```

The estimate is worth its caveat, and the errors run in both directions. It
handles `rhythm` because `y` is treated as a vowel letter. It undercounts
`being`, whose `ei` is one run of vowel letters spanning two nuclei, and it
overcounts `chocolate`, whose silent final `e` forms a run of its own. The
measure is adequate for equating two conditions on gross syllabic complexity and
it is not a substitute for a pronunciation dictionary.

### Computing dimensions on a pool

The bundled lexicon arrives with `n_density` and `old20` already present, but
neighbourhood measures depend on the reference set they are computed against, so
a lexicon derived elsewhere may need them recomputed. `add_neighbourhood()` does
that. Both measures are quadratic in the size of the reference set, which is why
the pipeline computes them on the pool alone, and why the reference should still
be the whole lexicon: a word's neighbours are its neighbours in the language,
not among the handful of words that survived your filters.

```{r neighbourhood}
sample_words <- lex[lex$word %in% c("cat", "dog", "house", "rhythm"), ]
add_neighbourhood(
  sample_words, reference = lex$word
)[, c("word", "n_density", "old20")]
```

`add_bigram_frequency()` works the same way, and `merge_norms()` joins any
word-keyed norm table, which is how dimensions the corpus does not carry
(concreteness, age of acquisition, behavioural measures from a megastudy) become
matchable. Once a column is present, the matcher treats it like any other.

```{r norms}
aoa <- data.frame(word = c("cat", "dog", "house"), aoa = c(3.1, 3.0, 3.4))
merged <- merge_norms(lex, aoa)
merged[merged$word %in% aoa$word, c("word", "frequency", "aoa")]
```

The join is case- and whitespace-insensitive on both sides, keeps the lexicon's
own row and column order, and refuses a norm column whose name already exists on
the lexicon. A silent rename would be the worst available outcome: it leaves the
dimension the design matches on under a name nothing looks for.

In a design you do not call this yourself. Name the tables in a `norms:` block and
the pipeline joins them before the pool is built, so a norm column can be filtered
on, matched on or spanned like any other dimension.

```yaml
norms:
  - path: norms/en_concreteness.csv
    on: word                  # optional join key, default 'word'
    columns: [concreteness]   # optional; default every column but the key

pool_filters:
  concreteness: [2.0, 5.0]
match_on: [length, concreteness]
```

lexsync ships no norm data: licences vary, and the citation is yours to honour.
What it does do is record what you joined. Every table appears in the materials
datasheet with its checksum and its per-column coverage, because a norm table can
supply the very variable a design manipulates, and a selection over columns of
unstated origin is not reproducible from the record that exists to make it so.
Coverage is recorded for a related reason: a word the table does not cover gets a
missing value, and the tolerance windows then drop it from the pool.

## Pool filters

A pool filter narrows the lexicon before any matching happens. It states what
counts as an eligible word at all, independently of condition, and it is the
right place for constraints you would apply to every item in the study.

```{r pool}
pool_filters <- list(length = c(3, 7), frequency = c(3.8, 7.0))
pool <- build_pool(lex, pool_filters)
c(lexicon = nrow(lex), pool = nrow(pool))
```

A numeric pair is read as an inclusive range. Any other vector is read as a set of
permitted values, so `list(pos = c("noun", "verb"))` would keep nouns and verbs
if the lexicon carried a `pos` column. A filter naming a column the lexicon does
not have is skipped silently, which keeps one schema usable across corpora that
differ in their optional columns, and which is worth knowing about when a filter
appears not to bite. Rows missing a value on a filtered column are always
dropped.

Filtering is a separate step from matching. `match_stimuli()` reads a design's
`conditions` and `match_on` fields but never its `pool_filters`, so passing a raw
lexicon where a pool was intended silently widens the selection, with no error to
warn you. The pipeline calls `build_pool()` first and the examples here do the
same.

## Supplied item pools

Narrowing a whole lexicon is the right way round when any word of the language
will do. It is the wrong way round when the candidate set is itself a research
decision: a list from a previous study, from a norming session, or one restricted
to a semantic category that no lexical filter can express. Such a list used to
have to masquerade as a lexicon, which meant inventing a `freq_zipf` column for
it, or else skip matching entirely by going in as `items.source: table`.

`items.source: pool` takes the list as it is.

```yaml
items:
  source: pool
  path: items/pool_en_concrete_nouns.csv
  lexicon: corpora/derived/en.csv
```

The list needs only a `word` column. Length and the syllable estimate are derived
from the form, and `items.lexicon` supplies the rest, frequency above all, by
lookup. A word the lexicon does not have raises a hard error, because the
tolerance windows drop rows with missing values silently and the pool would then
be smaller than it appears to be.

One subtlety would quietly invalidate the controls if it went the other way.
`n_density` and `old20` are properties of a word in its *language*, not among the
131 words of a supplied list, so they are computed against the lexicon's words. A
supplied pool with no lexicon falls back to itself, which is why one is given
above. `load_pool()` returns both parts of that, the pool and the reference word
list, if you call it directly outside a design.

Everything downstream is unchanged: the same conditions, the same `match_on`, the
same report, the same datasheet. What differs is only where the candidates came
from, and the datasheet records the list's checksum so that stays answerable
later. `design_en_supplied_pool.yaml` is a worked example.

## Relational designs, where the item is a pair

A priming study's theoretical variable is a property of the *pair*, such as
distributional similarity or associative strength, while its nuisance variables
are properties of each *member*: frequency, length, neighbourhood. A design could
previously have matching, or its own item table, but not both, so a relational
design had to be assembled by hand.

`items.members` promotes an ordinary item table to a pair-keyed one.

```yaml
items:
  source: table
  path: items/priming_pairs_en.csv
  members: [prime, target]          # the two word fields
  lexicon: corpora/derived/en.csv   # where each member's norms are looked up
  anchor_condition: related         # whose row represents the pair during selection
```

Each member's word-level norms are joined from the lexicon and prefixed, giving
`prime.frequency`, `target.length` and so on. The member name leads and the
dimension follows, for a reason specific to R: `prime.frequency` is safe because
`df$prime` still exact-matches the bare `prime` column, whereas
`frequency.prime` would be a partial-matching hazard. Those prefixed names are
usable wherever a dimension is: in `pool_filters`, in `match_on`, as a
`continuous` predictor or as a control.

`pair.lev` is the Levenshtein distance between the two members and `pair.overlap`
is `1 - pair.lev / length of the longer form`, a normalised orthographic
similarity. It scores edit operations, so it is not a count of shared material:
`abc` and `cba` share every letter but score 0.33. Orthographic overlap is the
standard confound control in a priming design, since a related pair that also
shares letters confounds semantic relatedness with orthographic similarity. The
pipeline adds both automatically when the two members are named `prime` and
`target`. Under any other member names, call `add_pair_overlap()` on the item
table yourself. A predictor that cannot be computed from the forms alone, such as
cosine similarity or forward associative strength, arrives as an extra column on
the item table and is used like any other.

```yaml
n_per_condition: 8                  # counts PAIRS: 8 pairs, 16 rows, 8 per list
pool_filters:
  target.length: [3, 7]
  target.frequency: [3.0, 6.0]
continuous:
  predictor: target.frequency
  controls: [target.length, pair.overlap]
match_on: [target.length, pair.overlap]
```

Two rules follow from the item being a pair. A filter applies to a pair only if
*every* one of its rows passes, and selection runs over one row per pair, the
`anchor_condition`, defaulting to the byte-first condition, with the result
re-expanded afterwards so every condition row of every chosen pair survives.
Neither is an optimisation. A filter on `prime.frequency` applied row by row would
keep a pair's related row and drop its unrelated one, leaving a set the
Latin-square counterbalancer cannot complete. Only a dimension that varies within
the pair can split it that way, which the prime-level and pair-level ones do and
the target-level ones do not, since the target is the same word in both rows.

A member name may not be one of `word`, `id`, `set`, `list`, `trial`, `condition`,
`replicate` or `item`, which the engines read directly, and a design using one is
refused. The datasheet gains a `relational` block recording the members, the pair
count, the member lexicon and its checksum, and the member-level dimensions
separately from the relational ones. `design_en_priming_continuous.yaml` is the
worked example, and its own comments go further into why each rule is there.

## Conditions, the anchor and matched sets

A condition is a named region of the pool, given by `define_by`, which is read
exactly as `pool_filters` is. The first condition listed is the anchor. Its items
are chosen by an even spread across the sorted candidate subpool, so the anchor
covers the whole of its region on the manipulated dimension and does not cluster
at one end. Every later condition is then matched to the anchor item by item, and
the `set` column records the pairing.

```{r design}
design <- list(
  name = "vignette_matching", language = "english", n_per_condition = 20,
  pool_filters = pool_filters,
  conditions = list(
    list(name = "high", define_by = list(frequency = c(5.2, 7.0))),
    list(name = "low",  define_by = list(frequency = c(3.8, 4.4)))
  ),
  match_on = list("length", "n_density", "old20")
)
stim <- match_stimuli(pool, design, schema)
head(
  stim[
    stim$set %in% 1:3,
    c("set", "condition", "word", "frequency", "length", "old20")
  ],
  6
)
```

Because the anchor is spread across its own subpool, the order of the conditions
is a design decision with consequences. The anchor's distribution on the
manipulated dimension is inherited from the pool, and every other condition's
distribution is whatever matching to that anchor produces.

## The four matching methods

The method is set with `matching.method` in a design, or left to the schema
default. All four minimise a distance in standardised space, computed by
z-scoring each `match_on` dimension against the whole pool, so that a letter and a
Zipf unit become comparable. They differ in what they minimise it over.

`standardised_euclidean`, the default, walks the anchor items in order and gives
each one its nearest unused counterpart, greedily. Because the anchor is selected
before any matching happens, it preserves the anchor's even spread exactly. Being
anchored, it generalises past two conditions, as `mahalanobis` does, while `joint`
and `optimal` score cross-condition pairs and are refused with any number of
conditions but two.

`joint` abandons the anchor. It scores every cross-condition pair, then takes the
best-matched pairs in ascending order of distance until it has enough, so an item
with no good counterpart is simply never selected. This equates the controls more
tightly than per-anchor matching when the manipulation is confounded with a
control, and the price is that the manipulated dimension's realised distribution
is no longer under your control. It requires exactly two conditions.

`mahalanobis` replaces the Euclidean metric with one weighted by the inverse of
the pool's correlation matrix, so that dimensions sharing variance are not
counted twice. It suits the case where the matched dimensions are strongly
intercorrelated, which is the normal case for orthographic measures.

```{r correlations}
round(cor(pool[, c("length", "n_density", "old20")]), 2)
```

Length, neighbourhood density and OLD20 are far from independent here. Under a
plain Euclidean metric, matching on all three counts the shared
length-neighbourhood variance three times over and lets an independent dimension
drift. Mahalanobis distance discounts that overlap (Rubin, 1980; Stuart, 2010).

`optimal` keeps the Euclidean metric but solves the pairing globally, minimising
the summed distance over all pairings at once (Gu & Rosenbaum, 1993; Hansen &
Klopfer, 2006). Greedy matching can spend the one good counterpart of a hard
item on an easy item that had alternatives. Optimal matching cannot. It requires
exactly two conditions and the `clue` package.

```{r methods}
methods <- c("standardised_euclidean", "joint", "mahalanobis")
if (requireNamespace("clue", quietly = TRUE)) methods <- c(methods, "optimal")

summarise_method <- function(method) {
  d <- design
  d$matching <- list(method = method)
  s <- match_stimuli(pool, d, schema)
  cmp <- match_report(
    s, c("frequency", "length", "n_density", "old20"), schema
  )$comparisons
  controls <- cmp[cmp$dimension != "frequency", ]
  high <- s$frequency[s$condition == "high"]
  low  <- s$frequency[s$condition == "low"]
  data.frame(
    method = method,
    worst_control_d = round(max(abs(controls$cohens_d)), 3),
    freq_d = round(cmp$cohens_d[cmp$dimension == "frequency"], 2),
    freq_raw_diff = round(mean(high) - mean(low), 2),
    freq_sd_high = round(sd(high), 2)
  )
}
do.call(rbind, lapply(methods, summarise_method))
```

The comparison rewards a careful reading, and it is a good illustration of why a
standardised difference should not be read as an effect size on its own. The two
globally pairing methods drive the control differences to zero, because a
free choice over all pairs can find counterparts that match exactly on all three
dimensions. Their standardised frequency separation is also the larger of the
four, yet their raw separation in Zipf units is the smaller. Both facts have one
cause: those methods take whichever pairs match best, and those pairs turn out to
cluster, which shrinks the SD of the high condition and so shrinks the pooled SD
that `freq_d` divides by. A narrower manipulation can raise `d`.

Whether that is a gain depends on the study. The anchor methods keep the high
condition spread across its whole band by construction, which supports a claim
about the band rather than about a cluster inside it, and it is variance you may
want if the manipulation is later modelled as continuous. The globally pairing
methods trade some of that spread for exact control of the nuisance dimensions.
These numbers describe this pool and this design. The point is how to make the
comparison, not which method wins.

### The parity caveat on two methods

The R and Python engines select byte-identical stimuli under
`standardised_euclidean` and `joint`. They do not guarantee it under
`mahalanobis` or `optimal`. Both of the latter route the distance through an
operation whose last bits are decided by the platform's linear-algebra
implementation, a covariance-matrix inverse in one case and an assignment solver
in the other, and once two candidate distances differ in their final bits the
tie-break that would otherwise make the engines agree never fires. `mahalanobis`
usually still agrees exactly, and `optimal` often selects an equally-optimal but
different set. Each run's datasheet records which case applies, under
`selection.cross_engine`. If cross-engine reproducibility is a requirement of
your project, the two default methods are the ones that carry the guarantee. The
reproducibility vignette sets out why the exception cannot simply be engineered
away.

## Tolerance windows

Matching by distance alone will always return the nearest available candidate,
however far away it is. The tolerance window is the guard against that. Before
the distance ranking, each later condition's candidates are filtered to those
lying within the anchor's mean plus or minus *k* standard deviations on every
matched dimension, with *k* set per dimension.

```{r tolerance}
str(schema$matching$tolerance_k)
```

Smaller *k* is stricter. The default of 2.0 on the orthographic dimensions is
permissive, admitting most of a normal distribution, and it is meant as a backstop
against a wild counterpart. The distance ranking within the window is what does
the work of control. A design tightens the window per dimension, and only the
named dimensions are overridden.

The window belongs to the anchored path, so it is what `standardised_euclidean`
and `mahalanobis` apply. `joint` and `optimal` score cross-condition pairs
directly and never read `tolerance_k` at all, which is worth knowing before
tightening a window and finding the selection unmoved.

```{r tolerance-override}
tight <- design
tight$matching <- list(tolerance_k = list(old20 = 0.25))
tight_stim <- match_stimuli(pool, tight, schema, verbose = TRUE)
round(match_report(tight_stim, "old20", schema)$comparisons$cohens_d, 3)
```

There is a behaviour here that will bite you if you do not know it. If a window
leaves a condition with fewer candidates than there are anchor items, the matcher
neither fails nor returns a short set. It relaxes that condition's window and
carries on. Tightening the window past what the pool can supply therefore looks
exactly like a window that was never set, and under the default policy
`verbose = TRUE` is the only thing that tells you which happened.

```{r tolerance-relaxed}
too_tight <- design
too_tight$matching <- list(tolerance_k = list(old20 = 0.05))
relaxed_stim <- match_stimuli(pool, too_tight, schema, verbose = TRUE)
```

Two details of the relaxation matter. It is per condition and it is
all-or-nothing: the fallback drops that condition's windows on *every* matched
dimension at once, not merely the one that proved unsatisfiable, so an
over-tightened `old20` also discards the `length` and `n_density` windows for
that condition. And it does not degrade the match as much as you might fear,
because the distance ranking still runs over the whole subpool and still prefers
near counterparts. Only when the relaxed subpool is still too small does the
matcher stop, with an error naming the condition and the counts.

Both behaviours are policies, and the schema's `matching` block names them.
`on_insufficient_tolerance` governs the fallback just described. Its default of
`relax` keeps the behaviour above, while `error` makes the matcher refuse and say
so, which is the honest setting when the window is itself the criterion,
as it is when reproducing a published study's stated windows. `shortfall`
governs what happens when fewer sets than `n_per_condition` can be selected at
all, and it defaults to `error` deliberately: the datasheet and the generated
methods paragraph state the requested n, so a set that silently shrank would
leave the record misstating the materials. A design that can live with fewer
sets says so with `shortfall: allow`, and either key may be overridden under
the design's own `matching:` block, exactly as `tolerance_k` was overridden
above.

Windows are the place to reproduce a published study's stated criteria. The
bundled `design_es_gender_repro.yaml` sets `frequency` to `0.111...`, which is
the mean plus or minus SD/9 used by González Alonso et al. (2025), alongside
`length` at the mean plus or minus 2 SD.

One asymmetry is worth knowing. Tolerance windows are derived from the anchor, so
they apply to the anchor-based methods, `standardised_euclidean` and
`mahalanobis`, and to continuous designs. `joint` and `optimal` never build an
anchor and never consult `tolerance_k`. Both instead cap each subpool to the 1200
candidates nearest the other condition's centroid before pairing, which bounds
the cost of scoring every pair.

## Continuous designs

A matched dichotomy is not the only way to study a graded property, and for a
property like frequency it is often not the best one. Splitting a continuum into
a high and a low band discards the variation within each band, spends statistical
power to do it, and leaves the estimate dependent on where the split was placed.
The selection can also manufacture its own artefacts: Liben-Nowell et al. (2019)
show that choosing a controlled subset from a large candidate set can produce an
apparent effect from a set where none exists, because the selection is itself a
search over the space of possible item sets. Kuperman (2015) makes the
constructive case, that a megastudy's variation can be used directly by sampling
across a predictor's range and modelling it.

A `continuous` block does that. Rather than conditions, it names one predictor to
span and the controls to hold still, and the two must agree with `match_on`.

```{r continuous}
cont_design <- list(
  name = "vignette_continuous", language = "english", n_per_condition = 60,
  pool_filters = pool_filters,
  continuous = list(
    predictor = "frequency",
    controls = c("length", "n_density", "old20")
  ),
  match_on = list("length", "n_density", "old20"),
  matching = list(
    tolerance_k = list(length = 1.5, n_density = 1.5, old20 = 1.5)
  )
)
cont <- select_continuous_stimuli(pool, cont_design, schema)
nrow(cont)
range(cont$frequency)
```

Selection runs in two passes. An even spread over the pool defines the control
windows from its own distribution, the pool is filtered to those windows, and a
second even spread over what survives is the selection. The result covers the
predictor's range while the controls stay near-constant, and therefore
near-uncorrelated with the predictor, which is what the regression needs.

```{r continuous-report}
cont_report <- match_report_continuous(
  cont, "frequency",
  c("length", "n_density", "old20"), schema
)
cont_report$comparisons
```

The report is read differently from a matched one. `predictor_span` is the
realised range of the predictor in the selected set and should be wide, since it
is the variation the regression has to work with. Each control's `pearson_r` is its correlation
with the predictor and should be near zero, since a control that still correlates
with the predictor is still a confound. Tightening `tolerance_k` pushes those
correlations down at the cost of a smaller eligible pool, and that is the dial to
turn when a control will not come loose.

The datasheet for a continuous design reports the span, the correlations and a
suggested regression model in which the controls enter as covariates.

## Reading the validation output

`match_report()` returns two data frames. The descriptives give per-condition n,
mean, SD, minimum, median and maximum for each dimension, which is the table that
belongs in a paper's materials section.

```{r report-desc}
report <- match_report(
  stim, c("frequency", "length", "n_density", "old20"), schema
)
knitr::kable(
  report$descriptives,
  caption = "Descriptive statistics per condition"
)
```

The comparisons contrast every other condition against the anchor on every
dimension, and they carry the four statistics that answer whether the matching
worked.

```{r report-comp}
knitr::kable(
  report$comparisons,
  caption = paste(
    "Realised control: standardised differences with 90% intervals,",
    "variance ratios and TOST verdicts"
  )
)
```

`cohens_d` is the standardised mean difference against the anchor, using the
pooled SD. On the manipulated dimension it should be large and on a controlled
dimension it should be near zero.

`d_ci_low` and `d_ci_high` bound it with a 90% confidence interval. The interval,
not the point estimate, is the primary summary of a controlled dimension. Its
upper limit is the largest imbalance still consistent with the stimuli you have,
and with few items that limit can be substantial even when the point estimate is
almost zero. Reporting the interval keeps the dependence on set size visible,
where a bare 'no significant difference' hides it. The width is a property of the
number of items, so a small set cannot demonstrate tight control however good its
point estimates look (Sassenhagen & Alday, 2016).

`var_ratio` is the condition's variance over the anchor's. Two conditions can
share a mean and still confound if one is far more spread out than the other,
which every mean-based statistic misses, so this is the distributional check. A
ratio near 1 is balanced and a common heuristic flags ratios outside roughly 0.5
to 2. See `?variance_ratio` for the sources behind that heuristic.

`tost_p` and `equivalent` are the two one-sided tests. A conventional
non-significant difference test is not evidence of equivalence, because failing
to reject a null is not the same as supporting it, and the failure is easiest to
achieve with few items and noisy measures. TOST inverts this by testing against
an explicit bound: the reported p is the larger of two one-sided p-values against
a difference of plus or minus `bound_d` pooled SDs, and a value below alpha
supports equivalence within that bound (Lakens, 2017). The bound and alpha are
the schema's, so the verdict is only as meaningful as the bound you chose.

```{r equivalence-settings}
str(schema$equivalence)
```

A bound of 0.5 SD is the schema default and it is a convention rather than a
finding. It is a permissive bound for a nuisance dimension in a design where that
dimension has a strong effect on the response, and a study that cares about a
particular confound should set a bound it can defend. The 90% interval is the
1 - 2 x alpha interval precisely so that it corresponds to the TOST decision at
alpha: the test declares equivalence exactly when the interval falls entirely
inside the bounds, which means the interval and the verdict are two readings of
one thing.

The individual statistics are exported and can be called directly on any two
vectors.

```{r direct-stats}
high <- stim$old20[stim$condition == "high"]
low  <- stim$old20[stim$condition == "low"]
str(cohens_d_ci(high, low))
str(tost_equiv(high, low, bound_d = 0.5))
round(variance_ratio(low, high), 3)
```

`balance_check()` is the last piece and answers a different question: whether the
levels of a column occur equally often. It returns human-readable warnings and an
empty vector when all is well, and the pipeline logs whatever it returns.

```{r balance}
balance_check(stim, "condition")
balance_check(rbind(stim, stim[1, ]), "condition")
```

## Pseudowords

Lexical decision needs non-words, and the difficulty is that they must be
orthographically plausible without being real. `build_lexdec_stimuli()` draws real
words by an even spread across the byte-ordered pool, then generates a
length-matched pseudoword for each.

```{r lexdec}
lexdec_pool <- build_pool(lex, list(length = c(4, 7), frequency = c(3.5, 6.0)))
stim_lexdec <- build_lexdec_stimuli(
  lexdec_pool, n = 8, reference_words = lex$word
)
head(stim_lexdec[, c("target", "condition", "length", "set")], 4)
```

The presented string is `target`, the conditions are `word` and `pseudoword`, and
`set` pairs each word with its non-word twin. `reference_words` should be the full
lexicon: it supplies the bigram statistics and the list of real forms a pseudoword
must avoid, and it defaults to the pool when omitted, which is rarely what you
want.

Two methods are available, chosen per design under `items.generation.method`.

`letter_substitution`, the default, changes as few letters as possible subject to
three constraints: every resulting bigram must be attested in the corpus, the form
must not be a real word or one already generated, and the length must be preserved
exactly. Single substitutions are searched first, then pairs as a fallback. Among
the legal candidates the most bigram-plausible wins, with a byte-order tie-break.

`subsyllabic` splits each word into onset, nucleus and coda constituents and swaps
whole constituents for attested alternatives of the same role and length, so the
pseudowords keep their syllabic structure. Codas and nuclei are changed before
onsets, which carry the most identifying orthography. Roughly two thirds of a
word's constituents are targeted, each swap preserves length, and a word with no
legal swap falls back to letter substitution, so every word yields a pseudoword.
It is a deterministic orthographic approximation of Wuggy (Keuleers & Brysbaert,
2010), trading Wuggy's phonological model for exact length matching and
cross-engine reproducibility.

```{r subsyllabic}
sub_stim <- build_lexdec_stimuli(
  lexdec_pool, n = 8, reference_words = lex$word, method = "subsyllabic"
)
# The real words come first in the frame, so pair each with its twin by `set`.
merge(
  sub_stim[sub_stim$condition == "word", c("set", "target")],
  sub_stim[sub_stim$condition == "pseudoword", c("set", "target")],
  by = "set", suffixes = c("_word", "_pseudoword")
)[1:4, ]
```

Both generators process base words in byte order, so the set of already-used
pseudowords evolves identically in the two engines and both select byte-identical
stimuli. Both are also orthographic models defined for a-z words:
`build_lexdec_stimuli()` filters the pool to those, and the subsyllabic segmenter
returns nothing for a word with an accent, a hyphen or a digit, which sends it
down the letter-substitution path. That is why lexical decision is demonstrated on
the English designs and never on the Chinese one.

`generate_pseudowords()` and `make_pseudoword()` are available directly if you
want the forms without the lexical-decision scaffolding around them.

```{r pseudowords}
generate_pseudowords(c("house", "table"), lex$word)
```

In a design this is `items.source: generate`, with the method under
`items.generation.method`. `design_en_lexdec.yaml` and `design_en_lexdec_wuggy.yaml`
are the worked examples of the two methods, and `design_en_lexdec_blocks.yaml`
adds practice and filler blocks to the same generated set.

## Items as a random factor

A matched set is one sample from the language, and an analysis that treats its
items as fixed generalises to those items rather than to the language they were
drawn from (Clark, 1973; Yarkoni, 2022). Crossed random effects are the standard
answer, and the datasheet suggests a model that includes them. `resample_stimuli()`
supports the stronger move of drawing several fully matched sets that share no
items, each one an independent sample, which lets a study run different item
samples across participant groups or show that an effect survives the change of
materials.

```{r resample}
reps <- resample_stimuli(pool, design, schema, n_sets = 2)
table(reps$replicate, reps$condition)
length(intersect(
  reps$word[reps$replicate == 1],
  reps$word[reps$replicate == 2]
))
```

Each replicate is matched independently against the pool with the earlier
replicates' items removed, so the sets are disjoint by construction. Later
replicates draw from a depleted pool, which means match quality can degrade
across replicates. The report will show it if it does.

## References

Clark, H. H. (1973). The language-as-fixed-effect fallacy: A critique of language
statistics in psychological research. *Journal of Verbal Learning and Verbal
Behavior*, *12*(4), 335–359. https://doi.org/10.1016/S0022-5371(73)80014-3

Coltheart, M., Davelaar, E., Jonasson, J. T., & Besner, D. (1977). Access to the
internal lexicon. In S. Dornic (Ed.), *Attention and Performance VI*
(pp. 535–555). Erlbaum.

González Alonso, J., Bernabeu, P., Silva, G., DeLuca, V., Poch, C., Ivanova, I., &
Rothman, J. (2025). Starting from the very beginning: Unraveling third language
(L3) development with longitudinal data from artificial language learning and
EEG. *International Journal of Multilingualism*, *22*(1), 119–142.
https://doi.org/10.1080/14790718.2024.2415993

Gu, X. S., & Rosenbaum, P. R. (1993). Comparison of multivariate matching
methods: Structures, distances, and algorithms. *Journal of Computational and
Graphical Statistics*, *2*(4), 405–420.
https://doi.org/10.1080/10618600.1993.10474623

Hansen, B. B., & Klopfer, S. O. (2006). Optimal full matching and related designs
via network flows. *Journal of Computational and Graphical Statistics*, *15*(3),
609–627. https://doi.org/10.1198/106186006X137047

Keuleers, E., & Brysbaert, M. (2010). Wuggy: A multilingual pseudoword
generator. *Behavior Research Methods*, *42*(3), 627–633.
https://doi.org/10.3758/BRM.42.3.627

Kuperman, V. (2015). Virtual experiments in megastudies: A case study of language
and emotion. *Quarterly Journal of Experimental Psychology*, *68*(8), 1693–1710.
https://doi.org/10.1080/17470218.2014.989865

Lakens, D. (2017). Equivalence tests: A practical primer for *t* tests,
correlations, and meta-analyses. *Social Psychological and Personality Science*,
*8*(4), 355–362. https://doi.org/10.1177/1948550617697177

Liben-Nowell, D., Strand, J., Sharp, A., Wexler, T., & Woods, K. (2019). The
danger of testing by selecting controlled subsets, with applications to
spoken-word recognition. *Journal of Cognition*, *2*(1), Article 2.
https://doi.org/10.5334/joc.51

Rubin, D. B. (1980). Bias reduction using Mahalanobis-metric matching.
*Biometrics*, *36*(2), 293–298. https://doi.org/10.2307/2529981

Sassenhagen, J., & Alday, P. M. (2016). A common misapplication of statistical
inference: Nuisance control with null-hypothesis significance tests. *Brain and
Language*, *162*, 42–45. https://doi.org/10.1016/j.bandl.2016.08.001

Stuart, E. A. (2010). Matching methods for causal inference: A review and a look
forward. *Statistical Science*, *25*(1), 1–21. https://doi.org/10.1214/09-STS313

van Heuven, W. J. B., Mandera, P., Keuleers, E., & Brysbaert, M. (2014).
SUBTLEX-UK: A new and improved word frequency database for British English.
*Quarterly Journal of Experimental Psychology*, *67*(6), 1176–1190.
https://doi.org/10.1080/17470218.2013.850521

Yarkoni, T. (2022). The generalizability crisis. *Behavioral and Brain Sciences*,
*45*, Article e1. https://doi.org/10.1017/S0140525X20001685

Yarkoni, T., Balota, D., & Yap, M. (2008). Moving beyond Coltheart's N: A new
measure of orthographic similarity. *Psychonomic Bulletin & Review*, *15*(5),
971–979. https://doi.org/10.3758/PBR.15.5.971