Nobody's data arrives ready. Items come in columns when models need rows, three waves live in three files, a scale needs reverse-scoring, timestamps are text, and forty people skipped an item. This volume is the unglamorous engine of every paper — and the part that separates people who can run an analysis from people who can run it on their own data.
9
days, 45–55 min
12
copy-ready recipes
1
cleaning script you keep for years
Day 845 minutes · concepts, almost no code
What counts as one row
Tidy data has three properties: each variable is a column, each observation is a row, each type of observational unit is its own table. The hard word is observation, and its meaning changes with your design. In a between-subjects study, one row is one person. In a repeated-measures study, one row is one person at one occasion. In a reaction-time task, one row is one trial. Your dataset needs to be tidy for the model you intend to fit, not tidy in the abstract.
Wide: one row per person
Columns bdi_t1, bdi_t2, bdi_t3. How data arrives from Qualtrics, and what a correlation matrix or a classical repeated-measures table wants.
Long: one row per observation
Columns id, time, bdi. What ggplot2, afex, lme4 and brms all require. This is the format you should store in.
Item-level: one row per response
Columns id, item, response. Needed for item-level psychometrics and IRT; usually pivoted from wide item grids.
library(tidyverse)
d <- read_csv("data/study1_raw.csv")
# 1. What does one row represent? Look, do not assume.
head(d, 3)
# 2. Is any person present more than once?
d |> count(id) |> filter(n > 1)
# 3. Are variable names encoding data? (bdi_t1 hides a "time" variable)
names(d)
# 4. Is one column holding two things? ("cbt_male", "2024-03-01 14:22")
d |> count(condition)#> # A tibble: 0 × 2 ← no duplicated ids, good
#> [1] "id" "condition" "age" "bdi_t1" "bdi_t2" "bdi_t3" "gad_t1" "gad_t2" "gad_t3"
#> # A tibble: 3 × 2
#> condition n
#> 1 cbt 84
#> 2 control 82
#> 3 waitlist 82
Those column names are the tell: bdi_t1 to bdi_t3 means "time" is a variable that has been smeared across headers. Tomorrow you pull it back into a column, and in doing so make every longitudinal model in this course possible.
Do this now · 20 minutes
Take your own dataset and write four comment lines answering the four questions above. Then write one sentence: "For the model I intend to fit, one row should be ___." Keep that sentence at the top of your cleaning script — it is the specification everything else serves.
Day 950 minutes · the two functions that unlock longitudinal work
pivot_longer, pivot_wider
Reshaping has a reputation for being hard. It is not — it is one function with four arguments you will memorise this week: which columns to gather, what the names become, what the values become, and how to split a compound name.
long <- d |>
pivot_longer(
cols = c(bdi_t1, bdi_t2, bdi_t3),
names_to = "time",
values_to = "bdi",
names_prefix = "bdi_t" # strips "bdi_t", leaving 1, 2, 3
) |>
mutate(time = as.integer(time))
long |> select(id, condition, time, bdi) |> head(6)
nrow(d); nrow(long) # 248 people become 744 observations#> # A tibble: 6 × 4
#> id condition time bdi
#> <chr> <chr> <int> <dbl>
#> 1 p001 cbt 1 31
#> 2 p001 cbt 2 22
#> 3 p001 cbt 3 18
#> 4 p002 control 1 26
#> 5 p002 control 2 25
#> 6 p002 control 3 27
#> [1] 248
#> [1] 744# Columns are bdi_t1, bdi_t2, bdi_t3, gad_t1, gad_t2, gad_t3:
# the names encode BOTH the measure and the time point.
long2 <- d |>
pivot_longer(
cols = matches("^(bdi|gad)_t[0-9]$"),
names_to = c("measure", "time"),
names_pattern = "(bdi|gad)_t([0-9])",
values_to = "score"
) |>
mutate(time = as.integer(time))
# ...and then one column per measure, keeping one row per person-occasion
analysis <- long2 |>
pivot_wider(names_from = measure, values_from = score)
head(analysis, 4)#> # A tibble: 4 × 6
#> id condition age time bdi gad
#> <chr> <chr> <dbl> <int> <dbl> <dbl>
#> 1 p001 cbt 29 1 31 14
#> 2 p001 cbt 29 2 22 11
#> 3 p001 cbt 29 3 18 9
#> 4 p002 control 41 1 26 17
That last object is the canonical analysis dataset for a longitudinal psychology study: one row per person per occasion, one column per measure, person-level variables repeated down the rows. Everything in volume 5 assumes this shape.
When pivot_wider warns about list-columns
The warning "values are not uniquely identified" means the rows you are widening are not unique on the identifying columns — usually a duplicated ID or a forgotten grouping variable. Do not reach for values_fn = mean to silence it; find the duplicates with count(id, time) |> filter(n > 1) first. Silencing that warning is how people average two different people together.
Do this now · 20 minutes
Pivot your own item grid or repeated measure into long format, check the row count equals people × occasions (minus genuine missings), then pivot it back to wide and verify you recover the original with all.equal(). Round-tripping proves you understand both directions.
Day 1050 minutes · where silent data loss happens
Joins, and auditing them before you trust them
Three waves in three files, a demographics sheet, a separate therapist log — joining is routine. It is also the most dangerous operation in this volume, because a mismatched ID does not raise an error. It quietly produces NA, or quietly duplicates rows, and your n changes without anyone noticing.
demo <- read_csv("data/demographics.csv") # id, age, gender, education
outcome <- read_csv("data/wave1.csv") # id, bdi, gad
# left_join keeps every row of the left table
d <- outcome |> left_join(demo, by = "id")
# THE AUDIT — three lines, every single time
nrow(outcome); nrow(d) # must be equal unless you intended growth
outcome |> anti_join(demo, by = "id") # who is in outcomes but missing demographics
demo |> anti_join(outcome, by = "id") # who has demographics but no outcome#> [1] 248
#> [1] 248
#> # A tibble: 3 × 3 ← three participants have no demographic record
#> id bdi gad
#> 1 p117 24 12
#> 2 p203 19 8
#> 3 p244 31 16
#> # A tibble: 12 × 4 ← twelve consented but never completed wave 1
Those two anti_join results are not an inconvenience — they are your attrition and data-management section. Report them. A CONSORT flow diagram is built from exactly these numbers.
outcome |> left_join(demo, by = "id") # keep all outcomes (the usual choice)
outcome |> inner_join(demo, by = "id") # complete cases on both — states your exclusion
outcome |> full_join(demo, by = "id") # keep everyone from both sides
outcome |> anti_join(demo, by = "id") # the mismatches only — an audit, not a dataset
# Different ID column names, and a common ID disaster
outcome |> left_join(demo, by = join_by(id == participant_id))
demo |> count(id) |> filter(n > 1) # duplicated keys multiply your rows
# Whitespace and case are the usual culprits behind "missing" matches
demo <- demo |> mutate(id = str_squish(str_to_lower(id)))#> # A tibble: 1 × 2
#> id n
#> 1 p088 2 ← this ONE duplicate would add 248 phantom rows if joined naively
Do this now · 20 minutes
Join two of your own files. Run the three-line audit and write the resulting numbers into a comment block: rows before, rows after, unmatched on each side. If the two sides do not match, find out why before doing anything else — that reason belongs in your methods section.
Day 1145 minutes · a statistics day disguised as a cleaning day
Factors decide what your coefficients mean
A factor's first level is the reference group, and every coefficient in your regression is a comparison against it. Level order is therefore a modelling decision that happens during cleaning, which is why so many people get a results table they cannot interpret.
library(forcats)
d <- d |>
mutate(
condition = factor(condition),
condition = fct_relevel(condition, "control"), # control becomes the baseline
condition = fct_recode(condition,
"CBT" = "cbt",
"Waitlist" = "waitlist"), # publication-ready labels
severity = fct_relevel(severity, "minimal", "mild", "moderate", "severe")
)
levels(d$condition)
levels(d$severity)
# Ordering a factor BY a number makes every plot readable at once
d |> mutate(condition = fct_reorder(condition, bdi, .fun = median)) |> levels()#> [1] "control" "CBT" "Waitlist"
#> [1] "minimal" "mild" "moderate" "severe"m1 <- lm(bdi ~ condition, data = d)
coef(m1) # each coefficient is "this condition minus control"
d2 <- d |> mutate(condition = fct_relevel(condition, "CBT"))
coef(lm(bdi ~ condition, data = d2)) # same data, same model, different questions#> (Intercept) conditionCBT conditionWaitlist
#> 21.049 -1.732 -1.888
#> (Intercept) conditioncontrol conditionWaitlist
#> 19.317 1.732 -0.156
Identical data, identical fit, different table. Nothing is wrong in either — but only one of them answers the question your preregistration asked. Set the reference level on purpose, and say in the paper which group it is.
Do this now · 15 minutes
Convert every categorical variable in your own data to a factor with deliberate level order and publication-ready labels, in one mutate() block. Print levels() for each and write a comment naming your reference group and why it is the right comparison.
Day 1245 minutes · messy human answers
Strings: cleaning what people actually typed
Free-text fields arrive as "Female", "female ", "F", "femal". Six stringr functions handle almost all of it, and the discipline is always the same: clean, then count, then look at what is left.
library(stringr)
d <- d |>
mutate(
gender_raw = gender,
gender = gender |> str_squish() |> str_to_lower(), # trim and normalise first
gender = case_when(
str_detect(gender, "^f") ~ "female",
str_detect(gender, "^m") ~ "male",
str_detect(gender, "non|nb|enby") ~ "non-binary",
.default = NA_character_
)
)
# ALWAYS look at what did not classify — never assume it is empty
d |> filter(is.na(gender)) |> count(gender_raw)#> # A tibble: 3 × 2
#> gender_raw n
#> 1 "prefer not to" 4
#> 2 "woman" 2
#> 3 "" 1
"woman" failed the ^f test. That is why you inspect the unclassified group rather than trusting your patterns — and why the raw column stays in the dataset. Never overwrite what a participant actually wrote.
str_detect(x, "anxiet") # contains a pattern
str_replace_all(x, ",", ".") # decimal commas in numeric text
str_extract(x, "[0-9]+") # pull the number out of "approx 12 sessions"
str_split_i(x, ";", 1) # first item of a multi-select answer
str_length(x) # response length: a data-quality flag
str_starts(x, "p0") # ID prefixes
# Open-text quality check: near-empty or suspiciously long answers
d |> mutate(n_char = str_length(notes)) |>
filter(n_char < 3 | n_char > 600) |>
select(id, n_char)#> # A tibble: 5 × 2
#> id n_char
#> 1 p034 1
#> 2 p088 2
#> 3 p156 812
Do this now · 15 minutes
Clean one free-text column in your own data with the trim → lowercase → classify → inspect-the-remainder sequence. Keep the raw column. Report how many responses did not classify and what you did with them.
Day 1345 minutes · needed for every diary study
Dates and times, which are never just text
Timestamps carry three research variables for free: elapsed time since baseline, time of day, and compliance. Parse them once at import and all three become available.
library(lubridate)
esm <- esm |>
mutate(
ts = ymd_hms(timestamp, tz = "Asia/Kolkata"), # parse ONCE
day = as_date(ts),
hour = hour(ts),
weekday = wday(ts, label = TRUE, week_start = 1),
beep_num = row_number(), .by = id
)
# Days since each person's first observation: the time variable for growth models
esm <- esm |>
mutate(t0 = min(ts), .by = id) |>
mutate(days_in = as.numeric(difftime(ts, t0, units = "days")))
# Compliance per person — a result, not an afterthought
esm |> summarise(beeps = n(),
span_days = round(max(days_in), 1),
compliance = round(n() / (7 * 5), 2), # 7 days x 5 beeps
.by = id) |> head(4)#> # A tibble: 4 × 4
#> id beeps span_days compliance
#> 1 p001 31 6.9 0.89
#> 2 p002 27 7.0 0.77
#> 3 p003 35 6.8 1
#> 4 p004 14 6.9 0.4
Two traps that ruin diary datasets
Time zones. If you do not set tz, R assumes UTC and your evening beeps move to the next day, which destroys day-level nesting. Ambiguous formats. 03/04/2025 is March or April depending on the platform's locale; use dmy() or mdy() explicitly rather than a generic parser, and check the minimum and maximum date afterwards.
Do this now · 15 minutes
Parse a date or timestamp column in your data with an explicit format and time zone. Derive elapsed time since each person's first record. Print the minimum and maximum date and check both are possible — impossible dates are the fastest way to catch a parsing error.
Day 1450 minutes · the day psychology data becomes psychology
Scoring scales: reverse-keying, means, and the rules you set in advance
Turning twenty items into one number involves at least four decisions: which items are reverse-keyed, mean or sum, how many missing items a person may have, and whether the scale is unidimensional at all. The first three are cleaning; the fourth is volume 6. Make all of them explicitly, in code, before you look at any result.
The rule "at most one missing item, then prorate the mean" is defensible and common. "Average whatever is there" is not, because someone who answered one item gets the same weight as someone who answered five. Whatever you choose, write it in the preregistration and report the number it excluded.
keys <- list(
agree = c("A1_r", "A2", "A3", "A4", "A5"),
consc = c("C1", "C2", "C3", "C4_r", "C5_r"),
extra = c("E1_r", "E2_r", "E3", "E4", "E5"),
neuro = c("N1", "N2", "N3", "N4", "N5"),
open = c("O1", "O2_r", "O3", "O4", "O5_r")
)
score_scale <- function(data, items, max_miss = 1) {
sub <- data[, items, drop = FALSE]
miss <- rowSums(is.na(sub))
out <- rowMeans(sub, na.rm = TRUE)
out[miss > max_miss] <- NA_real_
out
}
for (nm in names(keys)) scored[[nm]] <- score_scale(scored, keys[[nm]])
scored |> select(all_of(names(keys))) |> psych::describe() |> round(2)#> vars n mean sd median trimmed mad min max range skew kurtosis
#> agree 1 2782 4.65 0.89 4.80 4.72 0.89 1 6 5 -0.59 0.06
#> consc 2 2775 4.27 0.95 4.40 4.31 1.04 1 6 5 -0.37 -0.24
#> extra 3 2784 4.15 1.06 4.20 4.18 1.19 1 6 5 -0.28 -0.42
#> neuro 4 2778 3.18 1.20 3.20 3.17 1.19 1 6 5 0.03 -0.61
#> open 5 2780 4.59 0.81 4.60 4.62 0.89 1 6 5 -0.41 -0.08
You just wrote your first function — on day 14, because the alternative was writing the same five lines five times. That instinct ("this is repeating, make it a function") is what day 94 formalises.
Do this now · 20 minutes
Score every scale in your own data with an explicit reverse-key list, an explicit missingness rule, and a printed count of how many people the rule excluded. Save it as scripts/04-scoring.R. This file is the one your co-authors will ask for.
Day 1555 minutes · a methods-section day
Missing data: look first, decide second
Listwise deletion is a choice with consequences, not a neutral default — it costs power and, unless data are missing completely at random, biases your estimates. The decision depends on the pattern, so the first move is always to look.
library(naniar)
miss_summary <- d |> miss_var_summary() # per variable: n and % missing
miss_summary |> head(5)
n_case_complete(d); n_case_miss(d) # complete vs incomplete cases
vis_miss(d) # the pattern, visually
gg_miss_upset(d) # which COMBINATIONS co-occur
# Is missingness related to observed variables? (evidence against MCAR)
d |> mutate(miss_bdi = is.na(bdi_t3)) |>
summarise(m_baseline = mean(bdi_t1, na.rm = TRUE),
n = n(), .by = miss_bdi)#> # A tibble: 5 × 3
#> variable n_miss pct_miss
#> 1 notes 31 12.5
#> 2 bdi_t3 22 8.87
#> 3 gad_t3 19 7.66
#> 4 bdi_t2 7 2.82
#> 5 age 0 0
#> [1] 198
#> [1] 50
#> # A tibble: 2 × 3
#> miss_bdi m_baseline n
#> 1 FALSE 19.1 226
#> 2 TRUE 25.7 22 ← dropouts were more depressed at baseline
That last table is the whole lesson. The people who dropped out were more depressed to begin with, so complete-case analysis would report the treatment working better than it did. This is missing at random (MAR) — predictable from observed data — and it is exactly the case where multiple imputation or a likelihood-based model earns its keep.
library(mice)
# 1. Include your outcome and any variable predicting missingness
imp <- mice(d |> select(bdi_t1, bdi_t2, bdi_t3, age, condition),
m = 20, method = "pmm", seed = 2026, printFlag = FALSE)
# 2. Fit the model in EVERY imputed dataset
fits <- with(imp, lm(bdi_t3 ~ condition + bdi_t1))
# 3. Pool with Rubin's rules — never analyse one imputed set
pooled <- pool(fits)
summary(pooled, conf.int = TRUE) |> round(3)#> term estimate std.error statistic df p.value 2.5 % 97.5 %
#> 1 (Intercept) 7.412 1.503 4.932 191.442 0.000 4.447 10.377
#> 2 conditionCBT -3.918 0.988 -3.965 176.203 0.000 -5.868 -1.968
#> 3 conditionWaitlist -0.402 1.002 -0.401 182.771 0.689 -2.380 1.576
#> 4 bdi_t1 0.604 0.061 9.902 188.114 0.000 0.484 0.724
What to write in the paper
Percentage missing per variable, the pattern, whether missingness related to observed variables, your approach and why, and — if imputed — the number of imputations, the method, the variables in the imputation model, and that estimates were pooled with Rubin's rules. Four sentences. Reviewers in clinical psychology check for them specifically.
Do this now · 25 minutes
Produce a missingness table and plot for your own data, then test whether missingness on your main outcome relates to any baseline variable. Write the four-sentence missing-data paragraph now, while the numbers are in front of you.
Day 1650 minutes · the volume's deliverable
One cleaning script, and raw data you never touch again
The structure below is what your project looks like for the next four months. Numbered scripts, raw data treated as read-only, derived data written to its own folder, and a log of what every step cost you in cases.
r-course/
├── r-course.Rproj
├── data-raw/ # never edited, never overwritten, ideally read-only
│ ├── wave1.csv
│ └── demographics.csv
├── data/ # derived, reproducible, safe to delete
│ └── analysis.rds
├── scripts/
│ ├── 01-import.R
│ ├── 02-clean.R
│ ├── 03-score.R
│ └── 04-analysis.R
├── output/
│ ├── figures/
│ └── tables/
└── README.md # what this project is, and how to rerun it# ---- 02-clean.R --------------------------------------------------------
library(tidyverse); library(here)
raw <- read_csv(here("data-raw", "wave1.csv"), na = c("", "NA", "-99"))
log <- tibble(step = "imported", n = nrow(raw))
d <- raw |>
mutate(id = str_squish(str_to_lower(id)))
log <- add_row(log, step = "ids normalised", n = nrow(d))
d <- d |> filter(!duplicated(id))
log <- add_row(log, step = "duplicate ids removed", n = nrow(d))
d <- d |> filter(consent == 1, attention_check == 1)
log <- add_row(log, step = "consent + attention check", n = nrow(d))
d <- d |> filter(!is.na(bdi_t1))
log <- add_row(log, step = "baseline outcome present", n = nrow(d))
log |> mutate(lost = lag(n) - n) # your participant flow, computed not remembered
saveRDS(d, here("data", "analysis.rds"))
write_csv(log, here("output", "tables", "case_log.csv"))#> # A tibble: 5 × 3
#> step n lost
#> <chr> <int> <int>
#> 1 imported 263 NA
#> 2 ids normalised 263 0
#> 3 duplicate ids removed 261 2
#> 4 consent + attention check 251 10
#> 5 baseline outcome present 248 3
That table is your flow diagram, generated by the code rather than reconstructed from memory eight months later when a reviewer asks. It takes six extra lines and has saved more papers than any statistical technique in this course.
install.packages("renv")
renv::init() # creates a project-local library and renv.lock
renv::snapshot() # records exact versions of everything you used
sessionInfo() # paste the output into your manuscript's appendix#> * Lockfile written to '~/r-course/renv.lock'
#> R version 4.4.1 (2024-06-14), dplyr 1.1.4, lme4 1.1-35.5, lavaan 0.6-18
Do this now · 25 minutes
Restructure your project into the folder layout above, move the raw file into data-raw/ and mark it read-only at the operating-system level, and rewrite your cleaning as 02-clean.R with a case log. Run renv::init(). Restart R and run scripts 01 to 03 in order from cold.
Reference
Twelve recipes you will reuse constantly
{{ r.t }}
{{ r.c }}
Checkpoint
Six questions before volume 3
{{ quizCounter }}
{{ quizScore }}
{{ quizQ }}
{{ quizFb }}
Next: volume 3
Graphics is eight days on the grammar of graphics and on the specific figures behavioural journals expect — distributions rather than bar charts, within-subject lines, multi-panel layouts, and export at exact millimetre widths and 300 dpi.