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)173Como no Capítulo 8, em cada quadro Sua vez você ajusta um modelo e passa para checar um número extraído dele. Comece cada sessão com:
include("exercises/checar.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))
endContagens não são normais, proporções vivem em \([0,1]\), presenças são binárias. O modelo linear gaussiano trata todas como se fossem contínuas e ilimitadas, e falha de formas previsíveis: intervalos de confiança que incluem riquezas negativas, variância mal especificada, resíduos em leque.
O GLM resolve isso com duas escolhas: a distribuição da resposta e a função de ligação que conecta a média ao preditor linear.
| Resposta | Distribuição | Ligação | Em Julia |
|---|---|---|---|
| Contagem | Poisson | log | Poisson(), LogLink() |
| Contagem sobredispersa | Binomial negativa | log | negbin(...) |
| Presença/ausência | Bernoulli/Binomial | logit | Binomial(), LogitLink() |
| Proporção com denominador | Binomial | logit | Binomial(), LogitLink() |
| Contínua positiva assimétrica | Gama | 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
─────────────────────────────────────────────────────────────────────────────O processo gerador usou \(\log\mu = 0.9 + 0.28\log(\text{área}) - 0.011\,\text{canopy} + 0.03\,(\text{temp}-24.5) + \dots\)
DataFrame(
termo = coefnames(g1),
estimado = round.(coef(g1), digits = 4),
verdadeiro = [missing, 0.28, -0.011, 0.03],
)| Row | termo | estimado | verdadeiro |
|---|---|---|---|
| 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 |
Ajuste um GLM Poisson como o g1, incluindo hydro no fim da fórmula. Passe os coeficientes para checar("9.1", ...) e depois compare o AIC dele com o de g1.
Copie o g1 e acrescente + hydro. aic(g) dá o AIC; quanto menor, melhor.
g_hidro = glm(@formula(richness ~ log_area + canopy + temp_c + hydro),
ponds, Poisson(), LogLink())
coef(g_hidro)5-element Vector{Float64}:
-0.2955819444424689
0.3308987450570613
-0.011543093000739657
0.04752717672780732
-0.1895805246826413checar("9.1", coef(g_hidro))
(aic_g1 = aic(g1), aic_hidro = aic(g_hidro))✓ Correto!
(aic_g1 = 952.9021516391434, aic_hidro = 944.1070557908267)O AIC cai cerca de 9, então hydro merece o lugar. A simulação somou 0,22 à média log das lagoas permanentes; o modelo estima −0,19 para as temporárias, o mesmo efeito visto do outro lado.
Um coeficiente \(\beta\) na ligação log significa que um aumento de uma unidade no preditor multiplica a média por \(e^\beta\). Para canopy:
exp(coef(g1)[3]), (exp(coef(g1)[3]) - 1) * 100(0.9885114593513831, -1.148854064861693)Ou seja, cada ponto percentual de dossel a mais reduz a riqueza esperada em cerca de 1%.
Usando o modelo do 9.1, por qual fator a riqueza esperada muda em uma lagoa temporária em comparação com uma permanente, mantido o resto igual? Passe esse fator para checar("9.2", ...).
O coeficiente de hydro: temporary é o quinto. Na ligação log, um coeficiente \(\beta\) multiplica a média por \(e^\beta\).
fator_temp = exp(coef(g_hidro)[5])0.8273060956543652checar("9.2", fator_temp)✓ Correto!
Cerca de 0,83: espera-se que lagoas temporárias tenham cerca de 17% menos espécies que as permanentes.
A Poisson assume \(\text{Var}(Y) = \text{E}(Y)\). Dados ecológicos raramente cooperam.
# Resíduo de Pearson: (y - mu) / sqrt(Var(mu)). Na Poisson, Var(mu) = mu.
mu = predict(g1)
res_pearson = (ponds.richness .- mu) ./ sqrt.(mu)
gl = dof_residual(g1)
phi = sum(res_pearson .^ 2) / gl
(estatistica_dispersao = phi,
diagnostico = phi > 1.5 ? "sobredisperso" : "aceitável")(estatistica_dispersao = 1.952703419456855, diagnostico = "sobredisperso")Se \(\hat\phi \gg 1\), as alternativas são a binomial negativa ou um modelo misto com efeito aleatório em nível de observação.
ponds.csv foi gerado como um processo Poisson puro — não há variância extra-Poisson no código gerador. Ainda assim, \(\hat\phi \approx 1.95\). Por quê?
Porque o modelo gerador incluía um efeito de região que g1 omite. Heterogeneidade não modelada infla a estatística de dispersão exatamente como a sobredispersão genuína, e a estatística não distingue os dois casos.
Isso importa na prática: o reflexo de recorrer à binomial negativa sempre que \(\hat\phi > 1.5\) pode encobrir um preditor ausente. Verifique a estrutura do modelo antes de trocar a distribuição. O GLMM da Seção 9.6 trata da causa real aqui.
O quadro acima afirma que a sobredispersão de g1 vem de um preditor ausente. Teste isso: ajuste um GLM Poisson com log_area, canopy, temp_c, hydro e region, calcule a estatística de dispersão como no bloco de dispersão e passe-a para checar("9.3", ...).
Copie o bloco de dispersão, trocando g1 pelo seu modelo. predict(g) dá as médias ajustadas.
g_completo = glm(@formula(richness ~ log_area + canopy + temp_c + hydro + region),
ponds, Poisson(), LogLink())
mu_completo = predict(g_completo)
phi_completo = sum(((ponds.richness .- mu_completo) ./ sqrt.(mu_completo)) .^ 2) /
dof_residual(g_completo)0.9410041742922786checar("9.3", phi_completo)✓ Correto!
Cerca de 0,94: com os preditores do modelo gerador no lugar, a “sobredispersão” desaparece. Incluir só hydro (9.1) quase não a mudou; region era a peça que faltava.
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(
modelo = ["Poisson", "Binomial negativa"],
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 | modelo | aic | bic | loglik |
|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | |
| 1 | Poisson | 952.9 | 965.5 | -472.5 |
| 2 | Binomial negativa | 911.3 | 927.0 | -450.6 |
Quando o esforço amostral varia entre unidades, modele a taxa, não a contagem bruta, com um 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 |
Isso ajusta riqueza por unidade de área, fixando o coeficiente de log(area_m2) em 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(
termo = coefnames(g4),
coef = round.(coef(g4), digits = 4),
odds_ratio = round.(exp.(coef(g4)), digits = 3),
ic_inf = round.(exp.(confint(g4)[:, 1]), digits = 3),
ic_sup = round.(exp.(confint(g4)[:, 2]), digits = 3),
)| Row | termo | coef | odds_ratio | ic_inf | ic_sup |
|---|---|---|---|---|---|
| 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 |
A razão de chances para um ponto percentual de dossel fica perto de 1, o que é difícil de interpretar. Calcule, a partir de g4, a razão de chances para um aumento de 10 pontos percentuais e passe-a para checar("9.4", ...).
Dez pontos multiplicam as chances por \(e^{10\beta}\), não por \(10\,e^{\beta}\).
rc_10 = exp(10 * coef(g4)[2])1.1881878414800107checar("9.4", rc_10)✓ Correto!
Cerca de 1,19: cada 10 pontos a mais de dossel multiplicam as chances de encontrar a espécie-foco por 1,19, um aumento de 19%.
grade = DataFrame(
canopy = repeat(range(0, 100, length = 100), 2),
temp_c = fill(mean(ponds.temp_c), 200),
hydro = repeat(["permanent", "temporary"], inner = 100),
)
grade.p = predict(g4, grade)
fig = Figure(size = (700, 430))
ax = Axis(fig[1, 1], xlabel = "Cobertura de dossel (%)",
ylabel = "P(ocorrência)", title = "Espécie-foco")
for (h, cor) in zip(["permanent", "temporary"], [:steelblue, :darkorange])
sub = grade[grade.hydro .== h, :]
lines!(ax, sub.canopy, sub.p, color = cor, linewidth = 2.5, label = h)
end
scatter!(ax, ponds.canopy, Float64.(ponds.focal), color = (:grey40, 0.25))
axislegend(ax, position = :lt)
fig
Uma lagoa temporária tem 80% de dossel e água a 26 °C. Segundo o g4, qual a probabilidade de encontrar a espécie-foco nela? Passe-a para checar("9.5", ...).
Monte um DataFrame de uma linha com as colunas canopy, temp_c e hydro e passe-o para predict(g4, ...), como com grade acima.
lagoa_nova = DataFrame(canopy = [80.0], temp_c = [26.0], hydro = ["temporary"])
p_nova = predict(g4, lagoa_nova)1-element Vector{Union{Missing, Float64}}:
0.5246412509540388checar("9.5", p_nova)✓ Correto!
Cerca de 0,52. Em um GLM logístico, predict retorna probabilidades por padrão: a função de ligação é desfeita para você.
gm1 = fit(MixedModel,
@formula(richness ~ log_area + canopy + (1 | region)),
ponds, Poisson())Minimizing 103 Time: 0:00:00 ( 0.97 ms/it) objective: 211.67417092413703 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 |
Em R
glmer(richness ~ log_area + canopy + (1 | region),
data = ponds, family = poisson)A notação de fórmula é a mesma, (1 | grupo) incluído. Para o método de estimação usado e as opções disponíveis, consulte a documentação de MixedModels.jl (Bates et al. 2026): são detalhes que mudam entre versões, e não vale reproduzi-los aqui de memória.
Inclua hydro no gm1 e passe os efeitos fixos do novo modelo para checar("9.6", ...). Compare a estimativa de hydro com a do 9.1.
Acrescente + hydro antes do termo de efeito aleatório. fixef(gm) retorna os efeitos fixos.
gm_hidro = fit(MixedModel,
@formula(richness ~ log_area + canopy + hydro + (1 | region)),
ponds, Poisson(); progress = false)
fixef(gm_hidro)4-element Vector{Float64}:
0.6004949810030696
0.3694782152614212
-0.010471183105793884
-0.15990712087588418checar("9.6", fixef(gm_hidro))✓ Correto!
Este é o ponto em que ser honesto vale mais que ser entusiasta. A tabela abaixo é uma avaliação de agosto de 2026, feita consultando os registros públicos de pacotes; verifique antes de concluir que algo não existe.
| Precisa de | Em R | Em Julia hoje |
|---|---|---|
| Zero-inflação, hurdle | glmmTMB, pscl |
sem equivalente maduro |
| Regressão beta | betareg, glmmTMB |
sem equivalente maduro |
| Modelos aditivos (GAM) | mgcv |
sem equivalente maduro |
| Resíduos quantílicos simulados | DHARMa |
sem equivalente; implementável à mão |
| Médias marginais, contrastes | emmeans |
Effects.jl, cobertura parcial |
| Regressão ordinal | ordinal, MASS::polr |
cobertura parcial |
RCall.jlNão há vergonha em chamar R de dentro de Julia quando o pacote certo só existe lá. O fluxo continua reprodutível, e os dados atravessam a fronteira automaticamente.
using RCall
@rput ponds # envia o DataFrame para o R
R"""
library(mgcv)
m <- gam(richness ~ s(canopy) + s(log_area) + region,
family = poisson, data = ponds)
summary(m)
"""Para o repertório coberto por este capítulo — Poisson, binomial negativa, binomial, gama, offsets, GLMMs — Julia cobre o necessário, e este capítulo o demonstrou com código executado. Fora dele, RCall.jl é a ponte. Reconhecer a fronteira com clareza é mais útil ao leitor do que escondê-la.
g2. A binomial negativa resolveu a sobredispersão? Depois do 9.3, você ainda a escolheria?DHARMa para o g1: simule 250 conjuntos de dados a partir do modelo ajustado e, para cada lagoa, calcule a proporção de valores simulados abaixo do observado. Se o modelo estiver certo, essas proporções são uniformes. Elas são? Repita para o modelo do 9.3.