Análise de Controle de Plantas Daninhas
O controle de plantas daninhas com herbicidas é um experimento comum realizado na pesquisa de plantas daninhas. Geralmente, os experimentos são organizados em delineamento em blocos casualizados com vários tratamentos (herbicidas). Dias após o tratamento (DAT), uma nota de controle em percentagem (%) é dada às plantas daninhas. A percentagem no controle de plantas daninhas varia de 0% (plantas sem injúrias) e 100% (plantas completamente controladas). Portanto, antes de iniciar os experimentos, o pesquisador sabe que os valores de controle variam entre 0 a 100%.
Esta é uma análise estatística que não segue uma distribuição normal (gaussiana). A distribuição normal é uma premissa da ANOVA. Usaremos a distribuição beta (escala logit), que é uma família de distribuições de probabilidade contínuas definidas no intervalo 0 e 1. Para isso, as notas de controle de plantas daninhas devem ser alteradas para o intervalo 0 e 1. Por exemplo, se o controle de buva for 90%, ele deverá ser listado como 0,90 no seu arquivo de planilha. Além disso, a distribuição beta não aceita valores 0 e 1. Sendo assim, o controle de plantas daninhas de 0% e 100% deve ser alterado para 0,001 e 0,999, respectivamente.
Vamos começar
Baixe o programa R de acordo com seu sistema operacional: https://www.r-project.org
Baixe o programa RStudio de acordo com seu sistema operacional: https://rstudio.com
Crie um arquivo no R
- Clique em File -> New project… -> New directory -> New Project -> Salve o diretório de trabalho R onde quiser no seu laptop.
O arquivo salvo terá um projeto do R e você deve copiar e colar todos os seus dados (ex.,arquivo excel) no mesmo diretório de trabalho criado.
- Clique em Arquivo -> Novo arquivo -> script R (.R)
ou
- Clique em Arquivo -> Novo arquivo -> remarcação R … (.Rmd)
Embora eu goste de trabalhar no markdown R, o script R é mais fácil para iniciantes no Rstudio.
Instale os pacotes
Copie e cole os seguintes códigos no seu aquivo .R ou .Rmd.
install.packages(tidyverse)
install.packages(glmmTMB)
install.packages(lme4)
install.packages(lmerTest)
install.packages(emmeans)
install.packages(RCurl)
install.packages(car)
install.packages(kableExtra)
Execute todos os códigos abaixo clicando na opção Run no canto superior direito do seu aquivo .R ou .Rmd. Esses códigos instalam os pacotes necessários para analisar experimentos de controle de plantas daninhas (%). Depois de instalar esses pacotes, você não precisará instalá-los toda vez que abrir o Rstudio, a menos que você atualize o R ou o Rstudio.
Dados
library(tidyverse)
library(RCurl)
Os dados utilizados para o exercício são de um artigo científico de controle (%) de biótipo de Amaranthus tuberculatus (vou chamar de caruru aqui) com herbicidas em pós-emergência. Os nomes dos herbicidas foram alterados para A, B, C, D e E.
Copie cole os seguintes códigos abaixo e execute no R.
df_path <- getURL("https://raw.githubusercontent.com/maxwelco/General_codes/master/anova-beta/data/control.csv")
dado <- read_csv(df_path)
Após executar os codigos acima, os dados deve estar armazenado no nome dado.
dado
## # A tibble: 30 x 4
## year herbicide rep control
## <dbl> <chr> <dbl> <dbl>
## 1 2013 A 1 0.99
## 2 2013 A 2 0.97
## 3 2013 A 3 0.96
## 4 2013 B 1 0.99
## 5 2013 B 2 0.99
## 6 2013 B 3 0.97
## 7 2013 C 1 0.95
## 8 2013 C 2 0.9
## 9 2013 C 3 0.95
## 10 2013 D 1 0.5
## # … with 20 more rows
O aquivo dado possui o controle (control) de caruru com 5 trataments (herbicidas A, B, C, D, E), 2 anos (year, 2013 e 2014) e 3 blocos (rep).
Homogeneidade de variância
Testa a hipótese nula de que a homogeneidade de variância entre cada ano são iguais.
library(car)
Função Levene
leveneTest(control ~ herbicide, data = dado)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 0.5109 0.7282
## 25
Os resultados mostram o valor P-value = 0.7282, a hipótese nula é rejeitada, o que significa que as variações em cada ano são iguais.
Modelo
- Execute os pacotes abaixo
library(glmmTMB)
library(lme4)
library(lmerTest)
Generalized Linear Mixed Models using Template Model Builder de Mollie Brooks
Modelo misto
Familia beta com logit
O modelo (model) tem control (controle) como variável resposta, herbicida como efeitos fixos e year (ano) (junto com rep) como efeitos aleatórios.
Função glmmTMB
model <- glmmTMB(control ~ herbicide + (1|year/rep), beta_family(link = "logit"), data=dado)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
Anova
- Execute o pacote abaixo
Função Anova.glmmTMB
glmmTMB:::Anova.glmmTMB(model)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: control
## Chisq Df Pr(>Chisq)
## herbicide 212.35 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
O P-valor < 2.2e-16; existe evidência que é baixa a probabilidade de que uma diferença de alguma magnitude entre os herbicidas testados ocorre apenas por acaso, assumindo que a hipótese nula é verdadeira.
Avaliando o controle de caruru com herbicidas
- Execute o pacote abaixo
library(emmeans)
Gráficos de interação para médias marginais estimadas (emmip)
type=“response”, herbicida e intervalos de confiança de volta da escala logit
coord_flip(), inverte os eixos x e y para melhor visualização
Função emmip
emmip(model, ~herbicide, type="response") + coord_flip()
Esta figura ajuda a visualizar o efeito do herbicida no controle de caruru. É mais útil em experiemnts de fatorial.
Médias marginais estimadas (médias dos quadrados mínimos)
type=“response”, traz herbicida e intervalos de volta da escala logit
cont=“pairwise”, comparações entre cada tratamento com herbicida
adjust=“none”, diferença menos significativa de fisher (pode usar “tukey”)
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).
Função emmeans
lsmeans <- emmeans(model, ~ herbicide, cont="pairwise", adjust="none", type="response", alpha=0.05)
lsmeans
## $emmeans
## herbicide response SE df lower.CL upper.CL
## A 0.975 0.00793 22 0.952 0.987
## B 0.964 0.01011 22 0.936 0.980
## C 0.915 0.01618 22 0.874 0.943
## D 0.542 0.02999 22 0.479 0.603
## E 0.871 0.01975 22 0.824 0.907
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df t.ratio p.value
## A / B 1.488 0.6268 22 0.943 0.3559
## A / C 3.677 1.3992 22 3.422 0.0024
## A / D 33.295 11.6335 22 10.033 <.0001
## A / E 5.831 2.1417 22 4.801 0.0001
## B / C 2.472 0.8619 22 2.595 0.0165
## B / D 22.379 6.9838 22 9.960 <.0001
## B / E 3.919 1.3066 22 4.097 0.0005
## C / D 9.055 2.1694 22 9.196 <.0001
## C / E 1.586 0.4272 22 1.711 0.1011
## D / E 0.175 0.0373 22 -8.175 <.0001
##
## Tests are performed on the log odds ratio scale
Visualizar as médias de controle do caruru
comparisons=TRUE, comparações entre tratamentos com setas vermelhas
type=“response”, traz herbicida e intervalos de volta da escala logit
alpha=0.05, nível de significância a ser usado na construção de setas de comparação
adjust=“none”, diferença menos significativa de fisher (pode usar “tukey”)
Função plot
plot(lsmeans, ~ herbicide, comparisons =TRUE, type="response", alpha=0.05, adjust="none")
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. Veja os contrastes de lsmeans para verificar novamente.
Extrair e exibir informações em todas as comparações pareadas de médias marginais estimadas
alpha=0.05, valor numérico que fornece o nível de significância para as comparações
Letters=letters, exibir tratamentos de agrupamento de letras de acordo com comparações pareadas
adjust=“none”, diferença menos significativa de fisher (pode usar “tukey”)
reversed = TRUE, a ordem de uso das letras é invertida
Função CLD
cld <- cld(lsmeans$emmeans, alpha=0.05, Letters=letters, adjust="none", reversed = TRUE)
cld
## herbicide response SE df lower.CL upper.CL .group
## A 0.975 0.00793 22 0.952 0.987 a
## B 0.964 0.01011 22 0.936 0.980 a
## C 0.915 0.01618 22 0.874 0.943 b
## E 0.871 0.01975 22 0.824 0.907 b
## D 0.542 0.02999 22 0.479 0.603 c
##
## 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.05
Aviso
"A função e os métodos do CLD estão obsoletos. Os agrupamentos de letras (CLD) incentivam uma interpretação enganosa do teste de significância por meio de agrupamento visual cujas comparações têm P> alfa como se fossem iguais. No entanto, não provar que duas médias são diferentes não significa que as médias são iguais. Além disso, os CLDs fazem uma distinção difícil entre os valores de P quase iguais a alfa, mas em lados opostos."
Figura
Gráfico de pontos com intervalo de confiança
Como a função CLD está obsoleta, eu recomendo usar lsmeans para criar uma figura.
nd <- as.data.frame(lsmeans$emmeans)
Use nd para plotar uma figura.
ggplot function
ggplot(nd, aes(x=reorder(herbicide,response), y=response*100, color=herbicide)) +
geom_point(size=4) + ylim(0,100) + coord_flip() +
scale_color_manual(values=c("red", "blue", "green", "orange", "purple")) +
theme_bw() + labs(y="Controle de caruru (%)", x="Herbicidas") +
geom_linerange(aes(ymin = lower.CL*100, ymax = upper.CL*100), size=1.5) +
theme(axis.title = element_text(size=16),
axis.text = element_text(size=15),
legend.position = "none")
As médias de controle de caruru (%) e o intervalo de confiança de cada herbicida são exibidos em uma figura do gráfico ggplot.
Gráfico de barras com agrupamento de letras
Se você gosta de apresentar os resultados em gráfico de barras com letras. Execute os codigos abaixo.
ggplot(cld, aes(x=reorder(herbicide,response),
y=response*100, fill=herbicide, label=.group)) +
geom_bar(stat="identity") + ylim(0,105) +
scale_fill_manual(values=c("red", "blue", "green", "orange", "purple")) +
theme_bw() + labs(y="Waterhemp control (%)", x="Herbicides") +
geom_text(nudge_y = 7, size=8) +
theme(axis.title = element_text(size=16),
axis.text = element_text(size=15),
legend.position = "none") +
coord_flip()
Tabela
- Execute o pacote abaixo
library(kableExtra)
Este é uma tabela muito simples.
Função kable
cld %>%
kbl() %>%
kable_classic_2(full_width = F)
herbicide | response | SE | df | lower.CL | upper.CL | .group | |
---|---|---|---|---|---|---|---|
4 | A | 0.9752142 | 0.0079316 | 22 | 0.9522053 | 0.9872942 | a |
5 | B | 0.9635645 | 0.0101101 | 22 | 0.9357085 | 0.9796141 | a |
3 | C | 0.9145302 | 0.0161835 | 22 | 0.8744471 | 0.9426558 | b |
2 | E | 0.8709271 | 0.0197498 | 22 | 0.8241647 | 0.9066617 | b |
1 | D | 0.5416446 | 0.0299928 | 22 | 0.4791177 | 0.6028868 | c |
Essa postagem é a tradução do post em inglês One-Way ANOVA for Percentage Weed Control Data escrito por Maxwel C Oliveira e Rodrigo Werle no projeto Open Weed Science |