--- title: "Get started with rtprep" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Get started with rtprep} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", # dplyr is a Suggests: the vignette reads without it and runs with it eval = requireNamespace("dplyr", quietly = TRUE) ) ``` 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. ```{r setup, message = FALSE} library(rtprep) library(dplyr) ``` ## 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. ```{r data} head(rt_example) rt_example |> count(id, contam_rate, contaminant) |> filter(contaminant) ``` The 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. ```{r screen} screened <- rt_example |> mutate(rt_screen(rt, rule_sd(2.5)), .by = c(id, condition)) head(screened) ``` 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: ```{r rules} rule_cutoff(0.18, 3) rule_mad(2.5) rule_recursive("modified") rule_mixture("lognormal") ``` 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. ```{r keep} clean <- rt_example |> filter(rt_keep(rt, rule_sd(2.5), .by = list(id, condition))) nrow(clean) ``` The per-cell diagnostics, which `mutate()` would strip from the attribute, come from `screen_fits()`: ```{r 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) ``` ## Did 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: ```{r score} screened |> summarise( sensitivity = mean(!.keep[contaminant]), specificity = mean(.keep[!contaminant]), .by = id ) ``` The ±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: ```{r score-process} screened |> filter(contaminant) |> summarise(found = mean(!.keep), n = n(), .by = process) ``` The 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: ```{r estimate} 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) ``` `rt_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: ```{r error} 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) ``` Without 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 () 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. ```{r compare} 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 cmp$agreement ``` Agreement 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: ```{r guessing} 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) ``` Four 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()`: ```{r reporting-fits} rule <- rule_sd(2.5) rule 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) ) ``` Report the range across cells as well as the total. A criterion that removes `r sprintf("%.1f%%", 100 * sum(fits$n_dropped) / sum(fits$n_trials))` 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: ```{r reporting-reasons} rt_example |> mutate(rt_screen(rt, rule), .by = c(id, condition)) |> count(.rule, .reason) ``` Which gives a Methods sentence that can be written from the output rather than around it. `report_screening()` writes it: ```{r reporting-paragraph} scr <- rt_screen( rt_example$rt, rule, .by = list(participant = rt_example$id, condition = rt_example$condition) ) report_screening(scr) ``` 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: ```{r reporting-numbers} rep <- report_screening(scr) c(excluded = rep$n_excluded, screened = rep$n_screened, missing = rep$n_missing) rep$cells ``` Note 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.