Every response time analysis makes an exclusion decision before the
model is fitted, and most inherit it. The 200 ms and 3 s cutoffs come
from the last paper, the ±2.5 SD criterion from the field. What that
decision does to the parameters is invisible in the data you have,
because the contaminants are not labelled. rtprep addresses
both halves of that problem. It gives the screening rules one interface,
so that alternatives can be compared instead of assumed, and it
generates data with the contaminants labelled, so that a pipeline can be
tested instead of trusted.
This vignette walks the whole chain once, inside a dplyr pipeline: screen, filter, aggregate, estimate, and then check the result against the truth.
The data
rt_example holds four participants, two conditions, and
100 trials per cell. It was generated by r_contaminated(),
so every trial carries two columns a real data set cannot: whether it
was a contaminant, and which process produced it.
head(rt_example)
#> id condition trial rt response contaminant process true_drift contam_rate
#> 1 p1 easy 1 0.590 1 FALSE clean 1.615 0.02
#> 2 p1 easy 2 0.408 1 FALSE clean 1.615 0.02
#> 3 p1 easy 3 1.334 1 FALSE clean 1.615 0.02
#> 4 p1 easy 4 0.563 1 FALSE clean 1.615 0.02
#> 5 p1 easy 5 0.574 0 FALSE clean 1.615 0.02
#> 6 p1 easy 6 0.665 1 FALSE clean 1.615 0.02
rt_example |>
count(id, contam_rate, contaminant) |>
filter(contaminant)
#> id contam_rate contaminant n
#> 1 p1 0.02 TRUE 6
#> 2 p2 0.05 TRUE 10
#> 3 p3 0.10 TRUE 18
#> 4 p4 0.15 TRUE 26The participants contaminate at different rates, from 2% to 15%. That is deliberate. The question a preprocessing pipeline has to answer is not only whether it removes contaminants on average, but whether the error it leaves behind depends on how much each participant contaminated. We come back to that at the end.
Screening
A rule is a constructor that holds parameters;
rt_screen() applies it. Whatever the rule, the result is
one row per trial with the same four columns: the keep decision, the
probability that the trial came from the decision process, the rule’s
label, and the reason for a flag.
screened <- rt_example |>
mutate(rt_screen(rt, rule_sd(2.5)), .by = c(id, condition))
head(screened)
#> id condition trial rt response contaminant process true_drift contam_rate
#> 1 p1 easy 1 0.590 1 FALSE clean 1.615 0.02
#> 2 p1 easy 2 0.408 1 FALSE clean 1.615 0.02
#> 3 p1 easy 3 1.334 1 FALSE clean 1.615 0.02
#> 4 p1 easy 4 0.563 1 FALSE clean 1.615 0.02
#> 5 p1 easy 5 0.574 0 FALSE clean 1.615 0.02
#> 6 p1 easy 6 0.665 1 FALSE clean 1.615 0.02
#> .keep .prob .rule .reason
#> 1 TRUE 1 sd(2.5, mean, sd) <NA>
#> 2 TRUE 1 sd(2.5, mean, sd) <NA>
#> 3 FALSE 0 sd(2.5, mean, sd) too_slow
#> 4 TRUE 1 sd(2.5, mean, sd) <NA>
#> 5 TRUE 1 sd(2.5, mean, sd) <NA>
#> 6 TRUE 1 sd(2.5, mean, sd) <NA>The unnamed mutate() splices the four columns in, and
.by fits the rule within each participant and condition.
The rule could have been anything else in the roster with no change to
the lines around it. Absolute cutoffs, the median absolute deviation
criterion, the recursive criteria of Van Selst and Jolicoeur (1994), and
the contaminant mixture of Ratcliff and Tuerlinckx (2002) all come back
in this shape:
rule_cutoff(0.18, 3)
#> <rtprep rule> cutoff(0.18, 3)
#> Exclude trials outside [0.18, 3] seconds (bounds inclusive).
rule_mad(2.5)
#> <rtprep rule> sd(2.5, median, mad)
#> Exclude trials more than 2.5 x mad from the median, computed per group.
rule_recursive("modified")
#> <rtprep rule> recursive(modified)
#> Van Selst & Jolicoeur (1994) modified recursive criterion, iterated.
rule_mixture("lognormal")
#> <rtprep rule> mixture(lognormal)
#> Flag trials by their posterior probability under a uniform-contaminant / lognormal mixture.Filtering is one more verb. rt_keep() returns the keep
column alone and says once how many trials it dropped, so the exclusion
count sits next to the exclusion rather than being reconstructed later.
Pass .by to rt_keep() rather than to
filter(), as a list of the grouping columns, because inside
filter() there is no tidy selection and
c(id, condition) would concatenate them. The keep vector is
the same either way, and the count then covers the whole data set in one
line.
clean <- rt_example |>
filter(rt_keep(rt, rule_sd(2.5), .by = list(id, condition)))
#> sd(2.5, mean, sd): dropped 26 of 800 trials (3.2%)
nrow(clean)
#> [1] 774The per-cell diagnostics, which mutate() would strip
from the attribute, come from screen_fits():
rt_example |>
reframe(screen_fits(rt, rule_sd(2.5)), .by = c(id, condition)) |>
select(id, condition, n_trials, n_dropped, prop_dropped, lower, upper)
#> id condition n_trials n_dropped prop_dropped lower upper
#> 1 p1 easy 100 4 0.04 0.001490163 1.259977
#> 2 p1 hard 100 2 0.02 0.052523638 1.097803
#> 3 p2 easy 100 2 0.02 0.088801179 1.049293
#> 4 p2 hard 100 4 0.04 -0.160110325 1.472380
#> 5 p3 easy 100 4 0.04 -0.273151826 1.521747
#> 6 p3 hard 100 2 0.02 -0.167324763 1.449914
#> 7 p4 easy 100 3 0.03 -0.359058484 1.611060
#> 8 p4 hard 100 5 0.05 -0.092933241 1.310280Did the rule remove what you think it removed?
With real data the story ends here. With generated data it does not, and this is the point of generating any. The rule’s decisions can be set against the truth:
screened |>
summarise(
sensitivity = mean(!.keep[contaminant]),
specificity = mean(.keep[!contaminant]),
.by = id
)
#> id sensitivity specificity
#> 1 p1 0.3333333 0.9793814
#> 2 p2 0.2000000 0.9789474
#> 3 p3 0.3333333 1.0000000
#> 4 p4 0.1923077 0.9827586The ±2.5 SD criterion finds roughly a fifth to a third of each participant’s contaminants and keeps about 98% of the genuine trials. Which contaminants it finds depends on where they sit relative to the clean response time distribution, and splitting the hits by process shows that directly:
screened |>
filter(contaminant) |>
summarise(found = mean(!.keep), n = n(), .by = process)
#> process found n
#> 1 informationless 0.0500000 20
#> 2 delay 0.7368421 19
#> 3 leading_edge 0.0000000 21The delayed start-ups sit in the upper tail, where a symmetric criterion around the mean reaches them. The anticipations sit at the leading edge, which in a right-skewed distribution is well inside 2.5 standard deviations of the mean, so the same criterion removes none of them. The informationless responses overlap the clean core by construction, and a rule that reads only response times catches almost none of them, here 1 of 20. A rule’s hit rate is a property of the contaminant’s location, not of the rule alone.
Aggregation and estimation
Screening decides which trials survive; aggregation decides what the
survivors are summarised as. rt_summary() returns the
inputs the EZ-diffusion equations need, and ez_ddm()
inverts them into drift, bound, and non-decision time. Both take vectors
and return one row, so they fit the same grammar:
estimates <- clean |>
reframe(rt_summary(rt, response), .by = c(id, condition)) |>
mutate(ez_ddm(mean_rt, var_rt, n_upper / n_trials, n_trials))
estimates |>
select(id, condition, n_trials, drift, bound, ndt)
#> id condition n_trials drift bound ndt
#> 1 p1 easy 96 1.603160 1.0519220 0.3686653
#> 2 p1 hard 98 1.149562 0.9795129 0.3410422
#> 3 p2 easy 98 1.801162 1.0424950 0.3421352
#> 4 p2 hard 96 1.666531 1.2269532 0.3235057
#> 5 p3 easy 96 1.983975 1.2086318 0.3088695
#> 6 p3 hard 98 1.430030 1.1963321 0.3193125
#> 7 p4 easy 97 1.820633 1.4027536 0.2418490
#> 8 p4 hard 95 1.851720 1.1557183 0.3145477rt_summary() has a method argument for the
robust (median and IQR) and mixture-based moments, which act on the same
trials without removing any; ?rt_summary documents them,
and the aggregation article on the package website puts the routes side
by side.
Does the error track the contamination rate?
The participants differ in how much they contaminated, and the data carry the drift each cell was generated from. Setting the two side by side shows what a screen leaves behind:
truth <- rt_example |>
distinct(id, condition, true_drift, contam_rate)
no_screen <- rt_example |>
reframe(rt_summary(rt, response), .by = c(id, condition)) |>
mutate(ez_ddm(mean_rt, var_rt, n_upper / n_trials, n_trials)) |>
select(id, condition, drift_none = drift)
estimates |>
select(id, condition, drift_sd = drift) |>
left_join(no_screen, by = c("id", "condition")) |>
left_join(truth, by = c("id", "condition")) |>
mutate(
error_none = drift_none - true_drift,
error_sd = drift_sd - true_drift
) |>
select(id, condition, contam_rate, true_drift, error_none, error_sd) |>
arrange(contam_rate, condition)
#> id condition contam_rate true_drift error_none error_sd
#> 1 p1 hard 0.02 1.275 -0.24942309 -0.125437560
#> 2 p1 easy 0.02 1.615 -0.24452286 -0.011840246
#> 3 p2 hard 0.05 1.425 -0.04425564 0.241531380
#> 4 p2 easy 0.05 1.805 -0.18088400 -0.003838382
#> 5 p3 hard 0.10 1.575 -0.40770719 -0.144970326
#> 6 p3 easy 0.10 1.995 -0.58442080 -0.011024995
#> 7 p4 hard 0.15 1.800 -0.31071874 0.051719529
#> 8 p4 easy 0.15 2.280 -0.83510245 -0.459367027Without preprocessing every cell’s drift is underestimated, and the
two largest errors belong to the two participants who contaminated most.
Screening moves every cell, mostly toward the truth, but the heaviest
contaminator’s easy condition stays well below it. Four participants
with 100 trials per cell cannot separate the rate’s effect from sampling
error, and this vignette does not try to. Whether a participant’s error
tracks their own contamination rate needs more participants than this
vignette has. r_contaminated() generates them with the
truth attached, and the ground-truth article on the package website (https://www.gfrischkorn.org/rtprep/articles/ground-truth.html)
runs that check on data matched to a task of your own.
Comparing rules
Because the rules share one return shape, comparing them is one call.
screen_compare() applies every rule in a list and reports
how much each removed, how often each pair agrees, and how much their
excluded sets overlap.
cmp <- screen_compare(
rt_example$rt,
list(
cutoff = rule_cutoff(0.18, 3),
sd = rule_sd(2.5),
mad = rule_mad(2.5),
recursive = rule_recursive("modified"),
mixture = rule_mixture("lognormal")
),
.by = list(rt_example$id, rt_example$condition)
)
cmp
#> <rtprep comparison> 800 trials, 5 rules
#>
#> cutoff dropped 0.0%
#> sd dropped 3.2%
#> mad dropped 11.4%
#> recursive dropped 4.1%
#> mixture dropped 4.0%
#>
#> least agreement: cutoff vs mad, 88.6% of decisions (Jaccard 0.00)
cmp$agreement
#> rule_x rule_y agree jaccard n_only_x n_only_y n_both
#> 1 cutoff sd 0.96750 0.0000000 0 26 0
#> 2 cutoff mad 0.88625 0.0000000 0 91 0
#> 3 cutoff recursive 0.95875 0.0000000 0 33 0
#> 4 cutoff mixture 0.96000 0.0000000 0 32 0
#> 5 sd mad 0.91875 0.2857143 0 65 26
#> 6 sd recursive 0.99125 0.7878788 0 7 26
#> 7 sd mixture 0.99000 0.7575758 1 7 25
#> 8 mad recursive 0.92750 0.3626374 58 0 33
#> 9 mad mixture 0.92625 0.3516484 59 0 32
#> 10 recursive mixture 0.99875 0.9696970 1 0 32Agreement and the Jaccard overlap answer different questions, and diverge exactly where it matters: two rules that each drop 2% of trials and never the same ones agree on 96% of decisions and overlap not at all.
Checking the exclusions
One more check applies to real data, where the truth is not
available. If the fast trials a rule removed were guesses, their
accuracy should be at chance. check_guessing() tests that
with a Beta-Binomial Bayes factor. The ±2.5 SD criterion removed no fast
trial at all in this data set, so there is nothing for it to test there;
an absolute cutoff at 350 ms does remove some:
rt_example |>
mutate(rt_screen(rt, rule_cutoff(0.35, 3)), .by = c(id, condition)) |>
reframe(check_guessing(.keep, rt, response), .by = id) |>
select(id, n_tested, prop_upper, bf_01, bf_evidence)
#> id n_tested prop_upper bf_01 bf_evidence
#> 1 p1 4 0.5 1.875 anecdotal_for_guessing
#> 2 p2 5 0.6 1.875 anecdotal_for_guessing
#> 3 p3 5 0.4 1.875 anecdotal_for_guessing
#> 4 p4 5 0.6 1.875 anecdotal_for_guessingFour or five tested trials per participant give anecdotal evidence at
best, and that is the honest reading: the check needs a rule that
removes fast trials, and enough of them, before it can say anything.
Where n_tested is zero the test is silent, which is itself
informative about the rule.
Reporting what you did
A preprocessing step is part of the analysis, and a reader cannot repeat it from “outliers were removed”. Four things pin it down: which rule and at which setting, the grouping the criterion was computed within, how much it removed, and whether error trials went through the screen with the correct ones. All four are in the objects the code already produced, so none of them has to be typed from memory.
The rule and its setting print themselves, and the per-cell counts
come from screen_fits():
rule <- rule_sd(2.5)
rule
#> <rtprep rule> sd(2.5, mean, sd)
#> Exclude trials more than 2.5 x sd from the mean, computed per group.
fits <- rt_example |>
reframe(screen_fits(rt, rule), .by = c(id, condition))
fits |>
summarise(
cells = n(),
trials = sum(n_trials),
dropped = sum(n_dropped),
prop = sum(n_dropped) / sum(n_trials),
lowest_cell = min(prop_dropped),
highest_cell = max(prop_dropped)
)
#> cells trials dropped prop lowest_cell highest_cell
#> 1 8 800 26 0.0325 0.02 0.05Report the range across cells as well as the total. A criterion that removes 3.2% overall can be removing much more from one participant than another, and that spread is the thing this package exists to make visible.
The reasons say what the rule actually caught, which is worth checking before describing it:
rt_example |>
mutate(rt_screen(rt, rule), .by = c(id, condition)) |>
count(.rule, .reason)
#> .rule .reason n
#> 1 sd(2.5, mean, sd) too_slow 26
#> 2 sd(2.5, mean, sd) <NA> 774Which gives a Methods sentence that can be written from the output
rather than around it. report_screening() writes it:
scr <- rt_screen(
rt_example$rt, rule,
.by = list(participant = rt_example$id, condition = rt_example$condition)
)
report_screening(scr)
#> Response times were screened with a criterion of 2.5 standard
#> deviations around the mean, computed separately within each combination
#> of participant and condition (8 cells). This removed 26 of the 800
#> trials it screened (3.2%), between 2% and 5% per cell. All the
#> exclusions were slow trials. The criterion did not read accuracy, so
#> error trials passed through the screen on their response times alone.
#> The criterion is described by Miller (1991).
#>
#> [72 words; 1 reference; toBibtex(x) for BibTeX]Nothing in that paragraph was typed from memory, which is the point of generating it. Edit the rule above and the sentence follows; write the sentence by hand and it goes stale the first time the rule changes, with no warning and nothing to catch it.
The numbers it quotes stay reachable, for a sentence that has to be written differently:
rep <- report_screening(scr)
c(excluded = rep$n_excluded, screened = rep$n_screened, missing = rep$n_missing)
#> excluded screened missing
#> 26 800 0
rep$cells
#> n fitted min max
#> 1 8 8 0.02 0.05Note which denominator the percentage uses: the trials the criterion
actually saw. rt_example has no missing response times, so
here it is every trial. Where there are some, they never reached the
criterion, and counting them among its exclusions would overstate what
the rule did, so they get a sentence of their own instead.
The clause about accuracy is the one most often left out and the one
that most often changes the answer. The criterion here never looked at
response, so error trials passed through it on their
response times alone; a rule that does read accuracy, such as
rule_ewma() or
rule_mixture(use_accuracy = TRUE), makes the screen and the
dependent variable share information, and
report_screening() says so without being asked.
Where to go next
?rules documents every rule with the reference it
implements and the columns it adds to the fits table, and
?rules_compose covers rule_all(),
rule_any() and rule_then(), which combine
them. No single conventional rule reaches both ends of the distribution,
so combining a spread criterion with an accuracy control chart is often
better than tuning either. ?rt_summary covers the robust,
trimmed and mixture aggregation routes, ?adjust_accuracy
the accuracy correction that goes with the mixture route, and
?r_contaminated the three contaminant processes and how to
match the generator to a task of your own. ?rule_oracle
removes exactly the labelled contaminants, which is the ceiling any real
rule is read against.
?rtprep-glossary defines the terms the rest of the
documentation assumes. rule_custom() turns a function of
your own into a rule in one call, and ?extending gives the
full contract: what the function receives, and what it has to
return.