# (a) Lower lethal limit analysis # Calculating LT50 and LT90 from the lower lethal limit experiments # The lethal temperatures that would kill 50% and 90% of the population respectively library(dplyr) library(ggplot2) theme_bw.2<-theme_bw() + theme(legend.direction="horizontal", legend.position=c(0.5,1.047), legend.title=element_blank()) library(MASS) #load dataset (a_sayi_survival.csv) sayi<- read.csv(file.choose()) sayi sayi.may <- sayi%>% filter(round == "may") sayi.may attach(sayi.may) y <- cbind(alive, dead) mod.sayi.may <- glm(y ~ temperature, binomial) summary(mod.sayi.may) dose.p(mod.sayi.may, p=0.5) #temperature at which 50% of the population dies = -3.4 degrees C +/- 0.2 dose.p(mod.sayi.may, p=0.1) #temperature at which 90% of the population dies = -4.7 degrees C +/- 0.4 sayi.june <- sayi %>% filter(round == "june") sayi.june detach(sayi.may) attach(sayi.june) x <- cbind(alive, dead) mod.sayi.june <- glm(x ~ temperature, binomial) summary(mod.sayi.june) dose.p(mod.sayi.june, p=0.5) #temperature at which 50% of the population dies = -2.3 degrees C +/- 0.1 dose.p(mod.sayi.june, p=0.1) #temperature at which 90% of the population dies = -3.1 degrees C +/- 0.2 #load data (a_orbicollis_survival.csv) orbicollis <- read.csv(file.choose()) orbicollis orb.june <- orbicollis %>% filter(round == "june") orb.june detach(sayi.june) attach(orb.june) a <- cbind(alive, dead) mod.orb.june <- glm(a ~ temperature, binomial) summary(mod.orb.june) dose.p(mod.orb.june, p=0.5) #temperature at which 50% of the population dies = -2.9 degrees C +/- 0.1 dose.p(mod.orb.june, p=0.1) #temperature at which 90% of the population dies = -3.6 degrees C +/- 0.2 orb.july <- orbicollis %>% filter(round == "july") orb.july detach(orb.june) attach(orb.july) b<- cbind(alive, dead) mod.orb.july <- glm(b ~ temperature, binomial) summary(mod.orb.july) dose.p(mod.orb.july, p=0.5) #temperature at which 50% of the population dies = -2.4 degrees C +/- 0.1 dose.p(mod.orb.july, p=0.1) #temperature at which 90% of the population dies = -3.1 degrees C +/- 0.2 detach(orb.july) # create dataframe for results for plotting (below) prop.dead <- c(0.5, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9, 0.9) species <- c("sayi", "sayi", "orbicollis", "orbicollis", "sayi", "sayi", "orbicollis", "orbicollis") species.round <- c("sayi_may", "sayi_june", "orbicollis_june", "orbicollis_july", "sayi_may", "sayi_june", "orbicollis_june", "orbicollis_july") temp <- c(-3.393381, -2.349882, -2.88248, -2.412882, -4.688036, -3.089246, -3.591371, -3.105457) temp.SE <- c(0.2129695, 0.1241203, 0.1207549, 0.1216382, 0.4142359, 0.2437746, 0.2463814, 0.2308183) LTe <- data.frame(prop.dead, species.round, temp, temp.SE) ###################################################### # Lower lethal limit analysis # A logistic regression with logit link and binomial error distribution # (i.e. a probit analysis) is, a priori, most appropriate (e.g. Bürgi and Mills (2010), Sinclair (1997)). #load data (a_lethal_limit_rawdata.csv) dat1 <- read.csv(file.choose()) names(dat1) head(dat1) dat1$species.round <- as.factor(dat1$species.round) #compare all combinations of predictors across models library(MuMIn) library(heatmapFit) options(na.action = "na.fail") #prevent fitting models to different datasets bmodel <- glm(survival ~ temperature * species.round * mass.g, family = binomial(link="logit"), dat1) pred.bmodel <- predict(bmodel, type="response") heatmap.fit(bmodel$y, pred.bmodel, reps=1000) best <-dredge(bmodel, rank = "AICc") best subset(best, delta<2) #subset of models where delta AICc < 2 best.model <- glm(survival ~ temperature + species.round, family=binomial(link = "logit"), dat1) pred.best.model <- predict(best.model, type="response") heatmap.fit(best.model$y, pred.best.model, reps=1000) summary(best.model) ######################### summary(best.model) # intercept = orbicollis_july # sayi round 1 vs orb round 3 z=3.92, p<0.0001 # sayi round 2 vs orb round 3 z=-0.29, p=0.77 # orb round 2 vs orb round 3 z=2.22, p=0.027 dat1$species.round2[dat1$species.round == "orbicollis_june"]<-"a_orbicollis_june" dat1$species.round2[dat1$species.round == "orbicollis_july"]<-"orbicollis_july" dat1$species.round2[dat1$species.round == "sayi_june"]<-"sayi_june" dat1$species.round2[dat1$species.round == "sayi_may"]<-"sayi_may" dat1$species.round2 best.model2 <- glm(survival ~ temperature + species.round2, family = binomial(link="logit"), dat1) summary(best.model2) # intercept = orbicollis_june # sayi round 1 vs orb round 2 z=2.26, p=0.024 # sayi round 2 vs orb round 2 z=-2.47, p=0.014 dat1$species.round3[dat1$species.round == "orbicollis_june"]<-"orbicollis_june" dat1$species.round3[dat1$species.round == "orbicollis_july"]<-"orbicollis_july" dat1$species.round3[dat1$species.round == "sayi_june"]<-"asayi_june" dat1$species.round3[dat1$species.round == "sayi_may"]<-"sayi_may" dat1$species.round3 best.model3 <- glm(survival ~ temperature + species.round3, family = binomial(link="logit"), dat1) summary(best.model3) # intercept = sayi_june # sayi round 1 vs sayi round 2 z=4.08, p<0.0001 ############################### # Plot model results pframe <- data.frame(species.round = "sayi_may", temperature = seq(-7, 2, by = .01)) pp <- predict(best.model, newdata = pframe, se.fit = TRUE) linkinv <- family(best.model)$linkinv ## inverse-link function pframe$pred0 <- pp$fit pframe$pred <- linkinv(pp$fit) alpha <- 0.95 sc <- abs(qnorm((1-alpha)/2)) ## Normal approx. to likelihood pframe <- transform(pframe, lwr=linkinv(pred0-sc*pp$se.fit), upr=linkinv(pred0+sc*pp$se.fit)) pred.sayimay <- ggplot(pframe, aes(x=temperature, y=pred)) + ggtitle("Nicrophorus sayi in May") + xlab("Temperature (degrees C)") + ylab("Survival") + geom_ribbon(aes(ymin = lwr, ymax = upr), fill="grey80") + geom_line(aes(y=pred), size = 1) pred.sayimay + theme_bw.2 pframe1 <- data.frame(species.round = "sayi_june", temperature = seq(-7, 2, by = .01)) pp1 <- predict(best.model, newdata = pframe1, se.fit = TRUE) linkinv <- family(best.model)$linkinv ## inverse-link function pframe1$pred2 <- pp1$fit pframe1$pred1 <- linkinv(pp1$fit) alpha <- 0.95 sc <- abs(qnorm((1-alpha)/2)) ## Normal approx. to likelihood pframe1 <- transform(pframe1, lwr1=linkinv(pred2-sc*pp1$se.fit), upr1=linkinv(pred2+sc*pp1$se.fit)) pred.sayijune <- ggplot(pframe1, aes(x=temperature, y=pred1)) + ggtitle("Nicrophorus sayi in June") + xlab("Temperature (degrees C)") + ylab("Survival") + geom_ribbon(aes(ymin = lwr1, ymax = upr1), fill="grey80") + geom_line(aes(y=pred1), size = 1) pred.sayijune + theme_bw.2 pframe2 <- data.frame(species.round = "orbicollis_june", temperature = seq(-7, 2, by = .01)) pp2 <- predict(best.model, newdata = pframe2, se.fit = TRUE) linkinv <- family(best.model)$linkinv ## inverse-link function pframe2$pred4 <- pp2$fit pframe2$pred3 <- linkinv(pp2$fit) alpha <- 0.95 sc <- abs(qnorm((1-alpha)/2)) ## Normal approx. to likelihood pframe2 <- transform(pframe2, lwr2=linkinv(pred4-sc*pp2$se.fit), upr2=linkinv(pred4+sc*pp2$se.fit)) pred.orbicollisjune <- ggplot(pframe2, aes(x=temperature, y=pred3)) + ggtitle("Nicrophorus orbicollis in June") + xlab("Temperature (degrees C)") + ylab("Survival") + geom_ribbon(aes(ymin = lwr2, ymax = upr2), fill="grey80") + geom_line(aes(y=pred3), size = 1) pred.orbicollisjune + theme_bw.2 pframe3 <- data.frame(species.round = "orbicollis_july", temperature = seq(-7, 2, by = .01)) pp3 <- predict(best.model, newdata = pframe3, se.fit=TRUE) linkinv <- family(best.model)$linkinv ## inverse-link function pframe3$pred6 <- pp3$fit pframe3$pred5 <- linkinv(pp3$fit) alpha <- 0.95 sc <- abs(qnorm((1-alpha)/2)) ## Normal approx. to likelihood pframe3 <- transform(pframe3, lwr3=linkinv(pred6-sc*pp3$se.fit), upr3=linkinv(pred6+sc*pp3$se.fit)) pred.orbicollisjuly <- ggplot(pframe3, aes(x=temperature, y=pred5)) + ggtitle("Nicrophorus orbicollis in July") + xlab("Temperature (degrees C)") + ylab("Survival") + geom_ribbon(aes(ymin = lwr3, ymax = upr3), fill="grey80") + geom_line(aes(y=pred5), size = 1) pred.orbicollisjuly + theme_bw.2 # combined graph p <- ggplot() + #SAYI MAY (Orange) geom_ribbon(data=pframe,aes(x=temperature, ymin = lwr, ymax = upr), fill="grey80", alpha=0.3) + geom_line(data=pframe, aes(x=temperature, y=pred), colour="#E69F00", size=1)+ #SAYI JUNE (BLACK) geom_ribbon(data=pframe1,aes(x=temperature, ymin = lwr1, ymax = upr1), fill="grey80", alpha=0.3) + geom_line(data=pframe1, aes(x=temperature, y=pred1), colour="black", size=1)+ #ORBICOLLIS JUNE (BLUE) geom_line(data=pframe2, aes(x=temperature, y=pred3), colour="#0072B2", size=1)+ geom_ribbon(data=pframe2,aes(x=temperature, ymin = lwr2, ymax = upr2), fill="grey80", alpha=0.3) + #ORBICOLLIS JULY (GREEN) geom_line(data=pframe3, aes(x=temperature, y=pred5), colour="#009E73", size=1)+ geom_ribbon(data=pframe3,aes(x=temperature, ymin = lwr3, ymax = upr3), fill="grey80", alpha=0.3) + #TITLES, LABELS ggtitle("Lower Lethal Limit") + xlab("Temperature (degrees C)") + ylab("Survival") p+ theme_bw.2 # alternative figure LTe$prop.dead <- as.factor(LTe$prop.dead) #Re-order species_round LTe$species.round <- factor(LTe$species.round, levels = c("sayi_may", "sayi_june", "orbicollis_june", "orbicollis_july")) LTe50.fig <- ggplot(LTe, aes(x=species.round, y=temp, colour=species, shape=prop.dead)) + scale_color_manual(values=c('orange1','mediumpurple4')) + geom_errorbar(aes(ymin=temp-temp.SE, ymax=temp+temp.SE), width=.15) + ylim(-5.2, -2) + geom_point(size=4) LTe50.fig + theme_bw()+theme(legend.position = c(0.87, 0.22)) # figure 1b #####################################################