This glossary collects the key concepts, people, and ideas you will meet in this lesson. Use it as a reference while you work through the material, or as a review before assessments. Type in the search box to filter entries.

Introduction and Overview

An earlier lesson set up the structured workflow that takes raw data from collection to analysis-ready files. This lesson zooms into two specific phases of that pipeline that deserve their own treatment: cleaning the data (detecting errors, handling missing values, deciding how to deal with outliers) and producing descriptive analyses (the numerical summaries and visualisations that every analysis report eventually needs). The three content sections walk through this in order: data quality and types of error (this section), the cleaning workflow itself (a later section), and descriptive statistics and Table 1 (a later section).

The Data Pipeline

Before any statistical analysis can begin, data must travel through a pipeline: collection → entry → storage → analysis. Errors can be introduced at every stage, and understanding where problems arise is the first step toward ensuring high-quality data for epidemiologic research.

At the collection stage, errors may arise from poorly worded survey items, interviewer bias, instrument malfunction, or participant misunderstanding. During data entry, transcription mistakes and coding errors are common. Storage introduces risks of data corruption, version control problems, and format incompatibilities. Finally, at the analysis stage, incorrect variable coding, inappropriate transformations, and software errors can compromise results.

Types of Data Errors

Errors in epidemiologic datasets take many forms. Van den Broeck, Argeseanu Cunningham, Eeckels, & Herbst (2005) describe data cleaning as a three-stage process, screening, diagnosing, and editing, that should be planned into study design rather than improvised after collection. Click each card below to learn about common error types.

Screening means looking systematically for values that seem odd, such as an age of 220 or an empty cell where an answer was expected. Diagnosing means deciding what each odd value is: an error, a true but unusual value, or a missing-data code. Editing means acting on the diagnosis by correcting the value from the source document, setting it to missing, or leaving it unchanged, and recording the decision. Planning these steps before data collection allows the study team to build some of the checks into the data-entry system itself.

Missing Data

Missing data are ubiquitous in epidemiologic research. The mechanism by which data become missing has profound implications for the validity of subsequent analyses. Rubin (1976) identified three missing data mechanisms:

Data Dictionaries and Codebooks

A data dictionary (or codebook) is a structured document that describes every variable in a dataset: its name, label, type, valid range, coding scheme, and source. Maintaining a comprehensive data dictionary is essential for:

• Ensuring consistency across data entry operators
• Facilitating communication among research team members
• Supporting reproducibility when analyses are revisited months or years later
• Identifying and resolving coding discrepancies

Variable Types in Epidemiologic Data

Correctly classifying variables is fundamental to choosing appropriate cleaning strategies and analytic methods.

The Course Dataset and Its Codebook

The R activities in this lesson use phaa_survey.csv, a simulated cross-sectional survey of 800 adults that later lessons also use. The table below is a codebook excerpt for the variables this lesson works with most. A hard limit marks values that cannot be correct for this survey; values beyond it are treated as errors. A soft limit marks values that are possible but unusual; values between a soft and a hard limit are flagged for investigation, checked against the source where possible, and kept unless they are shown to be wrong. The limits are the rules applied in this lesson.

In this file a missing value is an empty cell, which R reads as NA (not available). Many real files instead store codes such as -9 (refused) or 999 (not asked); the codebook is where such codes are declared, and the next section shows how to convert them before any analysis. A full codebook would also record each variable's source, such as the questionnaire item or the measurement protocol; in this simulated file every variable comes from a single survey record per participant.

Getting Ready to Run R

Reading R Code: A Short Primer

R code is a sequence of instructions. Most lines either create an object (a named result stored in memory) or ask R to display something. The symbols in the table below appear throughout this lesson.

A First Look at the Raw File

# Read the raw file. Blank cells and the text "NA" become missing values (NA).
phaa <- read.csv("phaa_survey.csv", stringsAsFactors = FALSE,
                 na.strings = c("", "NA"))

dim(phaa)                 # number of rows (people) and columns (variables)
str(phaa)                 # each variable's type and its first few values
summary(phaa[, c("age", "systolic_bp", "bmi", "dep1")])
table(phaa$smoker, useNA = "ifany")   # counts for a text (character) variable
colSums(is.na(phaa))      # how many values are missing in each column

● Complete the knowledge check to continue.

Introduction and Overview

An earlier section named the kinds of errors that can enter a dataset. This section provides the cleaning workflow that addresses them: systematic verification, outlier identification, missing-data handling, transformations, recoding, and string cleaning. Each step has its own techniques and its own decisions to document.

Systematic Verification

Data cleaning should follow a systematic, documented process. The goal is not to “fix” data arbitrarily but to identify and resolve genuine errors while preserving the integrity of valid observations.

Identifying Outliers

Outliers are observations that lie far from the bulk of the data. They may represent genuine extreme values, data errors, or observations from a different population. The visual tools used to spot them, including box plots, stem-and-leaf displays, and the broader practice of looking at data before modelling, trace back to Tukey (1977). Several methods are used to identify them:

What the Visual Checks Reveal

Two kinds of problem are easier to see in a plot than in a table of numbers. A histogram (a bar chart of how many values fall in each range) shows gaps and separate clusters. When a decimal point is shifted during data entry, the shifted values form a small cluster at ten times the typical value, far from the main distribution, as in panel B below.

A scatter plot shows each person as a point with two measurements. A bivariate outlier is a combination of values that is implausible even though each value is ordinary on its own. In the simulated example below, a height of 191 cm and a weight of 46 kg each lie inside the boxplot whiskers, but together they imply a BMI of 12.6, which would be life-threatening. Only the scatter plot reveals the problem.

Outlier Rules in R

# Worked example: the eleven readings from the text
x <- c(96, 102, 106, 110, 112, 115, 118, 121, 125, 131, 184)
quantile(x, probs = c(0.25, 0.75))     # Q1 and Q3
IQR(x)                                  # Q3 minus Q1
boxplot.stats(x)$out                    # values beyond the fences

# The same steps on BMI in the course file (before any cleaning)
q     <- quantile(phaa$bmi, probs = c(0.25, 0.75), na.rm = TRUE)
iqr   <- IQR(phaa$bmi, na.rm = TRUE)
lower <- q[1] - 1.5 * iqr
upper <- q[2] + 1.5 * iqr
c(lower, upper)
flag <- which(phaa$bmi < lower | phaa$bmi > upper)
length(flag)                            # how many rows are flagged
phaa[flag, c("id", "bmi")]              # who they are, and their values

# z-scores: scale() computes (x - mean) / SD for every value
z <- as.numeric(scale(phaa$bmi))
phaa$id[which(abs(z) > 3)]              # IDs with |z| greater than 3
round(z[which(phaa$bmi == 31.6)], 2)   # z for the BMI of 31.6

# Masking: repeat without the impossible BMI of 95
bmi_ok <- phaa$bmi[phaa$bmi < 60]
round(c(mean = mean(phaa$bmi), sd = sd(phaa$bmi)), 2)
round(c(mean = mean(bmi_ok), sd = sd(bmi_ok)), 2)
round((31.6 - mean(bmi_ok)) / sd(bmi_ok), 2)

# Soft limits: count values that deserve a second look (nothing is changed)
sum(phaa$systolic_bp < 90 | phaa$systolic_bp > 180)
sum(phaa$bmi < 16 | phaa$bmi > 50)

Range Checks in Practice

The activity below applies the hard limits from the codebook to the course file. It is the editing stage of the screening, diagnosing and editing sequence: the values were screened in the first-look box, diagnosed as impossible in the outlier box, and are now replaced, with each change listed first.

# 1. Read the raw file. Treat blank cells as missing (na.strings).
phaa <- read.csv(file = "phaa_survey.csv",
                 stringsAsFactors = FALSE,
                 na.strings = c("", "NA"))
summary(phaa$age)                 # before the conversion in step 2

# 2. Make sure age and bmi are numeric (in a real export they often arrive
#    as character because a cell holds text such as "unknown"; here every
#    cell is a number, so no NAs are created at this step).
phaa$age <- as.numeric(phaa$age)
phaa$bmi <- as.numeric(phaa$bmi)
summary(phaa$age)                 # after the conversion

# 3a. Count and list the rows that break each hard limit BEFORE changing them
sum(phaa$age < 18 | phaa$age > 100, na.rm = TRUE)
phaa$id[which(phaa$age < 18 | phaa$age > 100)]
sum(phaa$systolic_bp < 80 | phaa$systolic_bp > 220, na.rm = TRUE)
phaa$id[which(phaa$systolic_bp < 80 | phaa$systolic_bp > 220)]
sum(phaa$bmi < 13 | phaa$bmi > 60, na.rm = TRUE)
phaa$id[which(phaa$bmi < 13 | phaa$bmi > 60)]

# 3b. Range checks: replace the impossible values with NA
phaa$age[phaa$age < 18 | phaa$age > 100] <- NA
phaa$systolic_bp[phaa$systolic_bp < 80 | phaa$systolic_bp > 220] <- NA
phaa$bmi[phaa$bmi < 13 | phaa$bmi > 60] <- NA
summary(phaa$age)                 # after the range check

# 4. Put the education levels in order (lowest to highest) so that
#    comparisons such as <= are meaningful
phaa$education <- factor(phaa$education,
  levels = c("Less than high school", "High school", "Some college",
             "Bachelor's", "Graduate degree"),
  ordered = TRUE)
levels(phaa$education)

Handling Missing Data

The choice of method for handling missing data depends on the assumed missing data mechanism and the proportion of missingness. Practical guidance on implementing multiple imputation via chained equations, including how to specify the imputation model, how many imputations to draw, and how to diagnose problems, is given by White, Royston, & Wood (2011).

Missing Data in R: How Much, and Who Is Missing?

# How much is missing, and where? (after the range checks above)
n_miss <- colSums(is.na(phaa))
n_miss[n_miss > 0]                       # only the columns with gaps
round(100 * n_miss[n_miss > 0] / nrow(phaa), 1)   # the same, as percentages
sum(complete.cases(phaa))                # rows with no gap in any column
nrow(na.omit(phaa))                      # na.omit() keeps only those rows

# Do people with and without a recorded income differ on what was observed?
inc_miss <- is.na(phaa$income)           # TRUE = income missing
table(inc_miss)
round(tapply(phaa$age, inc_miss, mean, na.rm = TRUE), 1)
round(tapply(phaa$systolic_bp, inc_miss, mean, na.rm = TRUE), 1)
round(prop.table(table(inc_miss, phaa$smoker), margin = 1) * 100, 1)

Reasoning About the Mechanism

The observed data can show that missingness is related to observed variables, which is evidence against MCAR. They cannot show whether missingness depends on the missing values themselves, because those values were never recorded. Whether MNAR is plausible is therefore judged from subject knowledge. Income is the standard example: refusal to answer income questions is often reported to be more common among people with high incomes, so an analyst should treat MNAR as a real possibility for income even when a comparison like the one above looks reassuring. By contrast, a blood pressure reading lost because a machine failed at random times is a case where MCAR is plausible.

Because the mechanism cannot be confirmed from the data, it is stated as an assumption in the methods section of a report. An example sentence is: "Income was missing for 35 of 800 participants (4.4%). We assumed that income was missing at random given age, gender, education and region, and handled it with multiple imputation; a sensitivity analysis that assumed non-responders had higher incomes than similar responders gave similar results." The sentence reports the amount missing, the assumed mechanism, the method used, and how the conclusion was checked.

set.seed(410)
income <- round(rlnorm(10000, meanlog = log(55), sdlog = 0.6))  # simulated, $1000s
mean(income)                            # the true mean in this population

# MCAR: everyone has the same 30% chance of leaving the question blank
inc_mcar <- income
inc_mcar[runif(10000) < 0.30] <- NA
mean(inc_mcar, na.rm = TRUE)

# MNAR: the top quarter of earners leave it blank far more often
p_miss <- ifelse(income > quantile(income, 0.75), 0.70, 0.15)
inc_mnar <- income
inc_mnar[runif(10000) < p_miss] <- NA
mean(inc_mnar, na.rm = TRUE)

# Standard error of a mean = SD / square root of n (smaller = more precise)
sd(income) / sqrt(10000)
sd(inc_mcar, na.rm = TRUE) / sqrt(sum(!is.na(inc_mcar)))

How Much Missing Data Is Too Much?

No single threshold applies to every study. As a rough guide, many applied texts suggest that when fewer than about 5% of the rows needed for an analysis have missing values and MCAR is plausible, complete-case analysis is unlikely to change the conclusions much. When more is missing, when missingness is related to observed variables, or when MNAR is plausible, multiple imputation under a stated MAR assumption, together with a sensitivity analysis, is the safer choice. The proportion that matters is the share of rows lost from the specific analysis, which can be much larger than the share missing on any one variable, as the 687 complete rows above show.

This guide also explains where complete-case analysis fits in this lesson. It is defensible when little is missing and MCAR is plausible, which is the situation in the scale activity in the next section, where na.omit() drops the 15 people (about 2%) who skipped one depression item before a factor analysis. It becomes hard to defend when many rows are lost or when the people with missing values differ from the rest.

Missing-Data Codes Such as -9 and 999

Many data systems cannot leave a cell empty, so they store a code such as -9 (refused) or 999 (not asked) instead. R treats these codes as real numbers. A single age of 999 among the 800 ages in the course file would raise the mean by about 1.2 years, and in a small sample, as in the example below, one code can move a mean by decades. The codebook lists the codes for each variable, and the analyst converts them to NA before any summary is calculated. The course file uses empty cells only, so the example below uses a short made-up vector. The operator %in% asks, for each value on its left, whether it appears in the list on its right.

age_raw <- c(34, -9, 51, 999, 47)    # -9 = refused, 999 = not asked
mean(age_raw)                        # wrong: the codes are averaged as if they were ages
age_raw[age_raw %in% c(-9, 999)] <- NA
age_raw
mean(age_raw, na.rm = TRUE)          # na.rm = TRUE leaves the NAs out
# The codes can also be converted while reading a file:
# mydata <- read.csv("mydata.csv", na.strings = c("", "NA", "-9", "999"))

library(mice)
mi_vars <- phaa[, c("age", "gender", "smoker", "bmi", "phys_act_min",
                    "social_support_score")]
mi_vars$gender <- factor(mi_vars$gender)
mi_vars$smoker <- factor(mi_vars$smoker)
imp <- mice(mi_vars, m = 5, seed = 2026, printFlag = FALSE)  # 5 completed datasets
fit <- with(imp, lm(phys_act_min ~ 1))   # lm(y ~ 1) estimates the mean of y
summary(pool(fit), conf.int = TRUE)      # pooled with Rubin's rules

# Complete-case estimate and its standard error, for comparison
mean(phaa$phys_act_min, na.rm = TRUE)
sd(phaa$phys_act_min, na.rm = TRUE) / sqrt(sum(!is.na(phaa$phys_act_min)))

Data Transformations

When continuous variables are heavily skewed, transformations can make the distribution more symmetric, stabilize variance, and improve the validity of parametric methods.

Parametric methods, such as the t-test and linear regression, assume that the data (or the errors of a model) follow a particular distribution, usually the bell-shaped normal distribution, which is fully described by its mean and standard deviation. When that assumption is badly wrong, their standard errors and confidence intervals can be misleading. Rank-based (non-parametric) methods, described in a later section, make fewer assumptions.

The log of zero does not exist (R returns -Inf, minus infinity), so a variable that contains zeros, such as the number of clinic visits, cannot be logged directly. The usual remedy is to add 1 to every value first, ln(x + 1), which turns 0 into ln(1) = 0 and keeps all values in the same order. The constant 1 is a convention; when most values are small, the choice of constant affects the results, and the report should state it. R has a dedicated function, log1p(), for ln(x + 1), and expm1() reverses it.

The reciprocal, 1/x, reverses the order of values: the largest value becomes the smallest. Analysts often use −1/x instead, so that larger original values still have larger transformed values and the direction of any association stays easy to read.

crp <- c(0.4, 0.8, 1.1, 1.5, 2.3, 3.0, 4.8, 9.5, 21.0)  # a right-skewed biomarker (mg/L)
mean(crp)                  # arithmetic mean, pulled up by the 21.0
median(crp)
log_crp <- log(crp)        # natural log of each value
round(log_crp, 2)
mean(log_crp)              # the mean on the log scale
exp(mean(log_crp))         # back-transformed: the geometric mean, in mg/L

visits <- c(0, 0, 1, 2, 5, 12)
log(visits)                # log(0) is -Inf, which cannot be analysed
log1p(visits)              # log1p(x) = log(x + 1), and log1p(0) = 0
expm1(log1p(visits))       # expm1() undoes log1p()

1 / c(2, 5, 10)            # the reciprocal reverses the order
-1 / c(2, 5, 10)           # a minus sign restores it

Recoding and Derived Variables

Recoding involves creating new variables from existing ones. Common examples include collapsing categories (e.g., combining “current smoker” and “former smoker” into “ever smoker”), categorizing continuous variables into clinically meaningful groups (e.g., BMI categories), and computing derived variables such as person-years of follow-up or age at onset.

In R, ifelse(condition, value_if_true, value_if_false) creates a new variable from a condition, and cut() turns a continuous variable into bands. In cut(), the breaks argument gives the cut-points, right = FALSE places each cut-point in the band above it (so a BMI of exactly 25.0 falls in the band 25 to 29.9), and labels names the bands. A cross-tabulation of the old and new variables, as in the box below, is the standard check that a recode did what was intended.

String Cleaning and Standardization

Free-text fields often contain inconsistencies: “Vancouver,” “vancouver,” “VANCOUVER,” and “Vancoouver” may all refer to the same city. Standardization strategies include converting to consistent case, trimming whitespace, applying regular expressions for pattern matching, and using lookup tables or fuzzy matching algorithms.

A regular expression is a short pattern that describes text. The pattern ^vanc means "text that starts with vanc" (the symbol ^ marks the start), and [0-9]+ means "one or more digits". The function grepl() asks whether a pattern occurs in each value, and gsub() replaces it. A lookup table is a two-column list that maps every known variant to its standard spelling. Fuzzy matching finds strings that differ by only a letter or two, which helps when the variants cannot all be listed in advance; its suggested matches should be checked by a person. The function duplicated() returns TRUE for the second and later copies of a value, so sum(duplicated(phaa$id)) counts repeated IDs.

city <- c("Vancouver", "vancouver ", "VANCOUVER", "Vancoouver", " Burnaby")
city <- trimws(tolower(city))    # remove stray spaces, then make lower case
city
city[grepl("^vanc", city)] <- "vancouver"   # regular expression: starts with "vanc"
table(city)

# Recoding with ifelse() and cut()
smoker01 <- ifelse(phaa$smoker == "Yes", 1, 0)   # 1 = smoker, 0 = non-smoker
table(phaa$smoker, smoker01)
bmi_group <- cut(phaa$bmi, breaks = c(0, 18.5, 25, 30, Inf), right = FALSE,
                 labels = c("Under 18.5", "18.5 to 24.9", "25 to 29.9", "30 or more"))
table(bmi_group, useNA = "ifany")

# Duplicates: is any ID present more than once?
sum(duplicated(phaa$id))

Documenting the Cleaning Process

Every data cleaning decision should be recorded in an audit trail. This includes what was changed, why, by whom, and when. Reproducibility requires that another analyst could replicate the cleaning process from the original raw data using your documentation, an expectation echoed in the Statistical Methods in Psychology Journals guidelines (Wilkinson & the Task Force on Statistical Inference, 1999).

In practice the audit trail has two parts. The cleaning script records every change as code, with a comment giving the reason, and it is always run from the untouched raw file. A cleaning log, a small table with one row per changed value, records the date, the analyst, the ID, the variable, the old and new values and the reason, in a form that a reviewer can read without running any code. The cleaned data are then saved under a new file name so that the raw file is never overwritten; in this lesson that happens in the last line of the scale activity in the next section, which writes phaa_survey_clean.csv.

# ---------------------------------------------------------------
# clean_phaa.R   Cleaning script for phaa_survey.csv
# Raw file: phaa_survey.csv (never edited)
# Output:   phaa_survey_clean.csv (written at the end of the script)
# ---------------------------------------------------------------
cleaning_log <- data.frame(
  date     = "2026-09-30",
  analyst  = "AB",
  id       = c("P0033", "P0612", "P0017", "P0245", "P0078"),
  variable = c("age", "age", "systolic_bp", "systolic_bp", "bmi"),
  old      = c(-3, 220, 300, 60, 95),
  new      = NA,
  reason   = c("below 18: impossible", "above 100: impossible",
               "above 220: implausible", "below 80: implausible",
               "above 60: implausible")
)
cleaning_log
write.csv(cleaning_log, "cleaning_log.csv", row.names = FALSE)

● Complete the reflection and knowledge check to continue.

Introduction and Overview

Earlier sections produced a clean, well-documented dataset. This section turns that dataset into the descriptive summaries that any analysis report opens with: means and medians, spreads and shapes, and frequency distributions, “Table 1” characteristics by exposure group, and the visualisations that go with them. Doing this carefully is what lets readers understand what your data look like before they trust your inferential models.

Measures of Central Tendency

Measures of central tendency summarize the “typical” value in a distribution. The choice depends on the variable type and distributional shape.

Measures of Spread

Spread (or dispersion) describes how much variability exists in the data around the central value.

The variance divides by n − 1 instead of n. The deviations are measured from the sample mean, which sits in the middle of the sample by construction, so they are slightly smaller on average than deviations from the true population mean would be. Dividing by n − 1, a slightly smaller number, corrects for this. The standard deviation is the square root of the variance, which brings the value back from squared units (squared mmHg, for instance) to the original units so that it can be reported beside the mean.

income5 <- c(30, 35, 40, 45, 50)        # five incomes, $1000s per year
mean(income5)
median(income5)
income6 <- c(income5, 500)              # add one very high earner
mean(income6)
median(income6)

# Standard deviation, step by step
sbp <- c(110, 116, 120, 124, 130)
dev <- sbp - mean(sbp)                  # distance of each value from the mean
dev
dev^2                                   # squared distances
sum(dev^2) / (length(sbp) - 1)          # the variance
sqrt(sum(dev^2) / (length(sbp) - 1))    # the standard deviation
sd(sbp)                                 # R's built-in function gives the same

Measures of Shape

Beyond central tendency and spread, the shape of a distribution provides critical information for choosing analytical methods.

Reading summary() and describe() Output

Two functions produce most numerical summaries in this course. Base R's summary() gives the minimum, first quartile, median, mean, third quartile and maximum, plus a count of missing values when there are any. The describe() function from the psych package gives a longer row for each variable, including the standard deviation, skew and excess kurtosis. The box below applies both to three variables that the activity below does not use.

summary(phaa$diastolic_bp)
library(psych)
describe(phaa[, c("diastolic_bp", "discrimination_score", "social_support_score")])

Frequency Distributions and Cross-Tabulations

For categorical variables, frequency distributions display the count and percentage of observations in each category. Cross-tabulations (contingency tables) display the joint distribution of two categorical variables and are the basis for computing measures of association such as odds ratios and risk ratios. Some readers will have calculated odds ratios in an earlier course; the worked example after the table builds them from the start.

# The outbreak table from the text
outbreak <- matrix(c(45, 15, 30, 60), nrow = 2,
                   dimnames = list(Exposure = c("Exposed", "Unexposed"),
                                   Group    = c("Case", "Control")))
outbreak
odds_cases    <- outbreak["Exposed", "Case"] / outbreak["Unexposed", "Case"]
odds_controls <- outbreak["Exposed", "Control"] / outbreak["Unexposed", "Control"]
odds_cases
odds_controls
odds_cases / odds_controls              # the odds ratio

# The same steps with the course data: hypertension by smoking status
tab <- table(Smoker = phaa$smoker, Hypertension = phaa$hypertension)
tab
odds_smokers     <- tab["Yes", "Yes"] / tab["Yes", "No"]
odds_non_smokers <- tab["No", "Yes"] / tab["No", "No"]
round(c(smokers = odds_smokers, non_smokers = odds_non_smokers,
        OR = odds_smokers / odds_non_smokers), 3)

Visualization Techniques

Appropriate visualization depends on the variable type and the question being asked. Exploratory data analysis emphasises looking at data first so structure, gaps, and anomalies are visible before any model is fit.

Histograms Boxplots Bar Charts Scatter & Line

Histograms display the distribution of a continuous variable by dividing values into bins and counting observations per bin. They reveal distributional shape, modality, gaps, and potential outliers. The choice of bin width affects interpretation: too few bins obscure patterns, too many create noise.

The figure below shows the same 850 simulated readings, drawn from two groups with different average blood pressures, using 5, 20 and 100 bins. In R the number of bins is suggested with the breaks argument, for example hist(phaa$systolic_bp, breaks = 20); R treats the number as a suggestion and chooses round cut-points near it.

Boxplots display the median, the IQR (the box), whiskers that extend to the most extreme observation lying within 1.5×IQR of the box, and individual points beyond the whiskers. They are ideal for comparing distributions across groups (e.g., blood pressure by treatment arm) and for quickly identifying asymmetry and outliers.

Bar charts display frequencies or proportions for categorical variables. Use them for comparing counts across groups. Avoid using bar charts for continuous data (use histograms instead) and avoid 3D bar charts, which distort proportions.

Scatter plots display the relationship between two continuous variables and can reveal linear or nonlinear associations, clusters, and outliers. Line graphs are appropriate for displaying trends over time (e.g., incidence rates by year). Both are essential in epidemiologic exploratory analysis.

# Categorical descriptives ---------------------------------------------------
table(phaa$gender)                            # frequencies
prop.table(table(phaa$gender))                # proportions
round(prop.table(table(phaa$gender)) * 100, 1)  # %

# Cross-tab: gender by smoker, row percentages
round(prop.table(table(phaa$gender, phaa$smoker), margin = 1) * 100, 1)

# Numeric descriptives --------------------------------------------------------
summary(phaa$age)
sd(phaa$age, na.rm = TRUE)
IQR(phaa$age, na.rm = TRUE)

# psych::describe() summarises many variables at once
library(psych)
describe(phaa[, c("age", "bmi", "systolic_bp", "phys_act_min")])

# Plots: histogram, boxplot stratified by exposure, bar chart, scatter plot ---
hist(phaa$systolic_bp, main = "Systolic BP", xlab = "mmHg")

boxplot(systolic_bp ~ smoker, data = phaa,
        main = "Systolic BP by smoking status", ylab = "mmHg")

barplot(table(phaa$gender), main = "Gender distribution")

plot(phaa$age, phaa$systolic_bp, main = "Systolic BP by age",
     xlab = "Age (years)", ylab = "Systolic BP (mmHg)")

Descriptive Statistics for Epidemiologic Data: “Table 1”

In epidemiologic publications, “Table 1” typically presents the baseline characteristics of the study population, often stratified by exposure or outcome status. How baseline characteristics and group comparisons are reported has been the subject of methodological scrutiny. Pocock, Assmann, Enos, & Kasten (2002) reviewed published trial reports and catalogued common problems with baseline comparisons, subgroup analyses and covariate adjustment. A well-constructed Table 1 includes:

• Demographics (age, sex, ethnicity) with appropriate summary statistics
• Key clinical or exposure variables
• Missing data counts for each variable
• A comparison column between groups, preferably standardised differences (explained below)

# Three small helper functions. function(v) { ... } defines a reusable recipe
# that is applied to the whole sample and then to each smoking group.
grp <- phaa$smoker                       # the grouping variable ("No", "Yes")

mean_sd <- function(v) {                 # "mean (SD)" and the standardised difference
  f <- function(x) sprintf("%.1f (%.1f)", mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
  m <- tapply(v, grp, mean, na.rm = TRUE)
  s <- tapply(v, grp, sd, na.rm = TRUE)
  d <- (m["Yes"] - m["No"]) / sqrt((s["Yes"]^2 + s["No"]^2) / 2)
  c(Overall = f(v), tapply(v, grp, f), SMD = sprintf("%.2f", d))
}
n_pct <- function(v, level) {            # "n (%)" and the standardised difference
  f <- function(x) sprintf("%d (%.1f%%)", sum(x == level, na.rm = TRUE),
                           100 * mean(x == level, na.rm = TRUE))
  p <- tapply(v == level, grp, mean, na.rm = TRUE)
  d <- (p["Yes"] - p["No"]) / sqrt((p["Yes"] * (1 - p["Yes"]) + p["No"] * (1 - p["No"])) / 2)
  c(Overall = f(v), tapply(v, grp, f), SMD = sprintf("%.2f", d))
}
n_miss <- function(v) {                  # number of missing values
  c(Overall = sum(is.na(v)), tapply(is.na(v), grp, sum), SMD = "")
}

table1 <- rbind(
  "N"                     = c(Overall = length(grp), table(grp), SMD = ""),
  "Age, mean (SD)"        = mean_sd(phaa$age),
  "  Age missing, n"      = n_miss(phaa$age),
  "Woman, n (%)"          = n_pct(phaa$gender, "Woman"),
  "Degree, n (%)"         = n_pct(phaa$education >= "Bachelor's", TRUE),
  "BMI, mean (SD)"        = mean_sd(phaa$bmi),
  "SBP, mean (SD)"        = mean_sd(phaa$systolic_bp),
  "Activity, mean (SD)"   = mean_sd(phaa$phys_act_min),
  "  Activity missing, n" = n_miss(phaa$phys_act_min),
  "Income missing, n (%)" = n_pct(is.na(phaa$income), TRUE)
)
noquote(table1)

Stratified Descriptive Analyses

Stratifying descriptive statistics by key exposure or covariate groups allows you to assess the distribution of potential confounders across exposure categories. This is a critical step before proceeding to multivariable modelling, as it helps identify imbalances that may require adjustment.

A confounder is a third variable that is related to both the exposure and the outcome and can therefore make the two appear more, or less, strongly associated than they truly are. Age is a common example, because it is linked to many exposures and to blood pressure. A large standardised difference in Table 1 marks a variable that differs between the exposure groups; whether it also affects the outcome is judged from subject knowledge and the later regression models. A variable that is itself affected by the exposure is a mediator and does not meet the definition of a confounder.

Building Scales: Internal Consistency and Factor Analysis

Many public-health surveys ask several items that together measure a latent construct (e.g., depression, anxiety, social support). Two questions need to be answered before you can use these items as a single scale variable in a regression: (1) internal consistency (do the items hang together?), assessed with Cronbach’s α; and (2) dimensionality (how many underlying factors do the items measure?), assessed with exploratory factor analysis (EFA). Once both checks pass, the items can be summed (or averaged) into a derived scale variable that you carry forward into later lessons.

Each depression item asks how often a symptom occurred, on a five-point response scale from 1 (Never) to 5 (Always), so a higher number always means more symptoms. Some questionnaires also include reverse-coded items, worded in the opposite direction (for example, "I felt hopeful about the future"), where a high answer means fewer symptoms. Such items must be flipped before alpha is calculated or a score is built; on a 1 to 5 scale the flipped value is 6 minus the original, so 5 becomes 1 and 1 becomes 5. The alpha() function prints a warning when an item correlates negatively with the rest of the scale, which is the usual sign of an item that has not been reversed. None of the course items is reverse-coded.

Cronbach's alpha can be understood as a measure of how well the items agree. If people who score high on one item tend to score high on the others, the items share a common source, the construct, and alpha is high. Alpha rises when the items are more strongly correlated with each other and when there are more of them. Under the assumptions behind it, an alpha of 0.90 indicates that about 90% of the variation in total scores reflects what the items have in common and about 10% reflects item-specific noise. Values of about 0.70 or higher are conventionally acceptable and 0.80 or higher good; values above about 0.95 can mean that some items are near-duplicates.

Two terms from factor analysis appear in the output. An eigenvalue measures how much of the shared variation among the items one factor accounts for; a large first eigenvalue followed by small ones suggests a single underlying factor. Parallel analysis compares each eigenvalue from the real data with the eigenvalue that random data of the same size would produce, and keeps factors, counting from the first, for as long as the real eigenvalue is larger. A factor loading is the correlation between an item and a factor; loadings of 0.40 or higher are conventionally treated as meaningful, and this lesson uses that threshold throughout.

R Activity: Cronbach’s α and exploratory factor analysis

Files for this activity phaa_survey.csvdataset Answer key (R script)revealed after you save your responses

The course dataset contains seven depression items (dep1 … dep7) and five anxiety items (anx1 … anx5). Below we (1) confirm internal consistency, (2) decide on the number of factors, and (3) build the derived scale scores we use for the rest of the course.

# 1. Internal consistency for each candidate scale ----------------------------
library(psych)
dep_items <- phaa[, c("dep1","dep2","dep3","dep4","dep5","dep6","dep7")]
anx_items <- phaa[, c("anx1","anx2","anx3","anx4","anx5")]

alpha(dep_items)             # raw_alpha >= 0.70 is acceptable, >= 0.80 is good
alpha(anx_items)
dep_alpha <- alpha(dep_items)            # store the result to read it at 3 decimals
round(dep_alpha$total$raw_alpha, 3)
round(dep_alpha$alpha.drop[, "raw_alpha", drop = FALSE], 3)  # alpha if an item is dropped

# 2. How many underlying factors? ------------------------------------------
fa.parallel(dep_items, fa = "fa")   # parallel analysis for factors

# 3. Inspect a one-factor solution for the depression items ----------------
factanal(x = na.omit(dep_items), factors = 1)
# na.omit(): factanal() needs complete rows and dep4 has 15 NAs, so these 15
# people are left out of this check (a complete-case analysis).
# With a single factor there is nothing to rotate, so no rotation is requested.
# Loadings >= 0.40 contribute meaningfully to the factor.

# 4. Two-factor solution combining dep + anx items -------------------------
combined <- na.omit(cbind(dep_items, anx_items))
factanal(x = combined, factors = 2, rotation = "varimax")
# dep1-7 should load on one factor; anx1-5 on the other.

# 5. Build derived scale variables for use in later lessons ------------------
# Score = mean of the answered items x number of items, calculated only when
# at most one item is missing (otherwise NA). A missing item is never
# counted as 0, which would understate the person's score.
n_dep <- rowSums(!is.na(dep_items))     # items answered by each person
n_anx <- rowSums(!is.na(anx_items))
phaa$dep_score <- ifelse(n_dep >= 6, rowMeans(dep_items, na.rm = TRUE) * 7, NA)
phaa$anx_score <- ifelse(n_anx >= 4, rowMeans(anx_items, na.rm = TRUE) * 5, NA)
summary(phaa$dep_score)

# 6. Save the cleaned, scale-augmented file under a new name -----------------
write.csv(phaa, "phaa_survey_clean.csv", row.names = FALSE)

Console output (depression items; the anxiety and two-factor output are left for the questions)

> alpha(dep_items) Reliability analysis Call: alpha(x = dep_items) raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r 0.9 0.9 0.89 0.57 9.4 0.0052 3 1.1 0.58 95% confidence boundaries lower alpha upper Feldt 0.89 0.9 0.91 Duhachek 0.89 0.9 0.91 Reliability if an item is dropped: raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r dep1 0.89 0.89 0.87 0.57 8.0 0.0061 0.00111 0.57 dep2 0.89 0.89 0.87 0.56 7.8 0.0062 0.00118 0.57 dep3 0.89 0.89 0.87 0.57 7.9 0.0062 0.00157 0.57 dep4 0.89 0.89 0.87 0.57 7.8 0.0062 0.00132 0.57 dep5 0.89 0.89 0.87 0.58 8.3 0.0059 0.00126 0.59 dep6 0.89 0.89 0.87 0.57 8.0 0.0061 0.00129 0.57 dep7 0.90 0.90 0.88 0.59 8.6 0.0057 0.00042 0.59 Item statistics n raw.r std.r r.cor r.drop mean sd dep1 800 0.80 0.80 0.76 0.72 3 1.4 dep2 800 0.82 0.82 0.78 0.74 3 1.4 dep3 800 0.81 0.81 0.77 0.73 3 1.4 dep4 785 0.81 0.81 0.78 0.74 3 1.4 dep5 800 0.78 0.78 0.73 0.69 3 1.4 dep6 800 0.80 0.80 0.75 0.72 3 1.4 dep7 800 0.75 0.75 0.69 0.65 3 1.4 ... (10 more lines not shown) > round(dep_alpha$total$raw_alpha, 3) [1] 0.904 > round(dep_alpha$alpha.drop[, "raw_alpha", drop = FALSE], 3) raw_alpha dep1 0.889 dep2 0.886 dep3 0.888 dep4 0.887 dep5 0.892 dep6 0.889 dep7 0.896 > fa.parallel(dep_items, fa = "fa") Parallel analysis suggests that the number of factors = 1 and the number of components = NA > factanal(x = na.omit(dep_items), factors = 1) Call: factanal(x = na.omit(dep_items), factors = 1) Uniquenesses: dep1 dep2 dep3 dep4 dep5 dep6 dep7 0.405 0.368 0.401 0.391 0.465 0.417 0.523 Loadings: Factor1 dep1 0.772 dep2 0.795 dep3 0.774 dep4 0.781 dep5 0.731 dep6 0.763 dep7 0.691 Factor1 SS loadings 4.030 Proportion Var 0.576 Test of the hypothesis that 1 factor is sufficient. The chi square statistic is 12.7 on 14 degrees of freedom. The p-value is 0.551

Output	What it shows and what to report
raw_alpha	Cronbach's alpha, printed to two decimals as 0.9. The line round(dep_alpha$total$raw_alpha, 3) gives 0.904, which is the value to report (α = 0.904).
std.alpha, average_r	Alpha calculated from standardised items (similar here because all items share the same response scale), and the average correlation between pairs of items (0.57).
95% confidence boundaries	The range of alpha values compatible with the data, 0.89 to 0.91.
Reliability if an item is dropped	Each row gives alpha with that item removed. A value above the overall 0.904 would mean the scale is better without the item. Every value here is lower (0.886 to 0.896 at three decimals), so every item contributes; dep7 contributes least, because removing it would lower alpha only to 0.896.
r.drop (item statistics)	The correlation of each item with the total of the other items, from 0.65 (dep7) to 0.74. Values below about 0.30 would mark an item that does not fit the scale.
fa.parallel()	The message reports one factor (the number of components is NA because fa = "fa" requests factors only). The annotated plot below shows why.
factanal() Loadings	All seven loadings lie between 0.691 (dep7) and 0.795 (dep2), well above 0.40. Uniquenesses are the share of each item that the factor does not explain (1 minus the squared loading).
Proportion Var and the chi-square test	The single factor explains 57.6% of the item variance. The test asks whether one factor is enough; its large p-value (0.551) gives no evidence that more factors are needed.

How to read this. An α of 0.904 indicates that the seven depression items measure essentially the same construct. The single factor explains 57.6% of the item variance and every loading is well above the 0.40 threshold (the lowest is 0.691, for dep7), so combining the items into dep_score is justified. Step 5 builds the score as the mean of the answered items multiplied by the number of items (a prorated sum), and only for people who answered at least six of the seven depression items (four of the five anxiety items). Using rowSums(dep_items, na.rm = TRUE) instead would count a skipped item as 0 and understate the scores of the 15 people who skipped dep4, by between 1.8 and 4.5 points (2.8 on average). The score is used as a continuous predictor in a later lesson and as a covariate in the logistic model for hypertension.

R Reflect on what you just ran

Use the questions below to interpret the output you produced. Look at your console / plot before answering.

1. What is the raw alpha for the depression items? For the anxiety items? Which scale (if either) clears the conventional 0.80 threshold for "good" internal consistency, and which item (look at the "Reliability if an item is dropped" table) contributes least to the depression scale?

Model answerRaw α for the seven depression items is 0.904 (good internal consistency), matching the console output above; the anxiety scale is a little lower at 0.855. Both scales clear the 0.80 threshold. In the "reliability if an item is dropped" table, the item that contributes least is the one whose alpha-if-dropped value is highest, equivalently the one with the lowest item-total correlation: dep7. Removing it would lower alpha the least (from 0.904 to 0.896), and its r.drop of 0.65 is the lowest (0.69 to 0.74 for the other six). Because every alpha-if-dropped value is below 0.904, every item, dep7 included, adds to the scale.

2. From fa.parallel(dep_items), how many factors have eigenvalues above the parallel-analysis line? Does the one-factor solution from factanal() show all seven loadings >= 0.40?

Model answerfa.parallel() shows exactly one factor above the parallel-analysis line for the depression items (first eigenvalue 4.02 against about 0.5 for random data, a value that varies slightly from run to run; the second is 0.05) and prints "the number of factors = 1", supporting a unidimensional structure. The one-factor solution from factanal() has all seven loadings ≥ 0.40 (they run from 0.691 for dep7 to 0.795 for dep2, and the factor explains 57.6% of the item variance), confirming that all items track the same latent construct.

3. In the two-factor combined solution, do the depression items load on a different factor than the anxiety items? Cite one specific loading from your output that supports your answer.

Model answerIn the two-factor combined solution, the depression and anxiety items load on different factors. A specific example: depression item dep3 loads 0.748 on Factor 1 and 0.205 on Factor 2, while anxiety item anx2 loads 0.124 on Factor 1 and 0.714 on Factor 2. The cross-loadings (no item loads above 0.21 on the other factor) and the loadings on each item's own factor (all at least 0.68, well above the 0.40 threshold) together support the interpretation that the two scales measure distinct, though related, latent constructs.

Save responses

Saved.

Normality Assessment

Assessing whether a variable follows a normal distribution informs the choice between parametric and non-parametric methods.

qqnorm(phaa$systolic_bp, main = "Systolic BP: normal Q-Q plot")
qqline(phaa$systolic_bp)          # the reference line
shapiro.test(phaa$systolic_bp)
shapiro.test(phaa$diastolic_bp)

From the Checks to a Decision

In practice the analyst combines the plot, the skew value and the purpose of the analysis. If the Q-Q plot is close to a straight line and |skew| is below about 1, the variable is summarised as mean (SD) and parametric methods are reasonable, even when a large sample makes the Shapiro-Wilk p-value small, as for systolic BP here. If the plot curves clearly or |skew| is above about 1, the variable is summarised as median (IQR), and the analyst either transforms it (as in an earlier section) or uses a rank-based method. Rank-based methods replace each value by its rank (1 for the smallest, 2 for the next, and so on), so a few extreme values cannot dominate the result. The table below pairs the common options.

● Complete the reflection and knowledge check to continue.

Bringing It All Together

This lesson followed a single dataset from raw to descriptive. The section on data quality named the errors (collection, entry, storage, and analysis-stage problems) and unpacked the missing-data mechanisms (MCAR, MAR, MNAR) that determine which cleaning approaches are valid. The section on cleaning strategies turned that vocabulary into a workflow: systematic range, type, and consistency checks; outlier triage; principled missing-data handling; transformations and derived variables; and a documentation trail that a peer reviewer could follow. The section on descriptive analyses closed the loop with the layer that every paper opens with: central tendency, spread, shape, frequency tables, visualisations, scale checks, and a Table 1 that stratifies on the primary exposure or outcome.

Descriptive analysis is a core stage of any analysis. It is the layer where data problems surface, where assumptions get tested, and where readers calibrate how much trust to place in the inferential models that follow. From here we can move into the regression machinery of the rest of the course knowing that the dataset under it has been examined and understood.

Final Knowledge Assessment

You must answer all 15 questions correctly (100%) and complete the final reflection to finish this lesson.