library(ggplot2) library(lme4) library(Matrix) library(MASS) library(lmerTest) library(emmeans) library(ggpubr) library(msme) library(olsrr) library(glmmTMB) library(DHARMa) library(MuMIn) library(effects) library(blmeco) #dataset including all data Expdata<-read.table("Prevalence for R.txt", header=T, sep="\t") #Converting to correct variables Expdata$colony<-factor(Expdata$colony) Expdata$treatment<-factor(Expdata$treatment) Expdata$infected<-factor(Expdata$infected, levels=c(0, 1), labels=c("uninfected", "infected")) #####optimal model m4<-glmer(Expdata$infected~ Expdata$treatment + Expdata$day_screened + Expdata$mass + Expdata$day_screened*Expdata$treatment + (1|Expdata$colony), family=binomial, data = Expdata) summary(m4) #comparing likelihood of infection across treatment groups, using tukey to correct for multiple comparisons emmeans(m4, list(pairwise ~ treatment), adjust = "tukey") #significant difference between all treatment comparisons #need to inverse log coefficients to find change in likelihood of infection in different treatments #inverse logit of means of each treatment exp(-1.12) #=0.326 for treatment group 1 #LCL exp(-1.7) #UCL exp(-0.54) exp(-3.33) #treatment group 2=0.03579 #LCL exp(-3.98) #UCL exp(-2.69) exp(-6.69) #treatment group 3=0.001243 #LCL exp(-8.93) #UCL exp(-4.45) #inverse logit of treatment coefficients #1-2 there's a 9.116% higher likelihood of being infected exp(2.21) #1-3 there's a 262.4341% higher likelihood exp(5.57) #2-3 there's a 28.789% higher likelihood exp(3.36) #m4 is best model so far, check assumptions residualsm4 <- simulateResiduals(m4) plot(residualsm4) #fits assumptions!!!!!!! :) :) :) #checked for overdispersion :):) dispersion_glmer(m4) #use drop1 to find if reduced model is a better fit drop1(m4, test="Chisq") #dropping mass and interactions increases the AIC and have sig different p values to m4 #therefore m4 is still best model. #to get chi square value for interaction term create model without interaction m4nointeraction<-glmer(Expdata$infected~ Expdata$treatment + Expdata$day_screened + Expdata$mass + (1|Expdata$colony), family=binomial, data = Expdata) summary(m4nointeraction) #chi square test between model with and without interaction anova(m4, m4nointeraction, test = "Chisq") #checking is random effect colony intercepts are different, looks like they are randoms <- ranef(m4) dotplot(randoms) #####Analyzing difference in R rate between treatment groups, glmer wouldn't work, couldn't get correct distribution #Using R number calculated (Colony count x proportion of infected in sample) Rdata3<-read.table("R rate by treatment group including dopey and doc.txt", header=T, sep="\t") #changing to correct variable type Rdata3$treatment<-factor(Rdata3$treatment, levels = c(1, 2, 3)) #log transforming R number logRnumber<-log(Rdata3$R_number) #check distribution of log r number hist(logRnumber) #more normal! #adding transformed r values to dataset Rdata4 Rdata4<-cbind(Rdata3, logRnumber) anova1<-aov(logRnumber ~ treatment, data=Rdata4) summary(anova1) plot(anova1) TukeyHSD(anova1) #significant difference between R number of 3-1 not others #try anova anova1<-aov(R_number ~ treatment, data=Rdata4) summary(anova1) #treatment has sig effect plot(anova1) TukeyHSD(anova1) #no significant difference #doesn't fit normality assumptions, try kruskall wallis kruskal.test(R_number ~ treatment, data = Rdata4) #no significant effect of treatment on R_number #try kruskall wallis on percentage prevalence after 7 days kruskal.test(percentage_infected ~ treatment, data = Rdata4) #Also not significant ##Repeating same analysis for R number without Doc and Dopey Rdata4<-read.table("R rate by treatment group not including dopey and doc.txt", header=T, sep="\t") #changing to correct variable type Rdata4$treatment<-factor(Rdata4$treatment, levels = c(1, 2, 3)) #log transforming R number logRnumber<-log(Rdata4$R_number) #check distribution of log r number hist(logRnumber) #more normal! #adding transformed r values to dataset Rdata4 Rdata4<-cbind(Rdata4, logRnumber) anova1<-aov(logRnumber ~ treatment, data=Rdata4) summary(anova1) plot(anova1) TukeyHSD(anova1) #significant difference between R number of 3-1 not others #try anova anova1<-aov(R_number ~ treatment, data=Rdata4) summary(anova1) #treatment has sig effect plot(anova1) TukeyHSD(anova1) #no significant difference #doesn't fit normality assumptions, try kruskall wallis kruskal.test(R_number ~ treatment, data = Rdata4) #no significant effect of treatment on R_number #try kruskall wallis on percentage prevalence after 7 days kruskal.test(percentage_infected ~ treatment, data = Rdata4) #Also not significant #######effect of treatment on colony growth colonycounts<-read.table("colony counts for R.txt", header=T, sep="\t") colonycounts$colony<-factor(colonycounts$colony) colonycounts$treatment<-factor(colonycounts$treatment) #distribution of counts hist(colonycounts$count) #glmer using poisson distribution countm1<-glmer(colonycounts$count~colonycounts$treatment+colonycounts$week+(1|colonycounts$colony), data = colonycounts, family = "poisson") summary(countm1) emmeans(countm1, list(pairwise ~ treatment), adjust = "tukey") dispersion_glmer(countm1) #doesn't fit assumptions? residualscountm1 <- simulateResiduals(countm1) plot(residualscountm1) #test dispersion of residuals testDispersion(residualscountm1) #p value not significant therefore not overdispersed i think, dpersion is 1.533 should be worries if over 2. #negative binomial glmer.nb, best model for colony counts countm3<-glmer.nb(colonycounts$count~colonycounts$treatment+colonycounts$week+colonycounts$treatment:colonycounts$week +(1|colonycounts$colony), data = colonycounts) summary(countm3) ##comparing model with and without treatment countm4<-glmer.nb(colonycounts$count~colonycounts$week+colonycounts$treatment:colonycounts$week + (1|colonycounts$colony), data = colonycounts) anova(countm4, countm3, test = "Chisq") ##shows treatment was not significant ##testing sig of interaction countm5<-glmer.nb(colonycounts$count~colonycounts$treatment+colonycounts$week +(1|colonycounts$colony), data = colonycounts) anova(countm5, countm3, test = "Chisq") ##interaction was significant emmeans(countm3, list(pairwise ~ treatment), adjust = "tukey") dispersion_glmer(countm3) #assumptions checking residualscountm3 <- simulateResiduals(countm3) plot(residualscountm3) #more patterns in residuals and more overdispersed than m1 ####Rerunning model m4 on subset of data, up until day 28 #dataset including all data prevalenceday28<-read.table("Prevalence for R up to day 28.txt", header=T, sep="\t") #Converting to correct variables prevalenceday28$colony<-factor(prevalenceday28$colony) prevalenceday28$treatment<-factor(prevalenceday28$treatment) prevalenceday28$infected<-factor(prevalenceday28$infected, levels=c(0, 1), labels=c("uninfected", "infected")) #rerun model 4 on days 7-28 m4day28<-glmer(prevalenceday28$infected~ prevalenceday28$treatment + prevalenceday28$day_screened + prevalenceday28$mass + prevalenceday28$day_screened:prevalenceday28$treatment + (1|prevalenceday28$colony), family=binomial, data = prevalenceday28) summary(m4day28) drop1(m4day28, test = "Chisq") #comparing models dropping terms m4day28notreatment <- glmer(prevalenceday28$infected~ prevalenceday28$day_screened + prevalenceday28$mass + prevalenceday28$day_screened:prevalenceday28$treatment + (1|prevalenceday28$colony), family=binomial, data = prevalenceday28) summary(m4day28notreatment) #difference in liklihood of infection still significant emmeans(m4day28, list(pairwise ~ treatment), adjust = "tukey") #getting chi sq statistic for model parts, running model without each term m4day28nointeraction<-glmer(prevalenceday28$infected~ prevalenceday28$treatment + prevalenceday28$day_screened + prevalenceday28$mass + (1|prevalenceday28$colony), family=binomial, data = prevalenceday28) #chisq test between models anova(m4day28notreatment, m4day28, test = "Chisq") anova(m4day28nointeraction, m4day28) ####treatment not significant anova(m4day28noday, m4day28, test = "Chisq") ###fails anova(m4day28nomass, m4day28, test = "Chisq") ####mass is significant anova(m4day28nointeraction, m4day28, test = "Chisq") drop1(m4day28, test = "Chisq") #gets chi sq value for interaction and mass drop1(m4day28nointeraction, test = "Chisq") #underdispersed, probs alright dispersion_glmer(m4day28) #checking residuals, best residuals I've ever seen :) residualsm4day28 <- simulateResiduals(m4day28) plot(residualsm4day28) ###logistic binary regression male prevalence maleprev<-read.table("males prevalence for r.txt", header=T, sep="\t") #Converting to correct variables maleprev$colony<-factor(maleprev$colony) maleprev$treatment<-factor(maleprev$treatment) maleprev$infection<-factor(maleprev$infection, levels=c(0, 1), labels=c("uninfected", "infected")) malem1 <- glm(maleprev$infection ~ maleprev$treatment + maleprev$mass + maleprev$colony, data = maleprev, family = binomial) summary(malem1) #treatment is not significant drop1(malem1, test="Chisq") anova(malem1) #model without treatment malem2<-glm(maleprev$infection ~ maleprev$mass + maleprev$colony, data = maleprev, family = binomial) anova(malem1, test = "Chisq") drop1(malem1, test = "F")