### RMBL Project: Carter Adamson 2019
### 30 July 2019
### Amy & Carter
### Updated Sept 2019 by Amy

#setwd("~/Dropbox/AMY/Students/RMBL Undergraduate Student Projects/Adamson Carter")

setwd("~/Dropbox/AMY/Research/OTCs & Pollination")

library(tidyverse)
library(glmmTMB)
library(car) 
library(bbmle) # AIC tables
library(MASS) # for glm.nb
library(pscl) # for zero-inflated nb models


# Ben Boelker's function for checking for overdispersion; use to get a sense of whether there might be overdispersion
overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}


####################
##                ##
##      DATA      ##
##                ##
####################

# (1) load pollinator observation data
chamber_data <- read.csv("Adamson_chamber_data.csv", header=T)
chamber_data$treatment <- as.factor(chamber_data$treatment)
str(chamber_data)

#(2) load pollen deposition data
pollen_data <- read.csv("pollen deposition data.csv", header=T)

# create a dataframe with means for each plot
mean_pollen_data <- pollen_data %>% 
  dplyr::group_by(species, treatment, plot) %>% 
  dplyr::summarize(mean_pollen_count = mean(pollen_count))


#######################
##                   ##
##     ANALYSIS      ##
##                   ##
#######################

#### (1) What is the effect of OTC warming chambers on the number of visitors? ####

#Since the data are already 0s and 1s, try a binomial model: i.e. "What is the effect of OTCs on the probabilty of a visitor being present?"
mod1 <- glmmTMB(visitors ~ treatment * plant_species + (1|day), data=chamber_data, family="binomial")
mod2 <- glmmTMB(visitors ~ treatment + plant_species + (1|day), data=chamber_data, family="binomial")
modn <- glmmTMB(visitors ~ 1 + (1|day), data=chamber_data, family="binomial")

#Using AIC comparison/Chi-square tests, model 2 is the "best" model
AICctab(mod1,mod2,modn)
summary(mod2)



##### (2) What is the effect of OTC warming chambers on the visitation rate? #####

## First try to use all of the data (all models appear to be overdispersed)
chamber_data$visit_rate <- chamber_data$visits/chamber_data$num_flowers
rate.mod1 <- glmmTMB(visit_rate ~ treatment * plant_species + (1|day), weights = num_flowers, data=chamber_data, family="binomial")
overdisp_fun(rate.mod1) 

rate.mod2 <- glmmTMB(visit_rate ~ treatment + plant_species + (1|day), weights = num_flowers, data=chamber_data, family="binomial")
overdisp_fun(rate.mod2)

rate.modn <- glmmTMB(visit_rate ~ 1 + (1|day), weights = num_flowers, data=chamber_data, family="binomial")
overdisp_fun(rate.modn)

## Next try a betabinomial distribution
rate.mod1a <- update(rate.mod1, family="betabinomial")
rate.mod2a <- update(rate.mod2, family="betabinomial")
rate.modna <- update(rate.modn, family="betabinomial")

#Beta-binomial fits much better, indicating overdispersion
AICctab(rate.mod1,rate.mod1a) 
AICctab(rate.mod2,rate.mod2a)
AICctab(rate.modn,rate.modna)
AICctab(rate.mod1a, rate.mod2a, rate.modna)
summary(rate.mod2a)


##### (2b) for sessions with at least one visit, do rates differ between treatments? #####
cha_data_noz <- chamber_data %>% 
  dplyr::filter(visits>0)

# First try to use all of the data (success)
rate.modz1 <- glmmTMB(visit_rate ~ treatment * plant_species + (1|day), weights = num_flowers, data=cha_data_noz, family="binomial")
overdisp_fun(rate.modz1)

rate.modz2 <- glmmTMB(visit_rate ~ treatment + plant_species + (1|day), weights = num_flowers, data=cha_data_noz, family="binomial")
overdisp_fun(rate.modz2)

rate.modzn <- glmmTMB(visit_rate ~ 1 + (1|day), weights = num_flowers, data=cha_data_noz, family="binomial")
overdisp_fun(rate.modzn)

AICctab(rate.modz1, rate.modz2, rate.modzn) 
summary(rate.modzn) 

# Check beta-binomial here too; binomial is better fit
rate.modz1a <- update(rate.modz1, family="betabinomial")
rate.modz2a <- update(rate.modz2, family="betabinomial")
rate.modzna <- update(rate.modzn, family="betabinomial")
AICctab(rate.modz1,rate.modz1a) 
AICctab(rate.modz2,rate.modz2a)
AICctab(rate.modzn,rate.modzna)



##### (3) What is the effect of OTC warming chambers on pollen deposition? #####

## Use data with an average # conspecific grains per stigma for each plot
pollen_mod1_m <- lm(mean_pollen_count ~ treatment * species, data = mean_pollen_data)
plot(fitted(pollen_mod1_m), resid(pollen_mod1_m))
shapiro.test(resid(pollen_mod1_m)) # using here because low sample size, so any departures from normality should be reliable

pollen_mod2_m <- lm(mean_pollen_count ~ treatment + species, data = mean_pollen_data)
plot(fitted(pollen_mod2_m), resid(pollen_mod2_m))
shapiro.test(resid(pollen_mod2_m))

pollen_modn_m <- lm(mean_pollen_count ~ 1, data = mean_pollen_data)
shapiro.test(resid(pollen_modn_m))

AICctab(pollen_mod1_m, pollen_mod2_m, pollen_modn_m)
summary(pollen_mod1_m)
Anova(pollen_mod1_m, type='III')



#### calculate percent change in response variables
perc_change_visitation <- chamber_data %>% 
  dplyr::group_by(plant_species, treatment) %>% 
  dplyr::summarize(mean_visitors = mean(visitors, na.rm=T), 
                   mean_visit_rate = mean(visit_rate, na.rm=T))

(0.296 -0.0556)/0.296  #DeNu visitors, 81%
(0.481-0.0571)/0.481 #PoPu visitors, 88%
(0.0550 -0.00452)/0.0550 #DeNu visit rate, 92%
(0.0661-0.00995)/0.0661 #PoPu visit rate, 85%

perc_change_pollen <- pollen_data %>% 
  dplyr::group_by(species, treatment) %>% 
  dplyr::summarize(overall_mean_pollen = mean(pollen_count))

(23.9 -6.464)/23.9  #DeNu, 73%



##################
##              ##
##    FIGURES   ##
##              ##
##################

## FIGURE 1

layout(matrix(c(1, 2,
                3, 4,
                5, 6), nrow=3, byrow=TRUE))

# visitors
# make separate dataframes for plotting Delphinium and Potentilla results
chamber_data_d <- chamber_data %>% 
  dplyr::filter(plant_species == 'Delphinium')

chamber_data_p <- chamber_data %>% 
  dplyr::filter(plant_species == 'Potentilla')

mean_data_d <- chamber_data_d %>% 
  dplyr::group_by(plant_species, treatment) %>% 
  dplyr::summarize(mean_visitors = mean(visitors))

mean_data_p <- chamber_data_p %>% 
  dplyr::group_by(plant_species, treatment) %>% 
  dplyr::summarize(mean_visitors = mean(visitors))


barplot(mean_visitors ~ treatment, data = mean_data_d,
        col=("slateblue1"), ylab="Mean probability of the presence of a visitor", 
        ylim=c(0, 0.5), space = 0.8)

barplot(mean_visitors ~ treatment, data = mean_data_p,
        col=("darkgoldenrod1"), ylab ="", ylim=c(0, 0.5), space = 0.8)


# visitation rates
boxplot(visit_rate ~ treatment, data=chamber_data_d, outline=F,
        boxwex=0.5, whisklty = 1, whisklwd=3, staplelty = 0, staplewex=0, add=FALSE, medlwd=5,
        col=("slateblue1"), boxcol=c("slateblue1"), whiskcol=c("slateblue1"), medcol=c("darkorchid4"), 
        xlab="", ylab="Visits per flower (per observation session per day)", ylim=c(0, 0.25))

boxplot(visit_rate ~ treatment, data=chamber_data_p, outline=F,
        boxwex=0.5, whisklty = 1, whisklwd=3, staplelty = 0, staplewex=0, add=FALSE, medlwd=5,
        col=("darkgoldenrod1"), boxcol=c("darkgoldenrod1"), whiskcol=c("darkgoldenrod1"), medcol=c("darkgoldenrod4"), 
        xlab="", ylab="", ylim=c(0, 0.25))


# pollen
mean_pollen_data_d <- mean_pollen_data %>% 
  dplyr::filter(species=="delphinium")
mean_pollen_data_p <- mean_pollen_data %>% 
  dplyr::filter(species=="potentilla")

boxplot(mean_pollen_count ~ treatment, data=mean_pollen_data_d, outline=F,
        boxwex=0.5, whisklty = 1, whisklwd=3, staplelty = 0, staplewex=0, add=FALSE, medlwd=5,
        col=("slateblue1"), boxcol=c("slateblue1"), whiskcol=c("slateblue1"), medcol=c("darkorchid4"), 
        xlab="", ylab="# Pollen grains per stigma", ylim=c(0, 35))

boxplot(mean_pollen_count ~ treatment, data=mean_pollen_data_p, outline=F,
        boxwex=0.5, whisklty = 1, whisklwd=3, staplelty = 0, staplewex=0, add=FALSE, medlwd=5,
        col=("darkgoldenrod1"), boxcol=c("darkgoldenrod1"), whiskcol=c("darkgoldenrod1"), medcol=c("darkgoldenrod4"), 
        xlab="", ylab="", ylim=c(0, 35))


## FIGURE 2

boxplot(visit_rate ~ treatment, data=cha_data_noz, outline=F,
        boxwex=0.3, whisklty = 1, whisklwd=3, staplelty = 0, staplewex=0, add=FALSE, medlwd=5,
        col=("grey70"), boxcol=c("grey70"), whiskcol=c("grey70"), medcol=c("grey40"), 
        xlab="", ylab="Visits per flower (visitors present)", ylim=c(0, 0.5))