###
## 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)