##############
###packages###
##############
library(moments)
library(ggplot2)
library(dplyr)
library(gridExtra)
library(xlsx)

#################
###Import data###
#################
all.bins <- read.csv("newdata.csv")
#log-transform mass
all.bins$log10.mass <- log10(all.bins$net.mass)
all.bins$log.mass <- log(all.bins$net.mass)


##normality tests
agostino.test(all.bins$log.mass)
anscombe.test(all.bins$log.mass)

##IQR summary
#create summary at genus level
genus_df <- subset(all.bins, all.bins$Genus != "")
genus_df$log_mass <- log10(genus_df$net.mass)
genus.IQR <- aggregate(genus_df$log_mass, list(genus_df$Genus), IQR)
genus.size <- aggregate(genus_df$log_mass, list(genus_df$Genus), length)
names(genus.IQR) <- c("Genus", "IQR")
names(genus.size) <- c("Genus", "BINs.weighed")
genus_summary <- merge(genus.size, genus.IQR,by = "Genus")
genus_summary <- subset(genus_summary, IQR != 0)
#median IQR
median(genus_summary$IQR) #=0.163


#ANOVAS
##One way
###family level
summary(aov(log.mass ~ Family, data=all.bins))
###family level with single data groups dropped
family.count <- aggregate(all.bins$log.mass, list(all.bins$Family), length)
family.count <- subset(family.count, family.count$x > 1)
all.bins.family <- subset(all.bins, all.bins$Family %in% family.count$Group.1)
summary(aov(log.mass ~ Family, data=all.bins.family))
###subfamily level
all.bins.subfamily = subset(all.bins, all.bins$Subfamily != "") #subset data
summary(aov(log.mass ~ Subfamily, data=all.bins.subfamily))
###subfamily level with single data groups dropped
subfamily.count <- aggregate(all.bins.subfamily$log.mass, list(all.bins.subfamily$Subfamily), length)
subfamily.count <- subset(subfamily.count, subfamily.count$x > 1)
all.bins.subfamily <- subset(all.bins.subfamily, all.bins.subfamily$Subfamily %in% subfamily.count$Group.1)
summary(aov(log.mass ~ Subfamily, data=all.bins.subfamily))
###genus level
all.bins.genus = subset(all.bins, all.bins$Genus != "")
genus.count <- aggregate(all.bins.genus$log.mass, list(all.bins.genus$Genus), length)
genus.count <- subset(genus.count, genus.count$x >= 2)
all.bins.genus <- subset(all.bins.genus, all.bins.genus$Genus %in% genus.count$Group.1)
summary(aov(log.mass ~ Genus, data=all.bins.genus))
###genus level, groups with more than 10 BINs
genus.count <- subset(genus.count, genus.count$x >= 10)
all.bins.genus <- subset(all.bins.genus, all.bins.genus$Genus %in% genus.count$Group.1)
summary(aov(log.mass ~ Genus, data=all.bins.genus))

###nested aov
#filter for genera with n>=2
all.bins.nested = subset(all.bins, all.bins$Genus != "")
genus.count <- aggregate(all.bins.nested$log.mass, list(all.bins.nested$Genus), length)
genus.count <- subset(genus.count, genus.count$x >= 2)
all.bins.nested <- subset(all.bins.nested, all.bins.nested$Genus %in% genus.count$Group.1)
#summarize by family counting genera
nested <- all.bins.nested %>% group_by(Family,Genus) %>%
  summarize(Count = n())
#filter for families with more than 1 genus
nested.count <- aggregate(nested$Genus, list(nested$Family), length)
nested.count <- subset(nested.count, nested.count$x > 1)
all.bins.nested <- subset(all.bins.nested, all.bins.nested$Family %in% nested.count$Group.1)
#run aov and display summary
log.fg.aov <- aov(log.mass ~ Family/Genus, data=all.bins.nested)
summary(log.fg.aov)

###families with n >= 100
family.count <- aggregate(all.bins$log.mass, list(all.bins$Family), length)
family.count <- subset(family.count, family.count$x >= 100)
large.nested <- subset(all.bins.nested, all.bins.nested$Family %in% family.count$Group.1)
#run aov and display summary
log.fg.aov100 <- aov(log.mass ~ Family/Genus, data=large.nested)
summary(log.fg.aov100)

###paired sample t-test for measurements repeated in 2014 and 2018
data <- read.xlsx("D:/data/BeetleBodyMass/comparison_2014-2018.xlsx", sheetName="Sheet1")
t.test(Weight.2018, Weight.2014, alternative='greater', conf.level=.95, paired=TRUE))

###linear regression of fresh vs. dry weights for the set of 73 specimens
FreshVsDry <- read.xlsx("D:/data/BeetleBodyMass/fresh_vs_dry.xlsx", sheetName="Sheet1")
lm(formula = Dried.Weight ~ Live.Weight, data = FreshVsDry)