How to read this chapter. Sections 1.1–1.5 set up why the book exists and what responsible practice requires. Sections 1.6–1.12 build the computational foundation everything else stands on. If you already write R daily, skim 1.6–1.11, but do work through §1.6.6 (Reproducibility), §1.6.7 (Data access), and §1.12 (Simulation), the book’s later chapters assume the conventions established there.

1 Learning objectives

After completing this chapter you will be able to:

Describe the defining characteristics of complex biomedical and health data, and identify which of them a given study exhibits.
Distinguish prediction, explanation, and causal inference, and state which one a given research question requires.
Apply FAIR and research-ethics principles to a concrete dataset, including computing a quantitative re-identification risk.
Configure a reproducible R environment: project structure, dependency pinning, seed discipline, and session recording.
Manipulate R’s core data structures, atomic vectors, lists, data frames, matrices, and arrays, and predict the result of type coercion.
Import, reshape, and export tabular data using tidy-data operations, and explain why long ↔︎ wide reshaping is invertible.
Simulate synthetic data matching the marginal characteristics of an observed dataset, and assess the fidelity of the simulation.
Evaluate a quantitative claim for provenance, deciding whether a published table, benchmark, or trend line supports the conclusion drawn from it.

Estimated time: 4–6 hours including exercises. Prerequisites: none. Some prior exposure to any programming language helps but is not assumed. Review the DSPA prerequisites and self-assessment pretest if you want to calibrate.

2 What this book is for

DSPA3 is the third edition of the textbook built around the HS650: Data Science and Predictive Analytics course at the University of Michigan (DSPA1 · DSPA2). It aims to give you a working understanding of the challenges, an appreciation of the opportunities, and a methodological foundation for designing, collecting, managing, processing, interrogating, analyzing, and interpreting complex health and biomedical data.

Vision, enable active learning by joining driving motivational challenges to mathematical foundations, computational statistics, and modern scientific inference.
Values, effective, reliable, reproducible, and transformative data-driven discovery in support of open science.
Strategic priorities, learners develop scientific intuition, computational skill, and data-wrangling ability sufficient to attack real biomedical and health data problems. Instructors receive documented, runnable code implementing everything from atomic data filters to end-to-end predictive analytics pipelines.

2.1 What is expected of you

You do not need to satisfy every prerequisite, master every mathematical foundation, or arrive as an experienced programmer.
You do not need to read linearly. Several suggested pathways through the topics map are sufficient for most readers.
The material will challenge you. When it does: engage with peers, search, read the referenced sources, and post in the discussion forum. The difficulty is the point, a skill that costs nothing to acquire is worth what it cost.
Some chapters will be too easy for you. Skim them, but attempt the practice problems anyway; they are calibrated to expose gaps that reading does not.
All code is available as an electronic supplement, and every example in this book is designed to run as printed.
These materials can be improved. Errors, typos, better code, interesting datasets, send them to DSPA.info @ umich.edu. Corrections are incorporated into subsequent revisions with attribution.

3 Driving problems

Every method in this book exists because someone needed to answer a question they could not answer by looking at the data directly. This section introduces the problems that recur throughout DSPA3.

The rightmost column is a forward map. It tells you where each dataset returns, so you can trace a single scientific problem across multiple analytic treatments, usually the fastest way to understand what a method actually buys you.

3.1 Neurodegeneration and dementia

Alzheimer’s disease. Can we relate observed clinical phenotypes to expected behavior? Can we prognosticate cognitive decline 3–12 months forward from imaging and clinical assessment? Can we recover meaningful subject clusters from harmonized multi-source data?

Source	Modalities & scale	Recurs in
ADNI · NACC	Clinical (demographics, cognitive assessment); imaging (sMRI, fMRI, DTI, PiB/FDG PET); genetics (Illumina SNP); chemical biomarkers (labs, proteomics). Cohort sizes typically $200 \le N \le 1200$	Dimensionality reduction, feature selection

A study of structural neuroimaging in AD illustrates the difficulty concretely: 808 ADNI subjects (200 AD, 383 MCI, 225 normal controls), sMRI parcellated with BrainParser, 80 neuroimaging biomarkers extracted by shape-analysis workflow, and 80 SNPs associated with those biomarkers via a pipeline implementation of PLINK. The correlations observed between genetic and imaging phenotypes suggest that AD/MCI/NC differences are best examined through joint models of morphometric, imaging, and genotypic data, not through any single modality.

Parkinson’s disease. Can we predict clinical diagnosis from available data, with and without the UPDRS assessment that currently defines the diagnosis? Can derived imaging and genetic biomarkers model progression well enough to support clinical decisions?

Source	Modalities & scale	Recurs in
PPMI	Demographics; clinical (physical, verbal learning, neurological, UPSIT olfactory, vitals, MDS-UPDRS, ADL, MoCA, GDS-15); sMRI; genetics (Illumina ImmunoChip, 196,524 variants; NeuroX, ~240,000 exonic variants; APOE e2/e3/e4). Three cohorts: de-novo PD $N=263$, SWEDD $N=40$, controls $N=127$. Visits at screening, baseline, 12, 24, 48 months	Simulation §1.12, dimensionality reduction, feature selection, longitudinal analysis

Amyotrophic lateral sclerosis. Which variables jointly predict ALS progression (ALSFRS, muscle function)? Can adverse events be predicted from phenotype plus 0–3 month clinical change?

Source	Modalities & scale	Recurs in
PRO-ACT	>100 clinical variables; 8,635 patients in the training set, 2,424 with valid gold-standard ALSFRS slopes. Demographics, ALSFRS, adverse events, `onset_delta`, `onset_site`, riluzole use	Feature selection, the primary ALS case study

3.2 Substance use and withdrawal

Is the Risk for Alcohol Withdrawal Syndrome (RAWS) screen valid and reliable for predicting withdrawal among adult medical inpatients? What AUDIT-C cutoff optimizes prediction? Should items be added or removed?

Source	Modalities & scale	Recurs in
MAWS / UMHS EHR · WHO Global Health Observatory	AUDIT-C scores (see Bush et al., 1998), including any current use (item 1), total AUDIT-C > 8, and prior withdrawal history. ~1,000 positive cases/year among ~10,000 adult inpatients	Classification

3.3 Genomics and high-throughput computing

Genetic forensics, the 2013–2016 Ebola outbreak. The HHMI BioInteractive disease-detective activity tracks Ebola virus spread through Sierra Leone using the fact that most of a genome is shared within a species, more so among related individuals. DNA profiling exploits the residual differences, particularly in noncoding regions, where variation is less exposed to selection. Short tandem repeats (STRs) are 2–5 base motifs repeated a variable number of times at defined loci, flanked by non-variable regions. An allele is named either by repeat count or by total length; the STR below is “6” (six GATA repeats) or equivalently 70 bp including flanks.

Next-generation sequencing. Whole-genome and exome sequencing carry the clues needed to identify genes behind Mendelian disorders, and population-genetic extensions push toward complex traits. The Graphical Pipeline for Computational Genomics (GPCG) implements alignment, quality control, SNP analysis, CNV identification, annotation, and visualization as composable workflow steps.

Neuroimaging-genetics. A computational infrastructure for high-throughput neuroimaging-genetics supports aggregation, harmonization, processing, and interpretation across imaging, genetics, clinical, and cognitive sources. Workflows are portable XML objects that carry execution instructions from a client to distributed pipeline servers, an early and durable example of computational provenance as a first-class artifact.

3.4 Interactive exploration

The SOCR Brain Visualization App preloads sMRI, ROI labels, and fiber tract models for a normal brain, and accepts drag-and-drop of your own stereotactic data.

3.5 Open data repositories

Datasets used throughout DSPA3 are drawn from open-access archives that follow FAIR principles (§1.5.1):

4 What makes health data hard

4.1 Beyond the “V” vocabulary

The familiar three-V framing, volume, velocity, variety, originates with Doug Laney’s 2001 META Group research note, 3D Data Management. Veracity was added later and popularized by IBM, producing the common “4 V’s.” The framing is widely repeated and worth knowing, but it was written for enterprise data management, and it does not tell a biomedical analyst what tool to reach for.

A constructive characterization is more actionable, because each dimension names a methodological gap and therefore a class of tool:

Dimension	What it means for your analysis	Where DSPA3 addresses it
Size	Data exceeds convenient memory or single-machine compute	Scaling & streaming
Complexity	Heterogeneous types in one study: images, sequences, text, tabular	Throughout
Incongruency	Sources disagree on units, coding, timing, or inclusion	Data handling
Multi-source	Joint modeling of elements never designed to be joined	Dimensionality reduction
Multi-scale	Macro → meso → micro observations of the same process	Longitudinal analysis
Time	Repeated, irregular, censored, or drifting measurement	Longitudinal analysis
Incomplete	Missingness that is itself informative	Missing data
Provenance	Where each value came from and what was done to it	§1.6.6, and every chapter

The last row is our addition. Provenance is not usually listed among the V’s, and its absence is arguably the most common cause of irreproducible health-data analysis: not that the statistics were wrong, but that nobody can reconstruct which version of which file was processed by which code.

4.2 Data science and predictive analytics

Data science is transdisciplinary, spanning theoretical, computational, experimental, and biosocial work, and deals with large volumes of complex, incongruent, dynamic, multi-source data. Its goal is algorithms, methods, tools, and services that ingest such data and produce semi-automated decision support: mining patterns, predicting outcomes, proposing groupings, computing signatures, and extracting actionable evidence. In practice this means data wrangling, harmonization and aggregation, exploratory and confirmatory analysis, predictive modeling, validation, and tuning.

Predictive analytics uses mathematical formulation, statistical computing, and software to represent, interrogate, and interpret data with the aim of forecasting: predicting patterns, or prognosticating process behavior inside or outside the observed range. Here process means the natural phenomenon under investigation, examined through proxy data.

The core task is to identify relationships, associations, arrangements, or motifs, across space, time, and features, that reduce the effective complexity of the data. From those, predictive analytics produces unknown outcomes, likelihood or parameter estimates, classification labels, and other aggregate or individualized forecasts. The assumptions of a technique determine its usability, bound its expected accuracy, and constrain what human action may responsibly follow from its output. The type of the target outcome, binary, polytomous, probability, scalar, vector, tensor, determines whether the appropriate strategy yields prediction, forecasting, labeling, likelihoods, grouping, or motif discovery.

4.3 Prediction is not explanation is not causation

This distinction is the most frequently blurred idea in applied data science, and blurring it produces confident, useless, and occasionally harmful results. Three different questions demand three different modeling strategies:

Goal	Question	Optimizes	Typical failure
Prediction	What will happen?	Out-of-sample accuracy	Works until the data-generating process shifts; unusable for intervention
Explanation	What is associated with the outcome, and how strongly?	Interpretable, well-specified structure	Confuses association with mechanism
Causal inference	What happens if I intervene?	Identifiability under stated assumptions	Assumes away confounding that is present

A model can be excellent at one and useless at the others. A model that predicts hospital readmission accurately from whether a discharge-planning note was filed tells you nothing about whether filing more notes would reduce readmissions. See Breiman’s Statistical Modeling: The Two Cultures (2001) and Shmueli’s To Explain or To Predict? (2010) for the canonical treatments.

Most of DSPA3 is about prediction, with explanation appearing in the regression and feature-selection chapters. When a method is being used for a purpose it was not designed for, we say so.

5 Responsible and reproducible practice

Technical skill is necessary and not sufficient. Analysts must also observe professional norms governing appropriate use, reproducibility, unbiased reporting, and both expected and unanticipated interpretation of their results. Five obligations organize the rest of this section: FAIR sharing, research ethics, honest accounting of benefits and harms, regulatory competence around sensitive data, and documented provenance.

5.1 FAIR resource sharing

Findable, Accessible, Interoperable, Reusable (Wilkinson et al., 2016) are the guiding principles for developing, deploying, and stewarding data, methods, tools, and services. FAIR is a property of the resource, not of the researcher’s intentions: a dataset with no persistent identifier is not findable no matter how much you wish it were.

5.2 Research ethics

Ethical data science demands responsible conduct and integrity across the whole investigation, planning, design, implementation, execution, and assessment. Being unaware of a norm is not an exemption from it.

5.3 Benefits and detriments of findings

Evidence-driven discovery generates questions as well as answers, and some are unexpected, unwelcome, or harmful. You are responsible for validating your results and for explicitly weighing the potential benefits against the probable detriments of both positive and negative findings.

5.4 Sensitive data: regulation and practice

Decisions about security, privacy, and confidentiality are not trivial, and the number of people, devices, algorithms, and services within arm’s length of raw data argues for layered protection. Under the U.S. HIPAA Privacy Rule, health data may be de-identified two ways:

Safe Harbor, remove 18 enumerated identifier categories (names, geographic subdivisions smaller than a state, all date elements finer than year, ages over 89, and so on).
Expert Determination, a qualified person applies statistical principles and documents that re-identification risk is very small.

Note what Safe Harbor is not: it is not a guarantee of anonymity, and it is not a statistical statement. It is a checklist. Expert determination requires actually computing something, which is the subject of §1.5.6.

Protection must persist across the data life cycle and may require preemptive, ongoing, and post-hoc analysis to find and patch vulnerabilities. There are genuine tradeoffs between the value of data and the risk of unrestrained inference from it; neither extreme is practical or sustainable.

5.5 Provenance and longevity

The digitalization of human experience and the promise of AI have produced enormous investment, excitement, and anxiety, along with the expectation that abundant information translates directly into insight, prediction, and cost-effective decisions. Recording the data, algorithmic, scientific, computational, and human factors behind a forecast is not bureaucracy, it is the part of the work that makes the forecast knowledge rather than output. §1.6.6 makes this operational.

5.6 Quantifying re-identification risk

Ethics sections are usually prose. This one is code, because re-identification risk is a computable quantity and treating it as a matter of good intentions is how it gets ignored.

Quasi-identifiers are variables that do not identify anyone on their own, but do so in combination. The classic triple is ZIP code, date of birth, and sex. $k$-anonymity (Sweeney, 2002) holds for a dataset when every record shares its quasi-identifier combination with at least $k-1$ others. When $k = 1$ for some record, that record is unique in the data and can be matched against any external source containing the same variables.

Let us measure it.

set.seed(2026)

n <- 5000
patients <- data.frame(
  zip5       = sample(48103:48109, n, replace = TRUE),
  birth_year = sample(1935:2005,   n, replace = TRUE),
  sex        = sample(c("F", "M"), n, replace = TRUE),
  # A sensitive attribute that de-identification is meant to protect:
  diagnosis  = sample(c("HC", "PD", "SWEDD"), n, replace = TRUE,
                      prob = c(0.5, 0.35, 0.15))
)

# A reusable k-anonymity report over any set of quasi-identifiers.
k_anonymity <- function(data, quasi_ids) {
  key        <- do.call(paste, c(data[quasi_ids], sep = "\r"))
  class_size <- table(key)
  sizes      <- as.integer(class_size[key])
  list(
    quasi_ids     = quasi_ids,
    k             = min(sizes),
    n_unique      = sum(sizes == 1L),
    pct_unique    = round(100 * mean(sizes == 1L), 2),
    pct_below_k5  = round(100 * mean(sizes < 5L),  2),
    n_classes     = length(class_size)
  )
}

str(k_anonymity(patients, c("zip5", "birth_year", "sex")))

#> List of 6
#>  $ quasi_ids   : chr [1:3] "zip5" "birth_year" "sex"
#>  $ k           : int 1
#>  $ n_unique    : int 38
#>  $ pct_unique  : num 0.76
#>  $ pct_below_k5: num 25.3
#>  $ n_classes   : int 985

Now watch what generalization buys you. Truncating ZIP to three digits and binning birth year into decades are both Safe-Harbor-flavored moves:

patients$zip3        <- patients$zip5 %/% 100
patients$birth_decade <- (patients$birth_year %/% 10) * 10

designs <- list(
  "zip5 + birth year + sex"   = c("zip5", "birth_year",   "sex"),
  "zip3 + birth year + sex"   = c("zip3", "birth_year",   "sex"),
  "zip3 + birth decade + sex" = c("zip3", "birth_decade", "sex"),
  "zip3 + birth decade"       = c("zip3", "birth_decade")
)

report <- do.call(rbind, lapply(names(designs), function(nm) {
  r <- k_anonymity(patients, designs[[nm]])
  data.frame(design = nm, k = r$k, unique_records = r$n_unique,
             pct_unique = r$pct_unique, pct_below_k5 = r$pct_below_k5,
             equivalence_classes = r$n_classes)
}))
report

Two lessons the numbers make unavoidable:

Generalization works, and it is not free. Each step upward reduces unique records and raises $k$, while destroying analytic resolution. You cannot study ZIP-level geographic variation in data generalized to three-digit ZIPs. Privacy and utility trade against each other along a curve, and someone chooses the operating point, make sure it is chosen deliberately and documented.
$k$-anonymity is a floor, not a guarantee. If every record in an equivalence class shares the same diagnosis, an attacker who locates the class learns the diagnosis without identifying the individual. That is homogeneity attack, and defending against it requires stronger notions, $\ell$-diversity, $t$-closeness, or differential privacy.

5.6.1 A cautionary tale, read carefully

The most-cited example in health privacy is the 1997 re-identification of Massachusetts Governor William Weld’s hospital record in a Group Insurance Commission dataset stripped of direct identifiers, linked against the Cambridge voter roll. The episode directly shaped the HIPAA de-identification provisions.

It is worth knowing that this account is contested. Barth-Jones (2012) re-examined the Cambridge demographics and argued that Weld was likely identifiable principally because he was a public figure whose hospitalization had been widely reported, not because the demographic linkage was reliable, and therefore that the attack does not generalize as cleanly as its popular retelling suggests.

We include both because the disagreement is instructive in two directions. Substantively, “how identifiable is de-identified data?” remains an open, quantitative question rather than a settled anecdote. Methodologically: read what you cite. The Barth-Jones paper is frequently cited as evidence for the Weld re-identification, which inverts its argument. Citations propagate faster than anyone re-checks them, and one of the more valuable habits this book can give you is the reflex to follow a reference to its source before repeating it, including references in this book.

5.7 When resources are misused

Beyond the guiding principles, some concrete failures:

Re-identification attempts and policy circumvention, see §1.5.7. Attempting to reverse de-identification, or to route around a data-use agreement, is a breach regardless of technical difficulty.
Uncalibrated reporting. Findings reported without reference to status-quo alternatives mislead. Ignoring placebo effects, methodological assumptions, conflicts of interest, or unmodeled randomization can materially change the apparent efficacy of a data-driven study.
Unintended misuse of legitimate access. In 2014 an Uber executive used the company’s internal “God View” tool to track a journalist’s location without consent, a technically impressive use of available data and a clear ethical breach. It contributed to a subsequent FTC settlement.
Gaming for personal gain. In 2015 a Morgan Stanley financial adviser downloaded account data on roughly 350,000 wealth-management clients, about 10% of that division’s clients, and portions surfaced online.
Statistical misuse, ranging from naive error to deliberate distortion.
Innocent error at scale. Selection bias, methodological drift, processing faults, unstable estimation, and misread metadata all produce wrong answers without anyone intending harm. At scale, innocent error and malicious manipulation are indistinguishable from the outside, which is the argument for provenance (§1.5.5).
Collecting the wrong data. Irrelevant volume adds cost, skews understanding, and distracts from real findings. More data is not a strategy.

6 The DSPA3 computing environment

6.1 Why R?

R is not the only reasonable choice, and this book’s use of it is a decision with tradeoffs rather than a verdict. The honest case:

What R is unusually good at

Statistical depth. Methods appear in R first. When a biostatistics paper publishes a method, the reference implementation is usually an R package.
Data-to-publication path. R Markdown and Quarto make a document that regenerates its own figures and numbers from source. This book is that artifact.
Domain coverage. Bioconductor alone provides thousands of curated packages for genomics, proteomics, and imaging.
Free and open. R is distributed under GPL-2/GPL-3, so anyone can inspect, audit, and extend the source. Note the actual copyleft condition: obligations attach when you distribute modified binaries or source, not when you modify privately. You may adapt R inside your institution indefinitely without publishing anything.
Interoperability. Native bridges to C/C++ (Rcpp), Python (reticulate), Java, databases, and web services.

What R is not good at

General-purpose software engineering. Large applications, services, and production infrastructure are better served elsewhere.
Idiomatic consistency. Three object systems (S3, S4, R5/R6), two major dialects (base and tidyverse), and decades of accumulated interfaces. R’s learning curve is steeper than the syntax suggests.
Raw single-threaded speed on tight numeric loops, absent vectorization or compiled extensions, though this matters far less than beginners assume (§1.6.2).

Let us verify a few claims rather than assert them:

R.version.string

#> [1] "R version 4.3.3 (2024-02-29 ucrt)"

paste("Platform:", R.version$platform)

#> [1] "Platform: x86_64-w64-mingw32"

At time of writing, the current release is R 4.6.1 (“Happy Hop”, 2026-06-24). If the line above reports something much older, see §1.6.4.

# Compute the CRAN package count rather than quoting a number that will be
# wrong by the time you read this.
cran <- available.packages(repos = "https://cloud.r-project.org")
cat("CRAN packages available today:", nrow(cran), "\n")

#> CRAN packages available today: 24437

The general rule, applied here and everywhere in DSPA3: if a number can be computed at render time, compute it. A hard-coded package count is wrong within months; a computed one never is. Watch for this pattern in the wild, when a textbook or paper quotes a figure without a date or a method, that is information about how much to trust it.

R is also polyglot. An .Rmd or .qmd document can execute chunks in many languages, which is what makes it viable as a protocol format rather than just an R notebook:

engines <- names(knitr::knit_engines$get())
cat("knitr supports", length(engines), "engines, including:\n")

#> knitr supports 52 engines, including:

print(head(sort(engines), 24))

#>  [1] "asis"      "asy"       "awk"       "bash"      "block"     "block2"   
#>  [7] "bslib"     "c"         "cat"       "cc"        "coffee"    "comment"  
#> [13] "css"       "ditaa"     "dot"       "embed"     "eviews"    "exec"     
#> [19] "fortran"   "fortran95" "gawk"      "go"        "groovy"    "haskell"

6.2 Reading a benchmark critically

You will encounter tables ranking languages by speed. Most are unusable, because a benchmark without a stated workload, machine, version, and implementation is not a measurement. Rather than reprint one, run one:

set.seed(7)
x <- rnorm(2e6)

# Three implementations of exactly the same computation.
loop_sum <- function(v) {
  total <- 0
  for (i in seq_along(v)) total <- total + v[i]
  total
}
vapply_sum    <- function(v) sum(vapply(v, identity, numeric(1)))
vectorized_sum <- function(v) sum(v)

timings <- rbind(
  `explicit for-loop` = system.time(loop_sum(x))[["elapsed"]],
  `vapply`            = system.time(vapply_sum(x))[["elapsed"]],
  `vectorized sum()`  = system.time(vectorized_sum(x))[["elapsed"]]
)
colnames(timings) <- "elapsed_seconds"
timings

#>                   elapsed_seconds
#> explicit for-loop            0.03
#> vapply                       0.69
#> vectorized sum()             0.00

cat("\nSlowdown vs. vectorized:\n")

#> 
#> Slowdown vs. vectorized:

round(timings[, 1] / timings["vectorized sum()", 1], 1)

#> explicit for-loop            vapply  vectorized sum() 
#>               Inf               Inf               NaN

The ratio between the slowest and fastest row is typically two to three orders of magnitude, larger than the gap most cross-language benchmarks report between R and a compiled language. The lesson is not that R is fast. It is that for most analytic code, the dominant performance factor is how you write it, not which language you write it in, and that a language-level speed ranking cannot tell you anything useful about your own code.

Three questions to ask of any benchmark you are shown:

What was measured? A tight scalar loop and a matrix decomposition rank languages in different orders.
Who ran it, when, on what? Version, hardware, BLAS implementation, and compiler flags all move results by large factors.
Is it reproducible? If you cannot re-run it, it is an assertion.

6.3 Measuring adoption: a construct-validity lesson

How popular is R? The question sounds empirical and is mostly a measurement problem. Two commonly used proxies, and what is wrong with each:

Proxy 1, literature mentions. The dataset below counts Google Scholar hits for R, SAS, and SPSS by publication year.

pubs <- read.csv(
  "https://umich.instructure.com/files/2361245/download?download_frd=1",
  header = TRUE
)
head(pubs)

library(ggplot2)
library(tidyr)

pubs_long <- pivot_longer(pubs, cols = c("R", "SAS", "SPSS"),
                          names_to = "software", values_to = "hits")

ggplot(pubs_long, aes(x = Year, y = hits, colour = software)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 1.6) +
  labs(
    title    = "Google Scholar mentions of statistical software, 1995-2015",
    subtitle = "Historical snapshot - see caveats below; not a measure of usage",
    x = "Year", y = "Scholar hits", colour = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

At least four things are wrong with treating this as evidence of adoption:

Mention ≠ use. A paper criticizing SPSS counts as an SPSS hit.
Not reproducible. Scholar hit counts change between queries and are not versioned; nobody can re-derive this series.
Ambiguity. “R” matches an enormous amount of non-R text.
It stops in 2015. Presenting an eleven-year-old series in the present tense is its own error.

We keep the figure because the caveats are the lesson. A plot can be correctly drawn from correctly-loaded data and still not support the claim placed underneath it.

Proxy 2, search interest. Google Trends is better in one respect (it is reproducible) and has a subtler flaw. Compare these two queries:

Searching the string "R" returns volume for a letter of the alphabet.
Searching the topic /m/0212jm returns volume for the entity R (programming language), as Google’s knowledge graph resolves it.

These measure different things, and the difference is not small.

# Not run: gtrendsR is rate-limited (HTTP 429) under repeated queries.
# Note the topic IDs - NOT the bare keyword strings.
library(gtrendsR)
trends <- gtrends(
  keyword = c("/m/0212jm",   # R (programming language)
              "/m/018fh1",   # SAS (software)
              "/m/02l0yf8",  # SPSS
              "/m/05z1_"),   # Python (programming language)
  gprop = "web", time = "all"
)[["interest_over_time"]]

# Pre-fetched snapshot (retrieved July 2023) so the chapter knits offline.
gt <- read.csv(
  "https://umich.instructure.com/files/31071103/download?download_frd=1",
  header = TRUE
)
str(gt, give.attr = FALSE)

#> 'data.frame':    235 obs. of  5 variables:
#>  $ Month : chr  "2004-01" "2004-02" "2004-03" "2004-04" ...
#>  $ R     : int  61 59 61 56 53 56 58 55 56 59 ...
#>  $ SAS   : int  11 12 11 11 12 12 12 12 12 12 ...
#>  $ Python: int  13 13 12 12 12 12 11 12 13 12 ...
#>  $ SPSS  : int  3 4 4 4 4 3 3 3 3 3 ...

gt$Month <- as.Date(paste0(gt$Month, "-01"))

gt_long <- pivot_longer(gt, cols = -Month,
                        names_to = "software", values_to = "rsv")

ggplot(gt_long, aes(x = Month, y = rsv, colour = software)) +
  geom_line(linewidth = 0.9) +
  labs(
    title    = "Google Trends relative search volume, 2004-2023",
    subtitle = "RSV is scaled to 100 at the series maximum - it is not a count",
    x = NULL, y = "Relative search volume", colour = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Construct validity is the question of whether a measurement measures the thing you named. It is not a statistical property you can test after the fact, it is a claim about the relationship between your operational definition and your concept, and it has to be argued. You will meet this again every time a model is trained on a proxy outcome: billing codes standing in for diagnosis, readmission standing in for treatment failure, engagement standing in for value. The failure mode is not a bad model; it is a good model of the wrong thing.

6.4 Installing R

R is free and runs on Windows, macOS, and Linux. From r-project.org:

Click CRAN and pick a mirror, e.g. the U-M mirror or the global cloud mirror, which routes automatically.
Choose your operating system, then base.
Download and install the current release (4.6.x or newer).

Verify the installation before continuing:

R.version.string
stopifnot(getRversion() >= "4.1.0")   # DSPA3 requires the native pipe |>

6.4.1 Keeping `R` current

Upgrade R roughly every 6–12 months and packages more often. The procedure differs by platform.

Platform	Upgrading `R`
Windows	`install.packages("installr"); installr::updateR()`, Windows only
macOS	Download the new `.pkg` from CRAN and install over the old version, or `brew install --cask r`
Linux	Distribution package manager, or the Posit r-builds project for side-by-side versions

Updating packages, on any platform:

update.packages(ask = FALSE, checkBuilt = TRUE)

checkBuilt = TRUE matters after a minor-version R upgrade: packages built against the previous version need recompilation, and skipping this produces confusing load failures weeks later.

6.5 Choosing an IDE

You can drive R from a terminal, but almost nobody does. Two mature options:

RStudio Desktop, the long-standing default. Free, open source, extremely well documented. Four panes:

Bottom-left, console. Where R actually evaluates.
Top-left, editor. Scripts and notebooks. Ctrl/Cmd + Enter sends the current line or selection to the console.
Top-right, environment / history. What is in memory; what you typed.
Bottom-right, files / plots / packages / help / viewer.

Positron, Posit’s newer data-science IDE, built on Code OSS. Polyglot by design (R and Python as equal citizens), integrated data explorer, native notebook editor, and package management pane. For a book that mixes R with reticulate-driven Python, Positron is the more natural fit.

Either works for every example in DSPA3. If you already know RStudio, stay; if you are starting fresh and expect to write Python too, try Positron.

6.6 Installing packages

Base R ships with a modest core; everything else is a package. Two mechanisms:

# Classic
install.packages(c("ggplot2", "dplyr", "tidyr"))

# Modern: pak resolves dependencies better, installs in parallel, and gives
# far clearer errors when a system library is missing.
install.packages("pak")
pak::pak(c("ggplot2", "dplyr", "tidyr"))

Installing puts a package on disk. Attaching makes its functions visible:

library(ggplot2)      # attach - errors immediately if not installed
ggplot2::geom_point   # qualified access - no attach needed

Prefer pkg::function() in scripts that others will run: it documents where each function came from and avoids masking surprises when two packages export the same name.

Pitfall. library(base) does nothing, base, stats, utils, graphics, grDevices, datasets, and methods are attached at startup. Calling library() on them implies otherwise.

6.6.1 The DSPA3 dependency budget

The full book touches roughly 150 packages. You do not need them all. This core set covers Chapters 1–4:

pak::pak(c(
  "ggplot2", "dplyr", "tidyr", "readr",      # data + graphics
  "knitr", "rmarkdown",                       # authoring
  "plotly", "DT",                             # interactivity
  "renv", "here",                             # reproducibility
  "Hmisc"                                     # descriptive summaries
))

Later chapters declare their own additional requirements at the top. Install per-chapter rather than all at once, it keeps the failure surface small.

6.7 Reproducibility

The book’s stated values include reproducible. That word has to cash out as mechanics, or it is decoration. Five practices, in increasing order of effort and payoff.

6.7.1 1. Use projects, not `setwd()`

setwd("/Users/you/Desktop/analysis") breaks the moment anyone else, or future you on another machine, runs the code. Instead, create an RStudio / Positron project (a directory with an .Rproj file) and reference files relative to its root:

library(here)
dat <- read.csv(here("data", "raw", "ppmi.csv"))

here() resolves against the project root regardless of where the session started.

6.7.2 2. Set seeds, always

Any code touching a pseudo-random number generator must be seeded, or its output changes on every run, including outputs your prose describes.

set.seed(42); round(rnorm(3), 4)

#> [1]  1.3710 -0.5647  0.3631

set.seed(42); round(rnorm(3), 4)   # identical

#> [1]  1.3710 -0.5647  0.3631

rnorm(3)                           # different: no reseed

#> [1]  0.6328626  0.4042683 -0.1061245

Convention for DSPA3: every chunk containing a random draw calls set.seed() first. Not once at the top of the file, per chunk. Chunks get re-run, cached, and reordered during editing; a single seed at the top does not survive that.

6.7.3 3. Never write to the working directory

Knitting a document should not leave files in your source tree.

out <- file.path(tempdir(), "example.csv")
write.csv(head(iris), out, row.names = FALSE)
file.exists(out)

#> [1] TRUE

head(read.csv(out), 3)

6.7.4 4. Pin your dependencies with `renv`

“It worked last year” is not reproducibility. renv records the exact version of every package in a project lockfile:

renv::init()      # create a project-local library + renv.lock
renv::snapshot()  # record current package versions into renv.lock
renv::restore()   # on another machine: reinstall exactly those versions

Commit renv.lock to version control. It is the difference between “install these packages” and “install these packages as they were.”

6.7.5 5. Record the session

Every DSPA3 chapter ends with its session information. Include it in anything you share:

sessionInfo()

Checkpoint 1. You have a colleague’s analysis script that runs on their laptop and fails on yours. Which of the five practices above, had they used it, would most likely have prevented the failure? (Answer: it depends on the error, but renv and here() between them eliminate the two most common causes, version drift and path assumptions.)

6.8 Accessing data reproducibly

DSPA3 loads most datasets over the network. That creates three problems: the chapter fails offline, it re-downloads on every knit, and it breaks permanently when a URL rots. One small helper solves all three, and we use it book-wide.

# ---- Cache-aware data loader ----------------------------------------------
dspa_cache_dir <- function() {
  d <- Sys.getenv("DSPA_CACHE", unset = file.path(tempdir(), "dspa-cache"))
  if (!dir.exists(d)) dir.create(d, recursive = TRUE)
  d
}

#' Download once, read many times.
#'
#' @param url     Remote location of the file.
#' @param name    Stable local filename used as the cache key.
#' @param reader  Function applied to the cached path (default `read.csv`).
#' @param refresh Force re-download even if cached.
dspa_read <- function(url, name, reader = utils::read.csv, ..., refresh = FALSE) {
  stopifnot(is.character(url),  length(url)  == 1L,
            is.character(name), length(name) == 1L)
  path <- file.path(dspa_cache_dir(), name)

  if (refresh || !file.exists(path)) {
    ok <- tryCatch({
      utils::download.file(url, destfile = path, mode = "wb", quiet = TRUE)
      file.exists(path) && file.size(path) > 0L
    }, error = function(e) FALSE, warning = function(w) FALSE)

    if (!ok) {
      if (file.exists(path)) unlink(path)
      stop("Could not retrieve '", name, "' from:\n  ", url,
           "\nIf you are offline, place the file at:\n  ", path,
           call. = FALSE)
    }
  }
  reader(path, ...)
}

#' Fingerprint a cached file so you can prove which version you analyzed.
dspa_fingerprint <- function(name) {
  path <- file.path(dspa_cache_dir(), name)
  if (!file.exists(path)) return(NA_character_)
  unname(tools::md5sum(path))
}

Using it:

ppmi <- dspa_read(
  url  = "https://umich.instructure.com/files/330397/download?download_frd=1",
  name = "05_PPMI_top_UPDRS_Integrated_LongFormat1.csv"
)
dim(ppmi)

#> [1] 1764   31

cat("MD5:", dspa_fingerprint("05_PPMI_top_UPDRS_Integrated_LongFormat1.csv"), "\n")

#> MD5: 03e27492af90cc68f0a8497483d183df

The fingerprint is the point. A results table is not reproducible unless you can demonstrate which bytes produced it. Record the checksum in your methods section alongside the URL and retrieval date. (tools::md5sum() needs no extra package; for anything security-sensitive prefer openssl::sha256().)

To persist the cache between sessions, point it somewhere durable:

# In ~/.Renviron  (then restart R)
DSPA_CACHE=~/dspa-data

6.9 Getting help

?lm                      # help for a known function
help("lm")               # equivalent
??"linear model"         # fuzzy search across installed packages
example(lm)              # run the documented examples
vignette(package = "dplyr")   # long-form tutorials
browseVignettes()             # everything installed
args(lm)                 # just the signature

Beyond the console:

CRAN Task Views, curated package lists by domain (Bayesian, Survival, MachineLearning, MedicalImaging…). The fastest route from “what should I use for X” to a shortlist.
R for Data Science, 2nd ed., free online; the standard modern introduction.
Advanced R, 2nd ed., for understanding why R behaves as it does.
UCLA OARC Statistical Methods and Data Analytics , worked examples across many methods.
Posit Community and Stack Overflow [r], post a reprex (minimal reproducible example) and you will usually get an answer within hours.

6.10 Working with an AI coding assistant

Most readers of this edition will use an LLM assistant while working through it. Pretending otherwise would be unhelpful, so here is the honest guidance.

What they are good at

Explaining unfamiliar syntax, error messages, and package APIs.
Producing a first draft of boilerplate, plotting scaffolds, reshaping, string manipulation.
Translating between idioms (“rewrite this base R using dplyr”).
Suggesting what to check when something is broken.

Where they fail, specifically in statistical work

Confidently wrong statistics. Assistants readily produce code that runs, returns a plausible number, and answers the wrong question, the wrong test for a design, an inappropriate correction, cross-validation that leaks. Nothing errors. You get a result.
Hallucinated APIs. Functions and arguments that do not exist, or existed in a prior version.
Stale idioms. Trained on years of code including deprecated patterns.
No access to your data. Assumptions about distributions, missingness, and coding that the assistant cannot verify and often does not flag.

A workable discipline

Never paste sensitive data into a general-purpose assistant. PHI, PII, and data under a use agreement do not leave your environment. This is a compliance requirement, not a preference.
Verify against documentation. If it suggests foo(bar = TRUE), run ?foo before believing it.
Ask for the reasoning, not just the code. “Why this test?” surfaces errors that reading the code does not.
Test on data with a known answer. Simulate a dataset where you know the truth (§1.12) and confirm the code recovers it. This single habit catches most silent statistical errors, from assistants and from yourself.
You own the output. “The AI wrote it” is not a methods section.

The skill this book teaches, knowing what a method assumes, what it estimates, and how it fails, is exactly the skill required to use an assistant safely. It gets more valuable as the tools improve, not less.

6.11 Quarto: authoring reproducible documents

Quarto is the successor to R Markdown: multi-language, multi-format, and maintained as the primary authoring path by Posit. Existing .Rmd files render under Quarto largely unmodified; new documents use .qmd.

Install from quarto.org/docs/get-started (RStudio and Positron bundle a copy). Create one via File → New File → Quarto Document, or the command palette (Ctrl/Cmd + Shift + P → “Quarto”).

Three features that matter for this book:

Inline code. Prose can compute. An inline R expression, written as `r nrow(iris)`, is evaluated at render time, so this sentence can state that the iris dataset has 150 rows without anyone ever typing that number. This is the mechanism behind the computed statistics throughout this chapter.

Multi-format output. One source renders to HTML, PDF, Word, or slides.

Multi-language chunks. {r}, {python}, {julia}, {ojs}, and many more, in a single document with shared state where the runtimes support it. Chapter 16’s Python integration depends on this.

Quarto also supports Observable JS for reactive, client-side interactive figures, useful when you want the reader to explore rather than watch.

7 R essentials

This section covers what you need to read and write the code in the rest of the book. It is deliberately short. Comprehensive function tables live in Appendix A; reach for them when you need a name, not when you are learning.

7.1 Assignment and the two pipes

x <- 2          # preferred assignment operator
y = 3           # also legal, but reserve `=` for function arguments
x + y

#> [1] 5

R gained a native pipe in version 4.1. DSPA3 uses it throughout:

c(4, 9, 16, 25) |> sqrt() |> sum()

#> [1] 14

# Equivalent, nested - harder to read as the chain grows
sum(sqrt(c(4, 9, 16, 25)))

#> [1] 14

You will encounter %>% (from magrittr, re-exported by dplyr) in older code and in many packages’ documentation. The two are near-equivalent for simple chains. |> requires no package and is slightly faster; %>% supports a . placeholder in more positions. DSPA3 uses |>.

R also gained shorthand lambda syntax in 4.1:

sapply(1:5, \(z) z^2)              # modern

#> [1]  1  4  9 16 25

sapply(1:5, function(z) z^2)       # classic - identical

#> [1]  1  4  9 16 25

Pitfall, T and F. TRUE and FALSE are reserved words. T and F are ordinary variables that happen to be initialized to them, and can be reassigned:

T                # TRUE

#> [1] TRUE

T <- FALSE       # perfectly legal R
T                # FALSE - and every `header = T` in your script is now wrong

#> [1] FALSE

rm(T)            # restore sanity

Always write TRUE/FALSE. Six extra keystrokes eliminate an entire category of silent bug.

7.2 Atomic vectors and coercion

The vector is R’s fundamental unit. There is no separate scalar type, a “single number” is a length-1 vector.

nums  <- c(1, 2, 3, 5, 6, 7, 10, 1, 4)
nums

#> [1]  1  2  3  5  6  7 10  1  4

length(nums)

#> [1] 9

typeof(nums)

#> [1] "double"

An atomic vector holds one type. Mixing types triggers coercion, along a fixed hierarchy: logical → integer → double → character.

c(TRUE, 1L, 2.5)              # -> double

#> [1] 1.0 1.0 2.5

c(TRUE, 1L, 2.5, "a")         # -> character; the numbers became strings

#> [1] "TRUE" "1"    "2.5"  "a"

typeof(c(TRUE, 1L, 2.5, "a"))

#> [1] "character"

This is the single most common source of confusing R bugs. One stray character value in a column silently converts the entire column to text, and every downstream numeric operation then either fails or, worse, succeeds on the wrong thing. You will see a full worked instance in §1.12.3.

Indexing is 1-based, unlike C, Python, or Java:

nums[1]                 # first element (not nums[0])

#> [1] 1

nums[-1]                # everything except the first

#> [1]  2  3  5  6  7 10  1  4

nums[2:4]               # a contiguous slice

#> [1] 2 3 5

nums[c(1, 4, 2)]        # arbitrary positions, in the order given

#> [1] 1 5 2

nums[nums > 4]          # logical subsetting - the R idiom

#> [1]  5  6  7 10

nums[nums > 3 & nums < 8]

#> [1] 5 6 7 4

Logical subsetting deserves emphasis: nums > 4 produces a logical vector of the same length, and indexing by it keeps the TRUE positions. Nearly all “filter the data” operations in R are this pattern.

7.3 Generating structured data

seq(1, 10, by = 2.5)          # by increment

#> [1] 1.0 3.5 6.0 8.5

seq(1, 10, length.out = 4)    # by count

#> [1]  1  4  7 10

seq_along(c("a", "b", "c"))   # safe 1..n - never fails on length 0

#> [1] 1 2 3

rep(c(1, 2, 3), times = 4)    # whole vector repeated

#>  [1] 1 2 3 1 2 3 1 2 3 1 2 3

rep(c(1, 2, 3), each  = 4)    # each element repeated

#>  [1] 1 1 1 1 2 2 2 2 3 3 3 3

1:5                           # the `:` shorthand

#> [1] 1 2 3 4 5

Pitfall. 1:length(x) gives c(1, 0) when x is empty, silently iterating twice over nothing. seq_along(x) gives integer(0), which is what you meant. Use seq_along() and seq_len().

7.4 Lists

A list can hold anything, including other lists, R’s general container.

l <- list(a = c(1, 2), b = "hi", c = -3 + 3i, d = list(nested = TRUE))
str(l)

#> List of 4
#>  $ a: num [1:2] 1 2
#>  $ b: chr "hi"
#>  $ c: cplx -3+3i
#>  $ d:List of 1
#>   ..$ nested: logi TRUE

Extraction uses three distinct operators, and confusing them is a rite of passage:

l["a"]      # single bracket -> a LIST containing element a

#> $a
#> [1] 1 2

l[["a"]]    # double bracket -> the element ITSELF

#> [1] 1 2

l$a         # $ -> same as [[ ]], with partial matching

#> [1] 1 2

l$a[[2]]    # index into the extracted vector

#> [1] 2

class(l["a"])     # list

#> [1] "list"

class(l[["a"]])   # numeric

#> [1] "numeric"

Mnemonic: a list is a train of boxcars. [ returns a shorter train; [[ returns the cargo.

7.5 Matrices, arrays, and tensors

A matrix is a vector with a dim attribute, which is why the coercion rules of §1.7.2 apply to it in full.

m <- matrix(1:12, nrow = 3, ncol = 4)
m

#>      [,1] [,2] [,3] [,4]
#> [1,]    1    4    7   10
#> [2,]    2    5    8   11
#> [3,]    3    6    9   12

dim(m)

#> [1] 3 4

m[2, 3]        # row 2, column 3

#> [1] 8

m[2, ]         # all of row 2

#> [1]  2  5  8 11

m[, 3]         # all of column 3

#> [1] 7 8 9

t(m)           # transpose

#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    4    5    6
#> [3,]    7    8    9
#> [4,]   10   11   12

Arrays generalize to any number of dimensions. In machine learning these are usually called tensors, indexed by their order (number of dimensions):

ar <- array(1:24, dim = c(3, 4, 2))    # two 3x4 slices
dim(ar)

#> [1] 3 4 2

ar[2, 3, 1]    # row 2, col 3, slice 1

#> [1] 8

ar[2,  , 1]    # row 2 of slice 1

#> [1]  2  5  8 11

Dimension names make code self-documenting:

dimnames(m) <- list(c("R1", "R2", "R3"), c("C1", "C2", "C3", "C4"))
m

#>    C1 C2 C3 C4
#> R1  1  4  7 10
#> R2  2  5  8 11
#> R3  3  6  9 12

m["R2", "C3"]

#> [1] 8

7.6 Data frames and tibbles

A data frame is a list of equal-length vectors presented as a table. It is the workhorse structure for statistics: columns may differ in type, rows are observations.

df <- data.frame(
  id     = 1:4,
  group  = c("a", "B", "C", "d"),
  score  = c(10.2, 11.7, 9.4, 13.1)
)
df

str(df)

#> 'data.frame':    4 obs. of  3 variables:
#>  $ id   : int  1 2 3 4
#>  $ group: chr  "a" "B" "C" "d"
#>  $ score: num  10.2 11.7 9.4 13.1

Recycling fills shorter columns, but only when lengths divide evenly:

data.frame(v = 1:4, n = c(10, 11))     # length 2 recycles into length 4

data.frame(v = 1:4, n = c(10, 11, 12)) # 3 does not divide 4 -> error

#> Error in `data.frame()`:
#> ! arguments imply differing number of rows: 4, 3

The second case is far more common in practice than the first, and the error message is your friend, silent recycling of mismatched data would be worse.

Tibbles (tibble/dplyr) are data frames with better printing, no partial matching, and no surprising type conversions. DSPA3 uses plain data frames where base R suffices and tibbles inside tidyverse pipelines; they are interchangeable for our purposes.

7.6.1 The `stringsAsFactors` story

Before R 4.0.0 (April 2020), data.frame() and read.table() silently converted character columns to factor. Previous recommendations to pass as.is = TRUE and stringsAsFactors = FALSE to prevent it have changed more recently to keep the character columns stay as characters.

class(data.frame(x = c("a", "b"))$x)   # "character" in R >= 4.0.0

#> [1] "character"

The old advice is now a no-op, and repeating it teaches a model of R that has been wrong for six years. Factors remain the right representation for genuinely categorical variables, you just now create them deliberately:

grp <- factor(c("PD", "HC", "SWEDD", "PD"),
              levels = c("HC", "PD", "SWEDD"))
grp

#> [1] PD    HC    SWEDD PD   
#> Levels: HC PD SWEDD

levels(grp)

#> [1] "HC"    "PD"    "SWEDD"

table(grp)

#> grp
#>    HC    PD SWEDD 
#>     1     2     1

as.integer(grp)     # the underlying integer codes

#> [1] 2 1 3 2

Setting levels explicitly is not optional in modeling work: the first level is the reference category in lm(), glm(), and most model formulas. Left to R, levels are alphabetical, so your reference group is whichever label sorts first, rarely what you intended.

Pitfall. as.numeric(factor) returns the codes, not the labels. To recover numbers stored as a factor: as.numeric(as.character(f)).

7.7 Missing values

R distinguishes three “nothing” values, and they are not interchangeable:

Value	Means	Length
`NA`	A value exists but is unknown	1
`NULL`	There is no value at all	0
`NaN`	Undefined numeric result (`0/0`)	1

v <- c(1, 3, NA, 9, 8)
is.na(v)

#> [1] FALSE FALSE  TRUE FALSE FALSE

sum(v)                  # NA propagates - this is a feature

#> [1] NA

sum(v, na.rm = TRUE)    # explicit removal

#> [1] 21

mean(v, na.rm = TRUE)

#> [1] 5.25

NA propagating through arithmetic is deliberate: it forces you to decide what missingness means rather than silently ignoring it. Most summary functions accept na.rm.

dfm <- data.frame(a = 1:5, b = c(1, 3, NA, 9, 8))
na.omit(dfm)                          # drop incomplete rows

complete.cases(dfm)                   # which rows are complete

#> [1]  TRUE  TRUE FALSE  TRUE  TRUE

colSums(is.na(dfm))                   # missingness per column

#> a b 
#> 0 1

Warning. na.omit() is complete-case analysis. It is unbiased only when data are missing completely at random, and it discards information even then. The missing-data chapter treats MCAR/MAR/MNAR and imputation properly; until then, at minimum always report how many rows you dropped.

7.8 Inspecting objects

str(iris)               # structure - the first thing to run on anything

#> 'data.frame':    150 obs. of  5 variables:
#>  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
#>  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
#>  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
#>  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
#>  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

dim(iris); nrow(iris); ncol(iris)

#> [1] 150   5

#> [1] 150

#> [1] 5

names(iris)

#> [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"

class(iris); typeof(iris)

#> [1] "data.frame"

#> [1] "list"

head(iris, 3)

summary(iris)

#>   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
#>  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
#>  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
#>  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
#>  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
#>  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
#>  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
#>        Species  
#>  setosa    :50  
#>  versicolor:50  
#>  virginica :50  
#>                 
#>                 
#>

str() deserves special mention: it is the single most useful diagnostic function in R, and running it before anything else prevents a large fraction of type-related confusion.

7.9 Control flow and vectorization

x <- 10
if (x > 10) "big" else "not big"

#> [1] "not big"

for (i in 1:3) cat("iteration", i, "\n")

#> iteration 1 
#> iteration 2 
#> iteration 3

i <- 0
while (i < 3) { i <- i + 1 }
i

#> [1] 3

Prefer vectorized alternatives. They are shorter, faster (§1.6.2), and less error-prone:

scores <- c(55, 91, 78, 43, 88)

# Loop version
grade_loop <- character(length(scores))
for (i in seq_along(scores)) {
  grade_loop[i] <- if (scores[i] >= 60) "pass" else "fail"
}

# Vectorized - preferred
grade_vec <- ifelse(scores >= 60, "pass", "fail")

identical(grade_loop, grade_vec)

#> [1] TRUE

grade_vec

#> [1] "fail" "pass" "pass" "fail" "pass"

7.10 Writing functions

standardize <- function(x, na.rm = TRUE) {
  if (!is.numeric(x)) stop("`x` must be numeric.", call. = FALSE)
  (x - mean(x, na.rm = na.rm)) / sd(x, na.rm = na.rm)
}

round(standardize(c(2, 4, 4, 4, 5, 5, 7, 9)), 3)

#> [1] -1.403 -0.468 -0.468 -0.468  0.000  0.000  0.935  1.871

standardize(c("a", "b"))

#> Error:
#> ! `x` must be numeric.

Three habits worth adopting immediately:

Validate inputs and fail with a message that says what went wrong.
Return one thing. R returns the last evaluated expression; explicit return() is for early exits.
Avoid depending on the global environment. Everything a function needs should arrive as an argument. Functions that reach outside themselves work until they are moved.

7.11 Iteration over structures

The *apply family applies a function across elements or margins.

m2 <- matrix(1:6, nrow = 2)
apply(m2, 1, sum)      # 1 = over rows

#> [1]  9 12

apply(m2, 2, sum)      # 2 = over columns

#> [1]  3  7 11

lapply(list(a = 1:5, b = c(9, 90, 999)), median)   # returns a list

#> $a
#> [1] 3
#> 
#> $b
#> [1] 90

sapply(list(a = 1:5, b = c(9, 90, 999)), median)   # simplifies if it can

#>  a  b 
#>  3 90

vapply(list(a = 1:5, b = c(9, 90, 999)), median, numeric(1))  # type-safe

#>  a  b 
#>  3 90

Prefer vapply() in code others will run: it declares the expected return type and fails loudly if reality disagrees, where sapply() will quietly hand back a list instead of a vector and break something three steps later.

v   <- c(1, 2, 4, 2, 2, 5, 6, 4, 7, 8, 8)
fac <- factor(rep(c("low", "mid", "high"), length.out = 11),
              levels = c("low", "mid", "high"))
tapply(v, fac, sum)     # split-apply-combine over a factor

#>  low  mid high 
#>   17   16   16

7.12 Frequencies and cross-tabulation

table(v)

#> v
#> 1 2 4 5 6 7 8 
#> 1 3 2 1 1 1 2

prop.table(table(v))    # relative frequencies - note the nested table()

#> v
#>          1          2          4          5          6          7          8 
#> 0.09090909 0.27272727 0.18181818 0.09090909 0.09090909 0.09090909 0.18181818

Pitfall. prop.table(v), passing the raw vector rather than table(v), divides each observation by the sum of all observations. It runs without complaint and returns something that is not a frequency distribution. The argument to prop.table() must already be a table.

Two-way tabulation and independence testing:

adm <- as.data.frame(UCBAdmissions)   # 1973 Berkeley admissions
head(adm)

xtabs(Freq ~ Gender + Admit, data = adm)

#>         Admit
#> Gender   Admitted Rejected
#>   Male       1198     1493
#>   Female      557     1278

summary(xtabs(Freq ~ ., data = adm))

#> Call: xtabs(formula = Freq ~ ., data = adm)
#> Number of cases in table: 4526 
#> Number of factors: 3 
#> Test for independence of all factors:
#>  Chisq = 2000.3, df = 16, p-value = 0

This dataset is the textbook illustration of Simpson’s paradox: aggregated over departments, admission rates appear to favor men; within nearly every individual department, they do not. See Bickel, Hammel & O’Connell (1975). Keep it in mind, the same reversal appears whenever a confounder is marginalized away, which is most of the time in observational health data.

7.13 Combining and reshaping data frames

d1 <- data.frame(id = c(1, 1, 7, 6, 8), b = c(1, 1, NA, 9, 8))
d2 <- data.frame(id = c(1, 1, 7, 6, 8), c = 1:5)
merge(d1, d2, by = "id")

rbind(head(iris, 2), tail(iris, 2))        # stack rows

cbind(head(iris[1:2], 3), extra = 1:3)     # append columns

unique(d1)                                  # drop duplicate rows

sort(v); rev(sort(v))                       # ordering

#>  [1] 1 2 2 2 4 4 5 6 7 8 8

#>  [1] 8 8 7 6 5 4 4 2 2 2 1

which.max(v); which.min(v); which(v == 4)   # positions

#> [1] 10

#> [1] 1

#> [1] 3 8

cut(v, breaks = c(0, 3, 6, 10))             # bin a numeric into a factor

#>  [1] (0,3]  (0,3]  (3,6]  (0,3]  (0,3]  (3,6]  (3,6]  (3,6]  (6,10] (6,10]
#> [11] (6,10]
#> Levels: (0,3] (3,6] (6,10]

aggregate(Sepal.Length ~ Species, data = iris, FUN = mean)

Note on cbind(). For data frames, cbind() dispatches to cbind.data.frame() and preserves column types. For matrices and vectors it produces a matrix and coerces everything to a common type. The same function name, two very different behaviors, §1.12.3 shows what happens when you get it wrong.

7.14 Sampling

set.seed(11)
sample(1:10, size = 5)                        # without replacement

#> [1] 10  2  8  1  7

sample(1:10, size = 5, replace = TRUE)        # with replacement

#> [1] 5 6 5 6 7

sample(c("PD", "HC"), 8, replace = TRUE, prob = c(0.7, 0.3))

#> [1] "PD" "PD" "PD" "PD" "PD" "PD" "PD" "PD"

set.seed(11)
train_idx <- sample(seq_len(nrow(iris)), size = 0.7 * nrow(iris))
train <- iris[train_idx, ]
test  <- iris[-train_idx, ]
c(train = nrow(train), test = nrow(test))

#> train  test 
#>   105    45

That last idiom, sample indices, then use -idx for the complement, is the basic train/test split, and it reappears in every supervised-learning chapter.

8 Data import and export

8.1 Reading tabular data

hdp <- dspa_read(
  url  = "https://umich.instructure.com/files/1628650/download?download_frd=1",
  name = "01_hdp.csv"
)
dim(hdp)

#> [1] 8525   27

str(hdp, list.len = 8)

#> 'data.frame':    8525 obs. of  27 variables:
#>  $ tumorsize   : num  68 64.7 51.6 86.4 53.4 ...
#>  $ co2         : num  1.53 1.68 1.53 1.45 1.57 ...
#>  $ pain        : int  4 2 6 3 3 4 3 3 4 5 ...
#>  $ wound       : int  4 3 3 3 4 5 4 3 4 4 ...
#>  $ mobility    : int  2 2 2 2 2 2 2 3 3 3 ...
#>  $ ntumors     : int  0 0 0 0 0 0 0 0 2 0 ...
#>  $ nmorphine   : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ remission   : int  0 0 0 0 0 0 0 0 0 0 ...
#>   [list output truncated]

The base readers and their defaults:

Function	Default separator	Notes
`read.table()`	any whitespace	most general; `header = FALSE` by default
`read.csv()`	`,`	`header = TRUE` by default
`read.delim()`	tab	`header = TRUE` by default
`read.fwf()`	fixed widths	supply `widths =`

Useful arguments: skip = (lines to ignore), comment.char = "" (stop treating # as a comment), na.strings = (which tokens mean missing), colClasses = (declare types up front, faster and prevents surprise coercion).

readr::read_csv() is faster on large files, never converts strings to factors, and reports a column-type summary you should actually read. rio dispatches on file extension across dozens of formats, which is convenient for exploration.

8.2 Native R formats

tmp_rds  <- file.path(tempdir(), "one_object.rds")
tmp_rda  <- file.path(tempdir(), "many_objects.RData")

x <- seq(1, 10, by = 0.5)
y <- list(a = 1, b = TRUE, c = "oops")

saveRDS(x, tmp_rds)              # ONE object; you name it on read
x2 <- readRDS(tmp_rds)
identical(x, x2)

#> [1] TRUE

save(x, y, file = tmp_rda)       # MANY objects, with their names baked in
load(tmp_rda)                    # restores `x` and `y` into the environment

Prefer saveRDS()/readRDS(). load() silently overwrites existing objects with whatever names the file happens to carry, an unpleasant surprise in a long session. readRDS() returns a value you assign yourself.

8.3 Built-in datasets

data("iris")
head(iris, 3)

data()                              # everything available
data(package = "datasets")          # from one package

8.4 Writing data

out_csv <- file.path(tempdir(), "iris_head.csv")
write.csv(head(iris), out_csv, row.names = FALSE)
readLines(out_csv, n = 2)

#> [1] "\"Sepal.Length\",\"Sepal.Width\",\"Petal.Length\",\"Petal.Width\",\"Species\""
#> [2] "5.1,3.5,1.4,0.2,\"setosa\""

Pitfall. write.table(x, file, sep = ", "), comma and space, is not CSV. Every field after the first acquires a leading space, so numeric columns come back as character on re-import. Use sep = ",", or just write.csv().

row.names = FALSE is almost always what you want; otherwise R writes an unnamed first column of row labels that becomes a mysterious X column when someone reads the file back.

8.5 Writing to the cloud

Web services distinguish the URL that displays a file from the URL that serves its bytes. Google Drive, for example, uses .../open?id=<ID> for viewing and .../uc?export=download&id=<ID> for download.

Do not build a workflow on this. Drive interposes a virus-scan interstitial for larger files, so the download URL returns an HTML page instead of your data, and read.csv() will cheerfully parse that HTML into a nonsensical data frame rather than erroring. For anything that must keep working, use a repository with a persistent identifier (Zenodo, Dataverse, OSF) and record the checksum (§1.6.7).

9 Tidy data and reshaping

9.1 The two shapes

Nearly all tabular data can be arranged two ways:

Wide, one row per subject; repeated measures spread across columns (score_t1, score_t2, score_t3).
Long, one row per subject-per-measurement; a time column and a value column.

Tidy data (Wickham, 2014) means each variable is a column, each observation a row, each type of observational unit a table. Long format is usually tidy; most statistical modeling and plotting functions expect it. Wide format is usually what people hand you.

Reshaping between them is bijective, no information is lost, which is worth demonstrating rather than asserting.

raw_wide <- read.table(header = TRUE, text = '
 CaseID Gender Age Condition1 Condition2
      1      M   5       13.0       10.5
      2      F   6       16.0       11.2
      3      F   8       10.0       18.3
      4      M   9        9.5       18.1
      5      M  10       12.1       19.0
')
raw_wide

library(tidyr)

data_long <- pivot_longer(
  raw_wide,
  cols      = c(Age, Condition1, Condition2),
  names_to  = "Feature",
  values_to = "Measurement"
)
head(data_long, 6)

data_back <- pivot_wider(
  data_long,
  names_from  = Feature,
  values_from = Measurement
)

# Reorder columns to match the original, then compare.
data_back <- as.data.frame(data_back)[, names(raw_wide)]
all.equal(raw_wide, data_back)

#> [1] TRUE

all.equal() confirming TRUE is the demonstration: the transformation is invertible.

Historical note. Note that some older utilities, e.g., reshape2::melt() and reshape2::dcast(), along with the entire reshape2 package, are now retired by their developers in favor of the newer tidyverse, e.g., tidyr package. Older code may still involve melt(), which now maps onto pivot_longer(). Base R’s reshape() also exists and does the same job with a more difficult interface. DSPA3 uses pivot_longer() / pivot_wider() throughout.

9.2 A real example

The SOCR Housing Price Index dataset holds state-by-year housing price index, unemployment rate, region, population, and percent change, natively long.

library(rvest)
library(dplyr)

page <- read_html(
  "https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_010309_HousingPriceIndex"
)
housing_long <- (page |> html_table())[[1]] |>
  mutate(
    Year    = as.integer(Year),
    HPI     = as.numeric(HPI),
    UR      = as.numeric(UR),
    Pop     = as.integer(Pop),
    Percent = as.numeric(Percent)
  )

dim(housing_long)

#> [1] 357   7

head(housing_long, 4)

housing_wide <- housing_long |>
  pivot_wider(
    names_from  = Year,
    values_from = c(HPI, UR, Region, Pop, Percent)
  )
dim(housing_wide)      # one row per state, year-tagged columns

#> [1] 51 36

housing_wide[1:3, 1:6]

housing_back <- housing_wide |>
  pivot_longer(
    cols      = -State,
    names_to  = c(".value", "Year"),
    names_sep = "_"
  ) |>
  mutate(Year = as.integer(Year))

# Same information, recovered.
dim(housing_back)

#> [1] 357   7

identical(sort(names(housing_long)), sort(names(housing_back)))

#> [1] TRUE

The ".value" sentinel in names_to tells pivot_longer() that part of each column name identifies which variable the value belongs to, and the rest is data. It is the mechanism that makes the round trip work when several variables were widened at once.

Pitfall. names_sep = "_" splits on the first underscore. If any of your value columns contains an underscore in its own name (birth_year), the split lands in the wrong place. Use names_pattern = with an explicit regular expression when names are not simple.

10 Mathematics, statistics, and distributions

Reference tables of every math and statistics function live in Appendix A. This section covers the ideas you need to read the rest of the book.

10.1 Descriptive statistics

x <- c(2, 4, 4, 4, 5, 5, 7, 9)

c(mean = mean(x), median = median(x), sd = sd(x), var = var(x),
  IQR = IQR(x), mad = mad(x))

#>     mean   median       sd      var      IQR      mad 
#> 5.000000 4.500000 2.138090 4.571429 1.500000 0.741300

quantile(x)                                  # 0, 25, 50, 75, 100%

#>   0%  25%  50%  75% 100% 
#>  2.0  4.0  4.5  5.5  9.0

quantile(x, probs = c(0.05, 0.95))

#>  5% 95% 
#> 2.7 8.3

range(x); diff(range(x))

#> [1] 2 9

#> [1] 7

scale(x)[, 1] |> round(3)                    # z-scores

#> [1] -1.403 -0.468 -0.468 -0.468  0.000  0.000  0.935  1.871

mean and sd summarize a distribution well only when it is roughly symmetric and light-tailed. For skewed or heavy-tailed data, which is most clinical measurement, median and IQR, or mad, describe the data more faithfully. Report both when you are unsure; the disagreement between them is itself informative.

10.2 Correlation and covariance

cor(iris$Sepal.Length, iris$Petal.Length)                       # Pearson

#> [1] 0.8717538

cor(iris$Sepal.Length, iris$Petal.Length, method = "spearman")  # rank-based

#> [1] 0.8818981

round(cor(iris[, 1:4]), 3)                                      # matrix

#>              Sepal.Length Sepal.Width Petal.Length Petal.Width
#> Sepal.Length        1.000      -0.118        0.872       0.818
#> Sepal.Width        -0.118       1.000       -0.428      -0.366
#> Petal.Length        0.872      -0.428        1.000       0.963
#> Petal.Width         0.818      -0.366        0.963       1.000

Pearson correlation measures linear association only. A perfect deterministic relationship can have a Pearson correlation near zero:

u <- seq(-3, 3, length.out = 200)
round(cor(u, u^2), 10)      # exactly determined, essentially uncorrelated

#> [1] 0

Always plot before trusting a correlation coefficient. Anscombe’s quartet (datasets::anscombe) exists precisely to make this point and is worth five minutes.

10.3 Probability distributions

R names distribution functions systematically. For any distribution dist, four functions exist:

Prefix	Returns	Normal example
`d`	density (or mass) $f(x)$	`dnorm(x)`
`p`	cumulative probability $F(x) = P(X \le x)$	`pnorm(x)`
`q`	quantile $F^{-1}(p)$	`qnorm(p)`
`r`	random draws	`rnorm(n)`

dnorm(0)                    # density of N(0,1) at 0 = 1/sqrt(2*pi)

#> [1] 0.3989423

pnorm(1.96)                 # P(Z <= 1.96)

#> [1] 0.9750021

qnorm(0.975)                # the 97.5th percentile -> 1.96

#> [1] 1.959964

set.seed(3); rnorm(4)

#> [1] -0.9619334 -0.2925257  0.2587882 -1.1521319

The same scheme covers the whole family:

set.seed(3)
c(binomial   = rbinom(1, size = 10, prob = 0.3),
  poisson    = rpois(1, lambda = 4),
  uniform    = round(runif(1), 3),
  exponential= round(rexp(1, rate = 2), 3),
  chisq      = round(rchisq(1, df = 3), 3),
  t          = round(rt(1, df = 10), 3))

#>    binomial     poisson     uniform exponential       chisq           t 
#>       2.000       6.000       0.385       0.502       2.551      -1.161

Also available: beta, gamma, weibull, cauchy, f, geom, hyper, logis, lnorm, nbinom, wilcox, signrank. Each takes its own parameters, check ?rgamma before assuming scale versus rate, a long-standing source of silently wrong simulations.

Further exploration: the Probability Distributome Project, the SOCR distribution calculators, and the bivariate/trivariate interactive calculators.

10.4 Hypothesis tests and model fitting

set.seed(5)
g1 <- rnorm(30, mean = 10, sd = 2)
g2 <- rnorm(30, mean = 11, sd = 2)

t.test(g1, g2)

#> 
#>  Welch Two Sample t-test
#> 
#> data:  g1 and g2
#> t = -2.4553, df = 57.81, p-value = 0.01711
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -2.3505675 -0.2391239
#> sample estimates:
#> mean of x mean of y 
#>  10.02260  11.31744

fit <- lm(Sepal.Length ~ Petal.Length + Species, data = iris)
summary(fit)$coefficients |> round(4)

#>                   Estimate Std. Error t value Pr(>|t|)
#> (Intercept)         3.6835     0.1061 34.7188        0
#> Petal.Length        0.9046     0.0648 13.9624        0
#> Speciesversicolor  -1.6010     0.1935 -8.2752        0
#> Speciesvirginica   -2.1177     0.2735 -7.7439        0

c(AIC = AIC(fit), BIC = BIC(fit), df_resid = df.residual(fit))

#>      AIC      BIC df_resid 
#> 106.2327 121.2859 146.0000

head(round(fitted(fit), 3), 4)

#>     1     2     3     4 
#> 4.950 4.950 4.859 5.040

head(round(residuals(fit), 3), 4)

#>      1      2      3      4 
#>  0.150 -0.050 -0.159 -0.440

Common generics that work across model classes: predict(), coef(), fitted(), residuals(), deviance(), logLik(), AIC(), anova(), confint(). This consistency is one of R’s real strengths, swapping lm() for glm() or randomForest() leaves most downstream code intact.

Other model fitters you will meet: glm() (generalized linear models), nls() (nonlinear least squares), loess() (local regression), aov() (ANOVA), optim() and nlm() (general optimization, covered in the optimization chapter).

11 A first look at visualization

Detailed treatment belongs to the visualization chapter; this is enough to read the figures in Chapters 2–4.

11.1 Three plotting systems

R has three graphics ecosystems, and DSPA3 uses each where it fits:

System	Strength	Use it for
base (`plot`, `hist`, `boxplot`)	Zero dependencies, instant	Quick diagnostic looks while working
`ggplot2`	Grammar of graphics; consistent, composable	Nearly all published figures
`plotly`	Interactive: zoom, hover, rotate	Genuinely interactive or 3-D figures

plotly figures embedded in a self_contained HTML document inflate file size substantially. Reserve it for cases where interaction adds information.

hist(iris$Sepal.Length, breaks = 15, col = "grey85",
     main = "Sepal length", xlab = "cm")

library(ggplot2)

ggplot(iris, aes(Sepal.Length, Petal.Length, colour = Species)) +
  geom_point(alpha = 0.75, size = 2) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.7) +
  labs(title = "Sepal vs. petal length by species",
       x = "Sepal length (cm)", y = "Petal length (cm)", colour = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

A ggplot2 figure is built from three required pieces, data, an aesthetic mapping (aes()), and at least one geom, plus optional scales, facets, coordinates, and themes. Everything composes with +. Once that structure is clear, unfamiliar plot types become easy to read.

11.2 Quantile-quantile plots

A Q-Q plot compares the quantiles of your data against a reference. If they agree, points fall on the identity line. It is the standard diagnostic for “does this look like the distribution I assumed?” and far more informative than a normality test, which mostly tells you about your sample size.

set.seed(2024)
n  <- 500
p  <- ppoints(n)                 # (i - 1/2)/n plotting positions

qq <- data.frame(
  theoretical = qnorm(p),
  Normal      = quantile(rnorm(n),          probs = p),
  Cauchy      = quantile(rcauchy(n),        probs = p),
  Exponential = quantile(rexp(n, rate = 1), probs = p)
)

qq_long <- tidyr::pivot_longer(qq, cols = -theoretical,
                               names_to = "sample", values_to = "quantile")

ggplot(qq_long, aes(theoretical, quantile)) +
  geom_point(size = 0.8, alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0, colour = "firebrick",
              linetype = "dashed") +
  facet_wrap(~ sample, scales = "free_y") +
  labs(title = "Q-Q plots against the standard normal",
       subtitle = "Dashed line = perfect agreement",
       x = "Theoretical N(0,1) quantiles", y = "Sample quantiles") +
  theme_minimal(base_size = 11)

Read the panels as three distinct diagnoses:

Normal, points hug the line; the assumption is reasonable.
Cauchy, extreme deviation at both tails, roughly linear in the middle: a heavy-tailed distribution. Note the y-axis range. The Cauchy has no finite mean, so any analysis of these data reporting a mean and standard error is reporting numbers that do not exist in the population.
Exponential, systematic curvature: skew, not tail weight.

The shape of the departure tells you what is wrong, which a p-value from a normality test cannot.

12 Simulation primer

12.1 Why simulate

Simulation is not a toy exercise. Four production uses:

Method validation. Generate data where you know the true parameter, then check whether your pipeline recovers it. This is the single most effective way to catch a silent statistical error, including one produced by an AI assistant (§1.6.8).
Power and design. Simulate the study you propose to run, at several sample sizes, and count how often you detect the effect.
Privacy-preserving sharing. Synthetic data matching the marginal and joint structure of a protected dataset can be shared for code development when the real data cannot leave its enclave.
Edge cases. Deliberately construct the pathological inputs, perfect separation, extreme collinearity, 90% missingness, that your code must survive.

12.2 Characterizing the observed data

Simulation starts by measuring what you are imitating. We use the PPMI-derived Parkinson’s dataset (SOCR case study).

ppmi <- dspa_read(
  url  = "https://umich.instructure.com/files/330397/download?download_frd=1",
  name = "05_PPMI_top_UPDRS_Integrated_LongFormat1.csv"
)

dim(ppmi)

#> [1] 1764   31

observed <- c(
  n      = sum(!is.na(ppmi$Age)),
  mean   = mean(ppmi$Age, na.rm = TRUE),
  sd     = sd(ppmi$Age,   na.rm = TRUE),
  median = median(ppmi$Age, na.rm = TRUE),
  min    = min(ppmi$Age,  na.rm = TRUE),
  max    = max(ppmi$Age,  na.rm = TRUE)
)
round(observed, 2)

#>       n    mean      sd  median     min     max 
#> 1764.00   61.07   10.27   62.15   31.18   83.03

# Fuller description (verbose output; run interactively)
Hmisc::describe(ppmi)

12.3 The coercion trap, demonstrated

Before building the simulation, a worked instance of the bug from §1.7.2, because it is the most common way a simulation silently produces garbage.

age <- c(62L, 71L, 55L)
dx  <- c("PD", "HC", "SWEDD")

# WRONG: cbind() on vectors builds a MATRIX, which must hold one type.
bad <- cbind(age, dx)
class(bad)

#> [1] "matrix" "array"

typeof(bad)          # "character" - the ages are now strings

#> [1] "character"

bad

#>      age  dx     
#> [1,] "62" "PD"   
#> [2,] "71" "HC"   
#> [3,] "55" "SWEDD"

mean(bad[, "age"])   # NA, with a warning: you cannot average text

#> [1] NA

# RIGHT: data.frame() keeps each column's type.
good <- data.frame(age = age, dx = dx)
str(good)

#> 'data.frame':    3 obs. of  2 variables:
#>  $ age: int  62 71 55
#>  $ dx : chr  "PD" "HC" "SWEDD"

mean(good$age)

#> [1] 62.66667

The failure is loud here because mean() complains. In a longer pipeline it is silent: the character column sorts lexicographically ("100" < "62"), plots as a factor, and joins on string equality. Run str() on every object you build. It costs one line and catches this class of bug immediately.

12.4 Building a synthetic longitudinal cohort

We now simulate a PPMI-like study: 282 subjects across three diagnostic groups, each measured at four visits, with imaging biomarkers, genetics, demographics, and longitudinally evolving clinical scores.

set.seed(1234)

n_subj <- 282
visits <- c(0, 6, 12, 18)          # months
group_n <- c(HC = 100, PD = 100, SWEDD = 82)
stopifnot(sum(group_n) == n_subj)

# --- Imaging biomarkers: Poisson counts, one column per ROI measure ---------
# Every biomarker generated here is carried into the final dataset.
roi_lambda <- c(
  L_caudate_area      =   600, L_caudate_volume      =   800,
  R_caudate_area      =   893, R_caudate_volume      =  1000,
  L_putamen_area      =   900, L_putamen_volume      =  1400,
  R_putamen_area      =  1300, R_putamen_volume      =  3000,
  L_hippocampus_area  =  1300, L_hippocampus_volume  =  3200,
  R_hippocampus_area  =  1500, R_hippocampus_volume  =  3800,
  cerebellum_area     = 16700, cerebellum_volume     = 14000,
  L_lingual_area      =  3300, L_lingual_volume      = 11000,
  R_lingual_area      =  3300, R_lingual_volume      = 12000,
  L_fusiform_area     =  3600, L_fusiform_volume     = 11000,
  R_fusiform_area     =  3300, R_fusiform_volume     = 10000
)

imaging <- as.data.frame(
  lapply(roi_lambda, \(lambda) rpois(n_subj, lambda))
)
dim(imaging)          # 282 x 22 - all 22 biomarkers retained

#> [1] 282  22

# Distinct seed per chunk: each chunk is then reproducible on its own, without
# depending on how many draws an earlier chunk happened to consume.
set.seed(1235)

# --- Subject-level demographics, diagnosis, genetics -----------------------
dx <- factor(rep(names(group_n), times = group_n),
             levels = c("HC", "PD", "SWEDD"))

# Genotype risk differs by group (0/1 minor-allele carrier indicator).
carrier_prob_chr12 <- c(HC = 0.30, PD = 0.60, SWEDD = 0.40)
carrier_prob_chr17 <- c(HC = 0.70, PD = 0.40, SWEDD = 0.50)

subjects <- data.frame(
  case_id             = sprintf("SIM%04d", seq_len(n_subj)),
  dx                  = dx,
  sex                 = factor(sample(c("F", "M"), n_subj, replace = TRUE),
                               levels = c("F", "M")),
  age                 = as.integer(rnorm(n_subj, mean = 62, sd = 10)),
  weight_kg           = as.integer(rnorm(n_subj, mean = 80, sd = 10)),
  chr12_rs34637584_GT = rbinom(n_subj, 1, carrier_prob_chr12[as.character(dx)]),
  chr17_rs11868035_GT = rbinom(n_subj, 1, carrier_prob_chr17[as.character(dx)])
)

# cbind() on two DATA FRAMES preserves types - contrast with Section 1.12.3.
subjects <- cbind(subjects, imaging)
str(subjects[, 1:8])

#> 'data.frame':    282 obs. of  8 variables:
#>  $ case_id            : chr  "SIM0001" "SIM0002" "SIM0003" "SIM0004" ...
#>  $ dx                 : Factor w/ 3 levels "HC","PD","SWEDD": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ sex                : Factor w/ 2 levels "F","M": 1 1 2 2 1 1 1 2 2 2 ...
#>  $ age                : int  55 83 64 29 63 66 68 78 59 65 ...
#>  $ weight_kg          : int  83 84 96 73 87 77 87 85 95 77 ...
#>  $ chr12_rs34637584_GT: int  1 1 1 1 1 0 0 0 1 0 ...
#>  $ chr17_rs11868035_GT: int  1 1 1 0 0 0 1 1 0 0 ...
#>  $ L_caudate_area     : int  570 607 608 582 612 585 624 584 588 557 ...

set.seed(1236)

# --- Expand to one row per subject-visit -----------------------------------
pd_sim <- subjects[rep(seq_len(n_subj), each = length(visits)), ]
pd_sim$visit_month <- rep(visits, times = n_subj)
rownames(pd_sim) <- NULL

# --- Clinical scores that actually evolve over time ------------------------
# Baseline severity and progression slope both depend on diagnostic group.
baseline_iii <- c(HC = 3,   PD = 22,  SWEDD = 12)
slope_iii    <- c(HC = 0.0, PD = 0.45, SWEDD = 0.15)   # points per month

g <- as.character(pd_sim$dx)
pd_sim$UPDRS_III <- pmax(0, round(
  baseline_iii[g] + slope_iii[g] * pd_sim$visit_month +
    rnorm(nrow(pd_sim), mean = 0, sd = 3)
))

baseline_ii <- c(HC = 1,   PD = 11,   SWEDD = 6)
slope_ii    <- c(HC = 0.0, PD = 0.20, SWEDD = 0.08)
pd_sim$UPDRS_II <- pmax(0, round(
  baseline_ii[g] + slope_ii[g] * pd_sim$visit_month +
    rnorm(nrow(pd_sim), mean = 0, sd = 2)
))

pd_sim$UPDRS_I <- rbinom(nrow(pd_sim), size = 2,
                         prob = c(HC = 0.15, PD = 0.35, SWEDD = 0.25)[g])

dim(pd_sim)

#> [1] 1128   33

str(pd_sim[, c("case_id", "dx", "sex", "age", "visit_month",
               "UPDRS_I", "UPDRS_II", "UPDRS_III")])

#> 'data.frame':    1128 obs. of  8 variables:
#>  $ case_id    : chr  "SIM0001" "SIM0001" "SIM0001" "SIM0001" ...
#>  $ dx         : Factor w/ 3 levels "HC","PD","SWEDD": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ sex        : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 2 2 ...
#>  $ age        : int  55 55 55 55 83 83 83 83 64 64 ...
#>  $ visit_month: num  0 6 12 18 0 6 12 18 0 6 ...
#>  $ UPDRS_I    : int  0 1 0 1 0 0 0 0 0 1 ...
#>  $ UPDRS_II   : num  2 0 2 0 1 0 3 0 3 0 ...
#>  $ UPDRS_III  : num  0 2 2 6 0 4 0 7 0 5 ...

Note what this simulation does that a simpler one would not: the clinical scores change over time, at group-specific rates. Note that repeating a single baseline value across all four visits produces data with zero within-subject variation. This is clearly useless for demonstrating any longitudinal method. When using synthetic data to test a mixed model or a GEE, the simulation must contain the structure the model is meant to find.

12.5 Assessing simulation fidelity

Generating data is easy. Generating data that resembles the target is the actual task, and it requires checking.

sim_age <- subjects$age          # subject-level, not row-duplicated

comparison <- rbind(
  observed  = c(n = length(na.omit(ppmi$Age)),
                mean = mean(ppmi$Age, na.rm = TRUE),
                sd   = sd(ppmi$Age,   na.rm = TRUE),
                q25  = unname(quantile(ppmi$Age, 0.25, na.rm = TRUE)),
                q75  = unname(quantile(ppmi$Age, 0.75, na.rm = TRUE))),
  simulated = c(n = length(sim_age), mean = mean(sim_age),
                sd = sd(sim_age),
                q25 = unname(quantile(sim_age, 0.25)),
                q75 = unname(quantile(sim_age, 0.75)))
)
round(comparison, 2)

#>              n  mean    sd   q25   q75
#> observed  1764 61.07 10.27 54.07 68.82
#> simulated  282 61.53  9.92 55.00 68.00

fidelity <- rbind(
  data.frame(age = as.numeric(na.omit(ppmi$Age)), source = "Observed (PPMI)"),
  data.frame(age = sim_age,                       source = "Simulated")
)

ggplot(fidelity, aes(age, fill = source, colour = source)) +
  geom_density(alpha = 0.35, linewidth = 0.8) +
  labs(title = "Simulation fidelity: age distribution",
       subtitle = "Marginal agreement only - joint structure is not shown here",
       x = "Age (years)", y = "Density", fill = NULL, colour = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

# A formal comparison of the two distributions.
suppressWarnings(ks.test(sim_age, as.numeric(na.omit(ppmi$Age))))

#> 
#>  Asymptotic two-sample Kolmogorov-Smirnov test
#> 
#> data:  sim_age and as.numeric(na.omit(ppmi$Age))
#> D = 0.05131, p-value = 0.544
#> alternative hypothesis: two-sided

Read that result carefully. A large p-value does not establish that the simulated data match the observed data; it establishes only that this test failed to detect a difference at this sample size. Failing to reject is not evidence of equivalence, a distinction worth internalizing now, because it recurs throughout model validation.

More importantly: matching marginals is not matching the joint distribution. Our simulation draws age independently of every biomarker, so any real correlation between age and hippocampal volume is absent by construction. Synthetic data built this way is suitable for testing that code runs and for teaching; it is not suitable for drawing scientific conclusions, and it should never be presented as though it were.

12.6 Exporting the simulated cohort

out_path <- file.path(tempdir(), "pd_simulated_cohort.csv")
write.csv(pd_sim, out_path, row.names = FALSE)

cat("Rows:",    nrow(pd_sim),
    "\nCols:",  ncol(pd_sim),
    "\nBytes:", file.size(out_path),
    "\nMD5:",   unname(tools::md5sum(out_path)), "\n")

#> Rows: 1128 
#> Cols: 33 
#> Bytes: 171348 
#> MD5: b94f78fb020db2d527719b2c4cf30729

Writing to tempdir() keeps the repository clean (§1.6.6). The checksum lets you prove later which version of the synthetic cohort produced a given result, the same provenance discipline applied to real data in §1.6.7.

13 Common pitfalls

A consolidated list of the traps introduced above, plus a few more. Most bugs you will hit in your first year of R are on this list.

#	Pitfall	Symptom	Fix
1	`T`/`F` instead of `TRUE`/`FALSE`	Code that worked yesterday now behaves oppositely	Always spell them out
2	`cbind()` on mixed-type vectors	Numbers become strings; `mean()` returns `NA`	Use `data.frame()`
3	`prop.table(v)` instead of `prop.table(table(v))`	Plausible numbers that are not proportions	Tabulate first
4	`1:length(x)` on an empty vector	Loop body runs twice on nothing	`seq_along(x)`
5	`as.numeric(factor)`	Returns 1,2,3… instead of the labels	`as.numeric(as.character(f))`
6	Unset factor `levels`	Reference category is alphabetical, not intended	Set `levels =` explicitly
7	Partial argument matching (`rnorm(n, m = 5)`)	Works until a new argument shadows it	Spell arguments in full
8	No `set.seed()`	Results change every run; prose contradicts output	Seed every stochastic chunk
9	`setwd()` in a script	Runs on one machine only	Projects + `here()`
10	Writing output to the working directory	Knitting litters the repo	`tempfile()` / `tempdir()`
11	`sep = ", "` in `write.table()`	Columns re-import as character	`sep = ","` or `write.csv()`
12	`load()` over an existing object	Silent overwrite of a variable you needed	`saveRDS()`/`readRDS()`
13	`na.omit()` without reporting	Silent, possibly biased, sample-size reduction	Report rows dropped; consider imputation
14	`library(base)`	Nothing; implies otherwise	Delete it
15	`=` vs `==`	Accidental assignment inside a condition	R errors on `if (x = 1)`, but check `subset()` calls
16	Floating-point equality	`0.1 + 0.2 == 0.3` is `FALSE`	`all.equal()` or `isTRUE(all.equal())`
17	`sapply()` returning a list	Downstream code breaks far from the cause	`vapply()` with an explicit template
18	Trusting a correlation without a plot	Nonlinear structure invisible	Plot first, always

0.1 + 0.2 == 0.3                     # FALSE - binary floating point

#> [1] FALSE

all.equal(0.1 + 0.2, 0.3)            # TRUE  - tolerance-aware

#> [1] TRUE

isTRUE(all.equal(0.1 + 0.2, 0.3))    # safe in an `if`

#> [1] TRUE

14 Practice problems

Attempt each before opening the solution. The point is the attempt.

14.1 Problem 1, Predict the type

Without running the code, predict typeof() for each expression. Then verify.

a <- c(1L, 2L, 3L)
b <- c(1L, 2.5)
d <- c(TRUE, FALSE, 1L)
e <- c(1, "2", TRUE)
f <- list(1L, "2", TRUE)
g <- c(list(1), 2)

Solution

a <- c(1L, 2L, 3L);     typeof(a)   # integer

#> [1] "integer"

b <- c(1L, 2.5);        typeof(b)   # double  - integer promotes to double

#> [1] "double"

d <- c(TRUE, FALSE, 1L);typeof(d)   # integer - logical promotes to integer

#> [1] "integer"

e <- c(1, "2", TRUE);   typeof(e)   # character - everything -> string

#> [1] "character"

f <- list(1L, "2", TRUE); typeof(f) # list - no coercion, lists are general

#> [1] "list"

g <- c(list(1), 2);     typeof(g)   # list - c() with a list gives a list

#> [1] "list"

The hierarchy is logical → integer → double → character, and any element of the higher type pulls the whole atomic vector up. Lists escape this entirely, which is why data frames (lists of columns) can mix types while matrices (dimensioned vectors) cannot.

14.2 Problem 2, Find a defensible generalization

Using the patients data from §1.5.6, find a generalization of the quasi-identifiers that achieves $k \ge 5$. Then check whether it is vulnerable to a homogeneity attack, that is, whether any equivalence class has a single distinct diagnosis.

Solution

set.seed(2026)
n <- 5000
patients <- data.frame(
  zip5       = sample(48103:48109, n, replace = TRUE),
  birth_year = sample(1935:2005,   n, replace = TRUE),
  sex        = sample(c("F", "M"), n, replace = TRUE),
  diagnosis  = sample(c("HC", "PD", "SWEDD"), n, replace = TRUE,
                      prob = c(0.5, 0.35, 0.15))
)
patients$zip3          <- patients$zip5 %/% 100
patients$birth_decade  <- (patients$birth_year %/% 10) * 10

candidates <- list(
  A = c("zip3", "birth_decade", "sex"),
  B = c("zip3", "birth_decade"),
  C = c("birth_decade", "sex")
)

for (nm in names(candidates)) {
  qi   <- candidates[[nm]]
  key  <- do.call(paste, c(patients[qi], sep = "\r"))
  size <- as.integer(table(key)[key])
  # Homogeneity: distinct diagnoses within each equivalence class
  ndx  <- tapply(patients$diagnosis, key, \(z) length(unique(z)))
  cat(sprintf("%s: k = %4d | min distinct diagnoses per class = %d\n",
              nm, min(size), as.integer(min(ndx))))
}

#> A: k =  155 | min distinct diagnoses per class = 3
#> B: k =  347 | min distinct diagnoses per class = 3
#> C: k =  155 | min distinct diagnoses per class = 3

Design A already reaches $k \ge 5$ here because the simulated population is uniform over a small set of ZIPs. Real populations are not uniform, and sparse cells appear wherever a rare combination exists, which is exactly where the vulnerable individuals are. Note also that even at $k \ge 5$, a class whose members all share one diagnosis leaks that diagnosis without identifying anyone. That is the gap $\ell$-diversity is designed to close.

14.3 Problem 3, Round-trip a messy table

The frame below is wide with an inconsistent naming scheme. Reshape it to long with separate measure and timepoint columns, then reshape it back and prove the round trip is lossless.

messy <- data.frame(
  subject   = c("s1", "s2", "s3"),
  bp_pre    = c(140, 135, 150),
  bp_post   = c(128, 130, 141),
  hr_pre    = c( 88,  76,  92),
  hr_post   = c( 74,  71,  80)
)

Solution

library(tidyr)

messy <- data.frame(
  subject = c("s1", "s2", "s3"),
  bp_pre  = c(140, 135, 150), bp_post  = c(128, 130, 141),
  hr_pre  = c( 88,  76,  92), hr_post  = c( 74,  71,  80)
)

tidy <- pivot_longer(
  messy,
  cols          = -subject,
  names_to      = c("measure", "timepoint"),
  names_sep     = "_",
  values_to     = "value"
)
tidy

back <- pivot_wider(
  tidy,
  names_from  = c(measure, timepoint),
  values_from = value,
  names_sep   = "_"
)
back <- as.data.frame(back)[, names(messy)]

all.equal(messy, back)

#> [1] TRUE

names_to = c("measure", "timepoint") with names_sep = "_" splits each column name into two variables at once. Had a measure been called resp_rate_pre, the split would land in the wrong place, that is when you switch to names_pattern = "^(.*)_(pre|post)$".

14.4 Problem 4, Estimate power by simulation

Write a function that simulates a two-group comparison ($n$ per group, true mean difference $\delta$, common SD $\sigma$), runs a two-sample $t$-test, and returns whether $p < 0.05$. Use it to estimate power at $n = 20, 40, 80$ for $\delta = 0.5$, $\sigma = 1$.

Solution

simulate_one <- function(n, delta, sigma = 1, alpha = 0.05) {
  g1 <- rnorm(n, mean = 0,     sd = sigma)
  g2 <- rnorm(n, mean = delta, sd = sigma)
  t.test(g1, g2)$p.value < alpha
}

estimate_power <- function(n, delta, sigma = 1, reps = 2000, seed = 99) {
  set.seed(seed)
  mean(replicate(reps, simulate_one(n, delta, sigma)))
}

powers <- vapply(c(20, 40, 80),
                 \(n) estimate_power(n, delta = 0.5),
                 numeric(1))
names(powers) <- paste0("n=", c(20, 40, 80))
round(powers, 3)

#>  n=20  n=40  n=80 
#> 0.342 0.594 0.872

# Cross-check against the closed-form calculation
vapply(c(20, 40, 80),
       \(n) power.t.test(n = n, delta = 0.5, sd = 1)$power,
       numeric(1)) |> round(3)

#> [1] 0.338 0.598 0.882

The two approaches agree to within Monte Carlo error, which validates the simulation. That agreement is the point: once you trust the machinery on a problem with a known answer, you can apply it to designs, clustered, unequal variance, non-normal, with dropout, where no closed form exists. This is exactly the “test on data with a known answer” discipline from §1.6.8.

14.5 Problem 5, Build a provenance record

Extend dspa_read() into dspa_read_provenance() that returns both the data and a metadata record: source URL, retrieval timestamp, file size, MD5 checksum, and dimensions. Explain why a Methods section should contain this.

Solution

dspa_read_provenance <- function(url, name, reader = utils::read.csv, ...) {
  dat  <- dspa_read(url, name, reader = reader, ...)
  path <- file.path(dspa_cache_dir(), name)

  list(
    data = dat,
    provenance = list(
      source_url     = url,
      local_name     = name,
      retrieved_at   = format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z"),
      bytes          = unname(file.size(path)),
      md5            = unname(tools::md5sum(path)),
      n_rows         = nrow(dat),
      n_cols         = ncol(dat),
      r_version      = R.version.string
    )
  )
}

rec <- dspa_read_provenance(
  "https://umich.instructure.com/files/330397/download?download_frd=1",
  "05_PPMI_top_UPDRS_Integrated_LongFormat1.csv"
)
str(rec$provenance)

#> List of 8
#>  $ source_url  : chr "https://umich.instructure.com/files/330397/download?download_frd=1"
#>  $ local_name  : chr "05_PPMI_top_UPDRS_Integrated_LongFormat1.csv"
#>  $ retrieved_at: chr "2026-08-25 12:27:58 EDT"
#>  $ bytes       : num 357211
#>  $ md5         : chr "03e27492af90cc68f0a8497483d183df"
#>  $ n_rows      : int 1764
#>  $ n_cols      : int 31
#>  $ r_version   : chr "R version 4.3.3 (2024-02-29 ucrt)"

Why it belongs in Methods: a URL alone identifies a location, not content. Files at stable URLs get corrected, re-exported, and re-uploaded. Absent a checksum, a reader who downloads “the same” dataset in three years and gets different numbers has no way to tell whether you made an error or the file changed. The checksum converts “we used the PPMI data” into a falsifiable claim. This is the mechanical form of the provenance obligation in §1.5.5.

14.6 Problem 6, Critique a claim

You are reviewing a manuscript containing this sentence:

“We selected Julia because it is 2,941 times faster than R (see Table 2).”

Table 2 has no source, no workload description, and no hardware details. Write a three-sentence reviewer comment. Then state what evidence would justify the sentence.

Solution

A defensible comment:

Table 2 reports speed ratios without specifying the benchmark workload, hardware, language versions, or BLAS configuration, so the reported factor cannot be reproduced or interpreted. Relative language performance is highly workload-dependent, vectorized R and optimized Julia often differ by less than an order of magnitude on array-heavy tasks, while tight scalar loops can differ by three, so a single scalar ratio is not informative about this study’s computations. Please replace the table with a timing of the actual analysis pipeline on the study hardware, with code and versions deposited so the comparison is reproducible.

Evidence that would justify it: a benchmark of the specific computation the paper performs, on stated hardware, with stated versions, with the implementation in both languages made available, and, because a single timing is a sample of size one, repeated enough times to report variability.

Compare with your own measurement in §1.6.2: the loop-versus- vectorized gap within R exceeded the gap most cross-language tables report between languages. Implementation dominates language choice for nearly all analytic code.

15 Checkpoint

Answer without scrolling back.

Name the three “nothing” values in R and one situation distinguishing each.
Your collaborator sends a script beginning setwd("C:/Users/pat/thesis"). Name two things that break and the fix for each.
k_anonymity() reports $k = 1$ for 340 of 5,000 records. What does that mean operationally, and what are two ways to address it?
Why does DSPA3 require set.seed() in every stochastic chunk rather than once at the top of a file?
A Q-Q plot against the normal shows points on the line in the middle and sharply off it at both ends. What does that indicate, and name one summary statistic that becomes untrustworthy.
Distinguish prediction, explanation, and causal inference in one sentence each, using the same clinical example for all three.

Answers

NA, a value exists but is unknown (an unrecorded lab result). NULL, there is no value at all (an absent list element; length 0). NaN, an undefined numeric result (0/0).
The path does not exist on any other machine, and the script’s behavior depends on the session’s starting directory. Fix: an RStudio/Positron project plus here() for all paths; pin dependencies with renv so package versions travel with the code.
340 records are unique on their quasi-identifiers and can in principle be matched against any external file containing the same variables. Address by generalizing (coarsen ZIP, bin age), suppressing those records, or moving to a formal framework such as differential privacy, and note that raising $k$ alone does not defeat a homogeneity attack.
Because chunks are re-run, cached, and reordered during editing. A single top-of-file seed only reproduces the output if every chunk executes exactly once, in order, in a fresh session, which is not how documents get written.
Heavy tails. The sample standard deviation (and therefore every standard error, confidence interval, and t-statistic derived from it) becomes unreliable; in the limiting Cauchy case the population mean does not exist.
Prediction: which patients will be readmitted within 30 days. Explanation: which factors are associated with readmission, and how strongly, adjusting for measured covariates. Causal inference: whether assigning a patient to the discharge-planning program would reduce their readmission probability.

16 Summary and what comes next

What this chapter established

Health data is hard for identifiable, nameable reasons, size, complexity, incongruency, multi-source, multi-scale, time, incompleteness, and provenance , and each names a class of tool the book will supply.
Prediction, explanation, and causal inference are different goals requiring different methods. Confusing them produces confident, useless results.
Ethics is measurable. Re-identification risk is a computation, not a sentiment, and privacy trades against utility along a curve someone must choose a point on.
Reproducibility is mechanics: projects, seeds, pinned dependencies, cached and checksummed data, recorded sessions.
R’s core structures are atomic vectors, lists, and the data frames built from them. Type coercion explains most confusing R behavior.
Tidy long/wide reshaping is invertible, and long format is what modeling and plotting functions expect.
Simulation validates methods, powers designs, and protects sensitive data, provided you check fidelity and are candid about what the simulation does not reproduce.
Quantitative claims require provenance. Compute what you can compute; cite what you cannot; read what you cite.

Where each thread continues

Thread from this chapter	Continues in
Data import, quality, missingness	Data handling
Visualization beyond first principles	Visual analytics & EDA
Matrices, linear models	Linear algebra & regression
What accuracy and cross-validation mean	Model assessment
Supervised prediction	Classification
High-dimensional structure	Dimensionality reduction
Which variables matter	Feature selection
Repeated measures over time	Longitudinal analysis
Fitting as optimization	Function optimization
Scale, streaming, deployment	Specialized topics

Before moving on, confirm you can: install and attach a package; create a project and read data through dspa_read(); explain what str() tells you about an object; reshape a table both directions; and simulate a dataset with a known answer and verify a method recovers it.

17 Chapter roadmap

Forward references throughout this chapter resolve here. Each entry names the DSPA3 chapter that develops the thread, so you can jump straight to a topic without reading linearly.

Data handling, quality, and missingness, import from every common format, assess data quality, and treat missingness properly (MCAR / MAR / MNAR, imputation, EM). Continues §1.7–§1.9.
Visual analytics and exploratory data analysis, the grammar of graphics in depth, composition/comparison/relationship/distribution displays, and interactive and 3-D visualization. Continues §1.11.
Linear algebra, matrix computing, and regression, matrix operations, eigenspectra, ordinary least squares, multiple regression, regression and model trees. Continues §1.7.5 and §1.10.4.
Model assessment, validation, and improvement, what accuracy, sensitivity, specificity, $\kappa$, ROC/AUC, and cross-validation actually mean, plus bootstrap estimation and hyperparameter tuning. Read this before the modeling chapters.
Supervised classification, $k$-NN, naive Bayes, decision trees, and rule learners, with clinical case studies.
Dimensionality reduction, PCA, ICA, factor analysis, SVD, t-SNE, and UMAP for high-dimensional imaging and genetics data.
Variable importance and feature selection, filter, wrapper, and embedded methods; regularization (ridge, LASSO); false-discovery-rate control and knockoff filtering.
Longitudinal and time-series analysis, ARIMA/ARIMAX, structural equation models, mixed models, GEE, and recurrent architectures for repeated measures. Continues §1.12.4.
Function optimization, gradient descent, Newton and quasi-Newton methods, constrained and Bayesian optimization. Every model fit in this book is an optimization problem underneath.
Specialized topics, specialized formats and databases, data streaming, parallel and out-of-core computing, and deployment.

18 Appendix A: R reference

Lookup material. Skim once; return when you need a name.

18.1 A.1 Mathematical functions

Function	Returns
`abs`, `sqrt`, `exp`, `log(x, base)`, `log10`, `log2`, `log1p`, `expm1`	Elementary functions
`sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `atan2`	Trigonometric
`ceiling`, `floor`, `round(x, n)`, `signif`, `trunc`	Rounding
`choose(n, k)`, `factorial(n)`, `gamma`, `lgamma`, `beta`	Combinatorial / special
`max`, `min`, `range`, `sum`, `prod`, `diff`	Aggregation
`cumsum`, `cumprod`, `cummax`, `cummin`	Cumulative
`pmax(x, y)`, `pmin(x, y)`	Element-wise extremes across vectors
`union`, `intersect`, `setdiff`, `setequal`, `is.element`	Set operations
`Re`, `Im`, `Mod`, `Arg`, `Conj`	Complex numbers
`fft`, `mvfft`, `convolve`, `filter`	Signal processing

Most aggregation functions accept na.rm = TRUE.

18.2 A.2 Statistical summaries

Function	Returns
`mean`, `median`, `weighted.mean(x, w)`	Central tendency
`var`, `sd`, `IQR`, `mad`	Spread
`quantile(x, probs)`, `fivenum`	Order statistics
`cor(x, y, method)`, `cov`	Association
`scale(x, center, scale)`	Standardization
`rank`, `order`, `sort`, `rev`	Ordering
`table`, `prop.table`, `xtabs`, `ftable`	Frequencies
`t.test`, `wilcox.test`, `chisq.test`, `fisher.test`	Common tests
`prop.test`, `binom.test`, `shapiro.test`, `ks.test`	More tests
`aov`, `anova`, `TukeyHSD`	ANOVA
`density(x)`, `ecdf(x)`	Distribution estimates

help.search("test") lists everything installed.

18.3 A.3 Matrix operations

Expression	Meaning
`t(x)`	Transpose
`diag(x)`	Diagonal (or build a diagonal matrix)
`A %*% B`	Matrix product
`A * B`	Element-wise product
`solve(A)`	Inverse
`solve(A, b)`	Solve $Ax = b$ (preferred over `solve(A) %*% b`)
`crossprod(A)`	$A^\top A$, faster than `t(A) %*% A`
`det(A)`, `qr(A)`, `chol(A)`, `svd(A)`, `eigen(A)`	Decompositions
`rowSums`, `colSums`, `rowMeans`, `colMeans`	Fast margins
`nrow`, `ncol`, `dim`, `dimnames`	Shape

A <- cbind(c(1, -1/5), c(-1/3, 1))
b <- c(1, 2)
solve(A, b)                      # solve the system

#> [1] 1.785714 2.357143

round(solve(A) %*% A, 12)        # inverse times original = identity

#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1

Full treatment in the linear algebra chapter.

18.4 A.4 Strings

Function	Purpose
`paste(..., sep)`, `paste0(...)`	Concatenate
`sprintf(fmt, ...)`	Formatted construction
`nchar`, `substr`, `substring`	Length, extraction
`strsplit(x, split, fixed)`	Split
`grep`, `grepl`, `regexpr`, `gregexpr`	Search
`sub`, `gsub`	Replace first / all
`toupper`, `tolower`, `trimws`	Case, whitespace
`startsWith`, `endsWith`	Prefix / suffix tests
`format`, `formatC`, `prettyNum`	Display formatting

a <- "When the going gets tough, the tough get going!"
substr(a, 10, 40)

#> [1] "going gets tough, the tough get"

strsplit("a.b.c", ".", fixed = TRUE)[[1]]

#> [1] "a" "b" "c"

gsub("tough", "TOUGH", a)

#> [1] "When the going gets TOUGH, the TOUGH get going!"

sprintf("%.2f%% complete (%d of %d)", 100 * 7 / 9, 7, 9)

#> [1] "77.78% complete (7 of 9)"

The stringr package offers a more consistent interface (str_detect, str_replace, str_extract) over the same regular-expression engine.

18.5 A.5 Dates and times

Date stores calendar dates; POSIXct stores date-times with time zones.

d <- as.Date("2026-08-23")
d + 30

#> [1] "2026-09-22"

format(d, "%A, %B %d, %Y")

#> [1] "Sunday, August 23, 2026"

seq(d, by = "month", length.out = 4)

#> [1] "2026-08-23" "2026-09-23" "2026-10-23" "2026-11-23"

difftime(as.Date("2026-12-25"), d, units = "days")

#> Time difference of 124 days

Code	Meaning	Code	Meaning
`%Y`	4-digit year	`%M`	Minute (00–59)
`%y`	2-digit year (avoid)	`%S`	Second (00–61, leap seconds)
`%m`	Month (01–12)	`%p`	AM/PM
`%b`, `%B`	Month name, abbrev./full	`%j`	Day of year (001–366)
`%d`	Day of month (01–31)	`%U`, `%W`	Week (Sunday-, Monday-based)
`%a`, `%A`	Weekday name, abbrev./full	`%z`	UTC offset (output)
`%H`, `%I`	Hour (00–23 / 01–12)	`%Z`	Time zone name (output)

See ?strptime. The lubridate package simplifies parsing (ymd(), mdy_hms()) and arithmetic across time zones and daylight-saving transitions, worth it as soon as dates become non-trivial.

18.6 A.6 Base graphics quick reference

High-level (start a plot): plot, hist, boxplot, barplot, pie, dotchart, matplot, pairs, coplot, mosaicplot, assocplot, fourfoldplot, stars, symbols, contour, filled.contour, image, persp, qqnorm, qqplot, curve, sunflowerplot, interaction.plot, ts.plot, plot.ts, termplot.

Low-level (add to an existing plot): points, lines, text, mtext, segments, arrows, abline, rect, polygon, legend, title, axis, rug, grid, box.

Common parameters (via par() or passed directly): main, sub, xlab, ylab, xlim, ylim, col, pch, lty, lwd, cex, font, las, bty, bg, adj, mar, mfrow, mfcol, pty, xaxt, yaxt, tck.

length(colors())          # named colors available

#> [1] 657

head(colors(), 8)

#> [1] "white"         "aliceblue"     "antiquewhite"  "antiquewhite1"
#> [5] "antiquewhite2" "antiquewhite3" "antiquewhite4" "aquamarine"

DSPA3 uses ggplot2 for nearly all published figures; base graphics remain excellent for fast diagnostic looks while you work.

19 Appendix B: Debugging

Programs that produce wrong answers usually contain logical errors, not syntax errors. Debugging is the systematic process of locating them.

An important distinction: some tools require an interactive session and do nothing useful inside a knitted document. The table marks which is which.

Tool	Purpose	Works in a knit?
`traceback()`	Show the call stack after an error	Interactive only
`tryCatch()`	Catch a condition and handle it	✅ Yes
`withCallingHandlers()`	Observe conditions without stopping	✅ Yes
`try()`	Continue past an error	✅ Yes
`stopifnot()`	Assert an invariant	✅ Yes
`browser()`	Pause at a point in a function	Interactive only
`debug()` / `debugonce()`	Step through a function	Interactive only
`trace()`	Inject code into an existing function	Interactive only
`options(error = recover)`	Post-mortem environment inspection	Interactive only

19.1 B.1 Reading an error

f <- function(x) x - g(x)
g <- function(y) y * h(y)
h <- function(z) {
  r <- log(z)          # log of a negative number -> NaN
  if (r < 10) r^2 else r^3
}

f(-1)

#> Error in `if (r < 10) ...`:
#> ! missing value where TRUE/FALSE needed

The error is not in f(). log(-1) returns NaN with a warning, then if (NaN < 10) fails because the condition is missing. The reported location and the actual defect are two calls apart, which is precisely why a call stack matters.

Interactively, run traceback() immediately after the error:

f(-1)
traceback()
#> 3: h(y)
#> 2: g(x)
#> 1: f(-1)

Read it bottom-up: f called g called h, and h is where it broke.

19.2 B.2 Catching conditions

safe_log <- function(z) {
  tryCatch(
    {
      r <- log(z)
      if (is.nan(r)) stop("log() produced NaN for input ", z, call. = FALSE)
      r
    },
    error   = function(e) { message("Handled error: ", conditionMessage(e)); NA_real_ },
    warning = function(w) { message("Handled warning: ", conditionMessage(w)); NA_real_ }
  )
}

safe_log(10)

#> [1] 2.302585

safe_log(-1)

#> [1] NA

Look closely at the second result. log(-1) emits a warning (“NaNs produced”) before our stop() is ever reached, so the warning handler fires and the error branch never runs. tryCatch() responds to the first matching condition and abandons the rest of the block, which is exactly the sort of control-flow detail that turns into a mysterious bug when a handler quietly swallows something you meant to see.

tryCatch() stops execution at the condition and runs your handler. When you want to observe a warning without aborting, use withCallingHandlers() plus invokeRestart("muffleWarning"):

noisy <- function() { warning("something is off"); "finished anyway" }

withCallingHandlers(
  noisy(),
  warning = function(w) {
    message("Logged: ", conditionMessage(w))
    invokeRestart("muffleWarning")
  }
)

#> [1] "finished anyway"

19.3 B.3 Failing fast

The cheapest debugging is the error you never have to find:

compute_bmi <- function(weight_kg, height_m) {
  stopifnot(
    is.numeric(weight_kg), is.numeric(height_m),
    length(weight_kg) == length(height_m),
    all(height_m > 0, na.rm = TRUE)
  )
  weight_kg / height_m^2
}

round(compute_bmi(c(70, 85), c(1.75, 1.80)), 2)

#> [1] 22.86 26.23

compute_bmi(c(70, 85), c(1.75, 0))

#> Error in `compute_bmi()`:
#> ! all(height_m > 0, na.rm = TRUE) is not TRUE

Assertions convert a wrong answer three steps downstream into an immediate, locatable failure. In a data pipeline they are worth writing liberally: check row counts after a join, check that proportions sum to one, check that no identifier is duplicated.

19.4 B.4 Interactive stepping

Run these in the console, not in a knit:

SS <- function(mu, x) {
  d  <- x - mu
  d2 <- d^2
  sum(d2)
}

set.seed(100); x <- rnorm(100)

debugonce(SS)    # step through the next call only
SS(1, x)

Inside the Browse[1]> prompt:

Key	Action
`n`	Execute the current line, print the next
`s`	Step into the function being called
`f`	Finish the current loop or function
`c`	Continue to the end (or next breakpoint)
`Q`	Quit the debugger
`where`	Print the call stack
`ls()`	List objects in the current frame

To pause at a chosen point rather than the start, insert browser(). Guard it so it never fires during a render:

h2 <- function(z) {
  if (interactive()) browser()
  r <- log(z)
  if (r < 10) r^2 else r^3
}

debug(f) flags a function permanently until undebug(f); debugonce(f) is almost always what you want, since a forgotten debug() flag makes later code behave mysteriously.

For post-mortem inspection of a failure deep in a call stack:

options(error = recover)   # choose a frame to inspect when an error occurs
# ... run the failing code ...
options(error = NULL)      # restore default behavior

19.5 B.5 A debugging protocol

Reproduce it, reliably, on the smallest input that still fails.
Read the actual message. Most R errors say what happened; the difficulty is that they say it about the symptom.
Locate it, traceback(), or bisect by commenting out.
Inspect state at the failure, str() on every input.
Form one hypothesis and test it. Changing three things at once means you will not know which one mattered.
Add an assertion so the bug cannot silently return.

Step 6 is the one people skip. A fixed bug without a regression check is a bug waiting to be reintroduced.

20 Appendix C: Further resources

Books (free online)

R for Data Science, 2nd ed., the standard modern introduction.
Advanced R, 2nd ed., R’s object systems, environments, metaprogramming.
R Packages, 2nd ed., when your analysis outgrows a script.
An Introduction to R , the official manual.

Tutorials and courses

Reference

CRAN Task Views, curated packages by domain.
Posit cheatsheets.
Bioconductor, genomics, proteomics, imaging.

SOCR / DSPA

Reproducibility

renv · here · targets for pipeline orchestration.
The Turing Way, a comprehensive handbook on reproducible research practice.

21 Session information

Every DSPA3 chapter records the environment that produced it. Include this in anything you share.

sessionInfo()

#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] dplyr_1.1.4   rvest_1.0.4   tidyr_1.3.1   ggplot2_4.0.1
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.6-5       gtable_0.3.6       jsonlite_1.8.9     compiler_4.3.3    
#>  [5] promises_1.3.2     Rcpp_1.0.14        tidyselect_1.2.1   xml2_1.3.6        
#>  [9] later_1.4.1        jquerylib_0.1.4    splines_4.3.3      scales_1.4.0      
#> [13] yaml_2.3.10        fastmap_1.2.0      lattice_0.22-6     R6_2.6.1          
#> [17] labeling_0.4.3     generics_0.1.3     curl_6.2.0         knitr_1.51        
#> [21] tibble_3.2.1       bslib_0.9.0        pillar_1.10.1      RColorBrewer_1.1-3
#> [25] rlang_1.1.5        websocket_1.4.1    cachem_1.1.0       xfun_0.52         
#> [29] sass_0.4.9         S7_0.2.1           otel_0.2.0         cli_3.6.3         
#> [33] mgcv_1.9-1         withr_3.0.2        magrittr_2.0.3     ps_1.9.0          
#> [37] processx_3.8.6     digest_0.6.37      grid_4.3.3         rstudioapi_0.18.0 
#> [41] nlme_3.1-165       chromote_0.4.0     lifecycle_1.0.5    vctrs_0.6.5       
#> [45] evaluate_1.0.3     glue_1.8.0         farver_2.1.2       rmarkdown_2.31    
#> [49] purrr_1.0.2        httr_1.4.7         tools_4.3.3        pkgconfig_2.0.3   
#> [53] htmltools_0.5.8.1

Prefix	Returns	Normal example
`d`	density (or mass) \(f(x)\)	`dnorm(x)`
`p`	cumulative probability \(F(x) = P(X \le x)\)	`pnorm(x)`
`q`	quantile \(F^{-1}(p)\)	`qnorm(p)`
`r`	random draws	`rnorm(n)`

DSPA3: Data Science and Predictive Analytics (UMich HS650)

Chapter 1: Foundations, Data Science, Ethics, and the R Environment

SOCR/MCAIM (Ivo Dinov)

August 25, 2026

1 Learning objectives

2 What this book is for

2.1 What is expected of you

3 Driving problems

3.1 Neurodegeneration and dementia

3.2 Substance use and withdrawal

3.3 Genomics and high-throughput computing

3.4 Interactive exploration

3.5 Open data repositories

4 What makes health data hard

4.1 Beyond the “V” vocabulary

4.2 Data science and predictive analytics

4.3 Prediction is not explanation is not causation

5 Responsible and reproducible practice

5.1 FAIR resource sharing

5.2 Research ethics

5.3 Benefits and detriments of findings

5.4 Sensitive data: regulation and practice

5.5 Provenance and longevity

5.6 Quantifying re-identification risk

5.6.1 A cautionary tale, read carefully

5.7 When resources are misused

6 The DSPA3 computing environment

6.1 Why R?

6.2 Reading a benchmark critically

6.3 Measuring adoption: a construct-validity lesson

6.4 Installing R

6.4.1 Keeping R current

6.5 Choosing an IDE

6.6 Installing packages

6.6.1 The DSPA3 dependency budget

6.7 Reproducibility

6.7.1 1. Use projects, not setwd()

6.7.2 2. Set seeds, always

6.7.3 3. Never write to the working directory

6.7.4 4. Pin your dependencies with renv

6.7.5 5. Record the session

6.8 Accessing data reproducibly

6.9 Getting help

6.10 Working with an AI coding assistant

6.11 Quarto: authoring reproducible documents

7 R essentials

7.1 Assignment and the two pipes

7.2 Atomic vectors and coercion

7.3 Generating structured data

7.4 Lists

7.5 Matrices, arrays, and tensors

7.6 Data frames and tibbles

7.6.1 The stringsAsFactors story

7.7 Missing values

7.8 Inspecting objects

7.9 Control flow and vectorization

7.10 Writing functions

7.11 Iteration over structures

7.12 Frequencies and cross-tabulation

7.13 Combining and reshaping data frames

7.14 Sampling

8 Data import and export

8.1 Reading tabular data

8.2 Native R formats

8.3 Built-in datasets

8.4 Writing data

8.5 Writing to the cloud

9 Tidy data and reshaping

9.1 The two shapes

9.2 A real example

10 Mathematics, statistics, and distributions

10.1 Descriptive statistics

10.2 Correlation and covariance

10.3 Probability distributions

10.4 Hypothesis tests and model fitting

11 A first look at visualization

11.1 Three plotting systems

11.2 Quantile-quantile plots

12 Simulation primer

12.1 Why simulate

6.4.1 Keeping `R` current

6.7.1 1. Use projects, not `setwd()`

6.7.4 4. Pin your dependencies with `renv`

7.6.1 The `stringsAsFactors` story