#1.Analise do poder enumero de Repetições ##1.1 Determinação de N-diferença padronizada pequena

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

#1.2 Determinação 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ção 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 do poder do teste- diferença 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.6 Estimativa do poder do teste- diferença padronizada media

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.7 Estimativa do poder do teste- diferença 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

Interpretação: diferencia padronizada pequena: tenho 7 % de rejeitar corretamente Ho diferencia padronizada média: 17% de rejeitar corretamente Ho diferencia padronizada média:40% de rejeitar corretamente Ho

#2.Analise exploratoria dos dados #2.5 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 ...

#2.6 Boxplot R-base

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

#2.7 Boxplot- usando pacote ‘ggplot2’

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

Interpretação: ploto minhas 32 observação de cromo contra os laboratorios. Quanto maior a caixa do boxplot, maior a variabilidade do dados. Para a distribuição ser simetrica, a mediana tem que estar no meio da caixa. Lab 2 e 4 apresentaram maior variabilidade. La 3 apresentou um valor outlander. Lab 1 apresentou menor variação,devido o tamanho da caixa.

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

OBS: quando adistribuição é simétrica: média e o desvio padrão são bons dados. Quando é assimetric: a medianda e MAD (descvio absoluto médio) é melhor opção. Variabilidade entre duas amostras para a distribuição simétrica: usar o coeficiente de variação (cv=devio padrão/média), expresso em porcentagem e para distribuição assimetrica: (CV=mediana/MAD). Skeamer= coeficiente assimetria, SE skeammer= erro padrão do coeficiente de assimetria, kurtosis=peso da cauda de distribuição.

#3. Estimação do modelo

mod1 <- aov(cromo ~ Lab, data = labs)

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

Interpretação: a esquerda do ~ é a variavel reposta e a direita é o fator de interesse.Grau de liberdade=df,quadrado médio (mean sq), soma de quadrado (sum sq),pr(>F)=pvalor

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

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

Interpretação: São sempre relativas ao fator de referencia, fator base= interceot. Primeiro faor como referencia lab1 contidas a ele. Os valores de lab2,lab3,lab4 é calculado da seuinte maneira lab1=lab2,…)

4. Diagnostico do modelo

##4.1 CV do experimento

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

##4.2 Analise dos residuos

erro do modelo antes da estimação:erro -\(\epsilon_{i}\) - erro da população

residuo-\(\hat{epsilon}_{i}\) - estimativa do erro da população

names(mod1)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"
#residuos
mod1$residuals
##       1       2       3       4       5       6       7       8       9 
##  2.7625 -1.8375 -1.3375 -0.7375  1.5625 -0.7375  0.4625 -0.1375  1.5125 
##      10      11      12      13      14      15      16      17      18 
##  2.9125  1.2125  0.6125  5.8125 -5.1875 -5.7875 -1.0875  0.5500 -4.6500 
##      19      20      21      22      23      24      25      26      27 
## -2.8500  0.0500  0.5500 -1.7500  6.9500  1.1500  4.1625  0.7625 -5.6375 
##      28      29      30      31      32 
##  2.4625 -5.6375 -0.4375  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
## [15] 11.0 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
## [29] 20.9 26.1 26.7 30.7
#valores previstos pelo modelo
mod1$fitted.values
##       1       2       3       4       5       6       7       8       9 
## 23.3375 23.3375 23.3375 23.3375 23.3375 23.3375 23.3375 23.3375 16.7875 
##      10      11      12      13      14      15      16      17      18 
## 16.7875 16.7875 16.7875 16.7875 16.7875 16.7875 16.7875 18.5500 18.5500 
##      19      20      21      22      23      24      25      26      27 
## 18.5500 18.5500 18.5500 18.5500 18.5500 18.5500 26.5375 26.5375 26.5375 
##      28      29      30      31      32 
## 26.5375 26.5375 26.5375 26.5375 26.5375

#4.3 Grafico dos residuos versus valores pelo modelo

plot(mod1,1)

#5.Teste Shapiro Wilks

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

#6.Teste 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

#7. Teste de comparação multiplas

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

#7.2 Teste Tukey

#funçao 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 grafico
plot(TukeyHSD(mod1))

#5. TUTORIAL

##5.1 Importanto e verificando os dados

turtles <- read.csv(file = "turtles.csv", header = TRUE) #importa os dados
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 ultimas 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

#Convertendo a variavel ‘temperatura’ para classe/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 Analise Exploratoria dos dados

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

OBS: Para acessar a paleta de cores do R base. digite ‘colours’ no console

Interpretação do grafico: as medias da temperatura de 15,20,25 nao vao dar estaticamente diferente -observando a posiçao das caixas,c mas a de 30 é diferente. A maior variabilidade de dados ocorreu para a tempertura de 15 graus.

#AULA 23/11

#Estimação do modelo

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

Interpretação: rejeitamos o Ho. Neste exemplo, o valor-p do teste da temperatura é 0.000090. Ou seja é muito improvavel de ter observado o valor F (ou os dados) se hipotese nula de que as quatro temperatura médias são iguais se fosse verdadeira

#Estimativas dos efeitos principais

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

Interpretação: são sempre relativa ao fator de referencia fator base= intercept (nivel de referencia), tomou como fator de referenciqa a temp 15, os valores de temp 20,25,30 é calculado da seguinte maneira temp20-temp15…

#Comparação multiplas

#Teste de Tukey 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 de Tukey de 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 de pacote Agricolae
plot(teste_tukey, main = "Teste de Tukey")

#Diagnostico: analise

#grafico dos residuos x valores previstos
plot(turtles.aov,1)

#Grafico quantil-quantil normal

#Grafico quantil-quantil normal
plot(turtles.aov, 2)

#Teste de Shapiro

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

Residuos tem variancia homogenica?

Grafico dos residuos padronizado

plot(turtles.aov, 3)

#Teste de Bartlet

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

#Complementos

#Alterando o nivel referencia

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 de Kruskal-wallis-ANOVA não parametrica

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