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)173Nos quadros Sua vez, você ajusta um modelo e passa para checar um número extraído dele — um valor-p, coeficientes, um intervalo. Todos usam as lagoas com dossel medido, como este capítulo. Comece cada sessão com:
include("exercises/checar.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))
endAntes dos modelos, o básico. HypothesisTests.jl cobre o repertório usual.
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.7448775533304762Em R
t.test(richness ~ hydro, data = ponds) usa Welch por padrão — o mesmo que UnequalVarianceTTest. Para o teste de variâncias iguais, use EqualVarianceTTest.
Compare body_size entre lagoas permanentes e temporárias com um teste t de Welch. Passe o valor-p para checar("8.1", ...).
Reaproveite perm e temp_ do bloco acima. pvalue(teste) extrai o valor-p.
t_corpo = UnequalVarianceTTest(perm.body_size, temp_.body_size)
pvalue(t_corpo)0.10269215805061108checar("8.1", pvalue(t_corpo))✓ Correto!
p ≈ 0,10: nenhuma evidência de diferença. Aqui podemos ir além do usual, porque na simulação o hidroperíodo realmente não afeta o tamanho corporal. Com dados reais, um teste não significativo não mostra que o efeito não existe — só que os dados não permitem dizer.
# normalidade e correlação
(
shapiro_like = ExactOneSampleKSTest(
(ponds.body_size .- mean(ponds.body_size)) ./ std(ponds.body_size),
Normal(0, 1)) |> pvalue,
correlacao = CorrelationTest(ponds.canopy, ponds.body_size) |> pvalue,
)(shapiro_like = 0.9802464726099859, correlacao = 4.331805547836094e-5)A resposta contínua do conjunto é body_size, o tamanho corporal médio.
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
─────────────────────────────────────────────────────────────────────────A interface de fórmulas de StatsModels.jl é a mesma do R: ~, +, *, :, -1 e transformações inline como log(x) funcionam como esperado.
(
r2 = r2(m1),
adj_r2 = adjr2(m1),
aic = aic(m1),
coefs = coef(m1),
intervalos = confint(m1),
)(r2 = 0.544209635943369, adj_r2 = 0.5415441952178917, aic = 801.3969591368274, coefs = [13.18368268706997, 0.9217361109462694], intervalos = [10.05281805410435 16.31454732003559; 0.7944032137607225 1.0490690081318161])Ajuste uma regressão simples de body_size em canopy e passe os dois coeficientes para checar("8.2", ...).
coef(m) retorna o intercepto e a inclinação, nessa ordem.
m_dossel = lm(@formula(body_size ~ canopy), ponds)
coef(m_dossel)2-element Vector{Float64}:
38.134435320100685
-0.04966562886720355checar("8.2", coef(m_dossel))✓ Correto!
A inclinação, cerca de −0,050, fica mais longe do verdadeiro −0,04 do que em m2 (próxima seção), porque este modelo deixa a temperatura de fora.
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
───────────────────────────────────────────────────────────────────────────O processo gerador (ver scripts/make_data.jl) usou \(38.0 + 0.9\,(\text{temp} - 24.5) - 0.04\,\text{canopy}\). Compare com as estimativas: essa é a vantagem de ensinar com dados simulados.
DataFrame(
termo = coefnames(m2),
estimado = round.(coef(m2), digits = 3),
verdadeiro = [38.0 - 0.9 * 24.5, 0.9, -0.04],
)| Row | termo | estimado | verdadeiro |
|---|---|---|---|
| String | Float64 | Float64 | |
| 1 | (Intercept) | 15.761 | 15.95 |
| 2 | temp_c | 0.901 | 0.9 |
| 3 | canopy | -0.042 | -0.04 |
Ajuste body_size ~ temp_c + canopy + hydro e calcule o intervalo de confiança de 95% para o coeficiente de temp_c. Ele contém o valor verdadeiro, 0,9? Passe o intervalo (dois números) para checar("8.3", ...).
confint(m) retorna uma linha por coeficiente; temp_c é a segunda, então confint(m)[2, :].
m_hidro = lm(@formula(body_size ~ temp_c + canopy + hydro), ponds)
ic_temp = confint(m_hidro)[2, :]2-element Vector{Float64}:
0.7785892869884607
1.0167713699307568checar("8.3", ic_temp)✓ Correto!
O intervalo, cerca de [0,78; 1,02], contém 0,9. Em cerca de um a cada 20 conjuntos de dados simulados, não conteria: é isso que “95%” significa.
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 trata String e CategoricalArray como fatores automaticamente, com codificação de tratamento (dummy) e o primeiro nível como referência — igual ao R. Para mudar a referência ou o esquema de contrastes:
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"Ajuste body_size ~ temp_c + region com a região Atlantic como nível de referência e passe os nomes dos coeficientes para checar("8.4", ...).
Copie o argumento contrasts do bloco acima e troque o nível.
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"checar("8.4", coefnames(m_atl))✓ Correto!
O nível de referência não tem coeficiente próprio: ele fica absorvido no intercepto, e as outras regiões são comparadas com ele. Sem contrasts, a referência é o primeiro nível em ordem alfabética, 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
───────────────────────────────────────────────────────────────────────────────────────modelos = ["temp" => m1,
"temp+canopy" => m2,
"temp+can+reg" => m3,
"temp*hydro" => m4]
DataFrame(
modelo = first.(modelos),
k = [length(coef(m)) for (_, m) in modelos],
aic = round.([aic(m) for (_, m) in modelos], digits = 1),
bic = round.([bic(m) for (_, m) in modelos], digits = 1),
adj_r2 = round.([adjr2(m) for (_, m) in modelos], digits = 3),
)| Row | modelo | 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 |
# teste F entre modelos aninhados, equivalente a anova() no 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 = "Valores ajustados", ylabel = "Resíduos",
title = "(a) Resíduos vs ajustados")
scatter!(ax1, fit_, res, color = (:steelblue, 0.6))
hlines!(ax1, 0, color = :grey40, linestyle = :dash)
ax2 = Axis(fig[1, 2], xlabel = "Quantis teóricos", ylabel = "Quantis amostrais",
title = "(b) QQ-plot")
# NOTA: Makie exporta `qqnorm` como recipe de gráfico, então os quantis são
# calculados explicitamente aqui em vez de virem de uma função auxiliar.
n = length(res)
teoricos = quantile.(Normal(0, 1), ((1:n) .- 0.5) ./ n)
scatter!(ax2, teoricos, sort(res), color = (:darkorange, 0.6))
ablines!(ax2, 0, std(res), color = :grey40, linestyle = :dash)
ax3 = Axis(fig[1, 3], xlabel = "Resíduos", ylabel = "Frequência",
title = "(c) Distribuição")
hist!(ax3, res, bins = 20, color = (:seagreen, 0.7))
fig
m3.
R tem plot(modelo), que devolve os quatro gráficos de diagnóstico de uma vez, e performance::check_model(), que vai muito além. Julia não tem equivalente consolidado; os diagnósticos acima são montados à mão. É uma lacuna genuína, e o tipo de coisa que um pacote da comunidade resolveria bem.
Se os resíduos de m3 forem aproximadamente normais, cerca de 5% deles devem estar a mais de dois desvios-padrão de zero. Calcule a proporção que está e passe-a para checar("8.5", ...).
Divida res (do bloco de diagnóstico) por std(res) e tire a média de abs.(z) .> 2 — a média de valores true/false é uma proporção.
z = res ./ std(res)
prop_grandes = mean(abs.(z) .> 2)0.057803468208092484checar("8.5", prop_grandes)✓ Correto!
Cerca de 5,8% — 10 de 173 resíduos, perto dos 5% esperados. É a versão numérica do QQ-plot acima.
MixedModels.jl (Bates et al. 2026) é mantido por Douglas Bates, primeiro autor do lme4 (Bates et al. 2015), e usa a mesma notação de fórmula — inclusive (1 | grupo).
Sobre desempenho: comparações entre os dois dependem fortemente do modelo, do tamanho dos dados e da estrutura dos efeitos aleatórios, e não conheço um benchmark publicado que eu possa citar aqui. Se o tempo de ajuste importa para o seu problema, meça no seu próprio caso.
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 |
(
efeitos_fixos = fixef(mm1),
desvios_random = VarCorr(mm1),
icc = first(mm1.σs.region)^2 /
(first(mm1.σs.region)^2 + mm1.σ^2),
)(efeitos_fixos = [15.79726348949173, 0.9002362213796093, -0.042569241640770514], desvios_random = Variance components:
Column Variance Std.Dev.
region (Intercept) 0.015537 0.124649
Residual 4.932481 2.220919
, icc = 0.0031401442761034163)Em R
lmer(body_size ~ temp_c + canopy + (1 | region), data = ponds)Idêntico, inclusive a notação (1 | grupo). MixedModels.jl também suporta (1 + x | grupo), (1 | a) + (1 | b) e termos aninhados.
Ajuste um modelo misto de body_size em temp_c com intercepto aleatório para region e passe os efeitos fixos para checar("8.6", ...). Depois olhe o VarCorr do seu modelo: quanto da variação a região explica?
Mesma estrutura de mm1, com uma fórmula menor. fixef(mm) retorna os efeitos fixos.
mm_temp = fit(MixedModel, @formula(body_size ~ temp_c + (1 | region)), ponds;
progress = false)
fixef(mm_temp)2-element Vector{Float64}:
13.183682687070387
0.9217361109462527checar("8.6", fixef(mm_temp))
VarCorr(mm_temp)✓ Correto!
| Column | Variance | Std.Dev | |
|---|---|---|---|
| region | (Intercept) | 0.00000 | 0.00000 |
| Residual | 5.81124 | 2.41065 |
A variância de região é estimada como exatamente zero — um ajuste singular. A simulação não tem efeito de região sobre o tamanho corporal, e três regiões trazem muito pouca informação sobre uma variância de qualquer forma. Com variância zero, o modelo misto se reduz a lm(@formula(body_size ~ temp_c), ponds): compare os coeficientes.
# intervalos por bootstrap paramétrico — rápido o bastante para ser rotina
using Random
boot = parametricbootstrap(MersenneTwister(42), 1000, mm1)
DataFrame(shortestcovint(boot))body_size ~ canopy separadamente para cada região e compare as inclinações com as de um único modelo com interação canopy * region. Por que elas coincidem?m4 tem interação significativa? A simulação não tem interação temp_c × hydro: se a sua fosse significativa, o que isso significaria?mm1 com um modelo que trata region como efeito fixo. Com três regiões, qual você reportaria, e por quê?