Habilidade competitiva
Aplicações continuas de mesmos fungicidas exercem presão de seleção sob as populações fúngicas presentes. A frequencias dos isolados integrantes da população vem sendo modificadas. Isolados com caracteristicas que permite sobreviver e se reproduzir nessa condição (“resistentes”) aumentam sua proporção.
O objetivo de este trabalho é verficar se retirando a pressão de seleção (uso continuo de um mesmo fungicida) a frequencia de isolados original da populaçao pode ser reestabelecida.
Para isso foram feitos tres experimentos (por duplicado) inoculando (concentração padrao )frutos de pessegos com um isolado resistente (SP09) e um isolado sensivel (PR09) ao longo de 5 transferencias.
Experimento 1 e 2: as unidades experimentales foram frutos de pessego frescos (variedade dorado). No exp1 foi avaliado unidades formadoras de colonia e no exp2 foi medida a germinação dos esporos dos conidios em meio de cultura com BDA com ou sem fungicida.
Experimento 3: as unidades experimentales foram pedaços de pessego enlatados. A variavel medida foi germinação dos esporos dos conidios em meio de cultura com BDA com ou sem fungicida.
Um objetivo secundario do trabalho foi testar se as diferentes metodologias conduzem a conclusões semelhantes.
Session set
pkg <- c("MASS",
"ggplot2","knitr","reshape2",
"dplyr","grid","gridExtra",
"Rmisc", "multcomp", "hnp")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## MASS ggplot2 knitr reshape2 dplyr grid gridExtra
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## Rmisc multcomp hnp
## TRUE TRUE TRUE#setwd("/home/epi/Dropbox/MyR/Análise otros/ISA/tese") # lab
setwd("/media/DATA/Dropbox/MyR/Análise otros/ISA/tese/aa habilidade")# DELL
opts_chunk$set(
warnings=FALSE,
cache=FALSE,
tidy=FALSE,
fig.width=8,
fig.height=6,
fig.align="center",
dpi=150)
options(width=90)
source('/media/DATA/Dropbox/MyR/Sources/boxplot.with.outlier.label.r', encoding='UTF-8')
#source("https://raw.githubusercontent.com/talgalili/R-code-snippets/master/boxplot.with.outlier.label.r")
trans.arcsine <- function(x){ asin(sign(x) * sqrt(abs(x))) }Unidades formadoras de colonias.
Importing dataset
dat1 = read.csv("fit4ufc.csv", colClasses = c(rep("integer", 15), rep("NULL", 4)),
dec=",", header=T, sep="\t", check.names=FALSE, na.strings=".")Manipulating & summarizing
long1 = melt(select(dat1, exp, rep, subrep, bda0, bda1, bda2, bda3,bda4,bda5),
id.vars = c("exp", "rep","subrep"),
variable.name="tran_bda", value.name="uf_bda") #; head(long1)
long2 = melt(select(dat1, exp, rep, subrep, f0, f1, f2, f3,f4,f5),
id.vars = c("exp", "rep","subrep"),
variable.name="tran_f", value.name="uf_f") #; head(long2)
long = cbind(long1, long2[,4:5])
long$tran_bda = gsub("bda", "", paste(long$tran_bda))
long$tran_f = gsub("f", "", paste(long$tran_f))
long = long[complete.cases(long),]
long = select(long, exp, rep, subrep, transf=tran_bda, uf_bda, uf_f)
long$resp = long$uf_f / long$uf_bda Explore data
da1 = aggregate(resp ~ exp * rep * transf, data = long, mean)
da1 = filter(da1, resp<1)
boxplot.with.outlier.label(da1$resp ~ da1$transf, row.names(da1), ylim = c(0,1.1))
head(da1); #str(da1)## exp rep transf resp
## 1 2 1 0 0.1947587
## 2 2 3 0 0.4555749
## 3 1 4 0 0.6724953
## 4 2 4 0 0.2805586
## 5 2 5 0 0.8136364
## 6 1 6 0 0.6155348for (i in 1:3) da1[,i] <- as.factor(da1[,i]) Opcoes de graficos
(sumda1 <- summarySEwithin(da1, measurevar="resp", withinvars="transf",na.rm=FALSE, conf.interval=.95))## transf N resp sd se ci
## 1 0 7 0.5626870 0.29062089 0.10984437 0.26877950
## 2 1 7 0.4633260 0.23830614 0.09007125 0.22039642
## 3 2 10 0.8654939 0.10870624 0.03437593 0.07776376
## 4 3 9 0.7629645 0.15910910 0.05303637 0.12230208
## 5 4 10 0.6248907 0.21184469 0.06699117 0.15154456
## 6 5 8 0.7628438 0.09562574 0.03380880 0.07994512p0<-ggplot(da1, aes(transf, resp)) + theme_bw() +
geom_point(aes(shape=factor(exp)), size=3) +
scale_shape_manual(values = c(1, 19), name = "Experimento") +
xlab("Transferencias") + ylab("Proporcao de resistentes") +
theme(legend.position = 'none') +
ggtitle("Observacoes por experimento")
p1<-ggplot(sumda1, aes(transf, resp)) + theme_bw() +
geom_point(size=3) +
geom_errorbar(aes(x = transf, ymin = resp-se, ymax = resp+se), colour="black", width=.1) +
geom_line() +
scale_y_continuous(limits= c(0.2,1)) +
xlab("Transferencias") + ylab("Proporcao de resistentes") +
ggtitle("Error bars = erro padrao")
p2<-ggplot(sumda1, aes(transf, resp)) + theme_bw() +
geom_point( size=3) +
geom_errorbar(aes(x = transf, ymin = resp-ci, ymax = resp+ci), colour="black", width=.1) +
scale_y_continuous(limits= c(0.2,1)) +
xlab("Transferencias") + ylab("Proporcao de resistentes") +
ggtitle("Error bars = confidence interval")
grid.arrange(p0,p1,p2, ncol=3)
ANOVA
mod1 = glm(resp ~ exp + transf, data=da1, family = quasibinomial)
#anova(mod.glm, test="Chisq") # binomial model
anova(mod1, test="F") # quasi-binomial model## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: resp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 50 11.2284
## exp 1 0.5620 49 10.6664 4.1474 0.0477438 *
## transf 5 4.2432 44 6.4232 6.2631 0.0001811 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1hnp1 = hnp(mod1, print.on=TRUE, plot=FALSE) ; plot(hnp1)## Quasi-binomial model
(glht1 = summary(glht(mod1, linfct=mcp(transf="Tukey")), test=adjusted(type="bonferroni")))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = resp ~ exp + transf, family = quasibinomial, data = da1)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 -0.21082 0.41060 -0.513 1.000000
## 2 - 0 == 0 1.67621 0.44564 3.761 0.002535 **
## 3 - 0 == 0 1.04819 0.41180 2.545 0.163728
## 4 - 0 == 0 0.35808 0.37583 0.953 1.000000
## 5 - 0 == 0 1.07783 0.42754 2.521 0.175532
## 2 - 1 == 0 1.88703 0.44602 4.231 0.000349 ***
## 3 - 1 == 0 1.25901 0.40396 3.117 0.027433 *
## 4 - 1 == 0 0.56890 0.37302 1.525 1.000000
## 5 - 1 == 0 1.28865 0.41600 3.098 0.029256 *
## 3 - 2 == 0 -0.62803 0.44956 -1.397 1.000000
## 4 - 2 == 0 -1.31814 0.41914 -3.145 0.024926 *
## 5 - 2 == 0 -0.59838 0.46242 -1.294 1.000000
## 4 - 3 == 0 -0.69011 0.37800 -1.826 1.000000
## 5 - 3 == 0 0.02964 0.42260 0.070 1.000000
## 5 - 4 == 0 0.71975 0.39272 1.833 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)let1 <- cld(glht1) ; let1 <- let1$mcletters$LettersGraficos
#a = with(da1, adjust(resp, by=transf, model=mod1))
(p_uf4 = ggplot(sumda1, aes(transf, resp)) + theme_bw() +
geom_point(size=3) +
geom_errorbar(aes(x = transf, ymin = resp-se, ymax = resp+se), colour="black", width=.1) +
geom_line() +
scale_y_continuous(limits= c(0.2,1)) +
xlab("Transferências") + ylab("Proporcão de isolados resistentes") +
ggtitle("UFC 4") +
annotate("text", x = 1:6, y = 1, label = let1)
)
Germinação (5)
dat2 = read.csv("fit5ger.csv",
#colClasses = c(rep("integer", 15), rep("NULL", 4)),
dec=",", header=T, sep="\t", check.names=FALSE, na.strings=".")
head(dat2)[,1:15]## exp rep subrep bda0 f0 bda1 f1 bda2 f2 bda3 f3 bda4 f4 bda5 f5
## 1 1 1 1 84 45 82 87 100 84 100 54 95 41 98 62
## 2 1 1 2 NA NA 86 75 100 91 100 87 95 54 92 35
## 3 1 1 3 91 40 86 70 97 85 100 74 95 45 100 50
## 4 1 2 1 87 25 76 80 100 87 100 85 87 38 94 55
## 5 1 2 2 74 43 83 65 100 89 100 88 82 38 100 61
## 6 1 2 3 84 35 74 61 100 87 91 87 84 23 96 75Manipulating & summarizing
long_bda2 = melt(select(dat2, exp, rep, subrep, bda0, bda1, bda2, bda3,bda4,bda5),
id.vars = c("exp", "rep","subrep"),
variable.name="tran_bda", value.name="ger_bda") #; head(long_bda2)
long_f2 = melt(select(dat2, exp, rep, subrep, f0, f1, f2, f3,f4,f5),
id.vars = c("exp", "rep","subrep"),
variable.name="tran_f", value.name="ger_f") #; head(long_f2)
long2 = cbind(long_bda2, long_f2[,4:5])
long2$tran_bda = gsub("bda", "", paste(long2$tran_bda))
long2$tran_f = gsub("f", "", paste(long2$tran_f))
long2 = long2[complete.cases(long2),]
long2 = select(long2, exp, rep, subrep, transf=tran_bda, ger_bda, ger_f)
long2$resp = long2$ger_f / long2$ger_bda Explore data
da2 = aggregate(resp ~ exp * rep * transf, data = long2, mean)
da2 = filter(da2, resp<1)
head(da2); #str(da2)## exp rep transf resp
## 1 1 1 0 0.4876374
## 2 2 1 0 0.6494900
## 3 1 2 0 0.4283680
## 4 2 2 0 0.7592593
## 5 1 3 0 0.4121087
## 6 2 3 0 0.7110656for (i in 1:3) da2[,i] <- as.factor(da2[,i])ANOVA
boxplot.with.outlier.label(da2$resp ~ da2$transf, row.names(da2))
mod2 = glm(resp ~ exp + transf, data=da2, family = quasibinomial)
#anova(mod.glm, test="Chisq") # binomial model
anova(mod2, test="F") # quasi-binomial model## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: resp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 61 7.6514
## exp 1 0.2356 60 7.4158 3.3243 0.0737 .
## transf 5 3.7963 55 3.6195 10.7113 3.212e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1hnp2 = hnp(mod2, print.on=TRUE, plot=FALSE) ; plot(hnp2)## Quasi-binomial model
(glht2 = summary(glht(mod2, linfct=mcp(transf="Tukey")), test=adjusted(type="bonferroni")))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = resp ~ exp + transf, family = quasibinomial, data = da2)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 1.74553 0.34412 5.072 5.89e-06 ***
## 2 - 0 == 0 1.25099 0.31183 4.012 0.000904 ***
## 3 - 0 == 0 1.20724 0.29749 4.058 0.000742 ***
## 4 - 0 == 0 0.21714 0.27214 0.798 1.000000
## 5 - 0 == 0 0.48456 0.27615 1.755 1.000000
## 2 - 1 == 0 -0.49454 0.34348 -1.440 1.000000
## 3 - 1 == 0 -0.53829 0.33091 -1.627 1.000000
## 4 - 1 == 0 -1.52839 0.30845 -4.955 1.08e-05 ***
## 5 - 1 == 0 -1.26097 0.31196 -4.042 0.000795 ***
## 3 - 2 == 0 -0.04375 0.29688 -0.147 1.000000
## 4 - 2 == 0 -1.03385 0.27188 -3.803 0.002148 **
## 5 - 2 == 0 -0.76642 0.27577 -2.779 0.081735 .
## 4 - 3 == 0 -0.99010 0.25540 -3.877 0.001589 **
## 5 - 3 == 0 -0.72267 0.25963 -2.783 0.080677 .
## 5 - 4 == 0 0.26742 0.23021 1.162 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)let2 <- cld(glht2) ; let2 <- let2$mcletters$LettersGrafico
(sumda2 <- summarySEwithin(da2, measurevar="resp", withinvars="transf",na.rm=FALSE, conf.interval=.95))## transf N resp sd se ci
## 1 0 6 0.5746548 0.16517634 0.06743296 0.17334194
## 2 1 10 0.8847246 0.07702567 0.02435765 0.05510084
## 3 2 10 0.8194327 0.06268197 0.01982178 0.04483998
## 4 3 12 0.8177941 0.11857000 0.03422821 0.07533578
## 5 4 12 0.6263508 0.19274806 0.05564157 0.12246628
## 6 5 12 0.6861935 0.12025940 0.03471590 0.07640917# a = with(da2, adjust(resp, by=transf, model=mod1) )
(p_ger5 = ggplot(sumda2, aes(transf, resp)) + theme_bw() +
geom_point(size=3) +
geom_errorbar(aes(x = transf, ymin = resp-se, ymax = resp+se), colour="black", width=.1) +
geom_line() +
scale_y_continuous(limits= c(0.2,1)) +
xlab("Transferências") + ylab("Proporcão de isolados resistentes") +
ggtitle("Germinaçao 5") +
annotate("text", x = 1:6, y = 1, label = let2)
)
Germinação (6)
dat3 = read.csv("fit6ger.csv", dec=",", header=T, sep="\t", check.names=FALSE, na.strings=".")
head(dat3)[,1:15]## exp rep subrep bda0 f0 bda1 f1 bda2 f2 bda3 f3 bda4 f4 bda5 f5
## 1 1 1 1 98 53 95 54 97 64 100 75 100 78 97 72
## 2 1 1 2 96 48 95 50 97 48 100 84 100 79 99 77
## 3 1 1 3 92 42 93 56 98 52 99 78 100 75 100 78
## 4 1 2 1 100 35 99 75 98 61 100 76 99 75 100 63
## 5 1 2 2 96 41 97 74 100 81 100 65 98 73 99 58
## 6 1 2 3 95 52 98 55 95 53 100 77 96 64 100 59Manipulating & summarizing
long1 = melt(select(dat3, exp, rep, subrep, bda0, bda1, bda2, bda3,bda4,bda5),
id.vars = c("exp", "rep","subrep"),
variable.name="tran_bda", value.name="ger_bda") # ; head(long1)
long2 = melt(select(dat3, exp, rep, subrep, f0, f1, f2, f3,f4,f5),
id.vars = c("exp", "rep","subrep"),
variable.name="tran_f", value.name="ger_f") # ; head(long2)
long3 = cbind(long1, long2[,4:5])
long3$tran_bda = gsub("bda", "", paste(long3$tran_bda))
long3$tran_f = gsub("f", "", paste(long3$tran_f))
long3 = long3[complete.cases(long3),]
long3 = select(long3, exp, rep, subrep, transf=tran_bda, ger_bda, ger_f)
long3$resp = long3$ger_f / long3$ger_bda Explore data
da3 = aggregate(resp ~ exp * rep * transf, data = long, mean)
da3 = filter(da3, resp<1)
head(da3); #str(da3)## exp rep transf resp
## 1 2 1 0 0.1947587
## 2 2 3 0 0.4555749
## 3 1 4 0 0.6724953
## 4 2 4 0 0.2805586
## 5 2 5 0 0.8136364
## 6 1 6 0 0.6155348for (i in 1:3) da3[,i] <- as.factor(da3[,i])
xtabs(~transf + exp, data=da3)## exp
## transf 1 2
## 0 2 5
## 1 5 2
## 2 4 6
## 3 5 4
## 4 5 5
## 5 5 3boxplot.with.outlier.label(da3$resp ~ da3$transf, row.names(da3))
ANOVA
mod3 = glm(resp ~ exp + transf, data=da3, family = quasibinomial)
#anova(mod.glm, test="Chisq") # binomial model
anova(mod3, test="F") # quasi-binomial model## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: resp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 50 11.2284
## exp 1 0.5620 49 10.6664 4.1474 0.0477438 *
## transf 5 4.2432 44 6.4232 6.2631 0.0001811 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1hnp3 = hnp(mod3, print.on=TRUE, plot=FALSE) ; plot(hnp3)## Quasi-binomial model
(glht3 = summary(glht(mod3, linfct=mcp(transf="Tukey")), test=adjusted(type="bonferroni")))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = resp ~ exp + transf, family = quasibinomial, data = da3)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 -0.21082 0.41060 -0.513 1.000000
## 2 - 0 == 0 1.67621 0.44564 3.761 0.002535 **
## 3 - 0 == 0 1.04819 0.41180 2.545 0.163728
## 4 - 0 == 0 0.35808 0.37583 0.953 1.000000
## 5 - 0 == 0 1.07783 0.42754 2.521 0.175532
## 2 - 1 == 0 1.88703 0.44602 4.231 0.000349 ***
## 3 - 1 == 0 1.25901 0.40396 3.117 0.027433 *
## 4 - 1 == 0 0.56890 0.37302 1.525 1.000000
## 5 - 1 == 0 1.28865 0.41600 3.098 0.029256 *
## 3 - 2 == 0 -0.62803 0.44956 -1.397 1.000000
## 4 - 2 == 0 -1.31814 0.41914 -3.145 0.024926 *
## 5 - 2 == 0 -0.59838 0.46242 -1.294 1.000000
## 4 - 3 == 0 -0.69011 0.37800 -1.826 1.000000
## 5 - 3 == 0 0.02964 0.42260 0.070 1.000000
## 5 - 4 == 0 0.71975 0.39272 1.833 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)let3 <- cld(glht3) ; let3 <- let3$mcletters$LettersGrafico
(sumda3 <- summarySEwithin(da3, measurevar="resp", withinvars="transf",na.rm=FALSE, conf.interval=.95))## transf N resp sd se ci
## 1 0 7 0.5626870 0.29062089 0.10984437 0.26877950
## 2 1 7 0.4633260 0.23830614 0.09007125 0.22039642
## 3 2 10 0.8654939 0.10870624 0.03437593 0.07776376
## 4 3 9 0.7629645 0.15910910 0.05303637 0.12230208
## 5 4 10 0.6248907 0.21184469 0.06699117 0.15154456
## 6 5 8 0.7628438 0.09562574 0.03380880 0.07994512(p_ger6 = ggplot(sumda3, aes(transf, resp)) + theme_bw() +
geom_point(size=3) +
geom_errorbar(aes(x = transf, ymin = resp-se, ymax = resp+se), colour="black", width=.1) +
geom_line() +
scale_y_continuous(limits= c(0.2,1)) +
xlab("Transferências") + ylab("Proporcão de isolados resistentes") +
ggtitle("Germinaçao 6") +
annotate("text", x = 1:6, y = 1, label = let3)
)
Todos os graficos
grid.arrange(p_uf4,p_ger5,p_ger6, ncol=3)