setwd("~/Desktop/R/RCallNumberSpecies")

###Power Analysis FOR BULBUL


RCallNumberBulbul = read.csv("RCallNumberBulbul.csv")
RCallNumberBulbul$Calls = as.numeric(RCallNumberBulbul$Calls)
RCallNumberBulbul$Proportion_Calls = as.numeric((RCallNumberBulbul$Proportion_Calls))

pf = sum(RCallNumberBulbul[1:6,4])/sum(RCallNumberBulbul$Calls)
pm = sum(RCallNumberBulbul[7:11,4])/sum(RCallNumberBulbul$Calls)
pa = sum(RCallNumberBulbul[12:14,4])/sum(RCallNumberBulbul$Calls)

muf = ((sum(RCallNumberBulbul[1:6,4])/6) - (sum(RCallNumberBulbul$Calls)/24))^2

mum = ((sum(RCallNumberBulbul[7:14,4])/8) - (sum(RCallNumberBulbul$Calls)/24))^2

mua = ((sum(RCallNumberBulbul[15:24,4])/10) - (sum(RCallNumberBulbul$Calls)/24))^2

test = lm(Calls~Habitat,RCallNumberBulbul)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##Bulbul histogram

library(ggplot2)
dodge = position_dodge(1)
RCallNumberBulbulHistogram = ggplot(data = RCallNumberBulbul, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2) +
  scale_fill_manual("Habitat", values = c("orange", "green", "blue"))
RCallNumberBulbulHistogram

## Summarizes data.
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)}
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm))}, measurevar)
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)}

##then create a new data frame 
RCallNumberBulbulsumm = summarySE(RCallNumberBulbul, measurevar = "Calls", groupvars=c("Species", "Habitat"))


hp4=ggplot(data=RCallNumberBulbulsumm, aes(x=Habitat, y=Calls, ymin=Calls-se, ymax=Calls+se, group = "Species", fill=factor(Habitat))) +
  geom_bar(stat="identity", position = "dodge")+
  geom_errorbar(width=0.1, position="dodge") +
  xlab("Red Bulbul") +
  theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())
  

hp4


test = aov(Calls~Habitat, RCallNumberBulbul)
summary(test)
TukeyHSD(test)


###FOR ZEBRA DOVE

RCallNumberZebraDove = read.csv("RCallNumberZebraDove.csv")
RCallNumberZebraDove$Calls = as.numeric(RCallNumberZebraDove$Calls)
RCallNumberZebraDove$Proportion_Calls = as.numeric((RCallNumberZebraDove$Proportion_Calls))

pf = sum(RCallNumberZebraDove[1:6,4])/sum(RCallNumberZebraDove$Calls)
pm = sum(RCallNumberZebraDove[7:14,4])/sum(RCallNumberZebraDove$Calls)
pa = sum(RCallNumberZebraDove[15:24,4])/sum(RCallNumberZebraDove$Calls)

muf = ((sum(RCallNumberZebraDove[1:6,4])/6) - (sum(RCallNumberZebraDove$Calls)/24))^2

mum = ((sum(RCallNumberZebraDove[7:14,4])/8) - (sum(RCallNumberZebraDove$Calls)/24))^2

mua = ((sum(RCallNumberZebraDove[15:24,4])/10) - (sum(RCallNumberZebraDove$Calls)/24))^2

test = lm(Calls~Habitat,RCallNumberZebraDove)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##Histogram for zebra dove

library(ggplot2)
dodge = position_dodge(1)
RCallNumberZebraDoveHistogram = ggplot(data = RCallNumberZebraDove, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2) 
RCallNumberZebraDoveHistogram

test = aov(Calls~Habitat, RCallNumberZebraDove)
summary(test)
TukeyHSD(test)


###FOR MYNA

RCallNumberMyna = read.csv("RCallNumberMyna.csv")
RCallNumberMyna$Calls = as.numeric(RCallNumberMyna$Calls)
RCallNumberMyna$Proportion_Calls = as.numeric(RCallNumberMyna$Proportion_Calls)

pf = sum(RCallNumberMyna[1:6,4])/sum(RCallNumberMyna$Calls)
pm = sum(RCallNumberMyna[7:14,4])/sum(RCallNumberMyna$Calls)
pa = sum(RCallNumberMyna[15:24,4])/sum(RCallNumberMyna$Calls)

muf = ((sum(RCallNumberMyna[1:6,4])/6) - (sum(RCallNumberMyna$Calls)/24))^2

mum = ((sum(RCallNumberMyna[7:14,4])/8) - (sum(RCallNumberMyna$Calls)/24))^2

mua = ((sum(RCallNumberMyna[15:24,4])/10) - (sum(RCallNumberMyna$Calls)/24))^2

test = lm(Calls~Habitat,RCallNumberMyna)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##HISTOGRAM FOR MYNA

library(ggplot2)
dodge = position_dodge(1)
RCallNumberMynaHistogram = ggplot(data = RCallNumberMyna, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2)
RCallNumberMynaHistogram

test = aov(Calls~Habitat, RCallNumberMyna)
summary(test)
TukeyHSD(test)


###FOR RED JUNGLE FOWL


RCallNumberRedJungleFowl = read.csv("RCallNumberRedJungleFowl.csv")
RCallNumberRedJungleFowl$Calls = as.numeric(RCallNumberRedJungleFowl$Calls)
RCallNumberRedJungleFowl$Proportion_Calls = as.numeric(RCallNumberRedJungleFowl$Proportion_Calls)

pf = sum(RCallNumberRedJungleFowl[1:6,4])/sum(RCallNumberRedJungleFowl$Calls)
pm = sum(RCallNumberRedJungleFowl[7:14,4])/sum(RCallNumberRedJungleFowl$Calls)
pa = sum(RCallNumberRedJungleFowl[15:24,4])/sum(RCallNumberRedJungleFowl$Calls)

muf = ((sum(RCallNumberRedJungleFowl[1:6,4])/6) - (sum(RCallNumberRedJungleFowl$Calls)/24))^2

mum = ((sum(RCallNumberRedJungleFowl[7:14,4])/8) - (sum(RCallNumberRedJungleFowl$Calls)/24))^2

mua = ((sum(RCallNumberRedJungleFowl[15:24,4])/10) - (sum(RCallNumberRedJungleFowl$Calls)/24))^2

test = lm(Calls~Habitat,RCallNumberRedJungleFowl)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##Histogram for Jungle Fowl

library(ggplot2)
dodge = position_dodge(1)
RCallNumberRedJungleFowlHistogram = ggplot(data = RCallNumberRedJungleFowl, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2)
RCallNumberRedJungleFowlHistogram

test = aov(Calls~Habitat, RCallNumberRedJungleFowl)
summary(test)
TukeyHSD(test)

###FOR SILVEREYE

RCallNumberSilvereye = read.csv("RCallNumberSilvereye.csv")
RCallNumberSilvereye$Calls = as.numeric(RCallNumberSilvereye$Calls)
RCallNumberSilvereye$Proportion_Calls = as.numeric(RCallNumberSilvereye$Proportion_Calls)

pf = sum(RCallNumberSilvereye[1:6,4])/sum(RCallNumberSilvereye$Calls)
pm = sum(RCallNumberSilvereye[7:14,4])/sum(RCallNumberSilvereye$Calls)
pa = sum(RCallNumberSilvereye[15:24,4])/sum(RCallNumberSilvereye$Calls)

muf = ((sum(RCallNumberSilvereye[1:6,4])/6) - (sum(RCallNumberSilvereye$Calls)/24))^2

mum = ((sum(RCallNumberSilvereye[7:14,4])/8) - (sum(RCallNumberSilvereye$Calls)/24))^2

mua = ((sum(RCallNumberSilvereye[15:24,4])/10) - (sum(RCallNumberSilvereye$Calls)/24))^2

test = lm(Calls~Habitat,RCallNumberSilvereye)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##Histogram for Silvereye

library(ggplot2)
dodge = position_dodge(1)
RCallNumberSilvereyeHistogram = ggplot(data = RCallNumberSilvereye, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2)
RCallNumberSilvereyeHistogram

test = aov(Calls~Habitat, RCallNumberSilvereye)
summary(test)
TukeyHSD(test)


###FOR KINGFISHER

RCallNumberKingfisher = read.csv("RCallNumberKingfisher.csv")
RCallNumberKingfisher$Calls = as.numeric(RCallNumberKingfisher$Calls)
RCallNumberKingfisher$Proportion_Calls = as.numeric(RCallNumberKingfisher$Proportion_Calls)

pf = sum(RCallNumberKingfisher[1:4,4])/sum(RCallNumberKingfisher$Calls)
pm = sum(RCallNumberKingfisher[5:7,4])/sum(RCallNumberKingfisher$Calls)
pa = sum(RCallNumberKingfisher[8:10,4])/sum(RCallNumberKingfisher$Calls)

muf = ((sum(RCallNumberKingfisher[1:4,4])/4) - (sum(RCallNumberKingfisher$Calls)/10))^2

mum = ((sum(RCallNumberKingfisher[5:7,4])/3) - (sum(RCallNumberKingfisher$Calls)/10))^2

mua = ((sum(RCallNumberKingfisher[8:10,4])/3) - (sum(RCallNumberKingfisher$Calls)/10))^2

test = lm(Calls~Habitat,RCallNumberKingfisher)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##Histogram for Kingfisher

library(ggplot2)
dodge = position_dodge(1)
RCallNumberKingfisherHistogram = ggplot(data = RCallNumberKingfisher, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2)
RCallNumberKingfisherHistogram

test = aov(Calls~Habitat, RCallNumberKingfisher)
summary(test)
TukeyHSD(test)


###FOR WAXBILL

RCallNumberWaxbill = read.csv("RCallNumberWaxbill.csv")
RCallNumberWaxbill$Calls = as.numeric(RCallNumberWaxbill$Calls)
RCallNumberWaxbill$Proportion_Calls = as.numeric(RCallNumberWaxbill$Proportion_Calls)

pf = sum(RCallNumberWaxbill[1:4,4])/sum(RCallNumberWaxbill$Calls)
pm = sum(RCallNumberWaxbill[5:7,4])/sum(RCallNumberWaxbill$Calls)
pa = sum(RCallNumberWaxbill[8:10,4])/sum(RCallNumberWaxbill$Calls)

muf = ((sum(RCallNumberWaxbill[1:4,4])/4) - (sum(RCallNumberWaxbill$Calls)/10))^2

mum = ((sum(RCallNumberWaxbill[5:7,4])/3) - (sum(RCallNumberWaxbill$Calls)/10))^2

mua = ((sum(RCallNumberWaxbill[8:10,4])/3) - (sum(RCallNumberWaxbill$Calls)/10))^2

test = lm(Calls~Habitat,RCallNumberWaxbill)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##Histogram for Waxbill

library(ggplot2)
dodge = position_dodge(1)
RCallNumberWaxbillHistogram = ggplot(data = RCallNumberWaxbill, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2)
RCallNumberWaxbillHistogram

test = aov(Calls~Habitat, RCallNumberWaxbill)
summary(test)
TukeyHSD(test)

###FOR FIRETAIL

RCallNumberFiretail = read.csv("RCallNumberFiretail.csv")
RCallNumberFiretail$Calls = as.numeric(RCallNumberFiretail$Calls)
RCallNumberFiretail$Proportion_Calls = as.numeric(RCallNumberFiretail$Proportion_Calls)

pf = sum(RCallNumberFiretail[1:4,4])/sum(RCallNumberFiretail$Calls)
pm = sum(RCallNumberFiretail[5:7,4])/sum(RCallNumberFiretail$Calls)
pa = sum(RCallNumberFiretail[8:10,4])/sum(RCallNumberFiretail$Calls)

muf = ((sum(RCallNumberFiretail[1:4,4])/4) - (sum(RCallNumberFiretail$Calls)/10))^2

mum = ((sum(RCallNumberFiretail[5:7,4])/3) - (sum(RCallNumberFiretail$Calls)/10))^2

mua = ((sum(RCallNumberFiretail[8:10,4])/3) - (sum(RCallNumberFiretail$Calls)/10))^2

test = lm(Calls~Habitat,RCallNumberFiretail)

sig = anova(test)["Residuals", "Mean Sq"]

f = sqrt(((pf*muf) + (pm*mum) + (pa*mua))/sig)

library(pwr)
pwr.anova.test(k =3, f = f, sig.level = 0.05, power = 0.8) 

##Histogram for Firetail

dodge = position_dodge(1)
RCallNumberFiretailHistogram = ggplot(data = RCallNumberFiretail, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2) +
  scale_fill_manual("Habitat", values = c("orange", "green", "blue"))
RCallNumberFiretailHistogram



test = aov(Calls~Habitat, RCallNumberFiretail)
summary(test)
TukeyHSD(test)


 ###For Fruit Dove
RCallNumberFruitDove = read.csv("RCallNumberFruitDove.csv")
RCallNumberFruitDove$Calls = as.numeric(RCallNumberFruitDove$Calls)
RCallNumberFruitDove$Proportion_Calls = as.numeric(RCallNumberFruitDove$Proportion_Calls)

##Fruit Dove Histogram 

dodge = position_dodge(1)
RCallNumberFruitDoveHistogram = ggplot(data = RCallNumberFruitDove, aes (x= Species, y= Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2) +
  scale_fill_manual("Habitat", values = c("orange", "green", "blue"))
RCallNumberFruitDoveHistogram

test = aov(Calls~Habitat, RCallNumberFruitDove)
summary(test)
TukeyHSD(test)
##COMBINED SPECIES FOR HISTOGRAM
setwd("~/Desktop/R/RCallNumberSpecies")

RCallNumberCombined <- rbind(RCallNumberBulbul, RCallNumberRedJungleFowl, RCallNumberZebraDove, RCallNumberMyna, RCallNumberSilvereye, RCallNumberWaxbill, RCallNumberFiretail, RCallNumberKingfisher, RCallNumberFruitDove)

##ANOVA TEST

test = aov(Calls~Habitat*Species, RCallNumberCombined)
summary(test)
TukeyHSD(test)

##HISTOGRAM


library(ggplot2)

## Summarizes data.
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)}
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm))}, measurevar)
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)}

##then create a new data frame 
RCallNumberCombinedsumm = summarySE(RCallNumberCombined, measurevar = "Calls", groupvars=c("Species", "Habitat"))


RCallNumberCombinedHistogram=ggplot(data=RCallNumberCombinedsumm, aes(Habitat, Calls, ymin=Calls-se, ymax=Calls+se, group = Species, fill=Habitat)) +
  geom_bar(stat="identity", position = dodge)+
  scale_fill_manual("Habitat", values = c("darkorange", "green4", "blue2")) +
  facet_wrap(~Species, scales="free", drop = TRUE) +
  geom_errorbar(width=0.1, position=dodge) +
  theme_bw() +
  theme(text= element_text(size= 20),
        axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.x=element_blank())+
  ggtitle("Average Number of Calls per Survey Replicate")
RCallNumberCombinedHistogram


##Correlations

##ElevationCorrelation
setwd("~/Desktop/R/RBirdCorrelation")

RElevationCorrelation = read.csv("RElevationCorrelation.csv")
library(ggplot2)

plot1 <- ggplot(RElevationCorrelation, aes(x=Elevation, y=Species_Richness)) +
  geom_point() +
  theme(text = element_text(size=20)) +
  geom_smooth(method="lm") + 
plot1


test = lm(Species_Richness~Elevation, RElevationCorrelation)
summary(test)


##CanopyCoverCorrelation

setwd("~/Desktop/R/RBirdCorrelation")

RCanopyCoverCorrelation = read.csv("RCanopyCoverCorrelation.csv")
library(ggplot2)

plot2 <- ggplot(RCanopyCoverCorrelation, aes(x=Canopy_Cover, y=Species_Richness)) +
  geom_point() +
  geom_smooth(method="lm") +
  theme(text = element_text(size=20))
plot2

test = lm(Species_Richness~Canopy_Cover, RCanopyCoverCorrelation)
summary(test)

##Invasive vs. Native

setwd("~/Desktop/R")
InvasiveNative = read.csv("InvasiveNative.csv")
library(ggplot2)
dodge=position_dodge(.7)

plot3 <- ggplot(data = InvasiveNative, aes(x = Species, y = Calls, group = Habitat, fill = Habitat)) +
  geom_bar(stat="identity")
plot3

InvasiveNative$Species = as.factor(InvasiveNative$Species)
test = aov(Calls~Species, InvasiveNative)
summary(test)
TukeyHSD(test)

##Grazing vs. Crop

setwd("~/Desktop/R")

GrazingCrop = read.csv("GrazingCrop.csv")
GrazingCrop$Calls = as.numeric(GrazingCrop$Calls)

GrazingCalls = subset(GrazingCrop, regexpr("Grazing",ignore.case = TRUE, GrazingCrop$Habitat) >0)
CropCalls = subset(GrazingCrop, regexpr("Crop", ignore.case = TRUE, GrazingCrop$Habitat)>0)

##FINDING P-VALUES FOR GRAZING VS. CROP CALLS
### FOR SHAPIRO.........GREATER THAN .05 USE T.TEST, 
### LESS THAN .05 USE WILCOX

##Subset Bulbul 
shapiro.test(GrazingCallsBulbul$Calls)
shapiro.test(CropCallsBulbul$Calls)

GrazingCallsBulbul = subset(GrazingCalls, regexpr("Red-vented Bulbul", ignore.case = TRUE, GrazingCalls$Species) >0)
CropCallsBulbul = subset (CropCalls, regexpr("Red-vented Bulbul", ignore.case = TRUE, CropCalls$Species) >0)

t.test(GrazingCallsBulbul$Calls, CropCallsBulbul$Calls, data=GrazingCrop)

##Subset Zebra Dove
shapiro.test(GrazingCallsZebraDove$Calls)
shapiro.test(CropCallsZebraDove$Calls)

GrazingCallsZebraDove= subset(GrazingCalls, regexpr("Zebra Dove", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsZebraDove = subset(CropCalls, regexpr("Zebra Dove", ignore.case = TRUE, CropCalls$Species) > 0)

t.test(GrazingCallsZebraDove$Calls, CropCallsZebraDove$Calls, data = GrazingCrop)

##Subset FruitDove
shapiro.test(GrazingCallsFruitDove$Calls)
shapiro.test(CropCallsFruitDove$Calls)

GrazingCallsFruitDove = subset(GrazingCalls, regexpr("Grey-green Fruit Dove", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsFruitDove = subset(CropCalls, regexpr("Grey-green Fruit Dove", ignore.case = TRUE, CropCalls$Species) > 0)

wilcox.test(GrazingCallsFruitDove$Calls, CropCallsFruitDove$Calls, data = GrazingCrop)

##Subset Jungle Fowl
shapiro.test(GrazingCallsRedJungleFowl$Calls)
shapiro.test(CropCallsRedJungleFowl$Calls)

GrazingCallsRedJungleFowl = subset(GrazingCalls, regexpr("Red Jungle Fowl", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsRedJungleFowl = subset(CropCalls, regexpr("Red Jungle Fowl", ignore.case = TRUE, CropCalls$Species) > 0)

t.test(GrazingCallsRedJungleFowl$Calls, CropCallsRedJungleFowl$Calls, data = GrazingCrop)
##Subset Kingfisher 
shapiro.test(GrazingCallsKingfisher$Calls)
shapiro.test(CropCallsKingfisher$Calls)

GrazingCallsKingfisher = subset(GrazingCalls, regexpr("Mo'orean Kingfisher", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsKingfisher = subset(CropCalls, regexpr("Mo'orean Kingfisher", ignore.case = TRUE, CropCalls$Species) > 0)

t.test(GrazingCallsKingfisher$Calls, CropCallsKingfisher$Calls, data = GrazingCrop)

##Subset for Common Waxbill
shapiro.test(GrazingCallsWaxbill$Calls)
shapiro.test(CropCallsKingfisher$Calls)

GrazingCallsWaxbill = subset(GrazingCalls, regexpr("Common Waxbill", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsWaxbill = subset(CropCalls, regexpr("Common Waxbill", ignore.case = TRUE, CropCalls$Species) > 0)

wilcox.test(GrazingCallsWaxbill$Calls, CropCallsWaxbill$Calls, data = GrazingCrop)

##Subset for Silvereye
shapiro.test(GrazingCallsSilvereye$Calls)
shapiro.test(CropCallsSilvereye$Calls)

GrazingCallsSilvereye = subset(GrazingCalls, regexpr("Silvereye", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsSilvereye = subset(CropCalls, regexpr("Silvereye", ignore.case = TRUE, CropCalls$Species) > 0)

t.test(GrazingCallsSilvereye$Calls, CropCallsSilvereye$Calls, data = GrazingCrop)

##Subset for Firetail
shapiro.test(GrazingCallsFiretail$Calls)
shapiro.test(CropCallsFiretail$Calls)

GrazingCallsFiretail = subset(GrazingCalls, regexpr("Firetail", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsFiretail = subset(CropCalls, regexpr("Firetail", ignore.case = TRUE, CropCalls$Species) > 0)

t.test(GrazingCallsFiretail$Calls, CropCallsFiretail$Calls, data = GrazingCrop)

##Subset for Myna
shapiro.test(GrazingCallsMyna$Calls)
shapiro.test(CropCallsMyna$Calls)

GrazingCallsMyna = subset(GrazingCalls, regexpr("Common Myna", ignore.case = TRUE, GrazingCalls$Species) > 0)
CropCallsMyna = subset(CropCalls, regexpr("Common Myna", ignore.case = TRUE, CropCalls$Species) > 0)

t.test(GrazingCallsMyna$Calls, CropCallsWaxbill$Calls, data = GrazingCrop)

## Summarizes data.
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)}
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm))}, measurevar)
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)}

##then create a new data frame 
GrazingCropsumm = summarySE(GrazingCrop, measurevar = "Calls", groupvars=c("Species", "Habitat"))

###IF this isn't work.. rerun "GrazingCrop" from above (wayy above..)

GrazingCropHistogram=ggplot(data=GrazingCropsumm, aes(Habitat, Calls, ymin=Calls-se, ymax=Calls+se, group = Species, fill=Habitat)) +
  geom_bar(stat="identity", position = dodge)+
  scale_fill_manual("Habitat", values = c("yellow", "orange")) +
  facet_wrap(~Species, scales="free", drop = TRUE) +
  geom_errorbar(width=0.1, position=dodge) +
  theme_bw() +
  theme(text= element_text(size= 20),
        axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.x=element_blank())+
  ggtitle("Average Number of Calls in Different Managed Land")
GrazingCropHistogram



##Don't use
library(ggplot2)
dodge = position_dodge(1)
GrazingCrop = ggplot(data = GrazingCrop, aes(x= Species, y= Calls, group = Habitat, fill=Habitat)) +
  geom_bar(stat="summary", fun.y="mean", position=dodge) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=dodge, width = .2) +
  scale_fill_manual("Habitat", values = c("yellow", "orange")) +
  theme(text = element_text(size=20)) +
  facet_wrap(~Species, scales="free") +
  theme_bw() +
  theme(text = element_text(size=20), 
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  ggtitle("Average Number of Calls in Grazing and Crop Lands")
GrazingCrop