Aula22

Notas

Para add pacotes sempre tem que digitar: library(nome do pacote) e depois rodar.

Pacotes add: library(ggplot2), library(summarytools), library(agricolae).

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

1.1 Determinação de N (Diferença pequena)

pwr.anova.test(k=4, f=0.1, sig.level = 0.05, power = 0.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

1.2 Determinação de N (Diferença média)

pwr.anova.test(k=4, f=0.25, sig.level = 0.05, power = 0.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ção de N (Diferença grande)

pwr.anova.test(k=4, f=0.4, sig.level = 0.05, power = 0.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 do Poder do Teste (Diferença pequena)

pwr.anova.test(k=4, n=8, f=0.1, sig.level = 0.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 do Poder do Teste (Diferença média)

pwr.anova.test(k=4, n=8, f=0.25, sig.level = 0.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 do Poder do Teste (Diferença grande)

pwr.anova.test(k=4, n=8, f=0.4, sig.level = 0.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 Importanto 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 ...

2.2 Gráfico 1

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

2.3 Gráfico 2

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

3. Estimação do Modelo

# aov: estima o modelo linear
mod1<- aov(cromo ~ Lab, data = labs)

# Exibe a tabela da 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

3.1 Estimativas dos efeitos principais \(\hat{\tau}_i\)

coef(mod1)

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

4. Diagnóstico do modelo

4.1 CV do experimento

library(agricolae)

cv.model(mod1)

## [1] 15.68622

4.2 Análises dos Resíduos

erro = erro da pop \(\epsilon_{ij}\) residuo =\(\hat{\epsilon}_{ij}\) = estimativa do erro da população a partir da amostra

# O que está contido no objeto mod1 ?
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

plot(mod1)

4.2.1 Teste de Shapiro-Wilks - Normalidade dos resíduos

shapiro.test(mod1$residuals)

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

4.2.2 Teste de Bartlett: variâncias homogêneas

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

4.3 Testes de Comparações Multiplas

4.3.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

4.3.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))

5. Tutorial - Executando a análise

5.1 Importanto e verificando os dados

turtles <- read.csv(file = "turtles.csv", header = TRUE) # importa o arquivo
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 # xibe 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

Convertendo a variável ‘temperatura’ para a classe ou tipo ‘factor’

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 ...

5.2 Análise Exploratória dos dados

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

Obs.: Para acessar a paleta de cores do R base, digite ‘colours()’ no console e clica enter.

5.3 Estimação do Modelo linear da ANAVA

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

5.3.1 Interpretação dos Efeitos Principais

coef(turtles.aov)

##   (Intercept) temperatura20 temperatura25 temperatura30 
##          58.4         -13.8          -9.2         -38.3

5.4 Teste de Tukey - Comparações Multiplas

### Teste Tukey 1 do Pacote Agricolae 

teste_tukey <- agricolae::HSD.test (turtles.aov,"temperatura", group = TRUE, console = FALSE)
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"

### Teste Tukey 2 do Pacote Agricolae 

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

### Teste de Tukey (grafico) do Pacote Agricolae)
plot(teste_tukey, main = "Teste de Tukey")

5.5 Diagnóstico: Análise dos Resíduos

5.5.1 Resíduos tem média zero

### Gráfico dos resíduos x valores previstos pelo modelo
plot(turtles.aov, 1)

5.5.2 Resíduos tem distribuição aproximadamente normal

### Gráfico quantil - quantil normal
plot(turtles.aov, 2)

### teste formal
shapiro.test(turtles.aov$residuals)

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

5.5.3 Resíduos tem variância homogênea

### Gráfico de escala-locação
plot(turtles.aov, 3)

### Teste formal
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

6. Complemento

6.1 Alterando o Nível de Referência

Comentário: (Ao invés de deixar o R buscar o primeiro nível automaticamente, eu posso informar o meu nível de referencia por.ex. nível 4 para ter como base.)

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

Teste Kruskal-Wallis - ANOVA não parametrica (Quando as hipóteses não são verificadas?)

### Teste formal Kruskal
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

Obs.: Parametrio: assumo a posição dos meus dados.Preciso estimar os parametros de “media” e “coeficiente” para determinar a normal por exemplo. Nao parametrico não assumo a posição dos meus dados. É quando é utilizado as localizações de onde os valores estão na tabela. ou seja, nao é considerado os valores de fato e sim suas posições.

Conclusão: o resultado de p-value do teste kruskal está muito abaixo de 0.05, portanto há pelo menos uma diferença entre os reusltados das medias das temperaturas. Logo, rejeita-se a hipotese normal (H0)