library(ggplot2)
require(ggeffects)
library(sciplot)
require(bbmle)
require(effects)
require(glmmTMB) # model at the bottom.
require(MASS)
require(emmeans)
require(lme4)
library(grid)



hss.index <- 1
lss.index <- 2
lv.index <- 3
wtc.index <- 4

plot.colors <- c(	"#FF0000", #HSS
					"#1122CC", #LSS
					"#449911", #LV
					"#993399" #WTC
			)
fill.colors <- c(	"#FF000033", #HSS
					"#1122CC33", #LSS
					"#44991133", #LV
					"#99339933" #WTC
			)


text.pars <- gpar(cex=1.5, fontface="bold")

#Setting up data set and sex comb teeth plots

sexcomb <- read.csv("../data/sex_comb_data_uptogen24.csv", header = T)

str(sexcomb) #check that treatment is a factor, generation is an interger and change replicate to a factor
sexcomb$rep <- factor(sexcomb$rep)
sexcomb$treatment <- relevel(sexcomb$treatment, "WTC")
#sexcomb <- na.omit(sexcomb)
sexcomb$individualid <- factor(paste(sexcomb$treatment, sexcomb$rep, sexcomb$gen, sexcomb$id, sep="."))

####################################
#Primary comb
#First: run models for primary comb using glmmTMB, only HSS and LSS


#NOTES HERE:
#The model runs just fine with only HSS and LSS treatments
primary_evolve_model_glmm <- glmmTMB(primary_teeth ~ gen*treatment 
    + (1|treatment:rep) + (1|individualid), 
    family = poisson, 
   #ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS",])

summary(primary_evolve_model_glmm)
plot(allEffects(primary_evolve_model_glmm))
#there is clearly no effect of generation


#However, when we add the WTC or LV treatments, we get warnings about model convergence
# --either non-positive-definite Hessian matrix, or
# --extreme or very small eigenvalues detected;
# but in both cases, the actual effect estimates are about the same (2.406 for the intercept, essentially 0 for the others,
#  but the std errors for those estimates are NA)
primary_evolve_model_glmm <- glmmTMB(primary_teeth ~ gen*treatment 
    + (1|treatment:rep) + (1|individualid), 
    family = poisson, 
   #ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS" | sexcomb$treatment=="LV" ,])

summary(primary_evolve_model_glmm)
plot(allEffects(primary_evolve_model_glmm))



#YET!
# These warnings go away, and we get standard error estimates, when we remove the random effects from the model
#  so the convergence problem probably has to do with difficulty in estimating random effects?
# (Is it too difficult to estimate the rep effect if there are only 2 reps for a treatment, or if 
#   some of the reps are missing data? e.g., there are only
#   2 reps for WTC, and the 3rd rep for LV is also missing some of the generations)
primary_evolve_model_glmm <- glmmTMB(primary_teeth ~ gen*treatment 
    + (1|individualid), 
    family = poisson, 
   #ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS" | sexcomb$treatment=="WTC" | sexcomb$treatment=="LV",])

summary(primary_evolve_model_glmm)
plot(allEffects(primary_evolve_model_glmm))

#with just the other random effect, the model works too
primary_evolve_model_glmm <- glmmTMB(primary_teeth ~ gen*treatment 
    + (1|treatment:rep), 
    family = poisson, 
   #ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS" | sexcomb$treatment=="WTC" | sexcomb$treatment=="LV",])

summary(primary_evolve_model_glmm)
plot(allEffects(primary_evolve_model_glmm))

#in any case the numbers are almost identical no matter which random effect we keep,
# and there's clearly no real effect of anything for
# primary comb tooth number






#Plot

as.data.frame(effect(c("gen", "treatment"), primary_evolve_model_glmm ))

new_dat <- data.frame(gen = c(1:24, 1:24), 
    treatment = c(rep("LSS", 24), rep("HSS", 24), rep("WTC", 24), rep("LV", 24)),
    individualid = NA, rep = NA)

#Do predictions based on the model fit
predictTrial_glmmTMB <- predict(primary_evolve_model_glmm, 
   new_dat, interval = "confidence",
   allow.new.levels=TRUE, type = "response", se.fit = TRUE)


length(predictTrial_glmmTMB$fit)
head(predictTrial_glmmTMB$fit)
cbind(new_dat, predictTrial_glmmTMB)

length(predictTrial_glmmTMB$se)

summary(primary_evolve_model_glmm)    
fixef(primary_evolve_model_glmm)

#Run similar model, but using lme4

primary_evolve_model_glmm_lme4 <- glmer(primary_teeth ~ gen*treatment 
    + (1|individualid), 
    family = poisson, 
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS" | sexcomb$treatment=="WTC" | sexcomb$treatment=="LV",])
##plot(allEffects(primary_evolve_model_glmm_lme4))

as.data.frame(Effect(c("treatment", "gen"), primary_evolve_model_glmm_lme4))
    
summary(primary_evolve_model_glmm_lme4)

#take-home here is that they give almost identical results, in spite of possible warnings about not converging

#Make plots

predictions <- ggemmeans(primary_evolve_model_glmm, terms=c("gen", "treatment")) # on link scale. transform back via exp(). Very similar to using effects


myplot.1 <- ggplot(predictions, aes(x, predicted, color=group)) + 
 geom_line() + geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha=0.2, color=NA) + 
 scale_color_manual(values=c(plot.colors[wtc.index], plot.colors[hss.index], plot.colors[lss.index], plot.colors[lv.index])) + 
 labs(x="Generation", y="Number of primary sex comb teeth") + 
 scale_fill_manual(values=c(fill.colors[wtc.index], fill.colors[hss.index], fill.colors[lss.index], fill.colors[lv.index])) +
 theme(legend.title=element_blank(), panel.background=element_blank(), panel.grid.major=element_blank(), panel.border=element_rect(color="black", fill=NA)) +
 theme(plot.margin=unit(c(1,1,1,1), "cm"))

print(myplot.1) 

####################################
#Primary comb - gaps
#First: run models for primary comb gaps using glmmTMB, only HSS and LSS



#The model runs just fine with only HSS and LSS treatments
sexcomb$primary_gaps_binom <- ifelse(sexcomb$primary_gaps >= 1, 1, 0)
gaps_evolve_model_glmm <- glmmTMB(primary_gaps_binom ~ gen*treatment 
    + (1|treatment:rep) + (1|individualid), 
    family = binomial(link="logit"), 
#   ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS",])

summary(gaps_evolve_model_glmm)
plot(allEffects(gaps_evolve_model_glmm))
#according to the model there is an effect of generation; the probability of having a gap defect in the primary
# comb decreases a bit each generation
#the interaction is not significant; however, it LOOKS very different on the plot -- much steeper for HSS
plogis(-1.92944-24*(0.05359)) #estimate decline in defects across 24 generations?


#When we add the LV treatment...
gaps_evolve_model_glmm <- glmmTMB(primary_gaps_binom ~ gen*treatment 
    + (1|treatment:rep) + (1|individualid), 
    family = binomial(link="logit"), 
#   ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS" | sexcomb$treatment=="LV" ,])

summary(gaps_evolve_model_glmm)
plot(allEffects(gaps_evolve_model_glmm))
#largely the result is the same, but the LV main effect is significant
#it appears that the LV treatment has a lower probability of gaps overall


#And we add the WTC treatment...
gaps_evolve_model_glmm <- glmmTMB(primary_gaps_binom ~ gen*treatment 
    + (1|treatment:rep) + (1|individualid), 
    family = binomial(link="logit"), 
#   ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS" | sexcomb$treatment=="LV" | sexcomb$treatment=="WTC",])

summary(gaps_evolve_model_glmm)
plot(allEffects(gaps_evolve_model_glmm))
#the results are different! the generation main effect is gone
#the HSS effect is significant (higher probability of having gaps overall) and the
# LSS effect is marginally significant (p=0.07)
#again, it LOOKS like probability of gaps decreases in the HSS treatment, but the interaction
# term is not significant






#Plot

as.data.frame(effect(c("gen", "treatment"), gaps_evolve_model_glmm ))

new_dat <- data.frame(gen = c(1:24, 1:24), 
    treatment = c(rep("LSS", 24), rep("HSS", 24), rep("WTC", 24), rep("LV", 24)),
    individualid = NA, rep = NA)

#Do predictions based on the model fit
predictTrial_glmmTMB <- predict(gaps_evolve_model_glmm, 
   new_dat, interval = "confidence",
   allow.new.levels=TRUE, type = "response", se.fit = TRUE)


length(predictTrial_glmmTMB$fit)
head(predictTrial_glmmTMB$fit)
cbind(new_dat, predictTrial_glmmTMB)

length(predictTrial_glmmTMB$se)

summary(gaps_evolve_model_glmm)    
fixef(gaps_evolve_model_glmm)

#Run similar model, but using lme4

gaps_evolve_model_glmm_lme4 <- glmer(primary_gaps_binom ~ gen*treatment 
    + (1|individualid), 
    family = binomial(link="logit"), 
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS" | sexcomb$treatment=="WTC" | sexcomb$treatment=="LV",])
plot(allEffects(gaps_evolve_model_glmm_lme4))

as.data.frame(Effect(c("treatment", "gen"), gaps_evolve_model_glmm_lme4))
    
summary(gaps_evolve_model_glmm_lme4)

#Make plots

predictions <- ggemmeans(gaps_evolve_model_glmm, terms=c("gen", "treatment")) # are these predictions on link scale?
#they look like they are already probabilities

myplot.2 <- ggplot(predictions, aes(x, predicted, color=group)) + 
 geom_line() + geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha=0.2, color=NA) + 
 scale_color_manual(values=c(plot.colors[wtc.index], plot.colors[hss.index], plot.colors[lss.index], plot.colors[lv.index])) + 
 labs(x="Generation", y="Proportion with primary comb defects") + 
 scale_fill_manual(values=c(fill.colors[wtc.index], fill.colors[hss.index], fill.colors[lss.index], fill.colors[lv.index])) + 
 theme(legend.title=element_blank(), panel.background=element_blank(), panel.grid.major=element_blank(), panel.border=element_rect(color="black", fill=NA)) +
 theme(plot.margin=unit(c(1,1,1,1), "cm"))
print(myplot.2)
 



####################################
#Ectopic comb
#First: run models for ectopic comb using glmmTMB, only HSS and LSS

secondary_evolve_model_glmm <- glmmTMB(secondary_teeth ~ gen*treatment 
    + (1|treatment:rep) + (1|individualid), 
    family = poisson, 
   #ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS",])

summary(secondary_evolve_model_glmm)
plot(allEffects(secondary_evolve_model_glmm))
#take-home message here is:
# there is an effect of generation (i.e., ectopic sex comb tooth number decreases)
# but slope does not differ between HSS and LSS



secondary_LV_evolve_model_glmm <- glmmTMB(secondary_teeth ~ gen + (1|individualid), 
    family = poisson, 
   #ziformula = ~gen, # no evidence of zero inflation
    data = sexcomb[sexcomb$treatment=="LV",])

summary(secondary_LV_evolve_model_glmm)
plot(allEffects(secondary_LV_evolve_model_glmm))
#take-home message here is:
# no effect of gen when LV is considered on its own


as.data.frame(effect(c("gen", "treatment"), secondary_evolve_model_glmm ))

new_dat <- data.frame(gen = c(1:24, 1:24), 
    treatment = c(rep("LSS", 24), rep("HSS", 24)),
    individualid = NA, rep = NA)

#Do predictions based on the model fit
predictTrial_glmmTMB <- predict(secondary_evolve_model_glmm, 
   new_dat, interval = "confidence",
   allow.new.levels=TRUE, type = "response", se.fit = TRUE)


length(predictTrial_glmmTMB$fit)
head(predictTrial_glmmTMB$fit)
cbind(new_dat, predictTrial_glmmTMB)

length(predictTrial_glmmTMB$se)

summary(secondary_evolve_model_glmm)    
fixef(secondary_evolve_model_glmm)

#Run similar model, but using lme4

secondary_evolve_model_glmm_lme4 <- glmer(secondary_teeth ~ gen*treatment 
    + (1|treatment:rep) + (1|individualid), 
    family = poisson, 
    data = sexcomb[sexcomb$treatment=="HSS" | sexcomb$treatment=="LSS",])
plot(allEffects(secondary_evolve_model_glmm_lme4))

as.data.frame(Effect(c("treatment", "gen"), secondary_evolve_model_glmm_lme4))
    
summary(secondary_evolve_model_glmm_lme4)

#take home here is that they give almost identical results, in spite of possible warnings about not converging


#Make plots

predictions <- ggemmeans(secondary_evolve_model_glmm, terms=c("gen", "treatment")) # on link scale. transform back via exp(). Very similar to using effects


myplot.3 <- ggplot(predictions, aes(x, predicted, color=group)) + 
 geom_line() + geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha=0.2, color=NA) + 
 scale_color_manual(values=c(plot.colors[hss.index], plot.colors[lss.index])) + 
 labs(x="Generation", y="Number of ectopic sex comb teeth") + 
 scale_fill_manual(values=c(fill.colors[hss.index], fill.colors[lss.index])) +
 theme(legend.title=element_blank(), panel.background=element_blank(), panel.grid.major=element_blank(), panel.border=element_rect(color="black", fill=NA)) +
 theme(plot.margin=unit(c(1,1,1,1), "cm"))

 

print(myplot.3)

####################################
#Make actual PDF 
pdf(file="../outputs/Figure5_EE.pdf", width=7, height=15)

vp.1 <- viewport(x=0, y=1, width=1, height=1/3, just=c("left", "top"))
print(myplot.1, vp=vp.1)
grid.text(label="A", x=0.01, y=0.98, just=c("left", "top"), gp=text.pars, vp=vp.1)

vp.2 <- viewport(x=0, y=2/3, width=1, height=1/3, just=c("left", "top"))
print(myplot.2, vp=vp.2)
grid.text(label="B", x=0.01, y=0.98, just=c("left", "top"), gp=text.pars, vp=vp.2)

vp.3 <- viewport(x=0, y=1/3, width=1, height=1/3, just=c("left", "top"))
print(myplot.3, vp=vp.3)
grid.text(label="C", x=0.01, y=0.98, just=c("left", "top"), gp=text.pars, vp=vp.3)

dev.off()