ANOVA
Maxwel Coura Oliveira
1 Importando os pacotes
Vamos importar alguns pacotes até o final dessa aula, porém vamos trabalhar na base tidyverse. O tidyverse é um pacote que contém vários pacotes que auxiliam na manipulação de dados e construção de gráficos. Um exemplo que vamos usar é o pipe (%>%) do pacote magrittr, que faz parte do tidyverse. Assim como o pacote dplyr para manipulaccção de dados.
library(tidyverse)2 Dados
Parte dos resultados de Confirmation and control of HPPD-inhibiting herbicide–resistant waterhemp (Amaranthus tuberculatus) in Nebraska de Oliveira et al. (2017). Nove herbicidas pré-emergentes aplicados em milho. Os estudos foi conduzido em campo em 2013 e 2014.
Vamos usar o pacote RCurl para baixar o pacote via GitHub.
library(RCurl)# Baixando via GitHub
df_path <- url("https://raw.githubusercontent.com/maxwelco/workshop-esalq/main/data/residual.csv")
# Lendo via read_csv
residual <- read_csv(df_path)
residual## # A tibble: 54 x 7
## weed year herbicide trt rep dat control
## <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 waterhemp 2014 SureStart 1 1 34 100
## 2 waterhemp 2014 SureStart 1 2 34 100
## 3 waterhemp 2014 SureStart 1 3 34 95
## 4 waterhemp 2014 Lumaz EZ + Harness 2 1 34 50
## 5 waterhemp 2014 Lumaz EZ + Harness 2 2 34 30
## 6 waterhemp 2014 Lumaz EZ + Harness 2 3 34 70
## 7 waterhemp 2014 Zidua_1 3 1 34 70
## 8 waterhemp 2014 Zidua_1 3 2 34 100
## 9 waterhemp 2014 Zidua_1 3 3 34 80
## 10 waterhemp 2014 Zidua_2 4 1 34 85
## # … with 44 more rows3 Manipulação dos dados
O resultado (control), % de controle de planta daninha (waterhemp) está em uma escala de 0 a 100, precisamos colocar em uma escala entre 0 e 1 (proporção)
residual1 <- residual %>%
mutate(control_prop = control / 100)
residual1## # A tibble: 54 x 8
## weed year herbicide trt rep dat control control_prop
## <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 waterhemp 2014 SureStart 1 1 34 100 1
## 2 waterhemp 2014 SureStart 1 2 34 100 1
## 3 waterhemp 2014 SureStart 1 3 34 95 0.95
## 4 waterhemp 2014 Lumaz EZ + Harness 2 1 34 50 0.5
## 5 waterhemp 2014 Lumaz EZ + Harness 2 2 34 30 0.3
## 6 waterhemp 2014 Lumaz EZ + Harness 2 3 34 70 0.7
## 7 waterhemp 2014 Zidua_1 3 1 34 70 0.7
## 8 waterhemp 2014 Zidua_1 3 2 34 100 1
## 9 waterhemp 2014 Zidua_1 3 3 34 80 0.8
## 10 waterhemp 2014 Zidua_2 4 1 34 85 0.85
## # … with 44 more rowsPronto! Os dados de controle estão na escala de proporção. Porém existem vários 0 e 1. Devemos eliminar o 0 e 1 e mudar para 0.001 e 0.999, respectivamente.
residual2 <- residual1 %>%
mutate(control_prop = case_when(
control_prop == 1.00 ~ 0.999,
control_prop == 0.00 ~ 0.001,
TRUE ~ control_prop
))
residual2## # A tibble: 54 x 8
## weed year herbicide trt rep dat control control_prop
## <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 waterhemp 2014 SureStart 1 1 34 100 0.999
## 2 waterhemp 2014 SureStart 1 2 34 100 0.999
## 3 waterhemp 2014 SureStart 1 3 34 95 0.95
## 4 waterhemp 2014 Lumaz EZ + Harness 2 1 34 50 0.5
## 5 waterhemp 2014 Lumaz EZ + Harness 2 2 34 30 0.3
## 6 waterhemp 2014 Lumaz EZ + Harness 2 3 34 70 0.7
## 7 waterhemp 2014 Zidua_1 3 1 34 70 0.7
## 8 waterhemp 2014 Zidua_1 3 2 34 100 0.999
## 9 waterhemp 2014 Zidua_1 3 3 34 80 0.8
## 10 waterhemp 2014 Zidua_2 4 1 34 85 0.85
## # … with 44 more rowsAgora sim, os dados estão prontos!
Mas antes um pequeno ajuste
residual3 <- residual2 %>%
mutate_if(is_character, as_factor) %>%
mutate_at(c("year", "trt", "rep"), as_factor)
residual3## # A tibble: 54 x 8
## weed year herbicide trt rep dat control control_prop
## <fct> <fct> <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 waterhemp 2014 SureStart 1 1 34 100 0.999
## 2 waterhemp 2014 SureStart 1 2 34 100 0.999
## 3 waterhemp 2014 SureStart 1 3 34 95 0.95
## 4 waterhemp 2014 Lumaz EZ + Harness 2 1 34 50 0.5
## 5 waterhemp 2014 Lumaz EZ + Harness 2 2 34 30 0.3
## 6 waterhemp 2014 Lumaz EZ + Harness 2 3 34 70 0.7
## 7 waterhemp 2014 Zidua_1 3 1 34 70 0.7
## 8 waterhemp 2014 Zidua_1 3 2 34 100 0.999
## 9 waterhemp 2014 Zidua_1 3 3 34 80 0.8
## 10 waterhemp 2014 Zidua_2 4 1 34 85 0.85
## # … with 44 more rowsO que fiz foi mudar as colunas de character (nomes) para fctr (fator) usando a função mutate_if e alguns colunas dbl (números) para fctr.
4 Visualização de dados
Sempre visualize seus dados antes de proceder com a análise estatística. Aqui vamos usar o pacote ggplot do grupo tidyverse. Você já importou o pacote ggplot via tidyverse.
4.1 Gráfico de pontos
Aqui plotei um gráfico bem simples para visualizar a dispersão dos dados de controle de planta daninha em dois anos.
ggplot(residual3, aes(x = herbicide, y = control_prop, color = year)) +
geom_point() +
coord_flip() +
ylim(0, 1)
4.2 Gráfico box-plot
Na figura abaixo adicionei um box-plot para melhorar a visualização da dispersão dos dados.
ggplot(residual3, aes(x = herbicide, y = control_prop, color = year)) +
geom_boxplot() + # veja que adicionei esse código para gerar o box plot
geom_point() +
coord_flip() +
ylim(0, 1)
Aqui mostrei uma pequena parte do que pode ser feito com a parte de visuazlição exploratória de dados. Existe n possibilidades. Com o ggplot é possivel fazer figuras fantásticas! Duvida? Veja isso aqui: TidyTuesday Bookmarks.
5 ANOVA
5.1 Premissas da ANOVA
Normalidade (Distribuição Gausseana)
Homogeneidade de variância
Independência
5.2 Normalidade
Percentagem são dados tecnicamente não normalizados e vamos usar beta distribution.
5.3 Homogeneidade de variâncias
Aqui vamos usar a função leveneTest para verificar a homogeneidade de variâncias entre os anos.
library(car)
leveneTest(control_prop ~ year, data = residual3)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 25.857 5.112e-06 ***
## 52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note que o P-valor é significativo, tecnicamente não poderia juntar os anos. Mas nessa análise vou juntar para facilitar.
- Definir juntar ou não dados de anos/locais distintos usando estatística é sempre uma discussão e não vamos entrar nesse mérito aqui.
5.4 Independência
Relacionado como foi conduzido o estudo e como coletou os dados.
5.5 Modelo
library(glmmTMB)modelo <- glmmTMB(control_prop ~ # variável resposta
herbicide + # variável fixa
(1|rep/year), # efeitos randomizados
family = beta_family(link = "logit"),
REML = FALSE,
data = residual3)5.5.1 Efeitos randomizados
5.6 Anova
glmmTMB:::Anova.glmmTMB(modelo)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: control_prop
## Chisq Df Pr(>Chisq)
## herbicide 61.851 8 2.018e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pronto! Executamos o modelo e a ANOVA. O resultado diz que herbicide é significativo Pr(>Chisq) = 2.018e-10 ***.
Temos que avaliar o desempenho dos herbicidas no controle da planta daninha!
6 Avaliando o controle da planta daninha com herbicidas
Vamos usar o pacote emmeans
library(emmeans)Primeira visualização usando a funccção emmip. Note que no código tem um atributo type="response", que faz a transformação dos dados de volta para a proporção (controle).
emmip(modelo, ~ herbicide, type="response") + coord_flip()
6.1 Médias marginais estimadas (médias dos quadrados mínimos)
A lsmeans forneceram os valores com o controle de plantas daninhas (prop), SE (erro padrão) e intervalos de confiança (lower.CL e upper.CL). Além disso, lsmeans fornece os contrastes aos pares entre tratamentos (herbicidas).
lsmeans <- emmeans(modelo, # modelo
~ herbicide, # efeito fixo
cont="pairwise", # tipo de comparação
type="response", # transforma os dados de volta ao original
alpha=0.05) # nivel
lsmeans## $emmeans
## herbicide response SE df lower.CL upper.CL
## SureStart 0.982 0.0114 42 0.937 0.995
## Lumaz EZ + Harness 0.760 0.0876 42 0.545 0.893
## Zidua_1 0.954 0.0254 42 0.866 0.985
## Zidua_2 0.901 0.0471 42 0.758 0.963
## Zidua_3 0.582 0.1137 42 0.352 0.782
## Anthem ATZ 0.808 0.0752 42 0.613 0.918
## Zidua + Sharpen + AAtrex 0.873 0.0564 42 0.711 0.950
## Verdict + Outlook 0.869 0.0577 42 0.704 0.948
## Corvus + AAtrex 0.840 0.0665 42 0.659 0.935
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df
## SureStart / (Lumaz EZ + Harness) 16.801 8.8571 42
## SureStart / Zidua_1 2.561 1.4099 42
## SureStart / Zidua_2 5.854 3.0977 42
## SureStart / Zidua_3 38.111 20.2834 42
## SureStart / Anthem ATZ 12.587 6.5982 42
## SureStart / (Zidua + Sharpen + AAtrex) 7.721 4.0773 42
## SureStart / (Verdict + Outlook) 8.036 4.2842 42
## SureStart / (Corvus + AAtrex) 10.075 5.3264 42
## (Lumaz EZ + Harness) / Zidua_1 0.152 0.0725 42
## (Lumaz EZ + Harness) / Zidua_2 0.348 0.1443 42
## (Lumaz EZ + Harness) / Zidua_3 2.268 0.7917 42
## (Lumaz EZ + Harness) / Anthem ATZ 0.749 0.2760 42
## (Lumaz EZ + Harness) / (Zidua + Sharpen + AAtrex) 0.460 0.1817 42
## (Lumaz EZ + Harness) / (Verdict + Outlook) 0.478 0.1869 42
## (Lumaz EZ + Harness) / (Corvus + AAtrex) 0.600 0.2279 42
## Zidua_1 / Zidua_2 2.286 1.1233 42
## Zidua_1 / Zidua_3 14.879 7.0578 42
## Zidua_1 / Anthem ATZ 4.914 2.3184 42
## Zidua_1 / (Zidua + Sharpen + AAtrex) 3.014 1.4545 42
## Zidua_1 / (Verdict + Outlook) 3.137 1.5466 42
## Zidua_1 / (Corvus + AAtrex) 3.933 1.9003 42
## Zidua_2 / Zidua_3 6.510 2.6663 42
## Zidua_2 / Anthem ATZ 2.150 0.8965 42
## Zidua_2 / (Zidua + Sharpen + AAtrex) 1.319 0.5715 42
## Zidua_2 / (Verdict + Outlook) 1.373 0.5951 42
## Zidua_2 / (Corvus + AAtrex) 1.721 0.7318 42
## Zidua_3 / Anthem ATZ 0.330 0.1180 42
## Zidua_3 / (Zidua + Sharpen + AAtrex) 0.203 0.0787 42
## Zidua_3 / (Verdict + Outlook) 0.211 0.0808 42
## Zidua_3 / (Corvus + AAtrex) 0.264 0.0980 42
## Anthem ATZ / (Zidua + Sharpen + AAtrex) 0.613 0.2433 42
## Anthem ATZ / (Verdict + Outlook) 0.638 0.2527 42
## Anthem ATZ / (Corvus + AAtrex) 0.800 0.3083 42
## (Zidua + Sharpen + AAtrex) / (Verdict + Outlook) 1.041 0.4334 42
## (Zidua + Sharpen + AAtrex) / (Corvus + AAtrex) 1.305 0.5319 42
## (Verdict + Outlook) / (Corvus + AAtrex) 1.254 0.5070 42
## t.ratio p.value
## 5.352 0.0001
## 1.709 0.7378
## 3.340 0.0419
## 6.840 <.0001
## 4.831 0.0006
## 3.870 0.0101
## 3.909 0.0091
## 4.369 0.0024
## -3.955 0.0080
## -2.546 0.2404
## 2.347 0.3390
## -0.784 0.9968
## -1.967 0.5736
## -1.887 0.6260
## -1.346 0.9112
## 1.682 0.7537
## 5.692 <.0001
## 3.375 0.0383
## 2.287 0.3728
## 2.319 0.3541
## 2.835 0.1356
## 4.574 0.0013
## 1.836 0.6591
## 0.639 0.9992
## 0.731 0.9980
## 1.277 0.9326
## -3.100 0.0750
## -4.109 0.0051
## -4.064 0.0058
## -3.588 0.0220
## -1.232 0.9444
## -1.133 0.9654
## -0.578 0.9996
## 0.096 1.0000
## 0.653 0.9991
## 0.559 0.9997
##
## P value adjustment: tukey method for comparing a family of 9 estimates
## Tests are performed on the log odds ratio scale6.2 Visualizar as médias de controle
plot(lsmeans, # lsmeans
~ herbicide, # efeito fixo
comparisons = TRUE, # quero comparar
type="response", # transforma os dados de volta ao original
alpha=0.05, # nivel
adjust="none") # ajuste, pode ser Tukey por exemplo - None = Fisher LSD
O ponto preto é a média do controle do herbicida e o roxo é o intervalo de confiança. A seta vermelha é uma comparação entre tratamentos; se uma seta vermelha não se sobrepõe, os tratamentos são diferentes.
#@ Separação por médias
Vamos usar o pacote multcomp
library(multcomp)cld <- cld(lsmeans$emmeans, # o uso do $emmeans filtra o resultado de lsmeans para as medias
alpha=0.05, # nivel
Letters=letters, # usar letras
adjust="none", # sem ajuste de medias
reversed = TRUE) # a ---> z
cld## herbicide response SE df lower.CL upper.CL .group
## SureStart 0.982 0.0114 42 0.937 0.995 a
## Zidua_1 0.954 0.0254 42 0.866 0.985 ab
## Zidua_2 0.901 0.0471 42 0.758 0.963 bc
## Zidua + Sharpen + AAtrex 0.873 0.0564 42 0.711 0.950 cd
## Verdict + Outlook 0.869 0.0577 42 0.704 0.948 cd
## Corvus + AAtrex 0.840 0.0665 42 0.659 0.935 cd
## Anthem ATZ 0.808 0.0752 42 0.613 0.918 cd
## Lumaz EZ + Harness 0.760 0.0876 42 0.545 0.893 d
## Zidua_3 0.582 0.1137 42 0.352 0.782 e
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05Pronto! Sua análise está feita. Existe n possibilidades de reportar esses dados. Você pode criar uma tabela ou usar figuras, segue algumas sugestões abaixo.
7 Reportando os dados
ggplot(cld, aes(x=reorder(herbicide, response), y= response * 100, color = herbicide)) +
geom_point(size=4) +
geom_linerange(aes(ymin = lower.CL*100, ymax = upper.CL*100), size=1.5) +
ylim(0,100) +
labs(y="Waterhemp control (%)", x="") +
theme_bw() +
theme(axis.title = element_text(size=16),
axis.text = element_text(size=15),
legend.position = "none") +
coord_flip() 
ggplot(cld, aes(x=herbicide, y=response * 100, fill=herbicide, label=.group)) +
geom_bar(stat="identity") +
ylim(0,110) +
labs(y="Waterhemp control (%)", x="Herbicides") +
geom_text(nudge_y = 7, nudge_x = 0, size = 8) +
theme_bw() +
theme(legend.position = "none")
7.1 Tabela
Vamos usar o pacote kableExtra
library(kableExtra)kable(cld) %>%
kable_styling()| herbicide | response | SE | df | lower.CL | upper.CL | .group | |
|---|---|---|---|---|---|---|---|
| 5 | SureStart | 0.9815083 | 0.0113984 | 42 | 0.9372839 | 0.9947233 | a |
| 2 | Zidua_1 | 0.9539639 | 0.0254327 | 42 | 0.8655886 | 0.9852244 | ab |
| 6 | Zidua_2 | 0.9006609 | 0.0471060 | 42 | 0.7580611 | 0.9632828 | bc |
| 9 | Zidua + Sharpen + AAtrex | 0.8730112 | 0.0564129 | 42 | 0.7111426 | 0.9504879 | cd |
| 8 | Verdict + Outlook | 0.8685050 | 0.0576837 | 42 | 0.7044304 | 0.9481974 | cd |
| 7 | Corvus + AAtrex | 0.8404720 | 0.0665060 | 42 | 0.6594246 | 0.9347929 | cd |
| 4 | Anthem ATZ | 0.8083161 | 0.0751900 | 42 | 0.6129593 | 0.9182229 | cd |
| 3 | Lumaz EZ + Harness | 0.7595689 | 0.0876114 | 42 | 0.5454170 | 0.8926849 | d |
| 1 | Zidua_3 | 0.5820680 | 0.1137236 | 42 | 0.3515661 | 0.7815461 | e |
8 Leitura complementar
Kniss et al. Stats4ag
Bates et al. Fitting Linear Mixed-Effects Models Using lme4