sensitivity_ake_b.Rmd 16.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
---
title: >
  Effect of fungicide sprays programs and pistachio hedging on
  sensitivity of *Alternaria alternata* to fluopyram, penthiopyrad and
  fluxapyroxad in *Pistachio orchard* of Tulare County, California
author: >
  [Paulo S. F. Lichtemberg](http://lattes.cnpq.br/8132272273348880)</br>
  Ryan D. Puckett</br>
  [Walmes M. Zeviani](http://www.leg.ufpr.br/doku.php/pessoais:walmes)</br>
  Connor G. Cunningham</br>
  Themis J. Michailides
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{Effect of fungicide sprays programs and pistachio hedging on sensitivity of Alternaria alternata}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Session Definition

```{r, message = FALSE, results = "hide"}
# https://github.com/walmes/wzRfun
# devtools::install_github("walmes/wzRfun")

library(lattice)
library(latticeExtra)
27
library(plyr)
28
library(wzRfun)
29 30 31 32
library(lme4)
library(lmerTest)
library(doBy)
library(multcomp)
33 34 35 36 37 38
```
```{r, eval = FALSE}
library(wzCoop)
```
```{r setup, include = FALSE}
source("config/setup.R")
Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
39 40 41 42
tbn_ <- captioner(prefix = "Table")
fgn_ <- captioner(prefix = "Figure")
tbl_ <- function(label) tbn_(label, display = "cite")
fgl_ <- function(label) fgn_(label, display = "cite")
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
```

# Experiment and Data Description

# Exploratory Analysis

```{r}
data(sensitivity_ake_b)
str(sensitivity_ake_b)

# Short object names are handy.
sen <- sensitivity_ake_b

# Divide diameters by 100 to convert to milimeters. Calculate mean
# diameter (dm).
sen <- within(sen, {
    d1 <- d1/100
    d2 <- d2/100
    dm <- (d1 + d2)/2
62
    ar <- pi * (dm/2)^2
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
})

# Population along years.
xtabs(~pop + yr, data = sen)

# Number of observations per factor combination.
ftable(xtabs(~yr + hed + tra + fun, data = sen))

# Number of isolates per factor combination.
with(sen,
     ftable(tapply(iso,
                   INDEX = list(yr, hed, tra, fun),
                   FUN = function(x) {
                       length(unique(x))
                   })))
78 79 80 81 82 83 84

# Isolates x fungicide missing cells.
xt <- xtabs(complete.cases(cbind(dm, dos)) ~ iso + fun,
            data = sen)
i <- xt == 0
sum(i)
xt[rowSums(i) > 0, ]
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
```

```{r}
xyplot(d2 ~ d1,
       data = sen,
       aspect = "iso",
       xlab = "First diameter (mm)",
       ylab = "Second diameter (mm)") +
    layer(panel.abline(a = 0, b = 1))

xyplot(d2 ~ d1 | iso,
       data = sen,
       aspect = "iso",
       groups = fun,
       strip = FALSE,
       as.table = TRUE,
       xlab = "First diameter (mm)",
       ylab = "Second diameter (mm)") +
    layer(panel.abline(a = 0, b = 1))
```

```{r, fig.cap = cap, fig.show = "hold", echo = -(1:2)}
cap <-
"Scatter plots of mean diamenter as function of
 dose grouped by *in vitro* fungicide in natual scale
 (top) and log-log scale (bottom)."
cap <- fgn_("dm-x-dos", cap)

# Natural scales.
xyplot(dm ~ dos | iso,
       data = sen,
       groups = fun,
       type = c("p", "a"),
       strip = FALSE,
       as.table = TRUE,
       xlab = "First diameter (mm)",
       ylab = "Second diameter (mm)")

# Log-log scales.
xyplot(dm ~ dos | iso,
       scales = list(log = TRUE),
       data = sen,
       groups = fun,
       type = c("p", "a"),
       strip = FALSE,
       as.table = TRUE,
       xlab = "First diameter (mm)",
       ylab = "Second diameter (mm)")
```

```{r, fig.cap = cap, fig.show = "hold", echo = -(1:2), fig.height = 12}
cap <-
"Scatter plot of mean diamenter as function of
 dose power 1/5 grouped by *in vitro* fungicide."
cap <- fgn_("dm-x-dos0.2", cap)

# 5th root of dose.
xyplot(dm ~ dos^0.2 | iso,
       data = sen,
       groups = fun,
       type = c("p", "a"),
       strip = FALSE,
       as.table = TRUE,
       xlab = "First diameter (mm)",
       ylab = "Second diameter (mm)")
```

Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
152 153
`r fgl_("dm-x-dos")` (top) shows that doses are very skewed.
`r fgl_("dm-x-dos")` (bottom) shows that in the log-log scale
154 155
there isn't a linear relation between mean diameter and fungicide dose.

Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
156
`r fgl_("dm-x-dos0.2")` shows that, under the transformed
157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
5th-root scale, doses levels are close to equally spaced. In fact, 0.2
was found by minimization of the variance of distance between doses
($\sigma^2$) a power transformation ($p$) of dose rescaled to a unit
interval as described by the steps below
$$
\begin{align*}
  z_i &= x_i^p\\
  u_i &= \frac{z_i - \min(z)}{\max(z) -\min(z)}, \text{then } u_i \in [0, 1] \\
  d_i &= u_{i+1} - u_{i}\\
  \bar{d} &= \sum_{i=1}^{k-1} d_i/k \\
 \sigma^2 &= \sum_{i=1}^{k-1} \frac{(d_i - \bar{d})^2}{k-2}
\end{align*}
$$
where $x$ are doses in natural scale, $z$ are doses power transformed,
$u$ are scaled to a unit interval, $d$ are diferences between doses,
$\bar{d}$ is the mean difference and $\sigma^2$ is the variance of
differences.

```{r}
# Unique fungicide dose levels.
x <- sort(unique(sen$dos))
x

# Variance of distance between doses scaled to a unit interval.
esp <- function(p) {
    u <- x^p
    u <- (u - min(u))
    u <- u/max(u)
    var(diff(u))
}

# Optimise de power parameter to the most equally spaced set.
op <- optimize(f = esp, interval = c(0, 1))
op$minimum

p <- seq(0, 1, by = 0.01)
v <- sapply(p, esp)
plot(log(v) ~ p, type = "o")
abline(v = op$minimum)
```

Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
198
So $x^{0.2}$ is the most equally spaced set obtained with a power
199 200 201 202 203
transformation. Equally spaced levels are beneficial beacause reduce
problems related to leverage.

# Half Effective Concentration (EC~50~) Estimation

Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
A cubic spline is function constructed of piecewise third-order
polynomials which pass through a set of $m + 1$ knots. These knots spans
the observed domain of the continous factor $x$, so the set of knots is
$$
  K = \{\xi_0, \xi_1, \ldots, \xi_m\}.
$$

A function $s(x)$ is a cubic spline if it is made of cubic polynomials
$s_{i}(x)$ in each interval $[x_{i-1}, x_{m}]$, $i = 1, \ldots, m$.

Those adjacent cubic pylinomials pieces must bind and be smooth at the
internal knots , so additional constrais are made to result in a
composite continuous smooth function. Requering continous derivatives,
we ensure that the resulting function is as smooth as possible.

For natural splines, two aditional boundary conditions are made
$$
  s^{''}_{1}(x) = 0, \quad s^{''}_{m}(x) = 0,
$$
that is, the pieces at borders aren't cubic but instead linear.

225 226 227 228 229 230 231 232
Natural cubic splines were used to estimate the half effective
concentration (EC~50~). A non linear model is usually applied in this
context but wasn't found a non linear model flexible enough to give a
good fit either a satisfactory convergence rate. So, despite splines
haven't a model equation, they are vey flexible and numerical
root-finding algorithms can e used to compute EG~50~ based on a linear
interpolated function on a predicted grid. Also, area under the
sensibility curve (AUSC) were computed by numerical integration under
Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
233
studied dose domain.
234

235 236 237 238 239 240 241
```{r}
sen$iso <- factor(sen$iso, levels = sort(unique(sen$iso)))
sen$ue <- with(sen, interaction(iso, fun, drop = TRUE))
sen$doz <- sen$dos^0.2

# A data frame without dose.
senu <- unique(subset(sen,
Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
242
                      select = c(ue, iso, fun, tra, hed, pop, yr, plot)))
243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258

# Splines.
library(splines)

# Function for fitting splines.
fitting <- function(x) {
    x <- na.omit(x)
    if (nrow(x) > 4) {
        n0 <- lm(dm ~ ns(doz, df = 3), data = x)
        yr <- range(x$dm)
        yfit <- predict(n0, newdata = pred)
        return(cbind(xfit = pred$doz, yfit = yfit))
    } else {
        NULL
    }
}
259

260 261
# Values of dose to predict diameter.
pred <- data.frame(doz = seq(0, max(sen$doz), length.out = 30))
262

263 264
# Appling to all experimental unit.
res <- ddply(sen, .(ue), .fun = fitting)
265

266 267 268
# Merge to pair `iso` and `fun`.
res <- merge(res, senu, by = "ue")
```
269 270

```{r}
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
# Estimatinf EC50 and area under curve.
get_ec50 <- function(x, y, interval = c(0, max(sen$doz))) {
    fa <- approxfun(x = x, y = y)
    int <- integrate(fa,
                     lower = 0,
                     upper = max(sen$doz))$value
    mxm <- optimize(fa, interval = interval, maximum = TRUE)
    ymid <- mxm[[2]]/2
    u <- try(uniroot(f = function(x) fa(x) - ymid,
                     interval = interval),
             silent = TRUE)
    if (class(u) == "try-error") {
        return(c(ec50 = NA, ev50 = NA, auc = int))
    }
    else {
        return(c(ec50 = u$root, ev50 = ymid, auc = int))
    }
288 289
}

290 291 292 293
# EC50 and AUC for each experimental unit.
ec <- ddply(res, .(ue), .fun = function(x) get_ec50(x$xfit, x$yfit))
str(ec)

294 295 296
# Proportion of not estimated EC50 and AUC.
cbind(AUC = sum(is.na(ec$auc)),
      EC50 = sum(is.na(ec$ev50)))/nrow(ec)
297 298 299 300 301 302 303 304 305 306

# Scatter plot matrix of estimated values.
splom(ec[, -1], type = c("p", "smooth"), col.line = 2)

# Correlation among variables.
cor(ec[, -1], use = "complete")

# Merge to pair `iso` and `fun`.
ec <- merge(ec, senu, by = "ue", all = TRUE)
str(ec)
307 308
```

309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351
```{r, fig.cap = cap, fig.show = "hold", echo = -(1:2), fig.height = 18}
cap <-
"Mean diameter as function of fungicide dose 5th root for each isolate.  Solid line is the natural cubic spline fitted for each fungicide in each isolate. Gray straight line indicates the EC~50~ on each curve."
cap <- fgn_("splines-fit", cap)

# View the results.
L <- split(ec, ec$iso)
xyplot(dm ~ dos^0.2 | iso,
       data = sen,
       groups = fun,
       cex = 0.4,
       strip = FALSE,
       as.table = TRUE,
       auto.key = list(columns = 3,
                       title = "In vitro fungicide",
                       cex.title = 1.1),
       ylim = c(0, NA),
       xlab = "In vitro fungicide dose 5th root",
       ylab = "Mean diameter (mm)") +
    as.layer(xyplot(yfit ~ xfit | iso,
                    groups = fun,
                    data = res,
                    type = "l")) +
    layer({
        with(L[[which.packet()]], {
            cl <- trellis.par.get()$superpose.symbol$col[as.integer(fun)]
            panel.segments(x0 = ec50,
                           y0 = ev50,
                           x1 = ec50,
                           y1 = 0,
                           col = "gray50")
            panel.segments(x0 = ec50,
                           y0 = ev50,
                           x1 = 0,
                           y1 = ev50,
                           col = "gray50")
            panel.points(x = ec50,
                         y = ev50,
                         pch = 19,
                         cex = 0.6,
                         col = cl)
        })
    })
352 353
```

354 355 356
****
# Analysis of Area Under Sensitivity Curve

357
```{r}
358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380
p <- xyplot(ec50 ~ fun | pop + hed,
            groups = tra,
            data = ec[!is.na(ec$ec50), ],
            type = c("p", "a"))
useOuterStrips(p)

p <- xyplot(auc ~ fun | pop + hed,
            groups = tra,
            data = ec[!is.na(ec$auc), ],
            type = c("p", "a"))
useOuterStrips(p)

p <- xyplot(auc ~ tra | pop + hed,
            groups = fun,
            data = ec[!is.na(ec$auc), ],
            type = c("p", "a"))
useOuterStrips(p)

p <- xyplot(auc ~ pop | tra + fun,
            groups = hed,
            data = ec[!is.na(ec$auc), ],
            type = c("p", "a"))
useOuterStrips(p)
381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398

#-----------------------------------------------------------------------
# Creates block and treatment cell factors.

ec$blk <- factor(as.integer(as.integer(substr(ec$plot, 0, 1)) > 2))
ec$cell <- with(ec, interaction(yr, blk, hed, tra, drop = TRUE))

# Number of isolates per cell combination.
ftable(xtabs(~pop + hed + tra, data = ec))/3

# ddply(ec,
#       ~yr + pop + tra + hed,
#       function(x) {
#           nlevels(droplevels(x$iso))
#       })

ec <- arrange(df = ec, yr, blk, hed, tra, iso, fun)
str(ec)
399 400
```

401 402
## 2015

Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
403 404 405
```{r}
#-----------------------------------------------------------------------

406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
ec15 <- subset(ec, yr == 2015)

# Mixed effects model.
m0 <- lmer(auc ~ blk + (1 | iso) + (pop + tra + hed + fun)^2,
           data = ec15,
           REML = FALSE)

# r <- residuals(m0)
# f <- fitted(m0)
# useOuterStrips(qqmath(~r | pop + tra, data = ec15))
# useOuterStrips(xyplot(r ~ f| pop + tra, data = ec15))

# Wald tests for the fixed effects.
anova(m0)

# A simpler model.
m1 <- update(m0, auc ~ blk + (1 | iso) + (pop + tra + fun))

# LRT between nested models.
anova(m1, m0)

# Parameter estimates.
summary(m1)

# Least squares means.
i <- c("pop", "tra", "fun")
L <- lapply(i,
       FUN = function(term){
           L <- LSmeans(m1, effect = term)
           rownames(L$K) <- L$grid[, 1]
           a <- apmc(L$K, m1, focus = term)
           names(a)[1] <- "level"
           a <- cbind(term = term, a)
           return(a)
       })
res <- ldply(L)
# str(res)

i <- c("Population", "In vivo fungicide", "In vitro fungicide")
```
```{r, fig.cap = cap, echo = -(1:2)}
cap <-
"Area under isolate sensitivity curve for levels of population, *in vivo* fungicide and *in vitro* fungicide. Pairs of means in a factor followed by the same letter are not statistically different at 5% significance level."
cap <- fgn_("auc-2015", cap)

resizePanels(
    segplot(level ~ lwr + upr | term,
            centers = fit,
            data = res,
            draw = FALSE,
            layout = c(1, NA),
            scales = list(y = list(relation = "free")),
            xlab = "Area under isolate sensitivity curve",
            ylab = "Levels of each factor",
            strip = strip.custom(factor.levels = i),
            cld = res$cld) +
    layer(panel.text(x = centers,
                     y = z,
                     labels = sprintf("%0.1f %s", centers, cld),
                     pos = 3)),
    h = sapply(L, nrow)
)
```
Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
469

470
## 2016
Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
471

472 473
```{r}
#-----------------------------------------------------------------------
Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
474

475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564
ec16 <- subset(ec, yr == 2016)

# Mixed effects model.
m0 <- lmer(auc ~ blk + (1 | iso) + (pop + tra + hed + fun)^2,
           data = ec16,
           REML = FALSE)

# r <- residuals(m0)
# f <- fitted(m0)
# useOuterStrips(qqmath(~r | pop + tra, data = ec15))
# useOuterStrips(xyplot(r ~ f| pop + tra, data = ec15))

# Wald tests for the fixed effects.
anova(m0)

# A simpler model.
m1 <- update(m0, auc ~ blk + (1 | iso) + pop * (tra + fun))

# LRT between nested models.
anova(m1, m0)

# Parameter estimates.
summary(m1)

# Least squares means.
res <- vector(mode = "list", length = 2)

L <- LSmeans(m1, effect = c("pop", "tra"))
g <- L$grid
L <- by(L$K, L$grid$pop, as.matrix)
L <- lapply(L, "rownames<-", levels(ec$tra))
L <- lapply(L, apmc, model = m1, focus = "tra")
res[[1]] <- ldply(L, .id = "pop")

L <- LSmeans(m1, effect = c("pop", "fun"))
g <- L$grid
L <- by(L$K, L$grid$pop, as.matrix)
L <- lapply(L, "rownames<-", levels(ec$fun))
L <- lapply(L, apmc, model = m1, focus = "fun")
res[[2]] <- ldply(L, .id = "pop")

L <- lapply(res,
            FUN = function(x) {
                x$by <- names(x)[2]
                names(x)[2] <- "term"
                return(x)
            })
res <- ldply(L)
res <- arrange(res, by, term, pop)

i <- c("In vivo fungicide", "In vitro fungicide")
p <- c(1, 2)
```
```{r, fig.cap = cap, echo = -(1:2)}
cap <-
"Area under isolate sensitivity curve for combination between population and *in vitro* fungicide (top) and population and *in vivo* fungicide (bottom). Pairs of means comparing fungicides (*in vivo* or *in vitro*) followed by the same letter are not statistically different at 5% significance level."
cap <- fgn_("auc-2016", cap)

resizePanels(
    segplot(term ~ lwr + upr | by,
            centers = fit,
            groups = pop,
            data = res,
            draw = FALSE,
            layout = c(1, NA),
            scales = list(y = list(relation = "free")),
            xlab = "Area under isolate sensitivity curve",
            ylab = "Levels of each factor",
            strip = strip.custom(factor.levels = i),
            key = list(title = "Population",
                       cex.title = 1.1,
                       columns = 2,
                       type = "b",
                       divide = 1,
                       lines = list(pch = p, lty = 1),
                       text = list(levels(res$pop))),
            cld = res$cld,
            panel = panel.groups.segplot,
            pch = p[as.integer(res$pop)],
            gap = 0.05) +
    layer(panel.text(x = centers[subscripts],
                     y = as.integer(z[subscripts]) +
                         centfac(groups[subscripts],
                                 space = gap),
                     labels = sprintf("%0.1f %s",
                                      centers[subscripts],
                                      cld[subscripts]),
                     pos = 3)),
    h = c(6, 8)
)
Walmes Marques Zeviani's avatar
Walmes Marques Zeviani committed
565 566
```

567 568 569 570 571 572 573 574 575 576 577 578
****
# Session information

```{r, echo=FALSE, results="hold"}
cat(format(Sys.time(), format = "%A, %d de %B de %Y, %H:%M"),
    "----------------------------------------", sep = "\n")
sessionInfo()
```

<!------------------------------------------- -->
[Paulo S. F. Lichtemberg]: http://lattes.cnpq.br/8132272273348880
[Walmes M. Zeviani]: http://www.leg.ufpr.br/doku.php/pessoais:walmes