Análise de Controle de Plantas Daninhas

Imagem: Ruthson Zimmerman

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

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)
herbicideresponseSEdflower.CLupper.CL.group
4A0.97521420.0079316220.95220530.9872942a
5B0.96356450.0101101220.93570850.9796141a
3C0.91453020.0161835220.87444710.9426558b
2E0.87092710.0197498220.82416470.9066617b
1D0.54164460.0299928220.47911770.6028868c
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
Pesquisador Associado

Trabalho para diminuir o impacto de plantas daninhas em ecossistemas. Sou professor, cientista de plantas daninhas e um entusiasta cientista de dados.