Descriptive Statistics and Tests
|
This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production. The reference material consulted while preparing these pages is listed in the section bibliography. |
Checked against R 4.6.1 (2026-06-24) on 2026-10-20; R 4.6.2 is scheduled for 2026-10-29. Almost everything on this page is base R (the base, stats, utils and datasets packages that ship with R), and every example runs. One add-on package is used, in a single example: broom 1.0.13, to turn the result of a test into a one-row data frame. The lintr 3.4.0 default linters were used only to check the scripts. Describing a sample and testing a statement about it are the first things almost every analysis does. R has a function for each step: mean(), median(), quantile() and summary() describe a vector, table() and xtabs() count, cor() measures association, aggregate() and tapply() do the same by group, and the classical tests (t.test(), wilcox.test(), chisq.test() and a dozen more) all return the same kind of object, an "htest", that prints in one familiar layout and can be taken apart like any list. This page walks through them in that order: the summaries first, then counts and correlations, then the "htest" object, then each test with its assumptions, its real output and the warnings it raises, then the correction for many tests, p.adjust(), and last the limits of what a p-value can tell you.
Every example is a script that was run with Rscript on R 4.6.1 with warnings printed as they happen (options(warn = 1)), and the lines that start with #> are what it printed. The scripts also pass the default lintr linters. The names books, loans, shelf_today, loan_days and the book() constructor come from the running library-shelf scenario described on the section index: a script that needs them starts with source("scenario.R"). The scenario data are synthetic and small, 16 books and 60 loans, so the numbers below show how the functions behave and not what libraries are like. Two other kinds of data are used where the shelf has too little to say: the datasets penguins (added to the datasets package in R 4.5.0), airquality and sleep, which are always available without loading anything. Scripts that draw random numbers call set.seed() first, so that their output can be reproduced; the random-number generators themselves are the subject of Distributions and Random Numbers. Where a message contains quotation marks, they are the typographic quotes R prints in a UTF-8 locale. Test results print with a tab before the title of the test (Welch Two Sample t-test and so on); that tab is part of the real output and is kept below.
The page leans on four neighbours. The data frames the examples summarise are described in Data Frames, the factors that table() and summary() count in Factors, the missing values that decide so many results in Missing and Special Values, and the formula y ~ x that the test functions accept in Formulas. Fitting models, which is what most of statistics becomes after the first descriptive pass, is in Statistical Models, and drawing the histograms and box plots that go with these numbers is in Base Graphics and Devices.
The data used on this page
Most examples need the length of a loan and the size of the book that was borrowed, so the first script builds one data frame with a row per loan. It tidies the genre labels as the Formulas page does, counts the days of each loan (a loan that is still out counts up to shelf_today), marks a loan as overdue when it lasted more than loan_days (21) days, and joins the books to the loans by ISBN. Two loans cannot be used: one has an impossible loan date, which is NA, and one was returned before it was taken out, whose negative length is replaced by NA because a negative number of days is an error and not a short loan. One loan has an ISBN that matches no book and drops out of the join. That leaves 59 loans of 15 different books; the 16th book, "Rivers of Glass", has an empty ISBN and so can match no loan. Every later script that needs these objects repeats the same lines at its top, so that each script runs alone.
source("scenario.R")
# The 16 books: tidy the genre labels, as the formulas page does.
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
# Days a loan lasted; a loan still out counts up to shelf_today.
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan: impossible
loans$overdue <- loans$days > loan_days
# One row per loan, with the book's pages, year and genre attached.
loan_books <- merge(loans, books, by = "isbn")
dim(loan_books)
#> [1] 59 12
str(loan_books[c("days", "overdue", "genre", "pages", "year")])
#> 'data.frame': 59 obs. of 5 variables:
#> $ days : int 10 19 31 25 58 23 33 4 31 7 ...
#> $ overdue: logi FALSE FALSE TRUE TRUE TRUE TRUE ...
#> $ genre : chr "fiction" "fiction" "fiction" "fiction" ...
#> $ pages : int 214 214 214 214 214 214 214 214 198 198 ...
#> $ year : int 2019 2019 2019 2019 2019 2019 2019 2019 2021 2021 ...
colSums(is.na(loan_books[c("days", "overdue", "pages", "year")]))
#> days overdue pages year
#> 2 2 10 3Source: R documentation: merge
The result has missing values in four columns (days and overdue 2 each, pages 10 and year 3, counted over loans), and the next sections show what each function does with them. The page counts of the books come from the CSV file with two NA, for "A Map of Salt" and "Ten Winters", and every loan of one of those books inherits the NA.
Means, medians and spread
The mean, trimming and missing values
mean() is the generic function for the arithmetic mean, and its default method has the signature mean(x, trim = 0, na.rm = FALSE, ...). With trim between 0 and 0.5 it drops that fraction of the observations from each end of the sorted data before averaging, so trim = 0.5 is the median; a value outside that range is taken as the nearest endpoint. The mean of a logical vector is the proportion of TRUE, which is the cheapest way to compute a rate, and the mean of dates is a date (see Dates and Times).
source("scenario.R")
pages <- books$pages
pages
#> [1] 387 432 214 198 301 NA 96 280 352 256 88 410 175 322 NA 244
mean(pages)
#> [1] NA
mean(pages, na.rm = TRUE)
#> [1] 268.2143
# trim drops that fraction of the observations from each end before averaging.
p <- sort(pages, na.last = NA) # na.last = NA removes the NAs
length(p)
#> [1] 14
mean(p, trim = 0.1)
#> [1] 269.5833
mean(p[-c(1, length(p))]) # the same: one value removed from each end
#> [1] 269.5833
mean(p, trim = 0.5) == median(p)
#> [1] TRUE
# A logical vector is a 0/1 vector: its mean is a proportion.
mean(books$year > 2010, na.rm = TRUE)
#> [1] 0.6666667
# Dates and difftimes have means too (see the Dates and Times page).
mean(loans$loan_date, na.rm = TRUE)
#> [1] "2025-09-20"Source: R documentation: mean
The first call returns NA, and that is the default on purpose. The mean of the 14 known page counts is 268.2 and their 10 percent trimmed mean is 269.6, but two books have no page count, and the two unknown values make the average of all 16 unknown too. na.rm = TRUE is a decision to answer a different question, the mean of the books for which the number is known, and it should be a visible one. The next section comes back to why the default is FALSE.
When the input is not a number
mean() of something that is not numeric or logical does not stop with an error. Its help page says that NA_real_ is returned, with a warning. The usual cause is a column that was read as text, so the first cure is to look at class() and str() of the column, and the second is to convert with as.numeric() and look at the values that did not convert. The other summaries are stricter or differ: var() coerces and warns, sd() of a factor is an error, and median() of a character vector sorts it and picks the middle string without complaint, which is rarely what was wanted.
# A tiny helper: the text of the first warning or error an expression raises.
message_of <- \(expr) {
tryCatch(
{
expr
"no condition"
},
warning = \(w) paste("warning:", conditionMessage(w)),
error = \(e) paste("error:", conditionMessage(e))
)
}
# mean() of something that is not numeric: a warning and NA, not an error.
message_of(mean("a"))
#> [1] "warning: argument is not numeric or logical: returning NA"
suppressWarnings(mean("a"))
#> [1] NA
message_of(mean(factor(c("a", "b"))))
#> [1] "warning: argument is not numeric or logical: returning NA"
# The functions that need a number are stricter or coerce.
message_of(var("a"))
#> [1] "warning: NAs introduced by coercion"
message_of(sd(factor(c("a", "b"))))
#> [1] "error: Calling var(x) on a factor x is defunct.\n Use something like 'all(duplicated(x)[-1L])' to test for a constant vector."
# The cure is to convert on purpose, and look at what did not convert.
x <- c("12", "7", "n/a", "30")
nums <- suppressWarnings(as.numeric(x))
nums
#> [1] 12 7 NA 30
mean(nums, na.rm = TRUE)
#> [1] 16.33333
x[is.na(nums)]
#> [1] "n/a"Source: R documentation: mean
The text of the first warning, argument is not numeric or logical: returning NA, is the one to search for when a script prints NA after reading a CSV file. The messages of the other functions on this page are collected in the section Messages and what they mean near the end.
Why missing values propagate
An Introduction to R states the rule in its section on missing values: in general, any operation on an NA becomes an NA, because if the specification of an operation is incomplete, the result cannot be known. The summary functions follow the rule by default, and each has an argument to say that the missing values should be dropped first. The argument is not the same everywhere, and neither is the default. The example uses airquality, the daily air quality in New York from May to September 1973 (datasets help page), in which the ozone column has 37 missing values among its 153 days.
# airquality (datasets): daily New York air quality, May to September 1973.
colSums(is.na(airquality))
#> Ozone Solar.R Wind Temp Month Day
#> 37 7 0 0 0 0
nrow(airquality)
#> [1] 153
# Missing values propagate: the mean of anything that includes an unknown is
# unknown. That is a deliberate default, and na.rm = TRUE is an explicit choice.
mean(airquality$Ozone)
#> [1] NA
mean(airquality$Ozone, na.rm = TRUE)
#> [1] 42.12931
sum(!is.na(airquality$Ozone)) # how many values that mean is based on
#> [1] 116
# The defaults differ by function: know which one you are calling.
median(airquality$Ozone)
#> [1] NA
sd(airquality$Ozone, na.rm = TRUE)
#> [1] 32.98788
tryCatch(quantile(airquality$Ozone), error = \(e) conditionMessage(e))
#> [1] "missing values and NaN's not allowed if 'na.rm' is FALSE"
range(airquality$Ozone, na.rm = TRUE)
#> [1] 1 168
summary(airquality$Ozone) # counts the NAs instead of failing
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
#> 1.00 18.00 31.50 42.13 63.25 168.00 37
table(is.na(airquality$Ozone))
#>
#> FALSE TRUE
#> 116 37
# na.rm in a grouped summary: the mean per month, with and without.
tapply(airquality$Ozone, airquality$Month, mean)
#> 5 6 7 8 9
#> NA NA NA NA NA
tapply(airquality$Ozone, airquality$Month, mean, na.rm = TRUE) |> round(1)
#> 5 6 7 8 9
#> 23.6 29.4 59.1 60.0 31.4
tapply(airquality$Ozone, airquality$Month, \(x) sum(!is.na(x)))
#> 5 6 7 8 9
#> 26 9 26 26 29
# Two columns, two different sets of rows: na.rm works column by column,
# complete cases work row by row.
mean(airquality$Ozone, na.rm = TRUE)
#> [1] 42.12931
mean(airquality$Ozone[complete.cases(airquality)])
#> [1] 42.0991
nrow(na.omit(airquality))
#> [1] 111
# Correlation: use = decides which rows enter.
cor(airquality$Ozone, airquality$Temp)
#> [1] NA
cor(airquality$Ozone, airquality$Temp, use = "complete.obs")
#> [1] 0.6983603
# Dropping is not neutral: the missing days are not a random sample.
tapply(is.na(airquality$Ozone), airquality$Month, sum)
#> 5 6 7 8 9
#> 5 21 5 5 1The pattern in the output is worth reading closely. mean(), median() and range() return NA unless na.rm = TRUE; quantile() goes further and stops with an error, missing values and NaN's not allowed if 'na.rm' is FALSE; summary() counts the missing values in an extra column instead of failing; and fivenum() drops them by default (na.rm = TRUE), unlike nearly every other function. tapply() has no such argument: extra arguments such as na.rm = TRUE are passed on to the function it applies. The formula interfaces of aggregate() and xtabs() have an na.action instead (the former omits incomplete rows, the latter keeps them), and cor() and cov() have use. Dropping is also not neutral: the table of missing days per month at the end shows that June has 21 missing ozone values among its 30 days and September only one. A mean computed with na.rm = TRUE is a mean of the days someone measured, which may not be a typical sample of the days that happened.
The two ways of dropping rows give different data. na.rm works column by column, so the mean of Ozone uses 116 values, whereas na.omit() and complete.cases() work row by row: only the 111 days with all six readings present remain, and the mean of ozone over those days is a different number (42.10 against 42.13). cor(..., use = "complete.obs") makes the same choice as na.omit(), and use = "pairwise.complete.obs" makes the choice of na.rm for each pair of columns. Which is right depends on why the values are missing, a question about the data and not about R; Missing and Special Values covers NA, NaN and the other special values.
Variance, standard deviation and range
var() and sd() use the denominator n - 1, not n. The help page of cor() (which also documents var() and cov()) says that this gives an unbiased estimator of the variance for independent identically distributed observations, and that these functions return NA when there is only one observation. The script computes the sum of squares by hand with both denominators to show which one R uses; the n - 1 version is what a textbook calls the sample variance and a calculator often labels s, and the n version is the variance of the data treated as the whole population. sd() is the square root of var(). range() returns the minimum and the maximum, not the width; diff(range(x)) is the width, and range() of data with a missing value is c(NA, NA) without na.rm.
source("scenario.R")
pages <- books$pages[!is.na(books$pages)]
median(pages)
#> [1] 268
var(pages)
#> [1] 11661.1
sd(pages)
#> [1] 107.9866
# var() and sd() divide by n - 1, not n.
n <- length(pages)
sum((pages - mean(pages))^2) / (n - 1)
#> [1] 11661.1
sum((pages - mean(pages))^2) / n
#> [1] 10828.17
sd(pages) == sqrt(var(pages))
#> [1] TRUE
# One value has no spread to estimate: the n - 1 denominator is 0.
var(pages[1])
#> [1] NA
sd(numeric())
#> [1] NA
# range() is the minimum and the maximum; diff() turns it into a width.
range(pages)
#> [1] 88 432
diff(range(pages))
#> [1] 344
range(books$pages) # the NAs make both ends NA
#> [1] NA NA
range(books$pages, na.rm = TRUE)
#> [1] 88 432
# Coefficient of variation: a unit-free spread.
sd(pages) / mean(pages)
#> [1] 0.4026131Running totals, weights and scaling
cumsum(), cumprod(), cummax() and cummin() return a vector of the same length. A missing value stops them: everything after the first NA is NA for cumsum(), so the usual repair is to replace the missing values by the neutral element first. A fifth member of the family is upcoming: the NEWS of R-devel (the future R 4.7.0) lists a new function cumvar() to calculate cumulative variances, which R 4.6 does not have (R-devel NEWS). weighted.mean(x, w) divides by the sum of the weights, which therefore need not add up to one; with na.rm = TRUE it removes a missing x together with its weight. The example weights the page count of each book by how many times the book was borrowed, which answers a different question from the plain mean: how long is the book that an average loan involves, as opposed to the average book. scale() centres each column on its mean and divides it by its standard deviation (the n - 1 one); it returns a matrix with the centres and scales stored as the attributes "scaled:center" and "scaled:scale", which is why sd(z) of the result is 1 and the mean is zero up to rounding.
source("scenario.R")
pages <- books$pages
cumsum(pages) # an NA poisons everything after it
#> [1] 387 819 1033 1231 1532 NA NA NA NA NA NA NA NA NA NA
#> [16] NA
cumsum(replace(pages, is.na(pages), 0L))
#> [1] 387 819 1033 1231 1532 1532 1628 1908 2260 2516 2604 3014 3189 3511 3511
#> [16] 3755
cummax(c(3, 1, 4, 1, 5, 9, 2, 6))
#> [1] 3 3 4 4 5 9 9 9
cummin(c(3, 1, 4, 1, 5, 9, 2, 6))
#> [1] 3 1 1 1 1 1 1 1
cumprod(1:6)
#> [1] 1 2 6 24 120 720
# Running total of loans per month, in date order.
month <- format(loans$loan_date, "%Y-%m")
per_month <- table(month)
head(per_month, 4)
#> month
#> 2025-01 2025-02 2025-03 2025-04
#> 4 2 3 3
head(cumsum(per_month), 4)
#> 2025-01 2025-02 2025-03 2025-04
#> 4 6 9 12
# weighted.mean(): the weights need not add to one. Here, how many pages does
# the average loan have, as opposed to the average book?
n_loans <- tabulate(match(loans$isbn, normalise_isbn(books$isbn)), nrow(books))
n_loans
#> [1] 7 3 8 6 1 3 3 3 0 2 3 4 2 4 7 3
mean(books$pages, na.rm = TRUE)
#> [1] 268.2143
weighted.mean(books$pages, w = n_loans, na.rm = TRUE)
#> [1] 267.7551
weighted.mean(c(1, NA, 3), w = c(3, 1, 1), na.rm = TRUE)
#> [1] 1.5
# scale() centres and divides by the standard deviation, column by column.
z <- scale(books$pages)
head(z, 3)
#> [,1]
#> [1,] 1.1000043
#> [2,] 1.5167227
#> [3,] -0.5020464
attr(z, "scaled:center")
#> [1] 268.2143
round(mean(z, na.rm = TRUE), 10)
#> [1] 0
sd(z, na.rm = TRUE)
#> [1] 1Source: R documentation: weighted.mean
Quantiles
quantile() and its nine types
A sample quantile is not defined uniquely, and R implements nine definitions of it, selected by the type argument of quantile(). The help page describes them following Hyndman and Fan (1996): types 1 to 3 are step functions of the probability (type 1 is the inverse of the empirical distribution function, type 2 averages at the steps, type 3 takes the nearest even order statistic), and types 4 to 9 interpolate linearly between consecutive order statistics, each with a different plotting position. The default is type 7, the definition used by S and by R before 2.0.0, and it has stayed the default; type 6 is the one used by Minitab and by SPSS, and the help page says that Hyndman and Fan recommended type 8. This matters whenever a number computed in R has to match a number from another program: the 25 percent point of the 14 known page counts is 202.00 in R and 192.25 in a program that uses type 6, and neither is wrong.
Other facts from the help page that show up in the examples. probs must lie in [0, 1], but values up to 2e-14 outside are moved to the endpoint. names = FALSE drops the labels, and digits (default 7, added in R 4.1.0) is the precision used to format them: before R 4.1.0 the names depended on getOption("digits"). Types 1 and 3 can be used with dates and ordered factors. fuzz protects against floating-point error when the position is computed; since R 4.5.0 it is used by every type, with a default of zero for type 7 so that results stay the same.
source("scenario.R")
pages <- books$pages
# Missing values are an error unless you ask for them to be dropped.
tryCatch(quantile(pages), error = \(e) conditionMessage(e))
#> [1] "missing values and NaN's not allowed if 'na.rm' is FALSE"
quantile(pages, na.rm = TRUE)
#> 0% 25% 50% 75% 100%
#> 88.0 202.0 268.0 344.5 432.0
p <- pages[!is.na(pages)]
sort(p)
#> [1] 88 96 175 198 214 244 256 280 301 322 352 387 410 432
# probs picks the points; names = FALSE drops the labels.
quantile(p, probs = c(0.1, 0.5, 0.9))
#> 10% 50% 90%
#> 119.7 268.0 403.1
quantile(p, probs = c(0.1, 0.5, 0.9), names = FALSE)
#> [1] 119.7 268.0 403.1
quantile(p, probs = 0.333333333)
#> 33.33333%
#> 224
quantile(p, probs = 0.333333333, digits = 3)
#> 33.3%
#> 224
# Default type 7 and the type 6 of Minitab and SPSS give different answers.
quantile(p, c(0.25, 0.75), type = 7)
#> 25% 75%
#> 202.0 344.5
quantile(p, c(0.25, 0.75), type = 6)
#> 25% 75%
#> 192.25 360.75
# The nine types in one table: rows are the types, columns the probabilities.
probs <- c(0.1, 0.25, 0.5, 0.9)
by_type <- t(vapply(1:9, \(type) quantile(p, probs, type = type), numeric(4)))
rownames(by_type) <- paste("type", 1:9)
round(by_type, 1)
#> 10% 25% 50% 90%
#> type 1 96.0 198.0 256 410.0
#> type 2 96.0 198.0 268 410.0
#> type 3 88.0 198.0 256 410.0
#> type 4 91.2 186.5 256 400.8
#> type 5 95.2 198.0 268 412.2
#> type 6 92.0 192.2 268 421.0
#> type 7 119.7 202.0 268 403.1
#> type 8 94.1 196.1 268 415.1
#> type 9 94.4 196.6 268 414.4Source: R documentation: quantile
The nine rows of the table agree on the median of 14 values only partly (256 for types 1, 3 and 4, and 268 for the others), because with an even number of values the median falls between the 7th and 8th value, 256 and 280, and the types disagree on which to take. The disagreement is largest in the tails, where the types differ most: at the 90 percent point the answers run from 400.8 to 421.0. With a few hundred values the types are nearly the same; with 14 they are not, and that is one more reason to say which type was used when a quantile is reported.
fivenum(), summary() and quantile()
Three functions return "the quartiles", and they do not agree. fivenum() returns Tukey’s five numbers: the minimum, the lower hinge, the median, the upper hinge and the maximum. The hinges are the medians of the lower and upper halves of the sorted data, and the median belongs to both halves when the number of values is odd. quantile() with its default type 7 interpolates at the position 1 + (n - 1) p, and summary() of a numeric vector calls quantile() (its argument quantile.type is 7 by default) and adds the mean. With 14 values the first hinge is the 4th value (198) and the quartile of type 7 is a quarter of the way from the 4th to the 5th (202); the upper hinge is the 11th value (352) and the upper quartile is 344.5. The differences are small and systematic, and they are the reason that boxplot() and fivenum() draw the same box, whereas quantile() can give a slightly different one. The example computes the hinges by hand from the halves, and shows the positions that types 7 and 6 use.
source("scenario.R")
p <- books$pages[!is.na(books$pages)]
n <- length(p)
sort(p)
#> [1] 88 96 175 198 214 244 256 280 301 322 352 387 410 432
# Tukey's five numbers: minimum, lower hinge, median, upper hinge, maximum.
fivenum(p)
#> [1] 88 198 268 352 432
quantile(p, na.rm = TRUE) # type 7: positions 1 + (n - 1) * p
#> 0% 25% 50% 75% 100%
#> 88.0 202.0 268.0 344.5 432.0
summary(p)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 88.0 202.0 268.0 268.2 344.5 432.0
# The hinges are medians of the lower and upper halves (including the median
# when n is odd). With n = 14 that is the 4th and 11th values.
sorted <- sort(p)
median(sorted[1:7])
#> [1] 198
median(sorted[8:14])
#> [1] 352
# Where each type puts the 25% point: a fractional position in sorted.
1 + (n - 1) * 0.25 # type 7
#> [1] 4.25
(n + 1) * 0.25 # type 6
#> [1] 3.75
quantile(p, 0.25, type = 7)
#> 25%
#> 202
quantile(p, 0.25, type = 6)
#> 25%
#> 192.25
# For n = 7 the type 7 quartiles happen to land on the hinges; type 6 does not.
q <- c(1, 3, 5, 7, 9, 11, 13)
fivenum(q)
#> [1] 1 4 7 10 13
quantile(q, type = 7)
#> 0% 25% 50% 75% 100%
#> 1 4 7 10 13
quantile(q, type = 6)
#> 0% 25% 50% 75% 100%
#> 1 3 7 11 13Source: R documentation: fivenum
Robust spread: IQR() and mad()
The standard deviation is pulled by long tails, and so is the mean. median(), IQR() and mad() are the robust companions. IQR(x, na.rm = FALSE, type = 7) is the distance between the quartiles; the help page notes that it uses quantile() and not Tukey’s hinges, and that dividing by 1.349 gives a normal-consistent estimate of the standard deviation. mad(x) is the median of the absolute deviations from the median, multiplied by 1.4826 (constant) so that it estimates the standard deviation for normal data; constant = 1 gives the raw median. The loan lengths of the shelf have a long right tail (a few loans of 58 to 84 days), so the mean (25.4 days) is higher than the median (23) and the standard deviation (17.1) is nearly twice the scaled MAD (8.9). The same data give the outliers boxplot.stats() would draw: values beyond 1.5 times the IQR from the quartiles, which are Tukey’s fences.
source("scenario.R")
back <- loans$return_date
back[is.na(back)] <- shelf_today
days <- as.integer(back - loans$loan_date)
days <- days[!is.na(days) & days >= 0L]
sort(days)
#> [1] 1 4 4 7 8 10 10 10 10 10 11 12 14 16 16 16 17 17 18 19 20 20 21 21 21
#> [26] 21 22 23 23 23 24 24 24 25 25 26 26 26 27 27 27 27 28 28 28 29 31 31 33 36
#> [51] 36 37 44 58 65 71 81 84
# Mean and standard deviation are pulled by the long loans; the median is not.
mean(days)
#> [1] 25.39655
median(days)
#> [1] 23
sd(days)
#> [1] 17.08281
# IQR() is the distance between the quartiles; type selects the quantile type.
IQR(days)
#> [1] 12
diff(quantile(days, c(0.25, 0.75)))
#> 75%
#> 12
IQR(days, type = 6)
#> [1] 12
# mad() is the median absolute deviation, scaled to match sd() for normal data.
mad(days)
#> [1] 8.8956
mad(days, constant = 1) # the raw median of |x - median(x)|
#> [1] 6
median(abs(days - median(days)))
#> [1] 6
# The trimmed mean sits between the mean and the median.
mean(days, trim = 0.1)
#> [1] 22.70833
# Tukey's fences: the rule boxplot() uses to draw points beyond the whiskers.
fences <- quantile(days, c(0.25, 0.75)) + c(-1.5, 1.5) * IQR(days)
fences
#> 25% 75%
#> -2 46
days[days < fences[[1L]] | days > fences[[2L]]]
#> [1] 84 81 71 65 58
boxplot.stats(days)$out
#> [1] 84 81 71 65 58Source: R documentation: IQR
summary()
One variable
summary() is a generic function, and its default method gives a different kind of summary for each type of vector: the six numbers (minimum, quartiles, mean and maximum) plus a count of NAs for a numeric vector, the counts of FALSE, TRUE and NAs for a logical vector, the counts of the levels of a factor and, since R 4.6.0, for a character vector the number of values, the number of unique ones, the number that are blank and the range of the string lengths. Printing rounds to four significant digits. Since R 3.4.0 the rounding is done by print() and format() and not by summary() itself, so summary(x)[["Mean"]] is the full value; the digits argument rounds the printed values, and the zdigits argument of the printing method, new in R 4.5.0 and described in its NEWS entry as somewhat experimental, controls the extra digits kept before zapsmall(). The result is a classed vector ("summaryDefault"), which you can index by name.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
summary(books$pages) # numeric: quartiles, mean, and a count of NAs
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
#> 88.0 202.0 268.0 268.2 344.5 432.0 2
summary(books$year > 2010) # logical: counts of FALSE, TRUE and NA
#> Mode FALSE TRUE NAs
#> logical 5 10 1
summary(books$title) # character: counts and string lengths (R 4.6.0)
#> Length N.unique N.blank Min.nchar Max.nchar
#> 16 16 0 10 21
summary(books$title, character.method = NULL) # the layout of R before 4.6.0
#> Length Class Mode
#> 16 character character
summary(factor(books$genre)) # factor: the counts of each level
#> classic fiction mystery non-fiction poetry
#> 2 4 2 3 2
#> science fiction
#> 3
summary(c(a = 1.234567, b = 20000.5)) # printed to four significant digits
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.235 5001.051 10000.867 10000.867 15000.684 20000.500
summary(books$pages, digits = 2)
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
#> 88 200 270 270 340 430 2
# The result of summary() on a number is a classed vector you can index.
s <- summary(books$pages)
class(s)
#> [1] "summaryDefault" "table"
names(s)
#> [1] "Min." "1st Qu." "Median" "Mean" "3rd Qu." "Max." "NAs"
s[["Mean"]]
#> [1] 268.2143
unclass(s)
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
#> 88.0000 202.0000 268.0000 268.2143 344.5000 432.0000 2.0000Source: R documentation: summary
The label of the missing values changed in R 4.6.0: the NEWS entry says that the default method and the factor method now label them NAs and not NA's. A script, or a report, that matched the old label with a regular expression stops matching. The character summary changed in the same release, and character.method = NULL brings back the older layout (Length, Class, Mode) while character.method = "factor" counts the values as a factor would, which is what the next section uses in a data frame.
A whole data frame
The method for data frames summarises each column and lays the results side by side as a character table of class "table". A numeric column shows its six numbers, a factor its counts (at most maxsum rows, with the rest lumped into (Other)), and a character column, in R 4.6.0 and later, the string summary. A column of text is described, not counted, unless character.method = "factor" is given or the column is converted with factor(), which is what read.csv() did to text columns before stringsAsFactors = FALSE became the default in R 4.0.0, so older tutorials show counts for them (see Data Frames).
source("scenario.R")
summary(books)
#> title author isbn year genre
#> Length :16 Length :16 Length :16 Min. :1813 Length :16
#> N.unique :16 N.unique :10 N.unique :16 1st Qu.:2006 N.unique : 9
#> N.blank : 0 N.blank : 0 N.blank : 1 Median :2016 N.blank : 0
#> Min.nchar:10 Min.nchar:10 Min.nchar: 0 Mean :1999 Min.nchar: 6
#> Max.nchar:21 Max.nchar:18 Max.nchar:17 3rd Qu.:2020 Max.nchar:15
#> Max. :2023
#> NAs :1
#> pages
#> Min. : 88.0
#> 1st Qu.:202.0
#> Median :268.0
#> Mean :268.2
#> 3rd Qu.:344.5
#> Max. :432.0
#> NAs :2
# summary() of a data frame is a character table of class "table".
s <- summary(books[c("year", "pages")])
class(s)
#> [1] "table"
dim(s)
#> [1] 7 2
dimnames(s)[[2L]]
#> [1] " year" " pages"
# A character column gets the string summary above. To count the values, ask
# for the factor method, or convert the column to a factor first.
books$genre <- tolower(trimws(books$genre))
summary(books[c("genre", "pages")], character.method = "factor")
#> genre pages
#> classic :2 Min. : 88.0
#> fiction :4 1st Qu.:202.0
#> mystery :2 Median :268.0
#> non-fiction :3 Mean :268.2
#> poetry :2 3rd Qu.:344.5
#> science fiction:3 Max. :432.0
#> NAs :2
books$genre <- factor(books$genre)
summary(books[c("genre", "pages")], maxsum = 4)
#> genre pages
#> fiction :4 Min. : 88.0
#> non-fiction :3 1st Qu.:202.0
#> science fiction:3 Median :268.0
#> (Other) :6 Mean :268.2
#> 3rd Qu.:344.5
#> Max. :432.0
#> NAs :2
# summary() of the date columns of loans.
summary(loans[c("loan_date", "return_date")])
#> loan_date return_date
#> Min. :2025-01-14 Min. :2025-01-22
#> 1st Qu.:2025-06-10 1st Qu.:2025-05-27
#> Median :2025-10-02 Median :2025-09-26
#> Mean :2025-09-20 Mean :2025-09-09
#> 3rd Qu.:2026-01-07 3rd Qu.:2025-12-12
#> Max. :2026-05-28 Max. :2026-04-19
#> NAs :1 NAs :7Source: R documentation: summary
The real output of summary(books) is worth keeping: the genre column shows nine distinct strings (N.unique : 9), where the tidy data have six, and the isbn column has one blank (N.blank : 1). The new character summary shows data problems that a numeric summary of the same table would not: the quickest way to find messy text is to read N.unique and N.blank.
Counting
table(): counts, missing values and names
table() cross-classifies one or more factors (or anything that can be turned into one, including numbers and strings) and counts the combinations. The result is an integer array of class "table" with named dimensions. It can be indexed by level name, sorted with sort(), and converted to a data frame with as.data.frame(), whose columns are Var1, … and Freq unless responseName or the names of the dimensions say otherwise. Three arguments matter most. useNA takes "no" (the default), "ifany" (add a column of NA counts only if there are some) or "always" (add it even when it is zero). exclude lists levels to remove; if it does not contain NA and useNA is not given, it implies useNA = "ifany", so table(x, exclude = NULL) counts the missing values, which is the idiom of older code. dnn names the dimensions, which makes the output readable and fills the empty header line over a one-way table.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
tab <- table(books$genre)
tab
#>
#> classic fiction mystery non-fiction poetry
#> 2 4 2 3 2
#> science fiction
#> 3
class(tab)
#> [1] "table"
names(tab)
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
tab[["fiction"]]
#> [1] 4
sort(tab, decreasing = TRUE)
#>
#> fiction non-fiction science fiction classic mystery
#> 4 3 3 2 2
#> poetry
#> 2
as.data.frame(tab)
#> Var1 Freq
#> 1 classic 2
#> 2 fiction 4
#> 3 mystery 2
#> 4 non-fiction 3
#> 5 poetry 2
#> 6 science fiction 3
as.data.frame(tab, responseName = "n_books", stringsAsFactors = FALSE)
#> Var1 n_books
#> 1 classic 2
#> 2 fiction 4
#> 3 mystery 2
#> 4 non-fiction 3
#> 5 poetry 2
#> 6 science fiction 3
# NAs are dropped by default; useNA says what to do with them.
table(loans$overdue)
#>
#> FALSE TRUE
#> 26 32
table(loans$overdue, useNA = "ifany")
#>
#> FALSE TRUE <NA>
#> 26 32 2
table(loans$overdue, useNA = "always")
#>
#> FALSE TRUE <NA>
#> 26 32 2
table(books$pages > 300, useNA = "ifany")
#>
#> FALSE TRUE <NA>
#> 8 6 2
table(c("a", "b"), useNA = "always") # an NA count of 0 is still shown
#>
#> a b <NA>
#> 1 1 0
# dnn names the dimensions; exclude removes levels (and then useNA is "no").
table(loans$overdue, loans$member != "", dnn = c("overdue", "member known"))
#> member known
#> overdue FALSE TRUE
#> FALSE 1 25
#> TRUE 0 32
table(c("x", "y", "z", "x"), exclude = "z")
#>
#> x y
#> 2 1
table(c(1, 2, NA, 2), exclude = NULL)
#>
#> 1 2 <NA>
#> 1 2 1Source: R documentation: table
R 3.4.0 made table() consistent here: the NEWS entry says that table(x, exclude = NULL) sets useNA = "ifany" and that excluded values are no longer counted among the NA. The argument deparse.level controls how the names of the dimensions are made from the arguments when dnn is not given: 1 (the default) uses an argument only if it is a symbol, which is why table(loan_books$genre, loan_books$overdue) has empty names and table(genre, overdue) would not.
Proportions and margins: proportions(), marginSums() and addmargins()
A table of counts is rarely the final answer; the proportions are. proportions(x, margin = NULL) divides by the grand total, or, with margin = 1, by the row totals, or with margin = 2 by the column totals, so that the margin you name is the one that sums to 1. Which margin to use is a question about what is being compared: row proportions answer "of the loans of each genre, which share was overdue", column proportions answer "of the overdue loans, which share was of each genre". R 4.0.0 added proportions() and marginSums() to replace the unfortunately named prop.table() and margin.table(), in the words of the NEWS entry; they are drop-in replacements that also accept the name of a dimension as margin, and the old names stay as aliases, so old code keeps working. addmargins() appends the totals as a new row and column (or only to the margins named), and rowSums() and colSums() also work on a two-way table.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
by_genre <- table(
loan_books$genre, loan_books$overdue,
dnn = c("genre", "overdue")
)
by_genre
#> overdue
#> genre FALSE TRUE
#> classic 2 3
#> fiction 10 12
#> mystery 2 2
#> non-fiction 2 3
#> poetry 2 4
#> science fiction 8 7
# proportions() (R 4.0.0) replaces the oddly named prop.table(); margin is the
# dimension that adds up to 1: 1 = rows, 2 = columns.
proportions(by_genre)
#> overdue
#> genre FALSE TRUE
#> classic 0.03508772 0.05263158
#> fiction 0.17543860 0.21052632
#> mystery 0.03508772 0.03508772
#> non-fiction 0.03508772 0.05263158
#> poetry 0.03508772 0.07017544
#> science fiction 0.14035088 0.12280702
round(proportions(by_genre, margin = 1), 2)
#> overdue
#> genre FALSE TRUE
#> classic 0.40 0.60
#> fiction 0.45 0.55
#> mystery 0.50 0.50
#> non-fiction 0.40 0.60
#> poetry 0.33 0.67
#> science fiction 0.53 0.47
round(proportions(by_genre, margin = 2), 2)
#> overdue
#> genre FALSE TRUE
#> classic 0.08 0.10
#> fiction 0.38 0.39
#> mystery 0.08 0.06
#> non-fiction 0.08 0.10
#> poetry 0.08 0.13
#> science fiction 0.31 0.23
identical(prop.table(by_genre, 1), proportions(by_genre, 1))
#> [1] TRUE
# marginSums() (R 4.0.0) replaces margin.table(): the totals themselves.
marginSums(by_genre, margin = 1)
#> genre
#> classic fiction mystery non-fiction poetry
#> 5 22 4 5 6
#> science fiction
#> 15
marginSums(by_genre)
#> [1] 57
rowSums(by_genre)
#> classic fiction mystery non-fiction poetry
#> 5 22 4 5 6
#> science fiction
#> 15
# Margins can be named, since the dimensions of a table are named.
proportions(by_genre, margin = "genre")["fiction", ]
#> FALSE TRUE
#> 0.4545455 0.5454545
# addmargins() appends the totals as a new row and column.
addmargins(by_genre)
#> overdue
#> genre FALSE TRUE Sum
#> classic 2 3 5
#> fiction 10 12 22
#> mystery 2 2 4
#> non-fiction 2 3 5
#> poetry 2 4 6
#> science fiction 8 7 15
#> Sum 26 31 57
addmargins(by_genre, margin = 1)
#> overdue
#> genre FALSE TRUE
#> classic 2 3
#> fiction 10 12
#> mystery 2 2
#> non-fiction 2 3
#> poetry 2 4
#> science fiction 8 7
#> Sum 26 31
addmargins(proportions(by_genre, 1), 2)
#> overdue
#> genre FALSE TRUE Sum
#> classic 0.4000000 0.6000000 1.0000000
#> fiction 0.4545455 0.5454545 1.0000000
#> mystery 0.5000000 0.5000000 1.0000000
#> non-fiction 0.4000000 0.6000000 1.0000000
#> poetry 0.3333333 0.6666667 1.0000000
#> science fiction 0.5333333 0.4666667 1.0000000Source: R documentation: proportions
xtabs(), ftable() and the formula interface
xtabs() builds the same tables from a formula and a data frame: the variables on the right-hand side classify, and an optional left-hand side is summed in each cell, which makes it a grouped sum as well as a count. Its result has class c("xtabs", "table"), and the summary() method for tables runs the chi-squared test of independence and reports, in one line, the same warning that chisq.test() gives when the expected counts are small. The default na.action of xtabs() is na.pass; the left-hand side NA are then summed as NA unless na.rm = TRUE, and print.xtabs() prints an NA as an empty cell (na.print = ""), which is why the first table of days per genre has blanks. Since R 4.4.0, xtabs(Freq ~ .) consistently uses na.pass and sums with na.rm = FALSE (NEWS), and addNA = TRUE (R 3.4.0) counts the missing values of the classifying factors as a level of their own. ftable() flattens a table with three or more dimensions for printing, with row.vars and col.vars choosing which variables go down the side, and has a formula interface of its own.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# xtabs() counts with a formula: the right-hand side names the classifiers.
xtabs(~genre, data = loan_books)
#> genre
#> classic fiction mystery non-fiction poetry
#> 6 23 4 5 6
#> science fiction
#> 15
xt <- xtabs(~ genre + overdue, data = loan_books)
xt
#> overdue
#> genre FALSE TRUE
#> classic 2 3
#> fiction 10 12
#> mystery 2 2
#> non-fiction 2 3
#> poetry 2 4
#> science fiction 8 7
class(xt)
#> [1] "xtabs" "table"
summary(xt) # the chi-squared test of independence, in one line
#> Call: xtabs(formula = ~genre + overdue, data = loan_books)
#> Number of cases in table: 57
#> Number of factors: 2
#> Test for independence of all factors:
#> Chisq = 0.8834, df = 5, p-value = 0.9714
#> Chi-squared approximation may be incorrect
names(dimnames(xt))
#> [1] "genre" "overdue"
# A left-hand side sums a variable instead of counting rows: days on loan.
# One NA makes a whole cell NA (it prints blank) unless na.rm = TRUE.
xtabs(days ~ genre, data = loan_books)
#> genre
#> classic fiction mystery non-fiction poetry
#> 83 106 171
#> science fiction
#> 392
xtabs(days ~ genre, data = loan_books, na.rm = TRUE)
#> genre
#> classic fiction mystery non-fiction poetry
#> 107 588 83 106 171
#> science fiction
#> 392
# Missing values: xtabs() drops incomplete rows unless addNA = TRUE.
sum(xtabs(~overdue, data = loan_books))
#> [1] 57
xtabs(~overdue, data = loan_books, addNA = TRUE)
#> overdue
#> FALSE TRUE <NA>
#> 26 31 2
# A data frame of counts and back again.
counts <- as.data.frame(xt)
head(counts, 4)
#> genre overdue Freq
#> 1 classic FALSE 2
#> 2 fiction FALSE 10
#> 3 mystery FALSE 2
#> 4 non-fiction FALSE 2
xtabs(Freq ~ genre + overdue, data = counts)
#> overdue
#> genre FALSE TRUE
#> classic 2 3
#> fiction 10 12
#> mystery 2 2
#> non-fiction 2 3
#> poetry 2 4
#> science fiction 8 7
# ftable() flattens a table with three or more dimensions for printing.
loan_books$long_book <- loan_books$pages > 300
xt3 <- xtabs(~ genre + overdue + long_book, data = loan_books)
dim(xt3)
#> [1] 6 2 2
ftable(xt3)
#> long_book FALSE TRUE
#> genre overdue
#> classic FALSE 0 2
#> TRUE 2 1
#> fiction FALSE 8 0
#> TRUE 8 0
#> mystery FALSE 1 1
#> TRUE 2 0
#> non-fiction FALSE 1 0
#> TRUE 1 0
#> poetry FALSE 2 0
#> TRUE 4 0
#> science fiction FALSE 0 8
#> TRUE 0 7
ftable(xt3, row.vars = c("long_book", "genre"))
#> overdue FALSE TRUE
#> long_book genre
#> FALSE classic 0 2
#> fiction 8 8
#> mystery 1 2
#> non-fiction 1 1
#> poetry 2 4
#> science fiction 0 0
#> TRUE classic 2 1
#> fiction 0 0
#> mystery 1 0
#> non-fiction 0 0
#> poetry 0 0
#> science fiction 8 7
ftable(overdue ~ long_book, data = loan_books)
#> overdue FALSE TRUE
#> long_book
#> FALSE 12 17
#> TRUE 11 8Source: R documentation: xtabs
Binning numbers: cut(), tabulate() and findInterval()
To tabulate a number, turn it into a factor of intervals. cut(x, breaks) takes either a vector of break points or a number of equal-width intervals, closes them on the right by default (right = TRUE, so (150,300]), and labels them. include.lowest = TRUE closes the first interval on its left edge too, which is needed when a break point equals the minimum, as the example with 88 pages shows; dig.lab sets how many digits the labels use. labels names the bins and ordered_result = TRUE makes them an ordered factor. Using quantile() as the breaks makes bins with about the same number of observations. tabulate(bin, nbins) is the low-level counter: it counts the positive integers up to nbins and nothing else, ignores zeros and negative numbers, and is a good deal faster than table() for a large vector of small integers. findInterval() gives the bin number directly, without making a factor.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# cut() turns a number into an ordered set of bins; table() counts them.
bins <- cut(books$pages, breaks = c(0, 150, 300, 450))
levels(bins)
#> [1] "(0,150]" "(150,300]" "(300,450]"
table(bins, useNA = "ifany")
#> bins
#> (0,150] (150,300] (300,450] <NA>
#> 2 6 6 2
# right = FALSE closes the bins on the left; labels name them.
size <- cut(
books$pages,
breaks = c(0, 150, 300, Inf), right = FALSE,
labels = c("short", "medium", "long")
)
table(size, useNA = "ifany")
#> size
#> short medium long <NA>
#> 2 6 6 2
table(size, books$genre)
#>
#> size classic fiction mystery non-fiction poetry science fiction
#> short 0 0 0 0 2 0
#> medium 1 3 1 1 0 0
#> long 1 0 1 1 0 3
# include.lowest keeps the left edge of the first bin; 88 pages is the minimum.
cut(88, breaks = c(88, 200, 300))
#> [1] <NA>
#> Levels: (88,200] (200,300]
cut(88, breaks = c(88, 200, 300), include.lowest = TRUE)
#> [1] [88,200]
#> Levels: [88,200] (200,300]
# breaks = n makes n equal-width bins over the range; dig.lab sets the labels.
table(cut(books$pages, breaks = 4, dig.lab = 4))
#>
#> (87.66,174] (174,260] (260,346] (346,432.3]
#> 2 5 3 4
table(cut(
books$pages,
breaks = quantile(books$pages, na.rm = TRUE), include.lowest = TRUE
))
#>
#> [88,202] (202,268] (268,344] (344,432]
#> 4 3 3 4
# tabulate() counts positive integers 1..nbins and is much faster than table().
tabulate(c(1, 2, 2, 5))
#> [1] 1 2 0 0 1
tabulate(c(1, 2, 2, 5), nbins = 3)
#> [1] 1 2 0
tabulate(as.integer(size), nbins = nlevels(size))
#> [1] 2 6 6
# Loans per ten-day block: 1 = days 0-9, 2 = days 10-19, and so on.
tabulate(loans$days[!is.na(loans$days)] %/% 10L + 1L)
#> [1] 5 15 26 6 1 1 1 1 2
# findInterval() gives the bin number without the factor.
findInterval(c(88, 214, 432), c(0, 150, 300, Inf))
#> [1] 1 2 3Source: R documentation: cut
Correlation and covariance
cor(x, y) returns Pearson’s correlation coefficient by default, method = "spearman" gives Spearman’s rho (the correlation of the ranks) and method = "kendall" gives Kendall’s tau, which is the tau-b variant when there are ties. The two rank-based measures are more robust and suit data that are not bivariate normal. cov() has the same arguments and returns the covariance, and cov2cor() rescales a covariance matrix into a correlation matrix. Given a matrix or a data frame as x, cor() returns the correlation matrix of its columns, and every column must be numeric: a data frame with a text column is an error, so subset the numeric columns first. The inputs must be numeric as is.numeric() sees them; logical values are also allowed, for historical compatibility.
The use argument decides what happens with missing values, and it is the argument most often left at its default by mistake: "everything" (the default) propagates NA; "all.obs" stops with an error if there is any; "complete.obs" drops every row with a missing value in any column (and is an error if no complete row is left); "na.or.complete" does the same but returns NA when no row is left; and "pairwise.complete.obs" uses, for each pair of columns, all the rows that are complete for that pair. The last choice can give a matrix that is not positive semi-definite, and the correlations in it are computed from different subsets of the rows. Abbreviations of the values of use and method are accepted. The loan lengths and the page counts have a correlation near zero, as the synthetic data intend: the example shows the values, not an effect.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loan_books <- merge(loans, books, by = "isbn")
# Pearson (the default) measures straight-line association.
cor(loan_books$pages, loan_books$days)
#> [1] NA
cor(loan_books$pages, loan_books$days, use = "complete.obs")
#> [1] 0.01723476
cor(loan_books$pages, loan_books$days, use = "pairwise.complete.obs")
#> [1] 0.01723476
cov(loan_books$pages, loan_books$days, use = "complete.obs")
#> [1] 30.03191
# Rank-based measures: Spearman and Kendall.
cor(
loan_books$pages, loan_books$days,
use = "complete.obs", method = "spearman"
)
#> [1] 0.07751257
cor(
loan_books$pages, loan_books$days,
use = "complete.obs", method = "kendall"
)
#> [1] 0.05721833
# On a data frame: the correlation matrix of the numeric columns.
num <- loan_books[c("days", "pages", "year")]
round(cor(num, use = "complete.obs"), 3)
#> days pages year
#> days 1.000 0.012 0.045
#> pages 0.012 1.000 -0.542
#> year 0.045 -0.542 1.000
round(cor(num, use = "pairwise.complete.obs"), 3)
#> days pages year
#> days 1.000 0.017 0.065
#> pages 0.017 1.000 -0.542
#> year 0.065 -0.542 1.000
round(cov(num, use = "complete.obs"), 1)
#> days pages year
#> days 286.5 20.9 39.8
#> pages 20.9 9940.3 -2815.8
#> year 39.8 -2815.8 2718.6
round(cov2cor(cov(num, use = "complete.obs")), 3)
#> days pages year
#> days 1.000 0.012 0.045
#> pages 0.012 1.000 -0.542
#> year 0.045 -0.542 1.000
# Rows used: complete.obs keeps only the rows with no NA at all.
sum(complete.cases(num))
#> [1] 45
sum(complete.cases(num[c("days", "pages")]))
#> [1] 48
nrow(num)
#> [1] 59
colSums(!is.na(num))
#> days pages year
#> 57 49 56The cells of a correlation matrix in the example differ between the two ways of dropping rows because the rows used differ: 45 rows have all three values, but 48 have both a length of loan and a page count. The messages of cor() are plain and are worth knowing; the second example collects them.
source("scenario.R")
message_of <- \(expr) {
tryCatch(
{
expr
"no condition"
},
warning = \(w) paste("warning:", conditionMessage(w)),
error = \(e) paste("error:", conditionMessage(e))
)
}
# Non-numeric input, a whole data frame with text columns, unequal lengths.
message_of(cor(books$title, books$pages))
#> [1] "error: 'x' must be numeric"
message_of(cor(books))
#> [1] "error: 'x' must be numeric"
message_of(cor(c(1, 2, 3), c(1, 2)))
#> [1] "error: incompatible dimensions"
# A constant has no variance, so the correlation is undefined: a warning, NA.
message_of(cor(c(1, 2, 3), c(5, 5, 5)))
#> [1] "warning: the standard deviation is zero"
suppressWarnings(cor(c(1, 2, 3), c(5, 5, 5)))
#> [1] NA
# A missing value: NA by default, an error with all.obs, a drop with use =.
cor(c(1, NA, 3, 4), c(2, 4, 7, 9))
#> [1] NA
message_of(cor(c(1, NA, 3, 4), c(2, 4, 7, 9), use = "all.obs"))
#> [1] "error: missing observations in cov/cor"
cor(c(1, NA, 3, 4), c(2, 4, 7, 9), use = "complete.obs")
#> [1] 0.9986254
message_of(cor(c(1, NA), c(NA, 4), use = "complete.obs"))
#> [1] "error: no complete element pairs"
# use is matched by its first letters, like method.
cor(1:5, c(2, 4, 5, 4, 5), method = "spear", use = "pair")
#> [1] 0.7378648A constant column has no variance and so no correlation; R returns NA with a warning, the standard deviation is zero. cor(books) fails with 'x' must be numeric because title, author, isbn and genre are text, and the message does not say which column is the culprit, so select the numeric ones with something like books[vapply(books, is.numeric, logical(1))]. The test for whether a correlation differs from zero is cor.test(), in the section on tests below.
Grouped summaries
A grouped summary applies a function to the pieces of a data frame defined by one or more grouping variables. Base R has four tools, which differ in what they return. aggregate() returns a data frame with one row per group and is the one to reach for first: its formula method (aggregate(days ~ genre, data, FUN)) omits the rows with missing values in any variable of the formula, and its default method takes a list of grouping vectors and leaves the missing-value decision to FUN. If FUN returns several numbers the result has a matrix column, which is handy and awkward (agg$days[, "mean"]). tapply(X, INDEX, FUN) returns an array indexed by the groups, a named vector for one grouping variable and a matrix for two, and gives NA for an empty combination. ave(x, g, FUN = mean) returns a vector the length of x, with the group value repeated on every row, which is the form wanted to centre or standardise within groups, to add a group total as a new column or to filter rows. by(data, INDICES, FUN) applies a function to the rows of a data frame, group by group, and prints in its own format. The apply family in general is described in The apply Family and Functionals, and split() plus vapply() (shown in the section on "htest" objects) is the building block under all of them.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# aggregate() with a formula: one row per group, rows with an NA are dropped.
aggregate(days ~ genre, data = loan_books, FUN = mean)
#> genre days
#> 1 classic 21.40000
#> 2 fiction 26.72727
#> 3 mystery 20.75000
#> 4 non-fiction 21.20000
#> 5 poetry 28.50000
#> 6 science fiction 26.13333
aggregate(cbind(days, pages) ~ genre, data = loan_books, FUN = median)
#> genre days pages
#> 1 classic 23.0 432
#> 2 fiction 22.0 214
#> 3 mystery 18.5 244
#> 4 non-fiction 21.5 175
#> 5 poetry 24.5 92
#> 6 science fiction 21.0 387
aggregate(days ~ overdue + genre, data = loan_books, FUN = length) |> head(4)
#> overdue genre days
#> 1 FALSE classic 2
#> 2 TRUE classic 3
#> 3 FALSE fiction 10
#> 4 TRUE fiction 12
# A function that returns several numbers gives a matrix column.
agg <- aggregate(
days ~ genre,
data = loan_books,
FUN = \(x) c(n = length(x), mean = mean(x), sd = sd(x))
)
agg
#> genre days.n days.mean days.sd
#> 1 classic 5.000000 21.400000 6.348228
#> 2 fiction 22.000000 26.727273 20.112024
#> 3 mystery 4.000000 20.750000 12.945398
#> 4 non-fiction 5.000000 21.200000 4.658326
#> 5 poetry 6.000000 28.500000 22.367387
#> 6 science fiction 15.000000 26.133333 17.983590
dim(agg)
#> [1] 6 2
agg$days[, "mean"]
#> [1] 21.40000 26.72727 20.75000 21.20000 28.50000 26.13333
# The non-formula form takes grouping vectors and drops nothing by itself, so
# FUN must cope with the NAs.
aggregate(
loan_books["days"],
by = list(genre = loan_books$genre),
FUN = mean, na.rm = TRUE
)
#> genre days
#> 1 classic 21.40000
#> 2 fiction 26.72727
#> 3 mystery 20.75000
#> 4 non-fiction 21.20000
#> 5 poetry 28.50000
#> 6 science fiction 26.13333
# tapply() returns an array indexed by the groups; by name, one value each.
tapply(loan_books$days, loan_books$genre, mean, na.rm = TRUE)
#> classic fiction mystery non-fiction poetry
#> 21.40000 26.72727 20.75000 21.20000 28.50000
#> science fiction
#> 26.13333
tapply(loan_books$days, list(loan_books$genre, loan_books$overdue), length)
#> FALSE TRUE
#> classic 2 3
#> fiction 10 12
#> mystery 2 2
#> non-fiction 2 3
#> poetry 2 4
#> science fiction 8 7
# ave() returns a vector the length of the data: the group value per row.
group_mean <- ave(
loan_books$days, loan_books$genre,
FUN = \(x) mean(x, na.rm = TRUE)
)
head(round(group_mean, 1), 3)
#> [1] 26.7 26.7 26.7
centred <- loan_books$days - group_mean # days relative to the genre average
round(tapply(centred, loan_books$genre, mean, na.rm = TRUE), 10)
#> classic fiction mystery non-fiction poetry
#> 0 0 0 0 0
#> science fiction
#> 0
# by() applies a function to the rows of a data frame, group by group.
by(
loan_books[c("days", "pages")], loan_books$genre,
\(d) round(colMeans(d, na.rm = TRUE), 1)
)
#> loan_books$genre: classic
#> days pages
#> 21.4 356.0
#> ------------------------------------------------------------
#> loan_books$genre: fiction
#> days pages
#> 26.7 213.2
#> ------------------------------------------------------------
#> loan_books$genre: mystery
#> days pages
#> 20.8 258.2
#> ------------------------------------------------------------
#> loan_books$genre: non-fiction
#> days pages
#> 21.2 175.0
#> ------------------------------------------------------------
#> loan_books$genre: poetry
#> days pages
#> 28.5 92.0
#> ------------------------------------------------------------
#> loan_books$genre: science fiction
#> days pages
#> 26.1 375.8Source: R documentation: aggregate
All of these are simple but base R only. The tidyverse packages (The tidyverse and dplyr) and data.table (data.table) have their own grouped-summary verbs, with a more regular syntax and, for large data, speed.
The htest object
Every classical test in the stats package returns a list of class "htest" (hypothesis test), with a fixed set of components that each test fills in as far as it applies. The help page of each test lists them in its Value section, and the print() method for the class, which is documented on the page of print.htest(), lays them out in the familiar block: a title (the method component), the data (data.name), the statistic and its degrees of freedom and p-value, the alternative hypothesis, the confidence interval and the sample estimates. The components of t.test() are statistic, parameter (the degrees of freedom), p.value, conf.int (with the confidence level stored in its "conf.level" attribute), estimate, null.value, stderr, alternative, method and data.name. stderr has been returned since R 3.6.0. Other tests have fewer components (shapiro.test() has no interval and no estimate) or more (chisq.test() also returns observed, expected, residuals and stdres). The fact that the object is a plain list is the useful property: $p.value is a number, $conf.int a vector of two, and str() shows the whole structure at once.
source("scenario.R")
# A test returns an object of class "htest": a plain list with a print method.
fit <- t.test(books$pages, mu = 250)
class(fit)
#> [1] "htest"
is.list(fit)
#> [1] TRUE
names(fit)
#> [1] "statistic" "parameter" "p.value" "conf.int" "estimate"
#> [6] "null.value" "stderr" "alternative" "method" "data.name"
fit
#>
#> One Sample t-test
#>
#> data: books$pages
#> t = 0.63111, df = 13, p-value = 0.5389
#> alternative hypothesis: true mean is not equal to 250
#> 95 percent confidence interval:
#> 205.8647 330.5639
#> sample estimates:
#> mean of x
#> 268.2143
#>
# Components, one by one.
fit$statistic
#> t
#> 0.6311118
fit$parameter
#> df
#> 13
fit$p.value
#> [1] 0.5389058
fit$estimate
#> mean of x
#> 268.2143
fit$conf.int
#> [1] 205.8647 330.5639
#> attr(,"conf.level")
#> [1] 0.95
attr(fit$conf.int, "conf.level")
#> [1] 0.95
fit$null.value
#> mean
#> 250
fit$stderr
#> [1] 28.86063
fit$alternative
#> [1] "two.sided"
fit$method
#> [1] "One Sample t-test"
fit$data.name
#> [1] "books$pages"
# str() shows the whole structure.
str(fit)
#> List of 10
#> $ statistic : Named num 0.631
#> ..- attr(*, "names")= chr "t"
#> $ parameter : Named num 13
#> ..- attr(*, "names")= chr "df"
#> $ p.value : num 0.539
#> $ conf.int : num [1:2] 206 331
#> ..- attr(*, "conf.level")= num 0.95
#> $ estimate : Named num 268
#> ..- attr(*, "names")= chr "mean of x"
#> $ null.value : Named num 250
#> ..- attr(*, "names")= chr "mean"
#> $ stderr : num 28.9
#> $ alternative: chr "two.sided"
#> $ method : chr "One Sample t-test"
#> $ data.name : chr "books$pages"
#> - attr(*, "class")= chr "htest"Source: R documentation: print.htest
That design is why code that calls a test and reports it can avoid parsing printed text: take what you need from the list. The second example does that. It pulls the p-value and an interval from one test, runs the same test on every genre with at least two page counts (a list of results, one per group, from split(), then vapply() to extract the p-values, because vapply() checks that every result is a single number), formats p-values for display with format.pval(), changes the number of digits of the printed block with print(x, digits = 3), and converts the result with broom. broom’s tidy() has a method for "htest" objects, and returns a one-row data frame whose columns are named estimate, statistic, p.value, parameter, conf.low, conf.high, method and alternative. It is what you want when the results of many tests must be stacked into one table; Statistical Models shows the same function on fitted models.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# A p-value or a confidence limit is just a number: pull it out by name.
fit <- t.test(books$pages, mu = 250)
fit$p.value < 0.05
#> [1] FALSE
round(fit$conf.int, 1)
#> [1] 205.9 330.6
#> attr(,"conf.level")
#> [1] 0.95
unname(fit$estimate)
#> [1] 268.2143
# Several tests at once: keep the results in a list and extract with vapply().
fits <- lapply(split(books$pages, books$genre), \(x) {
if (sum(!is.na(x)) < 2L) {
return(NULL)
}
t.test(x, mu = 250)
})
names(fits)
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
vapply(Filter(Negate(is.null), fits), \(f) f$p.value, numeric(1)) |> round(3)
#> classic fiction mystery non-fiction poetry
#> 0.396 0.255 0.575 0.904 0.016
#> science fiction
#> 0.043
# Rounding p-values for display: format.pval() and its eps argument.
format.pval(c(0.5, 0.012, 0.0000031, 1e-20), digits = 2)
#> [1] "0.500" "0.012" "3.1e-06" "< 2e-16"
format.pval(c(0.5, 0.012, 0.0000031, 1e-20), digits = 2, eps = 0.001)
#> [1] "0.500" "0.012" "<0.001" "<0.001"
# print.htest() prints with digits = getOption("digits"); change it per call.
print(fit, digits = 3)
#>
#> One Sample t-test
#>
#> data: books$pages
#> t = 0.6, df = 13, p-value = 0.5
#> alternative hypothesis: true mean is not equal to 250
#> 95 percent confidence interval:
#> 206 331
#> sample estimates:
#> mean of x
#> 268
#>
# broom turns an htest into a one-row data frame.
as.data.frame(broom::tidy(fit))
#> estimate statistic p.value parameter conf.low conf.high method
#> 1 268.2143 0.6311118 0.5389058 13 205.8647 330.5639 One Sample t-test
#> alternative
#> 1 two.sidedTwo details of format.pval(). Its eps argument is the threshold below which it prints <eps, and the default is .Machine$double.eps, which is why a very small p-value prints as < 2e-16; the 2.2e-16 of the printed test blocks has the same origin. A p-value of < 2.2e-16 is not zero and not a statement that the effect is huge: it only says that the number is smaller than the smallest difference the machine can distinguish from 1, and the honest summary of the result is the estimate and its interval. A second detail: since R 4.2.0 the result of ks.test() has class c("ks.test", "htest"), and the NEWS entry warns that code which tests the class with == must use inherits(x, "htest"); the same advice applies to every class vector in R.
Comparing means: t.test()
t.test() performs one-sample, two-sample and paired Student t tests. The signature of the default method is t.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, ...). With only x, it tests whether the mean is mu. With x and y, it tests whether the means differ by mu, and by default it uses the Welch (Satterthwaite) approximation, which does not assume that the two variances are equal, so the degrees of freedom are not whole numbers. var.equal = TRUE gives the classical pooled-variance Student test, which has whole degrees of freedom and slightly more power when the variances really are equal. paired = TRUE tests the mean of the differences x - y. Missing values are removed without a message, in pairs when paired = TRUE. alternative = "greater" is the alternative that x has the larger mean (in the one-sample case: that the mean is greater than mu), and a one-sided test returns a one-sided interval, with Inf as one limit. An Introduction to R first runs var.test() and then chooses between the Welch and the pooled test; the section on interpretation explains why the Welch test is the safer default and why a test of the variances first adds little.
One sample
The shelf asks whether the mean number of pages of the books is 250. The script prints the full block once, then takes single values from three variants: a one-sided alternative, a 99 percent interval, and the statistic computed by hand from its definition, which reproduces the printed t, the interval and the p-value.
source("scenario.R")
# H0: the mean length of a book is 250 pages. NAs are dropped by t.test().
t.test(books$pages, mu = 250)
#>
#> One Sample t-test
#>
#> data: books$pages
#> t = 0.63111, df = 13, p-value = 0.5389
#> alternative hypothesis: true mean is not equal to 250
#> 95 percent confidence interval:
#> 205.8647 330.5639
#> sample estimates:
#> mean of x
#> 268.2143
#>
# alternative: one-sided tests, and conf.level for the interval.
t.test(books$pages, mu = 250, alternative = "greater")
#>
#> One Sample t-test
#>
#> data: books$pages
#> t = 0.63111, df = 13, p-value = 0.2695
#> alternative hypothesis: true mean is greater than 250
#> 95 percent confidence interval:
#> 217.104 Inf
#> sample estimates:
#> mean of x
#> 268.2143
#>
t.test(books$pages, mu = 250, alternative = "less")$p.value
#> [1] 0.7305471
t.test(books$pages, mu = 250, conf.level = 0.99)$conf.int
#> [1] 181.2781 355.1505
#> attr(,"conf.level")
#> [1] 0.99
# By hand: t = (mean - mu) / (sd / sqrt(n)), and the interval from qt().
x <- books$pages[!is.na(books$pages)]
n <- length(x)
(mean(x) - 250) / (sd(x) / sqrt(n))
#> [1] 0.6311118
mean(x) + c(-1, 1) * qt(0.975, df = n - 1) * sd(x) / sqrt(n)
#> [1] 205.8647 330.5639
2 * pt(-abs((mean(x) - 250) / (sd(x) / sqrt(n))), df = n - 1)
#> [1] 0.5389058Source: R documentation: t.test
The statistic is (mean - mu) / (sd / sqrt(n)), the interval is the mean plus or minus the 0.975 quantile of the t distribution with n - 1 degrees of freedom times the standard error, and the two-sided p-value is twice the lower tail of the t distribution at minus the absolute statistic. The distributions are the subject of Distributions and Random Numbers. The p-value of 0.54 is large, and the interval is wide, from 206 to 331 pages: with 14 books the data are compatible with a mean well below and well above 250, which is a better summary of what 14 books can say than the statement that 250 was not rejected.
Two samples
Two independent samples are best shown with more data than the shelf has, so the next script uses penguins: the bill length of the Adelie and Chinstrap penguins (Palmer Station, Antarctica; the datasets help page gives the source). Welch is the default; var.equal = TRUE switches to the pooled test; the difference of the means is in the estimate component as two means, one per group, and the interval is the interval of the difference. The formula interface bill_len ~ species takes a response and a factor with exactly two levels; a factor with three levels is an error, with the message grouping factor must have exactly 2 levels, and subsetting the rows does not leave an unused level to confuse it, because the method calls factor() on the grouping variable.
# penguins is in the datasets package since R 4.5.0.
adelie <- penguins$bill_len[penguins$species == "Adelie"]
chinstrap <- penguins$bill_len[penguins$species == "Chinstrap"]
c(sum(!is.na(adelie)), sum(!is.na(chinstrap)))
#> [1] 151 68
# Two independent samples: Welch's test is the default (var.equal = FALSE).
welch <- t.test(adelie, chinstrap)
welch
#>
#> Welch Two Sample t-test
#>
#> data: adelie and chinstrap
#> t = -21.865, df = 106.97, p-value < 2.2e-16
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#> -10.952948 -9.131917
#> sample estimates:
#> mean of x mean of y
#> 38.79139 48.83382
#>
# var.equal = TRUE is the classical pooled-variance Student test.
pooled <- t.test(adelie, chinstrap, var.equal = TRUE)
pooled$method
#> [1] " Two Sample t-test"
pooled$parameter
#> df
#> 217
welch$parameter # Welch's degrees of freedom are not whole numbers
#> df
#> 106.967
c(welch = welch$p.value, pooled = pooled$p.value)
#> welch pooled
#> 2.900144e-41 2.011759e-62
# The estimate has one component per group; the difference is up to you.
welch$estimate
#> mean of x mean of y
#> 38.79139 48.83382
welch$estimate[[1L]] - welch$estimate[[2L]]
#> [1] -10.04243
welch$conf.int
#> [1] -10.952948 -9.131917
#> attr(,"conf.level")
#> [1] 0.95
# The formula interface: response ~ group, with a data frame. A group
# variable with more than two levels is an error, so subset first.
t.test(
bill_len ~ species,
data = penguins, subset = species != "Gentoo"
)$p.value
#> [1] 2.900144e-41
tryCatch(
t.test(bill_len ~ species, data = penguins),
error = \(e) conditionMessage(e)
)
#> [1] "grouping factor must have exactly 2 levels"Source: R documentation: t.test
The two p-values, 2.9e-41 and 2.0e-62, differ by many orders of magnitude and tell the same story. The Chinstrap bills are about 10 mm longer: the interval for Adelie minus Chinstrap runs from -10.95 to -9.13 mm (the printed estimates list the means in the order of the arguments). The effect size is what a reader needs, and the p-value is a footnote.
Paired samples
When each observation in one sample is matched with one in the other (the same patient under two drugs), the pairing carries information, and a paired test uses it. The classic example is sleep, the extra hours of sleep of ten patients under two drugs (Student, 1908; datasets help page). t.test(x, y, paired = TRUE) is the same test as the one-sample test of the differences. The formula interface for a paired test is Pair(x1, x2) ~ 1, introduced in R 4.0.0 when the formula method was extended to one-sample and paired tests. R 4.4.0 made passing paired = TRUE to the formula method an error, with the message cannot use 'paired' in formula method: an old call such as t.test(extra ~ group, data = sleep, paired = TRUE) used to be accepted and quietly paired the observations in their order of appearance, which is only right if the data are sorted consistently, so the change protects against a silent error.
# sleep: extra sleep of 10 patients under two drugs (datasets::sleep).
head(sleep, 3)
#> extra group ID
#> 1 0.7 1 1
#> 2 -1.6 1 2
#> 3 -0.2 1 3
x <- sleep$extra[sleep$group == 1]
y <- sleep$extra[sleep$group == 2]
# paired = TRUE tests the mean of the differences.
t.test(x, y, paired = TRUE)
#>
#> Paired t-test
#>
#> data: x and y
#> t = -4.0621, df = 9, p-value = 0.002833
#> alternative hypothesis: true mean difference is not equal to 0
#> 95 percent confidence interval:
#> -2.4598858 -0.7001142
#> sample estimates:
#> mean difference
#> -1.58
#>
# The same test, as a one-sample test of the differences.
t.test(y - x)$p.value
#> [1] 0.00283289
# Pair() (R 4.0.0) is the formula interface. paired = TRUE in a formula is an
# error since R 4.4.0.
wide <- data.frame(first = x, second = y)
t.test(Pair(first, second) ~ 1, data = wide)$p.value
#> [1] 0.00283289
tryCatch(
t.test(extra ~ group, data = sleep, paired = TRUE),
error = \(e) conditionMessage(e)
)
#> [1] "cannot use 'paired' in formula method"
# Ignoring the pairing gives a different answer: the Welch test is far less
# powerful here, because it throws away the information that the same person
# took both drugs.
t.test(x, y)$p.value
#> [1] 0.07939414Source: R documentation: t.test
The last line shows what the pairing buys. Ignoring it, the Welch test gives p = 0.079, not significant at 5 percent, whereas the paired test gives p = 0.0028. Pairing removes the variation between patients from the comparison. The converse mistake, treating two independent samples as paired, is as wrong, and the choice is a statement about how the data were collected, not about how the numbers look.
Rank-based tests: wilcox.test()
wilcox.test() performs the Wilcoxon signed rank test (one sample or paired) and the Wilcoxon rank sum test, also known as the Mann-Whitney test (two independent samples). They test a shift in location without assuming normality: the one-sample test asks whether the distribution is symmetric about mu, the two-sample test whether the two distributions differ by a location shift of mu. The help page cautions that the estimator for the difference that conf.int = TRUE returns is not the difference of the medians (a common misconception) but the median of the differences between a value of x and a value of y, the Hodges-Lehmann estimator, and that the one-sample estimate is the pseudomedian. The default method in R 4.6.1 is wilcox.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, exact = NULL, correct = TRUE, conf.int = FALSE, conf.level = 0.95, tol.root = 1e-4, digits.rank = Inf, ...). The formula interface, as for t.test(), takes a response and a two-level factor, or 1 on the right-hand side for the one-sample test, or Pair(x1, x2) ~ 1 for a paired one.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
two <- subset(loan_books, genre %in% c("fiction", "science fiction"))
fic <- two$days[two$genre == "fiction"]
sf <- two$days[two$genre == "science fiction"]
c(length(fic), length(sf))
#> [1] 23 15
sum(duplicated(c(fic, sf), incomparables = NA)) # tied loan lengths
#> [1] 9
# Two samples: the rank sum test (Mann-Whitney). The NAs are dropped.
wilcox.test(fic, sf)
#>
#> Wilcoxon rank sum exact test
#>
#> data: fic and sf
#> W = 168, p-value = 0.9331
#> alternative hypothesis: true location shift is not equal to 0
#>
wilcox.test(days ~ genre, data = two) # the same, with a formula
#>
#> Wilcoxon rank sum exact test
#>
#> data: days by genre
#> W = 168, p-value = 0.9331
#> alternative hypothesis: true location shift is not equal to 0
#>
# One sample: the signed rank test, against a median of loan_days (21).
wilcox.test(fic, mu = loan_days)
#>
#> Wilcoxon signed rank exact test
#>
#> data: fic
#> V = 149, p-value = 0.4675
#> alternative hypothesis: true location is not equal to 21
#>
# conf.int = TRUE adds the Hodges-Lehmann estimate and an interval.
w <- wilcox.test(fic, sf, conf.int = TRUE)
w$estimate
#> difference in location
#> 1
w$conf.int
#> [1] -9 9
#> attr(,"conf.level")
#> [1] 0.9558906Source: R documentation: wilcox.test
The two-sample test of the lengths of fiction and science fiction loans compares 22 and 15 loans with a known length (23 and 15 loans in all), and it gives an exact p-value although many loan lengths are the same whole number of days. This works since R 4.6.0.
Exact inference with ties in R 4.6
Before R 4.6.0, wilcox.test() could compute an exact p-value only when there were no ties and fewer than 50 observations. Real data nearly always have ties, because measurements are rounded: with ties, older versions printed a warning, cannot compute exact p-value with ties, and fell back to the normal approximation. The 4.6.0 NEWS entry reads that wilcox.test() can now perform exact (conditional) inference in case of ties. The help page of R 4.6.1 describes it: by default, if exact is not specified, an exact p-value is computed when the samples contain fewer than 50 finite values, and if there are ties the exact inference uses the conditional (permutation) distribution given the observed ranks, with an implementation of the Streitberg-Roehmel shift algorithm contributed by Torsten Hothorn. The message is therefore gone from wilcox.test(). It is not gone from R: cor.test() with method = "spearman" or "kendall" and the one-sample ks.test() still warn about ties (their texts are in the examples further down). An Introduction to R shows the old behaviour in its section 8.3, with the sample of ice-melting measurements of Rice (1995): the manual’s transcript ends with Warning message: Cannot compute exact p-value with ties in: wilcox.test(A, B) and a p-value of 0.007497. The container used for this page has only R 4.6.1, so the old warning cannot be reproduced here; what can be shown is that the new default gives a different p-value, that exact = FALSE gives the manual’s number back, and that the two samples of the manual have ties.
# The two samples of An Introduction to R, section 8.3 (latent heat of the
# fusion of ice, cal/gm; Rice 1995). Both samples contain ties.
a <- c(
79.98, 80.04, 80.02, 80.04, 80.03, 80.03, 80.04, 79.97,
80.05, 80.03, 80.02, 80.00, 80.02
)
b <- c(80.02, 79.94, 79.98, 79.97, 79.97, 80.03, 79.95, 79.97)
table(c(a, b))[table(c(a, b)) > 1]
#>
#> 79.97 79.98 80.02 80.03 80.04
#> 4 2 4 4 3
# The manual shows this call ending in a warning that an exact p-value cannot
# be computed with ties. R 4.6.1 computes it, conditionally on the ties.
wilcox.test(a, b)
#>
#> Wilcoxon rank sum exact test
#>
#> data: a and b
#> W = 89, p-value = 0.005435
#> alternative hypothesis: true location shift is not equal to 0
#>
# exact = FALSE reproduces the older behaviour: the normal approximation.
wilcox.test(a, b, exact = FALSE)
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: a and b
#> W = 89, p-value = 0.007497
#> alternative hypothesis: true location shift is not equal to 0
#>
# ks.test() with ties: exact in the two-sample case since R 4.2.0.
ks.test(a, b)
#>
#> Exact two-sample Kolmogorov-Smirnov test
#>
#> data: a and b
#> D = 0.59615, p-value = 0.03199
#> alternative hypothesis: two-sided
#>
# The tests that assume normality, for comparison.
t.test(a, b)$p.value
#> [1] 0.006939327
var.test(a, b)$p.value
#> [1] 0.393769The p-value changes from 0.0075 to 0.0054 when the exact conditional distribution replaces the normal approximation. It is the same data and the same test statistic, W = 89; what changed is how the p-value is found, and exact p-values are the reference that approximations try to match. The manual in this release still shows the pre-4.6.0 output, a reminder that tutorials, and manuals, lag behind the software.
What "conditional" means is easiest to see on a small example, where everything can be listed. Given the pooled ranks, with tied values sharing the average rank, the null hypothesis says that every choice of which observations belong to the first sample is equally likely. With 8 and 9 observations there are 24,310 such choices; counting those whose rank sum is as small as the observed one, or as large, gives the exact p-value. The script does that by brute force with combn() and gets the same number that wilcox.test() reports.
# Two small samples with ties: 2.2, 3.3, 4.4, 5.5 and 9.9 each occur more
# than once in the pooled data.
x <- c(1.1, 2.2, 2.2, 3.3, 4.4, 5.5, 5.5, 6.6)
y <- c(2.2, 3.3, 3.3, 4.4, 7.7, 8.8, 9.9, 9.9, 10.1)
table(c(x, y))[table(c(x, y)) > 1]
#>
#> 2.2 3.3 4.4 5.5 9.9
#> 3 3 2 2 2
# R 4.6.0 and later: exact inference conditional on the ties (the default for
# fewer than 50 observations, or exact = TRUE).
wilcox.test(x, y)
#>
#> Wilcoxon rank sum exact test
#>
#> data: x and y
#> W = 18.5, p-value = 0.09486
#> alternative hypothesis: true location shift is not equal to 0
#>
wilcox.test(x, y, exact = TRUE)$p.value
#> [1] 0.09485808
# exact = FALSE: the normal approximation, with a tie-corrected variance and,
# by default, a continuity correction.
wilcox.test(x, y, exact = FALSE)
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: x and y
#> W = 18.5, p-value = 0.09957
#> alternative hypothesis: true location shift is not equal to 0
#>
wilcox.test(x, y, exact = FALSE, correct = FALSE)$p.value
#> [1] 0.08999873
# What "conditional" means: keep the pooled ranks (ties get average ranks) and
# list every way of choosing which 8 of the 17 observations form the sample x.
r <- rank(c(x, y))
sums <- apply(combn(length(r), length(x)), 2L, \(i) sum(r[i]))
length(sums)
#> [1] 24310
observed <- sum(r[seq_along(x)])
lower <- mean(sums <= observed)
upper <- mean(sums >= observed)
min(1, 2 * min(lower, upper)) # two-sided: twice the smaller tail
#> [1] 0.09485808
wilcox.test(x, y)$p.value
#> [1] 0.09485808Source: R documentation: wilcox.test
Edgeworth-corrected p-values
When the normal approximation is used (exact = FALSE, or 50 or more observations) there is a second change in R 4.6.0: the NEWS entry says that wilcox.test() can optionally compute improved asymptotic p-values by including up to 3 correction terms of the Edgeworth series. The option is the existing argument correct: a logical value keeps its old meaning (TRUE, the default, is the continuity correction, and FALSE is none), and an integer from 0 to 3 gives the number of Edgeworth terms, where 0 is the same as TRUE. Anything else is an error, 'correct' must be an integer between 0 and 3. The help page states that the improved p-values are available when the normal approximation is used and there are no ties; with ties the integer is accepted and has no effect, as the second part of the script shows. R 4.6.1 fixed the p-value computation of wilcox.test(x, exact = FALSE, correct = k), which had lost a factor of dnorm(z), so p-values computed with an integer correct in R 4.6.0 should be recomputed.
The digits.rank argument, new in R 4.0.0 (NEWS), decides which values count as tied: if finite, the ranks are computed from signif(r, digits.rank). It exists because floating-point arithmetic splits ties that the data should have, a case of R FAQ 7.31: 0.4 - 0.3 and 0.3 - 0.2 are not exactly equal.
# Without ties, the normal approximation can be improved with Edgeworth terms:
# correct = k (an integer from 0 to 3) in place of TRUE or FALSE.
x <- c(1.1, 2.3, 2.7, 3.3, 4.4, 5.5, 5.9, 6.6)
y <- c(2.2, 3.1, 3.9, 4.8, 7.7, 8.8, 9.9, 10.5, 11.1)
p_exact <- wilcox.test(x, y)$p.value
p_normal <- vapply(
list(none = FALSE, continuity = TRUE, k0 = 0L, k1 = 1L, k2 = 2L, k3 = 3L),
\(k) wilcox.test(x, y, exact = FALSE, correct = k)$p.value,
numeric(1)
)
round(c(exact = p_exact, p_normal), 4)
#> exact none continuity k0 k1 k2 k3
#> 0.1139 0.1019 0.1124 0.1124 0.1139 0.1135 0.1141
# correct = TRUE and correct = 0 are the same thing: no Edgeworth terms, with
# the continuity correction. The method line says which was used.
wilcox.test(x, y, exact = FALSE, correct = 3)$method
#> [1] "Wilcoxon rank sum test with continuity correction"
tryCatch(
wilcox.test(x, y, exact = FALSE, correct = 4),
error = \(e) conditionMessage(e)
)
#> [1] "'correct' must be an integer between 0 and 3"
# With ties the Edgeworth terms are not available: every k gives the same
# continuity-corrected p-value.
x_tied <- c(1.1, 2.2, 2.2, 3.3, 4.4, 5.5, 5.5, 6.6)
y_tied <- c(2.2, 3.3, 3.3, 4.4, 7.7, 8.8, 9.9, 9.9, 10.1)
vapply(
0:3,
\(k) wilcox.test(x_tied, y_tied, exact = FALSE, correct = k)$p.value,
numeric(1)
) |> round(4)
#> [1] 0.0996 0.0996 0.0996 0.0996
# digits.rank decides which numbers count as tied. Floating-point noise can
# split a tie: these differences should all be 0.1 in absolute value.
a <- c(0.3, 0.4, 0.5, 0.9, 1.2, 0.35, 0.8, 0.7)
b <- c(0.2, 0.3, 0.4, 0.5, 0.6, 0.1, 0.9, 0.2)
d <- a - b
rank(abs(d))
#> [1] 2 4 2 6 8 5 2 7
rank(signif(abs(d), 7))
#> [1] 2.5 2.5 2.5 6.0 8.0 5.0 2.5 7.0
wilcox.test(d)$p.value
#> [1] 0.03125
wilcox.test(d, digits.rank = 7)$p.value
#> [1] 0.0390625
# Nothing left to test after the missing values are dropped; a bad mu.
tryCatch(wilcox.test(c(NA_real_, NA_real_)), error = \(e) conditionMessage(e))
#> [1] "not enough (non-missing) 'x' observations"
tryCatch(wilcox.test(1:5, 6:7, mu = "a"), error = \(e) conditionMessage(e))
#> [1] "'mu' must be a single number"Source: R documentation: wilcox.test
Counts: chisq.test() and fisher.test()
chisq.test() does two jobs. Given a vector of counts, it performs a goodness-of-fit test against the probabilities p (equal by default; rescale.p = TRUE rescales p if it does not sum to 1, otherwise that is an error). Given a matrix or two factors, it performs Pearson’s test of independence on the two-way table. For a 2 by 2 table, correct = TRUE (the default) applies Yates' continuity correction. The statistic is the sum of (observed - expected)^2 / expected, referred to a chi-squared distribution with (rows - 1) (columns - 1) degrees of freedom. The result holds the observed and expected counts, the Pearson residuals and the standardised residuals stdres, whose absolute values above about 2 point to the cells that depart from independence.
The approximation behind the p-value is poor when the expected counts are small, and R warns: Chi-squared approximation may be incorrect, when any expected count is below 5. R warns as soon as any expected count is below 5, which is the simple form of a common rule of thumb (Cochran’s rule: no expected count below 1 and at most a fifth of them below 5). There are two remedies. One is to simulate: simulate.p.value = TRUE computes a Monte Carlo p-value from B random tables with the same margins (2000 by default; the smallest possible p-value is 1/(B + 1)), which depends on the random seed, so call set.seed() first. The other is Fisher’s exact test, fisher.test(), which conditions on the margins and for a 2 by 2 table also returns the conditional maximum likelihood estimate of the odds ratio and its interval; for a larger table it is exact but can be slow. The shelf’s table of overdue loans by genre has expected counts as low as 1.8, so it triggers the warning. The script captures the warning text with withCallingHandlers() and still uses the result.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
by_genre <- table(
loan_books$genre, loan_books$overdue,
dnn = c("genre", "overdue")
)
# The test of independence on a 6 x 2 table. Some expected counts are below 5,
# so R warns; withCallingHandlers() keeps the warning text and the result.
res <- withCallingHandlers(
chisq.test(by_genre),
warning = \(w) {
cat("warning:", conditionMessage(w), "\n")
invokeRestart("muffleWarning")
}
)
#> warning: Chi-squared approximation may be incorrect
res
#>
#> Pearson's Chi-squared test
#>
#> data: by_genre
#> X-squared = 0.88335, df = 5, p-value = 0.9714
#>
round(res$expected, 2) # what the warning is about: cells below 5
#> overdue
#> genre FALSE TRUE
#> classic 2.28 2.72
#> fiction 10.04 11.96
#> mystery 1.82 2.18
#> non-fiction 2.28 2.72
#> poetry 2.74 3.26
#> science fiction 6.84 8.16
res$expected < 5
#> overdue
#> genre FALSE TRUE
#> classic TRUE TRUE
#> fiction FALSE FALSE
#> mystery TRUE TRUE
#> non-fiction TRUE TRUE
#> poetry TRUE TRUE
#> science fiction FALSE FALSE
round(res$stdres, 2) # standardised residuals: |value| > 2 is notable
#> overdue
#> genre FALSE TRUE
#> classic -0.26 0.26
#> fiction -0.02 0.02
#> mystery 0.18 -0.18
#> non-fiction -0.26 0.26
#> poetry -0.64 0.64
#> science fiction 0.70 -0.70
# Monte Carlo p-value: no approximation, but it depends on the random seed.
set.seed(282)
chisq.test(by_genre, simulate.p.value = TRUE, B = 10000)
#>
#> Pearson's Chi-squared test with simulated p-value (based on 10000
#> replicates)
#>
#> data: by_genre
#> X-squared = 0.88335, df = NA, p-value = 0.9733
#>
# Fisher's exact test on the same table (for r x c tables it can be slow).
fisher.test(by_genre)$p.value
#> [1] 0.9733191
# A 2 x 2 table: Yates' continuity correction is the default for chisq.test().
long <- table(
long_book = loan_books$pages > 300, overdue = loan_books$overdue
)
long
#> overdue
#> long_book FALSE TRUE
#> FALSE 12 17
#> TRUE 11 8
chisq.test(long)$p.value
#> [1] 0.409553
chisq.test(long, correct = FALSE)$p.value
#> [1] 0.2626761
fisher.test(long)
#>
#> Fisher's Exact Test for Count Data
#>
#> data: long
#> p-value = 0.3767
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#> 0.1347412 1.9247448
#> sample estimates:
#> odds ratio
#> 0.5206823
#>Source: R documentation: chisq.test
The three p-values for the 6 by 2 table, 0.971 from the approximation, 0.973 simulated and 0.973 from Fisher’s test, agree that these data give no evidence that the proportion of overdue loans depends on the genre. That is a statement about 57 loans, nearly all of fiction and science fiction, and not about libraries. For the 2 by 2 table of long and short books against overdue the corrected and uncorrected chi-squared tests give 0.41 and 0.26, and Fisher’s test 0.38; the continuity correction is conservative, and which one is quoted should not be decided by which one is smaller. summary() of a table or an xtabs object runs the same test of independence and prints the same warning, as the earlier example showed.
Normality: shapiro.test()
shapiro.test(x) performs the Shapiro-Wilk test of the null hypothesis that the sample comes from a normal distribution. It needs between 3 and 5000 non-missing values, and it stops with sample size must be between 3 and 5000 otherwise (and with all 'x' values are identical for a constant sample). The p-value is approximate: the help page says that the calculation is exact for n = 3, and uses approximations after that, said by Royston (1995) to be adequate for p-values below 0.1. The test is easy to run and easy to misuse. With few observations it has little power, so a large p-value says very little (the shelf’s 14 page counts give p = 0.83, and so would many non-normal samples of that size); with thousands of observations it rejects tiny departures from normality that do not matter for a t test (the 4000 draws from a t distribution with 10 degrees of freedom, which looks normal to the eye, reject at p = 2e-12). The question that matters is how far from normal the data are and how sensitive the method is to it, and a normal quantile plot (qqnorm() and qqline(), in Base Graphics and Devices) answers it better than a p-value. An Introduction to R runs shapiro.test() and ks.test() on the eruption data in its section 8.2, and notes that the distribution theory of the Kolmogorov-Smirnov test is not valid when the parameters of the normal distribution are estimated from the same sample.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# H0: the sample comes from a normal distribution. NAs are dropped.
shapiro.test(books$pages)
#>
#> Shapiro-Wilk normality test
#>
#> data: books$pages
#> W = 0.96678, p-value = 0.8309
#>
shapiro.test(loans$days) # loan lengths are right-skewed: strong evidence
#>
#> Shapiro-Wilk normality test
#>
#> data: loans$days
#> W = 0.81717, p-value = 5.256e-07
#>
shapiro.test(log(loans$days[loans$days > 0L]))$p.value
#> [1] 0.0001563376
# Bounds: 3 to 5000 non-missing values, and not all identical.
tryCatch(shapiro.test(c(1, 2)), error = \(e) conditionMessage(e))
#> [1] "sample size must be between 3 and 5000"
tryCatch(shapiro.test(rep(2, 10)), error = \(e) conditionMessage(e))
#> [1] "all 'x' values are identical"
set.seed(282)
tryCatch(shapiro.test(rnorm(5001)), error = \(e) conditionMessage(e))
#> [1] "sample size must be between 3 and 5000"
# With few data nothing is rejected; with many, tiny departures are.
set.seed(282)
shapiro.test(rnorm(30, mean = 270, sd = 100))$p.value
#> [1] 0.06891999
shapiro.test(rt(4000, df = 10))$p.value
#> [1] 2.353306e-12
# Penguin body mass: n = 342 and a small p-value, though the histogram is only
# mildly skewed. A normal quantile plot (qqnorm) shows how it departs.
shapiro.test(penguins$body_mass)$p.value
#> [1] 3.679039e-08Source: R documentation: shapiro.test
Proportions: prop.test() and binom.test()
prop.test(x, n, p = NULL, alternative, conf.level, correct = TRUE) tests proportions with a normal (chi-squared) approximation: given a number of successes and trials it tests one proportion against p, given vectors it tests that several proportions are equal, and given a two-column matrix of successes and failures it does the same. With one or two groups it applies Yates' continuity correction by default (correct = TRUE), and the confidence interval depends on the number of groups: for one group it is the interval for the proportion obtained by inverting the score test, for two groups it is an interval for the difference of the two proportions, both with a continuity correction where it applies, and for more groups there is none. binom.test(x, n, p = 0.5) is exact: it sums binomial probabilities, so it has no approximation to go wrong for small n, and its interval is the Clopper-Pearson one. For a large number of trials the two agree. Both need a count of successes, so a data frame with an overdue flag is summarised with sum() and sum(!is.na()) first; a missing flag is neither a success nor a failure, and counting it as a failure is a quiet error.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# One proportion: 32 of the 58 loans with a known outcome were overdue.
n_over <- sum(loans$overdue, na.rm = TRUE)
n_known <- sum(!is.na(loans$overdue))
c(n_over, n_known)
#> [1] 32 58
prop.test(n_over, n_known, p = 0.5)
#>
#> 1-sample proportions test with continuity correction
#>
#> data: n_over out of n_known, null probability 0.5
#> X-squared = 0.43103, df = 1, p-value = 0.5115
#> alternative hypothesis: true p is not equal to 0.5
#> 95 percent confidence interval:
#> 0.4162437 0.6803595
#> sample estimates:
#> p
#> 0.5517241
#>
prop.test(n_over, n_known, p = 0.5, correct = FALSE)$p.value
#> [1] 0.4307911
# binom.test() is exact: it sums binomial probabilities, no approximation.
binom.test(n_over, n_known, p = 0.5)
#>
#> Exact binomial test
#>
#> data: n_over and n_known
#> number of successes = 32, number of trials = 58, p-value = 0.5118
#> alternative hypothesis: true probability of success is not equal to 0.5
#> 95 percent confidence interval:
#> 0.4154233 0.6825624
#> sample estimates:
#> probability of success
#> 0.5517241
#>
# Two proportions: overdue rate of long and short books, from a 2 x 2 table.
long <- loan_books$pages > 300
ok <- !is.na(long) & !is.na(loan_books$overdue)
counts <- table(long[ok], loan_books$overdue[ok])
counts
#>
#> FALSE TRUE
#> FALSE 12 17
#> TRUE 11 8
prop.test(counts[, "TRUE"], rowSums(counts))
#>
#> 2-sample test for equality of proportions with continuity correction
#>
#> data: counts[, "TRUE"] out of rowSums(counts)
#> X-squared = 0.6801, df = 1, p-value = 0.4096
#> alternative hypothesis: two.sided
#> 95 percent confidence interval:
#> -0.1637395 0.4940480
#> sample estimates:
#> prop 1 prop 2
#> 0.5862069 0.4210526
#>
# The same with x = successes and n = trials per group; x and n as a matrix.
prop.test(c(17, 8), c(29, 19), correct = FALSE)$estimate
#> prop 1 prop 2
#> 0.5862069 0.4210526
prop.test(cbind(counts[, "TRUE"], counts[, "FALSE"]))$p.value
#> [1] 0.409553Source: R documentation: prop.test
The overdue rate is 32 out of 58 known loans, 55 percent, with a 95 percent interval from 42 to 68 percent: consistent with half, and with a majority. The exact and approximate p-values (0.512 and 0.512) agree to three decimals, and the two-sample comparison of long and short books has the same p-value, 0.41, as chisq.test() of the same table, because prop.test() for two groups is the chi-squared test of a 2 by 2 table.
Association, variances and distributions
cor.test(x, y) tests the hypothesis that a correlation is zero. For method = "pearson" it uses a t statistic with n - 2 degrees of freedom and returns a confidence interval from the Fisher z transformation. For Spearman’s rho and Kendall’s tau it uses exact p-values where it can and an approximation otherwise; with ties it cannot compute an exact p-value, and warns cannot compute exact p-value with ties. The formula interface is one-sided: ~ x + y. The loans of the shelf show a Pearson correlation of 0.017 between the number of pages and the length of a loan (p = 0.91, interval from -0.27 to 0.30): a correlation of this size cannot be told from zero with 48 pairs, and neither can a correlation of 0.25 be excluded.
var.test(x, y) compares the variances of two samples with an F test, and, like the F test in general, it is sensitive to non-normality, so treat it as a rough guide. Its formula interface is the same as that of t.test().
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# cor.test(): is the correlation different from zero? (pages and loan length)
cor.test(loan_books$pages, loan_books$days)
#>
#> Pearson's product-moment correlation
#>
#> data: loan_books$pages and loan_books$days
#> t = 0.11691, df = 46, p-value = 0.9074
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#> -0.2682134 0.2999008
#> sample estimates:
#> cor
#> 0.01723476
#>
# Rank correlation tests; ties make the exact p-value impossible.
res <- withCallingHandlers(
cor.test(loan_books$pages, loan_books$days, method = "spearman"),
warning = \(w) {
cat("warning:", conditionMessage(w), "\n")
invokeRestart("muffleWarning")
}
)
#> warning: cannot compute exact p-value with ties
res$estimate
#> rho
#> 0.07751257
res$p.value
#> [1] 0.6005178
cor.test(
loan_books$pages, loan_books$days,
method = "spearman", exact = FALSE
)$p.value
#> [1] 0.6005178
# The formula interface is one-sided: ~ x + y.
cor.test(~ pages + year, data = loan_books)$estimate
#> cor
#> -0.5416231
# var.test(): do two samples have the same variance? (an F test; sensitive to
# non-normality, so treat it as a rough guide only)
two <- subset(loan_books, genre %in% c("fiction", "science fiction"))
var.test(days ~ genre, data = two)
#>
#> F test to compare two variances
#>
#> data: days by genre
#> F = 1.2507, num df = 21, denom df = 14, p-value = 0.6779
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.4422365 3.2065290
#> sample estimates:
#> ratio of variances
#> 1.250716
#>Source: R documentation: cor.test
ks.test() is the Kolmogorov-Smirnov test, based on the largest distance between two cumulative distribution functions. The one-sample form compares a sample with a named distribution (the second argument is the name of the cumulative distribution function, here "pexp", with its parameters passed after it); the two-sample form compares two samples. Both assume a continuous distribution. R 4.2.0 made two changes (NEWS): the two-sample test gives exact p-values with ties (and Monte Carlo ones), and ks.test() has a formula interface, y ~ 1 for the one-sample and y ~ group for the two-sample test. The one-sample test still warns about ties, ties should not be present for the one-sample Kolmogorov-Smirnov test. kruskal.test() is the Kruskal-Wallis rank sum test, the rank-based extension of the Wilcoxon test to three or more groups; its input is a response and a grouping variable, or a list of samples. If all observations are in one group it stops with all observations are in the same group.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# ks.test(): one sample against a distribution (here exponential, with the
# rate we used to simulate), or two samples against each other.
set.seed(42)
x <- rexp(40, rate = 0.05)
ks.test(x, "pexp", rate = 0.05)
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: x
#> D = 0.12831, p-value = 0.4861
#> alternative hypothesis: two-sided
#>
class(ks.test(x, "pexp", rate = 0.05))
#> [1] "ks.test" "htest"
# The one-sample test assumes a continuous distribution: ties give a warning.
tied <- loans$days[!is.na(loans$days)]
tryCatch(
ks.test(tied, "pexp", rate = 1 / mean(tied)),
warning = \(w) conditionMessage(w)
)
#> [1] "ties should not be present for the one-sample Kolmogorov-Smirnov test"
# The two-sample test has been exact with ties since R 4.2.0, and takes a
# formula there.
two <- subset(loan_books, genre %in% c("fiction", "science fiction"))
ks.test(days ~ genre, data = two)
#>
#> Exact two-sample Kolmogorov-Smirnov test
#>
#> data: days by genre
#> D = 0.16364, p-value = 0.9167
#> alternative hypothesis: two-sided
#>
# kruskal.test(): the rank test for three or more groups.
kruskal.test(days ~ genre, data = loan_books)
#>
#> Kruskal-Wallis rank sum test
#>
#> data: days by genre
#> Kruskal-Wallis chi-squared = 0.39325, df = 5, p-value = 0.9955
#>
tryCatch(
kruskal.test(1:5, rep("a", 5)),
error = \(e) conditionMessage(e)
)
#> [1] "all observations are in the same group"Source: R documentation: ks.test
Several groups: aov(), TukeyHSD() and pairwise tests
When there are more than two groups, comparing them two at a time with t.test() multiplies the chances of a false positive, and the first step is a single overall test. aov(y ~ g, data) fits a one-way analysis of variance, and summary() of the fit prints the analysis of variance table with the F statistic and its p-value. The object has classes "aov" and "lm", because aov() is lm() underneath, and the extractor functions of the next page apply (Statistical Models). The test says that at least one mean differs, not which. For that, TukeyHSD() computes simultaneous confidence intervals for all the pairwise differences, with the 95 percent level applying to the whole family and its adjusted p-values labelled p adj. pairwise.t.test() compares each pair with t.test()-like statistics (with a pooled standard deviation by default, pool.sd = !paired) and adjusts the p-values with p.adjust(); its default is Holm’s method (p.adjust.method = p.adjust.methods, the first of which is "holm"). pairwise.wilcox.test() is the rank-based counterpart.
# One-way analysis of variance: does body mass differ among penguin species?
fit <- aov(body_mass ~ species, data = penguins)
summary(fit)
#> Df Sum Sq Mean Sq F value Pr(>F)
#> species 2 146864214 73432107 343.6 <2e-16 ***
#> Residuals 339 72443483 213698
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 2 observations deleted due to missingness
class(fit)
#> [1] "aov" "lm"
nobs(fit) # 2 of the 344 birds have no body mass
#> [1] 342
# aov() is lm() underneath: the coefficients are treatment contrasts.
coef(fit)
#> (Intercept) speciesChinstrap speciesGentoo
#> 3700.66225 32.42598 1375.35401
# Which pairs differ? TukeyHSD() controls the error over all three comparisons.
tuk <- TukeyHSD(fit)
tuk
#> Tukey multiple comparisons of means
#> 95% family-wise confidence level
#>
#> Fit: aov(formula = body_mass ~ species, data = penguins)
#>
#> $species
#> diff lwr upr p adj
#> Chinstrap-Adelie 32.42598 -126.5002 191.3522 0.8806666
#> Gentoo-Adelie 1375.35401 1243.1786 1507.5294 0.0000000
#> Gentoo-Chinstrap 1342.92802 1178.4810 1507.3750 0.0000000
#>
round(tuk$species, 3)
#> diff lwr upr p adj
#> Chinstrap-Adelie 32.426 -126.500 191.352 0.881
#> Gentoo-Adelie 1375.354 1243.179 1507.529 0.000
#> Gentoo-Chinstrap 1342.928 1178.481 1507.375 0.000
# pairwise.t.test() adjusts p-values with p.adjust(); holm is the default.
pairwise.t.test(penguins$body_mass, penguins$species)
#>
#> Pairwise comparisons using t tests with pooled SD
#>
#> data: penguins$body_mass and penguins$species
#>
#> Adelie Chinstrap
#> Chinstrap 0.63 -
#> Gentoo <2e-16 <2e-16
#>
#> P value adjustment method: holm
pairwise.t.test(
penguins$body_mass, penguins$species,
p.adjust.method = "bonferroni"
)
#>
#> Pairwise comparisons using t tests with pooled SD
#>
#> data: penguins$body_mass and penguins$species
#>
#> Adelie Chinstrap
#> Chinstrap 1 -
#> Gentoo <2e-16 <2e-16
#>
#> P value adjustment method: bonferroni
# The same question without the normality assumption: rank tests.
kruskal.test(body_mass ~ species, data = penguins)$p.value
#> [1] 5.609512e-48
pairwise.wilcox.test(penguins$body_mass, penguins$species)
#>
#> Pairwise comparisons using Wilcoxon rank sum test with continuity correction
#>
#> data: penguins$body_mass and penguins$species
#>
#> Adelie Chinstrap
#> Chinstrap 0.49 -
#> Gentoo <2e-16 <2e-16
#>
#> P value adjustment method: holmSource: R documentation: aov
Body mass differs by species (F = 343.6 on 2 and 339 degrees of freedom), and the pairwise comparisons say where: Gentoo is about 1.3 kg heavier than the other two, and Adelie and Chinstrap are not clearly different (adjusted p of 0.88 from Tukey and 0.63 from Holm). The "2 observations deleted due to missingness" in the table is a reminder of the two birds with no recorded mass. As before, when the groups are few and unbalanced the assumptions (similar variances, roughly normal residuals) matter more, and a rank test, kruskal.test(), is the cross-check.
Many tests: p.adjust()
If you run twenty tests of true null hypotheses at the 5 percent level, you expect one to come out significant by chance alone; if you then report only that one, the reported result is an artefact. p.adjust(p, method = p.adjust.methods, n = length(p)) corrects a vector of p-values for the number of tests. The method names are in the vector p.adjust.methods: "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr" (an alias of "BH") and "none". The first four control the family-wise error rate, the probability of at least one false positive. Bonferroni multiplies each p-value by the number of tests and is dominated by Holm’s method, which the help page says has no reason not to replace it, and which is valid under any dependence. "BH" (Benjamini and Hochberg) and "BY" control the false discovery rate, the expected share of false positives among the results declared significant, which is a weaker requirement and gives more power when many tests are run and some real effects are expected. n is the number of tests, and is larger than length(p) only if some p-values were not passed in. Missing p-values stay missing and are not counted unless n says so.
The result is in the same order as the input, so adjusted p-values can be put next to the raw ones. The example runs the t test of 250 pages in each genre (six p-values), and then a simulation of 20 tests on pure noise.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# Six one-sample t-tests (genre by genre, H0: mean pages = 250) give six
# p-values; only genres with at least two page counts can be tested.
by_genre <- split(books$pages, books$genre)
by_genre <- Filter(\(x) sum(!is.na(x)) >= 2L, by_genre)
p <- vapply(by_genre, \(x) t.test(x, mu = 250)$p.value, numeric(1))
round(p, 3)
#> classic fiction mystery non-fiction poetry
#> 0.396 0.255 0.575 0.904 0.016
#> science fiction
#> 0.043
# p.adjust() returns the adjusted p-values, in the same order.
p.adjust.methods
#> [1] "holm" "hochberg" "hommel" "bonferroni" "BH"
#> [6] "BY" "fdr" "none"
adjusted <- cbind(
raw = p,
bonferroni = p.adjust(p, method = "bonferroni"),
holm = p.adjust(p, method = "holm"),
BH = p.adjust(p, method = "BH")
)
round(adjusted, 3)
#> raw bonferroni holm BH
#> classic 0.396 1.000 1.000 0.594
#> fiction 0.255 1.000 1.000 0.510
#> mystery 0.575 1.000 1.000 0.689
#> non-fiction 0.904 1.000 1.000 0.904
#> poetry 0.016 0.097 0.097 0.097
#> science fiction 0.043 0.258 0.215 0.129
# After adjustment, how many are still below 0.05?
colSums(adjusted < 0.05)
#> raw bonferroni holm BH
#> 2 0 0 0
# n is the number of tests, if some p-values were not passed in.
p.adjust(p["poetry"], method = "holm", n = length(p))
#> poetry
#> 0.09668109
# Under a true null, 1 test in 20 gives p < 0.05 by chance alone.
set.seed(282)
null_p <- replicate(20, t.test(rnorm(15))$p.value)
sum(null_p < 0.05)
#> [1] 1
sum(p.adjust(null_p, "holm") < 0.05)
#> [1] 0
sum(p.adjust(null_p, "BH") < 0.05)
#> [1] 0
round(min(null_p), 3)
#> [1] 0.04
round(min(p.adjust(null_p, "holm")), 3)
#> [1] 0.805Source: R documentation: p.adjust
Two of the six genres look significant before adjustment (poetry p = 0.016 and science fiction p = 0.043), and none does after any of the three corrections (the smallest adjusted p-value is 0.097). The simulation shows the other side: with 20 tests of true nulls, one has p below 0.05 (the minimum is 0.04), and after the Holm correction its p-value is 0.8. Which method to use is decided before looking at the results. If one false positive is costly, for example a claim that goes into a report, use Holm; if you are screening many candidates and will follow up on the list, FDR control with BH is the usual choice. And in any case the number of tests is the number tried, not the number reported.
What a p-value does not say
The functions above produce numbers that are easy to misread. The limits, in order of how often they matter:
-
A p-value is the probability, if the null hypothesis and all the assumptions of the test are true, of data at least as extreme as the data observed. It is not the probability that the null hypothesis is true, not the probability that the result was produced by chance, and not the probability that the result will replicate. A p-value of 0.54 does not show that the mean is 250; it shows that these data are unsurprising if it were.
-
A small p-value says nothing about the size of the effect. With 300 penguins, a difference of a fraction of a millimetre is significant; with 14 books a difference of 100 pages is not. Report the estimate with its confidence interval, which has the units of the data, and use the p-value as a footnote.
t.test(),prop.test()andwilcox.test(conf.int = TRUE)return the interval. -
The tests assume that the observations are independent. The loans in the shelf data are not: there are 59 loans of 15 books by a handful of members, and loans of one book or by one member resemble each other. Treating every loan as an independent observation inflates the sample size and the degrees of freedom, so the test claims more precision than the data have; the example below shows how the answer changes when the unit of analysis becomes the book. The degrees of freedom fall from 32 to under 5, and the p-value moves from 0.93 to 0.61: pooling dependent observations misstates the uncertainty, and the resulting p-value can be wrong in either direction.
-
A non-significant result is not evidence of no effect. With 5 books per genre the power to detect a difference of 30 pages (with a standard deviation of 108) between two genres is about 6 percent, barely above the 5 percent that a test has when there is no difference at all.
power.t.test()works the problem in both directions: the sample size needed for a target power, or the power of a design. -
The data here are synthetic. The 16 books and 60 loans are generated to have realistic problems (missing values, spelling variants, a date that does not exist) and not realistic effects, so no result on this page, significant or not, is a finding about libraries or books.
-
Many tests, many chances: report how many were run, and correct for them (the previous section), or state that the analysis is exploratory.
-
Test assumptions are checked with diagnostics and subject knowledge, not with preliminary tests. A pre-test of normality or of equal variances before a t test makes the overall procedure a different procedure, with error rates nobody has derived; the safer defaults are the Welch test (no equal variance assumption) and a look at the data. The method still applies if the data are far from normal and the sample is large, because the t test is a test of a mean, and the mean of many observations is close to normal; it is then the outliers and the independence that need attention.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
back <- loans$return_date
back[is.na(back)] <- shelf_today
loans$days <- as.integer(back - loans$loan_date)
loans$days[loans$days < 0L] <- NA # a return before the loan is an error
loans$overdue <- loans$days > loan_days
loan_books <- merge(loans, books, by = "isbn")
# The loans of one book are not independent observations: they share the book.
# 59 loans, but only 15 distinct books.
c(loans = nrow(loan_books), books = length(unique(loan_books$isbn)))
#> loans books
#> 59 15
two <- subset(loan_books, genre %in% c("fiction", "science fiction"))
# Treating every loan as an independent observation...
t.test(days ~ genre, data = two)$p.value
#> [1] 0.9257102
t.test(days ~ genre, data = two)$parameter
#> df
#> 32.36671
# ...or averaging the loans of each book first, so that a book counts once.
per_book <- aggregate(days ~ isbn + genre, data = two, FUN = mean)
nrow(per_book)
#> [1] 7
res <- t.test(days ~ genre, data = per_book)
res$p.value
#> [1] 0.6091487
res$parameter
#> df
#> 4.713509
# Effect and uncertainty, not just the p-value: the difference in mean days
# and its 95% interval at the book level.
round(res$estimate, 1)
#> mean in group fiction mean in group science fiction
#> 29.6 25.8
round(res$conf.int, 1)
#> [1] -14.3 21.9
#> attr(,"conf.level")
#> [1] 0.95Source: R documentation: t.test
The unit of analysis changes the answer from p = 0.93 on the loans of the two genres to p = 0.61 on their 7 books, and the interval at the book level, from -14 to 22 days, is honest about how little 7 books can say.
Power and effect sizes
power.t.test() takes any four of n, delta, sd, sig.level and power, and solves for the fifth, in one of the three designs (type = "two.sample", "one.sample" or "paired", with alternative = "two.sided" or "one.sided"). n is the number in each group. The effect size used above, Cohen’s d, is the difference of the means in units of the pooled standard deviation; there is no function for it in base R, but the formula is short. Do the power calculation before collecting data: afterwards it only restates the p-value. The printed result is a "power.htest" object, printed by the second method on the page of print.htest(), and the NOTE: n is number in *each* group line is part of its output.
# Effect size and interval first: Adelie against Chinstrap bill length.
adelie <- penguins$bill_len[penguins$species == "Adelie"]
chinstrap <- penguins$bill_len[penguins$species == "Chinstrap"]
adelie <- adelie[!is.na(adelie)]
chinstrap <- chinstrap[!is.na(chinstrap)]
welch <- t.test(chinstrap, adelie)
round(welch$estimate[[1L]] - welch$estimate[[2L]], 2) # difference in mm
#> [1] 10.04
round(welch$conf.int, 2) # its 95% confidence interval
#> [1] 9.13 10.95
#> attr(,"conf.level")
#> [1] 0.95
# Cohen's d: the difference in units of the pooled standard deviation.
n1 <- length(adelie)
n2 <- length(chinstrap)
pooled_var <- ((n1 - 1) * var(adelie) + (n2 - 1) * var(chinstrap)) /
(n1 + n2 - 2)
pooled_sd <- sqrt(pooled_var)
round((mean(chinstrap) - mean(adelie)) / pooled_sd, 2)
#> [1] 3.48
# power.t.test(): give all but one of n, delta, sd, sig.level, power.
# How many books per group to detect a 30-page difference (sd 108) with 80%
# power? (108 is the standard deviation of the 14 known page counts.)
power.t.test(delta = 30, sd = 108, sig.level = 0.05, power = 0.8)
#>
#> Two-sample t test power calculation
#>
#> n = 204.4078
#> delta = 30
#> sd = 108
#> sig.level = 0.05
#> power = 0.8
#> alternative = two.sided
#>
#> NOTE: n is number in *each* group
#>
# The power of a comparison of two genres with 5 books each.
power.t.test(n = 5, delta = 30, sd = 108)$power
#> [1] 0.05815898
power.t.test(n = 5, delta = 100, sd = 108)$power
#> [1] 0.2520996
# One-sample and paired designs, and one-sided tests.
power.t.test(n = 14, delta = 50, sd = 108, type = "one.sample")$power
#> [1] 0.3612383
power.t.test(delta = 50, sd = 108, power = 0.8, type = "paired")$n |> ceiling()
#> [1] 39
power.t.test(delta = 50, sd = 108, power = 0.8, alternative = "one.sided")$n |>
ceiling()
#> [1] 59
# The smallest difference detectable with 80% power and 14 books per group.
power.t.test(n = 14, power = 0.8, sd = 108)$delta
#> [1] 118.8088Source: R documentation: power.t.test
Chinstrap bills are 10 mm longer than Adelie bills, a difference of 3.5 pooled standard deviations, which is a very large effect. To detect a difference of 30 pages with a standard deviation of 108 with 80 percent power one needs about 205 books in each genre, 14 per genre gives 80 percent power only for differences of about 119 pages, and 5 books per genre give a power of 6 percent for 30 pages and 25 percent for 100. These numbers are the honest way to say that a library catalogue of 16 books cannot answer the question.
Messages and what they mean
The messages below are what the functions of this page print. The script collects the text of the first warning or error of each call, so the wording is R 4.6.1’s. A message is not always a problem with the data. Chi-squared approximation may be incorrect is a property of the table, and cannot compute exact p-value with ties a property of the method; the others are mistakes in the call.
source("scenario.R")
# The first warning or error an expression raises, with its text.
message_of <- \(expr) {
tryCatch(
{
expr
"no condition"
},
warning = \(w) paste("warning:", conditionMessage(w)),
error = \(e) paste("error:", conditionMessage(e))
)
}
# The one- and two-sample t tests
message_of(t.test(5))
#> [1] "error: not enough 'x' observations"
message_of(t.test(c(NA_real_, 3)))
#> [1] "error: not enough 'x' observations"
message_of(t.test(c(4, 4, 4, 4)))
#> [1] "error: data are essentially constant"
message_of(t.test(1:5, 1:3, paired = TRUE))
#> [1] "error: not all arguments have the same length"
message_of(t.test(pages ~ genre, data = books))
#> [1] "error: grouping factor must have exactly 2 levels"
message_of(t.test("a"))
#> [1] "warning: argument is not numeric or logical: returning NA"
message_of(t.test(1:5, mu = c(1, 2)))
#> [1] "error: 'mu' must be a single number"
message_of(t.test(1:5, conf.level = 95))
#> [1] "error: 'conf.level' must be a single number between 0 and 1"
# Tables and counts
message_of(chisq.test(c(5, -1)))
#> [1] "error: all entries of 'x' must be nonnegative and finite"
message_of(chisq.test(matrix(c(5, 0, 5, 0), 2)))
#> [1] "warning: Chi-squared approximation may be incorrect"
message_of(chisq.test(c(3, 4), p = c(0.5, 0.2)))
#> [1] "error: probabilities must sum to 1."
message_of(chisq.test(table(books$genre)))
#> [1] "warning: Chi-squared approximation may be incorrect"
message_of(fisher.test(matrix(1:3, 1)))
#> [1] "error: 'x' must have at least 2 rows and columns"
message_of(table(1:3, 1:2))
#> [1] "error: all arguments must have the same length"
message_of(prop.test(11, 10))
#> [1] "error: elements of 'x' must not be greater than those of 'n'"
message_of(binom.test(2.5, 10))
#> [1] "error: 'x' must be nonnegative and integer"
message_of(binom.test(5, 10, p = 2))
#> [1] "error: 'p' must be a single number between 0 and 1"
# Correlation tests, ties and sample size
message_of(cor.test(1:2, 3:4))
#> [1] "error: not enough finite observations"
message_of(cor.test(c(1, 2, 2, 3), c(1, 3, 2, 4), method = "kendall"))
#> [1] "warning: cannot compute exact p-value with ties"
message_of(var.test(1, 1:5))
#> [1] "error: not enough 'x' observations"
message_of(ks.test(c(1, 2, 2, 3), "pnorm"))
#> [1] "warning: ties should not be present for the one-sample Kolmogorov-Smirnov test"
# Quantiles and adjusted p-values
message_of(quantile(1:5, 1.5))
#> [1] "error: 'probs' outside [0,1]"
message_of(quantile(1:5, type = 10))
#> [1] "error: 'type' must be an integer in 1..9"
message_of(p.adjust(c(0.01, 0.02), "nothing"))
#> [1] "error: 'arg' should be one of “holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”"
message_of(power.t.test(delta = 1))
#> [1] "error: exactly one of 'n', 'delta', 'sd', 'power', and 'sig.level' must be NULL"Source: R documentation: t.test
| Message | Meaning and cure |
|---|---|
|
|
|
|
|
the formula interface of |
|
R 4.4.0 and later; use |
|
|
|
|
|
a constant vector given to |
|
|
|
some expected counts are below 5; use |
|
|
|
|
|
|
|
|
|
|
|
|
Pitfalls
-
Forgetting that the default for missing values differs.
mean(),median()andrange()giveNA,quantile()stops,fivenum()silently drops,summary()counts,cor()givesNA, andt.test()andwilcox.test()drop silently. Know which function you call, and report the number of values used. -
Treating
sd()as the population standard deviation. It divides by n - 1; for a complete population multiply the variance by (n - 1) / n. -
Comparing quartiles from different programs. The nine quantile types differ;
fivenum(),quantile()andsummary()are three different answers, and a spreadsheet or Minitab may use a fourth. -
Using
table()on a number with many distinct values. Every distinct value is a level. Bin first withcut(). -
Reading a row-proportion table as a column-proportion table. The
marginargument ofproportions()is the one that sums to 1. -
Running
chisq.test()on percentages or on the raw data. It needs counts, or two factors that it tabulates itself. A table of percentages gives a meaningless statistic. -
The one-sided alternative chosen after seeing the data.
alternative = "greater"halves the p-value when the data go that way; decide before looking. -
Taking
exact = TRUEfor granted for large samples. The help page ofwilcox.test()warns that withexact = TRUEit can take a very long time to run if one sample is large (several thousands or more). The default switches to the normal approximation at 50 observations for that reason. -
Calling
wilcox.test()a test of medians. It tests a location shift under the assumption that the distributions have the same shape; with different shapes it also detects differences that are not in the medians. -
Rounding p-values to zero.
p-value < 2.2e-16is a bound, not a zero. Quote the estimate and the interval. -
Comparing floating-point numbers that should be ties. See R FAQ 7.31,
digits.rankinwilcox.test()andround()before ranking.
Version notes
The entries come from the NEWS files of the corresponding releases (R NEWS). Line numbers are not quoted; each statement is a NEWS entry or a statement of the help page of R 4.6.1.
| Version | Change |
|---|---|
R 3.4.0 |
|
R 3.6.0 |
|
R 4.0.0 |
|
R 4.0.3 |
|
R 4.1.0 |
The names of the result of |
R 4.2.0 |
|
R 4.4.0 |
Passing |
R 4.5.0 |
|
R 4.5.2 |
|
R 4.6.0 |
|
R 4.6.1 |
The p-value of |
Where older tutorials differ
Tutorials that predate 2020 or so, and the manual itself in places, show behaviour that has changed. None of the lines below is recommended code.
| Older material | Current behaviour | Why | Since |
|---|---|---|---|
|
an exact conditional p-value for fewer than 50 observations; |
the exact test with ties is implemented |
R 4.6.0 |
|
|
a more informative default |
R 4.6.0 |
|
|
the names were unfortunate, and the new functions accept named margins |
R 4.0.0 |
|
an error; |
the formula method silently assumed that the rows were already paired |
R 4.4.0 |
|
|
the argument says what it does |
R 3.4.0 |
|
exact two-sample p-value |
R computes it |
R 4.2.0 |
names of |
|
a global option should not change names |
R 4.1.0 |
|
|
|
always |
|
|
|
always |
|
look at the data ( |
pre-tests change the error rates of what follows |
always |
|
|
|
always |
Where to go next
-
Statistical Models for
lm()andglm(), the next step after a comparison of means: the formula, the summary, the diagnostics and the samebroom::tidy()on fitted models. -
Distributions and Random Numbers for the
d,p,qandrfunctions behind every p-value on this page,set.seed()andsample(). -
Base Graphics and Devices for
hist(),boxplot(),qqnorm()and the other pictures that should accompany these numbers. -
Data Frames and Factors for the objects that
summary(),aggregate(),table()andxtabs()work on. -
Missing and Special Values for
NA,NaN,Infand what each summary does with them. -
Formulas for the formula language used by
t.test(),aggregate(),xtabs()andkruskal.test(). -
The apply Family and Functionals and Vectorisation and Iteration for
vapply(),tapply()and running a test over many groups. -
Dates and Times for the dates behind the loan lengths used above.
-
The tidyverse and dplyr for
group_by()andsummarise(), the tidyverse way of the grouped summaries.
References
-
An Introduction to R: 2.5 Missing values
-
An Introduction to R: 8.2 Examining the distribution of a set of data
-
An Introduction to R: 8.3 One- and two-sample tests
-
R FAQ 7.31: Why doesn’t R think these numbers are equal?
-
R NEWS: Changes in R 4.0.0 to 4.6.1
-
R 3.x NEWS: Changes in R 3.4.0 and 3.6.0
-
R-devel NEWS: the upcoming
cumvar()