##ANTHROPOGENIC NOISE DECREASES ACTIVITY AND CALLING BEHAVIORS IN TWO SPECIES OF WILD MICE
##R version 3.2.2 (2015-03-09)

library(psych)
library(lattice)
library(Rmisc)
library(grid)
library(plyr)
library(ggplot2)
library(ggpubr)
library(nlme)
library(lme4)
library(MASS)
library(effsize)

#######Difference in Total Time Spent at the Focal Area by Deer Mice######
data<-read.csv ("Data_Activity.csv")

shapiro.test(data$Difference)

ggplot(data, aes(x=Treatment, y=Difference, fill=Treatment)) +
  geom_boxplot() +
  xlab("Treatment Type")+
  ylab("Difference in Total Time (minutes)")+
  scale_fill_manual(values = c("orange", "green"),labels=c("A" = "Anthropogenic Noise", "F" = "Natural Noise")) +
  ylim(-250,100)+   ### sets the y axis limit
  geom_jitter(position = position_jitter(0.1)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    #panel.border = element_blank(),
    axis.line = element_line(color = "black"),
    axis.line.x = element_line(color="black", size = 1),
    axis.line.y = element_line(color="black", size = 1),
    axis.title.x = element_text(colour="grey20",size=16,face="bold"),
    axis.title.y = element_text(colour="grey20",size=16,face="bold"), 
    axis.text.x = element_text(colour="grey20",size=12,face="bold"),
    axis.text.y = element_text(colour="grey20",size=12,face="bold")
  )

#Night Order
data$Night <- as.character(data$Night)
data$Night <- factor(data$Night, levels=unique(data$Night))

fit <- glmmPQL(DiffScale ~ Treatment + Night, 
               random = ~ 1|ID, family = "poisson", data=data)
summary(fit)

####Difference in Total Seeds Consumed by Both Species
data <- read.csv("Data_Foraging.csv")

shapiro.test(data$Difference)
shapiro.test(data$DiffScale)

#### night order
data$Night <- as.character(data$Night)
data$Night <- factor(data$Night, levels=unique(data$Night))

#### treatment order
data$Treatment <- as.character(data$Treatment)
data$Treatment <- factor(data$Treatment, levels=unique(data$Treatment))

fit <- glmmPQL(DiffScale ~ Treatment + Night, 
               random = ~ 1|ID, family = "poisson", data=data)
summary(fit)

tab<-aggregate(Difference ~ Treatment + Night ,data=data, FUN = "summary")
tab
sum = summarySE(data, 
                measurevar="Difference", 
                groupvars=c("Treatment", "Night"))
sum

#### Lefover husks 
data <- read.csv("Data_Husks.csv")

shapiro.test(data$Difference) 

data$Night <- as.character(data$Night)
data$Night <- factor(data$Night, levels=unique(data$Night))

data$Treatment <- as.character(data$Treatment)
data$Treatment <- factor(data$Treatment, levels=unique(data$Treatment))

fit <- glmmPQL(Difference ~ Treatment + Night, 
               random = ~ 1|ID, family = "poisson", data=data)
summary(fit)

tab<-aggregate(Difference ~ Treatment + Night ,data=data, FUN = "summary")
tab
sum = summarySE(data, 
                measurevar="Difference", 
                groupvars=c("Treatment", "Night"))
sum

####Difference in Latency to Enter the Focal Area and Time Spent in the Focal Area on the First Visit by Both Species
####Difference in Latency to Start Foraging and Time Spent Foraging by Both Species
data <- read.csv("Data_risk.csv")

shapiro.test(data$to) ###not normally distributed
shapiro.test(data$timeinFA) ###not normally distributed
shapiro.test(data$toGUD) ###not normally distributed
shapiro.test(data$timeinGUD)  ###normal distribution

data$Night <- as.character(data$Night)
data$Night <- factor(data$Night, levels=unique(data$Night))

data$Treatment <- as.character(data$Treatment)
data$Treatment <- factor(data$Treatment, levels=unique(data$Treatment))

fit <- glmmPQL(toScale ~ Treatment + Night, 
               random = ~ 1|ID, family = "poisson", data=data)
summary(fit)
tab<-aggregate(to ~ Treatment + Night ,data=data, FUN = "summary")
tab
sum = summarySE(data, 
                measurevar="to", 
                groupvars=c("Treatment", "Night"))
sum

fit <- glmmPQL(timeinFAScale ~ Treatment + Night, 
               random = ~ 1|ID, family = "poisson", data=data)
summary(fit)

tab<-aggregate(timeinFA ~ Treatment + Night ,data=data, FUN = "summary")
tab
sum = summarySE(data, 
                measurevar="timeinFA", 
                groupvars=c("Treatment", "Night"))
sum



fit <- glmmPQL(toGUDScale ~ Treatment + Night, 
               random = ~ 1|ID, family = "poisson", data=data)
summary(fit)

tab<-aggregate(toGUD ~ Treatment + Night ,data=data, FUN = "summary")
tab
sum = summarySE(data, 
                measurevar="toGUD", 
                groupvars=c("Treatment", "Night"))
sum

fit <- glmmPQL(timeinGUD ~ Treatment + Night, 
               random = ~ 1|ID, family = "gaussian", data=data)
summary(fit)

tab<-aggregate(timeinGUD ~ Treatment + Night ,data=data, FUN = "summary")
tab
sum = summarySE(data, 
                measurevar="timeinGUD", 
                groupvars=c("Treatment", "Night"))
sum


ggplot(data, aes(x=Treatment, y=toGUD, fill=Treatment)) +
  geom_boxplot() +
  xlab("Treatment Type")+
  ylab("Difference in latency to start foraging (seconds)")+
  scale_fill_manual(values = c("orange", "green"),labels=c("A" = "Anthropogenic Noise", "F" = "Natural Noise")) +
  ylim(-100,250)+   ### sets the y axis limit
  geom_jitter(position = position_jitter(0.1)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    #panel.border = element_blank(),
    axis.line = element_line(color = "black"),
    axis.line.x = element_line(color="black", size = 1),
    axis.line.y = element_line(color="black", size = 1),
    axis.title.x = element_text(colour="grey20",size=16,face="bold"),
    axis.title.y = element_text(colour="grey20",size=16,face="bold"), 
    axis.text.x = element_text(colour="grey20",size=12,face="bold"),
    axis.text.y = element_text(colour="grey20",size=12,face="bold")
  )

ggplot(data, aes(x=Treatment, y=timeinGUD, fill=Treatment)) +
  geom_boxplot() +
  xlab("Treatment Type")+
  ylab("Difference in foraging time (seconds)")+
  scale_fill_manual(values = c("orange", "green"),labels=c("A" = "Anthropogenic Noise", "F" = "Natural Noise")) +
  ylim(-400,200)+   ### sets the y axis limit
  geom_jitter(position = position_jitter(0.1)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    #panel.border = element_blank(),
    axis.line = element_line(color = "black"),
    axis.line.x = element_line(color="black", size = 1),
    axis.line.y = element_line(color="black", size = 1),
    axis.title.x = element_text(colour="grey20",size=16,face="bold"),
    axis.title.y = element_text(colour="grey20",size=16,face="bold"), 
    axis.text.x = element_text(colour="grey20",size=12,face="bold"),
    axis.text.y = element_text(colour="grey20",size=12,face="bold")
  )
#####Vocalization Production
data <- read.csv("Data_calls.csv")

### data not normally distributed use poisson distribution
shapiro.test(data$TotalCalls)

data$Night <- as.character(data$Night)
data$Night <- factor(data$Night, levels=unique(data$Night))

data$Treatment <- as.character(data$Treatment)
data$Treatment <- factor(data$Treatment, levels=unique(data$Treatment))

fit <- glmmPQL(TotalCalls ~ Treatment + Night, 
               random = ~ 1|ID, family = "poisson", data=data)
summary(fit)

tab<-aggregate(Calls ~ Treatment + Night ,data=data, FUN = "summary")
tab
sum = summarySE(data, 
                measurevar="TotalCalls", 
                groupvars=c("Treatment", "Night"))
sum