#(S2) SUPPLEMENTARY MATERIAL 2
#R/JAGS CODE FOR CONDUCTING ACTIVITY ANALYSIS


#Kayla Blincow & Brice Semmens
#3/9/2020

#IMPORTANT:::
#If you want to run the last half of this code without re-running the JAGS model
#load the FinalJAGSRun.RData file into the enivronment
#If you run the full code, it will take a long time

#Code Output: a figure that shows boxplot of activitiy metrics for each hour 
#in a 24 hour time period


#set working directory
setwd("C:/Users/kmbli/OneDrive - UC San Diego/Nassau Depth/Data")


rm(list = ls())

#load necessary libraries 
library(lubridate)
library(tidyverse)
library(rjags)
library(mcmcplots)
library(boot)
library(coda)
library(MCMCvis)


#If skipping the data prep/JAGS model run...
load("FinalJAGSRun.RData")
#then skip to line 212


####Data prep####
#read in the data
DNC <- read.csv("DNC_detections.csv", header = T)

#create datetime column that combines date and time
DNC$datetime <- paste(DNC$Date, DNC$Time)

#convert that object into a time format that R will recognize
DNC$datetime <- as.POSIXct(DNC$datetime, tz = "EST", 
                           format = "%m/%d/%Y %H:%M:%OS")

#make sure that worked okay..
summary(DNC)

#filter my data so that I am only looking at spawning months and home receivers
#non-spawning months
dnsp <- mutate(DNC, Month = month(datetime)) %>% 
  filter(Month > 2)

#home receiver
#tag 61's home receiver has two separate codes, so we will convert those to 
#the same so we make sure to get all the detections we need
dnsp$hydro <- as.numeric(dnsp$hydro)
dnsp$hydro[dnsp$hydro==104923] <- 5682

HR <- dnsp %>% 
  filter(hydro == 5221 | hydro == 5217 | hydro == 8018 | hydro == 5682 | 
           hydro == 5213) %>% 
  group_by(tag, hydro) %>% 
  mutate(counter = n()) %>% 
  group_by(tag) %>% 
  filter(counter == max(counter))

#check that worked
HR %>% group_by(tag, hydro) %>% 
  summarize(n())
#these hydrophone values match up with what I identified as home hydros for 
#each fish


#Now let's create a column that is the difference in time and one that's the 
#difference in depth between detections for each tag number
HR1 <- HR %>% group_by(tag) %>% 
  mutate(diff_dep = depth - lag(depth, default = depth[1]),
         diff_time = datetime - lag(datetime, default = datetime[1]))
#convert out of stupid tibble format
HR1 <- as.data.frame(HR1)
#note difftime is in seconds

#going to need a column that is just the hour, so let's make that
HR1 <- HR1 %>% 
  mutate(hour = hour(datetime))

#finally we need a column for the activity metric (absolute value of change in 
#depth divided by the corresponding change in time)
#filtered for detections that are within an hour of each other
hour <- HR1 %>% 
  filter(diff_time < 3600) %>% #cuz there are 3600 seconds in an hour...
  group_by(tag) %>% 
  mutate(act = abs(diff_dep)/as.numeric(diff_time))

hour <- as.data.frame(hour)


####let's start by looking at detections that occur within an hour of each other####

#need hour and tag to be a factor not a numeric
hour$hour <- as.factor(hour$hour)
hour$tag <- as.factor(hour$tag)

#convert NaNs to 0s--they are truly 0s (0 depth change/0 time change)
hour$act[is.na(hour$act)] <- 0

##################################################
# for JAG Delta GLMM model (note you could probably do this with https://rdrr.io/rforge/nwfscDeltaGLM/man/fitDeltaGLM.html)
##################################################

# Data we need to pass to JAGS:
# N.hrs - number of categorical hours in analysis (24)
N.hrs <- 24
# N.ind - number of individual fish for which we have observations
N.ind <- length(unique(hour$tag))
# nObs - number of total observations
nObs <- dim(hour)[1]
# nMoveObs - number of observations with positive movements
nMoveObs <- length(which(hour$act > 0))
# binaryObs - vector of binary movement observations in form [0,1]
binaryObs <- ceiling(hour$act / 10e10)
# moveObs - vector of all positive movements (logged before passing to JAGS)
moveObs<-log(hour$act[which(hour$act > 0)])
# hrsBinom - vector of which hour each binaryObs belongs to
hrsBinom<-hour$hour
levels(hrsBinom)<-rep(1:length(unique(hrsBinom)))
# hrsMove - vector of which hour each moveObs belongs to
hrsMove<-hour$hour[which(hour$act > 0)]
levels(hrsMove)<-rep(1:length(unique(hrsMove)))
# tagBinom - vector of which hour each binaryObs belongs to
tagBinom<-hour$tag
levels(tagBinom)<-rep(1:length(unique(tagBinom)))
# tagMove - vector of which hour each moveObs belongs to
tagMove<-hour$tag[which(hour$act > 0)]
levels(tagMove)<-rep(1:length(unique(tagMove)))


sink("deltaglmm.jags")
cat("
    model { 
    # Priors and constraints
    for (i in 1:N.hrs){
      mean.mu[i] ~ dunif(-10, 10)    # Prior for avg. move by individuals given movement
      mean.p[i] ~ dnorm(0, 0.304)  # Prior avg. prop of movement by individuals (logit space)
      
      mean.mu.sig[i] ~ dunif(0,10) # SD of random movement effect for each hr
      mean.mu.sig2[i] <- pow(mean.mu.sig[i], 2)       
      mean.mu.tau[i] <- pow(mean.mu.sig[i], -2)
      
      mean.p.sig[i] ~ dunif(0,10) # SD of random prob of movement effect for each hr
      mean.p.sig2[i] <- pow(mean.p.sig[i], 2)       
      mean.p.tau[i] <- pow(mean.p.sig[i], -2)
    }
    
    for (i in 1:N.hrs){
      for (j in 1:N.ind){
        Ydev[i,j] ~ dnorm(mean.mu[i],mean.mu.tau[i]) #expected mean movement for each ind in each hour
        logitpYdev[i,j] ~ dnorm(mean.p[i],mean.p.tau[i]) #expected prob of move for each ind in each hour in logit space
        logit(pYdev[i,j]) <-logitpYdev[i,j] #turn into proportion space (inverse of logit)
        resid.sigma[i,j] ~ dunif(0,10) #SD of movements for each ind in each hour
        resid.sigma2[i,j] <- pow(resid.sigma[i,j], 2)       
        resid.tau[i,j] <- pow(resid.sigma[i,j], -2)
      }
    }
    
    # Likelihood
    
    # likelihood for zero/nonzero movement
    for (i in 1:nObs){
      binaryObs[i]~dbern(pYdev[hrsBinom[i],tagBinom[i]])
    }
    
    #likelihood for positive movements
    for (i in 1:nMoveObs){
      moveObs[i]~dnorm(Ydev[hrsMove[i],tagMove[i]],resid.tau[hrsMove[i],tagMove[i]])
    }
    }
    ",fill = TRUE)
sink()


## Format data for JAGS
jags.data <- list(N.hrs=N.hrs,N.ind=N.ind,nObs=nObs,nMoveObs=nMoveObs,binaryObs=binaryObs,moveObs=moveObs,
                  hrsBinom=hrsBinom,hrsMove=hrsMove,tagBinom=tagBinom,tagMove=tagMove)


#run the models using rjags
dGlmm2 <- jags.model("deltaglmm.jags", data = jags.data, n.chains = 3,
                     n.adapt = 50000)

#set number of iterations and burn in
ni <- 300000
nb <- 50000

#pull samples from the model
samples <- coda.samples(dGlmm2,
             c('mean.mu', 'mean.p'), thin = 50,
             n.iter = ni)

#look at trace, density, and autocorrelation plots
mcmcplot(samples)

###Start here if you want to skip running the whole model###
#drop the burn in samples (with thinning we are left with 2000 samples)
samples_keep <- window(samples, start = nb, end = ni)


#look at Rhat values for the model parameters
MCMCsummary(samples_keep)

post.p <- inv.logit(MCMCchains(samples_keep, params = "mean.p"))
post.mu <- exp(MCMCchains(samples_keep, params = "mean.mu"))

meanps <- pivot_longer(as.data.frame(post.p), cols = starts_with("mean.p"), 
             names_to = "Hour", values_to = "PrMove")

meanmus <- pivot_longer(as.data.frame(post.mu), cols = starts_with("mean.mu"), 
                        names_to = "Hour", values_to = "MagMove")



#let's plot those estimates using the same format as my other figures
meanps$Hour <- factor(meanps$Hour, 
                     levels = c(
                       "mean.p[1]", "mean.p[2]", "mean.p[3]", "mean.p[4]", 
                       "mean.p[5]", "mean.p[6]", "mean.p[7]", "mean.p[8]",
                       "mean.p[9]", "mean.p[10]", "mean.p[11]", "mean.p[12]",
                       "mean.p[13]", "mean.p[14]", "mean.p[15]", "mean.p[16]",
                       "mean.p[17]", "mean.p[18]", "mean.p[19]", "mean.p[20]",
                       "mean.p[21]", "mean.p[22]", "mean.p[23]", "mean.p[24]"),
                     labels = c("1", "2", "3", "4", "5", "6", "7",
                                "8", "9", "10", "11", "12", "13", 
                                "14", "15", "16", "17", "18", "19",
                                "20", "21", "22", "23", "24"))

#manually calculate boxplot outlier cutoffs (values > 1.5*IQR and < 3*IQR)
a <- meanps %>%
  group_by(Hour) %>%
  summarize(iqr = IQR(PrMove),
            quantile_25 = quantile(PrMove, .25),
            quantile_75 = quantile(PrMove, .75)) %>%
  mutate(upperlim = quantile_75 + 3*iqr,
         lowerlim = quantile_25 - 3*iqr) %>%
  right_join(meanps)


#limit our data to just values within the outlier range
meanps2 <- filter(a, PrMove < upperlim & PrMove > lowerlim)

#plot
p2 <- ggplot(meanps2, aes(x=Hour, y=PrMove)) + theme_classic() + 
    labs(x = "Hour", y = "Mean Probability of Vertical Activity", size = 4) +
  theme(
    axis.title = element_text(color="black", size=15), 
    axis.text.x = element_text(color = "black", size = 13, margin = margin(r = 1)),
    axis.text.y = element_text(color = "black", size = 13),
    axis.ticks.length = unit(0.25, "cm")) +
  geom_rect(xmin = 4.5, xmax = 6.5, ymin = -0.05, ymax = 1.05, alpha = 0.1,
            fill = "gray90") +
  geom_rect(xmin = 18.5, xmax = 20.5, ymin = -0.05, ymax = 1.05, alpha = 0.1,
            fill = "gray90") +
  geom_rect(xmin = 0, xmax = 4.5, ymin = -0.05, ymax = 1.05, alpha = 0.1,
            fill = "gray69") +
  geom_rect(xmin = 20.5, xmax = 24.5, ymin = -0.05, ymax = 1.05, alpha = 0.1,
            fill = "gray69") +
  geom_boxplot() +
  geom_vline(xintercept = 4.5, linetype = "twodash") +
  geom_vline(xintercept = 6.5, linetype = "twodash") +
  geom_vline(xintercept = 18.5, linetype = "twodash") +
  geom_vline(xintercept = 20.5, linetype = "twodash")

 
png("Figure6.png", 
    units="in", 
    width=8, 
    height=4, 
    pointsize=12, 
    res=400)
suppressWarnings(print(p2))
dev.off()


#do this same plot for mean movement magnitude

#create appropriate hour factor labels for the plot
meanmus$Hour <- factor(meanmus$Hour, 
                      levels = c(
                        "mean.mu[1]", "mean.mu[2]", "mean.mu[3]", "mean.mu[4]", 
                        "mean.mu[5]", "mean.mu[6]", "mean.mu[7]", "mean.mu[8]",
                        "mean.mu[9]", "mean.mu[10]", "mean.mu[11]", "mean.mu[12]",
                        "mean.mu[13]", "mean.mu[14]", "mean.mu[15]", "mean.mu[16]",
                        "mean.mu[17]", "mean.mu[18]", "mean.mu[19]", "mean.mu[20]",
                        "mean.mu[21]", "mean.mu[22]", "mean.mu[23]", "mean.mu[24]"),
                      labels = c("1", "2", "3", "4", "5", "6", "7",
                                 "8", "9", "10", "11", "12", "13", 
                                 "14", "15", "16", "17", "18", "19",
                                 "20", "21", "22", "23", "24")
)



#manually calculate boxplot outlier cutoffs (values > 1.5*IQR and < 3*IQR)
b <- meanmus %>% 
  group_by(Hour) %>% 
  summarize(iqr = IQR(MagMove),
            quantile_25 = quantile(MagMove, .25),
            quantile_75 = quantile(MagMove, .75)) %>% 
  mutate(upperlim = quantile_75 + 3*iqr,
         lowerlim = quantile_25 - 3*iqr) %>% 
  right_join(meanmus)


#limit our data to just values within the outlier range
meanmus2 <- filter(b, MagMove < upperlim & MagMove > lowerlim)

#plot
p3 <- ggplot(meanmus2, aes(x=Hour, y=MagMove)) + theme_classic()+
  labs(x = "Hour", y = "Mean Magnitude of Vertical Activity", size = 4) +
  theme(
    axis.title = element_text(color="black", size=15), 
    axis.text.x = element_text(color = "black", size = 13, margin = margin(r = 1)),
    axis.text.y = element_text(color = "black", size = 13),
    axis.ticks.length = unit(0.25, "cm")) +
   geom_rect(xmin = 4.5, xmax = 6.5, ymin = -0.05, ymax = .05, alpha = 0.1,
             fill = "gray90") +
   geom_rect(xmin = 18.5, xmax = 20.5, ymin = -0.05, ymax = .05, alpha = 0.1,
             fill = "gray90") +
   geom_rect(xmin = 0, xmax = 4.5, ymin = -0.05, ymax = .05, alpha = 0.1,
             fill = "gray69") +
   geom_rect(xmin = 20.5, xmax = 24.5, ymin = -0.05, ymax = .05, alpha = 0.1,
             fill = "gray69") +
  geom_boxplot() +
  scale_y_continuous(limits = c(0, 0.0125))+
  geom_vline(xintercept = 4.5, linetype = "twodash") +
  geom_vline(xintercept = 6.5, linetype = "twodash") +
  geom_vline(xintercept = 18.5, linetype = "twodash") +
  geom_vline(xintercept = 20.5, linetype = "twodash")



png("Figure5.png", 
    units="in", 
    width=8, 
    height=4, 
    pointsize=12, 
    res=400)
suppressWarnings(print(p3))
dev.off()



#Done!