### ## Elastomere Tagging Experiment ## ## Chloe Fouilloux, Guillermo Garcia, and Bibiana Rojas ## #Code for Survival and Growth # #first let's load some packages# library(ggplot2) library(dplyr) library(lme4) library(lmerTest) library(survival) library(survminer) library(lsmeans) library(pbkrtest) library(coxme) #### #Now, let's read in our most recent updated data file. elastomere<- read.csv("Elastomer_Final.csv") View(elastomere) ############# ###subsets### ############# tagged <- elastomere[which(elastomere$Treatment=='Tag'), ] #Just tagged individuals tagged_na <- tagged[!is.na(tagged$Tag_Y_N), ] ##Tagged individuals with NAs dropped threemonths<- elastomere[ which(elastomere$Week < 13), ] #subset larval time tagged_na$Tag_0_1<- ifelse(tagged_na$Tag_Y_N == "Y", 1, 0) ## ################################## ###Addressing Reviewer Comments### ################################# #Comment (27) : #Did the same person observe the larvae at each every time-point? #Without knowing the answer, it is hard to know how confident we should #be in a false negative—were false negatives the result of observer error #or the mark physically migrating during development. table(elastomere$Observer) # CF EA NK # 315 69 114 315/(315+69+114) #63% ggplot(elastomere, aes(x = Observer))+ geom_bar(stat = "count")+ theme_bw() #Comment (32): Overstating the success in post-metamorphic individuals meta<- elastomere[ which(elastomere$Week > 13 & elastomere$Treatment == "Tag"), ] #subset individuals who survived past two weeks table(meta$ID) #11 individuals table(meta$ID, meta$Tag_Y_N, meta$Treatment == "Tag") #only 4 t<- ggplot(meta, aes(y = ID, x = Week, color = Tag_Y_N))+ geom_point(size = 3)+ ylab("Individual Tapdole Identity")+ xlab("Week of Development")+ labs(title = "Tag observation post-metamorphosis")+ theme_bw() t t + scale_color_manual(values=c("grey", "#0eb300"), name = "Tag observation", labels=c("No", "Yes", "") ) ############################################################ ######### Survival ####################################### ########################################################### #### Below two graphs good at showing death across control and tagged groups ggplot(elastomere, aes(x = Week, y = Weight, color = Dead))+ geom_point()+ facet_wrap(~Treatment) ggplot(elastomere, aes(x = Week, y = Weight, color = Dead))+ geom_point(aes(shape = Treatment))+ theme_bw() ##Binomial Response Plot to Tag Survival Across all time #Remember that this is skewed later in development because of all of our froglet death. ggplot(elastomere, aes(x = Week, y = Dead, color = Treatment))+ geom_point()+ geom_smooth(method = "glm", method.args = list(family = "binomial"),aes(fill = Treatment), se = T, fullrange = T)+ theme_bw() #If we subset this to just look at the first 3 months of development (larval) ##Binomial Response Plot to Tag Survival Across all time #This plot is to look at all of the invidividuals tested, see whose dead and #if I've entered any data incorrectly. ggplot(threemonths, aes(x = Week, y = Dead, color = Treatment))+ geom_point()+ geom_smooth(method = "glm", formula=y~(x), method.args = list(family = "binomial"),aes(fill = Treatment), se = F, fullrange = T)+ ylab("Larval death rate")+ geom_hline(yintercept = 0.1, linetype = "dashed", alpha = 0.3)+ theme_bw()+ scale_x_continuous(breaks = seq(0, 12, 4)) ############################################################ ######### Calculating Growth Rate ######################### ########################################################### ##Calculating average size at time of tagging weight_cont<- elastomere[ which(elastomere$Week == 0 & elastomere$Treatment== "Control"), ] weight_tag<- elastomere[ which(elastomere$Week == 0 & elastomere$Treatment== "Tag"), ] mean(weight_cont$Weight, na.rm = T) ##0.099 min(weight_cont$Weight, na.rm = T) #0.0307 max(weight_cont$Weight, na.rm = T) #0.18 mean(weight_tag$Weight, na.rm = T) ##0.12 min(weight_tag$Weight, na.rm = T) #0.0318 max(weight_tag$Weight, na.rm = T) #0.36 #se sd(weight_cont$Weight, na.rm=TRUE) / sqrt(length(weight_cont$Weight[!is.na(weight_cont$Weight)])) #0.0145 sd(weight_tag$Weight, na.rm=TRUE) / sqrt(length(weight_tag$Weight[!is.na(weight_tag$Weight)])) #0.0193 ####################### ### overall growth rate ### ####################### growth_rate = elastomere %>% # first sort by week for each individual tested arrange(ID, Week) %>% group_by(ID) %>% #this is to make sure that the difference is just within an individualso that #last week of one does not become the first week of the next mutate(Diff_week = Week - lag(Week), # Difference in time (just in case there are gaps) Diff_growth = Weight - lag(Weight), # Difference in weight between weeks Rate = (Diff_growth / Diff_week)/lag(Weight)) # EDIT (16/): Removed percent, gives growth rate per week g/week #Visualizing growth rate# ##We can also have the growth rate in mg per day if we so desire. #For now, I am going to keep g/week because that is the scale we actually measured on. growth_rate_day = elastomere %>% # first sort by week for each individual tested arrange(ID, Week) %>% group_by(ID) %>% #this is to make sure that the difference is just within an individualso that #last week of one does not become the first week of the next mutate(Diff_week = Week - lag(Week), # Difference in time (just in case there are gaps) ##g/week = week/ 7days = 1000mg/g = 142.85 mg/day Diff_growth = Weight - lag(Weight), # Difference in weight between weeks Rate = (Diff_growth / Diff_week)/lag(Weight)* (1000/7)) ####################### #Growth rate percent across development# ################### weekly<- ggplot(growth_rate, aes(color = Treatment, y = Rate, x = Week))+ geom_point()+ stat_smooth(method = "glm", fullrange = F, formula = y ~ log(x))+ ylab("Growth rate (g/week) across development")+ scale_x_continuous(breaks = seq(0, 24, 4))+ theme_bw()+ theme(legend.position = "top", panel.grid.minor.x = element_blank()) weekly+ scale_color_manual(values = c("black", "#0eb300")) ############################################################## ###### MODELLING GROWTH RATE ################################### ############################################################# #Looking at weight is great and all, but what we really want is growth RATE from ## week to week. #not possion or gamma because of growth rate can be negative. ######################### #In statistics code when comparing linear mixed models, #each of the linear models should be fit with maximum likelihood ### ########################################################## #Because our best models were LMMs we should compare model fits with MLs! weight.sub<-lmer(Rate ~ Treatment + (1|ID), data = growth_rate, REML = F) #changed to ML weight_week<- lmer(Rate ~ Treatment* Week + (1|ID), data = growth_rate, REML = F) weightweek1<- lmer(Rate~ Treatment+ Week + (1|ID), data = growth_rate, REML = F) #best weight.one.week<- lmer(Rate ~ 1+ Week + (1|ID), data = growth_rate, REML = F) #best weight.one<-lmer(Rate ~ 1 + (1|ID), data = growth_rate, REML = F) weight.lm<-lm(Rate ~ Treatment, data = growth_rate) weight.glm<-glm(Rate ~ Treatment, data = growth_rate, family = "quasi") #trying something a little different #All LMM models now fit with ML. Let's see what happens! #AIC Model Comparison Showdown aic<-AIC(weight.lm, weight.sub, weight.one, weight_week, weightweek1, weight.one.week, weight.glm) aic aic[order(aic$AIC),] #Fascinating!! So our lowest models here are weight.one.week and weightweek, #which are weighted practically identically. This is really interesting because #it means that the effect of treatment is not important to describing growth rate. #This is great because it shows that growth rate is not described by treatment, which #means that treatment doesn't affect tadpole growth rate. Yay! anova(weightweek1, weight.one.week) #As we can see the treatment term is not significant, so it could be omitted from the model. #However, for the sake of clarity when graphing and because the deltaAIC <2, I have decided #to keep the model with the Treatment term because it is biologically relevant. #Finally, we change the weightweek back to REML to get final model outputs (Following Zuur 2009) weightweek<- lmer(Rate~ Treatment+ Week + (1|ID), data = growth_rate, REML = T) summary(weightweek) #weight and week as addative effects ##Here we see that tagged individuals have an intercept that is -3.24 lower than controls. ##and that across weeks, tadpoles grow more slowly (-1.5) #Let's look at how this pans out in an anova' anova(weightweek, ddf="Kenward-Roger") ################ ## REPORTING LMM MODEL OUTPUT ##Show there is no difference in growth between the two treatments visually :) . #################################### #### PLOTTING CONFIDENCE LIMITS OF GROWTH RATE##### ################################# lsm<-lsmeans(weightweek, "Treatment") ##here we go! lsm with(summary(lsmeans(weightweek,spec="Treatment")),cbind(lower.CL,upper.CL)) plot(lsm) marginal = lsmeans(lsm,~ Treatment) CLD<- multcomp::cld(marginal, alpha=0.05, Letters=letters, adjust="tukey") CLD cld<-ggplot(CLD, aes(x= Treatment, y= lsmean, color = Treatment)) + ylim(0, 0.25)+ geom_point(shape = 9, size = 4, position = position_dodge(0.4))+ geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2, size = 0.7, position = position_dodge(0.4))+ geom_text( aes( label = .group), nudge_y = 0, nudge_x = 0.1, fontface = "bold", size = 5.5)+ ylab("Growth rate (g/week) across development")+ xlab("")+ theme_bw()+ theme(legend.position= "none", axis.text=element_text(size=12), axis.title.y = element_text(size = 12)) cld+ scale_color_manual(values = c("black", "#0eb300")) #CLs between both treatments overlap, which indicates that there is no difference #between groups. ############################################################## ###### TAG EFFECT ON SURVIVAL ################################### ############################################################# #Let's construct a survival curve! alive <- with(elastomere, Surv(Week, Dead)) head(alive,80) alive_fit <- survfit(Surv(Week, Dead) ~ 1, data=elastomere) summary(alive_fit, times = c(0, 4, 8, 12)) alive_trt_fit <- survfit(Surv(Week, Dead) ~ Treatment, data=elastomere) summary(alive_trt_fit) ## So, here we are seeing the effects of our tadpoles dying later in life ## because of our poor care protocol-- though everyone is dying, they ## are dying equally across treartments, meaning that our tag isn't the thing thats' ## causing them to die. SO thats. . . better than nothing? ## Now we can also just look at this for larval developement, ## which is the threemonths data set that I had subsetted earlier alive_fit_larval <- survfit(Surv(Week, Dead) ~ Treatment, data=threemonths) ?ggsurvplot() #Great! Now let's get graphing kids! ggsurvplot(alive_fit_larval, data = threemonths, palette = c("grey", "#0eb300"), linetype = "strata", conf.int = T, size = 1, legend.title = (""), legend.labs = c(" Control", "Tag"), legend = "top", font.legend = c(12), font.x = c(13), font.y = c(13), xlab= ("Time (weeks)"), xlim= c(0,12), ylim = c(50, 100), break.time.by = 4, fun = "pct", ggtheme= theme(panel.grid.minor.x = element_blank(), panel.grid = element_line(colour = "#f3f3f2"), panel.background = element_rect(fill = "white",colour = "black"), axis.text=element_text(size=13))) ############################################## ## Fitting a Cox Proportional Hazards Model## ############################################## cox <- coxph(Surv(Week, Dead) ~ Treatment, data = elastomere) cox summary(cox) cox1 <- coxph(Surv(Dead) ~ Treatment, data = elastomere) summary(cox1) #I realize that I should be taking into account the random effect of #tadpole ID because a single tadpole was weighed/measured multiples times. #let's run this with a random effect and see how it compares to the previous code. coxmix<- coxme(Surv(Week, Dead) ~ Treatment + (1 |ID), data = elastomere) summary(coxmix) AIC(cox, cox1, coxmix)