# Instalar os pacotes apenas se ainda não estiverem instalados:
# install.packages(c("ExpDes", "lmtest", "MASS", "readxl"))
# Carregar os pacotes utilizados no material:
library(ExpDes)
library(lmtest)
library(MASS)
library(readxl)
# Nível de significância adotado nos exemplos:
alfa <- 0.05Planejamento e Análise de Experimentos I
Bacharelado em Estatística - UFMT
Análise de experimentos no R
Este material apresenta a análise computacional dos delineamentos inteiramente casualizados, em blocos casualizados e em quadrado latino. Além da execução no R, serão discutidos o planejamento do experimento, a construção do modelo, a análise de variância, a verificação dos pressupostos, as transformações, os métodos não paramétricos e uma introdução aos experimentos fatoriais.
Conceitos fundamentais
Antes de analisar os dados, é necessário identificar os elementos do experimento.
- Unidade experimental: menor unidade que pode receber um tratamento de forma independente. Uma parcela, um animal, um vaso ou um lote podem ser unidades experimentais.
- Variável resposta: característica medida após a aplicação dos tratamentos, como altura, produtividade, teor ou tempo.
- Fator: variável controlada pelo pesquisador. Por exemplo, substrato, cultivar ou dose.
- Nível: modalidade ou valor particular de um fator.
- Tratamento: condição experimental aplicada à unidade. Em um experimento com um fator, os tratamentos são os níveis do fator; em um fatorial, são as combinações dos níveis dos fatores.
- Repetição: aplicação independente de um mesmo tratamento a mais de uma unidade experimental.
- Bloco: grupo de unidades semelhantes quanto a uma fonte conhecida de heterogeneidade.
- Erro experimental: variação entre unidades que receberam a mesma condição experimental e que não foi explicada pelo modelo.
Os três princípios básicos da experimentação são:
- repetição, que permite estimar o erro experimental e aumenta a precisão;
- casualização, que protege a inferência contra vieses sistemáticos;
- controle local, que agrupa unidades semelhantes e retira do erro fontes conhecidas de variação.
Casualização não é sinônimo de escolha aleatória da amostra. Ela se refere à atribuição aleatória dos tratamentos às unidades experimentais. A casualização deve ser realizada e registrada antes da instalação do experimento.
Organização dos dados
Para análise no R, recomenda-se o formato longo: cada linha representa uma unidade experimental e cada coluna representa uma variável. Identificadores de tratamento, bloco, linha e coluna devem ser armazenados como fatores; a resposta deve ser numérica.
| unidade | tratamento | bloco | resposta |
|---|---|---|---|
| 1 | A | 1 | 12,4 |
| 2 | B | 1 | 15,1 |
| 3 | A | 2 | 11,9 |
| 4 | B | 2 | 14,7 |
É importante conservar uma coluna que identifique inequivocamente cada unidade, além da ordem de execução, quando ela existir. Também devem ser mantidos os metadados de casualização, datas, operador, instrumentos e justificativas para qualquer valor ausente. Esses registros permitem distinguir variação experimental de erro de digitação ou falha operacional.
Delineamento inteiramente casualizado
O delineamento inteiramente casualizado (DIC) é apropriado quando as unidades experimentais podem ser consideradas suficientemente homogêneas. Todos os tratamentos concorrem, sem restrições, a todas as unidades experimentais.
Se há \(a\) tratamentos e \(r\) repetições por tratamento, o modelo estatístico é:
\[ Y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \]
em que:
- \(Y_{ij}\) é a observação do tratamento \(i\) na repetição \(j\);
- \(\mu\) é a média geral;
- \(\tau_i\) é o efeito do tratamento \(i\);
- \(\varepsilon_{ij}\) é o erro experimental;
- \(i=1,\ldots,a\) e \(j=1,\ldots,r\).
Para tornar os parâmetros identificáveis, é usual adotar:
\[ \sum_{i=1}^{a}\tau_i=0. \]
Os erros são considerados independentes, com média zero e variância constante:
\[ \varepsilon_{ij}\overset{iid}{\sim}N(0,\sigma^2). \]
A normalidade é necessária para a inferência exata do teste \(F\) em pequenas amostras. Ela se refere aos erros, avaliados por meio dos resíduos, e não necessariamente aos valores brutos da variável resposta.
Hipótese e decomposição da variabilidade
O teste global da ANOVA avalia:
\[ H_0:\tau_1=\tau_2=\cdots=\tau_a=0 \]
contra:
\[ H_1:\text{pelo menos um efeito de tratamento é diferente de zero}. \]
Em um DIC balanceado:
\[ SQ_{Total}=SQ_{Tratamentos}+SQ_{Erro}. \]
Definindo \(\bar{Y}_{i\cdot}\) como a média do tratamento \(i\) e \(\bar{Y}_{\cdot\cdot}\) como a média geral:
\[ SQ_{Tratamentos} =r\sum_{i=1}^{a} (\bar{Y}_{i\cdot}-\bar{Y}_{\cdot\cdot})^2, \]
\[ SQ_{Erro} =\sum_{i=1}^{a}\sum_{j=1}^{r} (Y_{ij}-\bar{Y}_{i\cdot})^2. \]
| Fonte de variação | Graus de liberdade | Quadrado médio | Estatística |
|---|---|---|---|
| Tratamentos | \(a-1\) | \(QM_{Trat}=SQ_{Trat}/(a-1)\) | \(F=QM_{Trat}/QM_E\) |
| Erro | \(a(r-1)\) | \(QM_E=SQ_E/[a(r-1)]\) | |
| Total | \(ar-1\) |
Um teste \(F\) significativo informa que nem todas as médias são iguais, mas não identifica quais diferem. Essa segunda pergunta requer comparações planejadas ou um procedimento de comparações múltiplas.
Contrastes definidos antes da coleta podem ser testados diretamente porque representam hipóteses científicas específicas. Para comparações exploratórias de todas as médias, deve-se controlar a multiplicidade. A exigência de um teste \(F\) global significativo antes do pós-teste é uma convenção frequente, mas o teste de Tukey já controla por si mesmo o erro familiar das comparações duas a duas. O essencial é definir a família de hipóteses e o procedimento antes de escolher os resultados que serão relatados.
Experimento 1: substratos para mudas de acácia
Um experimento foi conduzido no DIC, com quatro repetições, para avaliar a altura de plantas da espécie florestal acácia (Acacia holosericea) submetidas a diferentes substratos para formação de mudas.
| Substrato | Rep. 1 | Rep. 2 | Rep. 3 | Rep. 4 |
|---|---|---|---|---|
| CaC+EB | 18,2 | 16,4 | 20,5 | 19,4 |
| CaC+EA | 12,8 | 14,5 | 10,3 | 11,5 |
| HM | 21,8 | 19,3 | 24,5 | 22,6 |
| Cap+EA | 14,8 | 16,4 | 12,3 | 13,7 |
| Cap+EB | 13,0 | 11,2 | 15,4 | 12,6 |
| Plantmax | 16,0 | 14,2 | 18,1 | 16,4 |
O objetivo é realizar a análise de variância e comparar as médias pelo teste de Tukey, considerando \(\alpha=5\%\).
Uma casualização possível para as 24 unidades seria:
set.seed(2026)
plano_dic <- data.frame(
unidade = 1:24,
substrato = sample(rep(
c("CaC+EB", "CaC+EA", "HM", "Cap+EA", "Cap+EB", "Plantmax"),
each = 4
))
)
plano_dic unidade substrato
1 1 CaC+EB
2 2 CaC+EA
3 3 Cap+EA
4 4 Cap+EA
5 5 HM
6 6 CaC+EB
7 7 Cap+EA
8 8 CaC+EA
9 9 Plantmax
10 10 Cap+EB
11 11 HM
12 12 CaC+EB
13 13 CaC+EA
14 14 CaC+EB
15 15 Plantmax
16 16 Plantmax
17 17 Plantmax
18 18 HM
19 19 Cap+EB
20 20 Cap+EB
21 21 Cap+EA
22 22 Cap+EB
23 23 CaC+EA
24 24 HMO sorteio acima ilustra o planejamento. Os dados observados são criados a seguir.
Arquivo com os dados: dados_1.xlsx
# Leitura do arquivo e conferência da estrutura
dados_1 <- read_xlsx("dados_1.xlsx")
dados_1$substratos <- factor(dados_1$substratos)
str(dados_1)tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
$ repeticao : num [1:24] 1 2 3 4 1 2 3 4 1 2 ...
$ substratos: Factor w/ 6 levels "CaC+EA","CaC+EB",..: 2 2 2 2 1 1 1 1 5 5 ...
$ altura : num [1:24] 18.2 16.4 20.5 19.4 12.8 14.5 10.3 11.5 21.8 19.3 ...head(dados_1)# A tibble: 6 × 3
repeticao substratos altura
<dbl> <fct> <dbl>
1 1 CaC+EB 18.2
2 2 CaC+EB 16.4
3 3 CaC+EB 20.5
4 4 CaC+EB 19.4
5 1 CaC+EA 12.8
6 2 CaC+EA 14.5Antes da ANOVA, devem ser examinadas as observações e as estatísticas descritivas.
boxplot(
altura ~ substratos,
data = dados_1,
xlab = "Substrato",
ylab = "Altura (cm)",
col = "gray90"
)
aggregate(
altura ~ substratos,
data = dados_1,
FUN = function(x) c(
n = length(x),
media = mean(x),
desvio_padrao = sd(x),
mediana = median(x)
)
) substratos altura.n altura.media altura.desvio_padrao altura.mediana
1 CaC+EA 4.000000 12.275000 1.800694 12.150000
2 CaC+EB 4.000000 18.625000 1.755705 18.800000
3 Cap+EA 4.000000 14.300000 1.733974 14.250000
4 Cap+EB 4.000000 13.050000 1.746425 12.800000
5 HM 4.000000 22.050000 2.154840 22.200000
6 Plantmax 4.000000 16.175000 1.600781 16.200000Ajuste pela função lm()
# Modelo linear do DIC
lm_dic <- lm(altura ~ substratos, data = dados_1)
# Estimativas dos parâmetros na parametrização padrão do R
summary(lm_dic)
Call:
lm(formula = altura ~ substratos, data = dados_1)
Residuals:
Min 1Q Median 3Q Max
-2.7500 -1.0437 -0.1125 1.0500 2.4500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.2750 0.9034 13.587 6.66e-11 ***
substratosCaC+EB 6.3500 1.2776 4.970 9.90e-05 ***
substratosCap+EA 2.0250 1.2776 1.585 0.13038
substratosCap+EB 0.7750 1.2776 0.607 0.55169
substratosHM 9.7750 1.2776 7.651 4.60e-07 ***
substratosPlantmax 3.9000 1.2776 3.053 0.00685 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.807 on 18 degrees of freedom
Multiple R-squared: 0.8244, Adjusted R-squared: 0.7756
F-statistic: 16.9 on 5 and 18 DF, p-value: 3.03e-06# Análise de variância
anova_dic <- anova(lm_dic)
anova_dicAnalysis of Variance Table
Response: altura
Df Sum Sq Mean Sq F value Pr(>F)
substratos 5 275.817 55.163 16.898 3.03e-06 ***
Residuals 18 58.762 3.265
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Na tabela produzida por anova():
Dfapresenta os graus de liberdade;Sum Sqapresenta as somas de quadrados;Mean Sqapresenta os quadrados médios;F valueé a razão entre o quadrado médio do efeito e o quadrado médio residual;Pr(>F)é o valor-\(p\) do teste.
O valor-\(p\) é a probabilidade, sob \(H_0\) e sob os pressupostos do modelo, de obter uma estatística tão extrema quanto a observada. Ele não é a probabilidade de \(H_0\) ser verdadeira e não mede sozinho a importância prática do efeito.
O valor de \(F\) compara a variação explicada pelos substratos com a variação entre repetições submetidas ao mesmo substrato. Neste exemplo, o teste global rejeita \(H_0\) (\(F=16{,}90\); valor-\(p<0{,}001\)).
Uma medida descritiva da proporção da variabilidade total associada ao fator é:
\[ \eta^2=\frac{SQ_{Tratamentos}}{SQ_{Total}}. \]
eta2_dic <- anova_dic["substratos", "Sum Sq"] / sum(anova_dic[, "Sum Sq"])
eta2_dic[1] 0.8243691O tamanho de efeito não substitui a estimativa das diferenças e seus intervalos de confiança, mas complementa o valor-\(p\).
As médias ajustadas e seus intervalos de confiança também podem ser obtidos com o modelo. No DIC balanceado, elas coincidem com as médias aritméticas dos grupos.
novos_substratos <- data.frame(
substratos = factor(
levels(dados_1$substratos),
levels = levels(dados_1$substratos)
)
)
intervalos_medias_dic <- cbind(
novos_substratos,
predict(
lm_dic,
newdata = novos_substratos,
interval = "confidence",
level = 1 - alfa
)
)
intervalos_medias_dic substratos fit lwr upr
1 CaC+EA 12.275 10.37701 14.17299
2 CaC+EB 18.625 16.72701 20.52299
3 Cap+EA 14.300 12.40201 16.19799
4 Cap+EB 13.050 11.15201 14.94799
5 HM 22.050 20.15201 23.94799
6 Plantmax 16.175 14.27701 18.07299O coeficiente de variação experimental é uma medida relativa da dispersão residual:
\[ CV(\%)=100\frac{\sqrt{QM_E}}{\bar{Y}_{\cdot\cdot}}. \]
cv_dic <- 100 * sigma(lm_dic) / mean(dados_1$altura)
cv_dic[1] 11.237O CV deve ser comparado com referências da mesma variável, cultura, escala e condição experimental. Classificações universais como “baixo” ou “alto” podem ser enganosas. O CV também perde sentido quando a média está próxima de zero, pode mudar com a transformação e não deve ser comparado entre escalas distintas.
Comparações múltiplas
O teste de Tukey controla a probabilidade de pelo menos um erro do tipo I no conjunto de todas as comparações duas a duas.
aov_dic <- aov(altura ~ substratos, data = dados_1)
TukeyHSD(aov_dic, "substratos", conf.level = 1 - alfa) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = altura ~ substratos, data = dados_1)
$substratos
diff lwr upr p adj
CaC+EB-CaC+EA 6.350 2.2897054 10.4102946 0.0011850
Cap+EA-CaC+EA 2.025 -2.0352946 6.0852946 0.6178206
Cap+EB-CaC+EA 0.775 -3.2852946 4.8352946 0.9891641
HM-CaC+EA 9.775 5.7147054 13.8352946 0.0000060
Plantmax-CaC+EA 3.900 -0.1602946 7.9602946 0.0639232
Cap+EA-CaC+EB -4.325 -8.3852946 -0.2647054 0.0330224
Cap+EB-CaC+EB -5.575 -9.6352946 -1.5147054 0.0042751
HM-CaC+EB 3.425 -0.6352946 7.4852946 0.1281603
Plantmax-CaC+EB -2.450 -6.5102946 1.6102946 0.4235416
Cap+EB-Cap+EA -1.250 -5.3102946 2.8102946 0.9189774
HM-Cap+EA 7.750 3.6897054 11.8102946 0.0001241
Plantmax-Cap+EA 1.875 -2.1852946 5.9352946 0.6876671
HM-Cap+EB 9.000 4.9397054 13.0602946 0.0000185
Plantmax-Cap+EB 3.125 -0.9352946 7.1852946 0.1925714
Plantmax-HM -5.875 -9.9352946 -1.8147054 0.0025987O pacote ExpDes apresenta a ANOVA e o agrupamento de médias em uma única rotina.
crd(
treat = dados_1$substratos,
resp = dados_1$altura,
mcomp = "tukey"
)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 5 275.82 55.163 16.898 3.0304e-06
Residuals 18 58.76 3.265
Total 23 334.58
------------------------------------------------------------------------
CV = 11.24 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.1502373
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.9980527
According to the test of bartlett at 5% of significance, residuals can be considered homocedastic.
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a HM 22.05
ab CaC+EB 18.625
bc Plantmax 16.175
c Cap+EA 14.3
c Cap+EB 13.05
c CaC+EA 12.275
------------------------------------------------------------------------Quando houver interesse no método de agrupamento Scott–Knott:
crd(
treat = dados_1$substratos,
resp = dados_1$altura,
mcomp = "sk"
)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 5 275.82 55.163 16.898 3.0304e-06
Residuals 18 58.76 3.265
Total 23 334.58
------------------------------------------------------------------------
CV = 11.24 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.1502373
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.9980527
According to the test of bartlett at 5% of significance, residuals can be considered homocedastic.
------------------------------------------------------------------------
Scott-Knott test
------------------------------------------------------------------------
Groups Treatments Means
1 a HM 22.050
2 b CaC+EB 18.625
3 c Plantmax 16.175
4 c Cap+EA 14.300
5 c Cap+EB 13.050
6 c CaC+EA 12.275
------------------------------------------------------------------------Os procedimentos respondem a perguntas diferentes. Tukey testa todas as diferenças duas a duas; Scott–Knott procura formar grupos de médias. A escolha deve ser definida pelo objetivo da análise e, preferencialmente, antes de olhar os resultados.
Delineamento em blocos casualizados
O delineamento em blocos casualizados (DBC) é indicado quando existe uma fonte conhecida de heterogeneidade. Cada bloco contém todos os tratamentos e a casualização é feita separadamente dentro de cada bloco.
O modelo é:
\[ Y_{ij}=\mu+\tau_i+\beta_j+\varepsilon_{ij}, \]
em que \(\tau_i\) é o efeito do tratamento \(i\) e \(\beta_j\) é o efeito do bloco \(j\). Adotam-se, por exemplo, as restrições:
\[ \sum_{i=1}^{a}\tau_i=0 \qquad\text{e}\qquad \sum_{j=1}^{b}\beta_j=0. \]
O modelo é aditivo: não há termo tratamento \(\times\) bloco. Com uma observação por combinação, essa interação não pode ser separada do erro. Por isso, a resposta relativa dos tratamentos deve ser aproximadamente consistente entre os blocos.
Nos exemplos deste material, tratamentos representam níveis específicos de interesse e são tratados como efeitos fixos. Blocos podem ser considerados fixos quando o interesse se restringe aos blocos observados, ou aleatórios quando representam uma população de condições. Em um DBC completo e balanceado, o teste dos tratamentos baseado no erro intrablocos coincide nas formulações usuais de efeitos fixos ou aleatórios. Com dados ausentes, estruturas hierárquicas ou interesse em componentes de variância, um modelo misto deve representar explicitamente a natureza aleatória dos blocos.
Para o DBC completo e balanceado:
\[ SQ_{Total}=SQ_{Tratamentos}+SQ_{Blocos}+SQ_{Erro}. \]
| Fonte de variação | Graus de liberdade | Estatística \(F\) |
|---|---|---|
| Tratamentos | \(a-1\) | \(QM_{Trat}/QM_E\) |
| Blocos | \(b-1\) | \(QM_{Blocos}/QM_E\) |
| Erro | \((a-1)(b-1)\) | |
| Total | \(ab-1\) |
O teste dos blocos é secundário. O bloco é uma ferramenta de controle local e não o foco do estudo. A utilidade do bloqueamento deve ser avaliada também pela redução do quadrado médio do erro e pelo conhecimento do experimento.
Planejamento e casualização no DBC
Um bom bloco é internamente homogêneo e pode ser diferente dos demais blocos. Exemplos de critérios de blocagem incluem fertilidade do solo, declividade, peso inicial, lote de matéria-prima, dia, operador ou avaliador.
Com \(a\) tratamentos e \(b\) blocos, há \((a!)^b\) possíveis alocações quando uma permutação independente é sorteada em cada bloco.
set.seed(2026)
tratamentos_dbc <- c("A", "B", "C", "D")
plano_dbc <- do.call(
rbind,
lapply(1:6, function(j) {
data.frame(
bloco = factor(j),
parcela_no_bloco = 1:4,
tratamento = sample(tratamentos_dbc)
)
})
)
plano_dbc bloco parcela_no_bloco tratamento
1 1 1 A
2 1 2 D
3 1 3 C
4 1 4 B
5 2 1 A
6 2 2 C
7 2 3 D
8 2 4 B
9 3 1 D
10 3 2 A
11 3 3 C
12 3 4 B
13 4 1 B
14 4 2 D
15 4 3 C
16 4 4 A
17 5 1 C
18 5 2 A
19 5 3 B
20 5 4 D
21 6 1 D
22 6 2 B
23 6 3 C
24 6 4 AExperimento 2: épocas de corte da alfafa
O experimento avalia quatro épocas de corte, \(A\), \(B\), \(C\) e \(D\), na produtividade de matéria verde de alfafa. Foram formados seis blocos para controlar um gradiente de fertilidade associado à declividade da área.
| Bloco | \(A\) | \(B\) | \(C\) | \(D\) |
|---|---|---|---|---|
| I | 2,89 | 1,58 | 2,29 | 2,56 |
| II | 2,88 | 1,28 | 2,98 | 2,00 |
| III | 1,88 | 1,22 | 1,55 | 1,82 |
| IV | 2,90 | 1,21 | 1,95 | 2,20 |
| V | 2,20 | 1,30 | 1,15 | 1,33 |
| VI | 2,65 | 1,66 | 1,12 | 1,00 |
Os valores representam a produção de matéria verde em kg/parcela.
Arquivo com os dados: dados_2.xlsx
dados_2 <- read_xlsx("dados_2.xlsx")
dados_2$bloco <- factor(dados_2$bloco)
dados_2$epoca <- factor(dados_2$epoca)
str(dados_2)tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
$ bloco : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
$ epoca : Factor w/ 4 levels "A","B","C","D": 1 2 3 4 1 2 3 4 1 2 ...
$ produtividade: num [1:24] 2.89 1.58 2.29 2.56 2.88 1.28 2.98 2 1.88 1.22 ...aggregate(
produtividade ~ epoca,
data = dados_2,
FUN = function(x) c(media = mean(x), desvio_padrao = sd(x))
) epoca produtividade.media produtividade.desvio_padrao
1 A 2.5666667 0.4305191
2 B 1.3750000 0.1944994
3 C 1.8400000 0.7201111
4 D 1.8183333 0.5721684Ajuste do DBC
lm_dbc <- lm(produtividade ~ bloco + epoca, data = dados_2)
summary(lm_dbc)
Call:
lm(formula = produtividade ~ bloco + epoca, data = dados_2)
Residuals:
Min 1Q Median 3Q Max
-0.52583 -0.24000 -0.03125 0.23354 0.75500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.9967 0.2543 11.784 5.54e-09 ***
bloco2 -0.0450 0.2936 -0.153 0.880242
bloco3 -0.7125 0.2936 -2.426 0.028318 *
bloco4 -0.2650 0.2936 -0.902 0.381063
bloco5 -0.8350 0.2936 -2.844 0.012323 *
bloco6 -0.7225 0.2936 -2.461 0.026486 *
epocaB -1.1917 0.2397 -4.970 0.000168 ***
epocaC -0.7267 0.2397 -3.031 0.008425 **
epocaD -0.7483 0.2397 -3.121 0.007006 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4153 on 15 degrees of freedom
Multiple R-squared: 0.7341, Adjusted R-squared: 0.5923
F-statistic: 5.176 on 8 and 15 DF, p-value: 0.003083anova_dbc <- anova(lm_dbc)
anova_dbcAnalysis of Variance Table
Response: produtividade
Df Sum Sq Mean Sq F value Pr(>F)
bloco 5 2.7590 0.55179 3.1999 0.036559 *
epoca 3 4.3820 1.46068 8.4706 0.001572 **
Residuals 15 2.5866 0.17244
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O teste de tratamentos indica diferença entre as épocas (\(F=8{,}47\); valor-\(p=0{,}0016\)). O efeito de blocos também é detectável, o que é coerente com a presença do gradiente usado para formar os blocos.
aov_dbc <- aov(produtividade ~ bloco + epoca, data = dados_2)
TukeyHSD(aov_dbc, "epoca", conf.level = 1 - alfa) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = produtividade ~ bloco + epoca, data = dados_2)
$epoca
diff lwr upr p adj
B-A -1.19166667 -1.8826637 -0.50066968 0.0008626
C-A -0.72666667 -1.4176637 -0.03566968 0.0377333
D-A -0.74833333 -1.4393303 -0.05733634 0.0317423
C-B 0.46500000 -0.2259970 1.15599699 0.2537403
D-B 0.44333333 -0.2476637 1.13433032 0.2900919
D-C -0.02166667 -0.7126637 0.66933032 0.9997229rbd(
treat = dados_2$epoca,
block = dados_2$bloco,
resp = dados_2$produtividade,
mcomp = "tukey"
)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 3 4.3820 1.46068 8.4706 0.001572
Block 5 2.7590 0.55179 3.1999 0.036559
Residuals 15 2.5866 0.17244
Total 23 9.7276
------------------------------------------------------------------------
CV = 21.86 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.7947678
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.6274363
According to the test of oneillmathews at 5% of significance, the variances can be considered homocedastic.
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a A 2.566667
b C 1.84
b D 1.818333
b B 1.375
------------------------------------------------------------------------O teste de Duncan é apresentado apenas para mostrar a sintaxe disponível no pacote. Ele é menos conservador que Tukey e pode aumentar o risco de resultados falso-positivos quando muitas médias são comparadas.
rbd(
treat = dados_2$epoca,
block = dados_2$bloco,
resp = dados_2$produtividade,
mcomp = "duncan"
)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 3 4.3820 1.46068 8.4706 0.001572
Block 5 2.7590 0.55179 3.1999 0.036559
Residuals 15 2.5866 0.17244
Total 23 9.7276
------------------------------------------------------------------------
CV = 21.86 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.7947678
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.6274363
According to the test of oneillmathews at 5% of significance, the variances can be considered homocedastic.
------------------------------------------------------------------------
Duncan's test
------------------------------------------------------------------------
Groups Treatments Means
a A 2.566667
b C 1.84
b D 1.818333
b B 1.375
------------------------------------------------------------------------Eficiência do bloqueamento
Uma comparação descritiva consiste em confrontar o quadrado médio residual do modelo que ignora os blocos com o quadrado médio residual do DBC:
\[ ER(\%)=100\frac{QM_{E,\,DIC}}{QM_{E,\,DBC}}. \]
lm_sem_blocos <- lm(produtividade ~ epoca, data = dados_2)
qm_erro_sem_blocos <- anova(lm_sem_blocos)["Residuals", "Mean Sq"]
qm_erro_dbc <- anova_dbc["Residuals", "Mean Sq"]
eficiencia_relativa <- 100 * qm_erro_sem_blocos / qm_erro_dbc
eficiencia_relativa[1] 154.9969Neste conjunto, a razão é aproximadamente \(155\%\). Essa medida é uma descrição dos dados observados e não autoriza escolher o delineamento depois da coleta. O DBC deve ser definido no planejamento e analisado como foi casualizado.
Delineamento em quadrado latino
O delineamento em quadrado latino (DQL) controla duas fontes de heterogeneidade, representadas por linhas e colunas. Se há \(p\) tratamentos, o experimento deve ter \(p\) linhas, \(p\) colunas e \(p^2\) observações. Cada tratamento aparece uma vez em cada linha e uma vez em cada coluna.
O modelo aditivo é:
\[ Y_{ijk}=\mu+\rho_i+\kappa_j+\tau_k+\varepsilon_{ijk}, \]
em que \(\rho_i\) é o efeito da linha, \(\kappa_j\) é o efeito da coluna e \(\tau_k\) é o efeito do tratamento. Como o tratamento aplicado a cada célula é determinado pelo quadrado latino, linha e coluna localizam a observação.
Adotam-se as restrições:
\[ \sum_i\rho_i=0,\qquad \sum_j\kappa_j=0,\qquad \sum_k\tau_k=0. \]
| Fonte de variação | Graus de liberdade | Estatística \(F\) |
|---|---|---|
| Linhas | \(p-1\) | \(QM_L/QM_E\) |
| Colunas | \(p-1\) | \(QM_C/QM_E\) |
| Tratamentos | \(p-1\) | \(QM_{Trat}/QM_E\) |
| Erro | \((p-1)(p-2)\) | |
| Total | \(p^2-1\) |
O DQL fornece poucos graus de liberdade residuais quando \(p\) é pequeno. Além disso, o modelo supõe aditividade: não há graus de liberdade para estimar as interações tratamento \(\times\) linha ou tratamento \(\times\) coluna. Essa limitação deve ser considerada no planejamento e na interpretação.
Casualização no DQL
Não se sorteiam os tratamentos independentemente nas \(p^2\) células, pois isso destruiria a estrutura latina. Parte-se de um quadrado válido e realizam-se, de forma aleatória:
- uma permutação das linhas;
- uma permutação das colunas;
- uma permutação dos rótulos dos tratamentos;
- opcionalmente, uma transposição do quadrado.
set.seed(2026)
p <- 5
quadrado_base <- outer(0:(p - 1), 0:(p - 1), function(i, j) {
LETTERS[(i + j) %% p + 1]
})
ordem_linhas <- sample(1:p)
ordem_colunas <- sample(1:p)
novos_rotulos <- sample(LETTERS[1:p])
quadrado_dql <- quadrado_base[ordem_linhas, ordem_colunas]
quadrado_dql <- matrix(
novos_rotulos[match(quadrado_dql, LETTERS[1:p])],
nrow = p
)
quadrado_dql [,1] [,2] [,3] [,4] [,5]
[1,] "E" "C" "D" "A" "B"
[2,] "C" "A" "B" "D" "E"
[3,] "B" "E" "A" "C" "D"
[4,] "A" "D" "E" "B" "C"
[5,] "D" "B" "C" "E" "A" Experimento 3: formulações de um propulsor
Um engenheiro estuda cinco formulações de um propulsor. Os lotes de matéria-prima constituem as linhas e os operadores constituem as colunas. Cada formulação é avaliada uma vez em cada lote e uma vez por cada operador.
| Lote/Operador | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| 1 | A=24 | B=20 | C=19 | D=24 | E=24 |
| 2 | B=17 | C=24 | D=30 | E=27 | A=36 |
| 3 | C=18 | D=38 | E=26 | A=27 | B=21 |
| 4 | D=26 | E=31 | A=26 | B=23 | C=22 |
| 5 | E=22 | A=30 | B=20 | C=29 | D=31 |
Arquivo com os dados: dados_3.xlsx
dados_3 <- read_xlsx("dados_3.xlsx")
dados_3$lote <- factor(dados_3$lote)
dados_3$operador <- factor(dados_3$operador)
dados_3$formulacoes <- factor(dados_3$formulacoes)
str(dados_3)tibble [25 × 4] (S3: tbl_df/tbl/data.frame)
$ lote : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
$ operador : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
$ formulacoes: Factor w/ 5 levels "A","B","C","D",..: 1 2 3 4 5 2 3 4 5 1 ...
$ velocidade : num [1:25] 24 17 18 26 22 20 24 38 31 30 ...# Cada formulação deve aparecer uma vez em cada lote e em cada operador
xtabs(~ lote + formulacoes, data = dados_3) formulacoes
lote A B C D E
1 1 1 1 1 1
2 1 1 1 1 1
3 1 1 1 1 1
4 1 1 1 1 1
5 1 1 1 1 1xtabs(~ operador + formulacoes, data = dados_3) formulacoes
operador A B C D E
1 1 1 1 1 1
2 1 1 1 1 1
3 1 1 1 1 1
4 1 1 1 1 1
5 1 1 1 1 1Ajuste do DQL
lm_dql <- lm(
velocidade ~ lote + operador + formulacoes,
data = dados_3
)
summary(lm_dql)
Call:
lm(formula = velocidade ~ lote + operador + formulacoes, data = dados_3)
Residuals:
Min 1Q Median 3Q Max
-3.2 -1.2 -0.2 1.0 5.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.400 2.355 9.087 9.98e-07 ***
lote2 4.600 2.066 2.227 0.04586 *
lote3 3.800 2.066 1.840 0.09067 .
lote4 3.400 2.066 1.646 0.12568
lote5 4.200 2.066 2.033 0.06475 .
operador2 7.200 2.066 3.486 0.00450 **
operador3 2.800 2.066 1.356 0.20021
operador4 4.600 2.066 2.227 0.04586 *
operador5 5.400 2.066 2.614 0.02262 *
formulacoesB -8.400 2.066 -4.067 0.00156 **
formulacoesC -6.200 2.066 -3.002 0.01103 *
formulacoesD 1.200 2.066 0.581 0.57203
formulacoesE -2.600 2.066 -1.259 0.23207
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.266 on 12 degrees of freedom
Multiple R-squared: 0.8107, Adjusted R-squared: 0.6213
F-statistic: 4.281 on 12 and 12 DF, p-value: 0.008855anova_dql <- anova(lm_dql)
anova_dqlAnalysis of Variance Table
Response: velocidade
Df Sum Sq Mean Sq F value Pr(>F)
lote 4 68 17.000 1.5938 0.239059
operador 4 150 37.500 3.5156 0.040373 *
formulacoes 4 330 82.500 7.7344 0.002537 **
Residuals 12 128 10.667
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O efeito das formulações é significativo (\(F=7{,}73\); valor-\(p=0{,}0025\)), após controlar as diferenças entre lotes e entre operadores. Neste exemplo, operadores também explicam uma parcela relevante da variação, enquanto o efeito de lotes não é detectável ao nível de \(5\%\).
aov_dql <- aov(
velocidade ~ lote + operador + formulacoes,
data = dados_3
)
TukeyHSD(aov_dql, "formulacoes", conf.level = 1 - alfa) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = velocidade ~ lote + operador + formulacoes, data = dados_3)
$formulacoes
diff lwr upr p adj
B-A -8.4 -14.9839317 -1.8160683 0.0110827
C-A -6.2 -12.7839317 0.3839317 0.0684350
D-A 1.2 -5.3839317 7.7839317 0.9754380
E-A -2.6 -9.1839317 3.9839317 0.7194121
C-B 2.2 -4.3839317 8.7839317 0.8204614
D-B 9.6 3.0160683 16.1839317 0.0041583
E-B 5.8 -0.7839317 12.3839317 0.0944061
D-C 7.4 0.8160683 13.9839317 0.0254304
E-C 3.6 -2.9839317 10.1839317 0.4461852
E-D -3.8 -10.3839317 2.7839317 0.3966727latsd(
treat = dados_3$formulacoes,
row = dados_3$lote,
column = dados_3$operador,
resp = dados_3$velocidade,
mcomp = "tukey"
)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 4 330 82.500 7.7344 0.002537
Row 4 68 17.000 1.5938 0.239059
Column 4 150 37.500 3.5156 0.040373
Residuals 12 128 10.667
Total 24 676
------------------------------------------------------------------------
CV = 12.86 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.05810583
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a D 29.8
ab A 28.6
abc E 26
bc C 22.4
c B 20.2
------------------------------------------------------------------------Dados ausentes e desbalanceamento
As fórmulas clássicas apresentadas pressupõem dados completos e balanceados. Em um DIC, números diferentes de repetições ainda podem ser analisados por modelo linear, mas as médias têm precisões diferentes e a ordem dos termos pode importar em modelos com vários fatores não ortogonais.
No DBC, uma parcela perdida quebra a ortogonalidade entre tratamentos e blocos. No DQL, a perda de uma única observação é particularmente séria porque o delineamento já possui poucos graus de liberdade residuais. Nesses casos:
- não preencha o valor ausente pela média antes de ajustar o modelo;
- registre a causa da perda e verifique se ela pode depender do tratamento;
- use o modelo linear com os dados observados e métodos de somas de quadrados apropriados à hipótese;
- ajuste os graus de liberdade e as comparações de médias;
- não use rotinas que pressupõem balanceamento sem verificar seu comportamento;
- apresente uma análise de sensibilidade quando o mecanismo de ausência for duvidoso.
Em situações não ortogonais, os testes sequenciais de anova(lm(...)) dependem da ordem dos termos. A hipótese científica deve ser definida antes da escolha de somas de quadrados do tipo I, II ou III. O balanceamento planejado continua sendo a forma mais simples de preservar precisão e interpretação.
Análise de resíduos
A validade da ANOVA depende do modelo, do delineamento e do mecanismo de casualização. Os principais pressupostos sobre os erros são:
- média zero;
- independência;
- variância constante;
- distribuição aproximadamente normal para a inferência exata em pequenas amostras;
- aditividade e especificação correta do modelo.
O diagnóstico deve ser realizado depois do ajuste e com os resíduos. Testar a normalidade separadamente em cada tratamento, com apenas três ou quatro observações, tem pouca utilidade.
Tipos de resíduos
O resíduo bruto é:
\[ e_i=Y_i-\widehat{Y}_i. \]
O resíduo de Pearson no modelo linear homocedástico é:
\[ d_i=\frac{e_i}{\widehat{\sigma}}. \]
Como a variância de \(e_i\) depende da alavancagem \(h_{ii}\), o resíduo studentizado internamente é:
\[ r_i=\frac{e_i}{\widehat{\sigma}\sqrt{1-h_{ii}}}. \]
O resíduo studentizado externamente usa uma estimativa de \(\sigma\) calculada sem a observação \(i\) e é especialmente útil para investigar valores discrepantes.
diagnostico_dic <- data.frame(
observado = dados_1$altura,
ajustado = fitted(lm_dic),
residuo_bruto = residuals(lm_dic),
residuo_pearson = residuals(lm_dic) / sigma(lm_dic),
residuo_padronizado = rstandard(lm_dic),
residuo_studentizado = rstudent(lm_dic),
alavancagem = hatvalues(lm_dic),
distancia_cook = cooks.distance(lm_dic)
)
diagnostico_dic observado ajustado residuo_bruto residuo_pearson residuo_padronizado
1 18.2 18.625 -0.425 -0.23522044 -0.27160916
2 16.4 18.625 -2.225 -1.23144816 -1.42195386
3 20.5 18.625 1.875 1.03773722 1.19827572
4 19.4 18.625 0.775 0.42893138 0.49528730
5 12.8 12.275 0.525 0.29056642 0.33551720
6 14.5 12.275 2.225 1.23144816 1.42195386
7 10.3 12.275 -1.975 -1.09308320 -1.26218376
8 11.5 12.275 -0.775 -0.42893138 -0.49528730
9 21.8 22.050 -0.250 -0.13836496 -0.15977010
10 19.3 22.050 -2.750 -1.52201459 -1.75747106
11 24.5 22.050 2.450 1.35597663 1.56574695
12 22.6 22.050 0.550 0.30440292 0.35149421
13 14.8 14.300 0.500 0.27672992 0.31954019
14 16.4 14.300 2.100 1.16226568 1.34206881
15 12.3 14.300 -2.000 -1.10691970 -1.27816077
16 13.7 14.300 -0.600 -0.33207591 -0.38344823
17 13.0 13.050 -0.050 -0.02767299 -0.03195402
18 11.2 13.050 -1.850 -1.02390072 -1.18229871
19 15.4 13.050 2.350 1.30063065 1.50183891
20 12.6 13.050 -0.450 -0.24905693 -0.28758617
21 16.0 16.175 -0.175 -0.09685547 -0.11183907
22 14.2 16.175 -1.975 -1.09308320 -1.26218376
23 18.1 16.175 1.925 1.06541021 1.23022974
24 16.4 16.175 0.225 0.12452847 0.14379309
residuo_studentizado alavancagem distancia_cook
1 -0.26449923 0.25 4.098419e-03
2 -1.46672304 0.25 1.123307e-01
3 1.21393892 0.25 7.977026e-02
4 0.48464652 0.25 1.362831e-02
5 0.32708852 0.25 6.253989e-03
6 1.46672304 0.25 1.123307e-01
7 -1.28479515 0.25 8.850599e-02
8 -0.48464652 0.25 1.362831e-02
9 -0.15537884 0.25 1.418138e-03
10 -1.87652760 0.25 1.715947e-01
11 1.63720353 0.25 1.361980e-01
12 0.34276935 0.25 6.863788e-03
13 0.31142178 0.25 5.672552e-03
14 1.37485575 0.25 1.000638e-01
15 -1.30267064 0.25 9.076083e-02
16 -0.37417606 0.25 8.168475e-03
17 -0.03105461 0.25 5.672552e-05
18 -1.19638030 0.25 7.765724e-02
19 1.56057143 0.25 1.253067e-01
20 -0.28012783 0.25 4.594767e-03
21 -0.10872582 0.25 6.948876e-04
22 -1.28479515 0.25 8.850599e-02
23 1.24924032 0.25 8.408140e-02
24 0.13982209 0.25 1.148692e-03No DIC balanceado, as alavancagens são iguais dentro da estrutura do modelo. Em experimentos desbalanceados ou com covariáveis, elas podem variar de forma importante.
Sequência recomendada de diagnóstico
- Conferir a digitação, as unidades de medida e os níveis dos fatores.
- Verificar se o modelo reproduz o delineamento usado na casualização.
- Examinar resíduos versus valores ajustados.
- Examinar resíduos por tratamento e, quando aplicável, por bloco.
- Construir o gráfico quantil-quantil dos resíduos.
- Examinar resíduos na ordem real de execução ou coleta.
- Investigar observações influentes nos registros experimentais.
- Somente então considerar transformação, outro modelo ou método robusto.
Normalidade dos resíduos
O histograma pode auxiliar, mas é pouco informativo em amostras pequenas. O gráfico quantil-quantil costuma ser mais útil.
res_dic <- rstandard(lm_dic)
par(mfrow = c(1, 2))
hist(
res_dic,
main = "Resíduos padronizados",
xlab = "Resíduo",
col = "gray90"
)
qqnorm(res_dic, main = "Q-Q plot")
qqline(res_dic, col = "red", lwd = 2)
par(mfrow = c(1, 1))Uma verificação descritiva compara as frequências observadas em intervalos da normal padrão. Com poucos resíduos, desvios das proporções teóricas são naturais.
intervalos <- c("-1 a 1", "-1,64 a 1,64", "-1,96 a 1,96")
frequencias_observadas <- c(
mean(abs(res_dic) < 1),
mean(abs(res_dic) < 1.64),
mean(abs(res_dic) < 1.96)
)
data.frame(
intervalo = intervalos,
frequencia_observada = frequencias_observadas,
probabilidade_normal = c(0.68, 0.90, 0.95)
) intervalo frequencia_observada probabilidade_normal
1 -1 a 1 0.5000000 0.68
2 -1,64 a 1,64 0.9583333 0.90
3 -1,96 a 1,96 1.0000000 0.95O teste de Shapiro–Wilk avalia a compatibilidade com a normalidade:
shapiro.test(res_dic)
Shapiro-Wilk normality test
data: res_dic
W = 0.9384, p-value = 0.1502Não rejeitar \(H_0\) não prova normalidade. Em amostras pequenas, o teste tem baixa potência; em amostras grandes, desvios pouco relevantes podem ser detectados. Gráficos, delineamento, equilíbrio dos tamanhos amostrais e sensibilidade das conclusões devem ser considerados conjuntamente.
O teste de Kolmogorov–Smirnov usual não deve ser aplicado como se os parâmetros da normal fossem conhecidos quando média e variância foram estimadas com os próprios resíduos. Nessa situação, a distribuição de referência muda. Além disso, os resíduos não são uma amostra independente comum, pois satisfazem restrições impostas pelo ajuste.
Independência
Independência é principalmente uma consequência do planejamento. Repetições devem ser unidades experimentais distintas, e a casualização deve impedir que tratamentos fiquem confundidos com tempo, posição, operador ou lote.
Quando a ordem de execução é conhecida:
ordem_execucao <- seq_along(res_dic)
plot(
ordem_execucao,
res_dic,
type = "b",
xlab = "Ordem de execução registrada",
ylab = "Resíduo padronizado"
)
abline(h = 0, lty = 2)
No exemplo, a ordem no arquivo não foi informada como ordem real de execução; portanto, o gráfico é apenas uma demonstração de código. Testes como Durbin–Watson, Breusch–Godfrey e Ljung–Box só fazem sentido quando existe uma ordem temporal ou espacial relevante. Nenhum teste corrige pseudorrepetição.
Homogeneidade de variâncias
O primeiro diagnóstico é o gráfico de resíduos versus valores ajustados. Um formato de funil sugere que a variabilidade depende da média.
plot(
fitted(lm_dic),
res_dic,
xlab = "Valor ajustado",
ylab = "Resíduo padronizado"
)
abline(h = 0, lty = 2)
Também é útil comparar a dispersão em cada tratamento:
boxplot(
altura ~ substratos,
data = dados_1,
xlab = "Substrato",
ylab = "Altura (cm)",
col = "gray90"
)
O teste de Bartlett é eficiente sob normalidade, mas muito sensível a desvios desse pressuposto:
bartlett.test(altura ~ substratos, data = dados_1)
Bartlett test of homogeneity of variances
data: altura by substratos
Bartlett's K-squared = 0.27704, df = 5, p-value = 0.9981Uma versão do teste de Levene centrada nas medianas pode ser construída com funções básicas do R:
mediana_grupo <- ave(
dados_1$altura,
dados_1$substratos,
FUN = median
)
desvio_absoluto <- abs(dados_1$altura - mediana_grupo)
anova(lm(desvio_absoluto ~ dados_1$substratos))Analysis of Variance Table
Response: desvio_absoluto
Df Sum Sq Mean Sq F value Pr(>F)
dados_1$substratos 5 0.4737 0.09475 0.0896 0.9929
Residuals 18 19.0425 1.05792 ou usando a função leveneTest() do pacote car:
library(car)
leveneTest(altura ~ substratos, data = dados_1, center = median)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 0.0896 0.9929
18 O teste de Breusch–Pagan avalia se a variância está associada aos preditores do modelo:
bptest(lm_dic)
studentized Breusch-Pagan test
data: lm_dic
BP = 1.0259, df = 5, p-value = 0.9605Resultados de testes de variância devem ser interpretados com os gráficos. Em um DIC balanceado, a ANOVA costuma ser relativamente robusta a heterogeneidade moderada; a situação é mais preocupante com forte desigualdade e tamanhos amostrais diferentes.
Robustez da ANOVA
Robustez não significa ausência de pressupostos. Em delineamentos balanceados, o teste \(F\) costuma tolerar assimetria moderada e pequenas diferenças de variância, especialmente quando há número razoável de repetições. A preocupação aumenta quando ocorrem simultaneamente:
- poucas repetições;
- forte assimetria ou caudas pesadas;
- valores extremos influentes;
- razão muito grande entre as variâncias dos grupos;
- desbalanceamento associado às variâncias;
- dependência entre unidades.
Normalidade e homocedasticidade se referem à distribuição condicional dos erros. Não se deve abandonar automaticamente a ANOVA por um único teste de diagnóstico com valor-\(p\) ligeiramente inferior a \(0{,}05\). É mais informativo comparar a conclusão paramétrica com uma transformação, um método robusto ou um teste de randomização que preserve o delineamento.
Pontos discrepantes e influência
Um valor extremo não deve ser excluído automaticamente. Primeiro é necessário verificar erro de digitação, falha de equipamento, troca de tratamento ou uma resposta biologicamente legítima.
limite_cook <- 4 / nrow(dados_1)
subset(
diagnostico_dic,
abs(residuo_studentizado) > 2 | distancia_cook > limite_cook
) observado ajustado residuo_bruto residuo_pearson residuo_padronizado
10 19.3 22.05 -2.75 -1.522015 -1.757471
residuo_studentizado alavancagem distancia_cook
10 -1.876528 0.25 0.1715947Os limites \(|r_i|>2\) e \(D_i>4/n\) são regras de triagem, não testes automáticos de exclusão. Uma análise de sensibilidade, com e sem a observação justificada, pode mostrar se a conclusão depende de um único ponto.
Aditividade e especificação do modelo
Quando o tratamento é categórico, não existe uma exigência de “linearidade da resposta em relação ao tratamento” como haveria para uma covariável contínua. A questão correta é se a estrutura aditiva do modelo é adequada.
No DBC e no DQL, interações com blocos, linhas ou colunas não são estimadas com uma observação por célula. Perfis muito cruzados, padrões espaciais ou resíduos sistemáticos podem indicar não aditividade. A solução pode envolver outro delineamento, repetição dentro das células, covariáveis, modelos mistos ou uma transformação. Não se deve acrescentar termos que o experimento não permite estimar.
Diagnóstico dos três modelos
A função plot() aplicada a um modelo linear produz quatro gráficos úteis.
par(mfrow = c(2, 2))
plot(lm_dic)
par(mfrow = c(1, 1))O mesmo procedimento deve ser repetido para o DBC e o DQL:
# Exemplos para execução durante a aula:
# par(mfrow = c(2, 2))
# plot(lm_dbc)
# plot(lm_dql)
# par(mfrow = c(1, 1))Transformação de variáveis
Uma transformação substitui \(Y\) por \(g(Y)\) e ajusta o modelo nessa nova escala. Ela pode estabilizar variâncias, reduzir assimetria, tornar efeitos mais aditivos e diminuir a influência de valores extremos.
Transformar não corrige falta de independência, pseudorrepetição, tratamento aplicado à unidade errada, variável omitida ou delineamento inadequado. A escolha deve ter justificativa estatística e científica.
Relação entre média e variância
Suponha que o desvio-padrão seja proporcional a uma potência da média:
\[ \sigma_Y\propto\mu^\alpha. \]
Para uma transformação potência \(g(Y)=Y^\lambda\), o método delta fornece, aproximadamente:
\[ \operatorname{Var}[g(Y)] \approx [g'(\mu)]^2\operatorname{Var}(Y). \]
Como \(g'(\mu)=\lambda\mu^{\lambda-1}\), a variância transformada será aproximadamente constante quando:
\[ \lambda=1-\alpha. \]
Esse resultado explica por que algumas transformações aparecem com frequência.
| Situação aproximada | Relação média–variância | Transformação tradicional |
|---|---|---|
| variância constante | \(\sigma_Y\propto 1\) | nenhuma |
| contagens de Poisson | \(\sigma_Y\propto\mu^{1/2}\) | \(\sqrt{Y}\) ou \(\sqrt{Y+c}\) |
| erro multiplicativo/log-normal | \(\sigma_Y\propto\mu\) | \(\log(Y)\) |
| forte aumento da dispersão | \(\sigma_Y\propto\mu^{3/2}\) | \(1/\sqrt{Y}\) |
| proporção binomial | \(\operatorname{Var}(P)=p(1-p)/n\) | \(\arcsin(\sqrt{P})\) |
Para contagens pequenas, constantes como \(3/8\) podem reduzir viés na transformação raiz. Para proporções, modelos binomiais são frequentemente mais informativos que a transformação angular, pois respeitam o número de ensaios e permitem inferência diretamente na escala de probabilidade.
Família de Box–Cox
Para respostas estritamente positivas:
\[ Y^{(\lambda)}= \begin{cases} \dfrac{Y^\lambda-1}{\lambda}, & \lambda\neq 0,\\[6pt] \log(Y), & \lambda=0. \end{cases} \]
O parâmetro \(\lambda\) é escolhido pelo perfil da log-verossimilhança do modelo, incluindo o ajuste jacobiano da transformação. Valores usuais são:
- \(\lambda=1\): nenhuma transformação;
- \(\lambda=1/2\): raiz quadrada;
- \(\lambda=0\): logaritmo;
- \(\lambda=-1/2\): inverso da raiz;
- \(\lambda=-1\): inverso.
Não é necessário usar mecanicamente o valor exato que maximiza o perfil. Um valor simples dentro do intervalo de confiança costuma facilitar a comunicação. Depois da escolha, o diagnóstico deve ser repetido.
Experimento 4: dados adaptados para transformação
O exemplo do DIC foi modificado para criar assimetria e heterogeneidade mais evidentes. Essa adaptação é didática e não substitui os dados originais do Experimento 1.
Arquivo com os dados: dados_4.xlsx
dados_4 <- read_xlsx("dados_4.xlsx")
dados_4$substratos <- factor(dados_4$substratos)
lm_dic2 <- lm(altura ~ substratos, data = dados_4)
anova(lm_dic2)Analysis of Variance Table
Response: altura
Df Sum Sq Mean Sq F value Pr(>F)
substratos 5 275.83 55.167 3.2852 0.02783 *
Residuals 18 302.26 16.792
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow = c(1, 2))
plot(
fitted(lm_dic2),
rstandard(lm_dic2),
xlab = "Valor ajustado",
ylab = "Resíduo padronizado"
)
abline(h = 0, lty = 2)
qqnorm(rstandard(lm_dic2))
qqline(rstandard(lm_dic2), col = "red", lwd = 2)
par(mfrow = c(1, 1))
shapiro.test(rstandard(lm_dic2))
Shapiro-Wilk normality test
data: rstandard(lm_dic2)
W = 0.88827, p-value = 0.01224bartlett.test(altura ~ substratos, data = dados_4)
Bartlett test of homogeneity of variances
data: altura by substratos
Bartlett's K-squared = 14.889, df = 5, p-value = 0.01085Usando a função crd do pacote ExpDes:
crd(
treat = dados_4$substratos,
resp = dados_4$altura,
mcomp = "tukey"
)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 5 275.83 55.167 3.2852 0.027833
Residuals 18 302.26 16.792
Total 23 578.10
------------------------------------------------------------------------
CV = 25.55 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.01224471
WARNING: at 5% of significance, residuals can not be considered normal!
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.01084858
WARNING: at 5% of significance, residuals can not be considered homocedastic!
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a HM 22.05
ab CaC+EB 18.625
ab Plantmax 15.925
ab Cap+EA 14.3
ab Cap+EB 13.05
b CaC+EA 12.275
------------------------------------------------------------------------O exemplo apresenta evidências de não normalidade e de variâncias diferentes. O perfil de Box–Cox é obtido por:
perfil_boxcox <- boxcox(
lm_dic2,
lambda = seq(-3, 3, by = 0.05),
plotit = FALSE
)
lambda_estimado <- perfil_boxcox$x[which.max(perfil_boxcox$y)]
lambda_estimado[1] -1.1boxcox(lm_dic2, lambda = seq(-3, 3, by = 0.05))
O máximo do perfil fica próximo de \(-1{,}1\), mas o logaritmo é usado a seguir para mostrar uma transformação simples, interpretável e capaz de melhorar consideravelmente os diagnósticos. A decisão final deve considerar o perfil completo, o intervalo de confiança e o objetivo inferencial.
dados_4$altura_log <- log(dados_4$altura)
lm_dic_log <- lm(altura_log ~ substratos, data = dados_4)
anova(lm_dic_log)Analysis of Variance Table
Response: altura_log
Df Sum Sq Mean Sq F value Pr(>F)
substratos 5 0.84929 0.169858 3.4324 0.02365 *
Residuals 18 0.89075 0.049486
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow = c(1, 2))
plot(
fitted(lm_dic_log),
rstandard(lm_dic_log),
xlab = "Valor ajustado na escala log",
ylab = "Resíduo padronizado"
)
abline(h = 0, lty = 2)
qqnorm(rstandard(lm_dic_log))
qqline(rstandard(lm_dic_log), col = "red", lwd = 2)
par(mfrow = c(1, 1))
shapiro.test(rstandard(lm_dic_log))
Shapiro-Wilk normality test
data: rstandard(lm_dic_log)
W = 0.94999, p-value = 0.2708bartlett.test(altura_log ~ substratos, data = dados_4)
Bartlett test of homogeneity of variances
data: altura_log by substratos
Bartlett's K-squared = 9.2157, df = 5, p-value = 0.1008No exemplo, o valor-\(p\) de Shapiro–Wilk passa de aproximadamente \(0{,}012\) para \(0{,}271\), e o valor-\(p\) de Bartlett passa de aproximadamente \(0{,}011\) para \(0{,}101\). Isso não prova que o modelo transformado é verdadeiro, mas mostra uma melhoria conjunta dos diagnósticos.
crd(
treat = dados_4$substratos,
resp = dados_4$altura_log,
mcomp = "tukey"
)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 5 0.84929 0.169858 3.4324 0.023648
Residuals 18 0.89075 0.049486
Total 23 1.74004
------------------------------------------------------------------------
CV = 8.13 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.2708366
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.1007622
According to the test of bartlett at 5% of significance, residuals can be considered homocedastic.
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a HM 3.027798
ab CaC+EB 2.9211
ab Plantmax 2.749833
ab Cap+EA 2.654726
ab Cap+EB 2.562232
b CaC+EA 2.499521
------------------------------------------------------------------------Interpretação e retrotransformação
Quando o modelo é ajustado a \(\log(Y)\), uma diferença entre médias logarítmicas corresponde, após exponenciação, a uma razão de médias geométricas. Ela não é, em geral, uma diferença entre médias aritméticas na escala original.
medias_log <- aggregate(
altura_log ~ substratos,
data = dados_4,
FUN = mean
)
medias_log$media_geometrica <- exp(medias_log$altura_log)
medias_log substratos altura_log media_geometrica
1 CaC+EA 2.499521 12.17666
2 CaC+EB 2.921100 18.56170
3 Cap+EA 2.654726 14.22109
4 Cap+EB 2.562232 12.96472
5 HM 3.027798 20.65170
6 Plantmax 2.749833 15.64003Como \(E[\exp(X)]\neq\exp[E(X)]\), simplesmente exponenciar uma média logarítmica não estima automaticamente a média aritmética original. Para esse alvo podem ser necessárias correções de viés ou predições marginais apropriadas.
Transformar ou mudar de modelo?
- Para contagens, considere modelos Poisson ou binomial negativo.
- Para proporções com numerador e denominador conhecidos, considere o modelo binomial.
- Para respostas positivas contínuas e assimétricas, considere modelos Gama ou log-normal.
- Para variâncias diferentes conhecidas por grupo, considere modelos com estrutura de variância.
- Para dependência entre medidas, considere modelos mistos ou modelos de correlação.
Transformações e modelos generalizados não são rivais automáticos. A escolha depende da escala em que o efeito científico deve ser interpretado e do mecanismo probabilístico plausível para a resposta.
Testes não paramétricos
Métodos não paramétricos não são métodos sem pressupostos. Eles evitam uma forma paramétrica completa para a distribuição, mas ainda exigem independência ou pareamento correto, escala ordinal ou quantitativa e respeito ao delineamento.
Em testes por postos, observações ordenadas recebem postos \(1,\ldots,N\); empates recebem o posto médio. Transformações estritamente crescentes não alteram os postos.
Testes como Kruskal–Wallis e Friedman avaliam igualdade de distribuições ou de posições relativas. A interpretação estrita como comparação de medianas requer que as distribuições tenham aproximadamente a mesma forma e dispersão.
| Estrutura dos dados | Teste paramétrico | Alternativa por postos |
|---|---|---|
| dois grupos independentes | teste \(t\) | Mann–Whitney–Wilcoxon |
| dois grupos pareados | teste \(t\) pareado | Wilcoxon pareado |
| mais de dois grupos independentes/DIC | ANOVA de um fator | Kruskal–Wallis |
| mais de dois tratamentos em blocos/DBC | ANOVA em blocos | Friedman |
| DQL | ANOVA com linhas e colunas | não há substituto clássico único |
Nesta seção, os experimentos anteriores serão reanalisados por métodos não paramétricos mesmo quando os pressupostos da ANOVA forem razoáveis. A finalidade é comparar perguntas, resultados e limitações.
Reanálise do Experimento 1: Kruskal–Wallis no DIC
Para \(a\) tratamentos, \(N\) observações e soma de postos \(R_i\) no tratamento \(i\), a estatística sem a correção de empates é:
\[ H=\frac{12}{N(N+1)} \sum_{i=1}^{a}\frac{R_i^2}{n_i}-3(N+1). \]
Sob \(H_0\), usa-se aproximadamente \(H\sim\chi^2_{a-1}\). O R aplica a correção necessária para empates.
Com amostras muito pequenas ou muitos empates, a aproximação qui-quadrado pode ser imprecisa. Procedimentos exatos ou de permutação são preferíveis quando disponíveis e quando reproduzem corretamente a aleatorização.
kw_dados_1 <- kruskal.test(
altura ~ substratos,
data = dados_1
)
kw_dados_1
Kruskal-Wallis rank sum test
data: altura by substratos
Kruskal-Wallis chi-squared = 18.312, df = 5, p-value = 0.00258O teste também detecta diferenças entre os substratos (\(H=18{,}31\); valor-\(p=0{,}0026\)). Um tamanho de efeito descritivo é:
\[ \varepsilon_H^2= \frac{H-a+1}{N-a}. \]
n_kw <- nrow(dados_1)
k_kw <- nlevels(dados_1$substratos)
epsilon2_kw <- max(
0,
(as.numeric(kw_dados_1$statistic) - k_kw + 1) / (n_kw - k_kw)
)
epsilon2_kw[1] 0.739547Após um teste global significativo, podem ser realizadas comparações de Wilcoxon duas a duas com ajuste de Holm:
pairwise.wilcox.test(
x = dados_1$altura,
g = dados_1$substratos,
p.adjust.method = "holm",
exact = FALSE
)
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: dados_1$altura and dados_1$substratos
CaC+EA CaC+EB Cap+EA Cap+EB HM
CaC+EB 0.46 - - - -
Cap+EA 0.78 0.46 - - -
Cap+EB 0.94 0.46 0.94 - -
HM 0.46 0.56 0.46 0.46 -
Plantmax 0.48 0.49 0.78 0.48 0.46
P value adjustment method: holm Com apenas quatro observações por substrato, as comparações individuais têm baixa resolução e nenhuma permanece significativa após o ajuste de Holm, apesar do teste global significativo. Esse resultado não é contraditório: o teste global reúne evidência de todas as diferenças, enquanto cada comparação usa menos informação e sofre correção por multiplicidade.
Reanálise do Experimento 2: Friedman no DBC
O teste de Friedman ranqueia os \(a\) tratamentos dentro de cada bloco. Se \(R_{i\cdot}\) é a soma dos postos do tratamento \(i\) nos \(b\) blocos, a estatística sem empates é:
\[ Q=\frac{12}{ba(a+1)} \sum_{i=1}^{a}R_{i\cdot}^2-3b(a+1). \]
Sob \(H_0\), \(Q\) é comparada aproximadamente a \(\chi^2_{a-1}\).
O teste requer um bloco completo: todos os tratamentos devem aparecer em todos os blocos. Observações ausentes não podem ser simplesmente descartadas sem alterar a estrutura. Empates dentro de blocos recebem postos médios e exigem a correção incorporada pela implementação do R.
friedman_dados_2 <- friedman.test(
produtividade ~ epoca | bloco,
data = dados_2
)
friedman_dados_2
Friedman rank sum test
data: produtividade and epoca and bloco
Friedman chi-squared = 10.4, df = 3, p-value = 0.01545O teste detecta diferenças entre épocas (\(Q=10{,}40\); valor-\(p=0{,}0155\)), em concordância qualitativa com a ANOVA do DBC.
Os postos e suas médias podem ser examinados diretamente:
dados_2$posto_no_bloco <- ave(
dados_2$produtividade,
dados_2$bloco,
FUN = rank
)
aggregate(
posto_no_bloco ~ epoca,
data = dados_2,
FUN = function(x) c(soma = sum(x), media = mean(x))
) epoca posto_no_bloco.soma posto_no_bloco.media
1 A 23.000000 3.833333
2 B 9.000000 1.500000
3 C 13.000000 2.166667
4 D 15.000000 2.500000O coeficiente de concordância de Kendall associado ao teste é:
\[ W=\frac{Q}{b(a-1)},\qquad 0\leq W\leq 1. \]
b_friedman <- nlevels(dados_2$bloco)
a_friedman <- nlevels(dados_2$epoca)
kendall_w <- as.numeric(friedman_dados_2$statistic) /
(b_friedman * (a_friedman - 1))
kendall_w[1] 0.5777778Para comparações posteriores, cada par deve preservar o pareamento por bloco.
matriz_dbc <- xtabs(
produtividade ~ bloco + epoca,
data = dados_2
)
pares_epocas <- combn(colnames(matriz_dbc), 2, simplify = FALSE)
pos_friedman <- do.call(
rbind,
lapply(pares_epocas, function(par) {
teste <- wilcox.test(
matriz_dbc[, par[1]],
matriz_dbc[, par[2]],
paired = TRUE,
exact = FALSE
)
data.frame(
epoca_1 = par[1],
epoca_2 = par[2],
p_bruto = teste$p.value
)
})
)
pos_friedman$p_holm <- p.adjust(pos_friedman$p_bruto, method = "holm")
pos_friedman epoca_1 epoca_2 p_bruto p_holm
1 A B 0.03603169 0.2161901
2 A C 0.05917207 0.2366883
3 A D 0.03603169 0.2161901
4 B C 0.20841280 0.4266397
5 B D 0.14221324 0.4266397
6 C D 0.52823278 0.5282328Novamente, o número reduzido de blocos limita a potência das comparações pareadas com correção para multiplicidade.
Teste de randomização restrita no DBC
Outra alternativa usa diretamente a casualização. Sob a hipótese nula forte de nenhum efeito individual de tratamento, os rótulos podem ser permutados dentro de cada bloco. A estatística escolhida será o \(F\) de tratamentos.
Com \(B\) permutações Monte Carlo:
\[ \widehat p= \frac{1+\sum_{s=1}^{B}I(F_s\geq F_{obs})}{B+1}. \]
estatistica_f_dbc <- function(dados) {
ajuste <- lm(produtividade ~ bloco + epoca, data = dados)
anova(ajuste)["epoca", "F value"]
}
set.seed(2026)
B_dbc <- 4999
f_observado_dbc <- estatistica_f_dbc(dados_2)
f_permutado_dbc <- replicate(B_dbc, {
dados_perm <- dados_2
dados_perm$epoca <- ave(
as.character(dados_perm$epoca),
dados_perm$bloco,
FUN = sample
)
dados_perm$epoca <- factor(
dados_perm$epoca,
levels = levels(dados_2$epoca)
)
estatistica_f_dbc(dados_perm)
})
p_randomizacao_dbc <-
(1 + sum(f_permutado_dbc >= f_observado_dbc)) / (B_dbc + 1)
ep_monte_carlo_dbc <- sqrt(
p_randomizacao_dbc * (1 - p_randomizacao_dbc) / (B_dbc + 1)
)
c(
F_observado = f_observado_dbc,
p_monte_carlo = p_randomizacao_dbc,
erro_padrao_monte_carlo = ep_monte_carlo_dbc
) F_observado p_monte_carlo erro_padrao_monte_carlo
8.470588993 0.003200000 0.000798719 Permutar tratamentos livremente entre blocos seria incorreto, pois não reproduziria a casualização do DBC.
Reanálise do Experimento 3: uma alternativa para o DQL
Não existe um teste clássico por postos que seja, ao mesmo tempo, equivalente ao Friedman e controle as duas restrições do DQL. Kruskal–Wallis ignoraria lotes e operadores; Friedman controlaria somente uma dessas fontes.
Como análise didática de sensibilidade, será usado o procedimento de Freedman–Lane:
- ajustar o modelo reduzido com lote e operador;
- permutar os resíduos do modelo reduzido;
- somar os resíduos permutados aos valores ajustados;
- reajustar o modelo completo;
- comparar o \(F\) permutado de formulações com o observado.
Esse teste não exige normalidade, mas requer que os resíduos sejam permutáveis sob a hipótese nula. Ele é uma aproximação baseada no modelo e não deve ser confundido com a reconstrução exata de todas as casualizações possíveis do quadrado latino.
Permutabilidade ainda envolve condições: os erros devem ter distribuição comparável entre as células após o ajuste dos efeitos de linha e coluna. Forte heterogeneidade ou dependência também pode invalidar esse procedimento.
modelo_reduzido_dql <- lm(
velocidade ~ lote + operador,
data = dados_3
)
f_observado_dql <- anova(lm_dql)["formulacoes", "F value"]
set.seed(2026)
B_dql <- 4999
f_permutado_dql <- replicate(B_dql, {
dados_perm <- dados_3
dados_perm$velocidade <- fitted(modelo_reduzido_dql) +
sample(residuals(modelo_reduzido_dql))
ajuste_perm <- lm(
velocidade ~ lote + operador + formulacoes,
data = dados_perm
)
anova(ajuste_perm)["formulacoes", "F value"]
})
p_permutacao_dql <-
(1 + sum(f_permutado_dql >= f_observado_dql)) / (B_dql + 1)
ep_monte_carlo_dql <- sqrt(
p_permutacao_dql * (1 - p_permutacao_dql) / (B_dql + 1)
)
c(
F_observado = f_observado_dql,
p_monte_carlo = p_permutacao_dql,
erro_padrao_monte_carlo = ep_monte_carlo_dql
) F_observado p_monte_carlo erro_padrao_monte_carlo
7.734375000 0.003200000 0.000798719 O resultado Monte Carlo fica próximo do valor-\(p\) paramétrico e sustenta a conclusão de diferença entre formulações. O valor exato pode variar ligeiramente com a semente e com o número de permutações.
Reanálise do Experimento 4: Kruskal–Wallis nos dados assimétricos
O teste de Kruskal–Wallis pode ser aplicado aos dados brutos adaptados:
kw_dados_4 <- kruskal.test(
altura ~ substratos,
data = dados_4
)
kw_dados_4
Kruskal-Wallis rank sum test
data: altura by substratos
Kruskal-Wallis chi-squared = 10.602, df = 5, p-value = 0.05987Neste caso, a ANOVA bruta produz valor-\(p\approx0{,}028\), enquanto o teste de Kruskal–Wallis produz valor-\(p\approx0{,}060\). A diferença é didaticamente importante: testes paramétricos, transformações e testes por postos não respondem exatamente à mesma hipótese e atribuem pesos diferentes às magnitudes e aos valores extremos. Não se deve escolher o método apenas porque ele produz o menor valor-\(p\).
Comparação das estratégias
| Estratégia | Usa magnitudes? | Normalidade dos erros? | Respeita o delineamento? | Interpretação principal |
|---|---|---|---|---|
| ANOVA clássica | sim | para inferência exata pequena | sim | médias |
| ANOVA transformada | sim, em nova escala | aproximadamente | sim | médias transformadas |
| Kruskal–Wallis | somente ordem | não | somente grupos independentes | distribuições/postos |
| Friedman | postos intrablocos | não | DBC completo | distribuições/postos em blocos |
| randomização/permutação | depende da estatística | não | sim, se corretamente restrita | hipótese definida pela permutação |
| modelo generalizado | sim | não para a resposta bruta | sim, se bem especificado | média via função de ligação |
O método não paramétrico não é automaticamente superior quando a normalidade é questionável. O desenho experimental, a pergunta científica e a escala de interpretação continuam sendo centrais.
Introdução aos experimentos fatoriais
Um experimento fatorial estuda dois ou mais fatores simultaneamente. Ele não é, por si só, um delineamento: é uma estrutura de tratamentos. As combinações fatoriais podem ser distribuídas em DIC, DBC, quadrado latino ou outro delineamento adequado.
Em um fatorial completo, todas as combinações dos níveis dos fatores são observadas. Se o fator \(A\) tem \(a\) níveis e o fator \(B\) tem \(b\) níveis, existem:
\[ a\times b \]
combinações de tratamentos. Com \(r\) repetições por combinação, o número total de unidades é:
\[ N=abr. \]
Os fatores são cruzados porque cada nível de \(A\) ocorre com cada nível de \(B\). Isso é diferente de fatores aninhados, nos quais os níveis de um fator existem apenas dentro de um nível específico do outro.
Por que usar um fatorial?
- Estimar os efeitos de vários fatores em um único experimento.
- Usar as observações de todas as combinações para estimar efeitos principais.
- Investigar interação entre fatores.
- Produzir conclusões mais amplas do que experimentos separados de um fator.
- Identificar combinações de níveis que funcionam de maneira particular.
A principal informação adicional é a interação. Há interação quando o efeito de um fator depende do nível do outro.
Se \(\mu_{ij}\) é a média no nível \(i\) de \(A\) e no nível \(j\) de \(B\), a ausência de interação implica que diferenças entre níveis de \(A\) são constantes nos níveis de \(B\). Para dois níveis \(i\) e \(i'\) de \(A\):
\[ \mu_{ij}-\mu_{i'j} = \mu_{ij'}-\mu_{i'j'}. \]
Graficamente, perfis aproximadamente paralelos indicam pouca interação; linhas não paralelas sugerem interação. O gráfico é diagnóstico e descritivo: a decisão inferencial vem do modelo.
Fatorial \(3\times2\)
Considere dois fatores:
- fator \(A\): recipiente, com três níveis (\(R_1\), \(R_2\) e \(R_3\));
- fator \(B\): espécie, com dois níveis (\(E_1\) e \(E_2\)).
O experimento tem \(3\times2=6\) combinações de tratamentos. Com quatro repetições por combinação, são necessárias \(3\times2\times4=24\) unidades experimentais.
Em um DIC, o modelo é:
\[ Y_{ijk} = \mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\varepsilon_{ijk}, \]
em que:
- \(\alpha_i\) é o efeito do nível \(i\) do recipiente;
- \(\beta_j\) é o efeito do nível \(j\) da espécie;
- \((\alpha\beta)_{ij}\) é a interação recipiente \(\times\) espécie;
- \(k=1,\ldots,r\) identifica a repetição.
Restrições usuais são:
\[ \sum_i\alpha_i=0,\qquad \sum_j\beta_j=0, \]
\[ \sum_i(\alpha\beta)_{ij}=0 \quad\text{para cada }j, \qquad \sum_j(\alpha\beta)_{ij}=0 \quad\text{para cada }i. \]
Hipóteses e graus de liberdade
Para o fator \(A\):
\[ H_{0A}:\alpha_1=\alpha_2=\alpha_3=0. \]
Para o fator \(B\):
\[ H_{0B}:\beta_1=\beta_2=0. \]
Para a interação:
\[ H_{0AB}:(\alpha\beta)_{ij}=0 \quad\text{para todos }i,j. \]
| Fonte de variação | Graus de liberdade no \(3\times2\) | Estatística \(F\) |
|---|---|---|
| Recipiente (\(A\)) | \(3-1=2\) | \(QM_A/QM_E\) |
| Espécie (\(B\)) | \(2-1=1\) | \(QM_B/QM_E\) |
| \(A\times B\) | \((3-1)(2-1)=2\) | \(QM_{AB}/QM_E\) |
| Erro | \(3\cdot2\cdot(4-1)=18\) | |
| Total | \(24-1=23\) |
A soma de quadrados das seis combinações é decomposta ortogonalmente, no caso balanceado, em:
\[ SQ_{Tratamentos}=SQ_A+SQ_B+SQ_{AB}. \]
O teste da interação deve ser interpretado primeiro. Se a interação é importante, os efeitos principais são médias de comportamentos diferentes e podem esconder a relação científica relevante. Nesse caso, estudam-se os efeitos simples: o efeito de \(A\) dentro de cada nível de \(B\), ou o efeito de \(B\) dentro de cada nível de \(A\).
Planejamento e casualização do fatorial
As seis combinações devem ser tratadas como condições experimentais distintas e casualizadas às 24 unidades. Cada repetição precisa ser uma unidade independente; medir quatro folhas da mesma planta não cria quatro repetições do tratamento.
set.seed(2026)
combinacoes <- interaction(
rep(c("R1", "R2", "R3"), each = 2),
rep(c("E1", "E2"), times = 3),
sep = ":"
)
plano_fatorial <- data.frame(
unidade = 1:24,
tratamento = sample(rep(as.character(combinacoes), each = 4))
)
plano_fatorial unidade tratamento
1 1 R1:E1
2 2 R1:E2
3 3 R2:E2
4 4 R2:E2
5 5 R2:E1
6 6 R1:E1
7 7 R2:E2
8 8 R1:E2
9 9 R3:E2
10 10 R3:E1
11 11 R2:E1
12 12 R1:E1
13 13 R1:E2
14 14 R1:E1
15 15 R3:E2
16 16 R3:E2
17 17 R3:E2
18 18 R2:E1
19 19 R3:E1
20 20 R3:E1
21 21 R2:E2
22 22 R3:E1
23 23 R1:E2
24 24 R2:E1table(plano_fatorial$tratamento)
R1:E1 R1:E2 R2:E1 R2:E2 R3:E1 R3:E2
4 4 4 4 4 4 Exemplo: recipientes e espécies de eucalipto
O exemplo, adaptado do material de experimentação agrícola disponível no diretório, avalia três recipientes para produção de mudas e duas espécies de eucalipto. A resposta é a altura das mudas, em cm, aos 80 dias.
- \(R_1\): saco plástico pequeno;
- \(R_2\): saco plástico grande;
- \(R_3\): recipiente laminado;
- \(E_1\): Eucalyptus citriodora;
- \(E_2\): Eucalyptus grandis.
| Combinação | Rep. 1 | Rep. 2 | Rep. 3 | Rep. 4 |
|---|---|---|---|---|
| \(R_1E_1\) | 26,2 | 26,0 | 25,0 | 25,4 |
| \(R_1E_2\) | 24,8 | 24,6 | 26,7 | 25,2 |
| \(R_2E_1\) | 25,7 | 26,3 | 25,1 | 26,4 |
| \(R_2E_2\) | 19,6 | 21,1 | 19,0 | 18,6 |
| \(R_3E_1\) | 22,8 | 19,4 | 18,8 | 19,2 |
| \(R_3E_2\) | 19,8 | 21,4 | 22,8 | 21,3 |
recipiente <- factor(rep(c("R1", "R2", "R3"), each = 8))
especie <- factor(rep(rep(c("E1", "E2"), each = 4), times = 3))
repeticao <- factor(rep(1:4, times = 6))
altura_muda <- c(
26.2, 26.0, 25.0, 25.4,
24.8, 24.6, 26.7, 25.2,
25.7, 26.3, 25.1, 26.4,
19.6, 21.1, 19.0, 18.6,
22.8, 19.4, 18.8, 19.2,
19.8, 21.4, 22.8, 21.3
)
dados_fatorial <- data.frame(
repeticao,
recipiente,
especie,
altura_muda
)
writexl::write_xlsx(dados_fatorial, "dados_fatorial_3x2.xlsx")
head(dados_fatorial) repeticao recipiente especie altura_muda
1 1 R1 E1 26.2
2 2 R1 E1 26.0
3 3 R1 E1 25.0
4 4 R1 E1 25.4
5 1 R1 E2 24.8
6 2 R1 E2 24.6# Conferir as quatro repetições de cada combinação
xtabs(~ recipiente + especie, data = dados_fatorial) especie
recipiente E1 E2
R1 4 4
R2 4 4
R3 4 4Arquivo com os dados: dados_fatorial_3x2.xlsx
Exploração das médias e da interação
medias_fatoriais <- with(
dados_fatorial,
tapply(altura_muda, list(recipiente, especie), mean)
)
medias_fatoriais E1 E2
R1 25.650 25.325
R2 25.875 19.575
R3 20.050 21.325with(
dados_fatorial,
interaction.plot(
x.factor = recipiente,
trace.factor = especie,
response = altura_muda,
type = "b",
pch = c(16, 17),
lty = 1,
xlab = "Recipiente",
ylab = "Altura média (cm)",
trace.label = "Espécie"
)
)
Os perfis não são paralelos. Em \(R_1\), as espécies têm médias semelhantes; em \(R_2\), \(E_1\) apresenta altura muito maior que \(E_2\); em \(R_3\), a diferença muda de direção e é menor. Esse padrão sugere interação.
Ajuste do modelo fatorial no R
O operador * inclui os dois efeitos principais e a interação:
lm_fatorial <- lm(
altura_muda ~ recipiente * especie,
data = dados_fatorial
)
summary(lm_fatorial)
Call:
lm(formula = altura_muda ~ recipiente * especie, data = dados_fatorial)
Residuals:
Min 1Q Median 3Q Max
-1.5250 -0.6687 -0.1500 0.4500 2.7500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.6500 0.5663 45.294 < 2e-16 ***
recipienteR2 0.2250 0.8009 0.281 0.782
recipienteR3 -5.6000 0.8009 -6.992 1.58e-06 ***
especieE2 -0.3250 0.8009 -0.406 0.690
recipienteR2:especieE2 -5.9750 1.1326 -5.275 5.13e-05 ***
recipienteR3:especieE2 1.6000 1.1326 1.413 0.175
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.133 on 18 degrees of freedom
Multiple R-squared: 0.8838, Adjusted R-squared: 0.8516
F-statistic: 27.39 on 5 and 18 DF, p-value: 8.043e-08anova_fatorial <- anova(lm_fatorial)
anova_fatorialAnalysis of Variance Table
Response: altura_muda
Df Sum Sq Mean Sq F value Pr(>F)
recipiente 2 92.861 46.430 36.195 4.924e-07 ***
especie 1 19.082 19.082 14.875 0.001155 **
recipiente:especie 2 63.761 31.880 24.853 6.635e-06 ***
Residuals 18 23.090 1.283
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Os três testes são significativos, inclusive a interação (\(F=24{,}85\); valor-\(p<0{,}001\)). Portanto, a interpretação deve se concentrar nos efeitos simples, e não apenas nas médias marginais de recipientes e espécies.
Medidas de tamanho de efeito parcial podem ser calculadas por:
\[ \eta_p^2=\frac{SQ_{efeito}}{SQ_{efeito}+SQ_E}. \]
sq_erro_fatorial <- anova_fatorial["Residuals", "Sum Sq"]
eta2_parcial <- data.frame(
efeito = rownames(anova_fatorial)[1:3],
eta2_parcial = anova_fatorial[1:3, "Sum Sq"] /
(anova_fatorial[1:3, "Sum Sq"] + sq_erro_fatorial)
)
eta2_parcial efeito eta2_parcial
1 recipiente 0.8008639
2 especie 0.4524760
3 recipiente:especie 0.7341419Desdobramento da interação
Para manter o mesmo quadrado médio residual do experimento completo, os efeitos simples serão calculados usando \(QM_E\) e os 18 graus de liberdade residuais do modelo fatorial.
qm_erro_fatorial <- anova_fatorial["Residuals", "Mean Sq"]
gl_erro_fatorial <- anova_fatorial["Residuals", "Df"]
teste_efeito_simples <- function(dados, resposta, fator) {
y <- dados[[resposta]]
g <- droplevels(dados[[fator]])
medias <- tapply(y, g, mean)
tamanhos <- table(g)
sq <- sum(tamanhos * (medias - mean(y))^2)
gl <- nlevels(g) - 1
f <- (sq / gl) / qm_erro_fatorial
c(
SQ = sq,
gl = gl,
F = f,
p_valor = pf(f, gl, gl_erro_fatorial, lower.tail = FALSE)
)
}Efeito de espécie dentro de cada recipiente
efeito_especie_em_recipiente <- do.call(
rbind,
lapply(levels(dados_fatorial$recipiente), function(nivel) {
resultado <- teste_efeito_simples(
subset(dados_fatorial, recipiente == nivel),
resposta = "altura_muda",
fator = "especie"
)
data.frame(recipiente = nivel, t(resultado))
})
)
efeito_especie_em_recipiente$p_holm <- p.adjust(
efeito_especie_em_recipiente$p_valor,
method = "holm"
)
efeito_especie_em_recipiente recipiente SQ gl F p_valor p_holm
1 R1 0.21125 1 0.1646817 6.896656e-01 6.896656e-01
2 R2 79.38000 1 61.8813339 3.112001e-07 9.336003e-07
3 R3 3.25125 1 2.5345388 1.287883e-01 2.575765e-01A diferença entre espécies é pronunciada no recipiente \(R_2\). Em \(R_1\) e \(R_3\), as evidências são insuficientes para afirmar diferença entre espécies ao nível de \(5\%\).
Efeito de recipiente dentro de cada espécie
efeito_recipiente_em_especie <- do.call(
rbind,
lapply(levels(dados_fatorial$especie), function(nivel) {
resultado <- teste_efeito_simples(
subset(dados_fatorial, especie == nivel),
resposta = "altura_muda",
fator = "recipiente"
)
data.frame(especie = nivel, t(resultado))
})
)
efeito_recipiente_em_especie$p_holm <- p.adjust(
efeito_recipiente_em_especie$p_valor,
method = "holm"
)
efeito_recipiente_em_especie especie SQ gl F p_valor p_holm
1 E1 87.12167 2 33.95821 7.776197e-07 1.555239e-06
2 E2 69.50000 2 27.08965 3.730256e-06 3.730256e-06Há diferenças entre recipientes nas duas espécies. Para investigar quais combinações diferem, pode-se comparar as seis médias de tratamento usando o mesmo erro do experimento:
dados_fatorial$combinacao <- interaction(
dados_fatorial$recipiente,
dados_fatorial$especie,
sep = ":"
)
aov_combinacoes <- aov(
altura_muda ~ combinacao,
data = dados_fatorial
)
TukeyHSD(aov_combinacoes, "combinacao", conf.level = 1 - alfa) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = altura_muda ~ combinacao, data = dados_fatorial)
$combinacao
diff lwr upr p adj
R2:E1-R1:E1 0.225 -2.3201851 2.770185 0.9997185
R3:E1-R1:E1 -5.600 -8.1451851 -3.054815 0.0000204
R1:E2-R1:E1 -0.325 -2.8701851 2.220185 0.9983324
R2:E2-R1:E1 -6.075 -8.6201851 -3.529815 0.0000068
R3:E2-R1:E1 -4.325 -6.8701851 -1.779815 0.0004825
R3:E1-R2:E1 -5.825 -8.3701851 -3.279815 0.0000120
R1:E2-R2:E1 -0.550 -3.0951851 1.995185 0.9811892
R2:E2-R2:E1 -6.300 -8.8451851 -3.754815 0.0000041
R3:E2-R2:E1 -4.550 -7.0951851 -2.004815 0.0002705
R1:E2-R3:E1 5.275 2.7298149 7.820185 0.0000444
R2:E2-R3:E1 -0.475 -3.0201851 2.070185 0.9902110
R3:E2-R3:E1 1.275 -1.2701851 3.820185 0.6135909
R2:E2-R1:E2 -5.750 -8.2951851 -3.204815 0.0000143
R3:E2-R1:E2 -4.000 -6.5451851 -1.454815 0.0011258
R3:E2-R2:E2 1.750 -0.7951851 4.295185 0.2914242Essa comparação de todas as 15 diferenças é válida no DIC balanceado, mas pode ser mais ampla que a pergunta científica. Quando a interação é significativa, comparações planejadas de efeitos simples, com controle explícito de multiplicidade, costumam ser mais informativas.
Diagnóstico do modelo fatorial
Os mesmos pressupostos do modelo linear devem ser verificados.
par(mfrow = c(2, 2))
plot(lm_fatorial)
par(mfrow = c(1, 1))
shapiro.test(rstandard(lm_fatorial))
Shapiro-Wilk normality test
data: rstandard(lm_fatorial)
W = 0.9293, p-value = 0.09402Também é importante verificar a dispersão nas seis combinações. Com apenas quatro repetições por célula, testes formais têm pouca potência; o gráfico e o conhecimento do processo devem receber grande peso.
boxplot(
altura_muda ~ combinacao,
data = dados_fatorial,
xlab = "Combinação recipiente:espécie",
ylab = "Altura (cm)",
col = "gray90"
)
Como interpretar uma interação
Uma sequência segura é:
- examinar o teste e o gráfico da interação;
- se a interação for relevante, estimar efeitos simples;
- apresentar médias, diferenças, intervalos de confiança e unidades;
- controlar a multiplicidade das comparações não planejadas;
- evitar uma conclusão geral sobre um fator quando seu efeito muda com o outro.
Neste exemplo, não é adequado afirmar apenas que “\(R_1\) é o melhor recipiente” ou que “\(E_1\) é a melhor espécie”. A conclusão correta descreve as combinações: \(R_1\) produz alturas semelhantes para as duas espécies; \(R_2\) favorece fortemente \(E_1\); e \(R_3\) apresenta médias menores, com diferença pequena entre as espécies.
E se o fatorial for instalado em blocos?
Se houver uma fonte de heterogeneidade, as seis combinações podem ser casualizadas dentro de cada bloco. O modelo passa a incluir o efeito de bloco:
\[ Y_{ijk\ell} =\mu+\rho_\ell+\alpha_i+\beta_j +(\alpha\beta)_{ij}+\varepsilon_{ijk\ell}. \]
No R:
# Estrutura geral para um fatorial 3 x 2 em DBC:
# lm(resposta ~ bloco + fator_A * fator_B, data = dados)Os tratamentos continuam sendo as seis combinações. A casualização deve ocorrer separadamente em cada bloco, e cada bloco deve conter todas as combinações.
Síntese do procedimento de análise
Para qualquer experimento deste material:
- identifique unidade experimental, tratamentos, repetições e restrições de casualização;
- organize os dados no formato longo, uma linha por unidade experimental;
- declare tratamentos, blocos, linhas e colunas como fatores;
- faça a exploração gráfica e descritiva sem substituir a inferência;
- ajuste o modelo correspondente ao delineamento;
- interprete primeiro o teste global ou a interação;
- verifique resíduos, variâncias, independência e observações influentes;
- use comparações múltiplas coerentes com a pergunta;
- considere transformação, GLM ou teste não paramétrico quando justificado;
- relate estimativas, incerteza, unidades e relevância prática.
Erros frequentes
- Confundir subamostras com repetições independentes.
- Ignorar blocos, linhas ou colunas na análise.
- Testar normalidade dos dados brutos em vez dos resíduos.
- Escolher o pós-teste depois de observar qual fornece mais diferenças.
- Aplicar Kruskal–Wallis a dados blocados.
- Aplicar Friedman ao DQL ignorando uma das duas restrições.
- Transformar os dados sem repetir o diagnóstico.
- Retrotransformar médias sem esclarecer o alvo inferencial.
- Interpretar efeitos principais isoladamente quando há interação importante.
- Excluir outliers apenas porque ultrapassaram uma regra numérica.
- Relatar somente valor-\(p\), sem médias, diferenças ou unidades.
Modelo de relato
Um relato adequado deve indicar:
- delineamento e unidade experimental;
- fatores, níveis, número de repetições e método de casualização;
- modelo estatístico e estrutura de erro;
- diagnóstico dos pressupostos e eventuais adaptações;
- estatística, graus de liberdade e valor-\(p\);
- médias ou contrastes com medida de incerteza;
- método de ajuste para multiplicidade;
- escala de análise e forma de retrotransformação, quando aplicável;
- conclusão científica, sem transformar ausência de evidência em prova de igualdade.