k=Númerodeníveisdofator= laboratorios = 4 n=NúmerodeRepetiçõesdecada nível = 8

*= Itálico Ex.: Isaquinho Isaquinho **= Negrito ***=Itálico + Negrito Quantidade de “#” = tamanho da fonte (quanto mais menor a fonte)

1. Análise de Poder e Número de Repetições

1.1. Determinação de N para diferença padronizada pequena (f=0.1)

pwr.anova.test(k = 4, f = 0.1, sig.level = .05, power = .8)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 4
##               n = 273.5429
##               f = 0.1
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

Isso foi pra achar o n (número de repetições)! Sempre fica de fora d função o que se quer achar.

1.2. Determinaçao de N - diferença padronizada média

pwr.anova.test(k = 4, f = 0.25, sig.level = .05, power = .8)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 4
##               n = 44.59927
##               f = 0.25
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

1.3. Determinaçao de n - diferença padronizada grande

pwr.anova.test(k = 4, f = 0.4, sig.level = .05, power = .8)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 4
##               n = 18.04262
##               f = 0.4
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

1.4. Estimativa de Poder do Teste - Diferenca Padronizada Pequena

pwr.anova.test(k = 4, n=8, f = 0.1, sig.level = .05)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 4
##               n = 8
##               f = 0.1
##       sig.level = 0.05
##           power = 0.06694612
## 
## NOTE: n is number in each group

1.5. Estimativa de Poder do Teste - Diferenca Padronizada Média

pwr.anova.test(k = 4, n=8, f = 0.25, sig.level = .05)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 4
##               n = 8
##               f = 0.25
##       sig.level = 0.05
##           power = 0.1720053
## 
## NOTE: n is number in each group

1.6. Estimativa de Poder do Teste - Diferenca Padronizada Grande

pwr.anova.test(k = 4, n=8, f = 0.4, sig.level = .05)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 4
##               n = 8
##               f = 0.4
##       sig.level = 0.05
##           power = 0.3967438
## 
## NOTE: n is number in each group

2. Análise Exploratória dos Dados

2.1. Importando e verificando os dados

labs <- read.table("labs.txt", header = TRUE)
str(labs)
## 'data.frame':    32 obs. of  2 variables:
##  $ Lab  : Factor w/ 4 levels "Lab1","Lab2",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ cromo: num  26.1 21.5 22 22.6 24.9 22.6 23.8 23.2 18.3 19.7 ...
# Header = TRUE significa que as minhas variaveis tem o nome no cabeçalho

2.2. Boxplot - R Base

boxplot(cromo ~ Lab, data = labs, col = "lightgray")

2.3. Boxplot - Usando o pacote “ggplot2”

library(ggplot2)
ggplot(labs, aes(x=Lab, y=cromo)) + geom_boxplot()

2.4. Estatísticas Descritivas

library(summarytools)
## Registered S3 method overwritten by 'pryr':
##   method      from
##   print.bytes Rcpp
descr(labs$cromo, style = "rmarkdown")
## ### Descriptive Statistics  
## #### labs$cromo  
## **N:** 32  
## 
## |          &nbsp; |  cromo |
## |----------------:|-------:|
## |        **Mean** |  21.30 |
## |     **Std.Dev** |   5.04 |
## |         **Min** |  11.00 |
## |          **Q1** |  18.15 |
## |      **Median** |  21.20 |
## |          **Q3** |  25.20 |
## |         **Max** |  30.70 |
## |         **MAD** |   5.11 |
## |         **IQR** |   6.82 |
## |          **CV** |   0.24 |
## |    **Skewness** |  -0.04 |
## | **SE.Skewness** |   0.41 |
## |    **Kurtosis** |  -0.64 |
## |     **N.Valid** |  32.00 |
## |   **Pct.Valid** | 100.00 |

Std. Dev = Desvio Padrão Q3 - Quartil 3 MAD = Desvio Absoluto Médio IQR = Intervalo interquartil Skewness = simetria SE. Skewness = desvio da simetria CV = Coeficiente de variação - este serve para comparar distribuição de dados para amostras simétricas (duas amostras) que é quando a mediana está no meio do boxplot, o histograma está em forma de sino. CV = desvio padrão divido pela média x 100 = distribuição simetrica CV = mediana dividida “MAD” x 100

2.5. Estimativa do Modelo

mod1<- aov(cromo ~ Lab, data = labs)
summary(mod1)
##             Df Sum Sq Mean Sq F value     Pr(>F)    
## Lab          3  476.1  158.69   14.21 0.00000814 ***
## Residuals   28  312.7   11.17                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Exibe a tabela de ANAVA
summary(mod1)
##             Df Sum Sq Mean Sq F value     Pr(>F)    
## Lab          3  476.1  158.69   14.21 0.00000814 ***
## Residuals   28  312.7   11.17                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimativa dos Efeitos Principais \(\hat{\tau}_i\)

coef(mod1)
## (Intercept)     LabLab2     LabLab3     LabLab4 
##     23.3375     -6.5500     -4.7875      3.2000

3. Diagnóstico do Modelo

3.1. Cv do Experimento

library(agricolae)
cv.model(mod1)
## [1] 15.68622

3.2. Análise dos Resíduos

plot(mod1)

3.3. Teste de Shapiro-Wilks

shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.96141, p-value = 0.3001

Assim, o valor de p é menor que o nível de significância, aceito a hipótese nula de que o resíduo atende uma distribuição normal

3.4. Análise dos Resíduos

names(mod1)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"
# Resíduos
mod1$residuals
##       1       2       3       4       5       6       7       8       9      10 
##  2.7625 -1.8375 -1.3375 -0.7375  1.5625 -0.7375  0.4625 -0.1375  1.5125  2.9125 
##      11      12      13      14      15      16      17      18      19      20 
##  1.2125  0.6125  5.8125 -5.1875 -5.7875 -1.0875  0.5500 -4.6500 -2.8500  0.0500 
##      21      22      23      24      25      26      27      28      29      30 
##  0.5500 -1.7500  6.9500  1.1500  4.1625  0.7625 -5.6375  2.4625 -5.6375 -0.4375 
##      31      32 
##  0.1625  4.1625
# Dados da Amostra
labs$cromo
##  [1] 26.1 21.5 22.0 22.6 24.9 22.6 23.8 23.2 18.3 19.7 18.0 17.4 22.6 11.6 11.0
## [16] 15.7 19.1 13.9 15.7 18.6 19.1 16.8 25.5 19.7 30.7 27.3 20.9 29.0 20.9 26.1
## [31] 26.7 30.7
# Valores previstos pelo modelo
mod1$fitted.values
##       1       2       3       4       5       6       7       8       9      10 
## 23.3375 23.3375 23.3375 23.3375 23.3375 23.3375 23.3375 23.3375 16.7875 16.7875 
##      11      12      13      14      15      16      17      18      19      20 
## 16.7875 16.7875 16.7875 16.7875 16.7875 16.7875 18.5500 18.5500 18.5500 18.5500 
##      21      22      23      24      25      26      27      28      29      30 
## 18.5500 18.5500 18.5500 18.5500 26.5375 26.5375 26.5375 26.5375 26.5375 26.5375 
##      31      32 
## 26.5375 26.5375

3.5. Teste de Bartlett

bartlett.test(cromo ~ Lab, data = labs)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  cromo by Lab
## Bartlett's K-squared = 5.7637, df = 3, p-value = 0.1237

Como no teste de Barlett o valor de p está acima do nível de significância, aceito a hipótese nula de que os resíduos tem variância constante.

Como a ANAVA indicou diferenças entre as medições dos laboratórios, surgem outras perguntas. • Quais laboratórios apresentaram medições estatisticamente iguais e quais apresentaram medições estatisticamente iguais? • Necessitamos de procedimentos de comparações múltiplas de médias para responder estas questões, veremos dois; • Procedimentos de Tukey e de Fisher

3.6. Testes de Comparações Múltiplas

3.6.1. Teste de Fisher

LSD.test(mod1,"Lab", p.adj = "bon", console = TRUE)
## 
## Study: mod1 ~ "Lab"
## 
## LSD t Test for cromo 
## P value adjustment method: bonferroni 
## 
## Mean Square Error:  11.16665 
## 
## Lab,  means and individual ( 95 %) CI
## 
##        cromo      std r     LCL     UCL  Min  Max
## Lab1 23.3375 1.538030 8 20.9174 25.7576 21.5 26.1
## Lab2 16.7875 3.927717 8 14.3674 19.2076 11.0 22.6
## Lab3 18.5500 3.444250 8 16.1299 20.9701 13.9 25.5
## Lab4 26.5375 3.874435 8 24.1174 28.9576 20.9 30.7
## 
## Alpha: 0.05 ; DF Error: 28
## Critical Value of t: 2.838933 
## 
## Minimum Significant Difference: 4.743366 
## 
## Treatments with the same letter are not significantly different.
## 
##        cromo groups
## Lab4 26.5375      a
## Lab1 23.3375      a
## Lab3 18.5500      b
## Lab2 16.7875      b

3.6.2. Teste de Tukey

# Função Interna do R
TukeyHSD(mod1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = cromo ~ Lab, data = labs)
## 
## $Lab
##              diff        lwr        upr     p adj
## Lab2-Lab1 -6.5500 -11.111878 -1.9881216 0.0027656
## Lab3-Lab1 -4.7875  -9.349378 -0.2256216 0.0369829
## Lab4-Lab1  3.2000  -1.361878  7.7618784 0.2447531
## Lab3-Lab2  1.7625  -2.799378  6.3243784 0.7190930
## Lab4-Lab2  9.7500   5.188122 14.3118784 0.0000163
## Lab4-Lab3  7.9875   3.425622 12.5493784 0.0002808
# Exibição Via Gráfico
plot(TukeyHSD(mod1))

4. Tutorial

4.1. Importando os Dados

turtles <- read.csv(file = "turtles.csv", header = TRUE) 
str(turtles) #verifica a estrutura dos dados importados
## 'data.frame':    40 obs. of  2 variables:
##  $ temperatura: int  15 15 15 15 15 15 15 15 15 15 ...
##  $ dias       : int  37 43 45 54 56 65 62 73 74 75 ...
head(turtles) # exibe as primeiras linhas dos dados importados
##   temperatura dias
## 1          15   37
## 2          15   43
## 3          15   45
## 4          15   54
## 5          15   56
## 6          15   65
tail(turtles) # exibe as últimas linhas dos dados importados
##    temperatura dias
## 35          30   12
## 36          30   18
## 37          30   21
## 38          30   23
## 39          30   29
## 40          30   39
turtles # exibe todos os dados importados
##    temperatura dias
## 1           15   37
## 2           15   43
## 3           15   45
## 4           15   54
## 5           15   56
## 6           15   65
## 7           15   62
## 8           15   73
## 9           15   74
## 10          15   75
## 11          20   30
## 12          20   31
## 13          20   34
## 14          20   35
## 15          20   35
## 16          20   47
## 17          20   53
## 18          20   54
## 19          20   63
## 20          20   64
## 21          25   21
## 22          25   23
## 23          25   48
## 24          25   52
## 25          25   52
## 26          25   54
## 27          25   54
## 28          25   61
## 29          25   62
## 30          25   65
## 31          30   13
## 32          30   16
## 33          30   19
## 34          30   11
## 35          30   12
## 36          30   18
## 37          30   21
## 38          30   23
## 39          30   29
## 40          30   39
mean(turtles$dias)
## [1] 43.075
mean(turtles$temperatura) # Para acessar os frames usa-se o $
## [1] 22.5

4.2. Fator Variável de Temperatura

turtles$temperatura <- factor(turtles$temperatura) 
str(turtles)
## 'data.frame':    40 obs. of  2 variables:
##  $ temperatura: Factor w/ 4 levels "15","20","25",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ dias       : int  37 43 45 54 56 65 62 73 74 75 ...

4.3. Gráfico Boxplot Turtles

boxplot(dias ~ temperatura, data = turtles, col = "indianred")

### 4.4. Estimação do Modelo Linear Anava (Turtles)

turtles.aov <- aov(dias ~ temperatura, data = turtles) 
summary(turtles.aov)
##             Df Sum Sq Mean Sq F value      Pr(>F)    
## temperatura  3   8025  2675.2   15.98 0.000000908 ***
## Residuals   36   6027   167.4                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.5. Interpretação dos Efeitos Principais

coef(turtles.aov)
##   (Intercept) temperatura20 temperatura25 temperatura30 
##          58.4         -13.8          -9.2         -38.3

4.6. Teste de Tukey - Comparações Multiplas

# Teste Tukey Pacote Agricolae
teste_tukey <- agricolae::HSD.test (turtles.aov,"temperatura", group = TRUE, console = TRUE)
## 
## Study: turtles.aov ~ "temperatura"
## 
## HSD Test for dias 
## 
## Mean Square Error:  167.425 
## 
## temperatura,  means
## 
##    dias       std  r Min Max
## 15 58.4 13.696715 10  37  75
## 20 44.6 13.226237 10  30  64
## 25 49.2 15.266521 10  21  65
## 30 20.1  8.608136 10  11  39
## 
## Alpha: 0.05 ; DF Error: 36 
## Critical Value of Studentized Range: 3.808798 
## 
## Minimun Significant Difference: 15.58469 
## 
## Treatments with the same letter are not significantly different.
## 
##    dias groups
## 15 58.4      a
## 25 49.2      a
## 20 44.6      a
## 30 20.1      b
teste_tukey
## $statistics
##   MSerror Df   Mean       CV      MSD
##   167.425 36 43.075 30.03896 15.58469
## 
## $parameters
##    test      name.t ntr StudentizedRange alpha
##   Tukey temperatura   4         3.808798  0.05
## 
## $means
##    dias       std  r Min Max   Q25  Q50   Q75
## 15 58.4 13.696715 10  37  75 47.25 59.0 71.00
## 20 44.6 13.226237 10  30  64 34.25 41.0 53.75
## 25 49.2 15.266521 10  21  65 49.00 53.0 59.25
## 30 20.1  8.608136 10  11  39 13.75 18.5 22.50
## 
## $comparison
## NULL
## 
## $groups
##    dias groups
## 15 58.4      a
## 25 49.2      a
## 20 44.6      a
## 30 20.1      b
## 
## attr(,"class")
## [1] "group"
teste2 <- agricolae::HSD.test(turtles.aov,"temperatura",
                              group=FALSE)
                              print(teste2$comparison)
##         difference pvalue signif.        LCL      UCL
## 15 - 20       13.8 0.0983       .  -1.784689 29.38469
## 15 - 25        9.2 0.3970          -6.384689 24.78469
## 15 - 30       38.3 0.0000     ***  22.715311 53.88469
## 20 - 25       -4.6 0.8563         -20.184689 10.98469
## 20 - 30       24.5 0.0008     ***   8.915311 40.08469
## 25 - 30       29.1 0.0001     ***  13.515311 44.68469

4.7. Gráfico do Teste de Tukey

plot(teste_tukey, main = "Teste de Tukey")

### 4.8. Diagnóstico: Análise de Resíduos #### 4.8.1. Gráfico dos Resíduos x Valores Previstos do Modelo

plot(turtles.aov, 1)

4.8.2. Gráfico Quantil-Quantil Normal

plot(turtles.aov, 2)

#### 4.8.3. Teste de Shapiro-Wilks

shapiro.test(turtles.aov$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  turtles.aov$residuals
## W = 0.97008, p-value = 0.3619

Não rejeitamos a hipótese nula, visto que nível de significância do teste de Shapiro-Wilks apresentou valor p acima e o rsíduos possui distribuição normal com média = 0.

4.8.4. Gráfico de Variância Homôgenea

plot(turtles.aov, 3)

Não rejeitamos a hipótese nula, visto que a variância dos resíduos é aproximadamente constante.

4.9. Teste de Bartlett - Variância Homogenea

bartlett.test(dias ~ temperatura, data = turtles)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  dias by temperatura
## Bartlett's K-squared = 2.8101, df = 3, p-value = 0.4218

5. Complemento

5.1. Alterando o Nível de Referência

turtles2 <-within(turtles, temperatura <- relevel (temperatura, ref = 4))
turtle.ref <- aov(dias ~ temperatura, data = turtles2)
coef(turtle.ref)
##   (Intercept) temperatura15 temperatura20 temperatura25 
##          20.1          38.3          24.5          29.1

5.2. Teste de Kruskal-Wallis - ANOVA - Alternativa Não-paramétrica

kruskal.test(dias ~ temperatura, data = turtles)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  dias by temperatura
## Kruskal-Wallis chi-squared = 21.491, df = 3, p-value = 0.00008325