Planejamento e Análise de Experimentos I

Computacional

Author
Affiliation

Prof. Dr. Juliano Bortolini

Universidade Federal de Mato Grosso

Bacharelado em Estatística - UFMT

Análise de experimentos no R

# Instalar pacotes, se não estiver instalado:
# install.packages("ExpDes")
# install.packages("lmtest")
# install.packages("MASS")
# install.packages("readxl")

# Carregar os pacotes:
library(ExpDes)
library(lmtest)
library(MASS)
library(readxl)

Delineamento inteiramente casualizado

Modelo estatístico:

\[ Y_{ij} = \mu + \tau_i + \epsilon_{ij}, \] em que \(Y_{ij}\) é a observação referente ao tratamento \(i\) e repetição \(j\), \(i = 1, \cdots, I\), \(j = 1, \cdots, J\), \(\mu\) é uma constante comum a todas as observações, \(\tau_i\) é o efeito do tratamento \(i\) e \(\epsilon_{ij}\) é o erro experimental, que se supõe que seja independente e identicamente distribuído (normal) com média zero e variância \(\sigma^2\).

Experimento

  1. Um experimento foi conduzido no delineamento inteiramente casualizado (DIC), com quatro repetições, para avaliar a altura de plantas da espécie florestal acácia (Acacia holocericea) submetidas a diferentes substratos para formação de mudas. Os tratamentos e os valores de altura de planta (em cm) foram registrados na seguinte tabela:
Substratos 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

Realize a análise de variância com comparações múltiplas de Tukey, considerando um nível de significância de \(5\%\).

Arquivo com os dados: dados_1.xlsx

# Lendo os dados
dados_1 <- read_xlsx("dados_1.xlsx")
dados_1
# A tibble: 24 × 3
   repeticao substratos altura
       <dbl> <chr>       <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.5
 7         3 CaC+EA       10.3
 8         4 CaC+EA       11.5
 9         1 HM           21.8
10         2 HM           19.3
# ℹ 14 more rows
names(dados_1)
[1] "repeticao"  "substratos" "altura"    
# Anova usando a função lm do pacote stats
# A função lm é a função padrão para ajuste de modelos lineares no R

# A variável com os tratamentos precisa ser fator
dados_1$substratos <- as.factor(dados_1$substratos)

# definindo o modelo linear
lm_dic <- lm(altura ~ substratos, data = dados_1)

# estimativas dos parâmetros do modelo linear
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(lm_dic)
Analysis 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 ' ' 1
# ANOVA usando ExpDes
# teste Tukey é o padrão para comparações múltiplas
crd(treat = dados_1$substratos,
    resp = dados_1$altura)
------------------------------------------------------------------------
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 
------------------------------------------------------------------------

Note que há ambiguidade no teste de comparações múltipla usando o teste de Tukey. Uma outra opção para as comparações múltiplas é usar o teste Scott-Knott.

# teste 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
------------------------------------------------------------------------

Outras opções de teste de comparações múltiplas:

# ?crd
# mcomp
# Allows choosing the multiple comparison test; the default is the test of Tukey, however, the options are: the LSD test ('lsd'), the LSD test with Bonferroni protection ('lsdb'), the test of Duncan ('duncan'), the test of Student-Newman-Keuls ('snk'), the test of Scott-Knot ('sk'), the Calinski and Corsten test ('ccF') and bootstrap multiple comparison's test ('ccboot').

Delineamento em blocos casualizados

Modelo estatístico:

\[ Y_{ij} = \mu + \tau_i + \beta_j + \epsilon_{ij}, \] em que \(Y_{ij}\) é a observação referente ao tratamento \(i\) e bloco \(j\), \(i = 1, \cdots, I\), \(j = 1, \cdots, J\), \(\mu\) é uma constante comum a todas as observações, \(\tau_i\) é o efeito do tratamento \(i\), \(\beta_j\) é o efeito do bloco \(j\) e \(\epsilon_{ij}\) é o erro experimental, que se supõe que seja independente e identicamente distribuído (normal) com média zero e variância \(\sigma^2\).

Experimento

  1. Os dados da Tabela 1 referem-se a um ensaio sobre a influência de quatro épocas de corte na produtividade de matéria verde de uma variedade de alfafa. As épocas estudadas foram \(A\), \(B\), \(C\) e \(D\), sendo \(A\) a mais precoce e \(D\) a mais tardia. Foi utilizado o delineamento em Blocos Casualizados para controlar um possível gradiente de fertilidade do solo, já que a área experimental apresentava uma declividade de \(12\%\).

Tabela 1: Produções em \(\text{Kg/parcela}\) de matéria verde de alfafa.

Blocos \(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

Realize a análise de variância com teste Tukey para comparações múltiplas, ao nível de \(5\%\) de probabilidade.

Arquivo com os dados: dados_2.xlsx

# Lendo os dados
dados_2 <- read_xlsx("dados_2.xlsx")
dados_2
# A tibble: 24 × 3
   bloco epoca produtividade
   <dbl> <chr>         <dbl>
 1     1 A              2.89
 2     1 B              1.58
 3     1 C              2.29
 4     1 D              2.56
 5     2 A              2.88
 6     2 B              1.28
 7     2 C              2.98
 8     2 D              2   
 9     3 A              1.88
10     3 B              1.22
# ℹ 14 more rows
names(dados_2)
[1] "bloco"         "epoca"         "produtividade"
# Anova usando a função lm do pacote stats

# As variáveis com os blocos e tratamentos precisam ser fator
dados_2$bloco <- as.factor(dados_2$bloco)
dados_2$epoca <- as.factor(dados_2$epoca)

# definindo o modelo linear
lm_dbc <- lm(produtividade ~ bloco + epoca, data = dados_2)

# estimativas dos parâmetros do modelo linear
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.003083
# Análise de variância
anova(lm_dbc)
Analysis 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 ' ' 1
# ANOVA usando ExpDes
# teste Tukey é o padrão para comparações múltiplas
rbd(treat = dados_2$epoca,
    block = dados_2$bloco,
    resp = dados_2$produtividade)
------------------------------------------------------------------------
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 
------------------------------------------------------------------------
# Anova com o teste Duncan para comparações múltiplas
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 
------------------------------------------------------------------------

Delineamento em quadrado latino

Modelo estatístico:

\[ Y_{ijk} = \mu + \tau_i + L_j + C_k + \epsilon_{ijk}, \] em que \(Y_{ijk}\) é a observação referente ao tratamento \(i\) na linha \(j\) e coluna \(k\), \(i = 1, \cdots, p\), \(j = 1, \cdots, p\), \(k = 1, \cdots, p\), \(\mu\) é uma constante comum a todas as observações, \(\tau_i\) é o efeito do tratamento \(i\), \(L_j\) é o efeito da linha \(j\), \(C_k\) é o efeito da coluna \(k\) e \(\epsilon_{ijk}\) é o erro experimental, que se supõe que seja independente e identicamente distribuído (normal) com média zero e variância \(\sigma^2\).

Experimento

  1. Um engenheiro aeroespacial está estudando os efeitos de cinco formulações diferentes de um propulsor de foguete usado em sistemas de evacuação da tripulação.
  • Variável observada: velocidade de queima das formulações.
  • Cada formulação é misturada a partir de um lote de matéria-prima suficiente para cinco formulações.
  • As formulações são preparadas por diversos operadores.
  • Duas fontes de perturbação: lotes e operadores.
  • Delineamento apropriado: testar cada formulação apenas uma vez em cada lote de matéria-prima e cada operador preparar uma vez cada formulação.

Tabela 2: Velocidade de queima de diferentes formulações (A, B, C, D, E) para diferentes lotes de matéria-prima e operadores.

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

# Lendo os dados
dados_3 <- read_xlsx("dados_3.xlsx")
dados_3
# A tibble: 25 × 4
    lote operador formulacoes velocidade
   <dbl>    <dbl> <chr>            <dbl>
 1     1        1 A                   24
 2     2        1 B                   17
 3     3        1 C                   18
 4     4        1 D                   26
 5     5        1 E                   22
 6     1        2 B                   20
 7     2        2 C                   24
 8     3        2 D                   38
 9     4        2 E                   31
10     5        2 A                   30
# ℹ 15 more rows
names(dados_3)
[1] "lote"        "operador"    "formulacoes" "velocidade" 
# Anova usando a função lm do pacote stats

# As variáveis com as linhas, colunas e tratamentos precisam ser fator
dados_3$lote <- as.factor(dados_3$lote)
dados_3$operador <- as.factor(dados_3$operador)
dados_3$formulacoes <- as.factor(dados_3$formulacoes)

# definindo o modelo linear
lm_dql <- lm(velocidade ~ lote + operador + formulacoes, data = dados_3)

# estimativas dos parâmetros do modelo linear
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.008855
# Análise de variância
anova(lm_dql)
Analysis 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 ' ' 1
# ANOVA usando ExpDes
# teste Tukey é o padrão para comparações múltiplas
latsd(treat = dados_3$formulacoes,
    row = dados_3$lote,
    column = dados_3$operador,
    resp = dados_3$velocidade)
------------------------------------------------------------------------
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 
------------------------------------------------------------------------

Análise de resíduos

Os pressupostos a serem verificados em uma análise de experimentos (anova) são: resíduos normais, independentes e de variância constante (homocedasticidade). Também é desejável verificar se há pontos discrepantes e a linearidade da variável resposta.

O resíduo bruto, ou simplesmente resíduo, é definido por:

\[ \hat{\epsilon}_{ij} = Y_{ij} - \hat{Y}_{ij}, \]

em que \(Y_{ij}\) é o valor observado e \(\hat{Y}_{ij}\) é o valor estimado pelo modelo.

# Resíduo bruto
res_1 <- residuals(lm_dic)
res_1
     1      2      3      4      5      6      7      8      9     10     11 
-0.425 -2.225  1.875  0.775  0.525  2.225 -1.975 -0.775 -0.250 -2.750  2.450 
    12     13     14     15     16     17     18     19     20     21     22 
 0.550  0.500  2.100 -2.000 -0.600 -0.050 -1.850  2.350 -0.450 -0.175 -1.975 
    23     24 
 1.925  0.225 
dados_1$altura - fitted(lm_dic)
     1      2      3      4      5      6      7      8      9     10     11 
-0.425 -2.225  1.875  0.775  0.525  2.225 -1.975 -0.775 -0.250 -2.750  2.450 
    12     13     14     15     16     17     18     19     20     21     22 
 0.550  0.500  2.100 -2.000 -0.600 -0.050 -1.850  2.350 -0.450 -0.175 -1.975 
    23     24 
 1.925  0.225

O resíduo padronizado (resíduo de Pearson) é usado para facilitar a comparação entre diferentes observações:

\[ d_{ij} = \frac{\hat{\epsilon}_{ij}}{\sqrt{\hat{\sigma}^2}}, \] em que \(\hat{\sigma}^2 = QM_{Erro}\) é o quadrado médio do erro da análise de variância.

# Resíduo padronizado (ou Pearson)
anova(lm_dic) # ANOVA
Analysis 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 ' ' 1
sigma2 <- anova(lm_dic)[2,3] # quadrado médio do erro
sigma2
[1] 3.264583
res_2 <- residuals(lm_dic) / sqrt(sigma2)
res_2
          1           2           3           4           5           6 
-0.23522044 -1.23144816  1.03773722  0.42893138  0.29056642  1.23144816 
          7           8           9          10          11          12 
-1.09308320 -0.42893138 -0.13836496 -1.52201459  1.35597663  0.30440292 
         13          14          15          16          17          18 
 0.27672992  1.16226568 -1.10691970 -0.33207591 -0.02767299 -1.02390072 
         19          20          21          22          23          24 
 1.30063065 -0.24905693 -0.09685547 -1.09308320  1.06541021  0.12452847

É possível padronizar os resíduos considerando a variabilidade individual dos resíduos (alavancagem). Ele é calculado como:

\[ r_{ij} = \frac{\hat{\epsilon}_{ij}}{\sqrt{\hat{\sigma}^2(1 - h_{ii})}}, \] em que \(h_{ii}\) é o \(i\)-ésimo elemento da diagonal da matriz hat \(H\)

Esse resíduo possibilita uma identificação mais precisa de pontos influentes e discrepantes. Esse resíduo também é comumente chamado de resíduo padronizado. É mais comum chamar esse resíduo de padronizado ao invés do de Pearson! Outro nome possível para esse resíduos é studentizado internamente.

# Resíduo padronizado (studentizado internamente)
res_3 <- rstandard(lm_dic)
res_3
          1           2           3           4           5           6 
-0.27160916 -1.42195386  1.19827572  0.49528730  0.33551720  1.42195386 
          7           8           9          10          11          12 
-1.26218376 -0.49528730 -0.15977010 -1.75747106  1.56574695  0.35149421 
         13          14          15          16          17          18 
 0.31954019  1.34206881 -1.27816077 -0.38344823 -0.03195402 -1.18229871 
         19          20          21          22          23          24 
 1.50183891 -0.28758617 -0.11183907 -1.26218376  1.23022974  0.14379309 
# comparação entre os três tipos de resíduos para LM
data.frame(res_1, res_2, res_3)
    res_1       res_2       res_3
1  -0.425 -0.23522044 -0.27160916
2  -2.225 -1.23144816 -1.42195386
3   1.875  1.03773722  1.19827572
4   0.775  0.42893138  0.49528730
5   0.525  0.29056642  0.33551720
6   2.225  1.23144816  1.42195386
7  -1.975 -1.09308320 -1.26218376
8  -0.775 -0.42893138 -0.49528730
9  -0.250 -0.13836496 -0.15977010
10 -2.750 -1.52201459 -1.75747106
11  2.450  1.35597663  1.56574695
12  0.550  0.30440292  0.35149421
13  0.500  0.27672992  0.31954019
14  2.100  1.16226568  1.34206881
15 -2.000 -1.10691970 -1.27816077
16 -0.600 -0.33207591 -0.38344823
17 -0.050 -0.02767299 -0.03195402
18 -1.850 -1.02390072 -1.18229871
19  2.350  1.30063065  1.50183891
20 -0.450 -0.24905693 -0.28758617
21 -0.175 -0.09685547 -0.11183907
22 -1.975 -1.09308320 -1.26218376
23  1.925  1.06541021  1.23022974
24  0.225  0.12452847  0.14379309

Normalidade dos resíduos

  • Histograma dos resíduos padronizados (\(d_{ij}\) ou \(r_{ij}\)): Deve apresentar uma forma semelhante à curva normal padrão.
# histograma dos resíduos
hist(res_3)

  • Gráfico dos valores estimados versus resíduos padronizados: Deve apresentar pontos dispostos de forma aleatória, geralmente concentrados no intervalo \((-2, 2)\).
# Gráfico dos valores estimados versus resíduos padronizados
y_est <- fitted(lm_dic)
plot(y_est, res_3) # os pontos estão entre (-2, 2)?

  • Comparação de frequências: Consiste em comparar a frequência relativa dos resíduos padronizados nos intervalos \((-1, 1)\), \((-1,64, 1,64)\) e \((-1,96, 1,96)\) com as probabilidades teóricas de 68%, 90% e 95%, respectivamente, da distribuição normal padrão.
# comparação de frequências
intervalos <- c("-1 a 1", "-1.64 a 1.64", "-1.96 a 1.96")
frequencias_r <- c(mean(res_3 > -1 & res_3 < 1),
                 mean(res_3 > -1.64 & res_3 < 1.64),
                 mean(res_3 > -1.96 & res_3 < 1.96))
probabilidades <- c(0.68, 0.90, 0.95)
data.frame(intervalos, frequencias_r, probabilidades)
    intervalos frequencias_r probabilidades
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.95
  • Gráfico quantil-quantil (Q-Q plot): Com os resíduos ordenados versus os quantis da distribuição normal padronizada, o gráfico deve aproximar-se da reta bissetriz, demonstrando a aderência à normalidade.
# qq-plot
qqnorm(res_3)
qqline(res_3)

Além dessas abordagens visuais, a aplicação de testes de aderência (não-paramétricos) também pode ser realizada. Dentre esses, os mais conhecidos incluem:
Qui-quadrado, Kolmogorov-Smirnov, Jarque-Bera, Shapiro-Wilk, Anderson-Darling, Cramér-von Mises, D’Agostino-Pearson, Lilliefors, e Shapiro-Francia.

# Teste de Shapiro-Wilk
shapiro.test(res_3)

    Shapiro-Wilk normality test

data:  res_3
W = 0.9384, p-value = 0.1502
# Teste de Kolmogorov-Smirnov
ks.test(res_3, "pnorm", mean = 0, sd = 1)

    Exact one-sample Kolmogorov-Smirnov test

data:  res_3
D = 0.1346, p-value = 0.728
alternative hypothesis: two-sided

Independência dos Resíduos

Sem a informação sobre a ordem de execução dos experimentos, não é possível aplicar testes estatísticos para a independência que dependem de uma sequência temporal ou da coleta dos dados. Nesse caso, a inferência sobre a independência dos resíduos deve basear-se na robustez do desenho experimental - que, se bem realizado com aleatorização, repetição e controle local, sugere naturalmente a independência dos resíduos.

Alguns métodos para avaliar a independência incluem:

  • Gráfico dos valores estimados versus resíduos padronizados: Os pontos devem estar distribuídos aleatoriamente, sem a formação de padrões específicos.
# Gráfico dos valores estimados versus resíduos padronizados
plot(y_est, res_3)

  • Gráfico da ordem de coleta dos dados versus resíduos padronizados: Quando a ordem de coleta é conhecida, a ausência de padrões neste gráfico indica independência.

  • Função de autocorrelação: Permite avaliar a correlação entre os resíduos em diferentes defasagens.

  • Testes estatísticos: Entre os testes de independência, destaca-se o teste Durbin-Watson, além do teste de Breusch-Godfrey e Ljung-Box.

Homocedasticidade dos Resíduos

A homocedasticidade refere-se à igualdade das variâncias dos resíduos para todos os tratamentos. Essa verificação pode ser feita por:

  • Gráfico dos valores estimados versus resíduos padronizados: A ausência de padrão ou “funil” indica homogeneidade das variâncias.
# Gráfico dos valores estimados versus resíduos padronizados
plot(y_est, res_3)

  • Boxplot ou gráfico de dispersão dos valores observados \((Y_{ij}\) versus tratamentos: Visualmente, as dispersões entre os diferentes tratamentos devem ser semelhantes.
# Boxplot dos valores observados em função dos tratamentos
boxplot(dados_1$altura ~ dados_1$substratos)

  • Testes de homogeneidade de variâncias: Alguns testes utilizados incluem Bartlett, Breusch-Pagan, Levene, Samiuddin, O’Neill e Mathews, Layard, Park, White, Cochran, Hartley e Goldfeld-Quandt.
# Teste de Bartlett
bartlett.test(res_3 ~ dados_1$substratos)

    Bartlett test of homogeneity of variances

data:  res_3 by dados_1$substratos
Bartlett's K-squared = 0.27704, df = 5, p-value = 0.9981
# Teste de Breusch-Pagan
# library(lmtest)
bptest(lm_dic)

    studentized Breusch-Pagan test

data:  lm_dic
BP = 1.0259, df = 5, p-value = 0.9605

Pontos discrepantes (outliers)

Embora não seja um pressuposto obrigatório, a detecção de outliers é fundamental para identificar erros de coleta ou registros anômalos. Métodos simples para essa verificação são:

  • Boxplot dos resíduos padronizados para cada tratamento: Facilita a visualização de observações fora do padrão.
# Boxplot dos resíduos padronizados
boxplot(res_3 ~ dados_1$substratos)

  • Gráfico dos valores estimados versus resíduos padronizados: Observações que se afastam do intervalo \((-2, 2)\) devem ser analisadas como possíveis outliers.

Linearidade da Variável Resposta

A linearidade entre a variável resposta e os tratamentos é um pressuposto básico para muitos modelos estatísticos. Uma verificação simples pode ser realizada através de:

  • Gráfico de dispersão (ou boxplot) dos valores observados em função dos tratamentos: Se a variabilidade não for comum a todos os tratamentos, pode ser necessário transformar os dados para um melhor ajuste do modelo.
# Boxplot dos valores observados em função dos tratamentos
boxplot(dados_1$altura ~ dados_1$substratos)

Transformação de Variáveis

Quando os pressupostos do modelo não são atendidos, transformar os dados pode ajudar a melhorar a análise dos resíduos. A seguir, são apresentadas algumas transformações comuns:

Transformações Específicas

  • Dados de contagem (distribuição de Poisson):
    \(y = \sqrt{x}\) ou \(y = \sqrt{x} + 0,5\)

  • Dados com distribuição assimétrica (log-normal):
    \(y = \log(x)\) ou \(y = \log(x + c)\), onde \(0 < c \leq 1\) (comumente \(c = 1\) ou \(1/2\).

  • Dados de proporção (binomial):
    \(y = \arcsin(\sqrt{x})\), também conhecido como transformação arco-seno.

Transformações por Potência

Uma forma geral de transformação é a potência, definida como \(y = x^p\):

  • Quando \(p = 0\), utilizamos a transformação logarítmica: \(y = \log(x)\).
  • Quando \(p = 2\), obtemos: \(y = x^2\).
  • Quando \(p = \frac{1}{2}\), temos: \(y = \sqrt{x}\).
  • Quando \(p = -1\), a transformação é: \(y = \frac{1}{x}\).
  • Quando \(p = -\frac{1}{2}\), resulta em: \(y = \frac{1}{\sqrt{x}}\).

Transformação de Box-Cox

A transformação de Box-Cox é uma técnica poderosa para estabilizar a variância e aproximar a normalidade dos resíduos. Ela é definida por:

\[ y = \frac{x^p - 1}{p}, \quad \text{para } p \neq 0, \]

\[ y = \log(x), \quad \text{para } p = 0. \]

O parâmetro \(p\) pode ser estimado através da maximização da função:

\[ L(p) = -\frac{1}{2}\log\left[ QMErro(p) \right], \]

em que \(QMErro(p)\) representa o quadrado médio do erro obtido na análise de variância após a transformação de Box-Cox.

# transformação boxcox
# library(MASS)
boxcox(lm_dic, lambda = seq(-2, 2, 0.1))

Experimento

Arquivo com os dados: dados_4.xlsx

# Lendo os dados
# dados simulados
dados_4 <- read_xlsx("dados_4.xlsx")
dados_4
# A tibble: 24 × 3
   repeticao substratos altura
       <dbl> <chr>       <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.5
 7         3 CaC+EA       10.3
 8         4 CaC+EA       11.5
 9         1 HM           11.8
10         2 HM           19.3
# ℹ 14 more rows
names(dados_4)
[1] "repeticao"  "substratos" "altura"    
# Análise gráfica da variável resposta por tratamento
boxplot(dados_4$altura ~ dados_4$substratos)

# Ajuste do modelo
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 ' ' 1
# Resíduo padronizado
res_4 <- rstandard(lm_dic2)

# histograma dos resíduos
hist(res_4)

# Gráfico dos valores estimados versus resíduos padronizados
y_est <- fitted(lm_dic2)
plot(y_est, res_4) # os pontos estão entre (-2, 2)?

# qq-plot
qqnorm(res_4)
qqline(res_4)

# box plot dos resíduos por tratamento
boxplot(res_4 ~ dados_4$substratos)

# Teste de Shapiro-Wilk
shapiro.test(res_4)

    Shapiro-Wilk normality test

data:  res_4
W = 0.88827, p-value = 0.01224
# Teste de Kolmogorov-Smirnov
ks.test(res_4, "pnorm", mean = 0, sd = 1)

    Exact one-sample Kolmogorov-Smirnov test

data:  res_4
D = 0.16165, p-value = 0.5063
alternative hypothesis: two-sided
# teste de homogeneidade de variâncias
bartlett.test(res_4 ~ dados_4$substratos)

    Bartlett test of homogeneity of variances

data:  res_4 by dados_4$substratos
Bartlett's K-squared = 14.889, df = 5, p-value = 0.01085
# transformação boxcox
boxcox(lm_dic2, lambda = seq(-3, 3, 0.1))

dados_4$altura_2 <- log(dados_4$altura)
boxplot(dados_4$altura_2 ~ dados_4$substratos)

lm_dic3 <- lm(altura_2 ~ substratos, data = dados_4)
anova(lm_dic3)
Analysis of Variance Table

Response: altura_2
           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 ' ' 1
res_5 <- rstandard(lm_dic3)
bartlett.test(res_5 ~ dados_4$substratos)

    Bartlett test of homogeneity of variances

data:  res_5 by dados_4$substratos
Bartlett's K-squared = 9.2157, df = 5, p-value = 0.1008
shapiro.test(res_5)

    Shapiro-Wilk normality test

data:  res_5
W = 0.94999, p-value = 0.2708
# Anova e Tukey
crd(treat = dados_4$substratos,
    resp = dados_4$altura_2,
    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 
------------------------------------------------------------------------
# Scott-Knott
crd(treat = dados_4$substratos,
    resp = dados_4$altura_2,
    mcomp = "sk")
------------------------------------------------------------------------
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.
------------------------------------------------------------------------

Scott-Knott test
------------------------------------------------------------------------
  Groups Treatments    Means
1      a         HM 3.027798
2      a     CaC+EB 2.921100
3      b   Plantmax 2.749833
4      b     Cap+EA 2.654726
5      b     Cap+EB 2.562232
6      b     CaC+EA 2.499521
------------------------------------------------------------------------