Before the model, understand the sample
Before fitting a regression model, first understand the data.
Descriptive tables help us see sample size, missing values, category
frequencies, unusual values, skewed variables, and differences between
groups.
This tutorial uses real sample datasets supplied with R or
common R packages. No participant-level dataset is invented in
this document.
We use:
mtcars for categorical and continuous summaries;airquality for missing-data examples;ToothGrowth for grouped continuous summaries;Titanic for multi-way categorical tables;PlantGrowth for a compact grouped example;iris for simple categorical and continuous
summaries.
Setup
library(dplyr)
library(tidyr)
library(forcats)
library(knitr)
options(dplyr.summarise.inform = FALSE)
Install optional packages once from the R console if needed.
install.packages(c("dplyr", "tidyr", "forcats", "knitr", "gtsummary"))
Inspect the sample data first
We begin with mtcars, a built-in R dataset containing
measurements for 32 automobiles.
data(mtcars)
glimpse(mtcars)
## Rows: 32
## Columns: 11
## $ mpg <dbl> 21.0, 21.0, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2, 17.8,…
## $ cyl <dbl> 6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 8,…
## $ disp <dbl> 160.0, 160.0, 108.0, 258.0, 360.0, 225.0, 360.0, 146.7, 140.8, 16…
## $ hp <dbl> 110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, 180, 180…
## $ drat <dbl> 3.90, 3.90, 3.85, 3.08, 3.15, 2.76, 3.21, 3.69, 3.92, 3.92, 3.92,…
## $ wt <dbl> 2.620, 2.875, 2.320, 3.215, 3.440, 3.460, 3.570, 3.190, 3.150, 3.…
## $ qsec <dbl> 16.46, 17.02, 18.61, 19.44, 17.02, 20.22, 15.84, 20.00, 22.90, 18…
## $ vs <dbl> 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0,…
## $ am <dbl> 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,…
## $ gear <dbl> 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 4, 4, 4, 3, 3,…
## $ carb <dbl> 4, 4, 1, 1, 2, 1, 4, 2, 2, 4, 4, 3, 3, 3, 4, 4, 4, 1, 2, 1, 1, 2,…
dim(mtcars)
## [1] 32 11
summary(mtcars)
## mpg cyl disp hp
## Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
## 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
## Median :19.20 Median :6.000 Median :196.3 Median :123.0
## Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
## 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
## Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
## drat wt qsec vs
## Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
## 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
## Median :3.695 Median :3.325 Median :17.71 Median :0.0000
## Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
## 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
## Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
## am gear carb
## Min. :0.0000 Min. :3.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
## Median :0.0000 Median :4.000 Median :2.000
## Mean :0.4062 Mean :3.688 Mean :2.812
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :1.0000 Max. :5.000 Max. :8.000
colSums(is.na(mtcars))
## mpg cyl disp hp drat wt qsec vs am gear carb
## 0 0 0 0 0 0 0 0 0 0 0
Some variables stored as numbers actually represent categories. For
descriptive tables, it is helpful to create labelled factors.
mtcars2 <- mtcars |>
mutate(
transmission = factor(am, levels = c(0, 1), labels = c("Automatic", "Manual")),
cylinders = factor(cyl),
gears = factor(gear),
engine = factor(vs, levels = c(0, 1), labels = c("V-shaped", "Straight"))
)
glimpse(mtcars2)
## Rows: 32
## Columns: 15
## $ mpg <dbl> 21.0, 21.0, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.…
## $ cyl <dbl> 6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, …
## $ disp <dbl> 160.0, 160.0, 108.0, 258.0, 360.0, 225.0, 360.0, 146.7, 1…
## $ hp <dbl> 110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, …
## $ drat <dbl> 3.90, 3.90, 3.85, 3.08, 3.15, 2.76, 3.21, 3.69, 3.92, 3.9…
## $ wt <dbl> 2.620, 2.875, 2.320, 3.215, 3.440, 3.460, 3.570, 3.190, 3…
## $ qsec <dbl> 16.46, 17.02, 18.61, 19.44, 17.02, 20.22, 15.84, 20.00, 2…
## $ vs <dbl> 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, …
## $ am <dbl> 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, …
## $ gear <dbl> 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 4, 4, …
## $ carb <dbl> 4, 4, 1, 1, 2, 1, 4, 2, 2, 4, 4, 3, 3, 3, 4, 4, 4, 1, 2, …
## $ transmission <fct> Manual, Manual, Manual, Automatic, Automatic, Automatic, …
## $ cylinders <fct> 6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, …
## $ gears <fct> 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 4, 4, …
## $ engine <fct> V-shaped, V-shaped, Straight, Straight, V-shaped, Straigh…
One categorical variable: counts first, percentages second
The built-in iris dataset contains 150 flowers from
three species.
species_counts <- table(iris$Species)
species_counts
##
## setosa versicolor virginica
## 50 50 50
Relative frequencies are obtained with prop.table().
prop.table(species_counts)
##
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
100 * prop.table(species_counts)
##
## setosa versicolor virginica
## 33.33333 33.33333 33.33333
A tidy frequency table is often easier to export or present.
species_frequency <- as.data.frame(species_counts) |>
rename(species = Var1, count = Freq) |>
mutate(percent = count / sum(count) * 100)
species_frequency
## species count percent
## 1 setosa 50 33.33333
## 2 versicolor 50 33.33333
## 3 virginica 50 33.33333
species_frequency |>
mutate(percent_display = sprintf("%.1f%%", percent)) |>
select(species, count, percent_display) |>
kable(col.names = c("Species", "n", "%"),
caption = "Distribution of iris species")
Distribution of iris species
| setosa |
50 |
33.3% |
| versicolor |
50 |
33.3% |
| virginica |
50 |
33.3% |
The denominator is essential. A percentage becomes meaningful only
when we know what observations contribute to it.
Missing values and denominators
The built-in airquality dataset contains genuine missing
values in Ozone and Solar.R.
data(airquality)
colSums(is.na(airquality))
## Ozone Solar.R Wind Temp Month Day
## 37 7 0 0 0 0
Compare the total number of rows with the number of observed ozone
measurements.
nrow(airquality)
## [1] 153
sum(!is.na(airquality$Ozone))
## [1] 116
sum(is.na(airquality$Ozone))
## [1] 37
For a categorical demonstration, convert month to a labelled
factor.
airquality2 <- airquality |>
mutate(
Month = factor(
Month,
levels = 5:9,
labels = c("May", "June", "July", "August", "September")
),
ozone_available = factor(
if_else(is.na(Ozone), "Missing", "Observed"),
levels = c("Observed", "Missing")
)
)
table(airquality2$ozone_available)
##
## Observed Missing
## 116 37
prop.table(table(airquality2$ozone_available)) * 100
##
## Observed Missing
## 75.81699 24.18301
This makes the denominator explicit: all rows are included, and
missing ozone measurements form their own category.
Zero-count factor levels are different from missing values
A factor may contain a valid level that is not observed.
gear_extended <- factor(
mtcars$gear,
levels = c(3, 4, 5, 6),
labels = c("3 gears", "4 gears", "5 gears", "6 gears")
)
table(gear_extended, useNA = "ifany")
## gear_extended
## 3 gears 4 gears 5 gears 6 gears
## 15 12 5 0
The zero for 6 gears means that the category was defined
but no car in mtcars has that value. It is not the same as
NA.
fct_drop(gear_extended) |> table()
##
## 3 gears 4 gears 5 gears
## 15 12 5
Ordering categories intentionally
Categorical variables should be ordered in a way that helps
interpretation.
mtcars2 |>
count(cylinders, sort = TRUE)
## cylinders n
## 1 8 14
## 2 4 11
## 3 6 7
Frequency ordering can also be created with
fct_infreq().
mtcars2 |>
mutate(cylinders_frequency = fct_infreq(cylinders)) |>
count(cylinders_frequency, sort = FALSE)
## cylinders_frequency n
## 1 8 14
## 2 4 11
## 3 6 7
Frequency order can help exploratory displays, but substantive or
ordinal order should normally be preserved when it has scientific
meaning.
Two categorical variables: choose the denominator
Use transmission and number of cylinders from
mtcars.
transmission_cylinders <- table(
Transmission = mtcars2$transmission,
Cylinders = mtcars2$cylinders
)
transmission_cylinders
## Cylinders
## Transmission 4 6 8
## Automatic 3 4 12
## Manual 8 3 2
Row percentages
round(prop.table(transmission_cylinders, margin = 1) * 100, 1)
## Cylinders
## Transmission 4 6 8
## Automatic 15.8 21.1 63.2
## Manual 61.5 23.1 15.4
Each transmission row sums to 100%. This answers: among cars with a
particular transmission type, how are cylinder categories
distributed?
Column percentages
round(prop.table(transmission_cylinders, margin = 2) * 100, 1)
## Cylinders
## Transmission 4 6 8
## Automatic 27.3 57.1 85.7
## Manual 72.7 42.9 14.3
Each cylinder column sums to 100%. This answers: within a cylinder
category, what percentage of cars are automatic or manual?
Overall percentages
round(prop.table(transmission_cylinders) * 100, 1)
## Cylinders
## Transmission 4 6 8
## Automatic 9.4 12.5 37.5
## Manual 25.0 9.4 6.2
All cells together sum to 100%.
Counts and row percentages together
transmission_cylinders_tidy <- as.data.frame(transmission_cylinders) |>
as_tibble() |>
rename(count = Freq) |>
group_by(Transmission) |>
mutate(
row_total = sum(count),
row_percent = count / row_total * 100,
display = sprintf("%d (%.1f%%)", count, row_percent)
) |>
ungroup()
transmission_cylinders_tidy
## # A tibble: 6 × 6
## Transmission Cylinders count row_total row_percent display
## <fct> <fct> <int> <int> <dbl> <chr>
## 1 Automatic 4 3 19 15.8 3 (15.8%)
## 2 Manual 4 8 13 61.5 8 (61.5%)
## 3 Automatic 6 4 19 21.1 4 (21.1%)
## 4 Manual 6 3 13 23.1 3 (23.1%)
## 5 Automatic 8 12 19 63.2 12 (63.2%)
## 6 Manual 8 2 13 15.4 2 (15.4%)
transmission_cylinders_tidy |>
select(Transmission, Cylinders, display) |>
pivot_wider(names_from = Cylinders, values_from = display) |>
kable(caption = "Cylinder distribution within transmission type, n (row %)")
Cylinder distribution within transmission type, n (row
%)
| Automatic |
3 (15.8%) |
4 (21.1%) |
12 (63.2%) |
| Manual |
8 (61.5%) |
3 (23.1%) |
2 (15.4%) |
Totals
addmargins(transmission_cylinders)
## Cylinders
## Transmission 4 6 8 Sum
## Automatic 3 4 12 19
## Manual 8 3 2 13
## Sum 11 7 14 32
A third variable can change the story
The built-in Titanic dataset is a four-dimensional
contingency table containing class, sex, age group, survival and
frequencies.
Titanic
## , , Age = Child, Survived = No
##
## Sex
## Class Male Female
## 1st 0 0
## 2nd 0 0
## 3rd 35 17
## Crew 0 0
##
## , , Age = Adult, Survived = No
##
## Sex
## Class Male Female
## 1st 118 4
## 2nd 154 13
## 3rd 387 89
## Crew 670 3
##
## , , Age = Child, Survived = Yes
##
## Sex
## Class Male Female
## 1st 5 1
## 2nd 11 13
## 3rd 13 14
## Crew 0 0
##
## , , Age = Adult, Survived = Yes
##
## Sex
## Class Male Female
## 1st 57 140
## 2nd 14 80
## 3rd 75 76
## Crew 192 20
Collapse over passenger class and retain sex, age and survival.
titanic_three_way <- margin.table(
Titanic,
margin = c("Sex", "Age", "Survived")
)
titanic_three_way
## , , Survived = No
##
## Age
## Sex Child Adult
## Male 35 1329
## Female 17 109
##
## , , Survived = Yes
##
## Age
## Sex Child Adult
## Male 29 338
## Female 28 316
Separate two-way tables by age are often easier to understand.
for (age_group in dimnames(titanic_three_way)$Age) {
cat("\nAge stratum:", age_group, "\n")
print(titanic_three_way[, age_group, ])
print(round(prop.table(titanic_three_way[, age_group, ], margin = 1) * 100, 1))
}
##
## Age stratum: Child
## Survived
## Sex No Yes
## Male 35 29
## Female 17 28
## Survived
## Sex No Yes
## Male 54.7 45.3
## Female 37.8 62.2
##
## Age stratum: Adult
## Survived
## Sex No Yes
## Male 1329 338
## Female 109 316
## Survived
## Sex No Yes
## Male 79.7 20.3
## Female 25.6 74.4
Stratification can reveal patterns hidden in an overall table, but it
can also produce sparse cells. More dimensions do not automatically
create a better descriptive table.
Continuous variables: location, spread and availability
Use mpg from mtcars.
mpg <- mtcars$mpg
length(mpg)
## [1] 32
sum(!is.na(mpg))
## [1] 32
sum(is.na(mpg))
## [1] 0
mean(mpg)
## [1] 20.09062
median(mpg)
## [1] 19.2
sd(mpg)
## [1] 6.026948
var(mpg)
## [1] 36.3241
min(mpg)
## [1] 10.4
max(mpg)
## [1] 33.9
range(mpg)
## [1] 10.4 33.9
quantile(mpg, probs = c(0.10, 0.25, 0.50, 0.75, 0.90))
## 10% 25% 50% 75% 90%
## 14.340 15.425 19.200 22.800 30.090
summary(mpg)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.40 15.43 19.20 20.09 22.80 33.90
The mean describes the arithmetic centre, while the median describes
the middle ordered observation. SD describes variation around the mean.
IQR describes the width of the middle 50% of observations.
Skewness and influential values
The airquality$Ozone variable is useful for showing why
the mean and median may differ.
c(
observed = sum(!is.na(airquality$Ozone)),
mean = mean(airquality$Ozone, na.rm = TRUE),
median = median(airquality$Ozone, na.rm = TRUE),
sd = sd(airquality$Ozone, na.rm = TRUE),
iqr = IQR(airquality$Ozone, na.rm = TRUE)
)
## observed mean median sd iqr
## 116.00000 42.12931 31.50000 32.98788 45.25000
hist(
airquality$Ozone,
breaks = 12,
main = "Distribution of ozone measurements",
xlab = "Ozone"
)
boxplot(
airquality$Ozone,
horizontal = TRUE,
main = "Ozone boxplot",
xlab = "Ozone"
)
Mean and SD may be suitable for roughly symmetric variables. Median
and IQR are often informative when a distribution is skewed or contains
influential extreme values.
A reusable continuous-summary function
continuous_summary <- function(x) {
observed <- x[!is.na(x)]
if (length(observed) == 0) {
return(tibble(
n_total = length(x), n_observed = 0, n_missing = length(x),
mean = NA_real_, sd = NA_real_, median = NA_real_,
q1 = NA_real_, q3 = NA_real_, min = NA_real_, max = NA_real_
))
}
tibble(
n_total = length(x),
n_observed = length(observed),
n_missing = sum(is.na(x)),
mean = mean(observed),
sd = sd(observed),
median = median(observed),
q1 = unname(quantile(observed, 0.25)),
q3 = unname(quantile(observed, 0.75)),
min = min(observed),
max = max(observed)
)
}
continuous_summary(airquality$Ozone)
## # A tibble: 1 × 10
## n_total n_observed n_missing mean sd median q1 q3 min max
## <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int>
## 1 153 116 37 42.1 33.0 31.5 18 63.2 1 168
continuous_summary(mtcars$mpg)
## # A tibble: 1 × 10
## n_total n_observed n_missing mean sd median q1 q3 min max
## <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 32 32 0 20.1 6.03 19.2 15.4 22.8 10.4 33.9
Summarise several continuous variables without creating fake
observations.
continuous_variables <- c("mpg", "disp", "hp", "wt")
continuous_table <- lapply(
continuous_variables,
function(variable) {
continuous_summary(mtcars[[variable]]) |>
mutate(variable = variable, .before = 1)
}
) |>
bind_rows()
continuous_table
## # A tibble: 4 × 11
## variable n_total n_observed n_missing mean sd median q1 q3
## <chr> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mpg 32 32 0 20.1 6.03 19.2 15.4 22.8
## 2 disp 32 32 0 231. 124. 196. 121. 326
## 3 hp 32 32 0 147. 68.6 123 96.5 180
## 4 wt 32 32 0 3.22 0.978 3.32 2.58 3.61
## # ℹ 2 more variables: min <dbl>, max <dbl>
continuous_display <- continuous_table |>
mutate(
`Mean (SD)` = sprintf("%.1f (%.1f)", mean, sd),
`Median (IQR)` = sprintf("%.1f (%.1f-%.1f)", median, q1, q3),
`Min-Max` = sprintf("%.1f-%.1f", min, max)
) |>
select(
Variable = variable,
`Observed n` = n_observed,
`Missing n` = n_missing,
`Mean (SD)`,
`Median (IQR)`,
`Min-Max`
)
continuous_display |>
kable(caption = "Continuous-variable descriptive summary for mtcars")
Continuous-variable descriptive summary for mtcars
| mpg |
32 |
0 |
20.1 (6.0) |
19.2 (15.4-22.8) |
10.4-33.9 |
| disp |
32 |
0 |
230.7 (123.9) |
196.3 (120.8-326.0) |
71.1-472.0 |
| hp |
32 |
0 |
146.7 (68.6) |
123.0 (96.5-180.0) |
52.0-335.0 |
| wt |
32 |
0 |
3.2 (1.0) |
3.3 (2.6-3.6) |
1.5-5.4 |
Missingness changes the sample for a statistic
Most R summary functions need na.rm = TRUE when missing
values are present.
mean(airquality$Ozone)
## [1] NA
mean(airquality$Ozone, na.rm = TRUE)
## [1] 42.12931
Always report how many observations contributed to a statistic.
airquality |>
summarise(
across(
c(Ozone, Solar.R, Wind, Temp),
list(
observed = ~ sum(!is.na(.x)),
missing = ~ sum(is.na(.x)),
mean = ~ mean(.x, na.rm = TRUE)
),
.names = "{.col}_{.fn}"
)
)
## Ozone_observed Ozone_missing Ozone_mean Solar.R_observed Solar.R_missing
## 1 116 37 42.12931 146 7
## Solar.R_mean Wind_observed Wind_missing Wind_mean Temp_observed Temp_missing
## 1 185.9315 153 0 9.957516 153 0
## Temp_mean
## 1 77.88235
Continuous summaries across groups
ToothGrowth contains tooth length measurements for
guinea pigs receiving vitamin C through orange juice or ascorbic acid at
different doses.
data(ToothGrowth)
glimpse(ToothGrowth)
## Rows: 60
## Columns: 3
## $ len <dbl> 4.2, 11.5, 7.3, 5.8, 6.4, 10.0, 11.2, 11.2, 5.2, 7.0, 16.5, 16.5,…
## $ supp <fct> VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, V…
## $ dose <dbl> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, …
Summarise tooth length by supplement type.
tooth_by_supp <- ToothGrowth |>
summarise(
n = n(),
mean = mean(len),
sd = sd(len),
median = median(len),
q1 = quantile(len, 0.25),
q3 = quantile(len, 0.75),
.by = supp
)
tooth_by_supp
## supp n mean sd median q1 q3
## 1 VC 30 16.96333 8.266029 16.5 11.200 23.100
## 2 OJ 30 20.66333 6.605561 22.7 15.525 25.725
Add dose as another grouping variable.
ToothGrowth |>
mutate(dose = factor(dose)) |>
summarise(
n = n(),
mean = mean(len),
sd = sd(len),
median = median(len),
.by = c(supp, dose)
)
## supp dose n mean sd median
## 1 VC 0.5 10 7.98 2.746634 7.15
## 2 VC 1 10 16.77 2.515309 16.50
## 3 VC 2 10 26.14 4.797731 25.95
## 4 OJ 0.5 10 13.23 4.459709 12.25
## 5 OJ 1 10 22.70 3.910953 23.45
## 6 OJ 2 10 26.06 2.655058 25.95
These are descriptive differences. Statistical inference requires
additional methods and assumptions.
Efficient grouped summaries with across()
Use several continuous variables from mtcars and
summarise them by transmission.
grouped_continuous <- mtcars2 |>
summarise(
across(
c(mpg, hp, wt),
list(
n = ~ sum(!is.na(.x)),
mean = ~ mean(.x, na.rm = TRUE),
sd = ~ sd(.x, na.rm = TRUE),
median = ~ median(.x, na.rm = TRUE),
q1 = ~ quantile(.x, 0.25, na.rm = TRUE),
q3 = ~ quantile(.x, 0.75, na.rm = TRUE)
),
.names = "{.col}_{.fn}"
),
.by = transmission
)
grouped_continuous
## transmission mpg_n mpg_mean mpg_sd mpg_median mpg_q1 mpg_q3 hp_n hp_mean
## 1 Manual 13 24.39231 6.166504 22.8 21.00 30.4 13 126.8462
## 2 Automatic 19 17.14737 3.833966 17.3 14.95 19.2 19 160.2632
## hp_sd hp_median hp_q1 hp_q3 wt_n wt_mean wt_sd wt_median wt_q1
## 1 84.06232 109 66.0 113.0 13 2.411000 0.6169816 2.32 1.9350
## 2 53.90820 175 116.5 192.5 19 3.768895 0.7774001 3.52 3.4375
## wt_q3
## 1 2.7800
## 2 3.8425
A long format can be easier to format for reporting.
grouped_long <- mtcars2 |>
pivot_longer(
cols = c(mpg, hp, wt),
names_to = "variable",
values_to = "value"
) |>
summarise(
n = sum(!is.na(value)),
mean = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
median = median(value, na.rm = TRUE),
q1 = quantile(value, 0.25, na.rm = TRUE),
q3 = quantile(value, 0.75, na.rm = TRUE),
.by = c(transmission, variable)
)
grouped_long
## # A tibble: 6 × 8
## transmission variable n mean sd median q1 q3
## <fct> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Manual mpg 13 24.4 6.17 22.8 21 30.4
## 2 Manual hp 13 127. 84.1 109 66 113
## 3 Manual wt 13 2.41 0.617 2.32 1.94 2.78
## 4 Automatic mpg 19 17.1 3.83 17.3 15.0 19.2
## 5 Automatic hp 19 160. 53.9 175 116. 192.
## 6 Automatic wt 19 3.77 0.777 3.52 3.44 3.84
Variables have different units, so their raw SDs should not be
treated as directly comparable quantities.
Categorising a continuous variable
Sometimes a continuous variable must be grouped for reporting. With
mtcars, we can create fuel-economy bands.
mtcars_banded <- mtcars2 |>
mutate(
mpg_group = cut(
mpg,
breaks = c(-Inf, 15, 20, 25, Inf),
labels = c("15 or lower", ">15 to 20", ">20 to 25", "Above 25"),
right = TRUE
)
)
mtcars_banded |>
count(mpg_group, .drop = FALSE) |>
mutate(percent = n / sum(n) * 100)
## mpg_group n percent
## 1 15 or lower 6 18.75
## 2 >15 to 20 12 37.50
## 3 >20 to 25 8 25.00
## 4 Above 25 6 18.75
case_when() gives explicit control over boundaries.
mtcars_banded <- mtcars_banded |>
mutate(
mpg_group_explicit = case_when(
mpg <= 15 ~ "15 or lower",
mpg <= 20 ~ ">15 to 20",
mpg <= 25 ~ ">20 to 25",
TRUE ~ "Above 25"
),
mpg_group_explicit = factor(
mpg_group_explicit,
levels = c("15 or lower", ">15 to 20", ">20 to 25", "Above 25")
)
)
Categorisation loses information. Keep the original continuous
variable whenever later analysis benefits from its full variation.
Publication-style sample characteristics tables
A Table 1 commonly combines continuous and categorical variables.
Here transmission is used as the grouping variable.
Manual table components
Categorical component:
manual_categorical <- mtcars2 |>
select(transmission, cylinders, gears, engine) |>
pivot_longer(
cols = c(cylinders, gears, engine),
names_to = "variable",
values_to = "level"
) |>
count(transmission, variable, level, .drop = FALSE) |>
group_by(transmission, variable) |>
mutate(
denominator = sum(n),
percent = n / denominator * 100,
statistic = sprintf("%d (%.1f%%)", n, percent)
) |>
ungroup()
manual_categorical
## # A tibble: 42 × 7
## transmission variable level n denominator percent statistic
## <fct> <chr> <fct> <int> <int> <dbl> <chr>
## 1 Automatic cylinders 4 3 19 15.8 3 (15.8%)
## 2 Automatic cylinders 6 4 19 21.1 4 (21.1%)
## 3 Automatic cylinders 8 12 19 63.2 12 (63.2%)
## 4 Automatic cylinders 3 0 19 0 0 (0.0%)
## 5 Automatic cylinders 5 0 19 0 0 (0.0%)
## 6 Automatic cylinders V-shaped 0 19 0 0 (0.0%)
## 7 Automatic cylinders Straight 0 19 0 0 (0.0%)
## 8 Automatic engine 4 0 19 0 0 (0.0%)
## 9 Automatic engine 6 0 19 0 0 (0.0%)
## 10 Automatic engine 8 0 19 0 0 (0.0%)
## # ℹ 32 more rows
Continuous component:
manual_continuous <- mtcars2 |>
pivot_longer(
cols = c(mpg, hp, wt),
names_to = "variable",
values_to = "value"
) |>
summarise(
n = sum(!is.na(value)),
mean = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
median = median(value, na.rm = TRUE),
q1 = quantile(value, 0.25, na.rm = TRUE),
q3 = quantile(value, 0.75, na.rm = TRUE),
.by = c(transmission, variable)
) |>
mutate(
mean_sd = sprintf("%.1f (%.1f)", mean, sd),
median_iqr = sprintf("%.1f (%.1f-%.1f)", median, q1, q3)
)
manual_continuous
## # A tibble: 6 × 10
## transmission variable n mean sd median q1 q3 mean_sd
## <fct> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Manual mpg 13 24.4 6.17 22.8 21 30.4 24.4 (6.2)
## 2 Manual hp 13 127. 84.1 109 66 113 126.8 (84.1)
## 3 Manual wt 13 2.41 0.617 2.32 1.94 2.78 2.4 (0.6)
## 4 Automatic mpg 19 17.1 3.83 17.3 15.0 19.2 17.1 (3.8)
## 5 Automatic hp 19 160. 53.9 175 116. 192. 160.3 (53.9)
## 6 Automatic wt 19 3.77 0.777 3.52 3.44 3.84 3.8 (0.8)
## # ℹ 1 more variable: median_iqr <chr>
The same task with gtsummary
if (requireNamespace("gtsummary", quietly = TRUE)) {
table1 <- mtcars2 |>
select(mpg, cylinders, hp, wt, gears, engine, transmission) |>
gtsummary::tbl_summary(
by = transmission,
statistic = list(
gtsummary::all_continuous() ~ "{mean} ({sd})",
gtsummary::all_categorical() ~ "{n} ({p}%)"
),
missing = "ifany",
missing_text = "Missing",
label = list(
mpg ~ "Fuel economy, mpg",
cylinders ~ "Cylinders",
hp ~ "Gross horsepower",
wt ~ "Weight (1000 lbs)",
gears ~ "Forward gears",
engine ~ "Engine shape"
),
digits = list(
gtsummary::all_continuous() ~ 1,
gtsummary::all_categorical() ~ c(0, 1)
)
) |>
gtsummary::add_overall(last = FALSE)
table1
} else {
message("Install the gtsummary package to render the publication-style Table 1 example.")
}
| Characteristic |
Overall
N = 32 |
Automatic
N = 19 |
Manual
N = 13 |
| Fuel economy, mpg |
20.1 (6.0) |
17.1 (3.8) |
24.4 (6.2) |
| Cylinders |
|
|
|
| 4 |
11 (34.4%) |
3 (15.8%) |
8 (61.5%) |
| 6 |
7 (21.9%) |
4 (21.1%) |
3 (23.1%) |
| 8 |
14 (43.8%) |
12 (63.2%) |
2 (15.4%) |
| Gross horsepower |
146.7 (68.6) |
160.3 (53.9) |
126.8 (84.1) |
| Weight (1000 lbs) |
3.2 (1.0) |
3.8 (0.8) |
2.4 (0.6) |
| Forward gears |
|
|
|
| 3 |
15 (46.9%) |
15 (78.9%) |
0 (0.0%) |
| 4 |
12 (37.5%) |
4 (21.1%) |
8 (61.5%) |
| 5 |
5 (15.6%) |
0 (0.0%) |
5 (38.5%) |
| Engine shape |
|
|
|
| V-shaped |
18 (56.3%) |
12 (63.2%) |
6 (46.2%) |
| Straight |
14 (43.8%) |
7 (36.8%) |
7 (53.8%) |
P-values can be added when they serve the analytical purpose, but a
descriptive table does not require them automatically.
Missing data are part of the result
Use airquality to create a clear missingness table.
missing_long <- airquality |>
summarise(
across(everything(), ~ sum(is.na(.x)))
) |>
pivot_longer(
everything(),
names_to = "variable",
values_to = "missing_n"
) |>
mutate(
total_n = nrow(airquality),
observed_n = total_n - missing_n,
missing_percent = missing_n / total_n * 100
) |>
arrange(desc(missing_n), variable)
missing_long
## # A tibble: 6 × 5
## variable missing_n total_n observed_n missing_percent
## <chr> <int> <int> <int> <dbl>
## 1 Ozone 37 153 116 24.2
## 2 Solar.R 7 153 146 4.58
## 3 Day 0 153 153 0
## 4 Month 0 153 153 0
## 5 Temp 0 153 153 0
## 6 Wind 0 153 153 0
Compare missingness by month.
airquality2 |>
summarise(
observations = n(),
ozone_missing = sum(is.na(Ozone)),
solar_missing = sum(is.na(Solar.R)),
.by = Month
)
## Month observations ozone_missing solar_missing
## 1 May 31 5 4
## 2 June 30 21 0
## 3 July 31 5 0
## 4 August 31 5 3
## 5 September 30 1 0
If missingness differs across groups, observed-value summaries may
represent different subsets of the data.
Presentation that helps rather than distracts
A good table has clear labels, visible denominators, sensible
rounding, and no unnecessary precision.
transmission_display <- mtcars2 |>
count(transmission, .drop = FALSE) |>
mutate(
percent = n / sum(n) * 100,
`n (%)` = sprintf("%d (%.1f%%)", n, percent)
) |>
select(`Transmission` = transmission, `n (%)`)
transmission_display |>
kable(caption = "Transmission type in mtcars")
Transmission type in mtcars
| Automatic |
19 (59.4%) |
| Manual |
13 (40.6%) |
Keep calculations numeric until the presentation stage. Formatting
early makes later calculations more difficult.
Common mistakes
A percentage with no denominator
A percentage is incomplete unless readers know what observations
formed its denominator.
Row and column percentages are swapped
Row percentages and column percentages answer different questions.
Check which dimension sums to 100%.
Missing values disappear silently
Functions such as table() and numerical summaries can
exclude missing observations. Always show missing counts or observed
sample sizes when relevant.
A skewed variable is summarised only by mean and SD
Inspect distributions. Median and IQR can show a different and often
useful view of a skewed variable.
Standard deviation is confused with standard error
SD describes variation among observations. SE describes uncertainty
in an estimated mean. They answer different questions.
A zero-count category is called missing
A defined level with zero observations is not an NA
value.
Continuous variables are categorised without a reason
Cut-points simplify reporting but discard information and create
artificial boundaries.
Every possible statistic is included
More statistics do not necessarily make a better table. Select
summaries that fit the variable, research question and audience.
Descriptive differences are interpreted causally
A descriptive table reports observed differences. It does not by
itself establish causation, adjustment for confounding, or statistical
significance.
Complete worked example using only mtcars
Prepare labelled variables
analysis_data <- mtcars |>
mutate(
transmission = factor(am, c(0, 1), c("Automatic", "Manual")),
cylinders = factor(cyl),
gears = factor(gear),
engine = factor(vs, c(0, 1), c("V-shaped", "Straight"))
)
glimpse(analysis_data)
## Rows: 32
## Columns: 15
## $ mpg <dbl> 21.0, 21.0, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.…
## $ cyl <dbl> 6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, …
## $ disp <dbl> 160.0, 160.0, 108.0, 258.0, 360.0, 225.0, 360.0, 146.7, 1…
## $ hp <dbl> 110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, …
## $ drat <dbl> 3.90, 3.90, 3.85, 3.08, 3.15, 2.76, 3.21, 3.69, 3.92, 3.9…
## $ wt <dbl> 2.620, 2.875, 2.320, 3.215, 3.440, 3.460, 3.570, 3.190, 3…
## $ qsec <dbl> 16.46, 17.02, 18.61, 19.44, 17.02, 20.22, 15.84, 20.00, 2…
## $ vs <dbl> 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, …
## $ am <dbl> 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, …
## $ gear <dbl> 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 4, 4, …
## $ carb <dbl> 4, 4, 1, 1, 2, 1, 4, 2, 2, 4, 4, 3, 3, 3, 4, 4, 4, 1, 2, …
## $ transmission <fct> Manual, Manual, Manual, Automatic, Automatic, Automatic, …
## $ cylinders <fct> 6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, …
## $ gears <fct> 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 4, 4, …
## $ engine <fct> V-shaped, V-shaped, Straight, Straight, V-shaped, Straigh…
One-way categorical descriptions
one_way_variables <- c("transmission", "cylinders", "gears", "engine")
one_way_tables <- lapply(
one_way_variables,
function(variable) {
tab <- table(analysis_data[[variable]], useNA = "ifany")
tibble(
variable = variable,
level = names(tab),
count = as.integer(tab)
) |>
mutate(percent = count / sum(count) * 100)
}
) |>
bind_rows()
one_way_tables
## # A tibble: 10 × 4
## variable level count percent
## <chr> <chr> <int> <dbl>
## 1 transmission Automatic 19 59.4
## 2 transmission Manual 13 40.6
## 3 cylinders 4 11 34.4
## 4 cylinders 6 7 21.9
## 5 cylinders 8 14 43.8
## 6 gears 3 15 46.9
## 7 gears 4 12 37.5
## 8 gears 5 5 15.6
## 9 engine V-shaped 18 56.2
## 10 engine Straight 14 43.8
Main cross-tabulation
transmission_engine <- analysis_data |>
count(transmission, engine, .drop = FALSE) |>
group_by(transmission) |>
mutate(
denominator = sum(n),
percent = n / denominator * 100,
display = sprintf("%d (%.1f%%)", n, percent)
) |>
ungroup()
transmission_engine
## # A tibble: 4 × 6
## transmission engine n denominator percent display
## <fct> <fct> <int> <int> <dbl> <chr>
## 1 Automatic V-shaped 12 19 63.2 12 (63.2%)
## 2 Automatic Straight 7 19 36.8 7 (36.8%)
## 3 Manual V-shaped 6 13 46.2 6 (46.2%)
## 4 Manual Straight 7 13 53.8 7 (53.8%)
Continuous variables overall
worked_continuous <- c("mpg", "hp", "wt")
overall_continuous <- lapply(
worked_continuous,
function(variable) {
continuous_summary(analysis_data[[variable]]) |>
mutate(variable = variable, .before = 1)
}
) |>
bind_rows()
overall_continuous
## # A tibble: 3 × 11
## variable n_total n_observed n_missing mean sd median q1 q3 min
## <chr> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mpg 32 32 0 20.1 6.03 19.2 15.4 22.8 10.4
## 2 hp 32 32 0 147. 68.6 123 96.5 180 52
## 3 wt 32 32 0 3.22 0.978 3.32 2.58 3.61 1.51
## # ℹ 1 more variable: max <dbl>
Continuous variables by transmission
worked_by_transmission <- analysis_data |>
pivot_longer(
cols = all_of(worked_continuous),
names_to = "variable",
values_to = "value"
) |>
summarise(
n = sum(!is.na(value)),
missing = sum(is.na(value)),
mean = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
median = median(value, na.rm = TRUE),
q1 = quantile(value, 0.25, na.rm = TRUE),
q3 = quantile(value, 0.75, na.rm = TRUE),
.by = c(transmission, variable)
)
worked_by_transmission
## # A tibble: 6 × 9
## transmission variable n missing mean sd median q1 q3
## <fct> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Manual mpg 13 0 24.4 6.17 22.8 21 30.4
## 2 Manual hp 13 0 127. 84.1 109 66 113
## 3 Manual wt 13 0 2.41 0.617 2.32 1.94 2.78
## 4 Automatic mpg 19 0 17.1 3.83 17.3 15.0 19.2
## 5 Automatic hp 19 0 160. 53.9 175 116. 192.
## 6 Automatic wt 19 0 3.77 0.777 3.52 3.44 3.84
Compact mixed report table
continuous_report_rows <- worked_by_transmission |>
mutate(
statistic = sprintf("%.1f (%.1f)", mean, sd),
row_label = recode(
variable,
mpg = "Fuel economy, mpg: Mean (SD)",
hp = "Horsepower: Mean (SD)",
wt = "Weight, 1000 lbs: Mean (SD)"
)
) |>
select(row_label, transmission, statistic)
categorical_report_rows <- analysis_data |>
count(transmission, cylinders, .drop = FALSE) |>
group_by(transmission) |>
mutate(
denominator = sum(n),
statistic = sprintf("%d (%.1f%%)", n, n / denominator * 100),
row_label = paste("Cylinders:", cylinders)
) |>
ungroup() |>
select(row_label, transmission, statistic)
worked_report <- bind_rows(
continuous_report_rows,
categorical_report_rows
) |>
pivot_wider(names_from = transmission, values_from = statistic)
worked_report |>
kable(
col.names = c("Characteristic", "Automatic", "Manual"),
caption = "Selected characteristics of cars by transmission type"
)
Selected characteristics of cars by transmission type
| Fuel economy, mpg: Mean (SD) |
24.4 (6.2) |
17.1 (3.8) |
| Horsepower: Mean (SD) |
126.8 (84.1) |
160.3 (53.9) |
| Weight, 1000 lbs: Mean (SD) |
2.4 (0.6) |
3.8 (0.8) |
| Cylinders: 4 |
8 (61.5%) |
3 (15.8%) |
| Cylinders: 6 |
3 (23.1%) |
4 (21.1%) |
| Cylinders: 8 |
2 (15.4%) |
12 (63.2%) |
Consistency checks
stopifnot(nrow(analysis_data) == nrow(mtcars))
stopifnot(sum(one_way_tables$count[one_way_tables$variable == "transmission"]) == nrow(mtcars))
analysis_data |>
summarise(
cars = n(),
mpg_observed = sum(!is.na(mpg)),
hp_observed = sum(!is.na(hp)),
wt_observed = sum(!is.na(wt))
)
## cars mpg_observed hp_observed wt_observed
## 1 32 32 32 32
Practical reference
| How many observations are in each category? |
table(), count() |
Missing and zero-count levels |
| What percentage belongs to each category? |
prop.table(), n / sum(n) |
State the denominator |
| How are two categorical variables distributed together? |
table(x, y), xtabs() |
Row, column, or overall percentages |
| Do categories differ across strata? |
Multi-way table(), stratified tables |
Sparse cells and subgroup denominators |
| What is the typical continuous value? |
mean(), median() |
Distribution and scientific meaning |
| How variable are observations? |
sd(), IQR(), range(),
quantile() |
Units and influential extremes |
| How many values are observed or missing? |
sum(!is.na(x)), sum(is.na(x)) |
Variable-specific denominators |
| How does a continuous variable differ by group? |
.by or group_by() with
summarise() |
Descriptive, not automatically significant |
| How can several variables be summarised efficiently? |
across(), pivot_longer() |
Scales and readable output |
| How can continuous values be grouped? |
cut(), case_when() |
Boundaries and information loss |
| How can a simple table be presented? |
knitr::kable() |
Round at presentation stage |
| How can a mixed publication table be created? |
Manual dplyr workflow or gtsummary::tbl_summary() |
Labels, missingness and purpose |
Conclusion
A useful descriptive table begins with a question and a denominator.
Counts show how observations are distributed, while percentages make
those counts comparable only when their denominators are clear. Means
and SDs describe one aspect of continuous data; medians, IQRs, ranges,
percentiles and simple graphics can reveal features hidden by the
mean.
The examples in this tutorial use sample data already available in R
rather than fabricated participant data. That makes the code immediately
reproducible and gives students datasets they can load directly in a
fresh R session. Treat missingness as part of the result, distinguish
zero-count levels from missing observations, inspect distributions
before choosing summaries, and keep calculations numeric until the final
presentation stage.