Antes, limpar o “Ambiente Global” (Global Environment):
rm(list=ls(all=TRUE))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, usarando 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:
pesosTratamentos (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çõesMostrar a ordem das rações com função table:
table(rações)Esta na ordem alfabética porque as rações estão definidas como letras do alfabeto por isso não modifica a ordem. Colocar na ordem:
rações = factor(rações, levels = c('A',
'B',
'C',
'D'))Mostrar a ordem das rações em ordem:
table(rações)Criar o banco de dados com a variável resposta e o tratamento:
dados = data.frame(pesos,rações)
dadosAnálise descritiva:
Resumo númerico geral:
summary(dados)Resumo númerico por tratamento
tapply(pesos,rações,summary)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'))Boxplot 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")Com o histograma é possivel observar uma maior frequência nos valores 15 a 30.
Teste de Normalidade de todos os dados: Histograma:
ggplot(dados, aes(x=pesos))+
geom_histogram(bins = 6,
fill = "green4",
col= "magenta2")O histograma não está muito próximo da aparência de um sino, os últimos valores estão altos.
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)Os pontos seguem aproximadamente uma linha reta, porém um valor se encontra fora da linha.
Teste de Shapiro-Wilk para Normalidade: Hipotese nula: os dados de pesos possuem uma distribuição Normal; Hipotese alternativa: os dados de pesos nao possuem uma distribuição Normal.
shapiro.test(pesos)Como 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) 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_aov
summary(mod_aov) #produz a tabela da ANOVAComo 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
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
ajustadosGrá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")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)Os 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)Como 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)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)
TukeyHSD(mod_aov,
ordered = TRUE,
conf.level = 0.95)
res_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 nomesConsiderando 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 cultivar
df.residual(mod_aov) #Graus de liberdade do resíduo
glres = df.residual(mod_aov)
deviance(mod_aov)/glres #Quadrado médio do resíduo
QMRes = deviance(mod_aov)/glres
res_tukey2 = HSD.test(pesos,
rações,
glres,
QMRes,
alpha=0.05)
res_tukey2
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
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.
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 de 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) Tratamentos 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:
tratamentoMostrar a ordem dos tratamentos:
table(tratamento)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)Criar o objeto para representar os blocos. Usar a funcao “factor()”:
bloco = factor(rep(c("B1",
"B2",
"B3",
"B4"), times = 8))Criar objeto bloco e visualizar comprimento:
bloco
length(bloco) Organizar os objetos em um banco de dados (data frame):
dados = data.frame(bloco, tratamento, produção)
dadosEstrutura do banco de dados criado:
str(dados)Resumo numérico do banco de dados:
summary(dados)Análise descritva da variável resposta (produção) Análise por tratamento: Resumo númerico:
tapply(produção, tratamento, summary)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_tratMediana:
mediana_produção_trat = round(tapply(produção, tratamento, median),1)
mediana_produção_tratVariância:
var_produção_trat = round(tapply(produção, tratamento, var),1)
var_produção_tratDesvio padrão:
dp_produção_trat = round(sqrt(var_produção_trat),1)
dp_produção_tratCoeficiente 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"))dados |>
group_by(tratamento) |>
rstatix::get_summary_stats(produção, type = "full")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
mediana_produção_blocos = tapply(produção, bloco, median)
mediana_produção_blocos
var_produção_blocos = tapply(produção, bloco, var)
var_produção_blocosAo 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")O histograma não está com aparência próxima a de um sino, nos penultimos valores é possivel observar 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)")O histograma apresenta aparência próxima a de um sino, mesmo que os penúltimos valores possuem aumento na frequência.
Função “qqplot()” do pacote car:
qqPlot(produção, col="magenta1", pch=19)Os pontos seguem aproximadamente uma linha reta, 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)Como 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)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_aovTabela da ANOVA:
anova(mod_aov)Valor-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
residuosResumo numérico dos resíduos:
summary(residuos)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)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")Os 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)O 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)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
plot(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 TRAT
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_tukey2Representaçã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
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.
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))
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))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)
descA 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)
DescAAnalisando 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)
DescBAnalisando 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)O 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)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 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)Para 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)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)Graficamente 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)Como 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))O 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)
with(dados,fat2.dic(F1,F2,resp, mcomp="tukey"))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.