Base R data frame manipulation

Lesson 7 — R Programming

Lesson 7 of 13 Intermediate ~150 min

Learning objectives

Explain when a dependency-free solution is the right choice, not a fallback
Subset, sort, deduplicate and modify data frames with base R only
Group and summarise with aggregate(), tapply(), ave() and split()
Use merge() safely, and assert the cardinality it does not check
Reshape between long and wide with reshape()
Produce frequency tables with correct zero counts and denominators
Build a complete demographics summary without a single package
Translate confidently between base R and dplyr

Why base R only

Lessons 5 and 6 taught dplyr and tidyr, and for most analytical work they are the right tools. This lesson does the same work with no packages at all.

That is not an academic exercise. In pharmaceutical programming there are four situations where it is the correct choice:

1. A restricted environment. Many organisations maintain an approved package list. Getting dplyr and its dependency tree assessed, versioned and qualified is real work, and on some systems it has not been done. Base R ships with R itself and is covered by the R Foundation’s regulatory compliance documentation.

2. A package you are writing. Every entry in Imports is a dependency your users inherit and your validation covers. A standards package with Imports: stats, utils is markedly easier to qualify than one pulling in forty transitive dependencies. See Package development.

3. Legacy code. Study code written in 2015 is base R. You cannot maintain what you cannot read, and much of it will outlive the toolchain you prefer.

4. Understanding what dplyr does. dplyr is a well-designed layer over operations that base R performs explicitly. Seeing merge() multiply rows teaches you something that left_join() hides.

This is not an argument against the tidyverse

dplyr is clearer, safer in several specific ways documented below, and the right default when it is available. The argument is that “I cannot use dplyr here” should never mean “I cannot do this here”.

Throughout, examples use a small ADSL-shaped data frame:

adsl <- data.frame(
  USUBJID = c("001", "002", "003", "004", "005", "006"),
  SITEID  = c("001", "001", "002", "002", "003", "003"),
  ARM     = c("Placebo", "Drug A", "Placebo", "Drug A", "Drug A", "Placebo"),
  AGE     = c(45, 72, 38, 66, NA, 81),
  SEX     = c("F", "M", "F", "M", "F", "F"),
  WEIGHT  = c(62.5, 88.1, 55.0, 91.3, 68.9, 59.4),
  SAFFL   = c("Y", "Y", "N", "Y", "Y", "Y"),
  stringsAsFactors = FALSE
)

Subsetting rows

The `[` operator

adsl[adsl$AGE > 60, ]
adsl[adsl$SEX == "F" & adsl$SAFFL == "Y", ]
adsl[adsl$ARM %in% c("Placebo", "Drug A"), ]
adsl[1:3, ]
adsl[-1, ]

[ with NA in the condition creates phantom rows

Subject 005 has a missing age.

adsl[adsl$AGE > 60, ]
#>      USUBJID SITEID     ARM AGE  SEX WEIGHT SAFFL
#> 2        002    001  Drug A  72    M   88.1     Y
#> 4        004    002  Drug A  66    M   91.3     Y
#> NA      <NA>   <NA>    <NA>  NA <NA>     NA  <NA>   <- a row that does not exist
#> 6        006    003 Placebo  81    F   59.4     Y

adsl$AGE > 60 evaluates to TRUE, TRUE, FALSE, TRUE, NA, TRUE. Base [ treats the NA index as “a row I know nothing about” and materialises a row of NAs.

That phantom row then inflates nrow(), appears in counts, and propagates into joins. It is the single most dangerous difference between base subsetting and dplyr::filter().

Three fixes:

adsl[which(adsl$AGE > 60), ]      # which() discards NA — the idiomatic fix
subset(adsl, AGE > 60)            # subset() handles NA correctly
adsl[!is.na(adsl$AGE) & adsl$AGE > 60, ]   # explicit

Use which(). It is short, obvious once you know it, and correct.

`subset()`

subset(adsl, AGE > 60)
subset(adsl, SEX == "F" & SAFFL == "Y")
subset(adsl, AGE > 60, select = c(USUBJID, AGE, ARM))
subset(adsl, select = -c(WEIGHT, SITEID))
subset(adsl, AGE > 60, select = USUBJID:AGE)

subset() is convenient and handles NA correctly. One caution:

subset() is for interactive use

subset() uses non-standard evaluation — AGE is resolved inside the data frame. R’s own documentation warns against using it inside functions, because the column name is resolved at evaluation time in a way that is hard to control programmatically.

Inside a function, use [ with which():

filter_rows <- function(data, column, threshold) {
  data[which(data[[column]] > threshold), , drop = FALSE]
}

Note drop = FALSE — a single-column data frame would otherwise collapse to a vector.

Selecting columns

adsl[, c("USUBJID", "AGE", "ARM")]
adsl[c("USUBJID", "AGE")]                   # list-style, same result
adsl[, !names(adsl) %in% c("WEIGHT", "SITEID")]
adsl[, grepl("^S", names(adsl))]            # names starting with S
adsl[, sapply(adsl, is.numeric)]            # numeric columns only

adsl$WEIGHT <- NULL                         # delete a column
adsl[["SITEID"]] <- NULL

Reordering and renaming:

adsl <- adsl[, c("USUBJID", "ARM", "AGE", "SEX", "SAFFL")]

names(adsl)[names(adsl) == "ARM"] <- "TRT01P"      # rename one
names(adsl) <- toupper(names(adsl))                # rename all

# Move a column to the front
front <- "USUBJID"
adsl <- adsl[, c(front, setdiff(names(adsl), front))]

df[, "col"] returns a vector, df["col"] returns a data frame

class(adsl[, "AGE"])     #> [1] "numeric"
class(adsl["AGE"])       #> [1] "data.frame"
class(adsl[, c("AGE", "SEX")])   #> [1] "data.frame"

The return type of [, ] depends on how many columns you selected. In a function taking a column selection as an argument, that is a landmine: it works in testing with two columns and breaks in production with one.

adsl[, "AGE", drop = FALSE]      # always a data frame

Always pass drop = FALSE in a function. This is the inconsistency tibbles were created to remove.

Adding and modifying columns

adsl$AGE_MONTHS <- adsl$AGE * 12
adsl[["BMI"]] <- adsl$WEIGHT / (adsl$HEIGHT / 100)^2

# Conditional, two branches
adsl$AGEGR <- ifelse(adsl$AGE >= 65, ">=65", "<65")

# Conditional, several branches — nested ifelse is unreadable; use cut()
adsl$AGEGR1 <- cut(
  adsl$AGE,
  breaks = c(-Inf, 18, 65, 80, Inf),
  labels = c("<18", "18-64", "65-79", ">=80"),
  right  = FALSE                 # [lower, upper) — left-closed
)

cut() is the base-R answer to case_when() for numeric banding, and right = FALSE is almost always what a clinical specification means by “65–79”.

For non-numeric multi-branch logic, a lookup vector is cleaner than nested ifelse():

arm_labels <- c(PBO = "Placebo", A50 = "Drug A 50mg", A100 = "Drug A 100mg")
adsl$TRT01P <- unname(arm_labels[adsl$ARMCD])

arm_n <- c("Placebo" = 0, "Drug A 50mg" = 1, "Drug A 100mg" = 2)
adsl$TRT01PN <- unname(arm_n[adsl$TRT01P])

This is the named-vector-as-lookup pattern from lesson 3, and it replaces an entire case_when() block with one expression.

`transform()` and `within()`

adsl <- transform(adsl,
  AGE_MONTHS = AGE * 12,
  AGEGR      = ifelse(AGE >= 65, ">=65", "<65")
)

transform() cannot reference a column it creates in the same call. within() can:

adsl <- within(adsl, {
  HEIGHT_M <- HEIGHT / 100
  BMI      <- WEIGHT / HEIGHT_M^2        # uses HEIGHT_M from the line above
  BMI_CAT  <- cut(BMI, c(-Inf, 18.5, 25, 30, Inf),
                  labels = c("Under", "Normal", "Over", "Obese"))
})

within() adds columns in reverse order

names(within(data.frame(x = 1), { a <- 1; b <- 2 }))
#> [1] "x" "b" "a"

The block is an environment, and the order columns come back is not the order you wrote them. If column order matters — and in a submission dataset it does — reorder explicitly afterwards, or use transform() and sequential assignments.

Both functions also use non-standard evaluation, so the same caution applies as for subset(): fine in a script, awkward inside a function.

Sorting

adsl[order(adsl$AGE), ]
adsl[order(-adsl$AGE), ]                      # descending, numeric only
adsl[order(adsl$AGE, decreasing = TRUE), ]
adsl[order(adsl$ARM, adsl$AGE), ]             # multiple keys
adsl[order(adsl$ARM, -adsl$AGE), ]            # ARM ascending, AGE descending

For a descending character sort, negation does not work — use xtfrm() or rev() of the levels:

adsl[order(adsl$ARM, decreasing = TRUE), ]
adsl[order(-xtfrm(adsl$ARM), adsl$AGE), ]     # mixed directions on character

NA handling:

adsl[order(adsl$AGE), ]                       # NA last (default)
adsl[order(adsl$AGE, na.last = FALSE), ]      # NA first — SAS behaviour
adsl[order(adsl$AGE, na.last = NA), ]         # NA rows removed

SAS sorts missing first; R sorts it last

PROC SORT places missing values first because missing is smaller than any number. R’s order() defaults to na.last = TRUE.

In an independent double-programming comparison this produces different “first record per subject” selections wherever the sort key is missing — a difference that looks like a derivation error and is a convention difference. See Validation and testing.

Always reset row names after sorting, or the original positions persist and confuse later indexing:

adsl <- adsl[order(adsl$ARM, adsl$AGE), ]
rownames(adsl) <- NULL

Deduplicating and first/last per subject

A core clinical operation: one record per subject, or the earliest event per subject.

unique(adsl$ARM)                      #> distinct values of one column
unique(adsl[, c("ARM", "SEX")])       #> distinct combinations
duplicated(adsl$USUBJID)              #> logical, TRUE for repeats
sum(duplicated(adsl$USUBJID))         #> how many duplicates
anyDuplicated(adsl$USUBJID)           #> index of the first, or 0

First record per subject, equivalent to SAS if first.usubjid:

adae <- adae[order(adae$USUBJID, adae$ASTDT, adae$AESEQ), ]
first_ae <- adae[!duplicated(adae$USUBJID), ]

Last record per subject, SAS if last.usubjid:

last_ae <- adae[!duplicated(adae$USUBJID, fromLast = TRUE), ]

fromLast = TRUE is the whole trick and is widely under-known.

First per subject and parameter:

key <- paste(adlb$USUBJID, adlb$PARAMCD, sep = "|")
adlb <- adlb[order(adlb$USUBJID, adlb$PARAMCD, adlb$ADT), ]
baseline <- adlb[!duplicated(key), ]

Building a composite key with paste() and a separator that cannot appear in the data is the standard base-R approach to multi-column grouping.

The sort must be deterministic

duplicated() keeps whichever row happens to come first. If two records tie on every sort key, the one kept is arbitrary and may differ between runs, platforms or data extracts.

Always sort on enough keys to break every tie — typically down to --SEQ:

adae <- adae[order(adae$USUBJID, adae$ASTDT, adae$AETERM, adae$AESEQ), ]

This is the same requirement as deterministic --SEQ derivation in SDTM programming.

Grouping and summarising

Base R offers several tools. Which to use depends on the shape you want back.

`aggregate()`

The closest thing to group_by() |> summarise().

aggregate(AGE ~ ARM, data = adsl, FUN = mean)
#>       ARM      AGE
#> 1  Drug A 69.00000
#> 2 Placebo 54.66667

aggregate(AGE ~ ARM + SEX, data = adsl, FUN = mean)       # two grouping vars
aggregate(cbind(AGE, WEIGHT) ~ ARM, data = adsl, FUN = mean)   # two measures

Several statistics at once, returning a matrix column:

res <- aggregate(AGE ~ ARM, data = adsl, FUN = function(x) {
  c(n = length(x), mean = mean(x), sd = sd(x))
})

# The result needs flattening — a known awkwardness
res <- do.call(data.frame, res)
names(res) <- c("ARM", "n", "mean", "sd")
res
#>       ARM n     mean       sd
#> 1  Drug A 3 69.00000 13.07670
#> 2 Placebo 3 54.66667 19.13985

The formula interface silently drops missing values

nrow(adsl)                                          #> [1] 6
sum(aggregate(AGE ~ ARM, adsl, length)$AGE)         #> [1] 5

Subject 005 has a missing AGE and the formula interface applies na.action = na.omit by default — so the row vanishes and the counts no longer match the population.

For a clinical table this is wrong: the denominator must be the population, not the non-missing count. Two fixes:

aggregate(AGE ~ ARM, data = adsl, FUN = mean, na.action = na.pass)
aggregate(adsl["AGE"], by = list(ARM = adsl$ARM), FUN = mean, na.rm = TRUE)

The second form — by = with a named list — avoids the formula interface entirely and is safer for exactly this reason.

`tapply()`

One statistic, one grouping variable, returned as a named vector:

tapply(adsl$AGE, adsl$ARM, mean, na.rm = TRUE)
#>   Drug A  Placebo
#> 69.00000 54.66667

tapply(adsl$AGE, list(adsl$ARM, adsl$SEX), mean, na.rm = TRUE)   # a matrix

Convenient for a single number per group. Awkward when you want a data frame back.

`ave()` — group statistics broadcast to every row

The base-R equivalent of mutate(..., .by = group), and the most useful of these functions for clinical derivation.

adsl$ARM_MEAN_AGE <- ave(adsl$AGE, adsl$ARM, FUN = function(x) mean(x, na.rm = TRUE))

adsl[, c("USUBJID", "ARM", "AGE", "ARM_MEAN_AGE")]
#>   USUBJID     ARM AGE ARM_MEAN_AGE
#> 1     001 Placebo  45     54.66667
#> 2     002  Drug A  72     69.00000

The result is the same length as the input, so it can be assigned straight back as a column. This is how you carry a baseline value to every visit in base R:

adlb <- adlb[order(adlb$USUBJID, adlb$PARAMCD, adlb$AVISITN), ]

adlb$BASE <- ave(
  adlb$AVAL,
  adlb$USUBJID, adlb$PARAMCD,
  FUN = function(x) x[1]                 # first value in each group
)
adlb$CHG <- adlb$AVAL - adlb$BASE

ave() also gives you group sequence numbers and group sizes:

adae$AESEQ <- ave(seq_len(nrow(adae)), adae$USUBJID,
                  FUN = seq_along)                        # 1, 2, 3... per subject
adae$N_EVENTS <- ave(adae$USUBJID, adae$USUBJID, FUN = length)

`split()` plus `lapply()`

The most flexible approach: break the data frame into a list, operate on each piece, recombine.

by_arm <- split(adsl, adsl$ARM)
names(by_arm)                                  #> "Drug A" "Placebo"

summaries <- lapply(by_arm, function(d) {
  data.frame(
    ARM      = d$ARM[1],
    n        = nrow(d),
    n_age    = sum(!is.na(d$AGE)),
    mean_age = mean(d$AGE, na.rm = TRUE),
    sd_age   = sd(d$AGE, na.rm = TRUE),
    n_female = sum(d$SEX == "F"),
    stringsAsFactors = FALSE
  )
})

result <- do.call(rbind, summaries)
rownames(result) <- NULL
result
#>       ARM n n_age mean_age   sd_age n_female
#> 1  Drug A 3     2    69.00 4.242641        1
#> 2 Placebo 3     3    54.67 19.13985        3

This pattern — split, lapply, do.call(rbind, ...) — handles any summary aggregate() cannot express, and keeps n (the population) separate from n_age (the non-missing count), which is exactly the distinction a clinical table needs.

Split on several variables with interaction() or a list:

split(adsl, list(adsl$ARM, adsl$SEX))                  # all combinations
split(adsl, list(adsl$ARM, adsl$SEX), drop = TRUE)     # observed only

`rowsum()` and `by()`

rowsum(adsl[, c("AGE", "WEIGHT")], group = adsl$ARM)   # fast group sums

by(adsl, adsl$ARM, function(d) summary(d$AGE))          # prints per group

rowsum() is very fast for sums specifically. by() is mostly for printing.

Merging

merge() is base R’s join, and it needs care.

merge(adsl, adae, by = "USUBJID")                      # INNER join (default!)
merge(adsl, adae, by = "USUBJID", all.x = TRUE)        # left join
merge(adsl, adae, by = "USUBJID", all.y = TRUE)        # right join
merge(adsl, adae, by = "USUBJID", all = TRUE)          # full join

merge(adsl, sites, by.x = "SITEID", by.y = "site_number")   # different names
merge(adlb, ranges, by = c("PARAMCD", "SEX"))               # multiple keys
merge(a, b, by = "USUBJID", suffixes = c(".adsl", ".adae"))

Three things merge() does that left_join() does not

1. It defaults to an inner join. dplyr makes you name the join type. merge() silently drops unmatched rows unless you pass all.x = TRUE. A subject with no adverse events disappears from the result, and the row count quietly falls.

2. It reorders rows. merge() sorts by the key. left_join() preserves the order of the left table. If anything downstream depends on row order — a duplicated() call, a --SEQ derivation — the merge has silently changed it.

merge(adsl, adae, by = "USUBJID", all.x = TRUE, sort = FALSE)

sort = FALSE helps but does not guarantee the original order.

3. It does not check cardinality. This is the serious one. A merge on a non-unique key multiplies rows, and nothing warns you.

nrow(adsl)                                               #> [1] 6
m <- merge(adsl, adae, by = "USUBJID", all.x = TRUE)
nrow(m)                                                  #> [1] 23

dplyr::left_join() offers relationship = "many-to-one" and unmatched = "error". merge() offers nothing. You must assert it yourself:

safe_merge <- function(x, y, by, expect = c("many-to-one", "one-to-one")) {
  expect <- match.arg(expect)

  if (anyDuplicated(y[, by, drop = FALSE])) {
    stop("Right-hand side has duplicate keys on ", paste(by, collapse = ", "),
         " — the merge would multiply rows.", call. = FALSE)
  }
  if (expect == "one-to-one" && anyDuplicated(x[, by, drop = FALSE])) {
    stop("Left-hand side has duplicate keys; expected one-to-one.", call. = FALSE)
  }

  out <- merge(x, y, by = by, all.x = TRUE, sort = FALSE)

  if (nrow(out) != nrow(x)) {
    stop("Merge changed the row count: ", nrow(x), " -> ", nrow(out),
         call. = FALSE)
  }
  out
}

That function is worth having in every base-R study utility file. The row-count assertion alone catches the most common way a clinical analysis produces silently wrong numbers.

Diagnosing a merge:

setdiff(adsl$USUBJID, adae$USUBJID)      # in ADSL, no events
setdiff(adae$USUBJID, adsl$USUBJID)      # events, not in ADSL — a data error
intersect(adsl$USUBJID, adae$USUBJID)

# Which keys are duplicated on the right?
dup <- adae$USUBJID[duplicated(adae$USUBJID)]
unique(dup)

Merging several data frames:

Reduce(function(a, b) merge(a, b, by = "USUBJID", all.x = TRUE),
       list(adsl, exposure, disposition))

Stacking rows

rbind(df1, df2)                       # requires IDENTICAL column names
rbind(df1, df2, df3)
do.call(rbind, list_of_dataframes)

rbind() requires the same columns; bind_rows() does not

a <- data.frame(x = 1, y = 2)
b <- data.frame(x = 3)
rbind(a, b)
#> Error in rbind(deparse.level, ...) :
#>   numbers of columns of arguments do not match

dplyr::bind_rows() fills the gap with NA. In base R, align the columns first:

rbind_fill <- function(...) {
  dfs <- list(...)
  all_cols <- unique(unlist(lapply(dfs, names)))

  aligned <- lapply(dfs, function(d) {
    missing <- setdiff(all_cols, names(d))
    for (m in missing) d[[m]] <- NA
    d[, all_cols, drop = FALSE]
  })
  do.call(rbind, aligned)
}

rbind() also matches columns by name when given data frames, so a differently ordered second frame is handled correctly — unlike cbind().

cbind(df1, df2)          # by POSITION, no key matching, no checks

Never use cbind() to combine two datasets about the same subjects. It lines up row 1 with row 1 regardless of whether they describe the same person. Use merge() on a key.

Reshaping

reshape() is base R’s long/wide conversion. It is powerful and its interface is notoriously unfriendly — worth having a worked example to copy.

Wide to long:

wide <- data.frame(
  USUBJID = c("001", "002"),
  WEEK0   = c(120, 135),
  WEEK4   = c(118, 132),
  WEEK8   = c(115, 130)
)

long <- reshape(
  wide,
  direction = "long",
  varying   = c("WEEK0", "WEEK4", "WEEK8"),   # columns to stack
  v.names   = "SBP",                          # name of the value column
  timevar   = "WEEK",                         # name of the key column
  times     = c(0, 4, 8),                     # values for the key column
  idvar     = "USUBJID"
)

long <- long[order(long$USUBJID, long$WEEK), ]
rownames(long) <- NULL
long
#>   USUBJID WEEK SBP
#> 1     001    0 120
#> 2     001    4 118
#> 3     001    8 115
#> 4     002    0 135

Long to wide:

wide_again <- reshape(
  long,
  direction = "wide",
  idvar     = "USUBJID",
  timevar   = "WEEK",
  v.names   = "SBP"
)
names(wide_again)
#> [1] "USUBJID" "SBP.0"   "SBP.4"   "SBP.8"

Note the generated names use a dot separator. Clean them up:

names(wide_again) <- sub("^SBP\\.", "WEEK", names(wide_again))

reshape() leaves attributes behind

reshape() attaches a reshapeLong/reshapeWide attribute to its result, which it uses to reverse the operation. It also sets odd row names.

attributes(long)$reshapeLong        # the reshaping metadata
attr(long, "reshapeLong") <- NULL   # remove it
rownames(long) <- NULL

Harmless in analysis, but it will show up in a str() comparison during QC and confuse whoever is reading the diff.

Simpler alternatives for simple cases:

stack(wide[, c("WEEK0", "WEEK4", "WEEK8")])   # quick long format, loses the ID
unstack(long_df, values ~ ind)

t(as.matrix(df))                               # transpose a numeric frame

For anything beyond a single value column, reshape() is the only base option and the code above is the template.

Frequency tables

table(adsl$ARM)
#>  Drug A Placebo
#>       3       3

table(adsl$ARM, adsl$SEX)
table(adsl$ARM, adsl$SEX, dnn = c("Treatment", "Sex"))     # dimension names

table() drops NA and unobserved levels by default

table(adsl$AGEGR1)        # a category with no subjects simply does not appear
table(adsl$SEX)           # missing values excluded entirely

Both are wrong for a clinical table. A treatment arm with zero events must show 0, and missing values must be visible.

table(adsl$SEX, useNA = "ifany")      # show NA when present
table(adsl$SEX, useNA = "always")     # always show an NA row

# Zero counts: make it a factor with all expected levels
adsl$AGEGR1 <- factor(adsl$AGEGR1, levels = c("<18", "18-64", "65-79", ">=80"))
table(adsl$AGEGR1)
#>   <18 18-64 65-79  >=80
#>     0     2     2     2      <- the zero now appears

A blank cell reads as not evaluated; 0 reads as evaluated, none occurred. Declaring the factor levels is what guarantees the row exists.

Percentages and margins:

tab <- table(adsl$ARM, adsl$SEX)

prop.table(tab)                 # of the grand total
prop.table(tab, margin = 1)     # row percentages
prop.table(tab, margin = 2)     # column percentages
round(100 * prop.table(tab, 1), 1)

addmargins(tab)                 # adds row and column totals
margin.table(tab, 1)            # row totals only

xtabs() gives a formula interface and handles pre-aggregated counts:

xtabs(~ ARM + SEX, data = adsl)
xtabs(N ~ ARM + SEX, data = counts_df)     # weight by an existing count column
ftable(xtabs(~ ARM + SEX + AGEGR1, data = adsl))   # flatten 3+ dimensions

Converting a table to a data frame for further work:

as.data.frame(tab)
#>      Var1 Var2 Freq
#> 1  Drug A    F    1
#> 2 Placebo    F    3

as.data.frame(tab, responseName = "n")

Missing values

is.na(adsl$AGE)
sum(is.na(adsl$AGE))                  # count per column
colSums(is.na(adsl))                  # count for every column at once
complete.cases(adsl)                  # TRUE where the whole row is complete
adsl[complete.cases(adsl), ]          # listwise deletion
adsl[complete.cases(adsl[, c("AGE", "SEX")]), ]   # on selected columns
na.omit(adsl)                         # same, but attaches an attribute

colSums(is.na(df)) is the fastest way to see where the gaps are in a new dataset.

# Replace a missing-value sentinel
adsl$AGE[adsl$AGE %in% c(-1, -999)] <- NA

# Fill with a value
adsl$AGE[is.na(adsl$AGE)] <- median(adsl$AGE, na.rm = TRUE)

complete.cases() is an analysis decision

Dropping incomplete rows is listwise deletion. If missingness relates to the outcome — which in a trial it usually does — the result is biased.

Report what was removed rather than removing it silently:

n_before <- nrow(adsl)
analysis <- adsl[complete.cases(adsl[, c("AGE", "AVAL")]), ]
message(sprintf("Excluded %d of %d records (%.1f%%) for missing data",
                n_before - nrow(analysis), n_before,
                100 * (n_before - nrow(analysis)) / n_before))

A complete worked example

A demographics summary table, built with no packages at all.

#-------------------------------------------------------------------------------
# Demographics summary — base R only
# No package dependencies beyond stats and utils (both ship with R)
#-------------------------------------------------------------------------------

summarise_demographics <- function(data,
                                   arm_var   = "TRT01P",
                                   arm_levels = NULL,
                                   pop_filter = NULL) {

  stopifnot(is.data.frame(data), arm_var %in% names(data))

  # --- 1. Restrict to the analysis population ------------------------------
  if (!is.null(pop_filter)) {
    keep <- which(data[[pop_filter]] == "Y")       # which() drops NA
    data <- data[keep, , drop = FALSE]
  }
  if (nrow(data) == 0) stop("No records in the analysis population.")

  # --- 2. Fix the treatment order, including a Total column ----------------
  if (is.null(arm_levels)) arm_levels <- sort(unique(data[[arm_var]]))

  data$.arm <- factor(data[[arm_var]], levels = arm_levels)
  if (anyNA(data$.arm)) {
    stop("Values of ", arm_var, " outside the declared levels: ",
         paste(setdiff(unique(data[[arm_var]]), arm_levels), collapse = ", "))
  }

  # Duplicate every row under "Total" so summaries include it automatically
  total <- data
  total$.arm <- "Total"
  combined <- rbind(data, total)
  combined$.arm <- factor(combined$.arm, levels = c(arm_levels, "Total"))

  columns <- levels(combined$.arm)

  # --- 3. Big N: the population denominator --------------------------------
  big_n <- table(combined$.arm)

  # --- 4. Helpers -----------------------------------------------------------
  cont_rows <- function(var, label, digits = 1) {
    stats <- tapply(combined[[var]], combined$.arm, function(x) {
      x <- x[!is.na(x)]
      if (length(x) == 0) return(c(n = 0, mean = NA, sd = NA,
                                   median = NA, min = NA, max = NA))
      c(n = length(x), mean = mean(x), sd = sd(x),
        median = median(x), min = min(x), max = max(x))
    }, simplify = FALSE)

    fmt <- function(f) vapply(columns, function(a) {
      s <- stats[[a]]
      if (is.na(s["mean"])) "-" else f(s)
    }, character(1))

    rbind(
      c(label, rep("", length(columns))),
      c("  n",            fmt(\(s) sprintf("%d", s["n"]))),
      c("  Mean (SD)",    fmt(\(s) sprintf("%.*f (%.*f)", digits, s["mean"],
                                           digits + 1, s["sd"]))),
      c("  Median",       fmt(\(s) sprintf("%.*f", digits, s["median"]))),
      c("  Min, Max",     fmt(\(s) sprintf("%.0f, %.0f", s["min"], s["max"])))
    )
  }

  cat_rows <- function(var, label, levels_wanted = NULL) {
    v <- combined[[var]]
    if (is.null(levels_wanted)) levels_wanted <- sort(unique(v[!is.na(v)]))
    v <- factor(v, levels = levels_wanted)

    counts <- table(v, combined$.arm)            # factor levels -> zero rows kept

    body <- t(vapply(levels_wanted, function(lv) {
      vapply(columns, function(a) {
        n <- counts[lv, a]
        sprintf("%d (%.1f)", n, 100 * n / big_n[a])
      }, character(1))
    }, character(length(columns))))

    rbind(
      c(paste0(label, ", n (%)"), rep("", length(columns))),
      cbind(paste0("  ", levels_wanted), body)
    )
  }

  # --- 5. Assemble ----------------------------------------------------------
  out <- rbind(
    cont_rows("AGE", "Age (years)"),
    cat_rows("AGEGR1", "Age group", c("<65", "65-80", ">80")),
    cat_rows("SEX", "Sex", c("F", "M")),
    cont_rows("WEIGHT", "Weight (kg)")
  )

  out <- as.data.frame(out, stringsAsFactors = FALSE)
  names(out) <- c("Characteristic", columns)
  rownames(out) <- NULL

  attr(out, "big_n") <- big_n
  out
}

Using it:

tbl <- summarise_demographics(
  adsl,
  arm_var    = "TRT01P",
  arm_levels = c("Placebo", "Drug A 50mg", "Drug A 100mg"),
  pop_filter = "SAFFL"
)

print(tbl, right = FALSE)
#>   Characteristic     Placebo      Drug A 50mg   Drug A 100mg  Total
#> 1 Age (years)
#> 2   n                86           84            84            254
#> 3   Mean (SD)        75.2 (8.59)  75.7 (8.29)   74.4 (7.89)   75.1 (8.25)
#> 4   Median           76.0         77.5          76.0          76.0
#> 5   Min, Max         52, 89       51, 88        56, 88        51, 89
#> 6 Age group, n (%)
#> 7   <65              14 (16.3)    11 (13.1)     8 (9.5)       33 (13.0)
#> 8   65-80            42 (48.8)    47 (56.0)     45 (53.6)     134 (52.8)
#> 9   >80              30 (34.9)    26 (31.0)     31 (36.9)     87 (34.3)

attr(tbl, "big_n")
#>      Placebo  Drug A 50mg Drug A 100mg        Total
#>           86           84           84          254

Five techniques from this lesson are doing the work:

which() for the population filter, so a missing SAFFL cannot produce a phantom row
factor(levels = ) on both the arm and each categorical variable, so every expected row appears even at zero and an unexpected arm value is an error rather than a silent NA column
Row duplication for the Total column, so every summary includes it with the correct denominator
table() on factors rather than tapply() for counts, because factor levels guarantee the zero cells exist
Big N from table(combined$.arm), kept separate from the per-variable n — so percentages divide by the population, not by the non-missing count

No packages. It would run unchanged in an environment with nothing but base R installed.

Base R and dplyr side by side

Task	Base R	dplyr
Filter rows	`df[which(df$x > 1), ]`	`filter(df, x > 1)`
Select columns	`df[, c("a", "b")]`	`select(df, a, b)`
Drop a column	`df$a <- NULL`	`select(df, -a)`
Add a column	`df$b <- df$a * 2`	`mutate(df, b = a * 2)`
Add, referencing new	`within(df, {...})`	`mutate(df, ...)`
Rename	`names(df)[2] <- "new"`	`rename(df, new = old)`
Sort	`df[order(df$x), ]`	`arrange(df, x)`
Sort descending	`df[order(-df$x), ]`	`arrange(df, desc(x))`
Distinct rows	`unique(df)`	`distinct(df)`
First per group	`df[!duplicated(df$id), ]`	`slice_head(n = 1, by = id)`
Last per group	`df[!duplicated(df$id, fromLast = TRUE), ]`	`slice_tail(n = 1, by = id)`
Group summary	`aggregate(y ~ g, df, mean)`	`summarise(df, mean(y), .by = g)`
Group + keep rows	`ave(df$y, df$g, FUN = mean)`	`mutate(m = mean(y), .by = g)`
Count	`table(df$g)`	`count(df, g)`
Cross-tabulate	`table(df$a, df$b)`	`count(df, a, b)`
Left join	`merge(a, b, all.x = TRUE)`	`left_join(a, b)`
Inner join	`merge(a, b)`	`inner_join(a, b)`
Anti join	`a[!a$id %in% b$id, ]`	`anti_join(a, b)`
Stack rows	`rbind(a, b)`	`bind_rows(a, b)`
Wide to long	`reshape(direction = "long")`	`pivot_longer()`
Long to wide	`reshape(direction = "wide")`	`pivot_wider()`
Complete cases	`df[complete.cases(df), ]`	`drop_na(df)`
Column types	`sapply(df, class)`	`glimpse(df)`
Missing per column	`colSums(is.na(df))`	`summarise(across(everything(), ~ sum(is.na(.))))`
Row count	`nrow(df)`	`nrow(df)`

Where dplyr is genuinely safer

Not a style preference — these are correctness differences:

Behaviour	Base R	dplyr
`NA` in a filter condition	Creates a phantom row	Row excluded
Single-column selection	Returns a vector	Returns a data frame
Default join type	Inner — drops rows silently	Must be named
Join cardinality	Unchecked	`relationship =`
Unmatched join keys	Silent	`unmatched = "error"`
Row order after a join	Changed	Preserved
Partial name matching with `$`	Matches silently	Warns
Recycling in assignment	Permissive	Length 1 or n only

Every one of those has produced a wrong number in a real analysis. If dplyr is available, use it. If it is not, the mitigations in this lesson — which(), drop = FALSE, all.x = TRUE, and an explicit row-count assertion after every merge — close the gaps.

Performance

For study-sized data, base R and dplyr are both fast enough and the difference is irrelevant. Where it matters:

# Faster than aggregate for sums specifically
rowsum(df[, numeric_cols], group = df$id)

# Vectorised column operations beat apply
colMeans(df[, numeric_cols])            # not apply(df, 2, mean)
rowSums(df[, c("a", "b", "c")])         # not apply(df, 1, sum)

# Pre-allocate rather than growing
out <- vector("list", length(files))
for (i in seq_along(files)) out[[i]] <- read.csv(files[i])
combined <- do.call(rbind, out)

The single biggest base-R performance mistake is growing an object in a loop — result <- rbind(result, new_row) reallocates and copies the whole data frame every iteration. Collect into a list and rbind once at the end.

Common mistakes

Mistake	Consequence	Fix
`df[df$x > 1, ]` with `NA` present	Phantom `NA` row	`df[which(df$x > 1), ]`
`df[, "col"]` in a function	Returns a vector, not a frame	`drop = FALSE`
`merge()` without `all.x`	Inner join; rows silently dropped	`all.x = TRUE`
`merge()` without a row-count check	Silent row multiplication	Assert `nrow()` after
Relying on row order after `merge()`	Order changed	`sort = FALSE`, then re-sort
`aggregate(y ~ g, ...)` with missing `y`	Rows dropped; wrong denominator	`na.action = na.pass` or `by =`
`table()` without factor levels	Zero categories vanish	`factor(levels = )`
`table()` without `useNA`	Missing values invisible	`useNA = "ifany"`
`duplicated()` after a non-deterministic sort	Different record kept between runs	Sort on enough keys
Forgetting `na.last = FALSE` vs SAS	QC differences on missing keys	Agree the convention
`rbind()` on mismatched columns	Error	Align columns first
`cbind()` to combine datasets	Rows misaligned by position	`merge()` on a key
`within()` and column order	Columns in reverse order	Reorder explicitly
Growing a data frame in a loop	Quadratic time	Collect into a list, `rbind` once
`apply(df, 1, ...)`	Coerces every row to character	`ave()`, `rowSums()`, or `split`
Not resetting row names after subsetting	Confusing indices later	`rownames(df) <- NULL`

Exercise 7.1 — Baseline and change in base R

Given a long-format lab data frame with USUBJID, PARAMCD, AVISITN, ADT and AVAL, derive BASE (the last value at or before visit 0), CHG and PCHG, plus a baseline flag ABLFL — using no packages. Verify the result.

Show solution

derive_baseline_change <- function(adlb) {
  required <- c("USUBJID", "PARAMCD", "AVISITN", "ADT", "AVAL")
  missing  <- setdiff(required, names(adlb))
  if (length(missing) > 0) {
    stop("Missing required column(s): ", paste(missing, collapse = ", "),
         call. = FALSE)
  }

  n_in <- nrow(adlb)

  # --- 1. Deterministic sort ------------------------------------------------
  # Sort ASCENDING by date so the LAST pre-dose record is the latest one.
  adlb <- adlb[order(adlb$USUBJID, adlb$PARAMCD, adlb$ADT, adlb$AVISITN), ]
  rownames(adlb) <- NULL

  # --- 2. Flag the baseline record -----------------------------------------
  # Candidates: non-missing AVAL at or before visit 0
  is_pre <- !is.na(adlb$AVAL) & !is.na(adlb$AVISITN) & adlb$AVISITN <= 0

  key <- paste(adlb$USUBJID, adlb$PARAMCD, sep = "|")

  adlb$ABLFL <- NA_character_

  # Among the pre-dose rows only, the LAST one per key is baseline.
  # fromLast = TRUE marks all but the final occurrence as duplicated.
  pre_idx <- which(is_pre)
  if (length(pre_idx) > 0) {
    last_pre <- pre_idx[!duplicated(key[pre_idx], fromLast = TRUE)]
    adlb$ABLFL[last_pre] <- "Y"
  }

  # --- 3. Broadcast the baseline value to every row -----------------------
  base_lookup <- adlb$AVAL[!is.na(adlb$ABLFL)]
  names(base_lookup) <- key[!is.na(adlb$ABLFL)]

  adlb$BASE <- unname(base_lookup[key])

  # --- 4. Change and percent change ---------------------------------------
  adlb$CHG <- adlb$AVAL - adlb$BASE

  adlb$PCHG <- ifelse(
    !is.na(adlb$BASE) & adlb$BASE != 0,
    100 * (adlb$AVAL - adlb$BASE) / adlb$BASE,
    NA_real_
  )

  # --- 5. Postconditions ---------------------------------------------------
  stopifnot(
    "row count changed"            = nrow(adlb) == n_in,
    "more than one baseline flagged per key" =
      all(table(key[!is.na(adlb$ABLFL)]) == 1),
    "baseline flagged after visit 0" =
      all(adlb$AVISITN[!is.na(adlb$ABLFL)] <= 0),
    "PCHG contains infinity"       = !any(is.infinite(adlb$PCHG))
  )

  adlb
}

Verification on a worked case:

adlb <- data.frame(
  USUBJID = c("001","001","001","001", "002","002", "003"),
  PARAMCD = c("ALT","ALT","ALT","ALT", "ALT","ALT", "ALT"),
  AVISITN = c(  -1,    0,    4,   12,     0,    4,      4),
  ADT     = as.Date(c("2026-03-10","2026-03-14","2026-04-11","2026-06-06",
                      "2026-03-20","2026-04-17","2026-04-01")),
  AVAL    = c(  30,   32,   38,   45,     28,   35,     40),
  stringsAsFactors = FALSE
)

res <- derive_baseline_change(adlb)
res[, c("USUBJID","AVISITN","AVAL","ABLFL","BASE","CHG","PCHG")]
#>   USUBJID AVISITN AVAL ABLFL BASE CHG     PCHG
#> 1     001      -1   30  <NA>   32  -2  -6.2500
#> 2     001       0   32     Y   32   0   0.0000
#> 3     001       4   38  <NA>   32   6  18.7500
#> 4     001      12   45  <NA>   32  13  40.6250
#> 5     002       0   28     Y   28   0   0.0000
#> 6     002       4   35  <NA>   28   7  25.0000
#> 7     003       4   40  <NA>   NA  NA      NA

# Subjects with no baseline — expected, and worth reporting
no_base <- unique(res$USUBJID[is.na(res$BASE)])
if (length(no_base) > 0) {
  message(sprintf("%d subject(s) with no baseline: %s",
                  length(no_base), paste(no_base, collapse = ", ")))
}
#> 1 subject(s) with no baseline: 003

Six decisions worth explaining

Subject 001 has two pre-dose records (visits −1 and 0). The ascending sort plus duplicated(fromLast = TRUE) selects visit 0 — the latest qualifying record, which is the standard convention. Sorting descending and taking the first would also work; mixing the two is how a QC difference arises.

which(is_pre) then duplicated() on the subset. Flagging cannot be done on the full data frame, because duplicated() would consider post-baseline rows. Restricting to the candidate indices first, then mapping back, is the base-R way to express restrict_derivation() from admiral.

A named-vector lookup for BASE. base_lookup[key] broadcasts the baseline to every row in one expression — no merge, so no risk of row multiplication and no row reordering. This is the named-vector pattern from lesson 3, and it is faster than a merge as well as safer.

ave() would also work for the simple case:

adlb$BASE <- ave(adlb$AVAL, key, FUN = function(x) x[1])

but only if baseline is always the first row, which requires the pre-dose restriction to have already been applied by filtering. The lookup version handles the general case.

The BASE != 0 guard on PCHG. Without it, division by zero gives Inf, which then propagates through every mean and standard deviation to produce a summary table full of Inf. The final stopifnot asserts it did not happen.

Subject 003 has no pre-dose record, so BASE, CHG and PCHG are all NA. That is correct behaviour, and reporting the count is the right move — how those subjects are handled in the analysis is a statistical decision, not a programming one.

Exercise 7.2 — An adverse event incidence table in base R

Build a subject-incidence AE table: rows for “Any TEAE” and each preferred term, columns per treatment arm plus Total, showing n (%) with correct denominators and explicit zeros. Restrict to terms occurring in at least 5% of subjects in any arm. No packages.

Show solution

ae_incidence_table <- function(adae, adsl,
                               arm_levels,
                               threshold_pct = 5,
                               subject_var = "USUBJID",
                               arm_var     = "TRT01A",
                               term_var    = "AEDECOD") {

  # --- 1. Denominators from the POPULATION, not the AE data ---------------
  pop <- adsl[which(adsl$SAFFL == "Y"), , drop = FALSE]
  pop$.arm <- factor(pop[[arm_var]], levels = arm_levels)
  if (anyNA(pop$.arm)) {
    stop("Arm values outside the declared levels: ",
         paste(setdiff(unique(pop[[arm_var]]), arm_levels), collapse = ", "),
         call. = FALSE)
  }

  big_n <- c(table(pop$.arm), Total = nrow(pop))
  columns <- names(big_n)

  # --- 2. Subject incidence, not event counts -----------------------------
  teae <- adae[which(adae$TRTEMFL == "Y" &
                     adae[[subject_var]] %in% pop[[subject_var]]), , drop = FALSE]

  # One row per subject per term
  sub_term <- unique(teae[, c(subject_var, arm_var, term_var)])

  # One row per subject, for the "Any TEAE" row
  sub_any  <- unique(teae[, c(subject_var, arm_var)])

  # --- 3. Counts, with Total via row duplication --------------------------
  with_total <- function(d) {
    t2 <- d; t2[[arm_var]] <- "Total"
    out <- rbind(d, t2)
    out$.arm <- factor(out[[arm_var]], levels = columns)
    out
  }

  st  <- with_total(sub_term)
  any <- with_total(sub_any)

  # Factor on the term keeps every requested row, including zeros
  all_terms <- sort(unique(sub_term[[term_var]]))
  st$.term  <- factor(st[[term_var]], levels = all_terms)

  counts   <- table(st$.term, st$.arm)          # terms x arms
  any_counts <- table(any$.arm)

  # --- 4. Apply the >=5% threshold, excluding Total ------------------------
  arm_only <- setdiff(columns, "Total")
  pct <- sweep(counts[, arm_only, drop = FALSE], 2, big_n[arm_only], "/") * 100
  keep <- rownames(counts)[apply(pct, 1, max) >= threshold_pct]

  if (length(keep) == 0) {
    warning("No terms meet the ", threshold_pct, "% threshold.", call. = FALSE)
  }

  # --- 5. Order by descending Total frequency -----------------------------
  keep <- keep[order(-counts[keep, "Total"], keep)]

  # --- 6. Format ----------------------------------------------------------
  cell <- function(n, arm) sprintf("%d (%.1f)", n, 100 * n / big_n[arm])

  any_row <- c("Subjects with any TEAE",
               vapply(columns, function(a) cell(any_counts[a], a), character(1)))

  term_rows <- t(vapply(keep, function(tm) {
    vapply(columns, function(a) cell(counts[tm, a], a), character(1))
  }, character(length(columns))))

  out <- rbind(any_row, cbind(paste0("  ", keep), term_rows))

  out <- as.data.frame(out, stringsAsFactors = FALSE)
  names(out) <- c("Preferred Term", columns)
  rownames(out) <- NULL

  attr(out, "big_n") <- big_n
  attr(out, "n_terms_total") <- length(all_terms)
  attr(out, "n_terms_shown") <- length(keep)
  out
}

tbl <- ae_incidence_table(
  adae, adsl,
  arm_levels = c("Placebo", "Xanomeline Low Dose", "Xanomeline High Dose")
)

print(tbl, right = FALSE)
#>   Preferred Term              Placebo    Xanomeline Low Dose  ...  Total
#> 1 Subjects with any TEAE      69 (80.2)  77 (91.7)                 226 (89.0)
#> 2   APPLICATION SITE PRURITUS  3 (3.5)   27 (32.1)                  54 (21.3)
#> 3   APPLICATION SITE ERYTHEMA  3 (3.5)   24 (28.6)                  50 (19.7)
#> 4   DIARRHOEA                  9 (10.5)   5 (6.0)                   20 (7.9)

attr(tbl, "big_n")
#>              Placebo  Xanomeline Low Dose Xanomeline High Dose  Total
#>                   86                   84                   84    254

sprintf("Showing %d of %d terms", attr(tbl, "n_terms_shown"),
        attr(tbl, "n_terms_total"))
#> [1] "Showing 23 of 242 terms"

Seven points that make this correct

Denominators come from ADSL, not ADAE. big_n is computed from the safety population. Using table(adae$TRT01A) would count events, and the percentages would be meaningless. This is the single most common error in AE table programming.

unique() on subject + term is what makes it incidence. A subject with three headaches appears once. Without it, percentages can exceed 100%. If your ADAE carries the AOCCPFL flag, filter on that instead — the logic then lives in the dataset where it is tested once, per ADaM programming.

The “Any TEAE” row uses a separate unique() on subject and arm only. A subject with events in three different terms counts once in that row and three times across the term rows, so the “Any” count is genuinely less than the column sum. Summing the term rows to get the total is wrong.

factor() on the term keeps every level in table(), so a term with zero subjects in one arm yields 0 (0.0) rather than a missing cell.

The threshold excludes Total. arm_only drops it before computing the maximum, because the convention is “≥5% in any treatment group”. Including Total would admit terms that are common overall but rare in every arm.

sweep(..., 2, big_n[arm_only], "/") divides each column by its own denominator. counts / big_n would recycle the denominator down the columns instead of across them — a silent, wrong-numbers bug of exactly the kind lesson 3 warns about.

Ties broken by name. order(-counts[keep, "Total"], keep) sorts by descending frequency and then alphabetically, so the output is reproducible when two terms have the same count.

Base-R functions used: which, factor, table, unique, rbind, sweep, apply, vapply, order, sprintf, setdiff, t, cbind, as.data.frame, structure via attr. Nothing outside base.

Recap

Base R is the right choice in restricted environments, in packages, and for reading legacy code
df[df$x > 1, ] with NA present creates a phantom row — use which()
df[, "col"] returns a vector; pass drop = FALSE inside functions
merge() defaults to an inner join, reorders rows, and checks nothing — assert the row count
aggregate()’s formula interface drops missing values and corrupts denominators
ave() broadcasts a group statistic to every row — base R’s mutate(.by =)
duplicated(fromLast = TRUE) gives last-per-group; sort deterministically first
table() needs factor(levels = ) for zero counts and useNA for missing
A named vector is a lookup table — safer and faster than a merge
SAS sorts missing first, R sorts it last
Collect into a list and rbind once; never grow a data frame in a loop

Next: Dates, strings and factors.