using TidierData
using DataFrames, CSV
using GLM, StatsModels, StatsBase, Statistics, Distributions
using HypothesisTests
using CairoMakie
using MixedModels
CairoMakie.activate!(type = "png")
ponds = @chain CSV.read(joinpath("..", "data", "ponds.csv"), DataFrame) begin
@drop_missing(canopy)
@mutate(log_area = log(area_m2))
end
nrow(ponds)173In the Your turn boxes you fit a model and pass a number extracted from it — a p-value, coefficients, an interval — to check. All of them use the ponds with canopy recorded, as in this chapter. Start each session with:
include("exercises/check.jl")
using TidierData, DataFrames, CSV, GLM, StatsModels, Statistics
using HypothesisTests, MixedModels
ponds = @chain CSV.read(joinpath("data", "ponds.csv"), DataFrame) begin
@drop_missing(canopy)
@mutate(log_area = log(area_m2))
endBefore models, the basics. HypothesisTests.jl covers the usual repertoire.
perm = @chain ponds begin @filter(hydro == "permanent") end
temp_ = @chain ponds begin @filter(hydro == "temporary") end
UnequalVarianceTTest(perm.richness, temp_.richness)Two sample t-test (unequal variance)
------------------------------------
Population details:
parameter of interest: Mean difference
value under h_0: 0
point estimate: 1.01849
95% confidence interval: (-0.452, 2.489)
Test summary:
outcome with 95% confidence: fail to reject h_0
two-sided p-value: 0.1734
Details:
number of observations: [95,78]
t-statistic: 1.3673234271337709
degrees of freedom: 168.03764735213443
empirical standard error: 0.7448775533304762In R
t.test(richness ~ hydro, data = ponds) uses Welch by default — the same as UnequalVarianceTTest. For the equal-variance test, use EqualVarianceTTest.
Compare body_size between permanent and temporary ponds with a Welch t test. Pass its p-value to check("8.1", ...).
Reuse perm and temp_ from the block above. pvalue(test) extracts the p-value.
t_body = UnequalVarianceTTest(perm.body_size, temp_.body_size)
pvalue(t_body)0.10269215805061108check("8.1", pvalue(t_body))✓ Correct!
p ≈ 0.10: no evidence of a difference. Here we can go further than usual, because in the simulation hydroperiod really has no effect on body size. With real data, a non-significant test does not show that there is no effect — only that the data cannot tell.
# normality and correlation
(
ks_normal = ExactOneSampleKSTest(
(ponds.body_size .- mean(ponds.body_size)) ./ std(ponds.body_size),
Normal(0, 1)) |> pvalue,
correlation = CorrelationTest(ponds.canopy, ponds.body_size) |> pvalue,
)(ks_normal = 0.9802464726099859, correlation = 4.331805547836094e-5)The dataset’s continuous response is body_size, mean body size.
m1 = lm(@formula(body_size ~ temp_c), ponds)StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}
body_size ~ 1 + temp_c
Coefficients:
─────────────────────────────────────────────────────────────────────────
Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
─────────────────────────────────────────────────────────────────────────
(Intercept) 13.1837 1.5861 8.31 <1e-13 10.0528 16.3145
temp_c 0.921736 0.0645072 14.29 <1e-30 0.794403 1.04907
─────────────────────────────────────────────────────────────────────────The formula interface of StatsModels.jl matches R’s: ~, +, *, :, -1 and inline transformations such as log(x) all work as expected.
(
r2 = r2(m1),
adj_r2 = adjr2(m1),
aic = aic(m1),
coefs = coef(m1),
intervals = confint(m1),
)(r2 = 0.544209635943369, adj_r2 = 0.5415441952178917, aic = 801.3969591368274, coefs = [13.18368268706997, 0.9217361109462694], intervals = [10.05281805410435 16.31454732003559; 0.7944032137607225 1.0490690081318161])Fit a simple regression of body_size on canopy and pass both coefficients to check("8.2", ...).
coef(m) returns the intercept and the slope, in that order.
m_canopy = lm(@formula(body_size ~ canopy), ponds)
coef(m_canopy)2-element Vector{Float64}:
38.134435320100685
-0.04966562886720355check("8.2", coef(m_canopy))✓ Correct!
The slope, about −0.050, is further from the true −0.04 than in m2 (next section), because this model leaves temperature out.
m2 = lm(@formula(body_size ~ temp_c + canopy), ponds)StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}
body_size ~ 1 + temp_c + canopy
Coefficients:
───────────────────────────────────────────────────────────────────────────
Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept) 15.761 1.54229 10.22 <1e-18 12.7164 18.8055
temp_c 0.901451 0.059815 15.07 <1e-32 0.783375 1.01953
canopy -0.0423641 0.00777948 -5.45 <1e-06 -0.057721 -0.0270073
───────────────────────────────────────────────────────────────────────────The generating process (see scripts/make_data.jl) used \(38.0 + 0.9\,(\text{temp} - 24.5) - 0.04\,\text{canopy}\). Compare with the estimates: this is the payoff of teaching with simulated data.
DataFrame(
term = coefnames(m2),
estimated = round.(coef(m2), digits = 3),
truth = [38.0 - 0.9 * 24.5, 0.9, -0.04],
)| Row | term | estimated | truth |
|---|---|---|---|
| String | Float64 | Float64 | |
| 1 | (Intercept) | 15.761 | 15.95 |
| 2 | temp_c | 0.901 | 0.9 |
| 3 | canopy | -0.042 | -0.04 |
Fit body_size ~ temp_c + canopy + hydro and compute the 95% confidence interval for the temp_c coefficient. Does it contain the true value, 0.9? Pass the interval (two numbers) to check("8.3", ...).
confint(m) returns one row per coefficient; temp_c is the second row, so confint(m)[2, :].
m_hydro = lm(@formula(body_size ~ temp_c + canopy + hydro), ponds)
ci_temp = confint(m_hydro)[2, :]2-element Vector{Float64}:
0.7785892869884607
1.0167713699307568check("8.3", ci_temp)✓ Correct!
The interval, about [0.78, 1.02], contains 0.9. It would not in about one of every 20 simulated datasets: that is what “95%” means.
m3 = lm(@formula(body_size ~ temp_c + canopy + region), ponds)StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}
body_size ~ 1 + temp_c + canopy + region
Coefficients:
────────────────────────────────────────────────────────────────────────────────
Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
────────────────────────────────────────────────────────────────────────────────
(Intercept) 15.769 1.54634 10.20 <1e-18 12.7162 18.8218
temp_c 0.893539 0.0597492 14.95 <1e-31 0.775583 1.0115
canopy -0.0437397 0.00780336 -5.61 <1e-07 -0.059145 -0.0283344
region: Atlantic 0.664579 0.436092 1.52 0.1294 -0.196347 1.52551
region: Cerrado 0.018937 0.43474 0.04 0.9653 -0.83932 0.877194
────────────────────────────────────────────────────────────────────────────────Julia treats String and CategoricalArray columns as factors automatically, with treatment (dummy) coding and the first level as reference — same as R. To change the reference level or the contrast scheme:
m3b = lm(@formula(body_size ~ temp_c + region), ponds;
contrasts = Dict(:region => DummyCoding(base = "Cerrado")))
coefnames(m3b)4-element Vector{String}:
"(Intercept)"
"temp_c"
"region: Amazon"
"region: Atlantic"Fit body_size ~ temp_c + region with the Atlantic region as the reference level and pass the coefficient names to check("8.4", ...).
Copy the contrasts argument from the block above and change the level.
m_atl = lm(@formula(body_size ~ temp_c + region), ponds;
contrasts = Dict(:region => DummyCoding(base = "Atlantic")))
coefnames(m_atl)4-element Vector{String}:
"(Intercept)"
"temp_c"
"region: Amazon"
"region: Cerrado"check("8.4", coefnames(m_atl))✓ Correct!
The reference level has no coefficient of its own: it is absorbed by the intercept, and the other regions are compared with it. Without contrasts, the reference is the first level in alphabetical order, Amazon.
m4 = lm(@formula(body_size ~ temp_c * hydro), ponds)StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}
body_size ~ 1 + temp_c + hydro + temp_c & hydro
Coefficients:
───────────────────────────────────────────────────────────────────────────────────────
Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
───────────────────────────────────────────────────────────────────────────────────────
(Intercept) 14.8606 2.16225 6.87 <1e-09 10.5921 19.1291
temp_c 0.847037 0.0890498 9.51 <1e-16 0.671244 1.02283
hydro: temporary -3.33192 3.20923 -1.04 0.3006 -9.66726 3.00341
temp_c & hydro: temporary 0.147595 0.130228 1.13 0.2587 -0.109489 0.404679
───────────────────────────────────────────────────────────────────────────────────────models = ["temp" => m1,
"temp+canopy" => m2,
"temp+can+reg" => m3,
"temp*hydro" => m4]
DataFrame(
model = first.(models),
k = [length(coef(m)) for (_, m) in models],
aic = round.([aic(m) for (_, m) in models], digits = 1),
bic = round.([bic(m) for (_, m) in models], digits = 1),
adj_r2 = round.([adjr2(m) for (_, m) in models], digits = 3),
)| Row | model | k | aic | bic | adj_r2 |
|---|---|---|---|---|---|
| String | Int64 | Float64 | Float64 | Float64 | |
| 1 | temp | 2 | 801.4 | 810.9 | 0.542 |
| 2 | temp+canopy | 3 | 775.6 | 788.2 | 0.607 |
| 3 | temp+can+reg | 5 | 776.1 | 795.0 | 0.611 |
| 4 | temp*hydro | 4 | 803.5 | 819.3 | 0.541 |
# F test between nested models, equivalent to anova() in R
ftest(m1.model, m2.model, m3.model)F-test: 3 models fitted on 173 observations
─────────────────────────────────────────────────────────────────────
DOF ΔDOF SSR ΔSSR R² ΔR² F* p(>F)
─────────────────────────────────────────────────────────────────────
[1] 3 1005.3438 0.5442
[2] 4 1 856.0194 -149.3243 0.6119 0.0677 29.6549 <1e-06
[3] 6 2 839.0138 -17.0056 0.6196 0.0077 1.7026 0.1853
─────────────────────────────────────────────────────────────────────res = residuals(m3)
fit_ = predict(m3)
fig = Figure(size = (860, 350))
ax1 = Axis(fig[1, 1], xlabel = "Fitted values", ylabel = "Residuals",
title = "(a) Residuals vs fitted")
scatter!(ax1, fit_, res, color = (:steelblue, 0.6))
hlines!(ax1, 0, color = :grey40, linestyle = :dash)
ax2 = Axis(fig[1, 2], xlabel = "Theoretical quantiles", ylabel = "Sample quantiles",
title = "(b) QQ plot")
# NOTE: Makie exports `qqnorm` as a plot recipe, so the quantiles are
# computed explicitly here rather than taken from a helper.
n = length(res)
theoretical = quantile.(Normal(0, 1), ((1:n) .- 0.5) ./ n)
scatter!(ax2, theoretical, sort(res), color = (:darkorange, 0.6))
ablines!(ax2, 0, std(res), color = :grey40, linestyle = :dash)
ax3 = Axis(fig[1, 3], xlabel = "Residuals", ylabel = "Frequency",
title = "(c) Distribution")
hist!(ax3, res, bins = 20, color = (:seagreen, 0.7))
fig
m3.
R has plot(model), which returns the four diagnostic plots at once, and performance::check_model(), which goes much further. Julia has no consolidated equivalent; the diagnostics above are assembled by hand. This is a genuine gap, and exactly the kind of thing a community package would fix well.
If the residuals of m3 are roughly normal, about 5% of them should lie more than two standard deviations from zero. Compute the proportion that does and pass it to check("8.5", ...).
Divide res (from the diagnostics block) by std(res), then take the mean of abs.(z) .> 2 — the mean of true/false values is a proportion.
z = res ./ std(res)
prop_large = mean(abs.(z) .> 2)0.057803468208092484check("8.5", prop_large)✓ Correct!
About 5.8% — 10 of 173 residuals, close to the 5% expected. This is a numerical counterpart to the QQ plot above.
MixedModels.jl (Bates et al. 2026) is maintained by Douglas Bates, first author of lme4 (Bates et al. 2015), and uses the same formula notation — including (1 | group).
On performance: comparisons between the two depend heavily on the model, the size of the data and the random-effects structure, and I know of no published benchmark I can cite here. If fitting time matters for your problem, measure it on your own case.
mm1 = fit(MixedModel,
@formula(body_size ~ temp_c + canopy + (1 | region)),
ponds)| Est. | SE | z | p | σ_region | |
|---|---|---|---|---|---|
| (Intercept) | 15.7973 | 1.5290 | 10.33 | <1e-24 | 0.1246 |
| temp_c | 0.9002 | 0.0592 | 15.20 | <1e-51 | |
| canopy | -0.0426 | 0.0077 | -5.52 | <1e-07 | |
| Residual | 2.2209 |
(
fixed_effects = fixef(mm1),
variance = VarCorr(mm1),
icc = first(mm1.σs.region)^2 /
(first(mm1.σs.region)^2 + mm1.σ^2),
)(fixed_effects = [15.79726348949173, 0.9002362213796093, -0.042569241640770514], variance = Variance components:
Column Variance Std.Dev.
region (Intercept) 0.015537 0.124649
Residual 4.932481 2.220919
, icc = 0.0031401442761034163)In R
lmer(body_size ~ temp_c + canopy + (1 | region), data = ponds)Identical, down to the (1 | group) notation. MixedModels.jl also supports (1 + x | group), (1 | a) + (1 | b) and nested terms.
Fit a mixed model of body_size on temp_c with a random intercept for region and pass its fixed effects to check("8.6", ...). Then look at VarCorr of your model: how much variation does region account for?
Same structure as mm1, with a shorter formula. fixef(mm) returns the fixed effects.
mm_temp = fit(MixedModel, @formula(body_size ~ temp_c + (1 | region)), ponds;
progress = false)
fixef(mm_temp)2-element Vector{Float64}:
13.183682687070387
0.9217361109462527check("8.6", fixef(mm_temp))
VarCorr(mm_temp)✓ Correct!
| Column | Variance | Std.Dev | |
|---|---|---|---|
| region | (Intercept) | 0.00000 | 0.00000 |
| Residual | 5.81124 | 2.41065 |
The region variance is estimated as exactly zero — a singular fit. The simulation has no region effect on body size, and three regions give very little information about a variance anyway. With a zero variance, the mixed model reduces to lm(@formula(body_size ~ temp_c), ponds): compare the coefficients.
# parametric bootstrap intervals — fast enough to be routine
using Random
boot = parametricbootstrap(MersenneTwister(42), 1000, mm1)
DataFrame(shortestcovint(boot))body_size ~ canopy separately for each region and compare the slopes with those of a single model with a canopy * region interaction. Why do they agree?m4 have a significant interaction? The simulation has no temp_c × hydro interaction: if yours were significant, what would that mean?mm1 with a model that treats region as a fixed effect. With three regions, which would you report, and why?