library(dplyr)
theme_set(theme_bw())

bfsav<-read.csv("RevisedBiteforceData.csv")

bf<-bfsav
bf$rownum<-1:nrow(bf)

fit_hd <- lm((biteforce) ~ (hd), data = bf, na.action=na.exclude)
VALS_fit_hd<-cbind(fit_hd$model,res_hd=fit_hd$residuals)
summary(fit_hd)

fit_svl <- lm((biteforce) ~ (svl), data = bf, na.action=na.exclude)
VALS_fit_svl<-cbind(fit_svl$model,res_svl=fit_svl$residuals)
summary(fit_svl)

fit_hd_log <- lm(log(biteforce) ~ log(hd), data = bf, na.action=na.exclude)
VALS_fit_hd_log<-cbind(fit_hd_log$model,res_hd_log=fit_hd_log$residuals)
summary(fit_hd_log)

fit_svl_log <- lm(log(biteforce) ~ log(svl), data = bf, na.action=na.exclude)
VALS_fit_svl_log<-cbind(fit_svl_log$model,res_svl_log=fit_svl_log$residuals)
summary(fit_svl_log)

vals<-cbind(bf,
            r_hd=residuals(fit_hd),
            r_svl=residuals(fit_svl),
            r_hd_log=residuals(fit_hd_log),
            r_svl_log=residuals(fit_svl_log))

##We focus on log(SVL) because it has a higher adj R sq

#Compare distributions
for(i in c("r_hd","r_svl","r_hd_log","r_svl_log")){
  test1<-ks.test(vals[vals$toothtype=="Acrodont",c(i)],vals[vals$toothtype=="Pleurodont",c(i)],alternative=c("less"))
  print(paste(i,(test1$p.value)))
  }

#log
lab1<-"SVL (mm; log)"
lab2<-"Bite force (N; log)"
lab3<-"SVL-NBF"
labrs<-paste("A ","(Adj. R2 = ",round(summary(fit_svl_log)$adj.r.squared,digits = 2),")",sep="")
vals$xvals<-log(vals$svl)
vals$yvals<-log(vals$biteforce)
vals$res<-vals$r_svl_log

# #not log
lab1<-"SVL (mm)"
lab2<-"Bite force (N)"
lab3<-"SVL-NBF"
labrs<-paste("A ","(Adj. R2 = ",round(summary(fit_svl)$adj.r.squared,digits = 2),")",sep="")
vals$xvals<-(vals$svl)
vals$yvals<-(vals$biteforce)
vals$res<-vals$r_svl

library(patchwork)
library(ggplot2)

p1<-ggplot(vals, aes(x = xvals, y = yvals,fill=toothtype)) + 
  geom_smooth(method="lm",se=F,aes(group=1),color="grey",show.legend=F) + 
  geom_point(shape=21,show.legend=F,alpha=0.6,color="white",size=3)  + 
  ggtitle(labrs)+
  xlab(lab1) + ylab(lab2)

p2<-ggplot(vals,aes(x=toothtype,y=res,color=toothtype))+
  geom_boxplot(show.legend=F,fill="white")+
  geom_jitter(show.legend=F,alpha=0.2) +
  xlab("") + ylab(lab3) +
  ggtitle("B")

p3<-ggplot(vals,aes(x=toothtype,y=res,color=toothtype))+
  geom_boxplot(show.legend=F,fill="white")+
  geom_jitter(show.legend=F,alpha=0.2) +
  xlab("") + ylab(lab3) +
  ggtitle("C")+
  facet_grid(.~diet,scales="free",space="free_x")

p4<-ggplot(vals,aes(x=family,y=res,color=toothtype))+
  geom_boxplot(show.legend=F,fill="white")+
  geom_jitter(show.legend=F,alpha=0.2) +
  xlab("") + ylab(lab3) +
  ggtitle("D")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_grid(.~toothtype,space="free_x",drop=T,scales = "free_x") +
  theme(
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    panel.spacing = unit(0, "lines")
  )

#Seperator line
l<-ggplot()+geom_segment(aes(y=1,yend=1,x=0,xend=1),color="black")+theme_void()

#Plot as 7x7
q1<-(p1+p2)/l+plot_layout(heights=c(5,1),nrow=2)
q2<-p3/p4
q1/q2+plot_layout(nrow=3,heights=c(5,1,10))


####Compare distributions of diets
test_insect<-ks.test(vals[vals$toothtype=="Acrodont" & vals$diet=="Insectivorous","r_hd_log"],vals[vals$toothtype=="Pleurodont" & vals$diet=="Insectivorous","r_hd_log"],alternative=c("less"))
print(test_insect)

test_herb<-ks.test(vals[vals$toothtype=="Acrodont" & vals$diet=="Herbivorous","r_hd_log"],vals[vals$toothtype=="Pleurodont" & vals$diet=="Herbivorous","r_hd_log"],alternative=c("greater"))
print(test_herb)

####Tally of genera----
tal<-vals[,c("species","family","toothtype")]
tal<-unique(tal[,])

#Add missing dat
acrlist<-c("Sphenodontidae","	Chamaeleonidae","Agamidae","Trogonophidae")
tal$toothtype<-ifelse(tal$family=="Sphenodontidae","Acrodont",as.character(tal$toothtype))
tal$toothtype<-ifelse(tal$family=="Chamaeleonidae","Acrodont",as.character(tal$toothtype))
tal$toothtype<-ifelse(tal$family=="Agamidae","Acrodont",as.character(tal$toothtype))
tal$toothtype<-ifelse(tal$family=="Trogonophidae","Acrodont",as.character(tal$toothtype))
tal$toothtype[is.na(tal$toothtype)] <- "Pleurodont"

tal<-tal %>%
  group_by(family,toothtype) %>%
  summarise(famstudy=n())

#phlogenetic order
fam_ord<-c("Sphenodontidae", "Chamaeleonidae", "Agamidae", "Iguanidae", "Dactyloidae", "Crotaphytidae", "Liolaemidae", "Gekkonidae", "Phrynosomatidae", "Lacertidae", "Teiidae", "Cordylidae", "Scincidae", "Xenosauridae", "Anguidae", "Varanidae", "Trogonophidae")
tal$family <- ordered(tal$family, levels = c(fam_ord))

ggplot(tal)+geom_col(aes(x=family,y=famstudy,fill=toothtype)) + 
  xlab("Families") + ylab("Number of species in dataset") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  guides(fill=guide_legend(title="")) 


####Summary tables----
library(dplyr)
t1<- vals %>%
  group_by() %>%
  summarise(#Length
    minL=min(svl,na.rm=T),
    maxL=max(svl,na.rm=T),
    meanL=mean(svl,na.rm=T),
    sdL=sd(svl,na.rm=T),
    #Head
    minH=min(hd,na.rm=T),
    maxH=max(hd,na.rm=T),
    meanH=mean(hd,na.rm=T),
    sdH=sd(hd,na.rm=T),
    #Raw bite force
    minRBF=min(biteforce,na.rm=T),
    maxRBF=max(biteforce,na.rm=T),
    meanRBF=mean(biteforce,na.rm=T),
    sdRBF=sd(biteforce,na.rm=T),
    #Norm bite force SVL LOG
    minNL=min(r_svl_log,na.rm=T),
    maxNL=max(r_svl_log,na.rm=T),
    meanNL=mean(r_svl_log,na.rm=T),
    sdNL=sd(r_svl_log,na.rm=T),
    #Norm bite force HD LOG
    minNH=min(r_hd_log,na.rm=T),
    maxNH=max(r_hd_log,na.rm=T),
    meanNH=mean(r_hd_log,na.rm=T),
    sdNH=sd(r_hd_log,na.rm=T))

t2<- vals %>%
  group_by(toothtype) %>%
  summarise(#Length
    minL=min(svl,na.rm=T),
    maxL=max(svl,na.rm=T),
    meanL=mean(svl,na.rm=T),
    sdL=sd(svl,na.rm=T),
    #Head
    minH=min(hd,na.rm=T),
    maxH=max(hd,na.rm=T),
    meanH=mean(hd,na.rm=T),
    sdH=sd(hd,na.rm=T),
    #Raw bite force
    minRBF=min(biteforce,na.rm=T),
    maxRBF=max(biteforce,na.rm=T),
    meanRBF=mean(biteforce,na.rm=T),
    sdRBF=sd(biteforce,na.rm=T),
    #Norm bite force SVL LOG
    minNL=min(r_svl_log,na.rm=T),
    maxNL=max(r_svl_log,na.rm=T),
    meanNL=mean(r_svl_log,na.rm=T),
    sdNL=sd(r_svl_log,na.rm=T),
    #Norm bite force HD LOG
    minNH=min(r_hd_log,na.rm=T),
    maxNH=max(r_hd_log,na.rm=T),
    meanNH=mean(r_hd_log,na.rm=T),
    sdNH=sd(r_hd_log,na.rm=T))

nn3<-rbind(cbind(toothtype="All",t1),t2)
write.csv(nn3,"SummaryStats_200329.csv")


####Testing ANCOVA----


library(car)

#Variables are highly correlated
summary(lm(log(svl)~log(hd),data=vals))

stat_list<-c()
for(i in c("svl","hd")) {
  
  vals2<-vals
  vals2$chosen_var<-vals2[,c(i)]
  
  #Independence of variables
  anova_model <- aov(log(chosen_var) ~ toothtype, data = vals2)
  summary(anova_model)
  stat_pv<-summary(anova_model)[[1]][["Pr(>F)"]][[1]]
  #pval > 0.05 so Length and ToothType are independent
  
  #Homogeneity of variance
  leveneTest(log(biteforce)~toothtype,data=vals2)
  stat_lev<-leveneTest(log(biteforce)~toothtype,data=vals2)[1,3]
  # pval>0.05 so variances aren't different
  
  #fit ANCOVA model
  ancova_model <- aov(log(biteforce) ~ toothtype + log(chosen_var), data = vals2)
  #view summary of model
  Anova(ancova_model, type="III") 
  #pval for toothtype is small, indicating it has significant effect even after controlling for length
  stat_ano_F<-Anova(ancova_model, type="III")[c("toothtype"),c("F value")]
  stat_ano_P<-Anova(ancova_model, type="III")[c("toothtype"),c("Pr(>F)")]
  
  stat_list<-rbind(stat_list,c(i,stat_pv,stat_lev,stat_ano_F,stat_ano_P))
  
}

colnames(stat_list)<-c("Covariate","IndepenceOfVariables","DifferenceInVariance","AnovaFval","AnovaPval")