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
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
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
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
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
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
Interpretando: diferença padronizada pequena o poder de teste tem 7% de rejeitar corretamente H0.
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 ...
boxplot(cromo ~ Lab, data = labs, col = "lightpink")
library(ggplot2)
ggplot(labs, aes(x=Lab, y=cromo)) + geom_boxplot()
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 |
Observação: quando a distribuição é simétrica: a média e desvio padrão são bons dados. Quando é assimétrica: a mediana e MAD é melhor opção. Viabilidade entre duas amostras para distribuição simétrica.
# 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
Estimastivas dos efeitos principais \(\hat[\tau]_i\)
coef(mod1)
## (Intercept) LabLab2 LabLab3 LabLab4
## 23.3375 -6.5500 -4.7875 3.2000
library(agricolae)
cv.model(mod1)
## [1] 15.68622
erro = \(\epsilon_{ij}\) = erro da população
resíduo = \(\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, 1)
shapiro.test(mod1$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.96141, p-value = 0.3001
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
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
#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
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
Convertendo a variável ‘temperatura’ para a 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 ...
boxplot(dias ~ temperatura, data = turtles, col = "lightpink")
Obs: Para acessar a paleta de cores do R base, digite ‘colours()’ no console.
Interpretação: as médias de temperatura 15,20 e 25 não vão dar estatisticamente diferentes, - observando a posição das caixas - mas a de 30 é diferente . A maior variabilidade de dados ocorreu para a temperatura de 15.