#check what is in r console from last time ls() #remove everything rm(list=ls()) #check what working directory is set getwd() #setting to required one where data is based setwd("C:/Users/UEA/Documents/Dissertation and phd/d- data for phd/Neat data/life stage recovery") #reading in data survivalstage <- read.csv("heatwavesurvivallifestage2.csv", header = TRUE) # 10D data set ### checking for outliers/errors summary(survivalstage) # produces general (unsplit) range, quantiles, median, count and mean summary stats for each variable str(survivalstage) # checks the variable types # 'data.frame': 190 obs. of 10 variables: # $ Replicate : int 1 2 3 4 5 6 7 8 9 10 ... # $ Sex : Factor w/ 3 levels "female","male",..: 2 2 2 1 1 1 1 1 1 1 ... # $ Temperature.oC : int 30 30 30 30 30 30 30 30 30 30 ... # $ Lifestage : Factor w/ 5 levels "control","immature",..: 1 1 1 1 1 1 1 1 1 1 ... # $ Composite.treatment: Factor w/ 9 levels "30control","40immature",..: 1 1 1 1 1 1 1 1 1 1 ... combination of temperature and lifestage # $ N : int 20 20 20 20 14 20 20 20 20 20 ... number in group treated # $ Survivors : int 19 19 20 20 12 18 19 18 17 18 ... survivors in group # $ Deaths : int 1 1 0 0 2 2 1 2 3 2 ... deaths in group # $ Prop.survivors : num 0.95 0.95 1 1 0.86 0.9 0.95 0.9 0.85 0.9 ... # $ Batch : int 1 1 1 1 1 1 1 1 1 1 ... generation assayed str(survivalstage) is.na(survivalstage) # returns TRUE of x is missing # nothing missing survivalstage$Temperature.oC <- as.factor(survivalstage$Temperature.oC) levels(survivalstage$Composite.treatment) # "30control" "40immature" "40larvae" "40pupae" "41pupae" "42immature" "42larvae" "42mature" "42pupae" # reordering levels survivalstage$Composite.treatment <- factor(survivalstage$Composite.treatment,levels = c("30control","40larvae","42larvae","40pupae","41pupae","42pupae","40immature","42immature","42mature")) levels(survivalstage$Composite.treatment) levels(survivalstage$Temperature.oC) levels(survivalstage$Sex) # [1] "female" "male" survivalstage$Cbindsurvival<-cbind(survivalstage$Survivors, survivalstage$Deaths) # cbind response for model #! # complete separation will intefere with statistical testing, excluding from analysis survivalstageno4142pupae<-subset(survivalstage, !(Composite.treatment %in% c("41pupae","42pupae"))) survivalstageno4142pupae$Composite.treatment<-droplevels(survivalstageno4142pupae$Composite.treatment) names(survivalstageno4142pupae) survivalstageno41pupae<-subset(survivalstage, !(Composite.treatment %in% c("41pupae"))) survivalstageno41pupae$Composite.treatment<-droplevels(survivalstageno4142pupae$Composite.treatment) names(survivalstageno41pupae) ######################### male only life-stage ########## survivalstagemale <- subset(survivalstage, Sex %in% c("male")) names(survivalstagemale) ##################### ! PAPER FIG 3 PLOT ###################### names(survivalstage) str(survivalstage) library(ggplot2) temp <- expression(paste('Temperature (',degree,'C)')) #the temperature label with degrees sign # lightcoral # getting hex code indianred1 #205 106 106 higher = lighter rgb(255,106,106, maxColorValue = 255)#FF6A6A" # same as indian red rgb(255,140,140, maxColorValue = 255)#FF9696" #stage survival boxplot graphlifesur<-ggplot(survivalstageno41pupae, aes(x=Composite.treatment, y=Prop.survivors, fill= Composite.treatment)) + #change fill to colour is just lines and change 'scale_fill_manual' below to scale_color_manual geom_boxplot(notch=F, #change to F if want to get rid of notchs outlier.shape= NA, #shape of the outlier (hashtag out if dont want outliers marked) width=0.5, lwd=0.5, fatten=0.5, color="black", position=position_dodge(0.5)) + #size of the outlier (hashtag out if dont want outliers marked) stat_summary(fun.y="mean", geom= "point", size=3, position=position_dodge(1), color="black") + scale_fill_manual(values=c("white", "#FF9696", "firebrick2", "#FF9696", "firebrick2", "#FF9696", "firebrick2", "firebrick2"), # changes the colour of the bars name = temp, #adds in temperature label on the legend breaks = c("30control","40larvae","42larvae","40pupae","42pupae","40immature","42immature","42mature"), #the order listed in the legend label = c('Control\n \n30°C','Late \nlarvae\n40°C','Late \nlarvae\n42°C','Pupae\n \n40°C','Pupae\n \n42°C','Immature \nadult\n40°C','Immature \nadult\n42°C','Mature \nadult\n42°C')) + #how things are labeled in the lgend geom_jitter(shape=1, size=1.5, position=position_jitter(0.15)) + #so all the data points are not ontop of each other labs (x= "Heatwave treatment life stage group" , y="Proportion surviving heatwave treatment") + #adding title to the x axis and y axis scale_x_discrete(breaks = c("30control","40larvae","42larvae","40pupae","42pupae","40immature","42immature","42mature"), #the order of the variables on the x axis label = c("Control\n \n30°C",'Late \nlarvae\n40°C','Late \nlarvae\n42°C','Pupae\n \n40°C','Pupae\n \n42°C','Immature \nadult\n40°C','Immature \nadult\n42°C','Mature \nadult\n42°C')) + # the names on the x axis coord_cartesian(ylim=c(-0.02, 1.02)) + #set axis limits scale_y_continuous(breaks=seq(0, 1, 0.2), #ticks from 0 to 16000 and show number every 16000 expand = c(0, 0)) + #cuts the axis off at theme_classic() + #the theme of the whole plot theme( #legend.position="none", #get rid of the hashtag to get rid of legend panel.grid.major=element_blank(), #getting rid of majorgridlines panel.border=element_blank(), #getting rid of minorgridlines panel.grid.minor=element_blank(), axis.line.x=element_line(color="black", size = 1), axis.line.y=element_line(color="black", size = 1), axis.text.x=element_text(color="black", size=8), axis.text.y=element_text(color="black", size=12), axis.title.x=element_text(face = "bold", size=12, color="black", margin = margin(t = 10, r = 0, b = 0, l = 0)), axis.title.y=element_text(face = "bold", size=12, color="black", margin = margin(t = 0, r = 10, b = 0, l = 0)), legend.position="none", panel.background=element_blank(), plot.background=element_rect(fill="transparent", colour = NA)) setwd("C:/Users/UEA/Desktop") ggsave("fig3lifesur.png",width=6, height=4, dpi=300, bg = "transparent") setwd("C:/Users/UEA/Documents/Dissertation and phd/d- data for phd/Neat data/life stage recovery") #"white", "dark red", "firebrick2", "dark red", "firebrick2","firebrick2", "dark red", "firebrick2", "firebrick2" # "white", "red3", "firebrick2", "red3", "firebrick2","firebrick2", "red3", "firebrick2", "firebrick2" #stage survival boxplot with 41oC pupae graphlifesur<-ggplot(survivalstage, aes(x=Composite.treatment, y=Prop.survivors, fill= Composite.treatment)) + #change fill to colour is just lines and change 'scale_fill_manual' below to scale_color_manual geom_boxplot(notch=F, #change to F if want to get rid of notchs outlier.shape= NA, #shape of the outlier (hashtag out if dont want outliers marked) width=0.5, lwd=0.5, fatten=0.5, color="black", position=position_dodge(0.5)) + #size of the outlier (hashtag out if dont want outliers marked) stat_summary(fun.y="mean", geom= "point", size=3, position=position_dodge(1), color="black") + scale_fill_manual(values=c("white", "dark red", "firebrick2", "dark red", "firebrick2","firebrick2", "dark red", "firebrick2", "firebrick2"), # changes the colour of the bars name = temp, #adds in temperature label on the legend breaks = c("30control","40larvae","42larvae","40pupae","41pupae","42pupae","40immature","42immature","42mature"), #the order listed in the legend label = c('Control\n \n30°C','Late \nlarvae\n40°C','Late \nlarvae\n42°C','Pupae\n \n40°C','Pupae\n \n41°C','Pupae\n \n42°C','Immature \nadult\n40°C','Immature \nadult\n42°C','Mature \nadult\n42°C')) + #how things are labeled in the lgend geom_jitter(shape=1, size=1.5, position=position_jitter(0.15)) + #so all the data points are not ontop of each other labs (x= "Heatwave treatment life stage group" , y="Proportion surviving heatwave treatment") + #adding title to the x axis and y axis scale_x_discrete(breaks = c("30control","40larvae","42larvae","40pupae","41pupae","42pupae","40immature","42immature","42mature"), #the order of the variables on the x axis label = c("Control\n \n30°C",'Late \nlarvae\n40°C','Late \nlarvae\n42°C','Pupae\n \n40°C','Pupae\n \n41°C','Pupae\n \n42°C','Immature \nadult\n40°C','Immature \nadult\n42°C','Mature \nadult\n42°C')) + # the names on the x axis coord_cartesian(ylim=c(-0.02, 1.02)) + #set axis limits scale_y_continuous(breaks=seq(0, 1, 0.2), #ticks from 0 to 16000 and show number every 16000 expand = c(0, 0)) + #cuts the axis off at 0 theme_classic() + #the theme of the whole plot theme( #legend.position="none", #get rid of the hashtag to get rid of legend panel.grid.major=element_blank(), #getting rid of majorgridlines panel.border=element_blank(), #getting rid of minorgridlines panel.grid.minor=element_blank(), axis.line.x=element_line(color="black", size = 1), axis.line.y=element_line(color="black", size = 1), axis.text.x=element_text(color="black", size=8), axis.text.y=element_text(color="black", size=12), axis.title.x=element_text(face = "bold", size=12, color="black", margin = margin(t = 10, r = 0, b = 0, l = 0)), axis.title.y=element_text(face = "bold", size=12, color="black", margin = margin(t = 0, r = 10, b = 0, l = 0)), legend.position="none", panel.background=element_blank(), plot.background=element_rect(fill="transparent", colour = NA)) setwd("C:/Users/UEA/Desktop") ggsave("fig3lifesur.png",width=6, height=4, dpi=300, bg = "transparent") setwd("C:/Users/UEA/Documents/Dissertation and phd/d- data for phd/Neat data/life stage recovery") #"white", "dark red", "firebrick2", "dark red", "firebrick2","firebrick2", "dark red", "firebrick2", "firebrick2" # "white", "red3", "firebrick2", "red3", "firebrick2","firebrick2", "red3", "firebrick2", "firebrick2" names(survivalstagemale) #stage survival boxplot showing male subset of data little different to both sexes graphlifesurmale<-ggplot(survivalstagemale, aes(x=Composite.treatment, y=Prop.survivors, fill= Composite.treatment)) + #change fill to colour is just lines and change 'scale_fill_manual' below to scale_color_manual geom_boxplot(notch=F, #change to F if want to get rid of notchs outlier.shape= NA, #shape of the outlier (hashtag out if dont want outliers marked) width=0.5, lwd=0.5, fatten=0.5, color="black", position=position_dodge(0.5)) + #size of the outlier (hashtag out if dont want outliers marked) stat_summary(fun.y="mean", geom= "point", size=3, position=position_dodge(1), color="black") + scale_fill_manual(values=c("white", "dark red", "dark red", "firebrick2","firebrick2", "dark red", "firebrick2", "firebrick2"), # changes the colour of the bars name = temp, #adds in temperature label on the legend breaks = c("30control","40larvae","40pupae","41pupae","42pupae","40immature","42immature","42mature"), #the order listed in the legend label = c('Control\n \n30°C','Late \nlarvae\n40°C','Pupae\n \n40°C','Pupae\n \n41°C','Pupae\n \n42°C','Immature \nadult\n40°C','Immature \nadult\n42°C','Mature \nadult\n42°C')) + #how things are labeled in the lgend geom_jitter(shape=1, size=1.5, position=position_jitter(0.15)) + #so all the data points are not ontop of each other labs (x= "Heatwave treatment life stage group" , y="Proportion surviving heatwave treatment") + #adding title to the x axis and y axis scale_x_discrete(breaks = c("30control","40larvae","40pupae","41pupae","42pupae","40immature","42immature","42mature"), #the order of the variables on the x axis label = c("Control\n \n30°C",'Late \nlarvae\n40°C','Pupae\n \n40°C','Pupae\n \n41°C','Pupae\n \n42°C','Immature \nadult\n40°C','Immature \nadult\n42°C','Mature \nadult\n42°C')) + # the names on the x axis coord_cartesian(ylim=c(-0.02, 1.02)) + #set axis limits scale_y_continuous(breaks=seq(0, 1, 0.2), #ticks from 0 to 16000 and show number every 16000 expand = c(0, 0)) + #cuts the axis off at 0 theme_classic() + #the theme of the whole plot theme( #legend.position="none", #get rid of the hashtag to get rid of legend panel.grid.major=element_blank(), #getting rid of majorgridlines panel.border=element_blank(), #getting rid of minorgridlines panel.grid.minor=element_blank(), axis.line.x=element_line(color="black", size = 1), axis.line.y=element_line(color="black", size = 1), axis.text.x=element_text(color="black", size=8), axis.text.y=element_text(color="black", size=12), axis.title.x=element_text(face = "bold", size=12, color="black", margin = margin(t = 10, r = 0, b = 0, l = 0)), axis.title.y=element_text(face = "bold", size=12, color="black", margin = margin(t = 0, r = 10, b = 0, l = 0)), legend.position="none", panel.background=element_blank(), plot.background=element_rect(fill="transparent", colour = NA)) setwd("C:/Users/UEA/Desktop") ggsave("survmale.png",width=5.5, height=4, dpi=300, bg = "transparent") setwd("C:/Users/UEA/Documents/Dissertation and phd/d- data for phd/Neat data/life stage recovery") ############## ! DESCRIPTIVE STATS ########################### library(psych) #gives you vars n, mean, sd, median, trimmed, mad, min, max, range, skew, kurtosis, se describeBy(survivalstage$Prop.survivors, survivalstage$Composite.treatment,mat=TRUE) # item group1 vars n mean sd median trimmed mad min max range skew kurtosis se # X11 1 30control 1 68 0.91338235 0.06318913 0.900 0.9171429 0.074130 0.75 1.00 0.25 -0.5389232 -0.01418211 0.007662808 # X12 2 40larvae 1 9 0.97222222 0.03632416 1.000 0.9722222 0.000000 0.90 1.00 0.10 -0.7012161 -1.02688007 0.012108053 # X13 3 42larvae 1 38 0.06947368 0.09037302 0.050 0.0568750 0.074130 0.00 0.32 0.32 1.1699690 0.11949314 0.014660440 # X14 4 40pupae 1 17 0.61117647 0.21242300 0.600 0.6193333 0.148260 0.20 0.90 0.70 -0.3871881 -0.63386448 0.051520144 # X15 5 41pupae 1 3 0.05000000 0.00000000 0.050 0.0500000 0.000000 0.05 0.05 0.00 NaN NaN 0.000000000 # X16 6 42pupae 1 5 0.00000000 0.00000000 0.000 0.0000000 0.000000 0.00 0.00 0.00 NaN NaN 0.000000000 # X17 7 40immature 1 10 0.44400000 0.31542564 0.400 0.4300000 0.378063 0.05 0.95 0.90 0.3111712 -1.54716877 0.099746345 # X18 8 42immature 1 22 0.19000000 0.20239048 0.175 0.1638889 0.222390 0.00 0.88 0.88 1.6476541 3.44395835 0.043149794 # X19 9 42mature 1 18 0.60388889 0.26086144 0.700 0.6168750 0.252042 0.10 0.90 0.80 -0.5628777 -1.14143049 0.061485632 mean(survivalstageno41pupae$N) # 19.5 sd(survivalstageno41pupae$N) # 4.5 describeBy(survivalstage$N, survivalstage$Composite.treatment,mat=TRUE) ################################################################################################################################ NEW METHOD: USE GLM(M) WITH NON-GAUSSIAN ERROR STRUCTURE###################################################### library(car); library(MASS); library (lme4); library (nlme) library(glmmADMB)# glmmADMB() ############## MODEL SELECTION ########################### ########### Pupae present #### Binomial family error structures # As data is proportion bound at 0 and 1 fitting normal distibution does not give normal and homogenity of variance in residuals # Creating a global model globalmodbinom<-glm(Cbindsurvival ~ Composite.treatment, binomial(link = "logit"), data=survivalstage) globalmodbibnomLOG<-glm(Cbindsurvival ~ Composite.treatment, binomial(link = "log"), data=survivalstage) summary(globalmodbinom); summary(globalmodbibnomLOG); # No R^2, AIC given # AIC: 926.6, AIC: 926.6 # link change seem to do little pseudoR<-(globalmodbinom$null.deviance-globalmodbinom$deviance) / globalmodbinom$null.deviance # (thomas et al., 2015) pseudoR # 0.8109647 pseudoR<-(globalmodbibnomLOG$null.deviance-globalmodbibnomLOG$deviance) / globalmodbibnomLOG$null.deviance # (thomas et al., 2015) pseudoR # 0.8109647 # seems changing the link does nothing to R^2 or AIC AICc<-(-2*logLik(globalmodbinom))+((2*1*(1+1)/(190-1-1))); AICc # qAICc<-((-2*logLik(model1)/Theta)+((2*p*(p+1)/(n-p-1))); qAICc # AIC correcting for perameters(p) and sample size (n) # 314.3862 qAICc<-(-2*logLik(globalmodbinom)/2.884871)+((2*1*(1+1)/(75-1-1))); qAICc # 315.0083 ## Overdispersion check par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) theta<-globalmodbinom$deviance/globalmodbinom$df.residual; theta #dispersion perameter (thomas et al 2015) how much variation left unexplained after fitting distribution # theta = 2.884871, overdispersed is >1 is overdispersion. hist(survivalstage$Cbindsurvival, breaks= 50) # over 80% cases are 0, very positive skew. recommendation of using quasi binomial globalmodquasibin<-glm(Cbindsurvival ~ Composite.treatment, quasibinomial(link = "logit"), data=survivalstage) globalmodqausibibinmLOG<-glm(Cbindsurvival ~ Composite.treatment, quasibinomial(link = "log"), data=survivalstage) globalmodqausibibinmID<-glm(Cbindsurvival ~ Composite.treatment, quasibinomial(link = "identity"), data=survivalstage) # errpr summary(globalmodquasibin); summary(globalmodqausibibinmLOG); summary(globalmodqausibibinmID) pseudoR<-(globalmodquasibin$null.deviance-globalmodquasibin$deviance) / globalmodquasibin$null.deviance # (thomas et al., 2015) pseudoR<-(globalmodqausibibinmLOG$null.deviance-globalmodqausibibinmLOG$deviance) / globalmodqausibibinmLOG$null.deviance # (thomas et al., 2015) pseudoR<-(globalmodqausibibinmID$null.deviance-globalmodqausibibinmID$deviance) / globalmodqausibibinmID$null.deviance # (thomas et al., 2015) pseudoR # 0.6212537 # all produce the same output and pseudoR # 1) Errors normally distributed? - YES BUT NOT IMPROVED # binomial devresid<-resid(globalmodbinom, type = "deviance"); hist(devresid) shapiro.test(devresid);ks.test(devresid, pnorm) qqnorm(devresid,cex=1.8,pch=20); qqline(devresid,lty=2,lwd=2) par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) # quasibinomial devresid<-resid(globalmodquasibin, type = "deviance"); hist(devresid) shapiro.test(devresid);ks.test(devresid, pnorm) qqnorm(devresid,cex=1.8,pch=20); qqline(devresid,lty=2,lwd=2) par(mfrow=c(2,2)); plot(globalmodquasibin);par(mfrow=c(1,1)) # quasibinom Q-Q similar ok, devresid histogram pl and KS test similarly failed # 2) Homogenous/homoscedasticity variance of residuals - NO BUT NOT IMPROVED # binomial devresid<-resid(globalmodbinom, type = "deviance") plot(devresid ~ globalmodbinom$fitted.values, pch = 20, cex = 1, cex.lab = 1.5) fligner.test(devresid~survivalstage$Composite.treatment) par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) # quasibinomial devresid<-resid(globalmodquasibin, type = "deviance") plot(devresid ~ globalmodquasibin$fitted.values, pch = 20, cex = 1, cex.lab = 1.5) fligner.test(devresid~survivalstage$Composite.treatment) par(mfrow=c(2,2)); plot(globalmodquasibin);par(mfrow=c(1,1)) # similar wedging, test failed # 3) Independences of independent variables - NO PATTERNS WITH SAMPLE SIZE devresid<-resid(globalmodbinom, type = "deviance") plot(devresid ~ survivalstage$N) devresid<-resid(globalmodquasibin, type = "deviance") plot(devresid ~ survivalstage$N) # no real patterns # 4) No serial auto-correlation with time/space - No directional pattern over time #! need library(car) durbinWatsonTest(globalmodbinom) durbinWatsonTest(globalmodquasibin) # test failed devresid<-resid(globalmodbinom, type = "deviance") plot(devresid ~ survivalstage$Batch) plot(devresid ~ survivalstage$Replicate) devresid<-resid(globalmodquasibin, type = "deviance") plot(devresid ~ survivalstage$Batch) plot(devresid ~ survivalstage$Replicate) # 5) No bias by unduly influential datapoints - YES # binomial par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) influence<-influence.measures(globalmodbinom); summary(influence) # quasibinomial par(mfrow=c(2,2)); plot(globalmodquasibin);par(mfrow=c(1,1)) influence<-influence.measures(globalmodquasibin); summary(influence) # cooks distance of outliers lowered # 6) Independent variables measured without error - BEST OF ABILITY ########### Pupae absent #### Binomial family error structures # As data is proportion bound at 0 and 1 fitting normal distibution does not give normal and homogenity of variance in residuals # Creating a global model globalmodbinom<-glm(Cbindsurvival ~ Composite.treatment, binomial(link = "logit"), data=survivalstageno4142pupae) globalmodbibnomLOG<-glm(Cbindsurvival ~ Composite.treatment, binomial(link = "log"), data=survivalstageno4142pupae) summary(globalmodbinom); summary(globalmodbibnomLOG); # No R^2, AIC given # AIC: 916.76, AIC: 916.76 # link change seem to do little pseudoR<-(globalmodbinom$null.deviance-globalmodbinom$deviance) / globalmodbinom$null.deviance # (thomas et al., 2015) pseudoR # 0.7956907 pseudoR<-(globalmodbibnomLOG$null.deviance-globalmodbibnomLOG$deviance) / globalmodbibnomLOG$null.deviance # (thomas et al., 2015) pseudoR # 0.7956907 # seems changing the link does nothing to R^2 or AIC AICc<-(-2*logLik(globalmodbinom))+((2*1*(1+1)/(182-1-1))); AICc # qAICc<-((-2*logLik(model1)/Theta)+((2*p*(p+1)/(n-p-1))); qAICc # AIC correcting for perameters(p) and sample size (n) # 902.7785 qAICc<-(-2*logLik(globalmodbinom)/2.983736)+((2*1*(1+1)/(182-1-1))); qAICc # 302.5813 ## Overdispersion check par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) theta<-globalmodbinom$deviance/globalmodbinom$df.residual; theta #dispersion perameter (thomas et al 2015) how much variation left unexplained after fitting distribution # theta = 2.983736, overdispersed is >1 is overdispersion. hist(survivalstageno4142pupae$Cbindsurvival, breaks= 50) # over ~80% cases are 0, very positive skew. recommendation of using quasi binomial globalmodquasibin<-glm(Cbindsurvival ~ Temperature.oC, quasibinomial(link = "logit"), data=survivalstageno4142pupae) globalmodqausibibinmLOG<-glm(Cbindsurvival ~ Temperature.oC, quasibinomial(link = "log"), data=survivalstageno4142pupae) globalmodqausibibinmID<-glm(Cbindsurvival ~ Temperature.oC, quasibinomial(link = "identity"), data=survivalstageno4142pupae) summary(globalmodquasibin); summary(globalmodqausibibinmLOG); summary(globalmodqausibibinmID) pseudoR<-(globalmodquasibin$null.deviance-globalmodquasibin$deviance) / globalmodquasibin$null.deviance # (thomas et al., 2015) pseudoR<-(globalmodqausibibinmLOG$null.deviance-globalmodqausibibinmLOG$deviance) / globalmodqausibibinmLOG$null.deviance # (thomas et al., 2015) pseudoR<-(globalmodqausibibinmID$null.deviance-globalmodqausibibinmID$deviance) / globalmodqausibibinmID$null.deviance # (thomas et al., 2015) pseudoR # 0.6064984 # all produce the same output and pseudoR # 1) Errors normally distributed? - YES BUT NOT IMPROVED # binomial devresid<-resid(globalmodbinom, type = "deviance"); hist(devresid) shapiro.test(devresid);ks.test(devresid, pnorm) qqnorm(devresid,cex=1.8,pch=20); qqline(devresid,lty=2,lwd=2) par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) # quasibinomial devresid<-resid(globalmodquasibin, type = "deviance"); hist(devresid) shapiro.test(devresid);ks.test(devresid, pnorm) qqnorm(devresid,cex=1.8,pch=20); qqline(devresid,lty=2,lwd=2) par(mfrow=c(2,2)); plot(globalmodquasibin);par(mfrow=c(1,1)) # quasibinom Q-Q ok, devresid histogram ok and KS test similarly passed # 2) Homogenous/homoscedasticity variance of residuals - NOT NECASSARY BUT IMPROVED RELATIVE TO PUPAL INCLUSION # binomial devresid<-resid(globalmodbinom, type = "deviance") plot(devresid ~ globalmodbinom$fitted.values, pch = 20, cex = 1, cex.lab = 1.5) fligner.test(devresid~survivalstageno4142pupae$Composite.treatment) par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) # quasibinomial devresid<-resid(globalmodquasibin, type = "deviance") plot(devresid ~ globalmodquasibin$fitted.values, pch = 20, cex = 1, cex.lab = 1.5) fligner.test(devresid~survivalstageno4142pupae$Composite.treatment) par(mfrow=c(2,2)); plot(globalmodquasibin);par(mfrow=c(1,1)) # test failed but less wedging than with pupae in # 3) Independences of independent variables - NO PATTERNS WITH SAMPLE SIZE devresid<-resid(globalmodbinom, type = "deviance") plot(devresid ~ survivalstageno4142pupae$N) devresid<-resid(globalmodquasibin, type = "deviance") plot(devresid ~ survivalstageno4142pupae$N) # no pattern # 4) No serial auto-correlation with time/space - No directional pattern over time #! need library(car) durbinWatsonTest(globalmodquasibin) # test failed devresid<-resid(globalmodbinom, type = "deviance") plot(devresid ~ survivalstageno4142pupae$Batch) plot(devresid ~ survivalstageno4142pupae$Replicate) devresid<-resid(globalmodquasibin, type = "deviance") plot(devresid ~ survivalstageno4142pupae$Batch) plot(devresid ~ survivalstageno4142pupae$Replicate) # 5) No bias by unduly influential datapoints - YES # binomial par(mfrow=c(2,2)); plot(globalmodbinom);par(mfrow=c(1,1)) influence<-influence.measures(globalmodbinom); summary(influence) # quasibinomial par(mfrow=c(2,2)); plot(globalmodquasibin);par(mfrow=c(1,1)) influence<-influence.measures(globalmodquasibin); summary(influence) # much fewer outliers # 6) Independent variables measured without error - BEST OF ABILITY #PUPAE OR NO PUPAE? globalmodquasibinpup<-glm(Cbindsurvival ~ Composite.treatment, binomial(link = "logit"), data=survivalstage) globalmodquasibinnopup<-glm(Cbindsurvival ~ Composite.treatment, binomial(link = "logit"), data=survivalstageno4142pupae) summary(globalmodquasibinpup) # Estimate Std. Error z value Pr(>|z|) # (Intercept) 2.31448 0.09691 23.883 < 2e-16 *** # Composite.treatment40larvae 1.24087 0.46380 2.675 0.00746 ** # Composite.treatment42larvae -4.92001 0.17333 -28.385 < 2e-16 *** # Composite.treatment40pupae -1.86456 0.14764 -12.629 < 2e-16 *** # Composite.treatment41pupae -5.18616 0.60133 -8.625 < 2e-16 *** # Composite.treatment42pupae -20.91277 706.60709 -0.030 0.97639 # Composite.treatment40immature -2.44296 0.15668 -15.592 < 2e-16 *** # Composite.treatment42immature -3.80785 0.16176 -23.540 < 2e-16 *** # Composite.treatment42mature -1.97372 0.14919 -13.229 < 2e-16 *** summary(globalmodquasibinnopup) # Estimate Std. Error z value Pr(>|z|) # (Intercept) 2.31448 0.09691 23.883 < 2e-16 *** # Composite.treatment40larvae 1.24087 0.46379 2.675 0.00746 ** # Composite.treatment42larvae -4.92001 0.17333 -28.385 < 2e-16 *** # Composite.treatment40pupae -1.86456 0.14764 -12.629 < 2e-16 *** # Composite.treatment40immature -2.44296 0.15668 -15.592 < 2e-16 *** # Composite.treatment42immature -3.80785 0.16176 -23.540 < 2e-16 *** # Composite.treatment42mature -1.97372 0.14919 -13.229 < 2e-16 *** globalmodquasibinpup<-glm(Cbindsurvival ~ Composite.treatment, quasibinomial(link = "logit"), data=survivalstage) globalmodquasibinnopup<-glm(Cbindsurvival ~ Composite.treatment, quasibinomial(link = "logit"), data=survivalstageno4142pupae) summary(globalmodquasibinpup) # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 2.3145 0.1577 14.680 < 2e-16 *** # Composite.treatment40larvae 1.2409 0.7546 1.645 0.102 # Composite.treatment42larvae -4.9200 0.2820 -17.447 < 2e-16 *** # Composite.treatment40pupae -1.8646 0.2402 -7.763 5.92e-13 *** # Composite.treatment41pupae -5.1862 0.9783 -5.301 3.32e-07 *** # Composite.treatment42pupae -20.9128 1149.5853 -0.018 0.986 # Composite.treatment40immature -2.4430 0.2549 -9.584 < 2e-16 *** # Composite.treatment42immature -3.8078 0.2632 -14.469 < 2e-16 *** # Composite.treatment42mature -1.9737 0.2427 -8.132 6.51e-14 *** summary(globalmodquasibinnopup) # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 2.3145 0.1603 14.435 < 2e-16 *** # Composite.treatment40larvae 1.2409 0.7674 1.617 0.108 # Composite.treatment42larvae -4.9200 0.2868 -17.156 < 2e-16 *** # Composite.treatment40pupae -1.8646 0.2443 -7.633 1.42e-12 *** # Composite.treatment40immature -2.4430 0.2592 -9.424 < 2e-16 *** # Composite.treatment42immature -3.8078 0.2676 -14.227 < 2e-16 *** # Composite.treatment42mature -1.9737 0.2468 -7.996 1.68e-13 *** # ! pupae model doesnt make sense, says no difference between 42 pupae to controls; due to no variation # ! binomial models less reliable as says everything extremely difficent ############## ! MODEL SIGNIFICANCE ########################### ### MODEL REFINEMENT globalmodquasibin<-glm(Cbindsurvival ~ Composite.treatment, quasibinomial(link = "logit"), data=survivalstageno4142pupae); summary(globalmodquasibin) # RD 522.15 pseudoR<-(globalmodquasibin$null.deviance-globalmodquasibin$deviance) / globalmodquasibin$null.deviance; pseudoR # 0.7956907 nullmodquasibin<-glm(Cbindsurvival ~ 1, quasibinomial(link = "logit"), data=survivalstageno4142pupae); summary(nullmodquasibin) # RD 2555.7 pseudoR<-(nullmodquasibin$null.deviance-nullmodquasibin$deviance) / nullmodquasibin$null.deviance; pseudoR # 0 ### Minimum adequate model explains offspring significantly more than null model because # anova/lrtest(model, nullmodel) and aic functions do not work with quasibinomial anova(globalmodquasibin, test = "Chisq") drop1(globalmodquasibin, test = "Chisq") # Model: # Cbindsurvival ~ Composite.treatment # Df Deviance scaled dev. Pr(>Chi) # 522.15 # Composite.treatment 6 2555.70 742.84 < 2.2e-16 *** anova(globalmodquasibin, nullmodquasibin, test = "Chisq") # Model 1: Cbindsurvival ~ Composite.treatment # Model 2: Cbindsurvival ~ 1 # Resid. Df Resid. Dev Df Deviance Pr(>Chi) # 1 175 522.15 # 2 181 2555.70 -6 -2033.5 < 2.2e-16 *** # --- ############## ! MODEL POST HOC ########################### summary(globalmodquasibin) #glm(formula = Cbindsurvival ~ Composite.treatment, family = quasibinomial(link = "logit"), # data = survivalstageno4142pupae) # # Deviance Residuals: # Min 1Q Median 3Q Max # -4.5555 -0.9939 -0.1549 0.9779 5.0306 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 2.3145 0.1603 14.435 < 2e-16 *** # Composite.treatment40larvae 1.2409 0.7674 1.617 0.108 # Composite.treatment42larvae -4.9200 0.2868 -17.156 < 2e-16 *** # Composite.treatment40pupae -1.8646 0.2443 -7.633 1.42e-12 *** # Composite.treatment40immature -2.4430 0.2592 -9.424 < 2e-16 *** # Composite.treatment42immature -3.8078 0.2676 -14.227 < 2e-16 *** # Composite.treatment42mature -1.9737 0.2468 -7.996 1.68e-13 *** # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # (Dispersion parameter for quasibinomial family taken to be 2.737536) # # Null deviance: 2555.70 on 181 degrees of freedom # Residual deviance: 522.15 on 175 degrees of freedom # AIC: NA # estimate logit look ok # 30 exp(2.3145)/ (1+ exp(2.3145)) #0.9100708 1/(1+1/(exp(2.3145))) # 0.8972375 # crawley (2013) # 40 larvae exp(2.3145+1.2409)/ (1+ exp(2.3145+1.2409)) # 0.9722236 # 42 larvae exp(2.3145-4.9200)/ (1+ exp(2.3145-4.9200)) # 0.06878529 # 40 pupae exp(2.3145-1.8646)/ (1+ exp(2.3145-1.8646)) # 0.6106155 # 40 immature exp(2.3145-2.4430)/ (1+ exp(2.3145-2.4430)) # 0.4679191 # 42 immature exp(2.3145-3.8078)/ (1+ exp(2.3145-3.8078)) # 0.1834269 # 42 mature exp(2.3145-1.9737)/ (1+ exp(2.3145-1.9737)) # 0.5843848 describeBy(survivalstage$Prop.survivors, survivalstage$Composite.treatment,mat=TRUE) # item group1 vars n mean sd median trimmed mad min max range skew kurtosis se # X11 1 30 1 33 0.8972727 0.06587264 0.90 0.9022222 0.074130 0.75 1.00 0.25 -0.6151304 -0.0196779 0.01146695 # X12 2 larvae 1 9 0.1822222 0.10109127 0.21 0.1822222 0.059304 0.00 0.32 0.32 -0.5494953 -1.0709996 0.03369709 # X13 3 pupae 1 5 0.0000000 0.00000000 0.00 0.0000000 0.000000 0.00 0.00 0.00 NaN NaN 0.00000000 # X14 4 immature 1 17 0.2458824 0.19789703 0.20 0.2200000 0.148260 0.00 0.88 0.88 1.7212923 3.4720346 0.04799708 # X15 5 42 1 11 0.4727273 0.24120908 0.45 0.4833333 0.370650 0.10 0.75 0.65 -0.1808827 -1.6739800 0.07272727 library(lsmeans) lsmeans(globalmodquasibin, pairwise~Composite.treatment, adjust="tukey") # contrast estimate SE df z.ratio p.value # 30control - 40larvae -1.2408682 0.7673721 NA -1.617 0.6713 = # 30control - 42larvae 4.9200145 0.2867809 NA 17.156 <.0001 = # 30control - 40pupae 1.8645630 0.2442818 NA 7.633 <.0001 = # 30control - 40immature 2.4429582 0.2592333 NA 9.424 <.0001 = # 30control - 42immature 3.8078456 0.2676444 NA 14.227 <.0001 = # 30control - 42mature 1.9737204 0.2468442 NA 7.996 <.0001 = # 40larvae - 42larvae 6.1608827 0.7871998 NA 7.826 <.0001 # 40larvae - 40pupae 3.1054312 0.7727314 NA 4.019 0.0011 # 40larvae - 40immature 3.6838264 0.7775874 NA 4.738 <.0001 # 40larvae - 42immature 5.0487138 0.7804318 NA 6.469 <.0001 # 40larvae - 42mature 3.2145886 0.7735453 NA 4.156 0.0006 # 42larvae - 40pupae -3.0554515 0.3008275 NA -10.157 <.0001 # 42larvae - 40immature -2.4770563 0.3130902 NA -7.912 <.0001 # 42larvae - 42immature -1.1121689 0.3200892 NA -3.475 0.0092 # 42larvae - 42mature -2.9462941 0.3029119 NA -9.727 <.0001 # 40pupae - 40immature 0.5783952 0.2746922 NA 2.106 0.3492 # 40pupae - 42immature 1.9432826 0.2826436 NA 6.875 <.0001 # 40pupae - 42mature 0.1091574 0.2630322 NA 0.415 0.9996 # 40immature - 42immature 1.3648874 0.2956615 NA 4.616 0.0001 # 40immature - 42mature -0.4692378 0.2769734 NA -1.694 0.6200 # 42immature - 42mature -1.8341252 0.2848612 NA -6.439 <.0001 coefficients(globalmodquasibin) # (Intercept) Composite.treatment40larvae Composite.treatment42larvae Composite.treatment40pupae # 2.314480 1.240868 -4.920015 -1.864563 # Composite.treatment40immature Composite.treatment42immature Composite.treatment42mature # -2.442958 -3.807846 -1.973720 summary(globalmodquasibin) confint(globalmodquasibin) # 2.5 % 97.5 % # (Intercept) 2.012956621 2.643082 # Composite.treatment40larvae -0.009739147 3.175392 # Composite.treatment42larvae -5.509822725 -4.381877 # Composite.treatment40pupae -2.347423130 -1.387851 # Composite.treatment40immature -2.958229479 -1.940120 # Composite.treatment42immature -4.350367676 -3.299014 # Composite.treatment42mature -2.462003282 -1.492420 ############## ! ONE-SAMPLE COMPARISONS TO PUPAE ########################### ## comparing all treatments to 42 pupae (median being 0) # wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="30control"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 2346, p-value = 5.756e-13 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="40larvae"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 45, p-value = 0.007745 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="42larvae"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 231, p-value = 5.134e-05 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="40pupae"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 153, p-value = 0.0003149 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="41pupae"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 6, p-value = 0.1489 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="40immature"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 55, p-value = 0.005889 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="42immature"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 136, p-value = 0.0004627 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="42mature"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 171, p-value = 0.0002119 ## comparing all treatments to 41 pupae (median being 0.05) wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="30control"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 2346, p-value = 5.756e-13 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="40larvae"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 45, p-value = 0.007745 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="42larvae"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 236, p-value = 0.4418 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="40pupae"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 153, p-value = 0.0003149 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="42pupae"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 0, p-value = 0.03689 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="40immature"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 45, p-value = 0.009091 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="42immature"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 201, p-value = 0.002836 wilcox.test(survivalstage$Prop.survivors[survivalstage$Composite.treatment=="42mature"], mu = 0.05, exact = TRUE, conf.int = TRUE) # V = 171, p-value = 0.0002119 #################### males only check ################################## library(car); library(MASS); library (lme4); library (nlme) ###################### MODEL SELECITON ### MODEL REFINEMENT globalmodquasibin<-glm(Cbindsurvival ~ Temperature.oC, quasibinomial(link = "logit"), data=survivalstagemalenopupae); summary(globalmodquasibin) # RD 135.33 pseudoR<-(globalmodquasibin$null.deviance-globalmodquasibin$deviance) / globalmodquasibin$null.deviance; pseudoR # 0.809555 nullmodquasibin<-glm(Cbindsurvival ~ 1, quasibinomial(link = "logit"), data=survivalstagemalenopupae); summary(nullmodquasibin) # 754.52 pseudoR<-(nullmodquasibin$null.deviance-nullmodquasibin$deviance) / nullmodquasibin$null.deviance; pseudoR # 0 ### Minimum adequate model explains offspring significantly more than null model because # anova/lrtest(model, nullmodel) and aic functions do not work with quasibinomial anova(globalmodquasibin, test = "Chisq") ###################### MODEL SIGNIFICANCE drop1(globalmodquasibin, test = "Chisq") # Model: # Cbindsurvival ~ Temperature.oC # Df Deviance scaled dev. Pr(>Chi) # 70.84 # Temperature.oC 3 371.95 170.47 < 2.2e-16 *** anova(globalmodquasibin, nullmodquasibin, test = "Chisq") # Resid. Df Resid. Dev Df Deviance Pr(>Chi) # 1 37 70.84 # 2 40 371.95 -3 -301.12 < 2.2e-16 *** ###################### MODEL POST HOC summary(globalmodquasibin) # Call: # glm(formula = Cbindsurvival ~ Temperature.oC, family = quasibinomial(link = "logit"), # data = survivalstagemalenopupae) # # Deviance Residuals: # Min 1Q Median 3Q Max # -2.7614 -0.5739 0.2938 0.9253 3.5041 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 2.1102 0.2660 7.935 1.68e-09 *** # Temperature.oClarvae -3.6143 0.3717 -9.722 9.83e-12 *** # Temperature.oCimmature -3.0366 0.3354 -9.053 6.43e-11 *** # Temperature.oC42 -1.8383 0.3455 -5.321 5.21e-06 *** # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # # (Dispersion parameter for quasibinomial family taken to be 1.76638) # # Null deviance: 371.953 on 40 degrees of freedom # Residual deviance: 70.837 on 37 degrees of freedom # AIC: NA # # Number of Fisher Scoring iterations: 4 library(lsmeans) lsmeans(globalmodquasibin, pairwise~Temperature.oC, adjust="tukey") # $contrasts # contrast estimate SE df z.ratio p.value # 30 - larvae 3.6142906 0.3717460 NA 9.722 <.0001 # 30 - immature 3.0366221 0.3354449 NA 9.053 <.0001 # 30 - 42 1.8382795 0.3454819 NA 5.321 <.0001 # larvae - immature -0.5776685 0.3305433 NA -1.748 0.2990 # larvae - 42 -1.7760111 0.3407248 NA -5.212 <.0001 # immature - 42 -1.1983426 0.3007018 NA -3.985 0.0004 # comparisons of other lifestages to pupae wilcox.test(survivalstagemalenopupae$Prop.survivors[survivalstagemalenopupae$Temperature.oC=="30"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 91, p-value = 0.001524 wilcox.test(survivalstagemalenopupae$Prop.survivors[survivalstagemalenopupae$Temperature.oC=="larvae"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 36, p-value = 0.01403 wilcox.test(survivalstagemalenopupae$Prop.survivors[survivalstagemalenopupae$Temperature.oC=="immature"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 66, p-value = 0.003667 wilcox.test(survivalstagemalenopupae$Prop.survivors[survivalstagemalenopupae$Temperature.oC=="42"], mu = 0, exact = TRUE, conf.int = TRUE) # V = 36, p-value = 0.01403 ### in base # control hist(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "30"], main = list("Control", cex = 2), xlab = "survival", ylab ="Frequency", ylim = c(0,20), nclass = 10) # larvae hist(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "larvae"], col = "red", density = 30, angle = 180, border = "red", main = list("larvae", cex = 2), xlab = "survival", ylab ="Frequency", ylim = c(0,20), nclass = 10) # pupae hist(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "pupae"], col = "red", density = 30, angle = 180, border = "red", main = list("pupae", cex = 2), xlab = "survival", ylab ="Frequency", ylim = c(0,20), nclass = 10) # immature hist(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "immature"], col = "red", density = 30, angle = 180, border = "red", main = list("immature", cex = 2), xlab = "survival", ylab ="Frequency", ylim = c(0,20), nclass = 10) # mature hist(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "42"], col = "red", density = 30, angle = 180, border = "red", main = list("42", cex = 2), xlab = "survival", ylab ="Frequency", ylim = c(0,20), nclass = 10) # keep nclass = 10, keep scales default ###### plotting differences # base boxplots of data distribution grouped by temperature boxplot(survivalstage$Prop.survivors ~ survivalstage$Temperature.oC, ylab="survival", xlab="Temperature") # notice plot has automatically produced a scatterplot # if the mean.sperm.count is made as an integar ########### Normality - Failed, especially in one group with no variation shapiro.test (survivalstage$Prop.survivors[survivalstage$Temperature.oC == "30"]) # W = 0.9002, p-value = 0.00539 ks.test(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "30"], pnorm) # D = 0.77337, p-value < 2.2e-16 shapiro.test (survivalstage$Prop.survivors[survivalstage$Temperature.oC == "larvae"]) # W = 0.91828, p-value = 0.3782 ks.test(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "larvae"], pnorm) # D = 0.5, p-value = 0.02222 shapiro.test (survivalstage$Prop.survivors[survivalstage$Temperature.oC == "pupae"]) # Error ks.test(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "pupae"], pnorm) # D = 0.5, p-value = 0.1641 shapiro.test (survivalstage$Prop.survivors[survivalstage$Temperature.oC == "immature"]) # W = 0.80155, p-value = 0.00213 ks.test(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "immature"], pnorm) # D = 0.5, p-value = 0.0004069 shapiro.test(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "42"]) # W = 0.90239, p-value = 0.1977 ks.test(survivalstage$Prop.survivors[survivalstage$Temperature.oC == "42"], pnorm) # D = 0.53983, p-value = 0.003286 ########### Homogeneity of Variances - Failed, one group with no variation bartlett.test(survivalstage$Prop.survivors ~ survivalstage$Temperature.oC) # Bartlett's K-squared = Inf, df = 4, p-value < 2.2e-16 fligner.test(survivalstage$Prop.survivors ~ survivalstage$Temperature.oC) # Fligner-Killeen:med chi-squared = 27.264, df = 4, p-value = 1.758e-05 #! need library(car) leveneTest(survivalstage$Prop.survivors ~ survivalstage$Temperature.oC) # 4 7.6797 3.442e-05 *** ##### complete separation will intefere with statistical testing, excluding from analysis ########### survivalstagenopupae<-subset(survivalstage, !(Temperature.oC %in% c("pupae"))) survivalstagenopupae$Temperature.oC<-droplevels(survivalstagenopupae$Temperature.oC) str(survivalstagenopupae) ######################### male only stage library(psych) describeBy(survivalstagemale$Prop.survivors, survivalstagemale$Temperature.oC) # $`30` # vars n mean sd median trimmed mad min max range skew kurtosis se # X1 1 13 0.89 0.08 0.9 0.9 0.07 0.75 1 0.25 -0.6 -0.95 0.02 # # $larvae # vars n mean sd median trimmed mad min max range skew kurtosis se # X1 1 9 0.18 0.1 0.21 0.18 0.06 0 0.32 0.32 -0.55 -1.07 0.03 # # $pupae # vars n mean sd median trimmed mad min max range skew kurtosis se # X1 1 5 0 0 0 0 0 0 0 0 NaN NaN 0 # # $immature # vars n mean sd median trimmed mad min max range skew kurtosis se # X1 1 11 0.32 0.21 0.3 0.28 0.07 0.05 0.88 0.83 1.49 1.94 0.06 # # $`42` # vars n mean sd median trimmed mad min max range skew kurtosis se # X1 1 8 0.58 0.18 0.65 0.58 0.15 0.3 0.75 0.45 -0.41 -1.72 0.06 # Data are very unlikely to be fixed with transformations # As the data is both not normal and not homogenous in variance in groups there is debate over the best method so both tried #1) Homogeneity of variance not strict assumpation of Kurskall wallace kruskal.test(survivalstagenopupae$Prop.survivors ~ survivalstagenopupae$Temperature.oC) # Kruskal-Wallis chi-squared = 51.966, df = 3, p-value = 3.045e-11 # there is a significant difference between groups # Posthoc testing #! library(pgirmess) kruskalmc(survivalstagenopupae$Prop.survivors ~ survivalstagenopupae$Temperature.oC, probs = 0.05, cont=NULL) # Multiple comparison test after Kruskal-Wallis # p.value: 0.05 # Comparisons # obs.dif critical.dif difference # 30-larvae 39.631313 20.19069 TRUE # 30-immature 36.222816 16.02907 TRUE # 30-42 26.621212 18.69295 TRUE # larvae-immature 3.408497 22.13327 FALSE # larvae-42 13.010101 24.13249 FALSE # immature-42 9.601604 20.77605 FALSE # output seems unrealistic #! library(PMCMR) posthoc.kruskal.nemenyi.test(survivalstagenopupae$Prop.survivors ~ survivalstagenopupae$Temperature.oC, dist="Chisquare") #"Tukey" for no ties # 30 larvae immature # larvae 5.8e-06 - - # immature 8.2e-08 0.9828 - # 42 0.0026 0.5646 0.6828 dunn.test.control(x=survivalstagenopupae$Prop.survivors,g=survivalstagenopupae$Temperature.oC, p.adjust="bonferroni") dunn.test.control(x=survivalstagenopupae$Prop.survivors,g=survivalstagenopupae$Temperature.oC) # 30 # larvae 6.1e-07 # immature 6.5e-09 # 42 0.00049 # 30 # larvae 4.0e-07 # immature 6.5e-09 # 42 0.00016 posthoc.kruskal.dunn.test(x=survivalstagenopupae$Prop.survivors, g=survivalstagenopupae$Temperature.oC, p.adjust.method="bonferroni") posthoc.kruskal.dunn.test(x=survivalstagenopupae$Prop.survivors, g=survivalstagenopupae$Temperature.oC) # 30 larvae immature # larvae 1.2e-06 - - # immature 1.3e-08 1.00000 - # 42 0.00098 0.92056 1.00000 # 30 larvae immature # larvae 1.0e-06 - - # immature 1.3e-08 0.68344 - # 42 0.00065 0.46028 0.46028 #2) Convert data to ranks and submitt to a welch anova and tukey HSD Ruxton (2006) survivalstagenopupae$rank<-rank(survivalstagenopupae$Prop.survivors) # ranking all data inter-group by ascending count oneway.test(survivalstagenopupae$rank ~ survivalstagenopupae$Temperature.oC, var.equal=FALSE) # F = 71.407, num df = 3.000, denom df = 24.719, p-value = 2.583e-12 TukeyHSD(aov) # only works with aov not oneway, #! library(userfriendlyscience) posthocTGH(y=survivalstagenopupae$rank, x=survivalstagenopupae$Temperature.oC, method="games-howell") # use games-howell when different sample sizes # tukey for equal # n means variances # 30 33 54 99 # larvae 9 14 60 # immature 17 17 143 # 42 11 27 108 # # t df p # 30:larvae 12.78 16 4.4e-09 # 30:immature 10.73 28 1.4e-10 # 30:42 7.45 17 6.1e-06 # larvae:immature 0.88 23 8.2e-01 # larvae:42 3.21 18 2.3e-02 # immature:42 2.25 24 1.4e-01