Estimação da EC50 para isolados de Septoria
João Peterson F. Gardin e Walmes Zeviani
Download:
- Dados: septoriaEC50.txt.
- Script R: septoriaEC50.R.
Definições da sessão
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
## Instruções para instalação do wzRfun em:
## browseURL("https://github.com/walmes/wzRfun")
pkg <- c("lattice", "latticeExtra",
"nlme", "car",
"reshape", "plyr",
"wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## lattice latticeExtra nlme car reshape plyr wzRfun
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE##-----------------------------------------------------------------------------
## Configurações de cores da lattice.
colp <- brewer.pal(11, "Spectral")
colp <- colorRampPalette(colp, space="rgb")
mycol <- brewer.pal(8, "Dark2")
ps <- list(
box.dot=list(col="black", pch="|"),
box.rectangle=list(col="black", fill=c("gray70")),
box.umbrella=list(col="black", lty=1),
dot.symbol=list(col="black"),
dot.line=list(col="gray90", lty=1),
plot.symbol=list(col="black", cex=0.8),
plot.line=list(col="black", lty=1.2),
plot.polygon=list(col="gray80"),
superpose.line=list(col=mycol),
superpose.symbol=list(col=mycol),
superpose.polygon=list(col=mycol),
regions=list(col=colp(100)),
strip.background=list(col=c("gray90","gray50")),
strip.shingle=list(col=c("gray50","gray90"))
)
trellis.par.set(ps)
rm(list=c("colp", "mycol", "ps"))
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()## R version 3.4.0 (2017-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.80 plyr_1.8.4 reshape_0.8.6 car_2.1-4
## [5] nlme_3.1-131 latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-35
## [9] rmarkdown_1.5 knitr_1.15.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.10 compiler_3.4.0 nloptr_1.0.4 methods_3.4.0 tools_3.4.0
## [6] digest_0.6.12 lme4_1.1-12 evaluate_0.10 mgcv_1.8-17 Matrix_1.2-10
## [11] rpanel_1.1-3 yaml_2.1.14 parallel_3.4.0 SparseM_1.74 mvtnorm_1.0-5
## [16] stringr_1.2.0 MatrixModels_0.4-1 rprojroot_1.2 grid_3.4.0 nnet_7.3-12
## [21] survival_2.40-1 multcomp_1.4-6 TH.data_1.0-8 minqa_1.2.4 magrittr_1.5
## [26] codetools_0.2-15 backports_1.0.5 htmltools_0.3.6 MASS_7.3-47 splines_3.4.0
## [31] pbkrtest_0.4-6 quantreg_5.29 sandwich_2.3-4 stringi_1.1.2 doBy_4.5-15
## [36] zoo_1.7-14## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)Importação dos dados
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/septoriaEC50.txt"
path <- ifelse(Sys.info()["user"]=="walmes",
yes=basename(url), no=url)
## Importa os dados a partir do arquivo na web.
da <- read.table(file=path, header=TRUE, sep="\t",
colClasses=c("factor","numeric","integer")[c(1,1,2,3,1,3)])
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 1560 obs. of 6 variables:
## $ exper : Factor w/ 2 levels "I","II": 1 1 1 1 1 1 1 1 1 1 ...
## $ fung : Factor w/ 2 levels "Manzate","Tiofanato": 2 2 2 2 2 2 2 2 2 2 ...
## $ dose : num 0.1 0.1 0.1 0.1 0.1 1 1 1 1 1 ...
## $ rept : int 1 2 3 4 5 1 2 3 4 5 ...
## $ isol : Factor w/ 13 levels "475-1","475-2",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ germ100: int 95 94 95 96 97 92 90 93 92 94 ...## Tabela de frequência.
xtabs(~isol+dose, data=da)## dose
## isol 0 0.1 1 10 100 1000
## 475-1 20 20 20 20 20 20
## 475-2 20 20 20 20 20 20
## 475-3 20 20 20 20 20 20
## 475-4 20 20 20 20 20 20
## 475-5 20 20 20 20 20 20
## 475-6 20 20 20 20 20 20
## 475-7 20 20 20 20 20 20
## 475-8 20 20 20 20 20 20
## 475-9 20 20 20 20 20 20
## 476-1 20 20 20 20 20 20
## 476-2 20 20 20 20 20 20
## 476-3 20 20 20 20 20 20
## 476-4 20 20 20 20 20 20## Ordenar as linhas.
da <- arrange(da, exper, isol, rept, dose)
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(germ100~dose|isol, groups=fung, data=da, as.table=TRUE,
type=c("p","a"), jitter.x=TRUE, auto.key=TRUE,
scales=list(x=list(log=10)))
## Com log a dose zero desaparece pois log(0) não existe!
k <- 0.1
xyplot(germ100~(dose^k)|isol, groups=fung, data=da,
jitter.x=TRUE, auto.key=TRUE, as.table=TRUE,
type=c("p","a"))
## A escolha de qualquer potência é arbitraria. Será usada k=0.1 aqui
## porque apresenta uma disposição satisfatoriamente rugular entre as
## doses.
da$x <- da$dose^kEspecificação e estimação do modelo
O modelo considerado para germinação de cada isolado em função da dose para cada fungicida é representado por
\[ \begin{aligned} Y|x\, &\sim \text{Normal}(f(x), \sigma^2) \newline f(x)\, &= \theta_{a}\left(1-\left(\dfrac{x}{\theta_{l}}\right)^{\exp\{\theta_{c}\}}\right), \text{ se } x\leq \theta_{l}, \newline f(x)\, &= 0, \text{ se } x> \theta_{l}, \end{aligned} \]
em que \(Y\) representa a variável resposta germinação cujo valor médio é uma função \(f\) da dose \(x\) do fungicida descrita por um vetor de parâmetros desconhecidos e a variância é representada por \(\sigma^2\). No modelo não linear, os parâmetros têm o seguinte significado: * \(\theta_{a}\): germinação na ausência de fungicida (dose zero); * \(\theta_{l}\): dose de fungicida acima da qual a germinação é nula, por isso considera-se essa dose como sendo a limítrofe para o desenvolvimento; * \(\theta_{c}\): parâmetro de forma do modelo associado à curvatura.
A partir desse modelo, a dosse correpondente à uma germinação de 50% é obtida por: \[ \begin{aligned} EC_{50} &= \theta_{l} (1-50/\theta_{a})^{\exp\{-\theta_{c}\}}, \end{aligned} \]
que é a solução em \(x\) para a expressão do modelo quando \(f(x)=50\).
Após o ajuste do modelo, o erro-padrão, e por consequência, o intervalo de confiança para \(EC_{50}\) pode ser obtido via aplicação do método delta.
##-----------------------------------------------------------------------------
## Modelo não linear considerado para o caso.
L <- list(tha=80, thl=2.4, thc=1.3)
with(L,
curve(tha*(1-(x/thl)^exp(thc))*(x<=thl)+0, 0, 4,
ylim=c(0, 100),
xlab="x", ylab="f(x)"))
with(L,
{
d <- 0.2
ec50 <- thl*(1-50/tha)^exp(-thc)
abline(h=c(50, tha), v=c(ec50, thl), lty=2)
text(0, tha, label=expression(theta[a]), adj=c(-d,-d))
text(thl, 0, label=expression(theta[l]), adj=c(-d,-d))
text(ec50, 0, label=expression(EC[50]), adj=c(-d,-d))
})
## tha: germinação na dose 0;
## thl: dose acima da qual a germinação é nula;
## thc: fator relacionado à curvatura da função;
## O valor correspondente à uma germinação de 50% é dado por
## x_ec50 = thl*(1-50/tha)^exp(-thc)
##-----------------------------------------------------------------------------
## Aplicando o modelo à todos os casos, modelo de efeito fixo para os
## valores médios observados nas repetições.
db <- ddply(da, .(fung, exper, isol, x), summarise,
g100m=mean(germ100))
db$caso <- with(db, interaction(fung, isol))
xyplot(g100m~x|isol, groups=fung, data=db, as.table=TRUE)
n00 <- nlsList(g100m~tha*(1-(x/thl)^exp(thc))*(x<=thl)+0|caso,
start=c(tha=90, thl=2.1, thc=1.6),
data=db)
##-----------------------------------------------------------------------------
## Tabelas.
## Dos coeficientes.
coef(n00)## tha thl thc
## Manzate.475-1 82.81517 2.013078 2.242759
## Tiofanato.475-1 94.53497 2.144615 1.698787
## Manzate.475-2 93.48920 2.132101 1.367985
## Tiofanato.475-2 96.14151 2.131820 1.417472
## Manzate.475-3 82.25162 2.012664 1.969935
## Tiofanato.475-3 95.11726 2.217967 1.302250
## Manzate.475-4 85.09161 2.015707 1.691238
## Tiofanato.475-4 95.23327 2.183991 1.434912
## Manzate.475-5 86.55009 2.026815 1.659123
## Tiofanato.475-5 94.66575 2.210630 1.487843
## Manzate.475-6 89.64452 2.096351 1.093903
## Tiofanato.475-6 95.66668 2.229585 1.416346
## Manzate.475-7 85.33309 2.014153 1.705094
## Tiofanato.475-7 93.76357 2.139852 1.733826
## Manzate.475-8 92.74297 2.110189 1.318782
## Tiofanato.475-8 95.62572 2.157358 1.653316
## Manzate.475-9 94.27222 2.086466 1.765540
## Tiofanato.475-9 95.84372 2.158666 1.595519
## Manzate.476-1 92.52534 2.038571 1.962012
## Tiofanato.476-1 95.75385 2.170295 1.528655
## Manzate.476-2 91.05127 2.104647 1.035130
## Tiofanato.476-2 95.55963 2.019769 1.555194
## Manzate.476-3 95.15301 2.032306 1.703534
## Tiofanato.476-3 96.17787 2.022941 1.513091
## Manzate.476-4 93.29823 2.027958 1.949608
## Tiofanato.476-4 95.37629 2.066446 1.673697## Com os extremos do IC.
intervals(n00)## , , tha
##
## lower est. upper
## Manzate.475-1 80.83406 82.81517 84.79629
## Tiofanato.475-1 92.20976 94.53497 96.86018
## Manzate.475-2 90.76767 93.48920 96.21073
## Tiofanato.475-2 93.48694 96.14151 98.79607
## Manzate.475-3 80.14089 82.25162 84.36234
## Tiofanato.475-3 92.30464 95.11726 97.92988
## Manzate.475-4 82.75896 85.09161 87.42426
## Tiofanato.475-4 92.60183 95.23327 97.86470
## Manzate.475-5 84.18488 86.55009 88.91530
## Tiofanato.475-5 92.10271 94.66575 97.22879
## Manzate.475-6 86.54576 89.64452 92.74329
## Tiofanato.475-6 93.01061 95.66668 98.32275
## Manzate.475-7 83.01403 85.33309 87.65215
## Tiofanato.475-7 91.47180 93.76357 96.05534
## Manzate.475-8 89.95342 92.74297 95.53253
## Tiofanato.475-8 93.25447 95.62572 97.99698
## Manzate.475-9 92.00919 94.27222 96.53524
## Tiofanato.475-9 93.40966 95.84372 98.27778
## Manzate.476-1 90.40970 92.52534 94.64097
## Tiofanato.476-1 93.24142 95.75385 98.26628
## Manzate.476-2 87.87797 91.05127 94.22457
## Tiofanato.476-2 93.07902 95.55963 98.04025
## Manzate.476-3 92.83243 95.15301 97.47359
## Tiofanato.476-3 93.64638 96.17787 98.70936
## Manzate.476-4 91.17475 93.29823 95.42170
## Tiofanato.476-4 93.02604 95.37629 97.72654
##
## , , thl
##
## lower est. upper
## Manzate.475-1 2.001656 2.013078 2.024501
## Tiofanato.475-1 2.106007 2.144615 2.183223
## Manzate.475-2 2.093738 2.132101 2.170465
## Tiofanato.475-2 2.095180 2.131820 2.168460
## Manzate.475-3 1.998472 2.012664 2.026857
## Tiofanato.475-3 2.161375 2.217967 2.274560
## Manzate.475-4 1.997795 2.015707 2.033620
## Tiofanato.475-4 2.136005 2.183991 2.231977
## Manzate.475-5 2.007433 2.026815 2.046197
## Tiofanato.475-5 2.155377 2.210630 2.265883
## Manzate.475-6 2.057103 2.096351 2.135598
## Tiofanato.475-6 2.170131 2.229585 2.289040
## Manzate.475-7 1.996678 2.014153 2.031629
## Tiofanato.475-7 2.102031 2.139852 2.177673
## Manzate.475-8 2.074648 2.110189 2.145731
## Tiofanato.475-8 2.116089 2.157358 2.198628
## Manzate.475-9 2.060995 2.086466 2.111937
## Tiofanato.475-9 2.117277 2.158666 2.200054
## Manzate.476-1 2.022553 2.038571 2.054590
## Tiofanato.476-1 2.126021 2.170295 2.214568
## Manzate.476-2 2.063369 2.104647 2.145925
## Tiofanato.476-2 2.001352 2.019769 2.038187
## Manzate.476-3 2.014683 2.032306 2.049929
## Tiofanato.476-3 2.003632 2.022941 2.042250
## Manzate.476-4 2.013455 2.027958 2.042462
## Tiofanato.476-4 2.043858 2.066446 2.089034
##
## , , thc
##
## lower est. upper
## Manzate.475-1 2.0187474 2.242759 2.466770
## Tiofanato.475-1 1.5139656 1.698787 1.883608
## Manzate.475-2 1.2154607 1.367985 1.520509
## Tiofanato.475-2 1.2662216 1.417472 1.568722
## Manzate.475-3 1.7965410 1.969935 2.143329
## Tiofanato.475-3 1.1331886 1.302250 1.471311
## Manzate.475-4 1.5453459 1.691238 1.837131
## Tiofanato.475-4 1.2646439 1.434912 1.605180
## Manzate.475-5 1.5131713 1.659123 1.805075
## Tiofanato.475-5 1.3017423 1.487843 1.673943
## Manzate.475-6 0.9566831 1.093903 1.231122
## Tiofanato.475-6 1.2335730 1.416346 1.599119
## Manzate.475-7 1.5594662 1.705094 1.850721
## Tiofanato.475-7 1.5449628 1.733826 1.922690
## Manzate.475-8 1.1736204 1.318782 1.463943
## Tiofanato.475-8 1.4706768 1.653316 1.835955
## Manzate.475-9 1.5994373 1.765540 1.931643
## Tiofanato.475-9 1.4195621 1.595519 1.771476
## Manzate.476-1 1.7941375 1.962012 2.129886
## Tiofanato.476-1 1.3551198 1.528655 1.702191
## Manzate.476-2 0.9007962 1.035130 1.169464
## Tiofanato.476-2 1.4295036 1.555194 1.680885
## Manzate.476-3 1.5664191 1.703534 1.840649
## Tiofanato.476-3 1.3886622 1.513091 1.637519
## Manzate.476-4 1.7906756 1.949608 2.108540
## Tiofanato.476-4 1.5260339 1.673697 1.821359## Gráfico padrão dos IC.
plot(intervals(n00), as.table=TRUE)
##-----------------------------------------------------------------------------
## Calculando o valor da EC50 com intervalo de confiança baseado no
## método delta.
ec50 <- sapply(n00,
function(i){
dm <- deltaMethod(i,
g="thl*(1-50/tha)^exp(-thc)")
ic <- with(dm, c(lower=Estimate-1.96*SE,
est.=Estimate,
upper=Estimate+1.96*SE))
return(ic)
})
ec50 <- t(ec50)
##-----------------------------------------------------------------------------
## Intervalos de confiança dos demais parâmetros.
## IC para os parâmetros do modelo.
ic <- intervals(n00)
##-----------------------------------------------------------------------------
## Juntar EC50 aos outros, todos no mesmo array.
aux <- array(c(ic, ec50), dim=dim(ic)+c(0,0,1))
dn <- dimnames(ic)
dn[[3]] <- c(dn[[3]], "ec50")
dimnames(aux) <- dn
ic <- aux
ic <- adply(ic, .margins=c(1, 3))
names(ic)[1:2] <- c("caso", "param")
aux <- as.data.frame(stringsAsFactors=FALSE,
do.call(rbind, strsplit(as.character(ic$caso), "\\.")))
names(aux) <- c("fung", "isol")
ic <- cbind(ic, equallevels(aux, da))##-----------------------------------------------------------------------------
## Gráfico de segmentos para representar o IC 95%.
ic <- arrange(ic, param, isol, fung)
l <- levels(ic$fung)
key <- list(title="Fungicidas", cex.title=1.1,
columns=2, type="o", divide=1,
text=list(l),
lines=list(pch=1:length(l)))
segplot(isol~lower+upper|param, centers=est.,
data=ic, draw=FALSE, as.table=TRUE,
ylab="Isolados", xlab="Estimativa do parâmetro",
scales=list(x=list(relation="free")), key=key,
groups=fung, gap=0.12, pch=as.integer(ic$fung),
panel=panel.groups.segplot)
##-----------------------------------------------------------------------------
## Gráfico só do EC50 e na escala original.
ec50 <- droplevels(subset(ic, param=="ec50"))
i <- c("lower", "est.", "upper")
ec50[,i] <- ec50[,i]^(1/k)
segplot(isol~lower+upper, centers=est.,
data=ec50, draw=FALSE, as.table=TRUE,
ylab="Isolados",
xlab=expression(Estimativa~da~EC[50]),
scales=list(x=list(relation="free")), key=key,
groups=fung, gap=0.1, pch=as.integer(ic$fung),
panel=panel.groups.segplot)
##-----------------------------------------------------------------------------
## Gráfico dos valores preditos pelo modelo regressão não linear ajustado.
pred <- data.frame(x=seq(0, 2.5, length.out=50))
pdto <- lapply(n00,
function(i){
pred <- cbind(pred,
fit=predict(i, newdata=pred))
return(pred)
})
pdto <- ldply(pdto)
names(pdto)[1] <- "caso"
aux <- as.data.frame(stringsAsFactors=FALSE,
do.call(rbind, strsplit(as.character(pdto$caso), "\\.")))
names(aux) <- c("fung", "isol")
pdto <- cbind(pdto, equallevels(aux, da))
## xyplot(g100m~x|isol, groups=fung, data=db)
## xyplot(fit~x|isol, groups=fung, data=pdto, type="l")##-----------------------------------------------------------------------------
## Sobrepondo os gráficos e adicionando uma vertical da EC50.
L <- droplevels(
subset(ic, param=="ec50", select=c("isol", "fung", "est.")))
L <- dlply(L, .(isol), .fun=droplevels)
d <- c(0, 0.01, 1, 10, 100, 1000)
scales <- list(x=list(at=d^k, labels=d))
xyplot(g100m~x|isol, groups=fung, data=db,
layout=c(3,NA), as.table=TRUE,
## xlab=expression(Dose~de~fungicida^0.1),
xlab=expression(Dose~de~fungicida), scales=scales,
ylab="Germinação observada (%)",
xlim=extendrange(pred$x))+
as.layer(xyplot(fit~x|isol, groups=fung, data=pdto, type="l",
panel=function(...){
panel.xyplot(...)
wp <- which.packet()
col <- trellis.par.get(
"superpose.line")$col[1:nlevels(pdto$fung)]
panel.abline(v=L[[wp]]$est., lty=2,
col=col)
}))
##-----------------------------------------------------------------------------
## Tabela com os valores das estimativas e R2.
## R2 como o quadrado da correlação entre observados e ajustados.
## m0 <- lm(dist~speed, cars)
## r2 <- cor(x=fitted(m0)+residuals(m0),
## fitted(m0))^2
## r2; summary(m0)$r.squared
tab <- lapply(n00,
function(i){
coef <- coef(i)
r2 <- cor(x=fitted(i)+residuals(i),
fitted(i))^2
res <- as.data.frame(as.list(c(coef, R2=r2)))
return(res)
})
tab <- ldply(tab)
names(tab)[1] <- "caso"
aux <- as.data.frame(stringsAsFactors=FALSE,
do.call(rbind, strsplit(as.character(tab$caso), "\\.")))
names(aux) <- c("fung", "isol")
tab <- cbind(equallevels(aux, da),
subset(tab, select=-c(caso)))
tab## fung isol tha thl thc R2
## 1 Manzate 475-1 82.81517 2.013078 2.242759 0.9915465
## 2 Tiofanato 475-1 94.53497 2.144615 1.698787 0.9981540
## 3 Manzate 475-2 93.48920 2.132101 1.367985 0.9919622
## 4 Tiofanato 475-2 96.14151 2.131820 1.417472 0.9987158
## 5 Manzate 475-3 82.25162 2.012664 1.969935 0.9831894
## 6 Tiofanato 475-3 95.11726 2.217967 1.302250 0.9956300
## 7 Manzate 475-4 85.09161 2.015707 1.691238 0.9861063
## 8 Tiofanato 475-4 95.23327 2.183991 1.434912 0.9944946
## 9 Manzate 475-5 86.55009 2.026815 1.659123 0.9919060
## 10 Tiofanato 475-5 94.66575 2.210630 1.487843 0.9962558
## 11 Manzate 475-6 89.64452 2.096351 1.093903 0.9734795
## 12 Tiofanato 475-6 95.66668 2.229585 1.416346 0.9984303
## 13 Manzate 475-7 85.33309 2.014153 1.705094 0.9757344
## 14 Tiofanato 475-7 93.76357 2.139852 1.733826 0.9977691
## 15 Manzate 475-8 92.74297 2.110189 1.318782 0.9889762
## 16 Tiofanato 475-8 95.62572 2.157358 1.653316 0.9987181
## 17 Manzate 475-9 94.27222 2.086466 1.765540 0.9990232
## 18 Tiofanato 475-9 95.84372 2.158666 1.595519 0.9987143
## 19 Manzate 476-1 92.52534 2.038571 1.962012 0.9909503
## 20 Tiofanato 476-1 95.75385 2.170295 1.528655 0.9977367
## 21 Manzate 476-2 91.05127 2.104647 1.035130 0.9780223
## 22 Tiofanato 476-2 95.55963 2.019769 1.555194 0.9988407
## 23 Manzate 476-3 95.15301 2.032306 1.703534 0.9989225
## 24 Tiofanato 476-3 96.17787 2.022941 1.513091 0.9984596
## 25 Manzate 476-4 93.29823 2.027958 1.949608 0.9950014
## 26 Tiofanato 476-4 95.37629 2.066446 1.673697 0.9987166