Before any model, look at the data. Exploration has four jobs (Zuur et al. 2010): find errors and outliers, understand the structure of the data — what counts as a replicate —, see the distributions that models will have to describe, and detect collinearity between variables that a model will struggle to separate. Skipping it does not save time: the problems are still there, and they surface later as strange model results.
The simulated ponds of the previous chapters are too clean for this. This chapter uses real data, with all the problems real data have.
Like Chapter 10 and Chapter 11, this chapter is inspired by ecoR (Oliveira 2026), the R course wiki of the Instituto de Biociências, Universidade de São Paulo, produced and maintained by Alexandre Adalardo de Oliveira. Its sequence — checking the data frame, descriptive statistics, univariate graphs, relationships between variables, and the Anscombe quartet — follows ecoR’s exploratory-analysis module. The text, data, examples and exercises here are written from scratch for Julia; the authors of ecoR have not reviewed this book.
If you read Portuguese and want to go further, the original material is freely available (CC BY-SA 4.0):
ATLANTIC BIRD TRAITS (Rodrigues et al. 2019) compiles morphological measurements of birds from the Atlantic forests of South America: 72 483 records, gathered from museum collections, ringing projects and the literature between 1820 and 2018, by 179 authors. Its data paper states that there are no copyright or proprietary restrictions on its use, and asks that it be cited when used in publications or teaching — as here.
data/atlantic_bird_traits.csv keeps all records and 21 of the original 83 columns, with their original names and codes. Nothing is cleaned: finding the problems is the point. scripts/prepare_birds.jl documents how the subset was made from the original file, which is in the paper’s supporting information.
using TidierData
using DataFrames, CSV
using Statistics, StatsBase
using GLM, RDatasets, Distributions
using CairoMakie
CairoMakie.activate!(type = "png")
bird_file = joinpath("..", "data", "atlantic_bird_traits.csv")"../data/atlantic_bird_traits.csv"The answers are numbers computed from the bird data. Start each session with:
include("exercises/check.jl")
using TidierData, DataFrames, CSV, Statistics, StatsBase, GLM, Distributions
bird_file = joinpath("data", "atlantic_bird_traits.csv")The first section shows how to read the file; the later exercises assume you have read it that way, into birds, with the columns renamed.
Read the file as it comes and look at the type of each column:
raw = CSV.read(bird_file, DataFrame)
size(raw)(72483, 21)first(describe(raw, :eltype, :nmissing), 8)| Row | variable | eltype | nmissing |
|---|---|---|---|
| Symbol | DataType | Int64 | |
| 1 | ID_ABT | String7 | 0 |
| 2 | Order | String31 | 0 |
| 3 | Family | String31 | 0 |
| 4 | Binomial | String31 | 0 |
| 5 | Body_mass.g. | String7 | 0 |
| 6 | Wing_length.mm. | String7 | 0 |
| 7 | Tail_length.mm. | String7 | 0 |
| 8 | Tarsus_length.mm. | String7 | 0 |
Body mass, wing length, tail length: all read as text (String), and not one value is missing. The reason is in the file itself — missing values are written as the two letters NA, which CSV.jl does not treat as missing by default. One non-number in a column makes the whole column text. Tell CSV.read which code means missing:
birds = CSV.read(bird_file, DataFrame; missingstring = "NA")
first(describe(birds, :eltype, :nmissing), 8)| Row | variable | eltype | nmissing |
|---|---|---|---|
| Symbol | Type | Int64 | |
| 1 | ID_ABT | String7 | 0 |
| 2 | Order | String31 | 0 |
| 3 | Family | String31 | 0 |
| 4 | Binomial | String31 | 0 |
| 5 | Body_mass.g. | Union{Missing, Float64} | 11664 |
| 6 | Wing_length.mm. | Union{Missing, Float64} | 51134 |
| 7 | Tail_length.mm. | Union{Missing, Float64} | 26183 |
| 8 | Tarsus_length.mm. | Union{Missing, Float64} | 50742 |
In R
read.csv treats "NA" as missing by default, so in R this problem would pass unnoticed. Here it is visible: a numeric column that arrives as text is the first sign that the file uses a missing-value code you have not declared. Other common codes are "", "-", ".", -999 and 99999; missingstring also accepts a vector of them.
The column names mix units, dots and capitals, like Body_mass.g.. They are awkward to type and Tidier macros cannot use them directly. Rename them once, at the start:
rename!(birds,
"ID_ABT" => :id, "Order" => :order, "Family" => :family, "Binomial" => :species,
"Body_mass.g." => :mass_g, "Wing_length.mm." => :wing_mm,
"Tail_length.mm." => :tail_mm, "Tarsus_length.mm." => :tarsus_mm,
"Bill_length.mm." => :bill_mm, "Age" => :age, "Sex" => :sex, "Status" => :status,
"Recapture" => :recapture, "Ring" => :ring,
"Latitude_decimal_degrees" => :lat, "Longitude_decimal_degrees" => :lon,
"Altitude" => :altitude, "Year" => :year, "Main_researcher" => :source,
"Outside.range" => :outside_range, "Obs.spp" => :species_note)
names(birds)21-element Vector{String}:
"id"
"order"
"family"
"species"
"mass_g"
"wing_mm"
"tail_mm"
"tarsus_mm"
"bill_mm"
"age"
⋮
"recapture"
"ring"
"lat"
"lon"
"altitude"
"year"
"source"
"outside_range"
"species_note"How many records have a body mass? Pass the number to check("7.1", ...).
count(!ismissing, v) counts the non-missing values of a vector. It only means something once the file has been read with missingstring = "NA".
n_mass = count(!ismissing, birds.mass_g)60819check("7.1", n_mass)✓ Correct!
60 819 of 72 483: body mass is the most commonly recorded trait, and still missing for one record in six.
Categorical columns deserve a table of their levels before anything else:
(age = countmap(birds.age), sex = countmap(birds.sex))(age = Dict{Union{Missing, String15}, Int64}(missing => 4935, String15("Juvenile") => 8708, String15("Adult") => 52111, String15("Unknown") => 6727, String15("Nestling") => 2), sex = Dict{Union{Missing, String7}, Int64}(String7("Female") => 14586, missing => 3072, String7("Male") => 18420, String7("Unknown") => 36405))Parentheses with names make a named tuple: a tuple (Section 4.1) in which each value has a name. It returns in Section 10.1.
There are two kinds of missing here. missing means the information was never recorded; "Unknown" means someone looked at the bird and could not tell. They are not the same thing: a study of sex ratios should count the unknowns, but not the unrecorded.
Other columns disagree with themselves about spelling:
(recapture = countmap(birds.recapture), outside_range = countmap(birds.outside_range))(recapture = Dict{Union{Missing, String3}, Int64}(missing => 33682, String3("Yes") => 8556, String3("no") => 30245), outside_range = Dict{String3, Int64}("Yes" => 2725, "no" => 69758))"no" and "Yes": the capitalization differs between levels. Here it is harmless, because each level has only one spelling, but a single "yes" somewhere would create a third level that grouping would treat as different. Normalizing is cheap:
birds.recapture = passmissing(lowercase).(birds.recapture)
birds.outside_range = lowercase.(birds.outside_range)
unique(birds.recapture)3-element Vector{Union{Missing, String}}:
"no"
missing
"yes"passmissing(f) applies f to values and leaves missing alone.
Among records where sex was recorded (not missing), what proportion were sexed as "Male" or "Female"? Pass it to check("7.2", ...).
Keep the non-missing values with collect(skipmissing(birds.sex)), then take the mean of in(("Male", "Female")).(...).
sex_recorded = collect(skipmissing(birds.sex))
p_sexed = mean(in(("Male", "Female")).(sex_recorded))0.47551540822059907check("7.2", p_sexed)✓ Correct!
Less than half: many Atlantic forest passerines cannot be sexed in the hand. Any analysis by sex is an analysis of a non-random half of the data.
A row is a record, but records are not independent observations. They are grouped in several ways at once.
By species. The records are very unequally spread:
per_species = sort(combine(groupby(birds, :species), nrow => :n), :n, rev = true)
(names = nrow(per_species),
top_five = first(per_species, 5),
with_one_record = count(==(1), per_species.n))(names = 780, top_five = 5×2 DataFrame Row │ species n │ String31 Int64 ─────┼───────────────────────────────── 1 │ Chiroxiphia caudata 1857 2 │ Turdus albicollis 1686 3 │ Trichothraupis melanops 1493 4 │ Turdus leucomelas 1475 5 │ Basileuterus culicivorus 1419, with_one_record = 68)
780 names — some identified only to genus, such as Elaenia sp. — while the data paper reports 711 species. Before counting species, decide what counts as one. The ten most recorded species hold a fifth of all records; 68 names appear once.
By individual. A ringed bird caught again appears again:
ringed = dropmissing(birds[:, [:ring, :species]], :ring)
per_ring = combine(groupby(ringed, :ring), :species => (s -> length(unique(s))) => :n_species,
nrow => :n)
first(sort(per_ring[per_ring.n_species .> 1, :], :n, rev = true), 6)| Row | ring | n_species | n |
|---|---|---|---|
| String31 | Int64 | Int64 | |
| 1 | BRD | 25 | 31 |
| 2 | VRD | 24 | 27 |
| 3 | AZD | 21 | 26 |
| 4 | AZE | 22 | 24 |
| 5 | IND1 | 4 | 15 |
| 6 | BRE | 14 | 14 |
Some “rings” are shared by more than twenty species. They are not ring numbers but colour-band codes such as BRD or BLUE RIGHT, which many birds carry. A naive count of recaptures would take all those birds for one individual.
By source and method. Each record comes from one of 65 main researchers, and from live birds or museum skins, over two centuries:
combine(groupby(dropmissing(birds[:, [:status, :year]]), :status),
nrow => :n, :year => minimum => :first, :year => median => :median,
:year => maximum => :last)| Row | status | n | first | median | last |
|---|---|---|---|---|---|
| String7 | Int64 | Int64 | Float64 | Int64 | |
| 1 | live | 57862 | 1990 | 2010.0 | 2018 |
| 2 | museum | 11152 | 1820 | 2007.0 | 2017 |
Each of these groupings can make records from the same group more alike than records from different groups — the dependence that the mixed models of Section 8.7 account for. Exploring the structure tells you which grouping terms a model needs.
Count the individuals recorded more than once: rings that appear in more than one record of the same species, ignoring ring codes made only of letters (colour-band codes). Pass the number to check("7.3", ...).
Drop the letter-only codes with occursin(r"^[A-Za-z]+$", ring), which is true for them. Then group by ring and species, count rows, and count the groups with more than one.
real_rings = ringed[.!occursin.(r"^[A-Za-z]+$", ringed.ring), :]
per_bird = combine(groupby(real_rings, [:ring, :species]), nrow => :n)
n_recaptured = count(>(1), per_bird.n)4503check("7.3", n_recaptured)✓ Correct!
About 4 500 birds contribute two records or more, one of them up to 43. Treating their records as independent would overstate the sample size — and the colour codes would have added 27 false “individuals” to the count.
Body mass across all birds spans three orders of magnitude, from hummingbirds to large seabirds. On a linear scale almost everything piles up at the left; on a log scale the shape becomes visible:
masses = collect(skipmissing(birds.mass_g))
fig = Figure(size = (860, 340))
ax1 = Axis(fig[1, 1], xlabel = "Body mass (g)", ylabel = "Records", title = "(a) Linear")
hist!(ax1, masses, bins = 60, color = (:steelblue, 0.7))
ax2 = Axis(fig[1, 2], xlabel = "Body mass (g, log scale)", title = "(b) Logarithmic",
xscale = log10)
hist!(ax2, masses, bins = 10 .^ range(0, 3.6, length = 60), color = (:steelblue, 0.7))
fig
Sizes, areas and concentrations are bounded at zero and grow multiplicatively, so they are usually right-skewed and closer to symmetric on a log scale. That is why models of such variables so often use their logarithm.
Within one species the picture is different. For the white-necked thrush (Turdus albicollis), one of the most recorded species, compare the distribution with a normal curve of the same mean and standard deviation:
thrush = collect(skipmissing(birds.mass_g[birds.species .== "Turdus albicollis"]))
μ, σ = mean(thrush), std(thrush)
fig = Figure(size = (860, 340))
ax1 = Axis(fig[1, 1], xlabel = "Body mass (g)", ylabel = "Density", title = "(a) Histogram")
hist!(ax1, thrush, bins = 40, normalization = :pdf, color = (:seagreen, 0.6))
xs = range(minimum(thrush), maximum(thrush), length = 200)
lines!(ax1, xs, @.(exp(-(xs - μ)^2 / (2σ^2)) / (σ * sqrt(2π))), color = :black)
ax2 = Axis(fig[1, 2], xlabel = "Normal quantiles", ylabel = "Body mass (g)",
title = "(b) QQ plot")
n = length(thrush)
z = quantile.(Normal(), ((1:n) .- 0.5) ./ n)
scatter!(ax2, z, sort(thrush), color = (:seagreen, 0.5), markersize = 5)
ablines!(ax2, μ, σ, color = :black, linestyle = :dash)
fig
Close to normal in the middle, with a few birds far lighter or heavier than the rest — visible as the points that leave the line at both ends of the QQ plot. Those tails are the subject of the next section.
Compute the median and the mean body mass of Turdus albicollis and pass them, in that order, as a vector to check("7.4", ...). Which is less affected by the extreme birds?
thrush above already holds the masses. [median(x), mean(x)] builds the vector.
thrush_centre = [median(thrush), mean(thrush)]2-element Vector{Float64}:
66.0
66.74542818610924check("7.4", thrush_centre)✓ Correct!
66 and 66.7 g: nearly equal, as expected for a roughly symmetric distribution. The median ignores how extreme the tails are; the mean is pulled by them. With a few records at 26 g or 114 g, the mean moves slightly; with the gross errors of the next section, it would move much more.
An overall histogram cannot show an error in a 14 g bird recorded as 1 240 g: 1 240 g is a perfectly ordinary mass for a large bird. Outliers in trait data are outliers relative to the species. A simple, robust flag compares each mass with the median of its species:
weighed = dropmissing(birds[:, [:species, :mass_g, :age, :status, :source]], :mass_g)
weighed = transform(groupby(weighed, :species),
:mass_g => median => :species_median, nrow => :n_species)
weighed.ratio = weighed.mass_g ./ weighed.species_median
suspicious = sort(weighed[(weighed.n_species .>= 10) .& (weighed.ratio .> 3), :],
:ratio, rev = true)
first(suspicious[:, [:species, :mass_g, :species_median, :ratio, :age, :status]], 8)| Row | species | mass_g | species_median | ratio | age | status |
|---|---|---|---|---|---|---|
| String31 | Float64 | Float64 | Float64 | String15? | String7? | |
| 1 | Thlypopsis sordida | 1240.0 | 14.0 | 88.5714 | Unknown | museum |
| 2 | Thlypopsis sordida | 1070.0 | 14.0 | 76.4286 | Unknown | museum |
| 3 | Corythopis delalandi | 214.0 | 15.0 | 14.2667 | Adult | live |
| 4 | Tityra cayana | 940.0 | 73.0 | 12.8767 | Juvenile | museum |
| 5 | Tityra cayana | 840.0 | 73.0 | 11.5068 | Juvenile | museum |
| 6 | Platyrinchus mystaceus | 102.0 | 9.0 | 11.3333 | missing | missing |
| 7 | Thalurania glaucopis | 50.0 | 4.5 | 11.1111 | Adult | live |
| 8 | Phaethornis eurynome | 55.0 | 5.0 | 11.0 | Adult | live |
The top of the list shows several kinds of problem:
Now the other tail:
first(sort(weighed[(weighed.n_species .>= 10) .& (weighed.ratio .< 1/3), :], :ratio)[:,
[:species, :mass_g, :species_median, :ratio, :age, :status]], 6)| Row | species | mass_g | species_median | ratio | age | status |
|---|---|---|---|---|---|---|
| String31 | Float64 | Float64 | Float64 | String15? | String7? | |
| 1 | Crypturellus parvirostris | 8.5 | 185.0 | 0.0459459 | Unknown | museum |
| 2 | Aramides saracura | 22.5 | 346.0 | 0.0650289 | Unknown | museum |
| 3 | Aramides saracura | 22.5 | 346.0 | 0.0650289 | Juvenile | museum |
| 4 | Aramides saracura | 24.0 | 346.0 | 0.0693642 | Unknown | museum |
| 5 | Aramides saracura | 24.0 | 346.0 | 0.0693642 | Juvenile | museum |
| 6 | Asio clamator | 34.5 | 470.0 | 0.0734043 | Unknown | museum |
Here some of the “outliers” are real. A 24 g slaty-breasted wood-rail (Aramides saracura, median 346 g) marked as juvenile is plausibly a chick. Deleting it would remove a real bird; keeping it in an analysis of adult size would distort it.
A statistical rule can only say that a value is unusual, never that it is wrong. Flag suspicious values, look at them one by one, go back to the source where you can, and decide — and write down what you decided and why, in code, so that the decision can be reviewed. For many analyses the right choice is to restrict the data to what the question is about (adults, say) rather than to remove individual points.
The data also carry a quieter problem, precision. Most masses are whole numbers:
mean(weighed.mass_g .== round.(weighed.mass_g))0.721797464608099572% of masses are integers, which means many were weighed to the nearest gram. For a 3 g hummingbird, rounding to the gram is an error of up to 17%; for a thrush, under 1%. Precision that varies with the size of the animal is easy to miss and worth knowing before you compare small species.
Count the records whose mass is more than 3 times, or less than a third of, the median of their species, among species with at least 10 weighed records. Pass the number to check("7.5", ...).
weighed already has ratio and n_species. Combine the conditions with .& and .|, and count the trues with count.
n_suspicious = count((weighed.n_species .>= 10) .&
((weighed.ratio .> 3) .| (weighed.ratio .< 1/3)))291check("7.5", n_suspicious)✓ Correct!
291 of about 60 000 weighed records, 0.5%: few, but enough to wreck a species mean. The threshold of 3 is a choice, not a rule; with 2 you would flag more than twice as many.
Morphological traits go together: big birds have long wings, long tails, long legs and long bills. To see how strongly, summarize each species by the median of each trait — the median, because of the outliers above — and compare species:
traits = [:mass_g, :wing_mm, :tail_mm, :tarsus_mm, :bill_mm]
median_or_missing(x) = (v = collect(skipmissing(x)); isempty(v) ? missing : median(v))
by_species = dropmissing(combine(groupby(birds, :species),
[t => median_or_missing => t for t in traits]...))
log_traits = DataFrame(log.(Matrix(by_species[:, traits])), traits)
nrow(by_species)538round.(cor(Matrix(log_traits)); digits = 2)5×5 Matrix{Float64}:
1.0 0.93 0.71 0.81 0.72
0.93 1.0 0.77 0.72 0.64
0.71 0.77 1.0 0.56 0.5
0.81 0.72 0.56 1.0 0.64
0.72 0.64 0.5 0.64 1.0labels = ["log mass", "log wing", "log tail", "log tarsus", "log bill"]
fig = Figure(size = (820, 820))
for i in 1:5, j in 1:5
ax = Axis(fig[i, j], xlabel = i == 5 ? labels[j] : "", ylabel = j == 1 ? labels[i] : "",
xticklabelsvisible = i == 5, yticklabelsvisible = j == 1 && i != 1)
if i == j # diagonal: the distribution of each trait (counts)
hist!(ax, log_traits[:, i], bins = 25, color = (:grey60, 0.8))
else
scatter!(ax, log_traits[:, j], log_traits[:, i], markersize = 3,
color = (:steelblue, 0.5))
end
end
fig
Mass and wing length are almost a single variable across species. If both entered a model as predictors, the model could not tell their effects apart: the estimates would be unstable and their standard errors large. The variance inflation factor measures this. For each predictor, regress it on all the others; with \(R^2\) from that regression, \(\text{VIF} = 1 / (1 - R^2)\):
function vif(df, col)
others = setdiff(propertynames(df), [col])
f = Term(col) ~ sum(Term.(others))
return 1 / (1 - r2(lm(f, df)))
end
DataFrame(trait = traits, vif = round.([vif(log_traits, t) for t in traits]; digits = 1))| Row | trait | vif |
|---|---|---|
| Symbol | Float64 | |
| 1 | mass_g | 11.9 |
| 2 | wing_mm | 9.7 |
| 3 | tail_mm | 2.5 |
| 4 | tarsus_mm | 3.0 |
| 5 | bill_mm | 2.1 |
Zuur et al. (2010) use a cut-off of 3; other authors use 5 or 10. Mass and wing length are above all of them. The usual responses are to keep only one of each collinear group, to combine them (a size index), or to ask a question that does not need them separately.
What is the correlation between log body mass and log wing length across species, using the species medians above? Pass it to check("7.6", ...).
cor(x, y) on two columns of log_traits.
r_mass_wing = cor(log_traits.mass_g, log_traits.wing_mm)0.9326447165263585check("7.6", r_mass_wing)✓ Correct!
About 0.93. Across species, these two traits carry almost the same information. Within a species it is much weaker — about 0.25 among individual Turdus albicollis — because individuals of one species vary little in size: collinearity depends on the scale of the question.
Summaries are not a substitute for looking. Anscombe (1973) built four small data sets with the same means, variances, correlations and regression lines:
anscombe = dataset("datasets", "anscombe")
DataFrame(set = 1:4,
mean_x = [mean(anscombe[:, "X$i"]) for i in 1:4],
mean_y = [round(mean(anscombe[:, "Y$i"]); digits = 2) for i in 1:4],
cor_xy = [round(cor(anscombe[:, "X$i"], anscombe[:, "Y$i"]); digits = 3) for i in 1:4])| Row | set | mean_x | mean_y | cor_xy |
|---|---|---|---|---|
| Int64 | Float64 | Float64 | Float64 | |
| 1 | 1 | 9.0 | 7.5 | 0.816 |
| 2 | 2 | 9.0 | 7.5 | 0.816 |
| 3 | 3 | 9.0 | 7.5 | 0.816 |
| 4 | 4 | 9.0 | 7.5 | 0.817 |
fig = Figure(size = (760, 560))
for i in 1:4
x, y = anscombe[:, "X$i"], anscombe[:, "Y$i"]
ax = Axis(fig[(i - 1) ÷ 2 + 1, (i - 1) % 2 + 1], title = "Set $i",
limits = (3, 20, 3, 14))
scatter!(ax, x, y, color = :steelblue)
b = coef(lm(@formula(y ~ x), DataFrame(; x, y)))
ablines!(ax, b[1], b[2], color = :darkorange)
end
fig
The first looks like what a linear model assumes. The second is a curve. In the third one point pulls the line away from an otherwise perfect relationship. In the fourth the whole slope depends on a single observation. Only the plots show it.
lon, lat) and colour them by outside_range. Are the records flagged as outside the species’ range where you would expect? Are there points outside the Atlantic forest region altogether?"Unknown" sex.