##Analysis of Oocyst Prevalence### ##Take a look at raw data means & ci's## > aggregate( formula = O_Prev~Treatment+Block, + data = oocysts2, + FUN = mean ) Treatment Block O_Prev 1 0 1 0.4782685 2 1 1 0.7075087 3 0 2 0.4021973 4 1 2 0.6061574 library(lsr) aggregate( O_Prev~Treatment+Block, oocysts2, ciMean ) Treatment Block O_Prev.1 O_Prev.2 1 0 1 0.4554603 0.5010768 2 1 1 0.6828037 0.7322137 3 0 2 0.3748657 0.4295289 4 1 2 0.5815248 0.6307900 library(plyr) ###reshape dataset for sd's and se's### oplyr <- ddply(oocysts2, c("Block", "Treatment","Day"), summarise, + N = length(O_Prev), + mean = mean(O_Prev), + sd = sd(O_Prev), + se = sd / sqrt(N) + ) oplyr > oplyr2 <- ddply(oocysts2, c("Block", "Treatment"), summarise, + N = length(O_Prev), + mean = mean(O_Prev), + sd = sd(O_Prev), + se = sd / sqrt(N) + ) > oplyr2 ##make barplot## p <- ggplot(oplyr2, aes(fill=factor(Block), y=mean, x=Treatment)) p + geom_bar(position="dodge", stat="identity") dodge <- position_dodge(width=0.9) p + geom_bar(position=dodge,stat="identity") + geom_errorbar(limits, position=dodge, width=0.25) Call: glm(formula = cbind(infected, not.infected) ~ day + treatment + block, family = binomial(link = "logit"), data = obycup) Deviance Residuals: Min 1Q Median 3Q Max -3.1891 -0.3631 0.0035 0.3836 1.9426 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.51122 0.26704 1.914 0.05557 . day -0.02722 0.02847 -0.956 0.33899 treatment 0.81063 0.09705 8.353 < 2e-16 *** block -0.31152 0.09723 -3.204 0.00136 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 166.516 on 209 degrees of freedom Residual deviance: 86.349 on 206 degrees of freedom AIC: 629.19 Number of Fisher Scoring iterations: 3 ##"fortify" data to get residuals## oglm.5fort<-fortify(oglm.5) head(oglm.5fort) ##Plot residuals v. fitted## oresid<-ggplot(aes(x=.fitted,y=.resid),data=oglm.5fort)+geom_point()+geom_hline(yintercept=0)+geom_smooth(se=FALSE)+labs(x="Fitted Values",y="Residuals") oresid ##QQ Plot## ggplot(oglm.5fort,aes(sample=.stdresid))+stat_qq()+geom_abline() ##does not look good. #make new dataset with prevalence at cup level# #EU=elemental unit=Cup# oocysts$EU<-as.factor(paste(oocysts$Day,oocysts$Cup,oocysts$Treatment,oocysts$Block,sep="-")) oocysts2<-aggregate(cbind(Block,Treatment,Day,Cup,O_Prev,O_Intensity)~EU,data=oocysts,mean) oocysts2<-oocysts2[with(oocysts2,order(Block,Treatment,Cup,Day)),] head(oocysts2) EU Block Treatment Day Cup O_Prev 36 5-1-0-1 1 0 5 1 0.5000000 72 6-1-0-1 1 0 6 1 0.5000000 107 7-1-0-1 1 0 7 1 0.5000000 142 8-1-0-1 1 0 8 1 0.3750000 177 9-1-0-1 1 0 9 1 0.5000000 1 10-1-0-1 1 0 10 1 0.5454545 str(oocysts2) 'data.frame': 211 obs. of 6 variables: $ EU : Factor w/ 211 levels "10-1-0-1","10-1-0-2",..: 36 72 107 142 177 1 42 78 113 148 ... $ Block : num 1 1 1 1 1 1 1 1 1 1 ... $ Treatment: num 0 0 0 0 0 0 0 0 0 0 ... $ Day : num 5 6 7 8 9 10 5 6 7 8 ... $ Cup : num 1 1 1 1 1 1 2 2 2 2 ... $ O_Prev : num 0.5 0.5 0.5 0.375 0.5 ... ##Make initial boxplot## library(ggplot2) library(lattice) o<-ggplot(oocysts2,aes(y=O_Prev,x=Day)) o + geom_boxplot()+facet_grid(Treatment~Block) #violin plots# o + geom_violin()+facet_grid(Treatment~Block) ##calculating means and variance## mean() var() aggregate(O_Prev~Treatment,data=oocysts2,var) Treatment O_Prev 1 0 0.009708578 2 1 0.010753265 shapiro.test(oocysts2$O_Prev) Shapiro-Wilk normality test data: oocysts2$O_Prev W = 0.9758, p-value = 0.001055 shapiro.test(oocysts2$O_Prev[oocysts2$Treatment=="0"]) Shapiro-Wilk normality test data: oocysts2$O_Prev[oocysts2$Treatment == "0"] W = 0.9284, p-value = 3.143e-05 shapiro.test(oocysts2$O_Prev[oocysts2$Treatment=="1"]) Shapiro-Wilk normality test data: oocysts2$O_Prev[oocysts2$Treatment == "1"] W = 0.9505, p-value = 0.0005163 ###DATA ARE NOT NORMALLY DISTRIBUTED### ###SO, use Ansari-Bradley Test### ansari.test(O_Prev~Treatment,oocysts2) Ansari-Bradley test data: O_Prev by Treatment AB = 5601.5, p-value = 0.5879 alternative hypothesis: true ratio of scales is not equal to 1 ##VARIANCE IS EQUAL### ##Simple 2-sample t-test### t.test(O_Prev~Treatment,data=oocysts2,var.equal=TRUE) Two Sample t-test data: O_Prev by Treatment t = -15.3206, df = 209, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.2410336 -0.1860753 sample estimates: mean in group 0 mean in group 1 0.4376480 0.6512024 ##Using GLM## prevlm<-lm(formula=O_Prev~Treatment*Block*Day,method="qr",model=TRUE,qr=TRUE,data=oocysts2) summary(prevlm) prevlm2<-update(prevlm,.~.-Treatment:Block:Day) summary(prevlm2) prevlm3<-update(prevlm2,.~.-Block:Day) summary(prevlm3) prevlm4<-update(prevlm3,.~.-Treatment:Day) summary(prevlm4) prevlm5<-update(prevlm4,.~.-Treatment:Block) summary(prevlm5) prevlm6<-update(prevlm5,.~.-Day) summary(prevlm6) head(fortify(prevlm5)) prevlm5fort<-fortify(prevlm5) ##REsiduals plot## r1<-ggplot(aes(x=.fitted,y=.resid),data=prevlm5fort)+geom_point()+geom_hline(yintercept=0)+geom_smooth(se=FALSE)+labs(x="Fitted Values",y="Residuals") r1 ##QQ Plot## qq<-ggplot(prevlm5fort,aes(sample=.stdresid))+stat_qq()+geom_abline() qq ##AIC AIC(prevlm,prevlm2,prevlm3,prevlm4,prevlm5,prevlm6) df AIC prevlm 9 -401.3285 prevlm2 8 -403.2911 prevlm3 7 -405.2810 prevlm4 6 -407.2450 prevlm5 5 -408.2306 prevlm6 4 -407.1649 ##QUESTIONS## -Where is coefplot()? ##Intensity## #First take original dataset and select out all the 0's## intensity<- oocysts[ which(oocysts$O_Intensity > 0), ] ##then transform dataset to cup level instead of individual level### intensity$EU<-as.factor(paste(intensity$Day,intensity$Cup,intensity$Treatment,intensity$Block,sep="-")) intensity2<-aggregate(cbind(Block,Treatment,Day,Cup,O_Intensity)~EU,data=intensity,mean) intensity2<-intensity2[with(intensity2,order(Block,Treatment,Cup,Day)),] ##then look at means## aggregate( formula = O_Intensity~Treatment+Block, data = intensity2,FUN = mean ) ##ansari test - variance is equal between treatments, but not between blocks## ##linear model## intensem<-lm(formula=O_Intensity~Treatment*Block*Day,method="qr",data=intensity2) summary(intensem) Call: lm(formula = O_Intensity ~ Treatment * Block * Day, data = intensity2, method = "qr") Residuals: Min 1Q Median 3Q Max -1.40637 -0.38739 -0.09965 0.26437 2.92696 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.16501 0.89155 1.307 0.192787 Treatment 8.53531 1.25781 6.786 1.24e-10 *** Block 0.08145 0.55032 0.148 0.882484 Day -0.02783 0.11601 -0.240 0.810686 Treatment:Block -3.77963 0.77334 -4.887 2.07e-06 *** Treatment:Day -0.77645 0.16360 -4.746 3.91e-06 *** Block:Day 0.03846 0.07172 0.536 0.592321 Treatment:Block:Day 0.37392 0.10066 3.715 0.000263 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.623 on 203 degrees of freedom Multiple R-squared: 0.5863, Adjusted R-squared: 0.572 F-statistic: 41.1 on 7 and 203 DF, p-value: < 2.2e-16 ##max model is best model## ##fortify data and do residual, qq plots## intensemfort<-fortify(intensem) ri<-ggplot(aes(x=.fitted,y=.resid),data=intensemfort)+geom_point()+geom_hline(yintercept=0)+geom_smooth(se=FALSE)+labs(x="Fitted Values",y="Residuals") ri qq<-ggplot(intensemfort,aes(sample=.stdresid))+stat_qq()+geom_abline() qq ##create new dataset for graphing## library(plyr) intplyr <- ddply(intensity2, c("Block", "Treatment","Day"), summarise,N = length(O_Intensity),mean = mean(O_Intensity),sd = sd(O_Intensity),se = sd / sqrt(N)) ##create predictions from plyr'd dataset## predict(intensem,intplyr,se.fit=TRUE) #make barplot## p <- ggplot(intplyr, aes(fill=factor(Treatment), y=mean, x=Day)) p+geom_bar(position="dodge",stat="identity")+facet_wrap(~Block) dodge <- position_dodge(width=0.9) p + geom_bar(position=dodge,stat="identity") + facet_wrap(~Block)+geom_errorbar(limits, position=dodge, width=0.25) ##Prev vs. Intensity## #by cup# #Correlation# cor(oocysts2$O_Prev,oocysts2$O_Intensity) [1] 0.740873 Treatment<-as.factor(oocysts2$Treatment) class(Treatment) [1] "factor" levels(Treatment) [1] "0" "1" pi<-ggplot(oocysts2,aes(x=O_Intensity,y=O_Prev)) pi+geom_point(aes(color=Treatment))+facet_wrap(~Block) pi<-ggplot(oocysts2,aes(x=O_Intensity,y=O_Prev)) pi+geom_point(aes(color=Treatment))+facet_wrap(~Block)+stat_smooth(method="lm",formula=y~ns(x,2),color="black")+labs(y="Oocyst Intensity",x="Oocyst Prevalence") pi<-ggplot(oocysts2,aes(x=O_Intensity,y=O_Prev)) ##creates first layer## pi+geom_point(aes(color=Treatment))+facet_wrap(~Block) ##this layer separates by block, colors treatments## +stat_smooth(method="lm",formula=y~ns(x,2),color="black") ##this fits quadratic line &CI bands to points## +labs(y="Oocyst Intensity",x="Oocyst Prevalence") #relabels x and y axes# +scale_color_discrete(name="Food Treatment",labels=c("Low Food","High Food")) #edits legend to say detailed treatment names# ###Looking at EIP - how best to model it with making predictions?## #Facet grid model with spline fits/spline CI's## eip<-ggplot(EIP2,aes(x=Day,y=Spz.Prev)) eip+geom_point()+facet_grid(Treatment~Block)+stat_smooth(method="lm",formula=y~ns(x,3),color="black")+labs(y="Day Post Infection",x="Sporozoite Prevalence") ##Making plots in the style of Killeen et al 2014 paper## #First make normal plots## #make column titles lower case# names(gklow1) <- tolower(names(gklow1)) #rename columns# library(plyr) gklow1<-rename(gklow1, c("cumulative.survival.estimate"="survived", "cumulative.infectious.estimate"="infectious","survived.and.infectious"="both")) #melt into long-form# gkl1<-melt(gklow1,id.vars=c("treatment","block","day")) gk<-ggplot(gkl1, aes(x = day, y = value, colour = variable)) + + geom_line() + + ylab(label="Proportion Alive and/or Infectious") + + xlab("AGE (days)") + + scale_colour_manual(values=c("blue", "red","green")) gk ###EIP### ###model plot### op<-par(mfrow=c(1,2),mar=c(5,5,3,2)) asdf<-EIP2[which(EIP2$Block==1),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 1") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve(coef(m.xyzQRSh)[1]/(1+exp(-coef(m.xyzQRSh)[2]*(x-coef(m.xyzQRSh)[3]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xyzQRSh)[1]+coef(m.xyzQRSh)[4])/(1+exp(-(coef(m.xyzQRSh)[2]+coef(m.xyzQRSh)[5])*(x-(coef(m.xyzQRSh)[3]+coef(m.xyzQRSh)[6])))),from=9,to=16,add=TRUE,lwd=4) asdf<-EIP2[which(EIP2$Block==2),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 2") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve((coef(m.xyzQRSh)[1]+coef(m.xyzQRSh)[7])/(1+exp(-(coef(m.xyzQRSh)[2]+coef(m.xyzQRSh)[8])* (x-coef(m.xyzQRSh)[3]+coef(m.xyzQRSh)[9]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xyzQRSh)[1]+coef(m.xyzQRSh)[4]+coef(m.xyzQRSh)[7]+coef(m.xyzQRSh)[10])/ (1+exp(-(coef(m.xyzQRSh)[2]+coef(m.xyzQRSh)[5]+coef(m.xyzQRSh)[8])* (x-(coef(m.xyzQRSh)[3]+coef(m.xyzQRSh)[6]+coef(m.xyzQRSh)[9])))),from=9,to=16,add=TRUE,lwd=4) ###upper limit curves### op<-par(mfrow=c(1,2),mar=c(5,5,3,2)) asdf<-EIP2[which(EIP2$Block==1),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 1") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve(coef(m.xyzQRSh)[1]/(1+exp(-coef(m.xyzQRSh)[2]*(x-coef(m.xyzQRSh)[3]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xyzQRSh)[1]+coef(m.xyzQRSh)[4])/(1+exp(-(coef(m.xyzQRSh)[2]+coef(m.xyzQRSh)[5])*(x-(coef(m.xyzQRSh)[3]+coef(m.xyzQRSh)[6])))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xyzQRSh)[1]+.01941)/(1+exp(-(coef(m.xyzQRSh)[2].24982)*(x-(coef(m.xyzQRSh)[3]+.10508))),from=9,to=16,add=TRUE,lwd=4,col="darkslategray") curve((coef(m.xyzQRSh)[1]+.01941)+(coef(m.xyzQRSh)[4]+.02280)/(1+exp(-((coef(m.xyzQRSh)[2]+.24982)+(coef(m.xyzQRSh)[5]+1.11348))*(x-((coef(m.xyzQRSh)[3]+.10508)+(coef(m.xyzQRSh)[6]+.11026)))),from=9,to=16,add=TRUE,col="darkslategray",lwd=4) asdf<-EIP2[which(EIP2$Block==2),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 2") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve((coef(m.xyzQRSh)[1]+coef(m.xyzQRSh)[7])/(1+exp(-(coef(m.xyzQRSh)[2]+coef(m.xyzQRSh)[8])* (x-coef(m.xyzQRSh)[3]+coef(m.xyzQRSh)[9]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xyzQRSh)[1]+coef(m.xyzQRSh)[4]+coef(m.xyzQRSh)[7]+coef(m.xyzQRSh)[10])/ (1+exp(-(coef(m.xyzQRSh)[2]+coef(m.xyzQRSh)[5]+coef(m.xyzQRSh)[8])* (x-(coef(m.xyzQRSh)[3]+coef(m.xyzQRSh)[6]+coef(m.xyzQRSh)[9])))),from=9,to=16,add=TRUE,lwd=4) Formula: Spz.Prev ~ (b + x * Treatment + Q * Blocki + h * Treatment * Blocki)/(1 + exp(-(k + y * Treatment + R * Blocki) * (Day - (tm + z * Treatment + S * Blocki)))) Parameters: Estimate Std. Error t value Pr(>|t|) b 0.41890 0.01941 21.583 < 2e-16 *** k 1.29672 0.24982 5.191 4.17e-07 *** tm 12.00639 0.10508 114.260 < 2e-16 *** x 0.23298 0.02280 10.220 < 2e-16 *** y 2.71724 1.11348 2.440 0.015328 * z -1.15451 0.11026 -10.471 < 2e-16 *** Q 0.08209 0.02288 3.588 0.000396 *** R 1.02775 0.50697 2.027 0.043635 * S -0.14323 0.07913 -1.810 0.071414 . h -0.19238 0.02737 -7.028 1.75e-11 *** > confint(m.xyzQRSh,parm=c("b","k","tm","x","y","Q","R","S","h")) Waiting for profiling to be done... 2.5% 97.5% b 0.38276229 0.458076845 k 0.93902003 1.880870021 tm 11.80320134 12.222750197 x 0.18778079 0.276434136 y 1.14900576 NA Q 0.03640728 0.126589272 R 0.11736252 2.747269687 S -0.31407669 0.001282859 h -0.24580989 -0.138627799 ##m.xzQh### Formula: Spz.Prev ~ (b + x * Treatment + Q * Blocki + h * Treatment * Blocki)/(1 + exp(-(k) * (Day - (tm + z * Treatment)))) Parameters: Estimate Std. Error t value Pr(>|t|) b 0.38471 0.01547 24.862 < 2e-16 *** k 3.13975 0.38293 8.199 9.99e-15 *** tm 11.90023 0.06160 193.186 < 2e-16 *** x 0.26449 0.02035 12.998 < 2e-16 *** z -1.14715 0.07868 -14.580 < 2e-16 *** Q 0.11194 0.02026 5.526 7.73e-08 *** h -0.20943 0.02659 -7.877 8.32e-14 *** > confint(m.xzQh) Waiting for profiling to be done... 2.5% 97.5% b 0.35404027 0.4157437 k 2.37011366 4.9365534 tm 11.74544092 12.0268481 x 0.22413312 0.3048710 z -1.30355151 -0.9650701 Q 0.07206854 0.1518201 h -0.26182800 -0.1570810 op<-par(mfrow=c(1,2),mar=c(5,5,3,2)) asdf<-EIP2[which(EIP2$Block==1),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 1") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) + lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve(coef(m.xzQh)[1]/(1+exp(-coef(m.xzQh)[2]*(x-coef(m.xzQh)[3]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xzQh)[1]+coef(m.xzQh)[4])/(1+exp(-(coef(m.xzQh)[2]*(x-(coef(m.xzQh)[3]+coef(m.xzQh)[5]))))),from=9,to=16,add=TRUE,lwd=4) asdf<-EIP2[which(EIP2$Block==2),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 2") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) + lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve((coef(m.xzQh)[1]+coef(m.xzQh)[6])/(1+exp(-coef(m.xzQh)[2]*(x-coef(m.xzQh)[3]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xzQh)[1]+coef(m.xzQh)[4]+coef(m.xzQh)[6]+coef(m.xzQh)[7])/(1+exp(-(coef(m.xzQh)[2])*(x-(coef(m.xzQh)[3]+coef(m.xzQh)[5])))),from=9,to=16,add=TRUE,lwd=4) ##Plotting 95% CI's### op<-par(mfrow=c(1,1),mar=c(5,5,3,2)) asdf<-EIP2[which(EIP2$Block==1),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 1") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve(coef(m.xzQh)[1]/(1+exp(-coef(m.xzQh)[2]*(x-coef(m.xzQh)[3]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xzQh)[1]+coef(m.xzQh)[4])/(1+exp(-(coef(m.xzQh)[2]*(x-(coef(m.xzQh)[3]+coef(m.xzQh)[5]))))),from=9,to=16,add=TRUE,lwd=4) curve(.41574/(1+exp(-(3.13975)*(x-11.7544))),from=9,to=16,add=TRUE,lwd=4,col="darkgreen") curve((.41574+.30487)/(1+exp(-(3.13975)*(x-(11.7544-1.30355)))),from=9,to=16,add=TRUE,lwd=4,col="darkred") curve(.35404/(1+exp(-(3.13975)*(x-12.026848))),from=9,to=16,add=TRUE,lwd=4,col="darkgreen") curve((.35404+.22413)/(1+exp(-(3.13975)*(x-(12.026848-.9650701)))),from=9,to=16,add=TRUE,lwd=4,col="darkred") op<-par(mfrow=c(1,1),mar=c(5,5,3,2)) asdf<-EIP2[which(EIP2$Block==2),] plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 2") for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) lines(Spz.Prev~Day,data=a,col=Cup);points(Spz.Prev~Day,data=a,col=Cup)} curve((coef(m.xzQh)[1]+coef(m.xzQh)[6])/(1+exp(-coef(m.xzQh)[2]*(x-coef(m.xzQh)[3]))),from=9,to=16,add=TRUE,lwd=4) curve((coef(m.xzQh)[1]+coef(m.xzQh)[4]+coef(m.xzQh)[6]+coef(m.xzQh)[7])/(1+exp(-(coef(m.xzQh)[2])*(x-(coef(m.xzQh)[3]+coef(m.xzQh)[5])))),from=9,to=16,add=TRUE,lwd=4) curve((.41574+.15182)/(1+exp(-(3.13975)*(x-(11.7544)))),from=9,to=16,add=TRUE,lwd=4,col="darkgreen") curve((.41574+.3048710+.15182-.15708)/(1+exp(-(3.13975)*(x-(11.7544-.96507)))),from=9,to=16,add=TRUE,lwd=4,col="darkred") curve((.35404+.07206)/(1+exp(-(3.13975)*(x-(12.0268)))),from=9,to=16,add=TRUE,lwd=4,col="darkgreen") curve((.35404+.07206+.224133-.261828)/(1+exp(-(3.13975)*(x-(12.0268-1.30355)))),from=9,to=16,add=TRUE,lwd=4,col="darkred") EIP plot with points instead of lines op<-par(mfrow=c(1,1),mar=c(5,5,3,2)) > asdf<-EIP2[which(EIP2$Block==1),] > plot(0,ylim=c(0,1),xlim=c(9,16),xlab=tlab,ylab=plab,main = "Block 1") > for(i in levels(factor(asdf$mbgroup))){a<-subset(asdf,mbgroup==i) + points(Spz.Prev~Day,data=a,col=ifelse(Treatment==0,"red","blue"),pch=16)} > curve(coef(m.xzQh)[1]/(1+exp(-coef(m.xzQh)[2]*(x-coef(m.xzQh)[3]))),from=9,to=16,add=TRUE,lwd=4) > curve((coef(m.xzQh)[1]+coef(m.xzQh)[4])/(1+exp(-(coef(m.xzQh)[2]*(x-(coef(m.xzQh)[3]+coef(m.xzQh)[5]))))),from=9,to=16,add=TRUE,lwd=4)