1 Definições da sessão

#-----------------------------------------------------------------------
#                                            Prof. Dr. Walmes M. Zeviani
#                                leg.ufpr.br/~walmes · github.com/walmes
#                                        walmes@ufpr.br · @walmeszeviani
#                      Laboratory of Statistics and Geoinformation (LEG)
#                Department of Statistics · Federal University of Paraná
#                                       2020-jun-11 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
# Pacotes.

library(lme4)
library(lmerTest)
library(emmeans)

library(tidyverse)
library(readxl)

as_letters <- function(x) {
    xs <- strsplit(trimws(as.character(x)), "")
    xs <- lapply(xs, as.integer)
    m <- max(unlist(xs))
    sapply(xs, function(i) paste(rev(letters[m - i + 1]), collapse = ""))
}

2 Dados ao nível de inseto

#-----------------------------------------------------------------------
# Importação.

tb <- gdata::read.xls("DADOS 3 Experimento.xlsx", sheet = 1, skip = 1,
                      dec = ",") %>%
    as_tibble()
# tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 2, skip = 1)
str(tb)
## tibble [180 × 7] (S3: tbl_df/tbl/data.frame)
##  $ ensaio    : int [1:180] 1 1 1 1 1 1 1 1 1 1 ...
##  $ fonte     : int [1:180] 1 1 1 1 1 1 1 1 1 1 ...
##  $ tratamento: chr [1:180] "NAC" "NAC" "NAC" "NAC" ...
##  $ inseto    : int [1:180] 1 2 3 4 5 6 7 8 9 10 ...
##  $ qpcr      : num [1:180] 23.7 24.4 27.8 27.4 29.4 ...
##  $ teste     : int [1:180] 1 1 1 2 2 2 3 3 3 4 ...
##  $ morte     : int [1:180] 0 0 1 0 0 0 0 0 0 0 ...
tb <- tb %>%
    mutate_at(c("ensaio", "tratamento", "fonte", "teste"), "factor")

ggplot(data = tb,
       mapping = aes(x = teste, y = morte, color = tratamento)) +
    facet_wrap(facets = ~fonte, scale = "free_x") +
    geom_jitter(height = 0.1, width = 0)

ggplot(data = tb, #filter(tb, qpcr > 10000),
       mapping = aes(x = teste, y = qpcr, color = tratamento)) +
    facet_wrap(facets = ~fonte, scale = "free_x") +
    geom_point()

#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de qPCR.

# Dicotomizar.
tb$qpcr <- ifelse(tb$qpcr <= 30, 1, 0)

# Ajuste do modelo de efeitos aleatórios.
m0 <- glmer(qpcr ~ (1 | fonte/teste) + tratamento,
            data = tb,
            family = binomial)

# Estimativas dos componentes de variância.
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: qpcr ~ (1 | fonte/teste) + tratamento
##    Data: tb
## 
##      AIC      BIC   logLik deviance df.resid 
##    202.5    215.2    -97.2    194.5      176 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2520 -0.5235 -0.3708  0.4640  2.2812 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  teste:fonte (Intercept) 2.1430   1.4639  
##  fonte       (Intercept) 0.8158   0.9032  
## Number of obs: 180, groups:  teste:fonte, 61; fonte, 3
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    -1.9549     0.7110  -2.750  0.00597 **
## tratamentoNAC   1.5493     0.6084   2.547  0.01087 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tratamntNAC -0.515
# Médias com intervalo de confiança.
emm <- emmeans(m0, specs = ~tratamento, type = "response") %>%
    multcomp::cld()
emm
##  tratamento  prob     SE  df asymp.LCL asymp.UCL .group
##  CONTROLE   0.124 0.0772 Inf    0.0339     0.363  1    
##  NAC        0.400 0.1574 Inf    0.1556     0.707   2   
## 
## 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
emm <- emm %>%
    as.data.frame() %>%
    mutate(cld = as_letters(.group))

ggplot(data = emm,
       mapping = aes(x = prob, y = tratamento)) +
    geom_point() +
    geom_label(mapping = aes(label = sprintf("%0.2f %s", prob, cld)),
               vjust = -0.3) +
    geom_errorbarh(mapping = aes(xmin = asymp.LCL, xmax = asymp.UCL),
                   height = 0.1) +
    labs(x = "Resultado do qPCR",
         y = "Tratamento")

#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de mortalidade.

# Ajuste do modelo de efeitos aleatórios.
m0 <- glmer(morte ~ (1 | fonte/teste) + tratamento,
            family = binomial,
            data = tb)
## boundary (singular) fit: see ?isSingular
# Estimativas dos componentes de variância.
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: morte ~ (1 | fonte/teste) + tratamento
##    Data: tb
## 
##      AIC      BIC   logLik deviance df.resid 
##    133.8    146.6    -62.9    125.8      176 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0164 -0.2328 -0.1099 -0.1099  2.3549 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  teste:fonte (Intercept) 4.681    2.164   
##  fonte       (Intercept) 0.000    0.000   
## Number of obs: 180, groups:  teste:fonte, 61; fonte, 3
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    -2.1936     0.8011  -2.738  0.00618 **
## tratamentoNAC  -2.0551     0.8876  -2.315  0.02060 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tratamntNAC -0.231
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# Médias com intervalo de confiança.
emm <- emmeans(m0, specs = ~tratamento, type = "response") %>%
    multcomp::cld()
emm
##  tratamento   prob     SE  df asymp.LCL asymp.UCL .group
##  NAC        0.0141 0.0146 Inf   0.00182     0.101  1    
##  CONTROLE   0.1003 0.0723 Inf   0.02267     0.349   2   
## 
## 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
emm <- emm %>%
    as.data.frame() %>%
    mutate(cld = as_letters(.group))

ggplot(data = emm,
       mapping = aes(x = prob, y = tratamento)) +
    geom_point() +
    geom_label(mapping = aes(label = sprintf("%0.2f %s", prob, cld)),
               vjust = -0.3) +
    geom_errorbarh(mapping = aes(xmin = asymp.LCL, xmax = asymp.UCL),
                   height = 0.1) +
    labs(x = "Mortalidade",
         y = "Tratamento")

2.1 Dados ao nível de planta teste

#-----------------------------------------------------------------------
# Importação.

tb <- read_excel("DADOS 3 Experimento.xlsx", sheet = 2, skip = 1)
str(tb)
## tibble [180 × 6] (S3: tbl_df/tbl/data.frame)
##  $ ensaio    : num [1:180] 1 1 1 1 1 1 1 1 1 1 ...
##  $ fonte     : num [1:180] 1 1 1 1 1 1 1 1 1 1 ...
##  $ tratamento: chr [1:180] "NAC" "NAC" "NAC" "NAC" ...
##  $ inseto    : num [1:180] 1 2 3 4 5 6 7 8 9 10 ...
##  $ teste     : num [1:180] 1 1 1 2 2 2 3 3 3 4 ...
##  $ pcr       : num [1:180] 0 0 0 0 0 0 0 0 0 0 ...
# DANGER: por alguma razão os dados estão ao nível de inseto mas
# duplicados. Precisa criar a variável ensaio.
# xtabs(~fonte, data = tb)

tb <- tb %>%
    select(-inseto) %>%
    group_by(ensaio, tratamento, fonte, teste) %>%
    slice(1) %>%
    ungroup()

tb <- tb %>%
    mutate_at(c("ensaio", "tratamento", "fonte"), "factor")
str(tb)
## tibble [61 × 5] (S3: tbl_df/tbl/data.frame)
##  $ ensaio    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ fonte     : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ tratamento: Factor w/ 2 levels "CONTROLE","NAC": 1 1 1 1 1 1 1 1 1 1 ...
##  $ teste     : num [1:61] 11 12 13 14 15 16 17 18 19 20 ...
##  $ pcr       : num [1:61] 0 0 1 0 0 0 0 0 0 0 ...
ggplot(data = tb,
       mapping = aes(x = teste, y = pcr, color = tratamento)) +
    facet_wrap(facets = ~fonte, scale = "free_x") +
    geom_point()

# DANGER: um único resultado positivo no primeiro ensaio. Isso pode
# gerar problemas de estimação sem contar que a suposição de variância
# igual pro efeito aleatório de fonte e teste nos ensaios talvez não ser
# realista.

#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de PCR.

# Ajuste do modelo de efeitos aleatórios.
m0 <- glmer(pcr ~ (1 | fonte) + tratamento,
            family = binomial,
            data = tb)

# Estimativas dos componentes de variância.
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: pcr ~ (1 | fonte) + tratamento
##    Data: tb
## 
##      AIC      BIC   logLik deviance df.resid 
##     44.1     50.4    -19.1     38.1       58 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.4613 -0.3792 -0.3011 -0.1981  4.0079 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  fonte  (Intercept) 0.3379   0.5813  
## Number of obs: 61, groups:  fonte, 3
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    -1.9925     0.7078  -2.815  0.00488 **
## tratamentoNAC  -0.8371     0.9214  -0.908  0.36362   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tratamntNAC -0.424
# Médias com intervalo de confiança.
emm <- emmeans(m0, specs = ~tratamento, type = "response") %>%
    multcomp::cld()
emm
##  tratamento   prob     SE  df asymp.LCL asymp.UCL .group
##  NAC        0.0557 0.0470 Inf    0.0102     0.253  1    
##  CONTROLE   0.1200 0.0747 Inf    0.0329     0.353  1    
## 
## 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
emm <- emm %>%
    as.data.frame() %>%
    mutate(cld = as_letters(.group))

ggplot(data = emm,
       mapping = aes(x = prob, y = tratamento)) +
    geom_point() +
    geom_label(mapping = aes(label = sprintf("%0.2f %s", prob, cld)),
               vjust = -0.3) +
    geom_errorbarh(mapping = aes(xmin = asymp.LCL, xmax = asymp.UCL),
                   height = 0.1) +
    labs(x = "Resultado do PCR",
         y = "Tratamento")