Aula23nov
Camila Flávia
23/11/2019
# Pacotes que serão utlizados:
library(agricolae)
library(ggplot2)
library(asbio) # Teste Tukey de Aditividade## Loading required package: tcltklibrary(WRS2)
library(effects)## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.1. ANOVA com 2 Fatores (01 fator de Interesse e 1 Bloco)
1.1 Delineamento em Blocos Casualizados
1.1.1 Importando e verificando os Dados
# Importando os Dados
solucao <- read.table("solucao.txt", header = T)# Verificando os dados
str(solucao)## 'data.frame': 12 obs. of 3 variables:
## $ sol: Factor w/ 3 levels "S1","S2","S3": 1 1 1 1 2 2 2 2 3 3 ...
## $ dia: Factor w/ 4 levels "D1","D2","D3",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ y : int 13 22 18 39 16 24 17 44 5 4 ...summary(solucao)## sol dia y
## S1:4 D1:3 Min. : 1.00
## S2:4 D2:3 1st Qu.:11.00
## S3:4 D3:3 Median :17.50
## D4:3 Mean :18.75
## 3rd Qu.:22.50
## Max. :44.00# Boxplot: y x solução
ggplot(solucao, aes(x=sol, y= y, fill=sol)) + geom_boxplot()
Conclusão: no S1 e S2 possuem mesma media que diferem de S3. S3 provavelmente está produzindo menos bacterias em relação aos demais.
# Boxplot: y x dia
ggplot(solucao, aes(x=dia, y= y, fill=dia)) + geom_boxplot()
Conclusão: Dia 1 está produzindo menos bacterias do que os demais dias e que sua media está diferente das demais.
# Gráfico de Interação
interaction.plot(solucao$dia, solucao$sol, solucao$y, fixed = TRUE)
Conclusão: O professor se arrisca dizer que nao ha interação, qualitativamente dizendo.
# verificando de há interação Blocos x Níveis do Fator
# Teste de Atividade de Tukey
# Teste H0: Não há efeito interação
asbio::tukey.add.test(solucao$y, solucao$sol, solucao$dia)##
## Tukey's one df test for additivity
## F = 2.7732343 Denom df = 5 p-value = 0.1567331Conclusão: o resultado nao rejeita a hipotese nula H0. De fato não ha interação entre os dias e a solução. H0 = afirma que nao ha interação.
3. CV do Experimento
# CV do Pacote Agricolae
cv.model(mdic)## [1] 60.5163cv.model(mdbca)## [1] 15.675734. Análise Gráfica dos Resíduos
# Linearidade
plot(mdbca, which = 1) # Residudos x Valores previstos
# Análise Gráfica da Normalidade
plot(mdbca, which = 2) # grafico qqplot
Conclusão: Resultado satisfatório devido os valores estarem representados proximos da linha.
4.1 Testes Formais do Resíduos
#### Teste de Shapiro-Wilks
shapiro.test(mdbca$residuals)##
## Shapiro-Wilk normality test
##
## data: mdbca$residuals
## W = 0.93208, p-value = 0.4027Conclusão: Não rejeita a hipotese nula de distribuição normal.
#### Teste de Bartlett
# bartlett.test(y ~ interaction(sol, dia), data = solucao)5. Teste Tukey (Qual solução deveria ser recomendada?) use o Teste Tukey
### Teste de Tukey Agricolae
teste_tukey <- agricolae::HSD.test(mdbca,"sol", group = FALSE)
print(teste_tukey$comparison)## difference pvalue signif. LCL UCL
## S1 - S2 -2.25 0.5578 -8.626879 4.126879
## S1 - S3 15.00 0.0009 *** 8.623121 21.376879
## S2 - S3 17.25 0.0004 *** 10.873121 23.626879Conclusão: Pela avaliação do p-value (0.5578) ja se conclui que nao se rejeita a H0 devido o valor estar proximo de 0.05. Pela avaliação de LCL e UCl tambem, pois os valores estão entre o valor zero (ou seja, tem valor negativo a esquerda e positivo a direita)
A S3 é a melhor recomendada pois é a que apresenta menor crescimento bacteriano.
# Eficiencia do DBCA x DIC
qmr.di = 128.7
qmr.db = 8.6
qmr.di/qmr.db## [1] 14.96512Conclusão: O modelo DBCA (com bloco) é melhor do que o DIC (sem bloco) pois apresentou 14.96 menos residual.
6. Tutorial - Executando a Análise
# Importando e verificando os dados
sesseis <- read.csv(file = "sesseis.csv", header = TRUE)
str(sesseis)## 'data.frame': 60 obs. of 3 variables:
## $ riqueza : int 68 64 64 63 69 63 70 68 68 62 ...
## $ cobre : Factor w/ 3 levels "alta","baixa",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ orientacao: Factor w/ 2 levels "horizontal","vertical": 2 2 2 2 2 2 2 2 2 2 ...# Estimanção do modelo
sesseis.aov <- aov(riqueza ~ cobre*orientacao, data = sesseis)
summary(sesseis.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## cobre 2 3330 1665.0 192.53 < 0.0000000000000002 ***
## orientacao 1 240 240.0 27.75 0.000002463599 ***
## cobre:orientacao 2 571 285.4 33.00 0.000000000434 ***
## Residuals 54 467 8.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusão: As tres H0 são rejeitadas pois os resultados de Pr (F) estão muito abaixo de 0.05.
6.1 Teste de Tukey
TukeyHSD(sesseis.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = riqueza ~ cobre * orientacao, data = sesseis)
##
## $cobre
## diff lwr upr p adj
## baixa-alta 7.95 5.708825 10.19117 0
## nenhuma-alta 18.20 15.958825 20.44117 0
## nenhuma-baixa 10.25 8.008825 12.49117 0
##
## $orientacao
## diff lwr upr p adj
## vertical-horizontal -4 -5.522314 -2.477686 0.0000025
##
## $`cobre:orientacao`
## diff lwr upr
## baixa:horizontal-alta:horizontal 0.4 -3.48559529 4.2855953
## nenhuma:horizontal-alta:horizontal 14.2 10.31440471 18.0855953
## alta:vertical-alta:horizontal -11.7 -15.58559529 -7.8144047
## baixa:vertical-alta:horizontal 3.8 -0.08559529 7.6855953
## nenhuma:vertical-alta:horizontal 10.5 6.61440471 14.3855953
## nenhuma:horizontal-baixa:horizontal 13.8 9.91440471 17.6855953
## alta:vertical-baixa:horizontal -12.1 -15.98559529 -8.2144047
## baixa:vertical-baixa:horizontal 3.4 -0.48559529 7.2855953
## nenhuma:vertical-baixa:horizontal 10.1 6.21440471 13.9855953
## alta:vertical-nenhuma:horizontal -25.9 -29.78559529 -22.0144047
## baixa:vertical-nenhuma:horizontal -10.4 -14.28559529 -6.5144047
## nenhuma:vertical-nenhuma:horizontal -3.7 -7.58559529 0.1855953
## baixa:vertical-alta:vertical 15.5 11.61440471 19.3855953
## nenhuma:vertical-alta:vertical 22.2 18.31440471 26.0855953
## nenhuma:vertical-baixa:vertical 6.7 2.81440471 10.5855953
## p adj
## baixa:horizontal-alta:horizontal 0.9996312
## nenhuma:horizontal-alta:horizontal 0.0000000
## alta:vertical-alta:horizontal 0.0000000
## baixa:vertical-alta:horizontal 0.0587141
## nenhuma:vertical-alta:horizontal 0.0000000
## nenhuma:horizontal-baixa:horizontal 0.0000000
## alta:vertical-baixa:horizontal 0.0000000
## baixa:vertical-baixa:horizontal 0.1186003
## nenhuma:vertical-baixa:horizontal 0.0000000
## alta:vertical-nenhuma:horizontal 0.0000000
## baixa:vertical-nenhuma:horizontal 0.0000000
## nenhuma:vertical-nenhuma:horizontal 0.0705255
## baixa:vertical-alta:vertical 0.0000000
## nenhuma:vertical-alta:vertical 0.0000000
## nenhuma:vertical-baixa:vertical 0.0000651plot(TukeyHSD(sesseis.aov))
6.2 Verificação da Validade das Hipóteses do Modelo
## Normalidade dos resíduos
plot(sesseis.aov, 2)
O out-line que aparece no canto superior faz com que a H0 seja rejeitada. Mas os resultados se mantem em uma consição normal.
# Normalidade dos resíduos
shapiro.test(sesseis.aov$residuals)##
## Shapiro-Wilk normality test
##
## data: sesseis.aov$residuals
## W = 0.93677, p-value = 0.003894Conforme o valor de 0.03894, conforma a rejeição da H0 pois o valor está abaixo de 0.05.
# variância Homogênia
plot(sesseis.aov, 3)
fligner.test(riqueza ~ interaction(cobre,orientacao), data = sesseis)##
## Fligner-Killeen test of homogeneity of variances
##
## data: riqueza by interaction(cobre, orientacao)
## Fligner-Killeen:med chi-squared = 24.096, df = 5, p-value =
## 0.00020816.3 Alternativa Não-Parametrica
library(WRS2)t2way(riqueza ~ cobre*orientacao, data = sesseis)## Call:
## t2way(formula = riqueza ~ cobre * orientacao, data = sesseis)
##
## value p.value
## cobre 212.3893 0.001
## orientacao 19.3830 0.001
## cobre:orientacao 84.5326 0.001mcp2atm(riqueza ~ cobre*orientacao, data = sesseis)## Call:
## mcp2atm(formula = riqueza ~ cobre * orientacao, data = sesseis)
##
## psihat ci.lower ci.upper p-value
## cobre1 -15.00000 -19.51275 -10.48725 0.00000
## cobre2 -34.33333 -40.77738 -27.88929 0.00000
## cobre3 -19.33333 -25.29279 -13.37387 0.00000
## orientacao1 11.33333 5.94307 16.72360 0.00031
## cobre1:orientacao1 14.66667 10.15391 19.17942 0.00000
## cobre2:orientacao1 9.00000 2.55596 15.44404 0.00201
## cobre3:orientacao1 -5.66667 -11.62613 0.29279 0.022727. Gráfico dos Efeitos
library(effects)plot(allEffects(sesseis.aov))
X: Nivel de enriquecimento do cobre y: Riqueza media das especies