using TidierData
using DataFrames, CSV
using GLM, StatsModels, StatsBase, Statistics, Distributions
using MixedModels
using CairoMakie
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)173As in Chapter 8, you fit a model in each Your turn box and pass a number extracted from it to check. Start each session with:
include("exercises/check.jl")
using TidierData, DataFrames, CSV, GLM, StatsModels, Statistics, MixedModels
ponds = @chain CSV.read(joinpath("data", "ponds.csv"), DataFrame) begin
@drop_missing(canopy)
@mutate(log_area = log(area_m2))
endCounts are not normal, proportions live in \([0,1]\), presences are binary. The Gaussian linear model treats all of them as if they were continuous and unbounded, and fails in predictable ways: confidence intervals that include negative richness, misspecified variance, fan-shaped residuals.
The GLM fixes this with two choices: the distribution of the response and the link function connecting the mean to the linear predictor.
| Response | Distribution | Link | In Julia |
|---|---|---|---|
| Count | Poisson | log | Poisson(), LogLink() |
| Overdispersed count | Negative binomial | log | negbin(...) |
| Presence/absence | Bernoulli/Binomial | logit | Binomial(), LogitLink() |
| Proportion with denominator | Binomial | logit | Binomial(), LogitLink() |
| Positive skewed continuous | Gamma | log | Gamma(), LogLink() |
g1 = glm(@formula(richness ~ log_area + canopy + temp_c),
ponds, Poisson(), LogLink())StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Vector{Float64}, Poisson{Float64}, LogLink}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}
richness ~ 1 + log_area + canopy + temp_c
Coefficients:
─────────────────────────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|) Lower 95% Upper 95%
─────────────────────────────────────────────────────────────────────────────
(Intercept) -0.30001 0.289756 -1.04 0.3005 -0.867921 0.267902
log_area 0.331763 0.0300907 11.03 <1e-27 0.272786 0.390739
canopy -0.011555 0.0013362 -8.65 <1e-17 -0.0141739 -0.00893614
temp_c 0.0441958 0.00959497 4.61 <1e-05 0.02539 0.0630016
─────────────────────────────────────────────────────────────────────────────The generating process used \(\log\mu = 0.9 + 0.28\log(\text{area}) - 0.011\,\text{canopy} + 0.03\,(\text{temp}-24.5) + \dots\)
DataFrame(
term = coefnames(g1),
estimated = round.(coef(g1), digits = 4),
truth = [missing, 0.28, -0.011, 0.03],
)| Row | term | estimated | truth |
|---|---|---|---|
| String | Float64 | Float64? | |
| 1 | (Intercept) | -0.3 | missing |
| 2 | log_area | 0.3318 | 0.28 |
| 3 | canopy | -0.0116 | -0.011 |
| 4 | temp_c | 0.0442 | 0.03 |
Fit a Poisson GLM like g1, adding hydro at the end of the formula. Pass its coefficients to check("9.1", ...), then compare its AIC with that of g1.
Copy g1 and add + hydro. aic(g) gives the AIC; lower is better.
g_hydro = glm(@formula(richness ~ log_area + canopy + temp_c + hydro),
ponds, Poisson(), LogLink())
coef(g_hydro)5-element Vector{Float64}:
-0.2955819444424689
0.3308987450570613
-0.011543093000739657
0.04752717672780732
-0.1895805246826413check("9.1", coef(g_hydro))
(aic_g1 = aic(g1), aic_hydro = aic(g_hydro))✓ Correct!
(aic_g1 = 952.9021516391434, aic_hydro = 944.1070557908267)The AIC drops by about 9, so hydro earns its place. The simulation added 0.22 to the log mean of permanent ponds; the model estimates −0.19 for temporary ones, the same effect seen from the other side.
A coefficient \(\beta\) under a log link means that a one-unit increase in the predictor multiplies the mean by \(e^\beta\). For canopy:
exp(coef(g1)[3]), (exp(coef(g1)[3]) - 1) * 100(0.9885114593513831, -1.148854064861693)That is, each additional percentage point of canopy reduces expected richness by about 1%.
Using the model from 9.1, by what factor does expected richness change in a temporary pond compared with a permanent one, all else equal? Pass that factor to check("9.2", ...).
The hydro: temporary coefficient is the fifth. Under a log link, a coefficient \(\beta\) multiplies the mean by \(e^\beta\).
factor_temp = exp(coef(g_hydro)[5])0.8273060956543652check("9.2", factor_temp)✓ Correct!
About 0.83: temporary ponds are expected to hold about 17% fewer species than permanent ones.
The Poisson assumes \(\text{Var}(Y) = \text{E}(Y)\). Ecological data rarely cooperate.
# Pearson residual: (y - mu) / sqrt(Var(mu)). For the Poisson, Var(mu) = mu.
mu = predict(g1)
pearson_res = (ponds.richness .- mu) ./ sqrt.(mu)
df_res = dof_residual(g1)
phi = sum(pearson_res .^ 2) / df_res
(dispersion_statistic = phi,
verdict = phi > 1.5 ? "overdispersed" : "acceptable")(dispersion_statistic = 1.952703419456855, verdict = "overdispersed")If \(\hat\phi \gg 1\), the alternatives are the negative binomial or a mixed model with an observation-level random effect.
ponds.csv was generated as a pure Poisson process — there is no extra-Poisson variance in the data-generating code. Yet \(\hat\phi \approx 1.95\). Why?
Because the generating model included a region effect that g1 omits. Unmodelled heterogeneity inflates the dispersion statistic exactly as genuine overdispersion does, and the statistic cannot tell you which one you have.
This matters practically: the reflex of reaching for a negative binomial whenever \(\hat\phi > 1.5\) can paper over a missing predictor. Check your model structure before you change your distribution. The GLMM in Section 9.6 addresses the actual cause here.
The box above claims that the overdispersion of g1 comes from a missing predictor. Test it: fit a Poisson GLM with log_area, canopy, temp_c, hydro and region, compute its dispersion statistic as in the block above, and pass it to check("9.3", ...).
Copy the dispersion block, replacing g1 with your model. predict(g) gives the fitted means.
g_full = glm(@formula(richness ~ log_area + canopy + temp_c + hydro + region),
ponds, Poisson(), LogLink())
mu_full = predict(g_full)
phi_full = sum(((ponds.richness .- mu_full) ./ sqrt.(mu_full)) .^ 2) /
dof_residual(g_full)0.9410041742922786check("9.3", phi_full)✓ Correct!
About 0.94: with the predictors of the generating model in place, the “overdispersion” is gone. Adding hydro alone (9.1) barely moved it; region was the missing piece.
g2 = negbin(@formula(richness ~ log_area + canopy + temp_c), ponds, LogLink())StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Vector{Float64}, NegativeBinomial{Float64}, LogLink}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}
richness ~ 1 + log_area + canopy + temp_c
Coefficients:
─────────────────────────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|) Lower 95% Upper 95%
─────────────────────────────────────────────────────────────────────────────
(Intercept) -0.165087 0.424739 -0.39 0.6975 -0.997561 0.667387
log_area 0.334165 0.0432899 7.72 <1e-13 0.249319 0.419012
canopy -0.0120529 0.0019084 -6.32 <1e-09 -0.0157933 -0.00831254
temp_c 0.0391606 0.0141787 2.76 0.0057 0.0113708 0.0669503
─────────────────────────────────────────────────────────────────────────────DataFrame(
model = ["Poisson", "Negative binomial"],
aic = round.([aic(g1), aic(g2)], digits = 1),
bic = round.([bic(g1), bic(g2)], digits = 1),
loglik = round.([loglikelihood(g1), loglikelihood(g2)], digits = 1),
)| Row | model | aic | bic | loglik |
|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | |
| 1 | Poisson | 952.9 | 965.5 | -472.5 |
| 2 | Negative binomial | 911.3 | 927.0 | -450.6 |
When sampling effort varies across units, model the rate rather than the raw count, using an offset.
g3 = glm(@formula(richness ~ canopy + temp_c), ponds, Poisson(), LogLink();
offset = log.(ponds.area_m2))
coeftable(g3)| Coef. | Std. Error | z | Pr(>|z|) | Lower 95% | Upper 95% | |
|---|---|---|---|---|---|---|
| (Intercept) | -3.56359 | 0.244593 | -14.57 | <1e-47 | -4.04299 | -3.0842 |
| canopy | -0.0151224 | 0.001304 | -11.60 | <1e-30 | -0.0176782 | -0.0125666 |
| temp_c | 0.0310744 | 0.00949171 | 3.27 | 0.0011 | 0.012471 | 0.0496778 |
This models richness per unit area, fixing the coefficient of log(area_m2) at 1.
g4 = glm(@formula(focal ~ canopy + temp_c + hydro),
ponds, Binomial(), LogitLink())StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Vector{Float64}, Binomial{Float64}, LogitLink}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}
focal ~ 1 + canopy + temp_c + hydro
Coefficients:
────────────────────────────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|) Lower 95% Upper 95%
────────────────────────────────────────────────────────────────────────────────
(Intercept) -5.79176 1.6009 -3.62 0.0003 -8.92947 -2.65406
canopy 0.0172429 0.00759302 2.27 0.0232 0.0023609 0.032125
temp_c 0.209118 0.0616577 3.39 0.0007 0.088271 0.329965
hydro: temporary -0.926093 0.338929 -2.73 0.0063 -1.59038 -0.261805
────────────────────────────────────────────────────────────────────────────────DataFrame(
term = coefnames(g4),
coef = round.(coef(g4), digits = 4),
odds_ratio = round.(exp.(coef(g4)), digits = 3),
ci_lower = round.(exp.(confint(g4)[:, 1]), digits = 3),
ci_upper = round.(exp.(confint(g4)[:, 2]), digits = 3),
)| Row | term | coef | odds_ratio | ci_lower | ci_upper |
|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Float64 | |
| 1 | (Intercept) | -5.7918 | 0.003 | 0.0 | 0.07 |
| 2 | canopy | 0.0172 | 1.017 | 1.002 | 1.033 |
| 3 | temp_c | 0.2091 | 1.233 | 1.092 | 1.391 |
| 4 | hydro: temporary | -0.9261 | 0.396 | 0.204 | 0.77 |
The odds ratio for one percentage point of canopy is close to 1, which is hard to interpret. Compute the odds ratio for an increase of 10 percentage points, from g4, and pass it to check("9.4", ...).
Ten points multiply the odds by \(e^{10\beta}\), not by \(10\,e^{\beta}\).
or_10 = exp(10 * coef(g4)[2])1.1881878414800107check("9.4", or_10)✓ Correct!
About 1.19: each extra 10 points of canopy multiply the odds of finding the focal species by 1.19, an increase of 19%.
grid = DataFrame(
canopy = repeat(range(0, 100, length = 100), 2),
temp_c = fill(mean(ponds.temp_c), 200),
hydro = repeat(["permanent", "temporary"], inner = 100),
)
grid.p = predict(g4, grid)
fig = Figure(size = (700, 430))
ax = Axis(fig[1, 1], xlabel = "Canopy cover (%)",
ylabel = "P(occurrence)", title = "Focal species")
for (h, col) in zip(["permanent", "temporary"], [:steelblue, :darkorange])
sub = grid[grid.hydro .== h, :]
lines!(ax, sub.canopy, sub.p, color = col, linewidth = 2.5, label = h)
end
scatter!(ax, ponds.canopy, Float64.(ponds.focal), color = (:grey40, 0.25))
axislegend(ax, position = :lt)
fig
A temporary pond has 80% canopy and water at 26 °C. According to g4, what is the probability of finding the focal species there? Pass it to check("9.5", ...).
Build a one-row DataFrame with the columns canopy, temp_c and hydro, and pass it to predict(g4, ...), as with grid above.
new_pond = DataFrame(canopy = [80.0], temp_c = [26.0], hydro = ["temporary"])
p_new = predict(g4, new_pond)1-element Vector{Union{Missing, Float64}}:
0.5246412509540388check("9.5", p_new)✓ Correct!
About 0.52. predict returns probabilities for a logistic GLM by default: the link function is undone for you.
gm1 = fit(MixedModel,
@formula(richness ~ log_area + canopy + (1 | region)),
ponds, Poisson())Minimizing 133 Time: 0:00:00 ( 0.75 ms/it) objective: 211.6522702865303 Minimizing 281 Time: 0:00:00 ( 0.57 ms/it)
| Est. | SE | z | p | σ_region | |
|---|---|---|---|---|---|
| (Intercept) | 0.5193 | 0.2688 | 1.93 | 0.0533 | 0.3534 |
| log_area | 0.3702 | 0.0307 | 12.05 | <1e-32 | |
| canopy | -0.0103 | 0.0013 | -7.72 | <1e-13 |
In R
glmer(richness ~ log_area + canopy + (1 | region),
data = ponds, family = poisson)The formula notation is the same, (1 | group) included. For the estimation method used and the options available, consult the MixedModels.jl documentation (Bates et al. 2026): these are details that change between versions, and not worth reproducing here from memory.
Add hydro to gm1 and pass the fixed effects of the new model to check("9.6", ...). Compare the hydro estimate with the one from 9.1.
Add + hydro before the random-effect term. fixef(gm) returns the fixed effects.
gm_hydro = fit(MixedModel,
@formula(richness ~ log_area + canopy + hydro + (1 | region)),
ponds, Poisson(); progress = false)
fixef(gm_hydro)4-element Vector{Float64}:
0.6004949810030696
0.3694782152614212
-0.010471183105793884
-0.15990712087588418check("9.6", fixef(gm_hydro))✓ Correct!
This is the point where being honest is worth more than being enthusiastic. The table below is an assessment as of August 2026, made by checking the public package registries; check before concluding that something does not exist.
| You need | In R | In Julia today |
|---|---|---|
| Zero-inflation, hurdle | glmmTMB, pscl |
no mature equivalent |
| Beta regression | betareg, glmmTMB |
no mature equivalent |
| Additive models (GAM) | mgcv |
no mature equivalent |
| Simulated quantile residuals | DHARMa |
none; implementable by hand |
| Marginal means, contrasts | emmeans |
Effects.jl, partial coverage |
| Ordinal regression | ordinal, MASS::polr |
partial coverage |
RCall.jlThere is no shame in calling R from inside Julia when the right package only exists there. The workflow stays reproducible, and data crosses the boundary automatically.
using RCall
@rput ponds # send the DataFrame to R
R"""
library(mgcv)
m <- gam(richness ~ s(canopy) + s(log_area) + region,
family = poisson, data = ponds)
summary(m)
"""For the repertoire covered by this chapter — Poisson, negative binomial, binomial, gamma, offsets, GLMMs — Julia covers what you need, and this chapter has demonstrated it with executed code. Outside it, RCall.jl is the bridge. Naming the boundary clearly serves the reader better than hiding it.
g2. Did the negative binomial fix the overdispersion? After 9.3, would you still choose it?DHARMa-style quantile residuals for g1: simulate 250 datasets from the fitted model and, for each pond, compute the proportion of simulated values below the observed one. If the model is right, these proportions are uniform. Are they? Repeat for the model from 9.3.