Todos os testes do Capítulo 8 e do Capítulo 9 se apoiam em uma distribuição teórica: a t, a normal, a \(\chi^2\). Este capítulo constrói essas distribuições a partir dos próprios dados, embaralhando e reamostrando, e depois inverte a ideia: simula dados de um processo conhecido para perguntar com que frequência um estudo de certo tamanho encontraria um efeito.
Julia combina bem com isso. Reamostrar significa rodar o mesmo cálculo milhares de vezes, e em Julia o loop simples que expressa a ideia também é rápido (Seção 4.7); não é preciso contorcer o código para evitá-lo.
Como o Capítulo 10, este capítulo é inspirado no ecoR (Oliveira 2026), o wiki do curso de R do Instituto de Biociências da Universidade de São Paulo, produzido e mantido por Alexandre Adalardo de Oliveira. A sequência — permutação, depois bootstrap, depois simulação a partir de distribuições conhecidas, com um modelo nulo de coocorrência de espécies — segue o seu módulo de reamostragem. O texto, os exemplos e os exercícios daqui foram escritos do zero para Julia; os autores do ecoR não revisaram este livro.
Se quiser ir além, o material original está disponível gratuitamente (CC BY-SA 4.0):
using TidierData
using DataFrames, CSV
using Statistics, StatsBase, Random, Distributions
using GLM, HypothesisTests
using CairoMakie
CairoMakie.activate!(type = "png")
Random.seed!(2026)
ponds = @chain CSV.read(joinpath("..", "data", "ponds.csv"), DataFrame) begin
@drop_missing(canopy)
@mutate(log_area = log(area_m2))
end
nrow(ponds)173A maioria das respostas aqui vem de sorteios aleatórios; por isso a sua não vai bater com a da solução até a última casa, e checar aceita qualquer valor dentro da variação aleatória esperada. Use pelo menos tantas permutações ou simulações quanto o exercício pede: com poucas, a resposta fica ruidosa demais para passar. Comece cada sessão com:
include("exercises/checar.jl")
using TidierData, DataFrames, CSV, Statistics, StatsBase, Random, Distributions
using GLM, HypothesisTests
ponds = @chain CSV.read(joinpath("data", "ponds.csv"), DataFrame) begin
@drop_missing(canopy)
@mutate(log_area = log(area_m2))
endTrês funções fazem quase todo o trabalho deste capítulo. shuffle reordena um vetor; sample, do StatsBase.jl, sorteia elementos dele com ou sem reposição:
x = [1, 2, 3, 4, 5]
(embaralhado = shuffle(x),
sem_reposicao = sample(x, 3; replace = false),
com_reposicao = sample(x, 5)) # replace = true é o padrão(embaralhado = [1, 3, 4, 2, 5], sem_reposicao = [2, 4, 5], com_reposicao = [2, 5, 2, 5, 4])Os números aleatórios em Julia vêm de um gerador de números aleatórios (RNG). Fixar a semente, como o bloco inicial faz com Random.seed!(2026), faz um script produzir os mesmos sorteios toda vez que roda. Para ter mais controle, crie o seu próprio gerador e passe-o como primeiro argumento:
rng = Xoshiro(42)
sample(rng, x, 3)3-element Vector{Int64}:
4
3
3O manual do Julia avisa que uma versão minor nova (da 1.11 para a 1.12, por exemplo) pode mudar a sequência de números aleatórios gerada a partir de uma semente. A semente torna a sua análise reprodutível na sua instalação, não em todas as versões. Para reprodutibilidade exata, salve os dados aleatórios, ou fixe a versão do Julia e o Manifest.toml. É também por isso que os exercícios daqui aceitam uma faixa de respostas.
Na Seção 8.1, um teste t de Welch perguntou se o tamanho corporal difere entre lagoas permanentes e temporárias. Um teste de permutação faz a mesma pergunta sem supor nenhuma distribuição, em cinco passos:
eh_perm = ponds.hydro .== "permanent"
dif_medias(y, g) = mean(y[g]) - mean(y[.!g])
observado = dif_medias(ponds.body_size, eh_perm)
B = 9_999
difs_nulas = zeros(B)
rotulos = copy(eh_perm)
for b in 1:B
shuffle!(rotulos) # passos 3–4: embaralha, recalcula
difs_nulas[b] = dif_medias(ponds.body_size, rotulos)
end
extremos = count(abs.(difs_nulas) .>= abs(observado))
p_perm = (extremos + 1) / (B + 1)
(observado = observado, p_perm = p_perm)(observado = -0.8929446693657184, p_perm = 0.1073)Dois detalhes. abs torna o teste bilateral: uma diferença de −0,9 conta como tão extrema quanto uma de +0,9. E o + 1 no numerador e no denominador conta os dados observados como um dos arranjos possíveis, o que impede que o valor-p seja exatamente zero (Phipson e Smyth 2010).
fig = Figure(size = (700, 380))
ax = Axis(fig[1, 1], xlabel = "Diferença no tamanho corporal médio (mm), permanente − temporária",
ylabel = "Embaralhamentos")
hist!(ax, difs_nulas, bins = 50, color = (:grey60, 0.8))
vlines!(ax, [observado, -observado], color = :darkorange, linestyle = :dash, linewidth = 2)
fig
O valor-p, cerca de 0,10, fica perto do 0,103 do teste de Welch: com 173 lagoas e dados aproximadamente simétricos, a teoria e o embaralhamento concordam. Isso nem sempre vale com amostras pequenas ou dados assimétricos, e é aí que a permutação mostra seu valor. O HypothesisTests.jl empacota a mesma ideia:
lagoas_perm = ponds.body_size[eh_perm]
lagoas_temp = ponds.body_size[.!eh_perm]
pvalue(ApproximatePermutationTest(lagoas_perm, lagoas_temp, mean, 9_999))0.10461046104610461A cobertura de dossel é diferente nas lagoas onde a espécie-foco ocorre? Faça um teste de permutação bilateral para a diferença na média de canopy entre lagoas com focal == 1 e focal == 0, com 9 999 embaralhamentos, e passe o valor-p para checar("11.1", ...).
Reaproveite dif_medias e o loop acima, com ponds.canopy como resposta e ponds.focal .== 1 como grupos.
tem_focal = ponds.focal .== 1
obs_dossel = dif_medias(ponds.canopy, tem_focal)
grupos = copy(tem_focal)
nulo_dossel = [dif_medias(ponds.canopy, shuffle!(grupos)) for _ in 1:9_999]
p_dossel = (count(abs.(nulo_dossel) .>= abs(obs_dossel)) + 1) / (9_999 + 1)0.0353checar("11.1", p_dossel)✓ Correto!
Cerca de 0,04: lagoas com a espécie-foco têm cerca de 7 pontos a mais de dossel. A simulação pôs exatamente isso: a ocorrência aumenta com o dossel (scripts/make_data.jl). A compreensão faz o mesmo que o loop, em uma linha.
A mesma lógica testa uma inclinação. Se o tamanho corporal não depende do dossel, qualquer pareamento de tamanhos corporais com valores de dossel é tão provável quanto o observado; então embaralhar a resposta quebra a relação e mantém intactas as duas distribuições:
inclinacao(x, y) = cov(x, y) / var(x)
obs_inclinacao = inclinacao(ponds.canopy, ponds.body_size)
y = copy(ponds.body_size)
nulo_inclinacao = [inclinacao(ponds.canopy, shuffle!(y)) for _ in 1:9_999]
extremos = count(abs.(nulo_inclinacao) .>= abs(obs_inclinacao))
(inclinacao = obs_inclinacao, extremos = extremos, p = (extremos + 1) / (9_999 + 1))(inclinacao = -0.04966562886720383, extremos = 0, p = 0.0001)cov(x, y) / var(x) é a inclinação de mínimos quadrados, o mesmo número que o lm estima na Seção 8.2. Quase nenhum embaralhamento produz uma inclinação tão acentuada quanto a observada — a contagem está em extremos —, então o valor-p fica no menor valor que 9 999 embaralhamentos conseguem mostrar, 1/10 000, ou perto dele. Para dizer mais, seriam necessários mais embaralhamentos, mas “cerca de 1 em 10 000 ou menos” costuma ser tudo o que um resultado precisa.
Teste por permutação se o tamanho corporal está relacionado à área da lagoa: use a inclinação de body_size em log_area, bilateral, com 9 999 embaralhamentos. Passe o valor-p para checar("11.2", ...).
O bloco acima com ponds.log_area no lugar de ponds.canopy. Embaralhe uma cópia de body_size, não a própria coluna, para manter os dados intactos.
obs_area = inclinacao(ponds.log_area, ponds.body_size)
y = copy(ponds.body_size)
nulo_area = [inclinacao(ponds.log_area, shuffle!(y)) for _ in 1:9_999]
p_area = (count(abs.(nulo_area) .>= abs(obs_area)) + 1) / (9_999 + 1)0.2583checar("11.2", p_area)✓ Correto!
Cerca de 0,26. A simulação não tem efeito da área sobre o tamanho corporal, e o teste não encontra nenhum — corretamente, desta vez. Um teste em cada vinte encontraria um mesmo assim.
Um teste de permutação pergunta como os dados seriam se não houvesse efeito. O bootstrap (Efron 1979) faz outra pergunta: quanto uma estimativa mudaria se tivéssemos coletado outra amostra? Não podemos coletá-la, mas podemos imitá-la sorteando uma amostra nova do mesmo tamanho a partir dos próprios dados, com reposição. Algumas lagoas aparecem duas vezes, outras nenhuma, e a dispersão da estimativa ao longo de muitas dessas amostras mede a sua incerteza.
n = nrow(ponds)
medianas_boot = [median(sample(ponds.body_size, n)) for _ in 1:9_999]
(mediana = median(ponds.body_size),
ic_95 = quantile(medianas_boot, [0.025, 0.975]))(mediana = 35.56, ic_95 = [35.08, 36.54])Os quantis de 2,5% e 97,5% das medianas bootstrap dão um intervalo percentil de 95% — para uma estatística, a mediana, que não tem uma fórmula simples para o seu erro-padrão.
O bootstrap também confere a incerteza que o próprio modelo reporta. Reamostre lagoas inteiras (linhas), reajuste o modelo Poisson g1 da Seção 9.2 a cada vez e compare a dispersão do coeficiente de dossel com o erro-padrão que o modelo informa:
g1 = glm(@formula(richness ~ log_area + canopy + temp_c), ponds, Poisson(), LogLink())
dossel_boot = map(1:2_000) do _
linhas = sample(1:n, n)
coef(glm(@formula(richness ~ log_area + canopy + temp_c),
ponds[linhas, :], Poisson(), LogLink()))[3]
end
(ep_bootstrap = std(dossel_boot),
ep_modelo = stderror(g1)[3],
razao = std(dossel_boot) / stderror(g1)[3])(ep_bootstrap = 0.0018378945669366084, ep_modelo = 0.0013361976263705473, razao = 1.375466121676026)O erro-padrão bootstrap é cerca de 1,4 vez maior que o que o modelo Poisson informa — perto de \(\sqrt{1{,}95}\), a raiz da estatística de dispersão de g1 na Seção 9.3. Os erros-padrão do modelo supõem que a variância é igual à média; o bootstrap não supõe nada sobre ela, e enxerga a variação extra que o modelo deixa de fora.
do
map(1:2_000) do _ ... end é map(_ -> ..., 1:2_000) escrito como bloco, para funções anônimas com mais de uma linha. _ nomeia um argumento que você não usa.
Escreva boot_means(x, B), que retorna um vetor com as médias de B amostras bootstrap de x. Verifique com checar("11.3", boot_means).
Uma compreensão: [mean(sample(x, length(x))) for _ in 1:B]. Cada amostra bootstrap tem o mesmo tamanho de x e é sorteada com reposição.
boot_means(x, B) = [mean(sample(x, length(x))) for _ in 1:B]
checar("11.3", boot_means)✓ Correto!
Usando boot_means com 9 999 amostras, calcule um intervalo percentil de 95% para a riqueza média das lagoas da Amazônia. Passe-o como [inferior, superior] para checar("11.4", ...). Compare-o com o intervalo t usual.
ponds.richness[ponds.region .== "Amazon"] seleciona as lagoas da Amazônia. quantile(v, [0.025, 0.975]) retorna os dois extremos de uma vez.
amazonia = ponds.richness[ponds.region .== "Amazon"]
ic_boot = quantile(boot_means(amazonia, 9_999), [0.025, 0.975])2-element Vector{Float64}:
9.478260869565217
12.956521739130435checar("11.4", ic_boot)
confint(OneSampleTTest(amazonia))✓ Correto!
(9.384117221613428, 12.963708865343095)Os dois intervalos quase coincidem, de cerca de 9,4–9,5 a 13. O bootstrap começa um pouco mais alto: a riqueza é uma contagem assimétrica, e o bootstrap acompanha a assimetria, enquanto o intervalo t é simétrico por construção.
A reamostragem reaproveita os dados que você tem. A simulação cria dados que você ainda não tem, a partir de um processo que você especifica — e, como foi você que escolheu o processo, sabe a resposta certa. Essa é a lógica de todo este livro (Seção 5.1), e ela responde a uma pergunta que todo projeto de campo enfrenta antes de começar: quantas unidades amostrais bastam?
Suponha que lagoas permanentes tenham cerca de 25% mais espécies que as temporárias, como na simulação por trás de ponds.csv (\(e^{0{,}22} \approx 1{,}25\)). Esta função gera um levantamento de n lagoas a partir desse processo:
function simulate_ponds(n; p_temporary = 0.5, rng = Random.default_rng())
log_area = rand(rng, Normal(5.0, 0.9), n)
hydro = ifelse.(rand(rng, n) .< p_temporary, "temporary", "permanent")
mu = @. exp(0.9 + 0.28 * log_area + 0.22 * (hydro == "permanent"))
richness = rand.(rng, Poisson.(mu))
return DataFrame(; log_area, hydro, richness)
end
first(simulate_ponds(5), 5)| Row | log_area | hydro | richness |
|---|---|---|---|
| Float64 | String | Int64 | |
| 1 | 4.34568 | temporary | 13 |
| 2 | 6.30406 | permanent | 18 |
| 3 | 5.49257 | temporary | 12 |
| 4 | 4.17274 | permanent | 7 |
| 5 | 4.56875 | permanent | 9 |
Para um levantamento de n lagoas, ajuste o modelo e registre se ele detecta o efeito do hidroperíodo ao nível de 5%. A proporção de levantamentos simulados que detectam é o poder desse desenho:
function detects_hydro(n; kwargs...)
d = simulate_ponds(n; kwargs...)
m = glm(@formula(richness ~ log_area + hydro), d, Poisson(), LogLink())
return coeftable(m).cols[4][3] < 0.05 # valor-p do coeficiente de hydro
end
tamanhos = [20, 40, 60, 80, 100, 140, 180]
poder = [mean(detects_hydro(n) for _ in 1:500) for n in tamanhos]
DataFrame(n = tamanhos, poder = poder)| Row | n | poder |
|---|---|---|
| Int64 | Float64 | |
| 1 | 20 | 0.35 |
| 2 | 40 | 0.616 |
| 3 | 60 | 0.796 |
| 4 | 80 | 0.886 |
| 5 | 100 | 0.954 |
| 6 | 140 | 0.992 |
| 7 | 180 | 1.0 |
fig = Figure(size = (700, 380))
ax = Axis(fig[1, 1], xlabel = "Número de lagoas", ylabel = "Poder")
lines!(ax, tamanhos, poder, color = :steelblue, linewidth = 2.5)
scatter!(ax, tamanhos, poder, color = :steelblue)
hlines!(ax, 0.8, color = :grey40, linestyle = :dash)
ylims!(ax, 0, 1)
fig
Cerca de 60 lagoas atingem os 80% de poder usuais; com 20, dois levantamentos em cada três deixariam de ver um efeito que existe de verdade. E um resultado não significativo de um levantamento pequeno diz pouco: a maioria deles é não significativa exista ou não o efeito.
No campo, raramente se consegue o mesmo número de cada tipo de lagoa. Estime o poder com 60 lagoas quando só 20% são temporárias (p_temporary = 0.2), a partir de pelo menos 1 000 levantamentos simulados, e passe-o para checar("11.5", ...).
detects_hydro repassa seus argumentos nomeados para simulate_ponds (Seção 10.2), então detects_hydro(60; p_temporary = 0.2) já faz o que você precisa.
poder_desbalanceado = mean(detects_hydro(60; p_temporary = 0.2) for _ in 1:2_000)0.5885checar("11.5", poder_desbalanceado)✓ Correto!
Cerca de 0,60, abaixo dos 0,80 do desenho balanceado: as mesmas 60 lagoas dizem menos, porque a comparação se apoia em só umas 12 lagoas temporárias.
A permutação também é a base dos modelos nulos em ecologia de comunidades. A pergunta aqui é se as espécies coocorrem menos do que o esperado ao acaso, o que pode indicar competição ou exigências de hábitat diferentes. Os dados são uma matriz de presença e ausência, com uma linha por espécie e uma coluna por local.
O ponds.csv não tem listas de espécies; por isso esta seção usa uma pequena matriz simulada para o propósito: seis espécies em doze lagoas.
M = [1 1 1 0 1 1 0 0 0 0 0 1;
1 0 1 1 1 1 0 1 0 0 0 0;
0 1 1 1 1 1 0 0 0 1 0 0;
1 0 0 0 0 0 1 1 1 0 1 1;
0 0 1 0 0 0 1 1 0 1 1 1;
0 1 0 0 0 0 1 0 1 1 1 1]
(ocorrencias_especie = vec(sum(M, dims = 2)),
riqueza_lagoa = vec(sum(M, dims = 1)))(ocorrencias_especie = [6, 6, 6, 6, 6, 6], riqueza_lagoa = [3, 3, 4, 2, 3, 3, 3, 3, 2, 3, 3, 4])Uma estatística comum é o C-score (Stone e Roberts 1990). Para cada par de espécies, conte as “unidades em tabuleiro de xadrez”: com \(r_i\) e \(r_j\) o número de locais que cada uma ocupa e \(S\) o número de locais que compartilham, o par contribui com \((r_i - S)(r_j - S)\). O C-score é a média sobre todos os pares. Quanto menos as espécies compartilham locais, maior ele é.
Escreva cscore(M) para uma matriz com espécies nas linhas e locais nas colunas, e verifique com checar("11.6", cscore). O restante desta seção usa essa função.
r = vec(sum(M, dims = 2)) dá as ocorrências de cada espécie. Percorra os pares com for i in 1:s-1, j in i+1:s, onde s = size(M, 1); os locais compartilhados são sum(M[i, :] .& M[j, :]).
function cscore(M)
r = vec(sum(M, dims = 2))
s = size(M, 1)
total, pares = 0.0, 0
for i in 1:s-1, j in i+1:s
compartilhados = sum(M[i, :] .& M[j, :])
total += (r[i] - compartilhados) * (r[j] - compartilhados)
pares += 1
end
return total / pares
end
cscore(M)13.866666666666667checar("11.6", cscore)✓ Correto!
Um C-score de cerca de 13,9 é alto? Depende do que assembleias aleatórias destas espécies nestas lagoas dariam. Espécies comuns compartilham mais locais que as raras só por acaso, e lagoas ricas têm mais espécies que as pobres; então um modelo nulo justo mantém fixos tanto o número de ocorrências de cada espécie quanto a riqueza de cada lagoa, e aleatoriza apenas quais espécies vão para onde (Gotelli 2000).
O algoritmo de troca (swap) faz exatamente isso. Ele sorteia duas espécies e duas lagoas; se formarem um tabuleiro — [1 0; 0 1] ou [0 1; 1 0] —, ele o inverte. Os totais de linhas e colunas nunca mudam:
function swap!(M; attempts = 1_000)
nsp, nlocais = size(M)
for _ in 1:attempts
i, j = sample(1:nsp, 2; replace = false)
k, l = sample(1:nlocais, 2; replace = false)
if M[i, k] == M[j, l] && M[i, l] == M[j, k] && M[i, k] != M[i, l]
M[i, k], M[i, l] = M[i, l], M[i, k]
M[j, k], M[j, l] = M[j, l], M[j, k]
end
end
return M
end
c_nulos = [cscore(swap!(copy(M))) for _ in 1:999]
c_obs = cscore(M)
(observado = c_obs,
media_nula = mean(c_nulos),
p = (count(c_nulos .>= c_obs) + 1) / (999 + 1),
ses = (c_obs - mean(c_nulos)) / std(c_nulos))(observado = 13.866666666666667, media_nula = 12.801067734401064, p = 0.011, ses = 3.3547575925379842)fig = Figure(size = (700, 380))
ax = Axis(fig[1, 1], xlabel = "C-score", ylabel = "Assembleias nulas")
hist!(ax, c_nulos, bins = 25, color = (:grey60, 0.8))
vlines!(ax, c_obs, color = :darkorange, linestyle = :dash, linewidth = 2)
fig
O C-score observado fica na cauda superior: cerca de 1% das assembleias nulas são tão segregadas quanto ela. O teste é unilateral, porque a pergunta era sobre menos coocorrência que o esperado. O tamanho de efeito padronizado (ses), a distância até a média nula em desvios-padrão nulos, é a forma usual de comparar resultados como esse entre matrizes de tamanhos diferentes.
Esta é uma versão didática. Para uma análise real, o pacote vegan do R (nullmodel, oecosimu) implementa muitos algoritmos nulos com o cuidado que eles exigem — inclusive quantas trocas bastam para que as matrizes nulas sejam independentes, algo que esta versão trata por alto. A Seção 9.7.1 mostra como chamá-lo a partir de Julia.
median em vez de mean como estatística. Qual você reportaria para o tamanho corporal, e por quê?