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 = ""))
}#-----------------------------------------------------------------------
# Importação e inspeção.
# Sheet de 1 a 6.
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 1, skip = 1)
str(tb)
# NOTE: não existe tratamento. Apenas modelar para considerar o efeito
# dos fatores estruturais.
# · Observações ao nível de inseto.
# · Estrutura: 2 ensaios > 6 fontes/ensaio > 10 testes/fonte > 3
# insetos/teste.
# · Duas respostas.
# · Modelo: resposta ~ (1 | ensaio/fonte/teste).
xtabs(~ensaio + fonte, data = tb)
xtabs(~teste + fonte, data = tb)
ggplot(data = tb,
mapping = aes(x = teste, y = morte, color = ensaio)) +
facet_wrap(facets = ~fonte, scale = "free_x") +
geom_jitter(height = 0.1, width = 0)
#--------------------------------------------
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 2, skip = 1)
str(tb)
# NOTE: As plantas fontes receberam tratamento. O resto permanece
# igual. São plantas fontes independentes das do experimento anterior.
# · Observações ao nível de inseto.
# · Estrutura: 2 ensaios > 2 tratamentos/ensaio > 3 fontes/tratamento >
# 10 testes/fonte > 3 insetos/teste.
# · Duas respostas.
# · Modelo: resposta ~ (1 | ensaio/fonte/teste) + tratamento.
# · O tratamento foi casualizado nas fontes.
xtabs(~ensaio + tratamento, data = tb)
xtabs(~tratamento + fonte, data = tb)
xtabs(~teste + fonte, data = tb)
#--------------------------------------------
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 3, skip = 1)
str(tb)
# NOTE: Não existe tratamento.
# · Observações ao nível de planta teste.
# · Estrutura: 2 ensaios > 6 fontes/ensaio > 10 testes/fonte.
# · Uma respostas.
# · Modelo: resposta ~ (1 | ensaio/fonte).
xtabs(~teste + fonte, data = tb)
# DANGER: existe ainda a coluna `inseto` que não deveria e os dados
# estão ao nível de inseto. Será que estão repetidos.
# Verifica se os valores estão duplicados na dimensão do `inseto`.
tb %>%
group_by(fonte, teste) %>%
summarise_at("qpcr", c("mean", "var", "length")) %>%
print(n = Inf)
#--------------------------------------------
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 4, skip = 1)
str(tb)
# NOTE: Tratamento aplicado às plantas fonte.
# · Observações ao nível de planta teste.
# · Estrutura: 2 ensaios > 2 tratamentos/ensaio > 3 fontes/tratamento >
# 10 testes/fonte.
# · Uma respostas.
# · Modelo: resposta ~ (1 | ensaio/fonte) + tratamento.
# · O tratamento foi casualizado nas fontes.
xtabs(~ensaio + tratamento, data = tb)
xtabs(~tratamento + fonte, data = tb)
xtabs(~teste + fonte, data = tb)
# DANGER: existe ainda a coluna `inseto` que não deveria e os dados
# estão ao nível de inseto. Será que estão repetidos.
# Verifica se os valores estão duplicados na dimensão do `inseto`.
tb %>%
group_by(ensaio, tratamento, fonte, teste) %>%
summarise_at("qpcr", c("mean", "var", "length")) %>%
print(n = Inf)
#--------------------------------------------
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 5, skip = 1)
str(tb)
# NOTE: Não existe tratamento.
# · Observações ao nível de planta fonte.
# · Estrutura: 2 ensaios > 6 fontes/ensaio.
# · Duas respostas.
# · Modelo: resposta ~ (1 | ensaio).
xtabs(~ensaio + fonte, data = tb)
#--------------------------------------------
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 6, skip = 1)
str(tb)
# NOTE: Tratamento aplicado às plantas fonte.
# · Observações ao nível de planta fonte.
# · Estrutura: 2 ensaios > 2 tratamentos/ensaio > 3 fontes/tratamento.
# · Duas respostas.
# · Modelo: resposta ~ (1 | ensaio) + tratamento.
# · O tratamento foi casualizado nas fontes.
xtabs(~ensaio + tratamento, data = tb)2 Instante 0 · Sem tratamentos
Análise dos dados do instante 0 no qual não se aplicou tratamentos. O modelo irá considerar o efeito dos termos estruturais apenas para que sejam estimados os componentes de variância.
2.1 Dados ao nível de inseto
Respostas:
- Valor do QPCR dicotomizado em 30.
- Mortalidade do inseto 72 horas após exposição.
Estrutura experimental:
- 2 ensaios independentes.
- 6 plantas fonte por ensaio.
- 10 plantas teste por planta fonte.
- 3 insetos por planta teste.
- Total de 360 registros.
Modelo estatístico
\[ \begin{align*} y_{ijkl} &= \mu + \text{ENSAIO}_{i} + \text{FONTE}_{(i)j} + \text{TESTE}_{(ij)k} + \text{INSETO}_{(ijk)l} \\ \text{ENSAIO}_{i} &\sim \text{Normal}(0, \sigma^2_E) \\ \text{FONTE}_{(i)j} &\sim \text{Normal}(0, \sigma^2_F) \\ \text{TESTE}_{(ij)k} &\sim \text{Normal}(0, \sigma^2_T) \\ \text{INSETO}_{(ijk)l} &\sim \text{Normal}(0, \sigma^2_I) \end{align*} \]
#-----------------------------------------------------------------------
# Importação.
tb <- gdata::read.xls("DADOS 1 Experimento.xlsx", sheet = 1, skip = 1,
dec = ",") %>%
as.tibble()
# tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 1, skip = 1)
str(tb)## tibble [360 × 7] (S3: tbl_df/tbl/data.frame)
## $ ensaio : int [1:360] 1 1 1 1 1 1 1 1 1 1 ...
## $ tratamento: logi [1:360] NA NA NA NA NA NA ...
## $ fonte : int [1:360] 1 1 1 1 1 1 1 1 1 1 ...
## $ inseto : int [1:360] 1 2 3 4 5 6 7 8 9 10 ...
## $ qpcr : num [1:360] 25 29 32.6 32.2 31.9 ...
## $ teste : int [1:360] 1 1 1 2 2 2 3 3 3 4 ...
## $ morte : int [1:360] 1 1 1 0 0 0 0 0 0 1 ...tb <- tb %>%
mutate_at(c("ensaio", "fonte", "teste"), "factor")
# · Observações ao nível de inseto.
# · Estrutura: 2 ensaios > 6 fontes/ensaio > 10 testes/fonte > 3
# insetos/teste.
# · Duas respostas.
# · Modelo: resposta ~ (1 | ensaio/fonte/teste).
# xtabs(~ensaio + fonte, data = tb)
# xtabs(~teste + fonte, data = tb)
ggplot(data = tb,
mapping = aes(x = teste, y = morte, color = ensaio)) +
facet_wrap(facets = ~fonte, scale = "free_x") +
geom_jitter(height = 0.1, width = 0)
ggplot(data = tb,
mapping = aes(x = teste, y = qpcr, color = ensaio)) +
facet_wrap(facets = ~fonte, scale = "free_x") +
geom_point()
#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de qPCR.
# Dicotomizar em 30
tb$qpcr <- ifelse(tb$qpcr <= 30, 1, 0)
# Ajuste do modelo de efeitos aleatórios.
m0 <- glmer(qpcr ~ (1 | fonte/teste) + ensaio,
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) + ensaio
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 419.3 434.9 -205.7 411.3 356
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5540 -0.5962 -0.4736 0.8274 2.0969
##
## Random effects:
## Groups Name Variance Std.Dev.
## teste:fonte (Intercept) 0.01417 0.1190
## fonte (Intercept) 0.93387 0.9664
## Number of obs: 360, groups: teste:fonte, 120; fonte, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9951 0.4349 -2.288 0.02214 *
## ensaio2 1.9881 0.6286 3.163 0.00156 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## ensaio2 -0.703# qPCR média para cada ensaio.
predict(m0,
newdata = distinct(tb, ensaio),
re.form = ~0)## 1 2
## -0.9950984 0.9929715# Médias com intervalo de confiança.
emmeans(m0, specs = ~ensaio, type = "response")## ensaio prob SE df asymp.LCL asymp.UCL
## 1 0.27 0.0857 Inf 0.136 0.464
## 2 0.73 0.0882 Inf 0.529 0.866
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de mortalidade.
# Ajuste do modelo de efeitos aleatórios.
m0 <- glmer(morte ~ (1 | fonte/teste) + ensaio,
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: morte ~ (1 | fonte/teste) + ensaio
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 458.5 474.0 -225.2 450.5 356
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2641 -0.6454 -0.3996 0.6709 1.6716
##
## Random effects:
## Groups Name Variance Std.Dev.
## teste:fonte (Intercept) 1.2965 1.1387
## fonte (Intercept) 0.8861 0.9413
## Number of obs: 360, groups: teste:fonte, 120; fonte, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3959 0.4474 -0.885 0.376
## ensaio2 -0.4006 0.6387 -0.627 0.530
##
## Correlation of Fixed Effects:
## (Intr)
## ensaio2 -0.696# Mortalidade média para cada ensaio.
predict(m0,
newdata = distinct(tb, ensaio),
re.form = ~0,
type = "response")## 1 2
## 0.4022922 0.3107614# Médias com intervalo de confiança.
emmeans(m0, specs = ~ensaio, type = "response")## ensaio prob SE df asymp.LCL asymp.UCL
## 1 0.402 0.1076 Inf 0.219 0.618
## 2 0.311 0.0982 Inf 0.155 0.526
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale2.2 Dados ao nível de planta teste
#-----------------------------------------------------------------------
# Importação.
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 3, skip = 1)
str(tb)## tibble [360 × 5] (S3: tbl_df/tbl/data.frame)
## $ fonte : num [1:360] 1 1 1 1 1 1 1 1 1 1 ...
## $ tratamento: logi [1:360] NA NA NA NA NA NA ...
## $ inseto : num [1:360] 1 2 3 4 5 6 7 8 9 10 ...
## $ teste : num [1:360] 1 1 1 2 2 2 3 3 3 4 ...
## $ qpcr : num [1:360] 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(-tratamento, -inseto) %>%
mutate(ensaio = ifelse(fonte <= 6, 1, 2)) %>%
group_by(ensaio, fonte, teste) %>%
slice(1) %>%
ungroup()
tb <- tb %>%
mutate_at(c("ensaio", "fonte"), "factor")
str(tb)## tibble [120 × 4] (S3: tbl_df/tbl/data.frame)
## $ fonte : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ teste : num [1:120] 1 2 3 4 5 6 7 8 9 10 ...
## $ qpcr : num [1:120] 0 0 0 0 0 0 0 0 0 0 ...
## $ ensaio: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...# · Observações ao nível de plata teste.
# · Estrutura: 2 ensaios > 6 fontes/ensaio > 10 testes/fonte.
# · Uma respostas.
# · Modelo: resposta ~ (1 | fonte) + ensaio.
ggplot(data = tb,
mapping = aes(x = teste, y = qpcr, color = ensaio)) +
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(qpcr ~ (1 | ensaio:fonte) + ensaio,
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: qpcr ~ (1 | ensaio:fonte) + ensaio
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 88.4 96.7 -41.2 82.4 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0616 -0.1900 -0.1715 0.4851 5.2627
##
## Random effects:
## Groups Name Variance Std.Dev.
## ensaio:fonte (Intercept) 0.2193 0.4683
## Number of obs: 120, groups: ensaio:fonte, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.4642 0.7721 -4.487 7.22e-06 ***
## ensaio2 4.7127 0.8884 5.305 1.13e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## ensaio2 -0.902# Mortalidade média para cada ensaio.
predict(m0,
newdata = distinct(tb, ensaio),
re.form = ~0,
type = "response")## 1 2
## 0.03034915 0.77704960# Médias com intervalo de confiança.
emmeans(m0, specs = ~ensaio, type = "response")## ensaio prob SE df asymp.LCL asymp.UCL
## 1 0.0303 0.0227 Inf 0.00685 0.124
## 2 0.7770 0.0665 Inf 0.62160 0.881
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale2.3 Dados ao nível de planta fonte
#-----------------------------------------------------------------------
# Importação.
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 5, skip = 1)
str(tb)## tibble [12 × 5] (S3: tbl_df/tbl/data.frame)
## $ ensaio : num [1:12] 1 1 1 1 1 1 2 2 2 2 ...
## $ tratamento: logi [1:12] NA NA NA NA NA NA ...
## $ fonte : num [1:12] 1 2 3 4 5 6 7 8 9 10 ...
## $ qpcr : num [1:12] 334881 763663 1319380 365823 1064901 ...
## $ isolamento: num [1:12] 1062176 198020 496667 507895 43147 ...tb <- tb %>%
select(-tratamento) %>%
mutate_at(c("ensaio"), "factor")
str(tb)## tibble [12 × 4] (S3: tbl_df/tbl/data.frame)
## $ ensaio : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 2 2 2 2 ...
## $ fonte : num [1:12] 1 2 3 4 5 6 7 8 9 10 ...
## $ qpcr : num [1:12] 334881 763663 1319380 365823 1064901 ...
## $ isolamento: num [1:12] 1062176 198020 496667 507895 43147 ...# · Observações ao nível de plata fonte.
# · Estrutura: 2 ensaios > 6 fontes/ensaio.
# · Duas respostas.
# · Modelo: resposta ~ ensaio.
ggplot(data = tb,
mapping = aes(x = fonte, y = qpcr, color = ensaio)) +
geom_point() +
scale_y_log10()
ggplot(data = tb,
mapping = aes(x = fonte, y = isolamento, color = ensaio)) +
geom_point() +
scale_y_log10()
#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de qPCR.
# Ajuste do modelo de efeitos aleatórios.
m0 <- lm(log10(qpcr) ~ ensaio,
data = tb)
# Estimativas dos componentes de variância.
summary(m0)##
## Call:
## lm(formula = log10(qpcr) ~ ensaio, data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.80106 -0.26876 0.06969 0.38589 1.26861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.82244 0.44356 13.127 1.25e-07 ***
## ensaio2 0.07721 0.62729 0.123 0.904
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.087 on 10 degrees of freedom
## Multiple R-squared: 0.001513, Adjusted R-squared: -0.09834
## F-statistic: 0.01515 on 1 and 10 DF, p-value: 0.9045# qPCR média para cada ensaio.
predict(m0,
newdata = distinct(tb, ensaio))## 1 2
## 5.822440 5.899646# Médias com intervalo de confiança.
emmeans(m0, specs = ~ensaio)## ensaio emmean SE df lower.CL upper.CL
## 1 5.82 0.444 10 4.83 6.81
## 2 5.90 0.444 10 4.91 6.89
##
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de isolamento.
# Ajuste do modelo de efeitos aleatórios.
m0 <- lm(log10(isolamento) ~ ensaio,
data = tb)
# Estimativas dos componentes de variância.
summary(m0)##
## Call:
## lm(formula = log10(isolamento) ~ ensaio, data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8370 -0.1503 0.1145 0.2061 0.5543
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4719 0.1797 30.45 1.47e-09 ***
## ensaio2 1.9394 0.2542 7.63 6.13e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4019 on 8 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.8792, Adjusted R-squared: 0.8641
## F-statistic: 58.22 on 1 and 8 DF, p-value: 6.13e-05# qPCR média para cada ensaio.
predict(m0,
newdata = distinct(tb, ensaio))## 1 2
## 5.471939 7.411292# Médias com intervalo de confiança.
emmeans(m0, specs = ~ensaio)## ensaio emmean SE df lower.CL upper.CL
## 1 5.47 0.18 8 5.06 5.89
## 2 7.41 0.18 8 7.00 7.83
##
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.953 Instante 1 · Com tratamentos
3.1 Dados ao nível de inseto
#-----------------------------------------------------------------------
# Importação.
tb <- gdata::read.xls("DADOS 1 Experimento.xlsx", sheet = 2, skip = 1,
dec = ",") %>%
as.tibble()
# tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 2, skip = 1)
str(tb)## tibble [360 × 7] (S3: tbl_df/tbl/data.frame)
## $ ensaio : int [1:360] 1 1 1 1 1 1 1 1 1 1 ...
## $ tratamento: chr [1:360] "CONTROLE" "CONTROLE" "CONTROLE" "CONTROLE" ...
## $ fonte : int [1:360] 5 5 5 5 5 5 5 5 5 5 ...
## $ inseto : int [1:360] 361 362 363 364 365 366 367 368 369 370 ...
## $ qpcr : num [1:360] 35.6 33.6 31.4 34.4 34.2 ...
## $ teste : int [1:360] 1 1 1 2 2 2 3 3 3 4 ...
## $ morte : int [1:360] 0 0 1 0 1 1 0 0 0 0 ...tb <- tb %>%
mutate_at(c("ensaio", "tratamento", "fonte", "teste"), "factor")
# · Observações ao nível de inseto.
# · Estrutura: 2 ensaios > 2 tratamentos/ensaio > 6 fontes/ensaio > 10
# testes/fonte > 3 insetos/teste.
# · Duas respostas.
# · Modelo: resposta ~ (1 | ensaio/fonte/teste) + tratamento.
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 = ensaio)) +
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) + ensaio + 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) + ensaio + tratamento
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 391.6 411.0 -190.8 381.6 355
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0337 -0.6727 -0.2358 0.8175 4.1841
##
## Random effects:
## Groups Name Variance Std.Dev.
## teste:fonte (Intercept) 0.02708 0.1646
## fonte (Intercept) 0.75727 0.8702
## Number of obs: 360, groups: teste:fonte, 123; fonte, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4429 0.5590 -4.370 1.24e-05 ***
## ensaio2 2.0550 0.5954 3.451 0.000558 ***
## tratamentoNAC 1.2053 0.5899 2.043 0.041042 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ensai2
## ensaio2 -0.648
## tratamntNAC -0.632 0.124# 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.195 0.0677 Inf 0.0946 0.361 1
## NAC 0.448 0.0989 Inf 0.2701 0.640 2
##
## Results are averaged over the levels of: ensaio
## 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.05emm <- 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) + ensaio + 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: morte ~ (1 | fonte/teste) + ensaio + tratamento
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 337.2 356.6 -163.6 327.2 355
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1842 -0.4300 -0.3074 -0.1651 3.7007
##
## Random effects:
## Groups Name Variance Std.Dev.
## teste:fonte (Intercept) 0.6634 0.8145
## fonte (Intercept) 0.7278 0.8531
## Number of obs: 360, groups: teste:fonte, 123; fonte, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4894 0.5374 -2.772 0.00558 **
## ensaio2 -1.4808 0.6291 -2.354 0.01857 *
## tratamentoNAC 0.8137 0.6222 1.308 0.19097
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ensai2
## ensaio2 -0.441
## tratamntNAC -0.623 -0.052# 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.0971 0.0428 Inf 0.0397 0.219 1
## NAC 0.1953 0.0683 Inf 0.0938 0.363 1
##
## Results are averaged over the levels of: ensaio
## 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.05emm <- 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")
3.2 Dados ao nível de planta teste
#-----------------------------------------------------------------------
# Importação.
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 4, skip = 1)
str(tb)## tibble [360 × 6] (S3: tbl_df/tbl/data.frame)
## $ ensaio : num [1:360] 1 1 1 1 1 1 1 1 1 1 ...
## $ tratamento: chr [1:360] "CONTROLE" "CONTROLE" "CONTROLE" "CONTROLE" ...
## $ fonte : num [1:360] 5 5 5 5 5 5 5 5 5 5 ...
## $ inseto : num [1:360] 361 362 363 364 365 366 367 368 369 370 ...
## $ teste : num [1:360] 121 121 121 122 122 122 123 123 123 124 ...
## $ qpcr : num [1:360] 1 1 1 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 [120 × 5] (S3: tbl_df/tbl/data.frame)
## $ ensaio : Factor w/ 2 levels "1","2": 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 ...
## $ fonte : Factor w/ 12 levels "1","2","3","4",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ teste : num [1:120] 141 142 143 144 145 146 147 148 149 150 ...
## $ qpcr : num [1:120] 0 0 0 0 0 0 0 0 0 0 ...# · Observações ao nível de plata teste.
# · Estrutura: 2 ensaios > 2 tratamentos/ensaio >6 fontes/ensaio > 10
# testes/fonte.
# · Uma respostas.
# · Modelo: resposta ~ (1 | fonte) + ensaio + tratamento.
ggplot(data = tb,
mapping = aes(x = teste, y = qpcr, 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(qpcr ~ (1 | ensaio:fonte) + ensaio + 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: qpcr ~ (1 | ensaio:fonte) + ensaio + tratamento
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 56.0 67.2 -24.0 48.0 116
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6326 -0.1898 -0.1407 -0.1407 5.2675
##
## Random effects:
## Groups Name Variance Std.Dev.
## ensaio:fonte (Intercept) 1.652 1.285
## Number of obs: 120, groups: ensaio:fonte, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.4002 1.2960 -3.395 0.000686 ***
## ensaio2 0.7989 1.2526 0.638 0.523631
## tratamentoNAC 0.7989 1.2526 0.638 0.523624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ensai2
## ensaio2 -0.517
## tratamntNAC -0.517 -0.015# 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.0180 0.0196 Inf 0.00207 0.139 1
## NAC 0.0391 0.0394 Inf 0.00518 0.241 1
##
## Results are averaged over the levels of: ensaio
## 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.05emm <- 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 = "qPCR",
y = "Tratamento")
3.3 Dados ao nível de planta fonte
#-----------------------------------------------------------------------
# Importação.
tb <- read_excel("DADOS 1 Experimento.xlsx", sheet = 6, skip = 1)
str(tb)## tibble [12 × 5] (S3: tbl_df/tbl/data.frame)
## $ ensaio : num [1:12] 1 1 1 1 1 1 2 2 2 2 ...
## $ tratamento: chr [1:12] "NAC" "NAC" "NAC" "CONTROLE" ...
## $ fonte : num [1:12] 1 2 3 4 5 6 7 8 9 10 ...
## $ qpcr : num [1:12] 47658 120528 51049 648237 6150814 ...
## $ isolamento: num [1:12] 452586 NA NA 6300813 3349057 ...tb <- tb %>%
mutate_at(c("ensaio", "tratamento"), "factor")
str(tb)## tibble [12 × 5] (S3: tbl_df/tbl/data.frame)
## $ ensaio : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 2 2 2 2 ...
## $ tratamento: Factor w/ 2 levels "CONTROLE","NAC": 2 2 2 1 1 1 1 2 1 2 ...
## $ fonte : num [1:12] 1 2 3 4 5 6 7 8 9 10 ...
## $ qpcr : num [1:12] 47658 120528 51049 648237 6150814 ...
## $ isolamento: num [1:12] 452586 NA NA 6300813 3349057 ...# · Observações ao nível de plata fonte.
# · Estrutura: 2 ensaios > 6 fontes/ensaio.
# · Duas respostas.
# · Modelo: resposta ~ ensaio.
ggplot(data = tb,
mapping = aes(x = fonte, y = qpcr, color = tratamento)) +
geom_point() +
scale_y_log10()
ggplot(data = tb,
mapping = aes(x = fonte, y = isolamento, color = tratamento)) +
geom_point() +
scale_y_log10()
#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de qPCR.
# Ajuste do modelo de efeitos aleatórios.
m0 <- lm(log10(qpcr) ~ ensaio + tratamento,
data = tb)
# Estimativas dos componentes de variância.
summary(m0)##
## Call:
## lm(formula = log10(qpcr) ~ ensaio + tratamento, data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.89512 -0.22898 0.00398 0.39935 1.35643
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1431 0.4304 14.271 1.74e-07 ***
## ensaio2 -1.6103 0.4970 -3.240 0.0102 *
## tratamentoNAC -1.2701 0.4970 -2.555 0.0309 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8609 on 9 degrees of freedom
## Multiple R-squared: 0.6542, Adjusted R-squared: 0.5773
## F-statistic: 8.513 on 2 and 9 DF, p-value: 0.008409# Médias com intervalo de confiança.
emm <- emmeans(m0, specs = ~tratamento) %>%
multcomp::cld()## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimatesemm## tratamento emmean SE df lower.CL upper.CL .group
## NAC 4.07 0.351 9 3.27 4.86 1
## CONTROLE 5.34 0.351 9 4.54 6.13 2
##
## Results are averaged over the levels of: ensaio
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
## Note: contrasts are still on the log10 scale
## significance level used: alpha = 0.05emm <- emm %>%
as.data.frame() %>%
mutate(cld = as_letters(.group))
ggplot(data = emm,
mapping = aes(x = emmean, y = tratamento)) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.2f %s", emmean, cld)),
vjust = -0.3) +
geom_errorbarh(mapping = aes(xmin = lower.CL, xmax = upper.CL),
height = 0.1) +
labs(x = "Log10 do qPCR",
y = "Tratamento")
#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de isolamento.
# Ajuste do modelo de efeitos aleatórios.
m0 <- lm(log10(isolamento) ~ ensaio + tratamento,
data = tb)
# Estimativas dos componentes de variância.
summary(m0)##
## Call:
## lm(formula = log10(isolamento) ~ ensaio + tratamento, data = tb)
##
## Residuals:
## 1 4 5 6 8 9 11
## 0.01730 0.16458 -0.10990 -0.07198 0.23688 0.01730 -0.25418
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.6348 0.1104 60.115 4.59e-07 ***
## ensaio2 0.2655 0.1710 1.553 0.19539
## tratamentoNAC -0.9964 0.1710 -5.828 0.00432 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2035 on 4 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.8971, Adjusted R-squared: 0.8457
## F-statistic: 17.44 on 2 and 4 DF, p-value: 0.01058# Médias com intervalo de confiança.
emm <- emmeans(m0, specs = ~tratamento) %>%
multcomp::cld()## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimatesemm## tratamento emmean SE df lower.CL upper.CL .group
## NAC 5.77 0.121 4 5.44 6.11 1
## CONTROLE 6.77 0.110 4 6.46 7.07 2
##
## Results are averaged over the levels of: ensaio
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
## Note: contrasts are still on the log10 scale
## significance level used: alpha = 0.05emm <- emm %>%
as.data.frame() %>%
mutate(cld = as_letters(.group))
ggplot(data = emm,
mapping = aes(x = emmean, y = tratamento)) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.2f %s", emmean, cld)),
vjust = -0.3) +
geom_errorbarh(mapping = aes(xmin = lower.CL, xmax = upper.CL),
height = 0.1) +
labs(x = "Log10 do isolamento",
y = "Tratamento")
4 Análise conjunta dos instantes
OBS: fazer a análise conjunta dos instantes.
- Comparar o controle do instante 1 com os resultados do instante 0 para verificar efeito de instante.
- As mesmas 12 plantas do instante 0 foram reutilizadas no instante 1.
4.1 Análise conjunta ao nível de insetos
#-----------------------------------------------------------------------
# Importação e preparo.
# Importa as tabelas.
# tb <- list("A" = read_excel("DADOS 1 Experimento.xlsx",
# sheet = 1, skip = 1),
# "B" = read_excel("DADOS 1 Experimento.xlsx",
# sheet = 2, skip = 1))
tb <- list("A" = as.tibble(gdata::read.xls("DADOS 1 Experimento.xlsx",
sheet = 1, skip = 1, dec = ",")),
"B" = as.tibble(gdata::read.xls("DADOS 1 Experimento.xlsx",
sheet = 2, skip = 1, dec = ",")))
# Concatenha por linhas.
tb <- bind_rows(tb, .id = "instante")
# Coloca o nível controle para tratamento no instante 0.
tb <- tb %>%
replace_na(replace = list(tratamento = "CONTROLE"))
# Cria a condição experimental combinando níveis.
tb <- tb %>%
mutate(condicao = interaction(instante, tratamento, drop = TRUE),
qpcr = ifelse(qpcr <= 30, 1, 0))
# Converte fontes de variação para fator.
tb <- tb %>%
mutate_at(c("instante", "ensaio", "tratamento",
"condicao", "fonte", "teste"),
"factor")
# Tabela de frequência.
tb %>%
count(instante, ensaio, tratamento, condicao, fonte) %>%
arrange(instante, ensaio, tratamento, condicao, fonte) %>%
print(n = Inf)## # A tibble: 24 x 6
## instante ensaio tratamento condicao fonte n
## <fct> <fct> <fct> <fct> <fct> <int>
## 1 A 1 CONTROLE A.CONTROLE 1 30
## 2 A 1 CONTROLE A.CONTROLE 2 30
## 3 A 1 CONTROLE A.CONTROLE 3 30
## 4 A 1 CONTROLE A.CONTROLE 4 30
## 5 A 1 CONTROLE A.CONTROLE 5 30
## 6 A 1 CONTROLE A.CONTROLE 6 30
## 7 A 2 CONTROLE A.CONTROLE 7 30
## 8 A 2 CONTROLE A.CONTROLE 8 30
## 9 A 2 CONTROLE A.CONTROLE 9 30
## 10 A 2 CONTROLE A.CONTROLE 10 30
## 11 A 2 CONTROLE A.CONTROLE 11 30
## 12 A 2 CONTROLE A.CONTROLE 12 30
## 13 B 1 CONTROLE B.CONTROLE 4 30
## 14 B 1 CONTROLE B.CONTROLE 5 30
## 15 B 1 CONTROLE B.CONTROLE 6 30
## 16 B 1 NAC B.NAC 1 30
## 17 B 1 NAC B.NAC 2 30
## 18 B 1 NAC B.NAC 3 30
## 19 B 2 CONTROLE B.CONTROLE 9 30
## 20 B 2 CONTROLE B.CONTROLE 11 30
## 21 B 2 CONTROLE B.CONTROLE 12 30
## 22 B 2 NAC B.NAC 7 30
## 23 B 2 NAC B.NAC 8 30
## 24 B 2 NAC B.NAC 10 30ggplot(data = tb,
mapping = aes(x = teste,
y = qpcr,
color = condicao,
shape = interaction(instante, ensaio))) +
facet_wrap(facets = ~fonte, scale = "free_x") +
geom_point()
#-----------------------------------------------------------------------
# Análise para qPCR.
# Modelo que acomoda o termos estruturais com efeitos aleatórios.
m0 <- glmer(qpcr ~
(1 | instante:ensaio) +
(1 | instante:ensaio:fonte) +
(1 | instante:ensaio:fonte:teste) +
condicao,
data = tb,
family = binomial)
# Quadro de testes de Wald para os termos experimetais.
anova(m0)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## condicao 2 3.995 1.9975 1.9975# Resumo do modelo com estimativas dos parâmetros.
summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: qpcr ~ (1 | instante:ensaio) + (1 | instante:ensaio:fonte) +
## (1 | instante:ensaio:fonte:teste) + condicao
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 815.7 843.2 -401.8 803.7 714
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2208 -0.6658 -0.3783 0.8130 4.1486
##
## Random effects:
## Groups Name Variance Std.Dev.
## instante:ensaio:fonte:teste (Intercept) 0.01162 0.1078
## instante:ensaio:fonte (Intercept) 1.06433 1.0317
## instante:ensaio (Intercept) 0.83734 0.9151
## Number of obs: 720, groups:
## instante:ensaio:fonte:teste, 243; instante:ensaio:fonte, 24; instante:ensaio, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.003259 0.724359 0.004 0.996
## condicaoB.CONTROLE -1.421920 1.084899 -1.311 0.190
## condicaoB.NAC -0.195588 1.073811 -0.182 0.855
##
## Correlation of Fixed Effects:
## (Intr) cB.CON
## cB.CONTROLE -0.668
## condicB.NAC -0.674 0.809# Comparações múltiplas entre as condições experimentais.
emmeans(m0, spec = ~condicao, type = "response") %>%
multcomp::cld()## condicao prob SE df asymp.LCL asymp.UCL .group
## B.CONTROLE 0.195 0.127 Inf 0.0474 0.541 1
## B.NAC 0.452 0.197 Inf 0.1483 0.796 1
## A.CONTROLE 0.501 0.181 Inf 0.1952 0.806 1
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05#-----------------------------------------------------------------------
# Análise para mortalidade.
# Modelo que acomoda o termos estruturais com efeitos aleatórios.
m0 <- glmer(morte ~
(1 | instante:ensaio) +
(1 | instante:ensaio:fonte) +
(1 | instante:ensaio:fonte:teste) +
condicao,
data = tb,
family = binomial)
# Quadro de testes de Wald para os termos experimetais.
anova(m0)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## condicao 2 7.909 3.9545 3.9545# Resumo do modelo com estimativas dos parâmetros.
summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: morte ~ (1 | instante:ensaio) + (1 | instante:ensaio:fonte) +
## (1 | instante:ensaio:fonte:teste) + condicao
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 795.9 823.4 -392.0 783.9 714
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2856 -0.5776 -0.3174 0.6648 3.0036
##
## Random effects:
## Groups Name Variance Std.Dev.
## instante:ensaio:fonte:teste (Intercept) 1.00852 1.004
## instante:ensaio:fonte (Intercept) 1.06790 1.033
## instante:ensaio (Intercept) 0.01561 0.125
## Number of obs: 720, groups:
## instante:ensaio:fonte:teste, 243; instante:ensaio:fonte, 24; instante:ensaio, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5852 0.3528 -1.659 0.09717 .
## condicaoB.CONTROLE -1.7329 0.6402 -2.707 0.00679 **
## condicaoB.NAC -0.8666 0.6112 -1.418 0.15621
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cB.CON
## cB.CONTROLE -0.527
## condicB.NAC -0.564 0.344# Comparações múltiplas entre as condições experimentais.
emmeans(m0, spec = ~condicao, type = "response") %>%
multcomp::cld()## condicao prob SE df asymp.LCL asymp.UCL .group
## B.CONTROLE 0.0896 0.0444 Inf 0.0328 0.222 1
## B.NAC 0.1897 0.0776 Inf 0.0801 0.386 12
## A.CONTROLE 0.3577 0.0811 Inf 0.2181 0.527 2
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.054.2 Análise conjunta ao nível de plantas teste
#-----------------------------------------------------------------------
# Importação e preparo.
# Importa as tabelas.
tb <- list("A" = read_excel("DADOS 1 Experimento.xlsx",
sheet = 3, skip = 1),
"B" = read_excel("DADOS 1 Experimento.xlsx",
sheet = 4, skip = 1))
# Concatenha por linhas.
tb <- bind_rows(tb, .id = "instante")
# Coloca o nível controle para tratamento no instante 0.
tb <- tb %>%
replace_na(replace = list(tratamento = "CONTROLE"))
# Cria a condição experimental combinando níveis.
tb <- tb %>%
mutate(condicao = interaction(instante, tratamento, drop = TRUE),
ensaio = ifelse(fonte < 7, 1, 2))
# Converte fontes de variação para fator.
tb <- tb %>%
mutate_at(c("instante", "ensaio", "tratamento",
"condicao", "fonte"),
"factor")
# NOTE: Não tem variância no qpcr dentro de cada cela.
tb %>%
group_by(instante, ensaio, tratamento, condicao, fonte, teste) %>%
summarise_at("qpcr", "sd") %>%
pull(qpcr)## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [34] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [67] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [100] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [133] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [166] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [199] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [232] 0 0 0 0 0 0 0 0 0# DANGER! Os dados estão ao nível de inseto incorretamente.
tb <- tb %>%
group_by(instante, ensaio, tratamento, condicao, fonte, teste) %>%
slice(1) %>%
ungroup()
# Tabela de frequência.
tb %>%
count(instante, ensaio, tratamento, condicao, fonte) %>%
arrange(instante, ensaio, tratamento, condicao, fonte) %>%
print(n = Inf)## # A tibble: 24 x 6
## instante ensaio tratamento condicao fonte n
## <fct> <fct> <fct> <fct> <fct> <int>
## 1 A 1 CONTROLE A.CONTROLE 1 10
## 2 A 1 CONTROLE A.CONTROLE 2 10
## 3 A 1 CONTROLE A.CONTROLE 3 10
## 4 A 1 CONTROLE A.CONTROLE 4 10
## 5 A 1 CONTROLE A.CONTROLE 5 10
## 6 A 1 CONTROLE A.CONTROLE 6 10
## 7 A 2 CONTROLE A.CONTROLE 7 10
## 8 A 2 CONTROLE A.CONTROLE 8 10
## 9 A 2 CONTROLE A.CONTROLE 9 10
## 10 A 2 CONTROLE A.CONTROLE 10 10
## 11 A 2 CONTROLE A.CONTROLE 11 10
## 12 A 2 CONTROLE A.CONTROLE 12 10
## 13 B 1 CONTROLE B.CONTROLE 4 10
## 14 B 1 CONTROLE B.CONTROLE 5 10
## 15 B 1 CONTROLE B.CONTROLE 6 10
## 16 B 1 NAC B.NAC 1 10
## 17 B 1 NAC B.NAC 2 10
## 18 B 1 NAC B.NAC 3 10
## 19 B 2 CONTROLE B.CONTROLE 9 10
## 20 B 2 CONTROLE B.CONTROLE 11 10
## 21 B 2 CONTROLE B.CONTROLE 12 10
## 22 B 2 NAC B.NAC 7 10
## 23 B 2 NAC B.NAC 8 10
## 24 B 2 NAC B.NAC 10 10ggplot(data = tb,
mapping = aes(x = fonte,
y = qpcr,
color = condicao)) +
facet_wrap(facets = ~interaction(instante, ensaio), scale = "free_x") +
geom_jitter(height = 0.1, width = 0.05)
#-----------------------------------------------------------------------
# Análise para PCR.
# Modelo que acomoda o termos estruturais com efeitos aleatórios.
m0 <- glmer(qpcr ~
(1 | instante:ensaio) +
(1 | instante:ensaio:fonte) +
condicao,
data = tb,
family = binomial)
# Quadro de testes de Wald para os termos experimetais.
anova(m0)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## condicao 2 2.1709 1.0854 1.0854# Resumo do modelo com estimativas dos parâmetros.
summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: qpcr ~ (1 | instante:ensaio) + (1 | instante:ensaio:fonte) +
## condicao
## Data: tb
##
## AIC BIC logLik deviance df.resid
## 153.8 171.2 -71.9 143.8 235
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2630 -0.2446 -0.1901 -0.1306 5.3710
##
## Random effects:
## Groups Name Variance Std.Dev.
## instante:ensaio:fonte (Intercept) 0.8527 0.9234
## instante:ensaio (Intercept) 2.9841 1.7275
## Number of obs: 240, groups:
## instante:ensaio:fonte, 24; instante:ensaio, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.086 1.311 -0.828 0.407
## condicaoB.CONTROLE -2.815 2.006 -1.404 0.160
## condicaoB.NAC -1.893 1.937 -0.977 0.328
##
## Correlation of Fixed Effects:
## (Intr) cB.CON
## cB.CONTROLE -0.646
## condicB.NAC -0.667 0.845# Comparações múltiplas entre as condições experimentais.
emmeans(m0, spec = ~condicao, type = "response") %>%
multcomp::cld()## condicao prob SE df asymp.LCL asymp.UCL .group
## B.CONTROLE 0.0198 0.0297 Inf 0.00100 0.289 1
## B.NAC 0.0484 0.0664 Inf 0.00299 0.462 1
## A.CONTROLE 0.2523 0.2474 Inf 0.02518 0.815 1
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.054.3 Análise conjunta ao nível de plantas fonte
#-----------------------------------------------------------------------
# Importação e preparo.
# Importa as tabelas.
tb <- list("A" = read_excel("DADOS 1 Experimento.xlsx",
sheet = 5, skip = 1),
"B" = read_excel("DADOS 1 Experimento.xlsx",
sheet = 6, skip = 1))
# Concatenha por linhas.
tb <- bind_rows(tb, .id = "instante")
# Coloca o nível controle para tratamento no instante 0.
tb <- tb %>%
replace_na(replace = list(tratamento = "CONTROLE"))
# Cria a condição experimental combinando níveis.
tb <- tb %>%
mutate(condicao = interaction(instante, tratamento, drop = TRUE))
# Converte fontes de variação para fator.
tb <- tb %>%
mutate_at(c("instante", "ensaio", "tratamento",
"condicao", "fonte"),
"factor")
# Troca escala da variável resposta.
# tb <- tb %>%
# mutate(qpcr = qpcr/10000)
# Tabela de frequência.
tb %>%
count(instante, ensaio, tratamento, condicao) %>%
arrange(instante, ensaio, tratamento, condicao) %>%
print(n = Inf)## # A tibble: 6 x 5
## instante ensaio tratamento condicao n
## <fct> <fct> <fct> <fct> <int>
## 1 A 1 CONTROLE A.CONTROLE 6
## 2 A 2 CONTROLE A.CONTROLE 6
## 3 B 1 CONTROLE B.CONTROLE 3
## 4 B 1 NAC B.NAC 3
## 5 B 2 CONTROLE B.CONTROLE 3
## 6 B 2 NAC B.NAC 3ggplot(data = tb,
mapping = aes(x = fonte,
y = qpcr,
color = condicao)) +
facet_wrap(facets = ~interaction(instante, ensaio), scale = "free_x") +
geom_jitter(height = 0.1, width = 0.05) +
scale_y_log10()
ggplot(data = tb,
mapping = aes(x = fonte,
y = isolamento,
color = condicao)) +
facet_wrap(facets = ~interaction(instante, ensaio), scale = "free_x") +
geom_jitter(height = 0.1, width = 0.05) +
scale_y_log10()
#-----------------------------------------------------------------------
# Análise para qPCR.
# Modelo que acomoda o termos estruturais com efeitos aleatórios.
m0 <- lmer(log10(qpcr) ~
(1 | instante:ensaio) +
(1 | fonte) +
condicao,
data = tb)
# Quadro de testes de Wald para os termos experimetais.
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## condicao 6.8277 3.4138 2 2 3.5161 0.2214# Resumo do modelo com estimativas dos parâmetros.
summary(m0)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(qpcr) ~ (1 | instante:ensaio) + (1 | fonte) + condicao
## Data: tb
##
## REML criterion at convergence: 67.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.82794 -0.16644 0.04518 0.32266 1.29479
##
## Random effects:
## Groups Name Variance Std.Dev.
## fonte (Intercept) 0.001452 0.0381
## instante:ensaio (Intercept) 0.487744 0.6984
## Residual 0.970917 0.9854
## Number of obs: 24, groups: fonte, 12; instante:ensaio, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.8610 0.5700 1.9996 10.282 0.00933 **
## condicaoB.CONTROLE -0.5237 0.8547 2.5114 -0.613 0.59101
## condicaoB.NAC -1.7926 0.8547 2.5114 -2.097 0.14436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cB.CON
## cB.CONTROLE -0.667
## condicB.NAC -0.667 0.778# Comparações múltiplas entre as condições experimentais.
emmeans(m0, spec = ~condicao) %>%
multcomp::cld()## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## condicao emmean SE df lower.CL upper.CL .group
## B.NAC 4.07 0.651 3.1 2.03 6.10 1
## B.CONTROLE 5.34 0.651 3.1 3.30 7.37 1
## A.CONTROLE 5.86 0.570 2.0 3.41 8.31 1
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
## Note: contrasts are still on the log10 scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05#-----------------------------------------------------------------------
# Análise para isolamento.
# Modelo que acomoda o termos estruturais com efeitos aleatórios.
m0 <- lmer(log10(isolamento) ~
(1 | instante:ensaio) +
(1 | fonte) +
condicao,
data = tb)
# Quadro de testes de Wald para os termos experimetais.
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## condicao 1.279 0.63948 2 2 7.943 0.1118# Resumo do modelo com estimativas dos parâmetros.
summary(m0)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log10(isolamento) ~ (1 | instante:ensaio) + (1 | fonte) + condicao
## Data: tb
##
## REML criterion at convergence: 21.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1647 -0.2645 0.1571 0.4382 1.3768
##
## Random effects:
## Groups Name Variance Std.Dev.
## fonte (Intercept) 0.03713 0.1927
## instante:ensaio (Intercept) 0.94952 0.9744
## Residual 0.08051 0.2837
## Number of obs: 17, groups: fonte, 11; instante:ensaio, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.4382 0.6975 1.9905 9.231 0.0117 *
## condicaoB.CONTROLE 0.3738 0.9937 2.0500 0.376 0.7421
## condicaoB.NAC -0.7147 0.9962 2.0709 -0.717 0.5454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cB.CON
## cB.CONTROLE -0.698
## condicB.NAC -0.695 0.962# Comparações múltiplas entre as condições experimentais.
emmeans(m0, spec = ~condicao) %>%
multcomp::cld()## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## condicao emmean SE df lower.CL upper.CL .group
## B.NAC 5.72 0.727 2.19 2.84 8.61 1
## A.CONTROLE 6.44 0.698 1.97 3.39 9.49 12
## B.CONTROLE 6.81 0.720 2.14 3.89 9.73 2
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
## Note: contrasts are still on the log10 scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05