################################
## modeling Fossile et al ######

## load files for fish occurrence in sambaqui records and
## fish traits

data <- read.csv("sambaqui_catch.csv", header=T, sep=";", dec=",")
head(data)

## check dataset dimension

dim(data)


###########
## adjust  body size categories (from Quimbayo et al. 2021, for reef fish)
#	s1:  0 - 7 cm
#	s2:  7.1 - 15 cm
#	s3:  15.1 - 30 cm
#	s4:  30.1 - 50 cm
#	s5:  50.1 - 80 cm
#	s6:  > 80 cm


data$body_class <- data$Max_length
body_class <- as.numeric(data$body_class)
body_class <- ifelse(data$Max_length<=7,"s1", "s2")
body_class <- ifelse(data$Max_length>7.1 & data$Max_length<=15.0,"s2", body_class)
body_class <- ifelse(data$Max_length>15.1 & data$Max_length<=30.0,"s3", body_class)
body_class <- ifelse(data$Max_length>=30.1 & data$Max_length<=50.0,"s4", body_class)
body_class <- ifelse(data$Max_length>=50.1 & data$Max_length<=80.0,"s5", body_class)
body_class <- ifelse(data$Max_length>80,"s6", body_class)

length(body_class)

############
# paste column with body class in dataframe

data$body_class <- body_class
data$FG <- paste(data$Trophic_group,data$body_class)

###########
# check for correlation between body mass and maximum length

cor.test(data$Max_mass,data$Max_length)

# correlation is of 0.67, p-value < 2.2e-16 

###########

data1 <- data
data1$Non_ceramic <- ifelse(data1$Non_ceramic==1, data1$Trophic_level, data1$Non_ceramic)
data1$Ceramic <- ifelse(data1$Ceramic==1, data1$Trophic_level, data1$Ceramic)
data1$Modern <- ifelse(data1$Modern==1, data1$Trophic_level, data1$Modern)

head(data1)

###########

data2 <- data
data2$Non_ceramic <- ifelse(data2$Non_ceramic==1, data2$Max_length, data2$Non_ceramic)
data2$Ceramic <- ifelse(data2$Ceramic==1, data2$Max_length, data2$Ceramic)
data2$Modern <- ifelse(data2$Modern==1, data2$Max_length, data2$Modern)

head(data2)

###########
## adjust data for distribution of Body Mass
## across periods

data3 <- data
data3$Non_ceramic <- ifelse(data3$Non_ceramic==1, data3$Max_mass, data3$Non_ceramic)
data3$Ceramic <- ifelse(data3$Ceramic==1, data3$Max_mass, data3$Ceramic)
data3$Modern <- ifelse(data3$Modern==1, data3$Max_mass, data3$Modern)


###########
#Null model analysis

# load data 

# regional availability of species 

community_data <- read.csv ("sambaqui_data.csv",
                            h=T, sep=";")

# catch data, per period
catch_data <- read.csv ("sambaqui_catch.csv",
                        h=T,sep=";")
catch_data$Species <- gsub (" ","_", catch_data$Species)

# 1) match regional e local
# 24 species from the catch dataset are not in the regional community_data
(catch_data$Species [which(catch_data$Species %in% community_data$Species == F)])

# imput the pool
community_data<-rbind (community_data [,c("Species",
                                          "Trophic_level",
                                          "Max_length",
                                          "Max_mass", 
                                          "Trophic_group",
                                          "Non_ceramic",
                                          "Ceramic",
                                          "Modern")],# select cols in both datasets
                       catch_data [which(catch_data$Species %in% community_data$Species == F),])

# 2) sampling
# number of species in each period
spp_periods <- colSums (catch_data[,c("Non_ceramic",
                                      "Ceramic",
                                      "Modern")])
# set the minimum
min_spp_periods <- min (spp_periods)

# adjusting data ("," into ".", categorical to quantitative)
community_data$Max_length<-as.numeric(gsub (",",".",community_data$Max_length)) # max length
community_data$Max_mass<-as.numeric(gsub (",",".",community_data$Max_mass)) # max mass
community_data$Trophic_level<-as.numeric(gsub (",",".",community_data$Trophic_level)) # trophic level

# trophic group is not used

community_data<-community_data[,-which(colnames(community_data) == "Trophic_group")]

# ===================================================
# null model considering a random expectation for each period, separately
# ie. different number of species across periods
nsamples <- 1000
rdm_set <- lapply (spp_periods, function (i) 
  
  
  replicate (nsamples, # random composition
             sample (community_data$Species,
                     size = i, 
                     replace =F)))
#  each period
# and each random composition
# nested loop
community_data_rdm <- lapply (rdm_set, function (i) # for each period
  
  lapply (seq (1,ncol(i)), function (k) # for each one of the nsamples
    
    # get trait data
    community_data[which(community_data$Species %in% i[,k]),c("Trophic_level",
                                                              "Max_mass",
                                                              "Max_length")]
    
  )
)

# average values across random samples
# for each period
averaged_values<- lapply (community_data_rdm, function (i) # for each period
  
  Reduce("+",i)/length(i)) # get the average


# create a dataframe for analysis
data_analysis <- lapply (colnames(averaged_values[[1]]), function (i){
  
  data_analysis <- data.frame (
    
    rbind (
      
      # NUll  
      cbind ("values" = averaged_values [[1]][,i],
             "data" = "Non-ceramic",
             "null" = "yes"),
      cbind ("values" = averaged_values [[2]][,i],
             "data" = "Ceramic",
             "null" = "yes"),
      cbind ("values" = averaged_values [[3]][,i],
             "data" = "Modern",
             "null" = "yes"),
      # non-ceramic period
      cbind ("values" = community_data[which(community_data$Non_ceramic == 1),i],
             "data" = "Non-ceramic",
             "null" = "no"),
      # ceramic period
      cbind ("values" = community_data[which(community_data$Ceramic == 1),i],
             "data" = "Ceramic",
             "null" = "no"),
      # modern communities
      cbind ("values" = community_data[which(community_data$Modern == 1),i],
             "data" = "Modern",
             "null" = "no"))
  )
  # identify the data
  data_analysis <- cbind (data_analysis,
                          "Trait" = i)
  
  ; # return
  data_analysis
})

# melt to have a dataframe including all data
data_analysis<-do.call(rbind, data_analysis)
data_analysis$values<-as.numeric(data_analysis$values)
data_analysis$data<-factor(data_analysis$data,
                           levels = c("Non-ceramic",
                                      "Ceramic",
                                      "Modern"))
# plot
require(ggplot2)

ggplot(data_analysis[which(data_analysis$null == "no"),],  # subset (null dataset)
       aes (x = data, 
            y = log(values),
            fill=data)) +
  
  # facets
  facet_wrap(~Trait,ncol=3,scales = "free_y",
             labeller = as_labeller(c("Max_length" = "Maximum length (cm)", 
                                      "Max_mass" = "Maximum mass (g)",  
                                      "Trophic_level" = "Trophic level"))) +
  
  # violin plot for the null dataset
  geom_violin(size=1,
              col = "gray80",
              fill = "gray20",
              alpha = 0.5) +
  
  # boxplot for the observed dataset
  geom_boxplot(data = data_analysis[which(data_analysis$null == "yes"),],  # subset (observed dataset)
               aes (x = data, 
                    y = log(values), width = .3, alpha=0.9) +
                 
                 scale_fill_manual(values=c("#005F73","#EE9B00","#CA6702")) + 
                 
                 geom_point(data = data_analysis[which(data_analysis$null == "yes"),], shape = 21, position = position_jitterdodge(dodge.width = 0.3), alpha=0.5) +
                 
                 scale_fill_manual(values = c("#005F73","#EE9B00","#CA6702")) + 
                 
                 
                 #stat_summary(fun=mean,
                 #             geom="point",
                 #             shape=19, size=3) + 
                 
                 theme_gray() + xlab ("") + 
                 ylab (expression("Trait value")) +
                 theme (axis.title = element_text(size=12),
                        axis.text.x = element_blank(),
                        axis.ticks.x = element_blank(),
                        axis.text.y = element_text(size=10),
                        axis.ticks.y = element_blank(),
                        strip.background = element_blank(),
                        legend.title = element_blank())


###########
## test for significant differences between boxplots using ANOVA

## Apply ANOVA below
model1 = aov(log(box1$values)~box1$ind)
model2 = aov(log(box2$values)~box2$ind)
model3 = aov(log(box3$values)~box3$ind)


## Run post-hoc Tukey Test between groups
tk1 <- TukeyHSD(model1)
tk2 <- TukeyHSD(model2)
tk3 <- TukeyHSD(model3)

###########
## test for significant differences between observed and null models

# trophic level
trof <- subset(data_analysis, data_analysis$Trait=="Trophic_level", drop=T)
trof_nc <- subset(trof, trof$data=="Non-ceramic", drop=T)
trof_c <- subset(trof, trof$data=="Ceramic", drop=T)
trof_m <- subset(trof, trof$data=="Modern", drop=T)

pairwise.t.test(trof_m$values, trof_m$null, p.adjust="holm")
#p-value=0.00052 TROPHIC NC
#p-value=0.018 TROPHIC C
#p-value=0.69 TROPHIC M

# body mass
bmass <- subset(data_analysis, data_analysis$Trait=="Max_mass", drop=T)
mass_nc <- subset(bmass, bmass$data=="Non-ceramic", drop=T)
mass_c <- subset(bmass, bmass$data=="Ceramic", drop=T)
mass_m <- subset(bmass, bmass$data=="Modern", drop=T)

pairwise.t.test(mass_nc$values, mass_nc$null, p.adjust="holm", alternative="greater")
#p-value=0.15 BMASS NC
#p-value=0.14 BMASS C
#p-value=0.61 BMASS M

# body size
bsize <- subset(data_analysis, data_analysis$Trait=="Max_length", drop=T)
size_nc <- subset(bsize, bsize$data=="Non-ceramic", drop=T)
size_c <- subset(bsize, bsize$data=="Ceramic", drop=T)
size_m <- subset(bsize, bsize$data=="Modern", drop=T)

pairwise.t.test(size_m$values, size_m$null, p.adjust="holm")
#p-value=0.00079 SIZE NC
#p-value=0.0013 SIZE C
#p-value=0.52 SIZE M

###########
###########
# Barplot for proportion of trophic groups, body size classes and
# functional groups across studied periods 

fg_ceramic = table(data$FG,data$Ceramic)
fg_nceramic = table(data$FG,data$Non_ceramic)
fg_modern = table(data$FG,data$Modern)

fgs<- cbind(fg_ceramic,fg_nceramic,fg_modern)
colnames(fgs) <- paste(c("0","Ceramic","0","Non_Ceramic","0","Modern"))
fgs <- fgs[,c("Ceramic","Non_Ceramic","Modern")]

fgs_prop <- apply(fgs, 2, function(x){x/sum(x)})

###########
## set color to Functional Groups

color=as.data.frame(rownames(fgs_prop))
colnames(color)=paste("color")

color_fgs <- ifelse(color$color=="HERB s1","#DBE6DB","#BBCFBB")
color_fgs <- ifelse(color$color=="HERB s2","#BBCFBB",color_fgs)
color_fgs <- ifelse(color$color=="HERB s3","#ABC4AB",color_fgs)
color_fgs <- ifelse(color$color=="HERB s4","#94B394",color_fgs)
color_fgs <- ifelse(color$color=="HERB s5","#7CA27C",color_fgs)
color_fgs <- ifelse(color$color=="HERB s6","#668F66",color_fgs)
color_fgs <- ifelse(color$color=="INV s1","#FDD9E0",color_fgs)
color_fgs <- ifelse(color$color=="INV s2","#FBB1C0",color_fgs)
color_fgs <- ifelse(color$color=="INV s3","#FA9EB0",color_fgs)
color_fgs <- ifelse(color$color=="INV s4","#F87791",color_fgs)
color_fgs <- ifelse(color$color=="INV s5","#F65172",color_fgs)
color_fgs <- ifelse(color$color=="INV s6","#F42A52",color_fgs)
color_fgs <- ifelse(color$color=="MCAR s2","#D0B4CE",color_fgs)
color_fgs <- ifelse(color$color=="MCAR s3","#C09BBE",color_fgs)
color_fgs <- ifelse(color$color=="MCAR s4","#B082AD",color_fgs)
color_fgs <- ifelse(color$color=="MCAR s5","#A1689D",color_fgs)
color_fgs <- ifelse(color$color=="MCAR s6","#975E93",color_fgs)
color_fgs <- ifelse(color$color=="OMNI s1","#77B5E4",color_fgs)
color_fgs <- ifelse(color$color=="OMNI s2","#4499DA",color_fgs)
color_fgs <- ifelse(color$color=="OMNI s3","#2885CC",color_fgs)
color_fgs <- ifelse(color$color=="OMNI s4","#226FAA",color_fgs)
color_fgs <- ifelse(color$color=="OMNI s5","#1B5988",color_fgs)
color_fgs <- ifelse(color$color=="OMNI s6","#144366",color_fgs)
color_fgs <- ifelse(color$color=="PISC s3","#CA0902",color_fgs)
color_fgs <- ifelse(color$color=="PISC s4","#B70701",color_fgs)
color_fgs <- ifelse(color$color=="PISC s5","#8D0601",color_fgs)
color_fgs <- ifelse(color$color=="PISC s6","#790501",color_fgs) 
color_fgs <- ifelse(color$color=="PLANK s2","#F4E590",color_fgs)
color_fgs <- ifelse(color$color=="PLANK s3","#F0DC6A",color_fgs)
color_fgs <- ifelse(color$color=="PLANK s4","#EDD445",color_fgs)
color_fgs <- ifelse(color$color=="PLANK s5","#E9CB20",color_fgs)

###########

t_ceramic = table(data$Trophic_group,data$Ceramic)
t_nceramic = table(data$Trophic_group,data$Non_ceramic)
t_modern = table(data$Trophic_group,data$Modern)

tgroups<- cbind(t_ceramic,t_nceramic,t_modern)
colnames(tgroups) <- paste(c("0","Ceramic","0","Non.Ceramic","0","Modern"))
tgroups <- tgroups[,c("Ceramic","Non.Ceramic","Modern")]

tgroups_prop <- apply(tgroups, 2, function(x){x/sum(x)})
colSums(tgroups_prop) ## must equal 1 for each column

###########

b_ceramic = table(data$body_class,data$Ceramic)
b_nceramic = table(data$body_class,data$Non_ceramic)
b_modern = table(data$body_class,data$Modern)

bclass<- cbind(b_ceramic,b_nceramic,b_modern)
colnames(bclass) <- paste(c("0","Ceramic","0","Non_Ceramic","0","Modern"))
bclass <- bclass[,c("Ceramic","Non_Ceramic","Modern")]

bclass_prop <- apply(bclass, 2, function(x){x/sum(x)})
colSums(bclass_prop) ## must equal 1 for each column

#### build barplot for trophic groups only

par(mfrow=c(1,4))
#add body size class barplot
barplot(bclass_prop, col=c("#E8EFFF","#C3CDFF","#95A0D2","#6975A5","#3F4D79","#142850"),
        ylab="Body size class", names=c("NC","C","M"), cex.axis=1.1,
        border=NA, width=0.7, cex.lab=1.1, cex.names=1.1, las=1)
title("A",adj=0, line=1, cex.main=1.7)

#add trophic group plot
barplot(tgroups_prop, col=c("#668F66","#F42A52","#975E93","#144366","#790501","#E9CB20"),
        ylab="Trophic group proportion", names=c("NC","C","M"),
        border=NA, width=0.7, cex.lab=1.1, cex.axis=1.1, cex.names=1.1, las=1)
title("B",adj=0, line=1, cex.main=1.7)


#add FG proportion
barplot(fgs_prop, col=color_fgs, ylab="Functional group proportion", 
        names=c("NC","C","M"),
        border=NA, beside=F, width=0.7, cex.axis=1.1, cex.lab=1.1, cex.names=1.1, las=1)
title("C",adj=0, line=1, cex.main=1.7)

##add legends
# add invisible boxplot first
barplot(fgs_prop, col="white", width=0.7, ylab=NA, xlab=NA, axes=FALSE, border="white", axisnames=FALSE)
legend(-0.3,1, legend=c("Herb","Inv","Mcar","Omni","Pisc","Plank"),
       col=c("#668F66","#F42A52","#975E93","#144366","#790501","#E9CB20"),
       pch=15, cex=1.1, bty="n", pt.cex=1.4, title="Trophic groups", title.adj=0.4,
       x.intersp=0.6)
legend(0.3,0.6, legend=c("<7","7.1-15","15.1-30","30.1-50","50.1-80",">80 cm"),
       col=c("#E8EFFF","#C3CDFF","#95A0D2","#6975A5","#3F4D79","#142850"), 
       pch=15, cex=1.1, bty="n", pt.cex=1.4, title="Size class", 
       x.intersp=0.6, title.adj=0.2)

######################
######################
#####################
# Binomial tests for the proportion of body sizes
# and trophic levels across time periods

# check number of species within body size classes and trophic levels

bclass #species in body size classes
tgroups #species in trophic groups

colSums(bclass) #check species richness across periods. This number will be used in binomial tests
colSums(tgroups)

#Below, compare the number of species within a body size category in different periods.
#For example, test the difference in the proportion of large bodied species (size 6) 
#between the Ceramic period and the Modern period

#25 large-bodied relative to 34 species in the Ceramic period
#73 large-bodied relative to 317 species in the Modern period

#bsize

#Non-ceramic vs. ceramic
prop.test(x=c(13,3), n=c(62,34), p = NULL, alternative = "two.sided",
          correct = TRUE) #size4
prop.test(x=c(10,6), n=c(62,34), p = NULL, alternative = "two.sided",
          correct = TRUE) #size5
prop.test(x=c(37,25), n=c(62,34), p = NULL, alternative = "two.sided",
          correct = TRUE) #size6

#No significant differences detected for tests above

#Non-ceramic vs. modern
prop.test(x=c(2,14), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #size3
prop.test(x=c(13,23), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #size4
prop.test(x=c(10,13), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #size5
prop.test(x=c(37,43), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #size6

#Results were significant only for Size 3 p-value = 0.03743

#Ceramic vs. modern
prop.test(x=c(3,23), n=c(34,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #size4
prop.test(x=c(6,13), n=c(34,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #size5
prop.test(x=c(25,43), n=c(34,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #size6

#Results were significant only for Size 6 p-value = 0.009833

#######################
#trophic groups

#Non-ceramic vs. ceramic
prop.test(x=c(3,2), n=c(62,34), p = NULL, alternative = "two.sided",
          correct = TRUE) #HERB
prop.test(x=c(25,11), n=c(62,34), p = NULL, alternative = "two.sided",
          correct = TRUE) #INV
prop.test(x=c(24,12), n=c(62,34), p = NULL, alternative = "two.sided",
          correct = TRUE) #MCAR
prop.test(x=c(5,4), n=c(62,34), p = NULL, alternative = "two.sided",
          correct = TRUE) #PISC

#No significant differences detected for tests above

#Non-ceramic vs. modern
prop.test(x=c(3,4), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #HERB
prop.test(x=c(25,42), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #INV
prop.test(x=c(24,34), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #MCAR
prop.test(x=c(5,4), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #PISC
prop.test(x=c(0,4), n=c(62,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #PLANK

#No significant differences detected for tests above


#Ceramic vs. modern
prop.test(x=c(2,4), n=c(34,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #HERB
prop.test(x=c(11,42), n=c(34,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #INV
prop.test(x=c(12,34), n=c(34,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #MCAR
prop.test(x=c(4,4), n=c(34,94), p = NULL, alternative = "two.sided",
          correct = TRUE) #PISC

#No significant differences detected for tests above