rm(list=ls(all=TRUE))Trabalho de Experimentação
Delineamento Inteiramente Casualizado (DIC)
Antes, limpar o “Ambiente Global” (Global Environment):
Primeiramente é necessário carregar os pacotes que serão utilizados:
library("tidyverse")
library("ggplot2")
library("PMCMRplus")
library("car")
library("agricolae")
library("rstatix")
library("emmeans")
library("ScottKnott")
library("multcomp")Dados: Um experimento de alimentação de porcos, usando quatro rações (A, B, C,D), cada uma fornecida a cinco animais escolhidos ao acaso. Os aumentos de peso observados, em quilogramas, constam a seguir.
Entrada de dados da variável resposta:
pesos = c(35, 19, 31, 15,
30, 40, 35, 46,
41, 33, 39, 27,
20, 29, 45, 27,
12, 13, 28, 30)Criar objeto com variável pesos:
pesos [1] 35 19 31 15 30 40 35 46 41 33 39 27 20 29 45 27 12 13 28 30Tratamentos (rações):
rações = factor(c(rep('A',5),
rep('B',5),
rep('C',5),
rep('D',5)))Criar objeto para tratamentos (rações):
rações [1] A A A A A B B B B B C C C C C D D D D D
Levels: A B C DMostrar a ordem das rações com função table:
table(rações)rações
A B C D
5 5 5 5 Criar o banco de dados com a variável resposta e o tratamento:
dados = data.frame(pesos,rações)
dados pesos rações
1 35 A
2 19 A
3 31 A
4 15 A
5 30 A
6 40 B
7 35 B
8 46 B
9 41 B
10 33 B
11 39 C
12 27 C
13 20 C
14 29 C
15 45 C
16 27 D
17 12 D
18 13 D
19 28 D
20 30 DAnálise descritiva:
Resumo númerico geral:
summary(dados) pesos rações
Min. :12.00 A:5
1st Qu.:25.25 B:5
Median :30.00 C:5
Mean :29.75 D:5
3rd Qu.:36.00
Max. :46.00 Resumo númerico por tratamento
tapply(pesos,rações,summary)$A
Min. 1st Qu. Median Mean 3rd Qu. Max.
15 19 30 26 31 35
$B
Min. 1st Qu. Median Mean 3rd Qu. Max.
33 35 40 39 41 46
$C
Min. 1st Qu. Median Mean 3rd Qu. Max.
20 27 29 32 39 45
$D
Min. 1st Qu. Median Mean 3rd Qu. Max.
12 13 27 22 28 30 Resumo númerico com desvio padrão e coeficiente de variação:
nrep = tapply(pesos,rações, length)
media = tapply(pesos,rações, mean)
variancia = tapply(pesos,rações, var)
desvpad = tapply(pesos,rações, sd)
coef_var = round(100*(tapply(pesos,rações,sd)/tapply(pesos,rações, mean)), 1)
data.frame(Repeticoes = nrep,
Media = media,
Variancia = variancia,
DesvioPadrao =desvpad,
CoefVariacao = coef_var,
row.names = c('A','B','C','D')) Repeticoes Media Variancia DesvioPadrao CoefVariacao
A 5 26 73.0 8.544004 32.9
B 5 39 26.5 5.147815 13.2
C 5 32 99.0 9.949874 31.1
D 5 22 76.5 8.746428 39.8Boxplot para visualizar melhor os dados e comparar os tratamentos:
car::Boxplot(pesos~rações,
las=1,
col="lightblue", xlab="",
ylab=expression("pesos"(Kg)))
points(media,col="red", pch=8)
É possível observar que as caixas das rações D e A se encontram mais abaixo quando comparadas com as demais, o valor máximo da D esta abaixo do valor mínimo da caixa B. As caixas A, B e C estão assimétricos a equerda.
Criar histograma para todos os dados:
hist(pesos, col = "red", col.axis = "blue", col.lab = "green",
xlab = "Pesos (kg)", ylab="Frequência",
main = "Frêquencia de pesos dos porcos")
O histograma não está muito próximo da aparência de um sino, sendo possivel observar uma maior frequência nos valores 25 a 35.
Teste de Normalidade de todos os dados: Histograma:
ggplot(dados, aes(x=pesos))+
geom_histogram(bins = 6,
fill = "green4",
col= "magenta2")
Graficamente podemso observar que o histograma não apresenta uma aparência próxima de um sino.
Boxplot:
ggplot(dados)+
geom_boxplot(aes(y = pesos),
fill = "deepskyblue3",
col= "turquoise1")
QQ-plot:
ggplot(dados, aes(sample = pesos)) +
stat_qq() +
stat_qq_line() +
labs(y = "pesos (kg)")
qqPlot(pesos, col="turquoise1", pch=19)
[1] 17 18Os pontos seguem aproximadamente uma linha reta, porém um valor se encontra fora da linha.
Teste de Shapiro-Wilk para Normalidade: Hipótese nula: os dados de pesos possuem uma distribuição Normal; Hipótese alternativa: os dados de pesos não possuem uma distribuição Normal.
shapiro.test(pesos)
Shapiro-Wilk normality test
data: pesos
W = 0.95835, p-value = 0.5115Como o p-valor é aproximadamente 0,51 e é maior do que 0,05, a hipótese nula não é rejeitada, logo, há evidências que comprovem distribuição normal dos dados dos pesos.
Teste da homogeneidade de variâncias: Hipotese nula: As rações possuem variâncias homogêneas para pesos dos porcos; Hipotese alternativa: As variâncias dos pesos dos porcos em cada ração não são homogêneas.
leveneTest(pesos ~ rações) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 0.3324 0.802
16 Como o p-valor é maior que 0,05, não há evidências para rejeitar a hipótese nula, logo, há evidências que comprovem que as são variâncias são homogêneas.
Tabela da análise de variâncias: Testar as hipóteses: Hip. nula: a média dos pesos é a mesma para todas as rações (tratamentos); Hip. alternativa: em pelo menos uma ração a média dos pesos se diferem.
Ajuste do modelo aos dados:
mod_aov = aov(pesos ~ rações, data = dados)
mod_aovCall:
aov(formula = pesos ~ rações, data = dados)
Terms:
rações Residuals
Sum of Squares 823.75 1100.00
Deg. of Freedom 3 16
Residual standard error: 8.291562
Estimated effects may be unbalancedsummary(mod_aov) #produz a tabela da ANOVA Df Sum Sq Mean Sq F value Pr(>F)
rações 3 823.8 274.58 3.994 0.0267 *
Residuals 16 1100.0 68.75
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como valor-p é menor do que 0,05, não há evidências para aceitar a hipótese nula, logo, em pelo menos uma ração a média dos pesos se diferem.
Contrução do gráfico “resíduos padronizados x valores ajustados”: Resíduos ordinários:
residuos = mod_aov$residuals
residuos 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
9 -7 5 -11 4 1 -4 7 2 -6 7 -5 -12 -3 13 5 -10 -9 6 8 dados = data.frame(dados, residuos)Resíduos X Tratamento (rações):
ggplot(data = dados,
aes(x = rações, y = residuos)) +
geom_jitter(width = 0.05,
color="violetred") +
labs(y = "Residuos ordinarios",
x = "rações") +
stat_summary(show.legend = T,
fun = "mean",
geom = "crossbar",
width = 0.25,
linewidth = 0.25,
color = "turquoise1")
Valores ajustados (teóricos) do modelo:
ajustados = mod_aov$fitted.values
ajustados 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
26 26 26 26 26 39 39 39 39 39 32 32 32 32 32 22 22 22 22 22 Gráfico dos Resíduos X valores ajustados:
plot(ajustados,
residuos,
main = "Residuos X Valores ajustados")
Resíduos padronizados: (verificar se existe algum valor resíduo maior que 3 ou menor que -3):
residuos_p = rstandard(mod_aov)Gráfico Resíduos Padronizados x Valores ajustados:
plot(ajustados,
residuos_p,
xlab = "valores ajustados",
ylab = "resíduos padronizados",
ylim=c(-3.5,3.5),
pch=16,
col="green2",
las=1)
abline(h=seq(-3,3,1), lty=3, col="plum1")
abline(h=c(-3,3), lty=4, col="magenta")
title("Residuos Padronizados X Valores ajustados")
Os resíduos estão bem ajustados porque estão dispersos e próximos do zero, não havendo valores discrepantes que ultrapassem a linha tracejada no entanto, é possível observar um padrão entre os valores.
Gráfico para normalidade dos resíduos:
hist(residuos)
ggplot(dados, aes(x=residuos))+
geom_histogram(bins = 5,
fill = "green",
col= "maroon4")+
labs(y = "Frequencia absoluta")
O histograma apresenta uma aparência próxima a de um sino, indicando uma possível normalidade dos resíduos.
qqnorm(residuos,
ylab = "Residuos",
main = "Grafico Normal de Probabilidade dos Residuos",
pch=19,
col="darkorchid3")
qqline(residuos,
col="turquoise1") 
ggplot(dados, aes(sample = residuos)) +
stat_qq() +
stat_qq_line() +
labs(y = "Residuos ordinarios")
qqPlot(residuos, col="deeppink1", pch=19)
[1] 15 13Os resíduos parecem seguir uma distribuição normal, os pontos seguem aproximadamente uma linha reta.
Teste de hipótese para a normalidade dos resíduos com shapiro: Hipótese nula: os resíduos tem distribuição normal; Hipótese alternativa: os resíduos não tem distribuição normal:
shapiro.test(residuos)
Shapiro-Wilk normality test
data: residuos
W = 0.93874, p-value = 0.227Como valor-p é maior do que 0,05, não há evidências para rejeitar a hipótese nula, logo, os resíduos tem distribuição normal.
Teste de Homogeneidade de variâncias para os resíduos: Ho: As variâncias são homogêneas; Ha: As variâncias não são homogêneas.
leveneTest(residuos ~ rações)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 0.3324 0.802
16 Como o valor-p é maior que 0,05, não há evidências para rejeitar a hipótese nula, logo, as variâncias são homogêneas.
Teste de comparações multiplas de Tukey: Teste de Tukey
TukeyHSD(mod_aov) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = pesos ~ rações, data = dados)
$rações
diff lwr upr p adj
B-A 13 -2.003315 28.003315 0.1018285
C-A 6 -9.003315 21.003315 0.6687032
D-A -4 -19.003315 11.003315 0.8698923
C-B -7 -22.003315 8.003315 0.5553529
D-B -17 -32.003315 -1.996685 0.0237354
D-C -10 -25.003315 5.003315 0.2640642TukeyHSD(mod_aov,
ordered = TRUE,
conf.level = 0.95) Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered
Fit: aov(formula = pesos ~ rações, data = dados)
$rações
diff lwr upr p adj
A-D 4 -11.003315 19.00331 0.8698923
C-D 10 -5.003315 25.00331 0.2640642
B-D 17 1.996685 32.00331 0.0237354
C-A 6 -9.003315 21.00331 0.6687032
B-A 13 -2.003315 28.00331 0.1018285
B-C 7 -8.003315 22.00331 0.5553529res_tukey1 = TukeyHSD(mod_aov,
ordered = TRUE,
conf.level = 0.95)
plot(res_tukey1)
Para melhorar a visualização e embelezar o gráfico. define-se o padrão que já existe:
def.par = par() Criar um padrão para este gráfico:
par(mai=c(1.3, #margem inferior
1.4, #margem da esquerda
1.1, #margem do topo
.5)) #margem da direita
plot(res_tukey1,
col.axis = "slategray4",
col = "cyan4",
las=1, #muda a posicao dos nomes do eixo Y
cex.axis = 0.75) #tamanho menor dos nomes
Considerando o valor-p, notamos que as rações B e D são menores que o nível de significância adotado, ou seja, é menor do que 0,05. Ao analisar o gráfico nota-se que o valor zero não está contido no intervalo de confiança das duas rações.
Comparando as médias:
xtabs(~rações) #Repetições por nível de cultivarrações
A B C D
5 5 5 5 df.residual(mod_aov) #Graus de liberdade do resíduo [1] 16glres = df.residual(mod_aov)
deviance(mod_aov)/glres #Quadrado médio do resíduo[1] 68.75QMRes = deviance(mod_aov)/glres
res_tukey2 = HSD.test(pesos,
rações,
glres,
QMRes,
alpha=0.05)
res_tukey2$statistics
MSerror Df Mean CV MSD
68.75 16 29.75 27.8708 15.00331
$parameters
test name.t ntr StudentizedRange alpha
Tukey rações 4 4.046093 0.05
$means
pesos std r Min Max Q25 Q50 Q75
A 26 8.544004 5 15 35 19 30 31
B 39 5.147815 5 33 46 35 40 41
C 32 9.949874 5 20 45 27 29 39
D 22 8.746428 5 12 30 13 27 28
$comparison
NULL
$groups
pesos groups
B 39 a
C 32 ab
A 26 ab
D 22 b
attr(,"class")
[1] "group"plot(res_tukey2,
cex.axis=0.85,
ylab = " média de pesos (kg)",
xlab = "rações",
las=1)
Ao comparar as médias em ordem decrescente da média maior para menor, temos ração B, C, A e por fim D. A ração B possui média maior porém sua média se difere significamente somente da ração D. E a D só se difere da B.
Apresentação final dos resultados. Função “cld()” pertence ao pacote “multcomp”:
res_medias = emmeans(mod_aov, spec = ~rações) %>%
cld(Letters = letters, sort = TRUE, reverse = TRUE) %>%
as.data.frame()
res_medias rações emmean SE df lower.CL upper.CL .group
B 39 3.708099 16 31.13918 46.86082 a
C 32 3.708099 16 24.13918 39.86082 ab
A 26 3.708099 16 18.13918 33.86082 ab
D 22 3.708099 16 14.13918 29.86082 b
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. ggplot(data = res_medias,
mapping = aes(y = reorder(rações, emmean))) +
geom_errorbarh(mapping = aes(xmin = lower.CL, xmax = upper.CL),
height = 0.06) +
labs(y = "rações",
x = "média pesos (kg)") +
geom_point(mapping = aes(x = emmean)) +
geom_text(mapping = aes(x = emmean,
label = sprintf("%0.1f %s",
emmean,
trimws(.group))),
vjust = -1)
Por fim, ao comparar os experimentos das rações conclui-se que pelo menos uma média é diferente, e de acordo com os testes a ração B se difere significativamente da ração B. Ao analisar as médias a ração B gerou mais aumento de peso nos porcos.
Delineamento em Bloscos Casualizados (DBC)
Dados: Em um experimento de competição de variedades de batatinha, feito pelo Engenheiro Agronomo Oscar A. GAray em Balcare, Argentina, em blocos casualizados, as produção obtidas, em t/ha, foram as seguintes.
Foram utilizadas quatro repetições e oito tratamentos: T1: Kennebec T2: Huinkul T3: S. Rafaela T4: Buena Vista T5: B 25-50 E T6: B 1-52 T7: B 116-51 T8: B 72-53 A
Entrada de dados para variável resposta:
produção = c(9.2, 13.4, 11.0, 9.2,
21.1, 27.0, 26.4, 25.7,
22.6, 29.9, 24.2, 25.1,
15.4, 11.9, 10.1, 12.3,
12.7, 18.0, 18.2, 17.1,
20.0, 21.1, 20.0, 28.0,
23.1, 24.2, 26.4, 16.3,
18.0, 24.6, 24.0, 24.6)Número de unidades experimentais:
length(produção) [1] 32Tratamentos: criar a variável qualitativa para representar os tratamentos:
tratamento = factor(rep(c("Kennebec",
"Huinkul",
"S. Rafaela",
"Buena Vista",
"B 25-50 E",
"B 1-52",
"B 116-51",
"B 72-53 A"), each=4))Criar objeto com variável tratamentos:
tratamento [1] Kennebec Kennebec Kennebec Kennebec Huinkul Huinkul
[7] Huinkul Huinkul S. Rafaela S. Rafaela S. Rafaela S. Rafaela
[13] Buena Vista Buena Vista Buena Vista Buena Vista B 25-50 E B 25-50 E
[19] B 25-50 E B 25-50 E B 1-52 B 1-52 B 1-52 B 1-52
[25] B 116-51 B 116-51 B 116-51 B 116-51 B 72-53 A B 72-53 A
[31] B 72-53 A B 72-53 A
8 Levels: B 1-52 B 116-51 B 25-50 E B 72-53 A Buena Vista Huinkul ... S. RafaelaMostrar a ordem dos tratamentos:
table(tratamento)tratamento
B 1-52 B 116-51 B 25-50 E B 72-53 A Buena Vista Huinkul
4 4 4 4 4 4
Kennebec S. Rafaela
4 4 Não estão em ordem, colocar na ordem:
tratamento = factor(tratamento, c("Kennebec",
"Huinkul",
"S. Rafaela",
"Buena Vista",
"B 25-50 E",
"B 1-52",
"B 116-51",
"B 72-53 A"))table(tratamento)tratamento
Kennebec Huinkul S. Rafaela Buena Vista B 25-50 E B 1-52
4 4 4 4 4 4
B 116-51 B 72-53 A
4 4 Criar o objeto para representar os blocos. Usar a função “factor()”:
bloco = factor(rep(c("B1",
"B2",
"B3",
"B4"), times = 8))Visualizar o comprimento do objeto bloco:
bloco [1] B1 B2 B3 B4 B1 B2 B3 B4 B1 B2 B3 B4 B1 B2 B3 B4 B1 B2 B3 B4 B1 B2 B3 B4 B1
[26] B2 B3 B4 B1 B2 B3 B4
Levels: B1 B2 B3 B4length(bloco) [1] 32Organizar os objetos em um banco de dados (data frame):
dados = data.frame(bloco, tratamento, produção)
dados bloco tratamento produção
1 B1 Kennebec 9.2
2 B2 Kennebec 13.4
3 B3 Kennebec 11.0
4 B4 Kennebec 9.2
5 B1 Huinkul 21.1
6 B2 Huinkul 27.0
7 B3 Huinkul 26.4
8 B4 Huinkul 25.7
9 B1 S. Rafaela 22.6
10 B2 S. Rafaela 29.9
11 B3 S. Rafaela 24.2
12 B4 S. Rafaela 25.1
13 B1 Buena Vista 15.4
14 B2 Buena Vista 11.9
15 B3 Buena Vista 10.1
16 B4 Buena Vista 12.3
17 B1 B 25-50 E 12.7
18 B2 B 25-50 E 18.0
19 B3 B 25-50 E 18.2
20 B4 B 25-50 E 17.1
21 B1 B 1-52 20.0
22 B2 B 1-52 21.1
23 B3 B 1-52 20.0
24 B4 B 1-52 28.0
25 B1 B 116-51 23.1
26 B2 B 116-51 24.2
27 B3 B 116-51 26.4
28 B4 B 116-51 16.3
29 B1 B 72-53 A 18.0
30 B2 B 72-53 A 24.6
31 B3 B 72-53 A 24.0
32 B4 B 72-53 A 24.6Estrutura do banco de dados criado:
str(dados)'data.frame': 32 obs. of 3 variables:
$ bloco : Factor w/ 4 levels "B1","B2","B3",..: 1 2 3 4 1 2 3 4 1 2 ...
$ tratamento: Factor w/ 8 levels "Kennebec","Huinkul",..: 1 1 1 1 2 2 2 2 3 3 ...
$ produção : num 9.2 13.4 11 9.2 21.1 27 26.4 25.7 22.6 29.9 ...Resumo numérico do banco de dados:
summary(dados) bloco tratamento produção
B1:8 Kennebec :4 Min. : 9.20
B2:8 Huinkul :4 1st Qu.:14.90
B3:8 S. Rafaela :4 Median :20.55
B4:8 Buena Vista:4 Mean :19.71
B 25-50 E :4 3rd Qu.:24.60
B 1-52 :4 Max. :29.90
(Other) :8 Análise descritva da variável resposta (produção) Análise por tratamento: Resumo númerico:
tapply(produção, tratamento, summary)$Kennebec
Min. 1st Qu. Median Mean 3rd Qu. Max.
9.2 9.2 10.1 10.7 11.6 13.4
$Huinkul
Min. 1st Qu. Median Mean 3rd Qu. Max.
21.10 24.55 26.05 25.05 26.55 27.00
$`S. Rafaela`
Min. 1st Qu. Median Mean 3rd Qu. Max.
22.60 23.80 24.65 25.45 26.30 29.90
$`Buena Vista`
Min. 1st Qu. Median Mean 3rd Qu. Max.
10.10 11.45 12.10 12.43 13.07 15.40
$`B 25-50 E`
Min. 1st Qu. Median Mean 3rd Qu. Max.
12.70 16.00 17.55 16.50 18.05 18.20
$`B 1-52`
Min. 1st Qu. Median Mean 3rd Qu. Max.
20.00 20.00 20.55 22.27 22.82 28.00
$`B 116-51`
Min. 1st Qu. Median Mean 3rd Qu. Max.
16.30 21.40 23.65 22.50 24.75 26.40
$`B 72-53 A`
Min. 1st Qu. Median Mean 3rd Qu. Max.
18.0 22.5 24.3 22.8 24.6 24.6 Número de repetições
nrep = tapply(produção, tratamento, length)Média:
medias_produção_trat = round(tapply(produção, tratamento, mean),1)
medias_produção_trat Kennebec Huinkul S. Rafaela Buena Vista B 25-50 E B 1-52
10.7 25.0 25.4 12.4 16.5 22.3
B 116-51 B 72-53 A
22.5 22.8 Mediana:
mediana_produção_trat = round(tapply(produção, tratamento, median),1)
mediana_produção_trat Kennebec Huinkul S. Rafaela Buena Vista B 25-50 E B 1-52
10.1 26.0 24.6 12.1 17.6 20.6
B 116-51 B 72-53 A
23.6 24.3 Variância:
var_produção_trat = round(tapply(produção, tratamento, var),1)
var_produção_trat Kennebec Huinkul S. Rafaela Buena Vista B 25-50 E B 1-52
4.0 7.2 9.9 4.8 6.6 14.8
B 116-51 B 72-53 A
19.0 10.3 Desvio padrão:
dp_produção_trat = round(sqrt(var_produção_trat),1)
dp_produção_trat Kennebec Huinkul S. Rafaela Buena Vista B 25-50 E B 1-52
2.0 2.7 3.1 2.2 2.6 3.8
B 116-51 B 72-53 A
4.4 3.2 Coeficiente de variação:
coef_var = round(100*(dp_produção_trat/medias_produção_trat),1)
data.frame(Repeticoes = nrep,
Media = medias_produção_trat,
Mediana = mediana_produção_trat,
Variancia = var_produção_trat,
DesvioPadrao =dp_produção_trat,
CoefVariacao = coef_var,
row.names = c("Kennebec",
"Huinkul",
"S. Rafaela",
"Buena Vista",
"B 25-50 E",
"B 1-52",
"B 116-51",
"B 72-53 A")) Repeticoes Media Mediana Variancia DesvioPadrao CoefVariacao
Kennebec 4 10.7 10.1 4.0 2.0 18.7
Huinkul 4 25.0 26.0 7.2 2.7 10.8
S. Rafaela 4 25.4 24.6 9.9 3.1 12.2
Buena Vista 4 12.4 12.1 4.8 2.2 17.7
B 25-50 E 4 16.5 17.6 6.6 2.6 15.8
B 1-52 4 22.3 20.6 14.8 3.8 17.0
B 116-51 4 22.5 23.6 19.0 4.4 19.6
B 72-53 A 4 22.8 24.3 10.3 3.2 14.0dados |>
group_by(tratamento) |>
rstatix::get_summary_stats(produção, type = "full")# A tibble: 8 × 14
tratamento variable n min max median q1 q3 iqr mad mean
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Kennebec produção 4 9.2 13.4 10.1 9.2 11.6 2.4 1.33 10.7
2 Huinkul produção 4 21.1 27 26.0 24.6 26.6 2 0.964 25.0
3 S. Rafaela produção 4 22.6 29.9 24.6 23.8 26.3 2.5 1.85 25.4
4 Buena Vista produção 4 10.1 15.4 12.1 11.4 13.1 1.62 1.63 12.4
5 B 25-50 E produção 4 12.7 18.2 17.6 16 18.0 2.05 0.815 16.5
6 B 1-52 produção 4 20 28 20.6 20 22.8 2.82 0.815 22.3
7 B 116-51 produção 4 16.3 26.4 23.6 21.4 24.8 3.35 2.45 22.5
8 B 72-53 A produção 4 18 24.6 24.3 22.5 24.6 2.1 0.445 22.8
# ℹ 3 more variables: sd <dbl>, se <dbl>, ci <dbl>Como as unidades experimentais são muitas, é melhor visualizar as medidas descritivas em gráfico para comparar os tratamentos.
Criar boxplot para comparar os tratamentos:
car::Boxplot(produção~tratamento,
col = "turquoise1",
border = "orchid4",
ylab = "produção (t/ha)",
main = c("produção (t/ha)",
"variedades(batatinhas)"),
cex.main=0.75)
Com o boxplot é possível observar que os tratamentos “Kennebec”, “Buena Vista” e “B 25-50 E” estão com as ciaxas mais abaixo, indicando menores valores, o tratamento “S. Rafaela” é o que o maior valor e esta com a caixa mais acima, é possível observar que as caixas não possuem uma simetria, ocorrendo a presença de caixas tanto assimetricas positivas quanto negativas.
Criar histograma para todos os dados:
hist(produção, col = "red", col.axis = "blue", col.lab = "green",
xlab = "produção (t/ha)", ylab="Frequência",
main = "frêquencia da produção de batatinha")
O histograma apresenta uma distribuição dos valores com maior frequências entre 20 a 25 e menor entrre 5 a 10.
Análise descritiva dos blocos:
medias_produção_blocos = tapply(produção, bloco, mean)
medias_produção_blocos B1 B2 B3 B4
17.7625 21.2625 20.0375 19.7875 mediana_produção_blocos = tapply(produção, bloco, median)
mediana_produção_blocos B1 B2 B3 B4
19.00 22.65 22.00 20.85 var_produção_blocos = tapply(produção, bloco, var)
var_produção_blocos B1 B2 B3 B4
24.65982 41.06268 42.54268 48.76125 Ao analisar os blocos nota-se maior média no bloco 2 e menor no bloco 1, sendo que ocorreu maior variação no bloco 4.
Teste de Normalidade da variável resposta: Histograma:
ggplot(dados, aes(x=produção))+
geom_histogram(bins = 6,
fill = "magenta1",
col= "cyan1")
Graficamente, o histograma apresenta aparência próxima a de um sino, memso que os penúltimos valores possuam um aumento na frequência.
Teste de Normalidade com QQ-normal: QQ-plot:
ggplot(dados, aes(sample = produção)) +
stat_qq() +
stat_qq_line() +
labs(y = "produção (t/ha)")
Os pontos seguem aproximadamente uma linha reta, indicando normalidade dos valores.
Função “qqplot()” do pacote car:
qqPlot(produção, col="magenta1", pch=19)
[1] 1 4Todos os pontos estão dentro dos limites do envelope simulado, indicando normalidade dos valores.
Teste de Shapiro-Wilk para Normalidade: Hipótese nula: os dados de produção possuem uma distribuição Normal; Hipótese alternativa: os dados de produção não possuem uma distribuição Normal.
shapiro.test(produção)
Shapiro-Wilk normality test
data: produção
W = 0.94276, p-value = 0.0897Como o p-valor é maior do que 0,05, a hipótese nula não é rejeitada, logo, há evidências que comprovem distribuição normal dos dados de produção.
Teste de homogeneidade de variâncias dos tratamentos com Teste de Levene: Ho: As variâncias são homogêneas; Ha: As variâncias não são homogêneas.
leveneTest(produção~tratamento)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 7 0.14 0.9939
24 O valor-p foi maior que 0,05, portanto não há evidências para rejeitar a hipótese nula.
Tabela análise de variância do modelo: Ajuste do modelo Testar as hipóteses: Hip. nula: a produção média é o mesmo para todos os tratamentos; Hip. alternativa: em pelo menos um tratamento a produção média difere.
mod_aov = aov(produção ~ bloco + tratamento)
mod_aovCall:
aov(formula = produção ~ bloco + tratamento)
Terms:
bloco tratamento Residuals
Sum of Squares 50.530 919.720 179.465
Deg. of Freedom 3 7 21
Residual standard error: 2.923346
Estimated effects may be unbalancedTabela da ANOVA:
anova(mod_aov)Analysis of Variance Table
Response: produção
Df Sum Sq Mean Sq F value Pr(>F)
bloco 3 50.53 16.843 1.9709 0.1493
tratamento 7 919.72 131.389 15.3744 5.723e-07 ***
Residuals 21 179.46 8.546
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Valor-p < 0,05 então há evidências para rejeitar a hipótese nula, portanto, em pelo menos um tratamento a produção média se difere.
Contrução do gráfico “resíduos padronizados x valores ajustados” Resíduos ordinários:
residuos = mod_aov$residuals
residuos 1 2 3 4 5 6 7 8 9 10 11
0.450 1.150 -0.025 -1.575 -2.000 0.400 1.025 0.575 -0.900 2.900 -1.575
12 13 14 15 16 17 18 19 20 21 22
-0.425 4.925 -2.075 -2.650 -0.200 -1.850 -0.050 1.375 0.525 -0.325 -2.725
23 24 25 26 27 28 29 30 31 32
-2.600 5.650 2.550 0.150 3.575 -6.275 -2.850 0.250 0.875 1.725 Resumo numérico dos resíduos:
summary(residuos) Min. 1st Qu. Median Mean 3rd Qu. Max.
-6.2750 -1.6438 0.0625 0.0000 1.0562 5.6500 Resíduos X Tratamento:
ggplot(data = dados,
aes(x = tratamento, y = produção)) +
geom_jitter(width = 0.05,
color="purple3") +
labs(y = "Residuos ordinarios",
x = "Tratamento") +
stat_summary(show.legend = T,
fun = "mean",
geom = "crossbar",
width = 0.25,
linewidth = 0.25,
color = "blue")
Resumo numérico dos resíduos:
summary(residuos) Min. 1st Qu. Median Mean 3rd Qu. Max.
-6.2750 -1.6438 0.0625 0.0000 1.0562 5.6500 Valores ajustados (teóricos):
ajustados = mod_aov$fitted.valuesGráfico resíduos X valores ajustados:
plot(ajustados, #eixo X
residuos, #eixo Y
ylab="residuos",
xlab="valores ajustados",
main="Residuos X Valores ajustados",
pch = 16,
col="darkorange2")
residuos_p = rstandard(mod_aov)Gráfico resíduos padronizados com valores ajustados:
plot(ajustados,
residuos_p,
ylab="residuos padronizados",
xlab="valores ajustados",
ylim = c(-3.5,3.5),
main="Residuos padronizados X valores ajustados",
pch = 16,
col="darkviolet")
abline(h=c(-3,3),
col="darkgreen",
lty = "dotted")
Os resíduos estão aleatoriamente distribuídos em torno do valor zero e entre -3 e +3, não havendo valores discrepantes que ultrapassem a linha tracejada.
Gráfico para verificar a normalidade dos Resíduos:
qqnorm(residuos,
ylab="Residuos",
main="Grafico Normalidade dos Residuos",
pch = 16,
col="dodgerblue3")
qqline(residuos)
QQ-plot:
ggplot(dados, aes(sample = residuos)) +
stat_qq() +
stat_qq_line() +
labs(y = "Residuos ordinarios")
qqPlot(residuos, pch = 16, col = "maroon")
[1] 28 24Os resíduos parecem seguir uma distribuição normal, os pontos seguem aproximadamente uma linha reta, no entanto existe um ponto ultrapassando o limite.
Teste para normalidade dos resíduos com shapiro-wilk: Hipótese nula: os resíduos tem distribuição normal; Hipótese alternativa: os resíduos não tem distribuição normal:
shapiro.test(residuos)
Shapiro-Wilk normality test
data: residuos
W = 0.96895, p-value = 0.471O valor-p é maior do que 0,05, portanto não há evidências para rejeitar a hipótese nula de que os resíduos tem distribuição normal.
Teste de homogeneidade de variâncias dos resíduos:
leveneTest(residuos~tratamento)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 7 0.8304 0.5725
24 Como o valor-p é maior que 0,05, não há evidências para rejeitar a hipótese nula, logo, as variâncias são homogêneas. Teste de Tukey:
res_tukey1 = TukeyHSD(mod_aov, "tratamento", ord=T)
res_tukey1 Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered
Fit: aov(formula = produção ~ bloco + tratamento)
$tratamento
diff lwr upr p adj
Buena Vista-Kennebec 1.725 -5.2084131 8.658413 0.9888092
B 25-50 E-Kennebec 5.800 -1.1334131 12.733413 0.1463333
B 1-52-Kennebec 11.575 4.6415869 18.508413 0.0003326
B 116-51-Kennebec 11.800 4.8665869 18.733413 0.0002605
B 72-53 A-Kennebec 12.100 5.1665869 19.033413 0.0001884
Huinkul-Kennebec 14.350 7.4165869 21.283413 0.0000176
S. Rafaela-Kennebec 14.750 7.8165869 21.683413 0.0000117
B 25-50 E-Buena Vista 4.075 -2.8584131 11.008413 0.5219209
B 1-52-Buena Vista 9.850 2.9165869 16.783413 0.0022078
B 116-51-Buena Vista 10.075 3.1415869 17.008413 0.0017232
B 72-53 A-Buena Vista 10.375 3.4415869 17.308413 0.0012384
Huinkul-Buena Vista 12.625 5.6915869 19.558413 0.0001072
S. Rafaela-Buena Vista 13.025 6.0915869 19.958413 0.0000701
B 1-52-B 25-50 E 5.775 -1.1584131 12.708413 0.1495926
B 116-51-B 25-50 E 6.000 -0.9334131 12.933413 0.1223292
B 72-53 A-B 25-50 E 6.300 -0.6334131 13.233413 0.0926538
Huinkul-B 25-50 E 8.550 1.6165869 15.483413 0.0091664
S. Rafaela-B 25-50 E 8.950 2.0165869 15.883413 0.0059314
B 116-51-B 1-52 0.225 -6.7084131 7.158413 1.0000000
B 72-53 A-B 1-52 0.525 -6.4084131 7.458413 0.9999952
Huinkul-B 1-52 2.775 -4.1584131 9.708413 0.8721456
S. Rafaela-B 1-52 3.175 -3.7584131 10.108413 0.7803373
B 72-53 A-B 116-51 0.300 -6.6334131 7.233413 0.9999999
Huinkul-B 116-51 2.550 -4.3834131 9.483413 0.9124298
S. Rafaela-B 116-51 2.950 -3.9834131 9.883413 0.8349125
Huinkul-B 72-53 A 2.250 -4.6834131 9.183413 0.9523508
S. Rafaela-B 72-53 A 2.650 -4.2834131 9.583413 0.8956159
S. Rafaela-Huinkul 0.400 -6.5334131 7.333413 0.9999993plot(res_tukey1,
cex.axis=0.6,
las=1)
De 28 grupos pelos menos 12 médias se constrastam, para analisar melhor os dados é interessante fazer o teste de Ducan.
Testes de comparações múltiplas:
require(agricolae)
xtabs(~tratamento) #Repetições por nível de TRATtratamento
Kennebec Huinkul S. Rafaela Buena Vista B 25-50 E B 1-52
4 4 4 4 4 4
B 116-51 B 72-53 A
4 4 glres = mod_aov$df.residual #Graus de liberdade dos resíduos
QMRes = deviance(mod_aov)/glres #Quadrado médio do ResíduoTeste de Tukey:
res_tukey2= HSD.test(produção,
tratamento,
glres,
QMRes,
alpha=0.05)
res_tukey2$statistics
MSerror Df Mean CV MSD
8.545952 21 19.7125 14.82991 6.933413
$parameters
test name.t ntr StudentizedRange alpha
Tukey tratamento 8 4.743477 0.05
$means
produção std r Min Max Q25 Q50 Q75
B 1-52 22.275 3.851731 4 20.0 28.0 20.00 20.55 22.825
B 116-51 22.500 4.355074 4 16.3 26.4 21.40 23.65 24.750
B 25-50 E 16.500 2.578113 4 12.7 18.2 16.00 17.55 18.050
B 72-53 A 22.800 3.212476 4 18.0 24.6 22.50 24.30 24.600
Buena Vista 12.425 2.202082 4 10.1 15.4 11.45 12.10 13.075
Huinkul 25.050 2.686385 4 21.1 27.0 24.55 26.05 26.550
Kennebec 10.700 1.989975 4 9.2 13.4 9.20 10.10 11.600
S. Rafaela 25.450 3.141656 4 22.6 29.9 23.80 24.65 26.300
$comparison
NULL
$groups
produção groups
S. Rafaela 25.450 a
Huinkul 25.050 a
B 72-53 A 22.800 ab
B 116-51 22.500 ab
B 1-52 22.275 ab
B 25-50 E 16.500 bc
Buena Vista 12.425 c
Kennebec 10.700 c
attr(,"class")
[1] "group"Representação gráfica:
plot(res_tukey2,
xlab = "Tratamento",
ylab = "produção (t/ha)",
cex.axis=0.85,
las=1)
Como as médias estão em ordem decrescente o tratamento S. Rafaela obteve maior média de produtividade de batatinha, seguido por Huinkul. Os tratamentos Kennebeck e Buena Vista foram os que tiveram médias menores. Os outros tratamentos tiveram médias próximas das primeiras e últimas médias.
Apresentação final dos resultados:
res_medias = emmeans(mod_aov, spec = ~tratamento) %>%
multcomp:: cld(Letters = letters, sort = TRUE, reverse = TRUE) %>%
as.data.frame()
res_medias tratamento emmean SE df lower.CL upper.CL .group
S. Rafaela 25.450 1.461673 21 22.410284 28.48972 a
Huinkul 25.050 1.461673 21 22.010284 28.08972 a
B 72-53 A 22.800 1.461673 21 19.760284 25.83972 ab
B 116-51 22.500 1.461673 21 19.460284 25.53972 ab
B 1-52 22.275 1.461673 21 19.235284 25.31472 ab
B 25-50 E 16.500 1.461673 21 13.460284 19.53972 bc
Buena Vista 12.425 1.461673 21 9.385284 15.46472 c
Kennebec 10.700 1.461673 21 7.660284 13.73972 c
Results are averaged over the levels of: bloco
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 8 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. ggplot(data = res_medias,
mapping = aes(y = reorder(tratamento, emmean))) +
geom_errorbarh(mapping = aes(xmin = lower.CL, xmax = upper.CL),
height = 0.06) +
labs(y = "Tratamento",
x = "produção medio (g)") +
geom_point(mapping = aes(x = emmean)) +
geom_text(mapping = aes(x = emmean,
label = sprintf("%0.1f %s",
emmean,
trimws(.group))),
vjust = -1)
Ao comparar a produtividade das variedades conclui-se que existe um número grande de médias contrastando, e de acordo com os testes as variedades S. Rafaela e Huinkul não se diferem significamente na produtividade com médias altas. E como se trata de competição de produtividade as variedades Kennebeck e Buena Vista foram as que menos se destacaram.
Experimento Fatorial
Dados: Considerando um experimento inteiramento casualizado, com cinco repetições no esquema fatorial 4 x 2, para testar os efeitos da aplicação de dicloroisocianurato de sódio (DUP) em 4 épocas diferentes de aplicação, em soja inoculada com Rhizobium e não inoculada, sobre o número de nódulos.
Inoculação: NI (Soja não inoculada); IN (Soja inoculada).
Épocas: Plantio; V1+15; V3+15 e R1+15.
Há interesse em testar hipóteses sobre os efeitos do fator época: H0: o número de nódulos será igual para todas as épocas; Ha: o número de nódulos se difere em pelo menos uma época.
Há interesse em testar hipóteses sobre os efeitos do fator inoculação: H0: o número de nódulos será igual para as duas inoculações; Ha: o número de nódulos se difere em pelo menos uma inoculação.
Há interesse em testar hipótes se há interação entre os fatores: H0: não há intereção para todas as combinações de fatores; Ha: há interação em pelo menos uma combinação de fatores.
Entrada de dados:
NN=c(339,332,163,230,300,
163,172,123,083,161,
171,069,095,046,079,
335,235,217,174,222,
284,136,225,098,110,
082,038,092,053,046,
196,252,346,468,258,
032,038,063,048,160)Criar um objeto com os fatores de inoculação:
(Inoculacao=rep(c("IN","NI"),e=20)) [1] "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN" "IN"
[16] "IN" "IN" "IN" "IN" "IN" "NI" "NI" "NI" "NI" "NI" "NI" "NI" "NI" "NI" "NI"
[31] "NI" "NI" "NI" "NI" "NI" "NI" "NI" "NI" "NI" "NI"Criar um objeto com o fatores de época com 5 repetições:
(epoca=rep(c("Plantio","V1+15","V3+15","R1+15"),e=5,2)) [1] "Plantio" "Plantio" "Plantio" "Plantio" "Plantio" "V1+15" "V1+15"
[8] "V1+15" "V1+15" "V1+15" "V3+15" "V3+15" "V3+15" "V3+15"
[15] "V3+15" "R1+15" "R1+15" "R1+15" "R1+15" "R1+15" "Plantio"
[22] "Plantio" "Plantio" "Plantio" "Plantio" "V1+15" "V1+15" "V1+15"
[29] "V1+15" "V1+15" "V3+15" "V3+15" "V3+15" "V3+15" "V3+15"
[36] "R1+15" "R1+15" "R1+15" "R1+15" "R1+15" Renomeando o objeto inoculação e transformando em fator:
F1=as.factor(Inoculacao)Renomeando o objeto época e transformando em fator:
F2=as.factor(epoca)Criando um objeto com os tratamentos, cada fator de inoculação, combinado com cada fator de época:
Trat=paste(F1,F2)Criando um banco de dados com os tratamentos e seus respectivos números de nódulos:
dados=data.frame(F1,F2,resp=NN)
X="";Y="Número de nódulos"Estatística descritiva dos números de nódulos:
Média:
Média = with(dados, mean(resp))Desvio padrão:
Desvio = with(dados, sd(resp))Variância:
Variancia = with(dados, var(resp))Coeficiente de variação:
CV = Desvio / Média * 100Criando uma tabela com os resultados das medidas descritivas dos números de nódulos:
desc = cbind(Média, Desvio, Variancia, CV)
desc Média Desvio Variancia CV
[1,] 168.35 106.8336 11413.41 63.45921A partir das medidas descritivas de todos os valores de número de nódulos encontramos uma média de 168,35, com um desvio padrão de 106,8336 e uma variância de 11413,41, resultando em um coeficiente de variação aparentemente alto, sendo de aproximadamente 63%.
Medidas descritivas por inoculação:
Média:
MediaA = with(dados, tapply(resp, F1, mean))Variância:
VarianciaA = with(dados, tapply(resp, F1, var))Desvio padrão:
DesvioA = with(dados, tapply(resp, F1, sd))Coeficiente de variação:
CVA = DesvioA / MediaA * 100Criando uma tabela com os resultados das medidas descritivas por inoculação:
DescA = cbind(MediaA, VarianciaA, DesvioA, CVA)
DescA MediaA VarianciaA DesvioA CVA
IN 185.45 8229.208 90.71498 48.91614
NI 151.25 14582.724 120.75895 79.84063Analisando as medidas descritivas por inoculação, podemos perceber que a média do número de nódulos das plantas de soja inoculadas, é maior que a média do número de nódulos das plantas não inoculadas, apresentando também uma menor variância e menor desvio padrão, e consequentemente, um menor coeficiente de variação.
Medidas descritivas por época de aplicação:
Média:
MediaB = with(dados, tapply(resp, F2, mean))Variância:
VarianciaB = with(dados, tapply(resp, F2, var))Desvio padrão:
DesvioB = with(dados, tapply(resp, F2, sd))Coeficiente de variação:
CVB = DesvioB / MediaB * 100Criando uma tabela com os resultados das medidas descritivas por época de aplicação:
DescB = cbind(MediaB, VarianciaB, DesvioB, CVB)
DescB MediaB VarianciaB DesvioB CVB
Plantio 221.7 8287.344 91.03485 41.06218
R1+15 152.4 10686.933 103.37762 67.83309
V1+15 101.3 2559.122 50.58777 49.93857
V3+15 198.0 18507.556 136.04248 68.70832Analisando as medidas descritivas por época de aplicação podemos perceber que na época de plantio e V3 as plantas de soja apresentam as maiores médias de nódulos. Já a variância é menor da época V1. E o desvio padrão e o coeficiente de variação são maiores nas épocas V3 e R1.
Análise descritiva por gráficos: criar boxplot para comparar os fatores de inoculação:
par(bty='l', mai=c(1, 1, .2, .2))
par(cex=0.7)
caixas=with(dados, car::Boxplot(resp ~ F1, vertical=T,las=1, col='Lightgreen',
xlab=X, ylab=Y))
mediab=with(dados,tapply(resp, F1, mean))
points(mediab, pch='+', cex=1.5, col='red')
title("Comparação dos Fatores de Inoculação")
Analisando o box plot que compara o número de nódulos por inoculação, podemos perceber graficamente que a média do número de nódulos das plantas de soja inoculadas, é maior que a média do número de nódulos das plantas não inoculadas. E por apresentar uma menor amplitude dos dados podemos notar uma menor variação dos valores nas plantas inoculadas.
Criar boxplot para comparar os fatores de época:
par(bty='l', mai=c(1, 1, .2, .2))
par(cex=0.7)
caixas=with(dados, car::Boxplot(resp ~ F2, vertical=T,las=1, col='Lightblue',
xlab=X, ylab=Y))
mediab=with(dados,tapply(resp, F2, mean))
points(mediab, pch='+', cex=1.5, col='red')
title("Comparação dos Fatores de Época de Aplicação")
Analisando graficamente o boxplot por época de aplicação podemos perceber que, pela amplitude das caixas, na época V1 tem-se a menor variação dos dados, e na época de plantio e V3 as plantas de soja apresentam as maiores médias de nódulos.
Criar box plot para comparar os dois fatores combinados:
par(bty='l', mai=c(1, 1, .2, .2))
par(cex=0.7)
caixas=with(dados, car::Boxplot(resp ~ F1*F2, vertical=T,las=1, col='Lightpink',
xlab=X, ylab=Y))
title("Comparação dos dois fatores combinados")
Fazendo análise gráfica, percebemos que nos tratamentos “IN.R1”, “NI.R1”, “IN.V1”, “NI.V1” e “IN.V3”, a variação dos dados é menor quando comparado aos tratamentos “IN.Plantio”, “NI.Plantio” e NI.V3”. Analizando as medianas, notamos que são maiores nos tratamentos “IN.Plantio” e NI.V3”.
Criar um histograma para todos os dados:
hist(dados$resp, col = "lightgreen", col.axis = "blue", col.lab = "red",
xlab = "Número de nódulos", ylab="Frequência",
main = "Histograma da frequência de número de nódulos")
Analisando graficamente, os dados não aparentam ter distribuição normal.
Fazer o teste de normalidade de todos os dados testando as hipóteses: H0: Os dados seguem uma distribuição normal. Ha: Os dados não seguem uma distribuição normal.
shapiro.test(dados$resp)
Shapiro-Wilk normality test
data: dados$resp
W = 0.93251, p-value = 0.01946O valor-p obtido foi 0,01946, sendo menor que 0,05, dessa forma, a hipótese nula é rejeitada, logo, há evidências que comprovem que os dados não tem distribuição normal para o nível de significância de 5%.
Fazer o teste de homogeneidade de variâncias dos dados por tratamentos H0: As variâncias são homogêneas; Ha: As variâncias não são homogêneas.
bartlett.test(resp ~ Trat, data = dados)
Bartlett test of homogeneity of variances
data: resp by Trat
Bartlett's K-squared = 9.8754, df = 7, p-value = 0.1957O valor-p obtido foi 0,1957, sendo maior que 0,05, dessa forma, a hipótese nula não deve ser rejeitada, havendo evidências que comprovem que os dados apresentam homogeneidade de variâncias dos dados por tratamentos para o nível de significância de 5%.
Análise de variância Hipóteses do fator 1 (inoculação): H0: o número de nódulos será igual para as duas inoculações; Ha: o número de nódulos se difere entre as inoculações.
Hipóteses do fator 2 (épocas): H0: o número de nódulos será igual para todas as épocas; Ha: o número de nódulos se difere em pelo menos uma época.
Hipótes da interação entre os fatores: H0: Todas as combinações entre os níveis do fator 1 e do fator 2 têm o mesmo efeito; Ha: Pelo menos duas combinações entre os níveis do fator 1 e do fator 2 têm efeitos diferentes.
mod = with(dados, aov(resp~F1*F2))
anova(mod)Analysis of Variance Table
Response: resp
Df Sum Sq Mean Sq F value Pr(>F)
F1 1 11696 11696 2.7579 0.106542
F2 3 84754 28251 6.6615 0.001272 **
F1:F2 3 212960 70987 16.7382 1.02e-06 ***
Residuals 32 135712 4241
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para o Fator 1 (F1), o resultado não é significativo, indicando que estatisticamente o número de nódulos será igual na soja inoculada com Rhizobium e na soja não inoculada;
Analisando a Fator 2 (F2), percebemos que o resultado foi significativo para 0,01, indicando que, para esse nível de significância, o número de nódulos se difere em pelo menos uma época de aplicação.
Fazendo a análise da interação entre os fatores, obtivemos um resultado muito significativo, em que concluimos que, com a significância de 0,001, pelo menos uma combinação entre os níveis do fator 1 e do fator 2 têm efeitos diferentes das demais.
Extrair do modelo os valores ajustados e criar um objeto para armazenar os valores ajustados (teóricos):
ajustados = fitted(mod) Extrair do modelo os residuos e criar um objeto para armazenar os residuos:
residuos = residuals(mod)
summary(residuos) Min. 1st Qu. Median Mean 3rd Qu. Max.
-109.80 -37.85 -13.80 0.00 30.25 164.00 Residuos padronizados:
residuosp = rstandard(mod)Grafico dos residuos padronizados (eixo Y) e valores ajustados (eixo X):
plot(ajustados,
residuosp,
xlab= "Valores ajustados",
ylab= "Resíduos padronizados",
ylim=c(-3.5, 3.5),
pch = 16,
col="red",
main = "Resíduos padronizados X Valores ajustados")
abline(h =c(-3, 3), lty = 2, col = "blue")
Com a análise do gráfico, é possível observar que não possue pontos com valores discrepantes que ultrapassem a linha tracejada nos valore 3 e -3. Porém, podemos perceber um determinado padrão entre os valores.
Gráfico para verificar a Normalidade dos resíduos:
car:: qqPlot(residuos, col="red", pch = 16)
[1] 34 21Graficamente os dados apresentam distribuição normal, todos os pontos estão dentro do envelope, apenas dois pontos estão mais a margem.
Teste de normalidade dos residuos: Shapiro Wilk H0: os residuos tem distribuição normal; Ha: os resíduos não tem distribuição normal.
shapiro.test(residuos)
Shapiro-Wilk normality test
data: residuos
W = 0.96809, p-value = 0.3125Como o valor-p encontrado é 0,3125, sendo maior que 0,05, então não rejeitamos a hipótese nula, indicando que os resíduos tem distribuição normal.
Comparando com o resultado encontrado no teste de Normalidade de todos os dados, em que o valor-p obtido foi 0,01946, rejeitando a hipótese nula, encontramos uma diferença ao comparar com a normalidade dos resíduos, visto que, os resíduos apresentaram distribuição normal.
Teste de homogeneidade de variâncias dos resíduos por tratamentos: H0: As variâncias são homogêneas; Ha: As variâncias não são homogêneas.
O usando o teste de Bartlett:
with(dados, bartlett.test(mod$residuals~Trat))
Bartlett test of homogeneity of variances
data: mod$residuals by Trat
Bartlett's K-squared = 9.8754, df = 7, p-value = 0.1957O valor é igual ao teste de homogeneidade de variâncias, onde o valor-p obtido foi 0,1957, sendo maior que 0,05, dessa forma, a hipótese nula não deve ser rejeitada, havendo evidências que comprovem que os dados apresentam homogeneidade de variâncias dos resíduos por tratamentos.
Procedendo com os desdobramentos adequados: Observeando primeiramente as interações:
Interação entre o fator 1:
with(dados, interaction.plot(F2, F1, resp, las=1, col=1:6, bty='l',
xlab='', ylab='CBM', trace.label="FATOR1"))
Nesse caso existe uma interação devido a diferença na direção da resposta.
Interação entre o fator 2
with(dados, interaction.plot(F1, F2, resp, las=1, col=1:6, bty='l',
xlab='', ylab='CBM', trace.label="FATOR2"))
Na época V3 existe uma interação devido a diferença na direção das demais respostas. Já nas épocas plantio, R1 e V1, existe uma interação devido à diferença na grandeza da resposta.
Teste de Comparações:
library(ExpDes.pt)
Attaching package: 'ExpDes.pt'The following object is masked from 'package:MASS':
ginvThe following objects are masked from 'package:agricolae':
lastC, order.group, tapply.statwith(dados,fat2.dic(F1,F2,resp, mcomp="tukey"))------------------------------------------------------------------------
Legenda:
FATOR 1: F1
FATOR 2: F2
------------------------------------------------------------------------
Quadro da analise de variancia
------------------------------------------------------------------------
GL SQ QM Fc Pr>Fc
F1 1 11696 2 2.7579 0.106542
F2 3 84754 3 6.6615 0.001272
F1*F2 3 212960 5 16.7382 0.000001
Residuo 32 135712 4
Total 39 445123 1
------------------------------------------------------------------------
CV = 38.68 %
------------------------------------------------------------------------
Teste de normalidade dos residuos (Shapiro-Wilk)
valor-p: 0.3125183
De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.
------------------------------------------------------------------------
Interacao significativa: desdobrando a interacao
------------------------------------------------------------------------
Desdobrando F1 dentro de cada nivel de F2
------------------------------------------------------------------------
------------------------------------------------------------------------
Quadro da analise de variancia
------------------------------------------------------------------------
GL SQ QM Fc Pr.Fc
F2 3 84754.5 28251.50 6.6615 0.0013
F1:F2 Plantio 1 26112.1 26112.10 6.1571 0.0185
F1:F2 R1+15 1 70896.4 70896.40 16.7169 3e-04
F1:F2 V1+15 1 15288.1 15288.10 3.6048 0.0667
F1:F2 V3+15 1 112360.0 112360.00 26.4938 0
Residuo 32 135712.0 4241.00
Total 39 445123.1 11413.41
------------------------------------------------------------------------
F1 dentro do nivel Plantio de F2
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a 1 272.8
b 2 170.6
------------------------------------------------------------------------
F1 dentro do nivel R1+15 de F2
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a 1 236.6
b 2 68.2
------------------------------------------------------------------------
F1 dentro do nivel V1+15 de F2
De acordo com o teste F, as medias desse fator sao estatisticamente iguais.
------------------------------------------------------------------------
Niveis Medias
1 1 140.4
2 2 62.2
------------------------------------------------------------------------
F1 dentro do nivel V3+15 de F2
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a 2 304
b 1 92
------------------------------------------------------------------------
Desdobrando F2 dentro de cada nivel de F1
------------------------------------------------------------------------
------------------------------------------------------------------------
Quadro da analise de variancia
------------------------------------------------------------------------
GL SQ QM Fc Pr.Fc
F1 1 11696.4 11696.40 2.7579 0.1065
F2:F1 IN 3 105043.8 35014.58 8.2562 3e-04
F2:F1 NI 3 192671.0 64223.65 15.1435 0
Residuo 32 135712.0 4241.00
Total 39 445123.1 11413.41
------------------------------------------------------------------------
F2 dentro do nivel IN de F1
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a 1 272.8
ab 2 236.6
bc 3 140.4
c 4 92
------------------------------------------------------------------------
F2 dentro do nivel NI de F1
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a 4 304
b 1 170.6
b 2 68.2
b 3 62.2
------------------------------------------------------------------------Analisando o quadro da ánalise de variância para o desdobramento do fator 1 dentro dos níveis de fator 2 (épocas) podemos observar que apenas a interação de F1 com a época de aplicação V1+15 apresentou um valor de F maior que 0,05, não sendo estatísticamente significativa.
Ao analisar o teste de Tukey para cada desdobramento individual de F1 dentro dos níveis de F2, observamos que o tratamento inoculado apresentou as maiores médias para três épocas e, na época V3+15 a média para o tratamento não inculado foi maior. Sendo assim, a combinação entre a época V3+15 com a não inoculação apresentou o meior número médio de nódulos.
O desdobramento das épocas dentro do fator inoculação: Dentro do tratamento com sementes inoculadas as aplicações nas épocas “Plantio” e “R1+15” apresentaram os maiores números médios de nódulos e não diferiram estatisticamente entre si. Para o tratamento com sementes não inoculadas apenas o número médio de nódulos para a aplicação na época “V3+15” foi estatisticamente diferente, sendo a maior média dentre as quatro épocas.