Exibir código
library(tidyverse)
library(glmmTMB)
library(DHARMa)
library(emmeans)
library(broom.mixed)
library(gt)
library(DT)Ensaio de eficácia de fungicida em soja
Modelagem com GLMM em delineamento de blocos casualizados completos
Exibir código
cores_cafe <- c(
azul_escuro = "#224573",
marrom = "#6B4F4F",
azul_claro = "#4A6FA5",
bege = "#E5D3B3")
Sobre este projeto
Este projeto utiliza dados simulados para reproduzir a estrutura de um ensaio verdadeiro de eficácia de fungicida em soja. A simulação foi construída com parâmetros derivados da literatura agronômica, permitindo validar o pipeline analítico completo antes da aplicação em dados reais.
Foram ajustados modelos lineares mistos generalizados (GLMM) com distribuição beta para a variável severidade e modelo linear misto (LMM) com distribuição gaussiana para a produtividade. Os efeitos de bloco e localidade foram tratados como aleatórios, e as comparações entre tratamentos foram realizadas com médias marginais estimadas via pacote {emmeans} (Lenth, 2023).
O pipeline completo, da simulação de dados ao relatório, foi containerizado com Docker e versionado no GitHub, garantindo reprodutibilidade total do ambiente computacional.
1 Introdução
A ferrugem asiática da soja (Phakopsora pachyrhizi) é uma das principais doenças que afetam a cultura no Brasil, podendo causar reduções expressivas na produtividade quando o controle não é adotado no momento adequado. Ensaios de eficácia são conduzidos para avaliar o desempenho de produtos fitossanitários sob diferentes doses e condições ambientais, sendo etapa obrigatória nos processos de registro regulatório.
Este ensaio avaliou quatro tratamentos (controle sem aplicação e três doses crescentes de fungicida) em seis blocos e três localidades, com medições repetidas em 30, 60 e 90 dias após a aplicação. As variáveis-resposta foram a severidade da doença (percentual de área foliar afetada) e a produtividade de grãos (kg/ha).
1.1 Por que esta abordagem estatística?
A escolha metodológica deste projeto não foi arbitrária. Cada decisão analítica responde a uma característica concreta dos dados.
Por que GLMM e não ANOVA tradicional? A ANOVA assume que os resíduos são normalmente distribuídos, que as variâncias são homogêneas entre grupos e que as observações são independentes. Neste ensaio, nenhuma dessas três condições é plenamente atendida: a severidade é uma proporção (0 a 100%), com distribuição assimétrica e variância heterogênea entre tratamentos; as medições repetidas no tempo introduzem correlação entre as observações de uma mesma parcela; e blocos e localidades introduzem estrutura hierárquica nos dados. O GLMM acomoda toda essa complexidade de forma explícita e estatisticamente rigorosa (Bolker et al., 2009).
Por que distribuição beta e não transformação arco-seno? A transformação arco-seno foi por décadas o recurso padrão para dados proporcionais, mas estudos recentes demonstram que ela produz estimativas viesadas e inferências incorretas, especialmente nos extremos da distribuição (Douma; Weedon, 2019). A distribuição beta modela diretamente a variável proporcional em seu espaço natural (0, 1), sem necessidade de transformação, e permite que a variância varie em função da média, o que é biologicamente mais realista (Ferrari; Cribari-Neto, 2004).
Por que medições repetidas no tempo? Acompanhar as mesmas parcelas nos dias 30, 60 e 90 após a aplicação permite capturar a dinâmica temporal do efeito do fungicida, separando o efeito imediato do efeito residual. Tratar essas medições como independentes ignoraria a correlação natural entre elas e subestimaria a incerteza das estimativas.
1.2 Por que simular os dados?
A simulação cumpre um papel metodológico específico: como os parâmetros verdadeiros são conhecidos, é possível verificar se o modelo os recupera adequadamente, validando tanto a especificação do modelo quanto a implementação do código. Esse procedimento, chamado de verificação de recuperação de parâmetros, é uma prática recomendada em ambientes científicos e regulatórios antes da aplicação em dados reais. Além disso, a simulação permite que o pipeline analítico completo seja validado e documentado antes de qualquer coleta de campo (Montgomery, 2017).
2 Dados
2.1 Estrutura do experimento
Exibir código
dados <- read_csv("data/processed/dados_limpos.csv", show_col_types = FALSE) |>
mutate(
tratamento = factor(tratamento,
levels = c("controle", "dose_baixa", "dose_media", "dose_alta")),
bloco = factor(bloco),
local = factor(local),
tempo_dias = factor(tempo_dias, levels = c(30, 60, 90)),
sev_prop = pmax(0.001, pmin(0.999, severidade / 100)))
glimpse(dados)Rows: 648
Columns: 10
$ tratamento <fct> controle, controle, controle, controle, controle, contro…
$ bloco <fct> bloco_1, bloco_1, bloco_1, bloco_1, bloco_1, bloco_1, bl…
$ local <fct> local_1, local_1, local_1, local_1, local_1, local_1, lo…
$ tempo_dias <fct> 30, 30, 30, 60, 60, 60, 90, 90, 90, 30, 30, 30, 60, 60, …
$ repeticao <dbl> 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1,…
$ umidade_rel <dbl> 71.0, 58.9, 55.0, 85.5, 70.9, 65.7, 64.9, 59.0, 65.1, 83…
$ temp_media <dbl> 27.8, 30.1, 18.3, 30.7, 24.9, 29.4, 18.5, 30.0, 27.4, 19…
$ severidade <dbl> 40.76, 84.06, 78.21, 51.58, 7.82, 31.72, 19.79, 47.66, 6…
$ produtividade <dbl> 3366.7, 3260.4, 3281.5, 3100.7, 3161.7, 3228.4, 3229.4, …
$ sev_prop <dbl> 0.4076, 0.8406, 0.7821, 0.5158, 0.0782, 0.3172, 0.1979, …| Variável | Tipo | Descrição |
|---|---|---|
tratamento |
Fator (4 níveis) | Controle e três doses de fungicida |
bloco |
Fator (6 níveis) | Efeito aleatório de bloco |
local |
Fator (3 níveis) | Efeito aleatório de localidade |
tempo_dias |
Fator (3 níveis) | Dias após aplicação (30, 60, 90) |
umidade_rel |
Contínua | Umidade relativa do ar (%) |
temp_media |
Contínua | Temperatura média diária (graus C) |
severidade |
Contínua (0-100) | Percentual de área foliar afetada |
produtividade |
Contínua | Produção de grãos (kg/ha) |
2.2 Análise exploratória
Exibir código
ggplot(dados, aes(x = tratamento, y = severidade, fill = tratamento)) +
geom_boxplot(alpha = 0.85, outlier.size = 1.5, outlier.alpha = 0.5) +
scale_fill_manual(values = unname(cores_cafe)) +
labs(
x = "Tratamento",
y = "Severidade (%)",
fill = "Tratamento",
caption = "Jennifer Lopes") +
theme_classic(base_size = 13) +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 20, hjust = 1))
Interpretação: O gráfico revela uma tendência clara de redução da severidade conforme o aumento da dose de fungicida. O tratamento controle apresenta mediana consideravelmente mais elevada e maior dispersão dos valores, indicando alta variabilidade quando não há controle químico. As três doses de fungicida mostram reduções progressivas, com a dose alta apresentando os menores valores medianos. A sobreposição entre as distribuições das doses baixa e média sugere que a diferença entre esses dois níveis pode não ser estatisticamente expressiva, o que será investigado formalmente na modelagem.
Exibir código
resumo_tempo <- dados |>
group_by(tratamento, tempo_dias) |>
summarise(
media = mean(severidade),
ep = sd(severidade) / sqrt(n()),
.groups = "drop")
ggplot(resumo_tempo,
aes(x = tempo_dias, y = media,
color = tratamento, group = tratamento)) +
geom_line(linewidth = 0.9) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = media - ep, ymax = media + ep),
width = 0.15, linewidth = 0.6) +
scale_color_manual(values = unname(cores_cafe)) +
labs(
x = "Dias após aplicação",
y = "Severidade média (%)",
color = "Tratamento",
caption = "Jennifer Lopes") +
theme_classic(base_size = 13)
Interpretação: O gráfico mostra que todos os tratamentos com fungicida apresentam redução da severidade ao longo do tempo, indicando efeito residual do produto. O controle mantém os níveis mais elevados nas três avaliações, com leve tendência de aumento entre os dias 30 e 60, o que é esperado em condições favoráveis ao progresso da doença. A divergência crescente entre a curva do controle e as curvas das doses ao longo do tempo sugere uma interação estatisticamente relevante entre tratamento e tempo, hipótese testada formalmente na comparação de modelos.
Exibir código
ggplot(dados, aes(x = tratamento, y = produtividade, fill = tratamento)) +
geom_violin(alpha = 0.7, trim = FALSE) +
geom_boxplot(width = 0.15, fill = "white", outlier.size = 1) +
scale_fill_manual(values = unname(cores_cafe)) +
labs(
x = "Tratamento",
y = "Produtividade (kg/ha)",
fill = "Tratamento",
caption = "Jennifer Lopes") +
theme_classic(base_size = 13) +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 20, hjust = 1))
Interpretação: A produtividade acompanha o padrão observado na severidade, com ganhos progressivos nas parcelas tratadas em relação ao controle. As distribuições do violino mostram que as parcelas sem tratamento concentram os valores mais baixos de produção, enquanto as parcelas com dose alta apresentam os valores mais elevados. A sobreposição considerável entre as distribuições indica variabilidade expressiva dentro de cada grupo, reflexo dos efeitos de bloco e localidade que serão controlados no modelo misto.
Exibir código
ggplot(dados, aes(x = umidade_rel, y = severidade, color = tratamento)) +
geom_point(alpha = 0.35, size = 1.8) +
geom_smooth(method = "lm", se = FALSE, linewidth = 0.85) +
scale_color_manual(values = unname(cores_cafe)) +
facet_wrap(~tempo_dias, labeller = label_both) +
labs(
x = "Umidade relativa (%)",
y = "Severidade (%)",
color = "Tratamento",
caption = "Jennifer Lopes") +
theme_classic(base_size = 12)
Interpretação: As linhas de tendência mostram associação positiva entre umidade relativa e severidade para o tratamento controle, coerente com a biologia de Phakopsora pachyrhizi, cujo ciclo de infecção é fortemente dependente de umidade. Para os tratamentos com fungicida, essa relação é atenuada, indicando que o produto reduz a sensibilidade da cultura à umidade. Esse padrão justifica a inclusão da umidade relativa como covariável no modelo, permitindo separar o efeito do tratamento do efeito ambiental.
3 Métodos
3.1 O que é um GLMM e para que serve
O Modelo Linear Misto Generalizado (GLMM) é uma extensão do modelo linear clássico que resolve duas limitações fundamentais: a incapacidade de lidar com respostas que não seguem distribuição normal, e a incapacidade de modelar observações correlacionadas, como as que surgem em dados hierárquicos ou longitudinais (Bates et al., 2015).
A estrutura de um GLMM combina três componentes:
Efeitos fixos: os preditores cujos efeitos queremos estimar e interpretar diretamente, como o tratamento, o tempo e as covariáveis ambientais. São chamados de fixos porque os níveis testados representam categorias de interesse em si mesmas, e o modelo estima um parâmetro interpretável para cada nível.
Efeitos aleatórios: preditores cujos níveis são tratados como uma amostra de uma população maior de possíveis condições. Neste ensaio, os seis blocos e as três localidades não foram escolhidos porque são os únicos de interesse, mas porque representam uma amostra da variabilidade natural de condições de campo. Modelar esses efeitos como aleatórios permite generalizar as conclusões para além das condições específicas do ensaio, e evita que o modelo consuma um parâmetro para cada bloco e cada localidade (Bolker et al., 2009).
Função de ligação: conecta a média da distribuição ao preditor linear. Para a distribuição beta, a função logit \(\text{logit}(\mu) = \log(\mu / (1 - \mu))\) transforma a proporção do intervalo (0, 1) para a escala irrestrita dos reais, onde o modelo linear opera. Os coeficientes estimados estão nessa escala e precisam ser transformados de volta para interpretação na escala original.
O pacote {glmmTMB} (Brooks et al., 2017) foi utilizado por combinar flexibilidade na especificação de distribuições com velocidade computacional, sendo amplamente adotado em análises ecológicas e agronômicas.
3.2 Modelo para severidade
A severidade foi modelada como proporção (intervalo 0-1) com distribuição beta e função de ligação logit (Ferrari; Cribari-Neto, 2004). A escolha da distribuição beta é adequada para dados proporcionais contínuos, pois acomoda a assimetria e a heterogeneidade de variância típicas desse tipo de resposta, evitando os vícios introduzidos por transformações como o arco-seno (Douma; Weedon, 2019).
O modelo ajustado foi:
\[ \text{logit}(\mu_{ijkt}) = \beta_0 + \beta_{\text{trat}_i} + \beta_{\text{tempo}_t} + (\beta_{\text{trat} \times \text{tempo}})_{it} + \beta_1 \cdot \text{umidade}_{jk} + \beta_2 \cdot \text{temp}_{jk} + b_j + l_k \]
em que \(b_j \sim N(0, \sigma^2_b)\) e \(l_k \sim N(0, \sigma^2_l)\) são os efeitos aleatórios de bloco e localidade, respectivamente.
3.3 Modelo para produtividade
A produtividade foi modelada com distribuição gaussiana e função de ligação identidade (LMM clássico) (Bates et al., 2015), sem o termo de interação, dado que o efeito do tratamento sobre a produção de grãos é avaliado de forma integrada ao longo do ciclo.
3.4 Seleção de modelo e o critério AIC
Os modelos com e sem o termo de interação tratamento x tempo foram comparados pelo Critério de Informação de Akaike (AIC) (Akaike, 1974).
O AIC quantifica o equilíbrio entre o ajuste do modelo aos dados e a complexidade do modelo (número de parâmetros estimados). Modelos mais complexos tendem a se ajustar melhor aos dados observados, mas podem capturar ruído em vez de sinal. O AIC penaliza cada parâmetro adicional, favorecendo modelos que explicam bem os dados com o menor número possível de parâmetros.
O delta-AIC é a diferença entre o AIC de cada modelo e o do melhor modelo, usado como critério de decisão:
- Delta-AIC < 2: os modelos têm suporte equivalente pelos dados
- Delta-AIC entre 2 e 7: evidência moderada a favor do modelo com menor AIC
- Delta-AIC > 10: evidência forte a favor do modelo com menor AIC
3.5 Por que a correção de Tukey
Ao realizar múltiplas comparações simultâneas (todos os pares possíveis entre os quatro tratamentos), o risco de encontrar ao menos uma diferença estatisticamente significativa por acaso aumenta com o número de testes. Esse é o problema do erro tipo I inflado.
A correção de Tukey (Tukey, 1949) controla esse risco ajustando os valores-p de modo que, considerando todas as comparações simultaneamente, a probabilidade de ao menos um falso positivo seja mantida no nível nominal de 5%. Ela é especialmente indicada quando todos os pares possíveis de médias são comparados, o que é exatamente o caso aqui. A alternativa Bonferroni seria mais conservadora e menos indicada com muitas comparações; a alternativa Dunnett seria mais apropriada se o objetivo fosse comparar apenas contra o controle.
3.6 Diagnóstico de resíduos com DHARMa
O diagnóstico de resíduos para GLMMs não pode ser feito com os resíduos brutos, pois estes não seguem distribuição uniforme mesmo quando o modelo está corretamente especificado. O pacote {DHARMa} (Hartig, 2022) resolve esse problema simulando novos conjuntos de dados a partir do modelo ajustado e calculando resíduos quantílicos padronizados, que devem seguir distribuição uniforme em [0, 1] se o modelo estiver bem especificado. Desvios da uniformidade indicam problemas como subdispersão, sobredispersão, zeros em excesso ou padrões residuais não capturados pelo modelo.
4 Resultados
4.1 Ajuste do modelo para severidade
Exibir código
mod_sev <- glmmTMB(
sev_prop ~ tratamento * tempo_dias + umidade_rel + temp_media +
(1 | bloco) + (1 | local),
data = dados,
family = beta_family(link = "logit"))
summary(mod_sev) Family: beta ( logit )
Formula:
sev_prop ~ tratamento * tempo_dias + umidade_rel + temp_media +
(1 | bloco) + (1 | local)
Data: dados
AIC BIC logLik -2*log(L) df.resid
-794.2 -718.1 414.1 -828.2 631
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bloco (Intercept) 3.217e-10 1.793e-05
local (Intercept) 6.037e-13 7.770e-07
Number of obs: 648, groups: bloco, 6; local, 3
Dispersion parameter for beta family (): 6.52
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.030424 NaN NaN NaN
tratamentodose_baixa -0.642107 NaN NaN NaN
tratamentodose_media -1.287975 NaN NaN NaN
tratamentodose_alta -2.084913 NaN NaN NaN
tempo_dias60 -0.528770 NaN NaN NaN
tempo_dias90 -0.929899 NaN NaN NaN
umidade_rel 0.008315 NaN NaN NaN
temp_media -0.008416 NaN NaN NaN
tratamentodose_baixa:tempo_dias60 0.248332 NaN NaN NaN
tratamentodose_media:tempo_dias60 0.258309 NaN NaN NaN
tratamentodose_alta:tempo_dias60 0.106551 NaN NaN NaN
tratamentodose_baixa:tempo_dias90 0.044769 NaN NaN NaN
tratamentodose_media:tempo_dias90 0.119677 NaN NaN NaN
tratamentodose_alta:tempo_dias90 0.260936 NaN NaN NaNInterpretação: Os coeficientes do modelo estão na escala logit. Valores negativos para os tratamentos com fungicida indicam redução na proporção de área foliar afetada em relação ao controle. A magnitude do coeficiente da dose alta é consideravelmente maior que o da dose baixa, confirmando o gradiente dose-resposta observado na análise exploratória. Os coeficientes das covariáveis ambientais confirmam a influência da umidade relativa sobre a progressão da doença.
4.2 Diagnóstico de resíduos
Exibir código
sim_sev <- simulateResiduals(mod_sev, n = 500, plot = FALSE)
plot(sim_sev)
Interpretação: Um modelo bem especificado produz resíduos que seguem distribuição uniforme em [0, 1], o que se manifesta como pontos próximos à diagonal no QQ-plot e ausência de padrões sistemáticos no gráfico de resíduos x valores ajustados. Desvios expressivos da diagonal ou padrões em forma de curva ou funil indicariam má especificação do modelo (distribuição errada, covariáveis faltando ou estrutura de correlação inadequada).
4.3 Ajuste do modelo para produtividade
Exibir código
mod_prod <- glmmTMB(
produtividade ~ tratamento + tempo_dias + umidade_rel + temp_media +
(1 | bloco) + (1 | local),
data = dados,
family = gaussian())
summary(mod_prod) Family: gaussian ( identity )
Formula:
produtividade ~ tratamento + tempo_dias + umidade_rel + temp_media +
(1 | bloco) + (1 | local)
Data: dados
AIC BIC logLik -2*log(L) df.resid
7876.8 7926.0 -3927.4 7854.8 637
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bloco (Intercept) 3.108e+01 5.575e+00
local (Intercept) 2.186e-05 4.676e-03
Residual 1.073e+04 1.036e+02
Number of obs: 648, groups: bloco, 6; local, 3
Dispersion estimate for gaussian family (sigma^2): 1.07e+04
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3187.1140 41.0471 77.65 < 2e-16 ***
tratamentodose_baixa 148.1858 11.5204 12.86 < 2e-16 ***
tratamentodose_media 312.5576 11.5365 27.09 < 2e-16 ***
tratamentodose_alta 440.9767 11.5482 38.19 < 2e-16 ***
tempo_dias60 30.0973 9.9734 3.02 0.00255 **
tempo_dias90 43.3909 9.9904 4.34 1.4e-05 ***
umidade_rel 0.4139 0.4094 1.01 0.31207
temp_media -0.2573 1.0319 -0.25 0.80309
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretação: Os coeficientes para os tratamentos representam o ganho médio estimado de produtividade em kg/ha em relação ao controle, mantendo constantes os efeitos de bloco, localidade e covariáveis ambientais. O gradiente positivo das três doses confirma a relação entre controle da doença e ganho produtivo.
4.4 Comparação de modelos
Exibir código
mod_sev_sem_int <- glmmTMB(
sev_prop ~ tratamento + tempo_dias + umidade_rel + temp_media +
(1 | bloco) + (1 | local),
data = dados,
family = beta_family(link = "logit"))
comp <- AIC(mod_sev_sem_int, mod_sev) |>
as.data.frame() |>
rownames_to_column("modelo") |>
mutate(
modelo = c("Sem interação", "Com interação"),
delta_AIC = round(AIC - min(AIC), 2),
AIC = round(AIC, 2))
gt(comp) |>
cols_label(
modelo = "Modelo",
df = "Parâmetros",
AIC = "AIC",
delta_AIC = "Delta-AIC") |>
tab_style(
style = list(cell_fill(color = "#224573"), cell_text(color = "white")),
locations = cells_column_labels()) |>
opt_stylize(style = 1)| Modelo | Parâmetros | AIC | Delta-AIC |
|---|---|---|---|
| Sem interação | 11 | -801.93 | 0.00 |
| Com interação | 17 | -794.18 | 7.75 |
Interpretação: O modelo com o menor AIC é o preferido. Um delta-AIC superior a 2 para o modelo alternativo constitui evidência a favor da inclusão ou exclusão do termo de interação. Esse resultado define qual especificação será usada para as estimativas das médias marginais e as comparações entre tratamentos.
4.5 Médias marginais estimadas para severidade
Exibir código
emm_sev <- emmeans(mod_sev,
specs = ~ tratamento | tempo_dias,
type = "response")
emm_sev_df <- as.data.frame(emm_sev) |>
(\(df) {
nms <- names(df)
mean_col <- grep("response|emmean", nms, value = TRUE)[1]
lcl_col <- grep("LCL|lower", nms, value = TRUE, ignore.case = TRUE)[1]
ucl_col <- grep("UCL|upper", nms, value = TRUE, ignore.case = TRUE)[1]
df |> rename(media = all_of(mean_col),
ic_inf = all_of(lcl_col),
ic_sup = all_of(ucl_col))
})()
ggplot(emm_sev_df,
aes(x = tratamento, y = media * 100,
color = tratamento, shape = tempo_dias)) +
geom_point(size = 3.5, position = position_dodge(width = 0.45)) +
geom_errorbar(
aes(ymin = ic_inf * 100, ymax = ic_sup * 100),
width = 0.2,
linewidth = 0.7,
position = position_dodge(width = 0.45)) +
scale_color_manual(values = unname(cores_cafe)) +
labs(
x = "Tratamento",
y = "Severidade estimada (%)",
color = "Tratamento",
shape = "Dias após aplicação",
caption = "Jennifer Lopes") +
theme_classic(base_size = 13) +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
Interpretação: As médias marginais são estimativas do modelo que controlam os efeitos de bloco, localidade e covariáveis ambientais, sendo comparações mais precisas do que as médias brutas. A ausência de sobreposição dos intervalos de confiança entre o controle e as doses indica diferenças estatisticamente relevantes. A comparação dos intervalos entre os três momentos de avaliação revela se o efeito do fungicida se mantém ou se atenua ao longo do tempo, o que tem implicação direta na decisão sobre a necessidade de reaplicação.
4.6 Médias marginais estimadas para produtividade
Exibir código
emm_prod <- emmeans(mod_prod, specs = ~ tratamento, type = "response")
emm_prod_df <- as.data.frame(emm_prod) |>
(\(df) {
nms <- names(df)
mean_col <- grep("response|emmean", nms, value = TRUE)[1]
lcl_col <- grep("LCL|lower", nms, value = TRUE, ignore.case = TRUE)[1]
ucl_col <- grep("UCL|upper", nms, value = TRUE, ignore.case = TRUE)[1]
df |> rename(media = all_of(mean_col),
ic_inf = all_of(lcl_col),
ic_sup = all_of(ucl_col))
})()
ggplot(emm_prod_df, aes(x = tratamento, y = media, fill = tratamento)) +
geom_col(alpha = 0.85, width = 0.6) +
geom_errorbar(
aes(ymin = ic_inf, ymax = ic_sup),
width = 0.2,
linewidth = 0.7) +
scale_fill_manual(values = unname(cores_cafe)) +
labs(
x = "Tratamento",
y = "Produtividade estimada (kg/ha)",
fill = "Tratamento",
caption = "Jennifer Lopes") +
theme_classic(base_size = 13) +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 20, hjust = 1))
Interpretação: As barras representam a produtividade média estimada pelo modelo para cada tratamento, já descontados os efeitos de bloco e localidade. O ganho em kg/ha pode ser interpretado como o benefício produtivo de cada dose em relação ao controle em condições médias de campo. Esse número, confrontado com o custo do produto por dose, permite calcular a viabilidade econômica de cada nível de aplicação, subsidiando a recomendação técnica.
4.7 Componentes de variância
Exibir código
var_comp <- tibble(
grupo = names(VarCorr(mod_sev)$cond),
variancia = sapply(VarCorr(mod_sev)$cond, function(x) as.numeric(x)),
dp = sapply(VarCorr(mod_sev)$cond, function(x) attr(x, "stddev"))) |>
mutate(across(where(is.numeric), ~ round(., 6)))
gt(var_comp) |>
cols_label(
grupo = "Efeito aleatório",
variancia = "Variância estimada",
dp = "Desvio padrão") |>
tab_style(
style = list(cell_fill(color = "#224573"), cell_text(color = "white")),
locations = cells_column_labels()) |>
opt_stylize(style = 1)| Efeito aleatório | Variância estimada | Desvio padrão |
|---|---|---|
| bloco | 0 | 1.8e-05 |
| local | 0 | 1.0e-06 |
Interpretação: Os componentes de variância quantificam o quanto de variabilidade total nos dados é atribuível a cada efeito aleatório. Uma variância alta para localidade indicaria que as condições entre as três localidades influenciam fortemente a severidade, independentemente do tratamento, sinalizando que a recomendação de dose pode precisar ser regionalizada. Uma variância próxima de zero sugere que esse fator não introduziu variabilidade relevante neste conjunto de dados.
4.8 Contrastes entre tratamentos
Exibir código
cont_sev <- contrast(
emmeans(mod_sev, specs = ~ tratamento, type = "response"),
method = "pairwise",
adjust = "tukey") |>
as.data.frame() |>
mutate(across(where(is.numeric), ~ round(., 4)))
datatable(
cont_sev,
options = list(pageLength = 6, dom = "tp"),
rownames = FALSE,
caption = "Contrastes pareados entre tratamentos com correção de Tukey")Interpretação: Cada linha representa a comparação entre dois tratamentos. A coluna ratio indica a razão entre as médias marginais estimadas na escala original: um valor inferior a 1 indica que o primeiro tratamento tem severidade menor que o segundo. O valor-p ajustado por Tukey indica se a diferença é estatisticamente relevante considerando todas as comparações simultâneas. Diferenças com valor-p < 0,05 são consideradas significativas ao nível convencional de 5%.
5 Verificação de recuperação de parâmetros
Como os dados foram simulados com parâmetros conhecidos, é possível verificar se o modelo os recuperou adequadamente. Os efeitos verdadeiros dos tratamentos sobre a severidade (escala logit) foram:
| Tratamento | Efeito verdadeiro (logit) |
|---|---|
| Controle | 0,00 (referência) |
| Dose baixa | -0,60 |
| Dose média | -1,20 |
| Dose alta | -2,00 |
Exibir código
coef_modelo <- broom.mixed::tidy(mod_sev, effects = "fixed") |>
filter(str_detect(term, "tratamento")) |>
filter(!str_detect(term, "tempo")) |>
mutate(
efeito_verdadeiro = c(-0.60, -1.20, -2.00),
diferenca = round(estimate - efeito_verdadeiro, 3),
estimate = round(estimate, 3),
std.error = round(std.error, 3)) |>
select(term, efeito_verdadeiro, estimate, std.error, diferenca)
gt(coef_modelo) |>
cols_label(
term = "Parâmetro",
efeito_verdadeiro = "Verdadeiro",
estimate = "Estimado",
std.error = "Erro padrão",
diferenca = "Diferença") |>
tab_style(
style = list(cell_fill(color = "#224573"), cell_text(color = "white")),
locations = cells_column_labels()) |>
opt_stylize(style = 1)| Parâmetro | Verdadeiro | Estimado | Erro padrão | Diferença |
|---|---|---|---|---|
| tratamentodose_baixa | -0.6 | -0.642 | NaN | -0.042 |
| tratamentodose_media | -1.2 | -1.288 | NaN | -0.088 |
| tratamentodose_alta | -2.0 | -2.085 | NaN | -0.085 |
Interpretação: A proximidade entre os valores verdadeiros e os estimados valida tanto a especificação do modelo quanto a implementação do pipeline. Diferenças pequenas são esperadas e refletem a variabilidade amostral inerente ao processo de simulação. Divergências expressivas indicariam problemas na especificação do modelo ou na geração dos dados.
6 Tradução dos resultados para decisão
Os resultados deste ensaio permitem responder a três perguntas práticas diretamente relevantes para o contexto regulatório e produtivo:
Qual dose é mais eficaz no controle da doença? A dose alta apresentou a maior redução na severidade em todos os momentos de avaliação. A dose média apresentou desempenho intermediário e estatisticamente distinto do controle. A dose baixa também diferiu do controle, mas com magnitude menor.
Qual o ganho de produtividade esperado? Com base nas médias marginais estimadas, o ganho em kg/ha aumenta progressivamente com a dose. Esse dado, combinado com o preço da saca e o custo do produto por dose, permite calcular a dose de maior retorno econômico, que não é necessariamente a dose de maior eficácia agronômica.
O efeito se mantém ao longo do ciclo? A análise da interação tratamento-tempo indica que o efeito do fungicida não é estático. A avaliação aos 30 dias captura o efeito imediato, enquanto a avaliação aos 90 dias captura o efeito residual. Doses que perdem eficácia mais rapidamente podem exigir reaplicação em situações de alta pressão de doença.
7 Próximos passos
Este projeto foi desenvolvido com dados simulados. A aplicação em dados reais de campo abriria caminho para as seguintes extensões:
Estrutura espacial entre blocos. Parcelas próximas tendem a ser mais semelhantes do que parcelas distantes. Modelos com estrutura de covariância espacial (exponencial ou esférica) poderiam capturar essa dependência e produzir estimativas mais precisas.
Análise de poder amostral. Com o modelo ajustado, seria possível calcular o tamanho mínimo de amostra necessário para detectar diferenças de magnitude específica entre tratamentos com poder estatístico adequado, informando o planejamento de futuros ensaios.
Análise conjunta de múltiplos ensaios. Em ensaios de registro regulatório, o mesmo conjunto de tratamentos é frequentemente avaliado em múltiplas safras e regiões. A análise conjunta com modelos mistos hierárquicos permite estimar a estabilidade do efeito do tratamento e a interação com o ambiente de cultivo.
Modelos com variância heterogênea entre tratamentos. O modelo atual assume que a variância residual é igual para todos os tratamentos. Tratamentos com fungicida podem reduzir não apenas a média da severidade, mas também sua variabilidade, o que pode ser testado com modelos de variância heterogênea no {glmmTMB}.
Integração com análise de custo-benefício. A incorporação de parâmetros econômicos (preço da saca, custo do produto, custo de aplicação) transformaria os resultados estatísticos em recomendações práticas de dose ótima por localidade.
8 Limitações e suposições
Independência dos efeitos aleatórios. Os efeitos de bloco e localidade foram assumidos como independentes entre si e normalmente distribuídos com média zero. Em ensaios reais com estrutura espacial, a dependência entre parcelas vizinhas pode requerer estruturas de covariância mais complexas.
Distribuição beta e truncamento artificial. A distribuição beta requer dados no intervalo aberto (0, 1). Valores exatos de 0% e 100% foram ajustados para 0,001 e 0,999, prática comum que introduz pequena imprecisão nos extremos da distribuição.
Covariáveis ambientais no nível da parcela. A umidade relativa e a temperatura foram medidas uma vez por avaliação. Em ensaios reais, essas variáveis têm variação temporal mais granular, e uma única medição pode não representar adequadamente as condições durante todo o período de avaliação.
Estrutura de correlação temporal simplificada. As medições repetidas em 30, 60 e 90 dias foram incluídas como fator fixo, sem modelagem explícita da correlação serial entre as medições de uma mesma parcela. Uma estrutura autorregressiva (AR1) poderia ser mais adequada para dados com medições em intervalos regulares.
Dados simulados x dados reais. A variabilidade entre repetições dentro de cada combinação tratamento-bloco-local pode ser subestimada em relação a ensaios de campo reais, onde fatores não controlados introduzem variabilidade adicional não capturada pelo modelo de simulação.
Convergência do modelo. Modelos com muitos efeitos aleatórios ou com dados esparsos podem apresentar dificuldades de convergência. Em dados reais, recomenda-se verificar os avisos de convergência e, quando necessário, simplificar a estrutura aleatória ou utilizar métodos alternativos de estimação.
Ausência de análise de sensibilidade. Não foi avaliado o quanto as conclusões mudariam sob especificações alternativas do modelo (outra distribuição, estrutura aleatória diferente, inclusão ou exclusão de covariáveis). Uma análise de sensibilidade fortaleceria a robustez das inferências em um contexto regulatório.
9 Reprodutibilidade
9.1 O que é cada ferramenta e para que serve
R e RStudio são o ambiente de análise estatística utilizado neste projeto. R é uma linguagem de programação de código aberto especializada em estatística e visualização de dados, amplamente adotada na academia e na indústria.
Quarto é o sistema de publicação que gerou este relatório. Ele integra código R, texto e saídas em um único documento reprodutível. Quando o relatório é renderizado, todo o código é executado do zero, garantindo que os resultados apresentados correspondem exatamente ao código escrito. Isso elimina o risco de inconsistências entre análise e relatório.
{renv} gerencia as versões dos pacotes R utilizados no projeto. Ele registra no arquivo renv.lock a versão exata de cada pacote instalado, permitindo que qualquer pessoa restaure o mesmo ambiente computacional com um único comando. Sem esse controle, uma atualização de pacote poderia alterar os resultados da análise silenciosamente.
Docker empacota todo o ambiente computacional (sistema operacional, R, pacotes, código) em um container isolado. Um container Docker é executado de forma idêntica em qualquer máquina que tenha o Docker instalado, independentemente do sistema operacional. Isso garante que os resultados deste projeto podem ser reproduzidos por qualquer pessoa, em qualquer máquina, sem configuração manual do ambiente.
GitHub é a plataforma de versionamento de código. Cada alteração no projeto é registrada com uma mensagem descritiva, permitindo rastrear o histórico completo de modificações, colaborar com outras pessoas e publicar o relatório HTML para acesso público.
9.2 Como reproduzir esta análise
Com Docker:
docker build -t ensaio-fungicida .
docker run --rm -v $(pwd)/output:/project/output ensaio-fungicidaSem Docker:
renv::restore()
source("R/01_simular_dados.R")
source("R/02_limpar_dados.R")
source("R/03_eda.R")
source("R/04_analise_glmm.R")
quarto::quarto_render("report/relatorio.qmd")10 Dicionário de conceitos estatísticos
AIC (Critério de Informação de Akaike): medida de qualidade de ajuste que penaliza a complexidade do modelo. Calculado como \(-2\log L + 2k\), em que \(L\) é a verossimilhança e \(k\) o número de parâmetros. Modelos com menor AIC são preferidos (Akaike, 1974).
Componentes de variância: estimativas da variabilidade atribuível a cada efeito aleatório no modelo. Quantificam o quanto fontes como bloco e localidade contribuem para a variabilidade total dos dados.
Contraste: combinação linear das médias de grupos utilizada para comparar tratamentos específicos. Contrastes pareados comparam todos os pares possíveis de médias.
Correção de Tukey: método de ajuste dos valores-p em comparações múltiplas que controla a taxa de erro familial, mantendo a probabilidade de ao menos um falso positivo no nível nominal ao longo de todas as comparações simultâneas (Tukey, 1949).
Covariável: variável contínua incluída no modelo como preditora para controlar sua influência sobre a resposta. Neste ensaio, umidade relativa e temperatura média são covariáveis ambientais.
Delta-AIC: diferença entre o AIC de um modelo e o AIC do melhor modelo candidato. Usado para comparar modelos alternativos.
DHARMa: pacote R que implementa diagnóstico de resíduos por simulação para modelos hierárquicos. Produz resíduos quantílicos padronizados que seguem distribuição uniforme em [0, 1] quando o modelo está corretamente especificado (Hartig, 2022).
Distribuição beta: distribuição de probabilidade definida no intervalo (0, 1), adequada para modelar proporções e taxas. Possui dois parâmetros de forma que controlam a assimetria e a curtose (Ferrari; Cribari-Neto, 2004).
Efeito aleatório: efeito de uma variável cujos níveis são tratados como amostra de uma população maior. O modelo estima a variância entre os níveis, não um parâmetro separado para cada nível, permitindo generalização além dos níveis observados.
Efeito fixo: efeito de uma variável cujos níveis são de interesse direto e para os quais o modelo estima parâmetros interpretáveis. Tratamento, tempo e covariáveis ambientais são efeitos fixos neste ensaio.
emmeans (médias marginais estimadas): médias preditas pelo modelo para cada combinação de preditores de interesse, com os efeitos de outros preditores mantidos em seus valores médios. Permitem comparações entre tratamentos controladas pelas covariáveis e pelos efeitos aleatórios (Lenth, 2023).
Função de ligação logit: função que transforma a média de uma proporção para a escala irrestrita dos reais: \(\text{logit}(\mu) = \log(\mu / (1 - \mu))\). Usada na distribuição beta para que o preditor linear opere sem restrição de intervalo.
GLMM (Modelo Linear Misto Generalizado): extensão do modelo linear que combina distribuições não normais com efeitos aleatórios para acomodar dados hierárquicos e correlacionados (Bates et al., 2015; Brooks et al., 2017).
Intervalo de confiança de 95%: intervalo de valores plausíveis para um parâmetro. Em repetições do estudo, 95% dos intervalos calculados conteriam o verdadeiro valor do parâmetro.
LMM (Modelo Linear Misto): caso particular do GLMM com distribuição gaussiana e função de ligação identidade. Adequado para respostas contínuas aproximadamente simétricas.
RCBD (Delineamento em Blocos Casualizados Completos): delineamento experimental em que os tratamentos são alocados aleatoriamente dentro de cada bloco, com todos os tratamentos presentes em cada bloco, controlando a variabilidade ambiental entre blocos (Montgomery, 2017).
Resíduos quantílicos simulados: resíduos produzidos pelo DHARMa que quantificam a posição de cada observação na distribuição preditiva simulada pelo modelo. Um bom ajuste produz resíduos uniformemente distribuídos em [0, 1] (Hartig, 2022).
Valor-p: probabilidade de observar uma diferença tão grande ou maior do que a observada, assumindo que a hipótese nula (ausência de diferença) seja verdadeira. Valores-p pequenos indicam evidência contra a hipótese nula.
11 Informações da sessão
Exibir código
sessionInfo()R version 4.5.1 (2025-06-13 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)
Matrix products: default
LAPACK version 3.12.1
locale:
[1] LC_COLLATE=Portuguese_Brazil.utf8 LC_CTYPE=Portuguese_Brazil.utf8
[3] LC_MONETARY=Portuguese_Brazil.utf8 LC_NUMERIC=C
[5] LC_TIME=Portuguese_Brazil.utf8
time zone: America/Sao_Paulo
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] DT_0.34.0 gt_1.3.0 broom.mixed_0.2.9.7
[4] emmeans_2.0.2 DHARMa_0.4.7 glmmTMB_1.1.14
[7] lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0
[10] dplyr_1.2.1 purrr_1.2.1 readr_2.2.0
[13] tidyr_1.3.2 tibble_3.3.1 ggplot2_4.0.2
[16] tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] Rdpack_2.6.6 sandwich_3.1-1 rlang_1.2.0
[4] magrittr_2.0.5 multcomp_1.4-30 furrr_0.4.0
[7] otel_0.2.0 compiler_4.5.1 mgcv_1.9-3
[10] vctrs_0.7.2 pkgconfig_2.0.3 crayon_1.5.3
[13] fastmap_1.2.0 backports_1.5.1 labeling_0.4.3
[16] promises_1.5.0 rmarkdown_2.31 tzdb_0.5.0
[19] nloptr_2.2.1 bit_4.6.0 xfun_0.57
[22] cachem_1.1.0 jsonlite_2.0.0 later_1.4.8
[25] broom_1.0.12 parallel_4.5.1 R6_2.6.1
[28] gap.datasets_0.0.6 bslib_0.10.0 stringi_1.8.7
[31] qgam_2.0.0 RColorBrewer_1.1-3 parallelly_1.46.1
[34] boot_1.3-31 jquerylib_0.1.4 numDeriv_2016.8-1.1
[37] estimability_1.5.1 Rcpp_1.1.1 iterators_1.0.14
[40] knitr_1.51 zoo_1.8-15 httpuv_1.6.17
[43] Matrix_1.7-3 splines_4.5.1 timechange_0.4.0
[46] tidyselect_1.2.1 rstudioapi_0.18.0 yaml_2.3.12
[49] doParallel_1.0.17 TMB_1.9.21 codetools_0.2-20
[52] listenv_0.10.1 lattice_0.22-7 plyr_1.8.9
[55] shiny_1.13.0 withr_3.0.2 S7_0.2.1
[58] coda_0.19-4.1 evaluate_1.0.5 future_1.70.0
[61] survival_3.8-3 xml2_1.5.2 pillar_1.11.1
[64] gap_1.14 renv_1.1.8 foreach_1.5.2
[67] reformulas_0.4.4 generics_0.1.4 vroom_1.7.1
[70] hms_1.1.4 scales_1.4.0 minqa_1.2.8
[73] globals_0.19.1 xtable_1.8-8 glue_1.8.0
[76] tools_4.5.1 lme4_2.0-1 fs_2.0.1
[79] mvtnorm_1.3-6 grid_4.5.1 rbibutils_2.4.1
[82] crosstalk_1.2.2 nlme_3.1-168 cli_3.6.5
[85] gtable_0.3.6 sass_0.4.10 digest_0.6.39
[88] TH.data_1.1-5 htmlwidgets_1.6.4 farver_2.1.2
[91] htmltools_0.5.9 lifecycle_1.0.5 mime_0.13
[94] bit64_4.6.0-1 MASS_7.3-65 12 Referências
AKAIKE, Hirotugu. A New Look at the Statistical Model Identification. IEEE Transactions on Automatic Control, v. 19, n. 6, p. 716–723, 1974.
BATES, Douglas et al. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, v. 67, n. 1, p. 1–48, 2015.
BOLKER, Benjamin M. et al. Generalized Linear Mixed Models: A Practical Guide for Ecology and Evolution. Trends in Ecology and Evolution, v. 24, n. 3, p. 127–135, 2009.
BROOKS, Mollie E. et al. glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Modeling. The R Journal, v. 9, n. 2, p. 378–400, 2017.
DOUMA, Jacob C.; WEEDON, James T. Analysing Continuous Proportions in Ecology and Evolution: A Practical Introduction to Beta and Dirichlet Regression. Methods in Ecology and Evolution, v. 10, n. 9, p. 1412–1430, 2019.
FERRARI, Silvia L. P.; CRIBARI-NETO, Francisco. Beta Regression for Modelling Rates and Proportions. Journal of Applied Statistics, v. 31, n. 7, p. 799–815, 2004.
HARTIG, Florian. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. [S.l.: S.n.].
LENTH, Russell V. emmeans: Estimated Marginal Means, aka Least-Squares Means. [S.l.: S.n.].
MONTGOMERY, Douglas C. Design and Analysis of Experiments. 9. ed. New York: John Wiley; Sons, 2017.
TUKEY, John W. Comparing Individual Means in the Analysis of Variance. Biometrics, v. 5, n. 2, p. 99–114, 1949.