###############ANALYSIS OF VARIANCE####
###upload of  "tabglm"

tabglm$class<-as.factor(tabglm$class)  
tabglm$sector<-as.factor(tabglm$sector) 
tabglm$site<-as.factor(tabglm$site)
tabglm$habitat<-as.factor(tabglm$habitat)


####sea urchin density under commercial size (< 5cm diameter test) ################################################

tabtot<-subset(tabglm, class=="undersize")
summary(tabtot)

library(nortest)
pearson.test (tabtot$den)###normale

#####anova non bilanciata x habitat dei settori
#####nested fixed factors senza random factor  ###DEF
library(lme4)


####con  habitat (OK)
fit1<-glm (den~sector+habitat, data=tabtot)
summary(fit1)
anova(fit1, test="F") ###The functions summary and anova are used to obtain and print a summary and analysis of variance table of the results.

mcheck<-function (obj, ... ) {   
  rs<-obj$resid
  fv<-obj$fitted
  par(mfrow=c(1,2))
  plot(fv,rs,xlab="Fitted values",ylab="Residuals")
  abline(h=0, lty=2)
  qqnorm(rs,xlab="Normal scores",ylab="Ordered residuals")
  qqline(rs,lty=2)
  par(mfrow=c(1,1))
  invisible(NULL)}

mcheck(fit1)####ok- supplementary material


###variability of recruit density  #####################
tabrec<-subset(tabglm, class=="a")
summary(tabrec)
pearson.test (tabrec$den)###non-normal distribution 


###Percentage of zeros
100*sum(tabrec$den == 0) / nrow(tabrec)###34% very high
# Based on our experience, you should raise your eyebrow
# if you have more than 20% of zeros. 


summary(tabrec)
plot(tabrec$den1)
dotchart(tabrec$den1, ylab = "Order of observations", xlab = "Density x plot", main = "rec density")
hist(tabrec$den1)####
####no outliers but zero inflation

####binomiale negativa zero inflation
install.packages("glmmADMB", repos="http://www.math.mcmaster.ca/bolker/R",type="source")

install.packages("glmmADMB", repos=c("http://glmmadmb.r-forge.r-project.org/repos", getOption("repos")))
library(glmmADMB)

###omit sector 5 --> perfectly multicollinear between sector and habitat 
tabrec2<-subset(tabrec, settore!="5")
summary(tabrec2)

####glm
fit2b<-glmmadmb(den1~settore+habitat, family="nbinom", zeroInflation=T, data=tabrec2)
summary(fit2b)###supplementary material
mcheck(fit2b)###supplementary material
p<-aov(fit2b)
summary(p)####summary result table


#####variability of middle-size sea urchin density
tabmed<-subset(tabglm,class=="young" )
summary(tabmed)
pearson.test (tabmed$den)###normal distribution 

fitmed<-glm(den~settore+habitat, data=tabmed) ###
summary(fitmed)##supplementary material
anova(fitmed, test="F")####summary result table 
mcheck(fitmed)####supplementary material 


##########################################################################


#####CORRELATION RECRUIT CLASS 0-1cm AND BOTTOM CURRENT SPEED

#### upload tabrec
###



tabrec$id<-as.factor(tabrec$id)
tabrec$habitat<-as.factor(tabrec$habitat)  
tabrec$subarea<-as.factor(tabrec$subarea) 
tabrec$stazione<-as.factor(tabrec$stazione) 
tabrec$zona<-as.factor(tabrec$zona) 
tabrec$sito<-as.factor(tabrec$sito)
tabrec$prof<-as.factor(tabrec$prof)
#################






#####correlation rec class a and bottom current speed
summary(tabrec)

library(nortest)
###variabile risposta non normale
shapiro.test(tabrec$a) ####non normale
pearson.test(tabrec$corrente) ####normale

#######grafico  regressione

p1<-ggplot(tabrec, aes(x=corrente, y=a))  
p1<-p1+ geom_point(aes(color=subarea), size=2)
prec<-p1+ geom_smooth(method="lm", se=FALSE, fullrange=TRUE, color="black", span=0.5, lty=2)+
  labs( x="current speed (m/s)", y = "recruit density")+
  scale_y_continuous(limits=c(0, 10))+ scale_x_continuous(limits=c(0.03,0.17))+
  theme(panel.background = element_rect(fill = "white", colour = "grey50"))
prec


#####REGRESSIONE SPEARMAN TEST (REC NO SEGUONO DISTRIBUZIONE NORMALE)
cor.test(tabrec$corrente, tabrec$a,method="spearman") ####uso spearman perch? den reclute non ? normale (pearson)
###p-value = 0.002932
###correlation coefficient (Cor.coeff = -0.4)

####predazione

####tabA




tabpredef$class<-as.factor(tabpredef$class)  
tabpredef$zona<-as.factor(tabpredef$zona) 
tabpredef$anno<-as.factor(tabpredef$anno) 
tabpredef$settore<-as.factor(tabpredef$settore) 
tabpredef$sito<-as.factor(tabpredef$sito)

summary(tabpredef)


############PESCI PREDATORI#####

###pesci predatori

####densit?? totpred
totpred<-summarySE(tabpred, measurevar="pred_den", groupvars=c("sector"))
totpred###


p1a<-ggplot(totpred, aes(x=sector, y=pred_den)) + 
  geom_errorbar(aes(ymin=pred_den, ymax=pred_den+se),width=.2,position=position_dodge(.9))+
  geom_bar(position=position_dodge(),stat="identity") 

p1a<-p1a + labs( x="sector", y = "abbondance", title="totpred")+ scale_fill_manual(values=cbPalette) + theme_classic()+ylim(0,100)

p1a



####################################CORRELATION MIDDLE-SIZE SEA URCHINS AND PREDATORY FIH


#####correlation class CDE-predatory fish #####

tabpregen <- subset(tabpred, class=="cde")

p5<-ggplot(tabpregen, aes(x=pred_den, y=urchin_den))  
p5<-p5+ geom_point(aes(color=sector), size=2) 
p6<-p5+ geom_smooth(method="lm", se=FALSE, fullrange=TRUE, color='black', span=0.5, lty=2)+
  labs( x="fish predators", y = "urchins")+
  scale_y_continuous(limits=c(0, 10))+ scale_x_continuous(limits=c(0, 200))+
  theme(panel.background = element_rect(fill = "white", colour = "grey50"))
p6

library(nortest)
shapiro.test(tabpregen$pred_den) ####non normale
pearson.test(tabpregen$pred_den)  ###normale###non significativo x il terzo decimale

shapiro.test(tabpregen$urchin_den)####normale

cor.test(tabpregen$urchin_den,tabpregen$pred_den,method="pearson")###p-value = 0.04268


###########################################################
####correlation sparidi-classe c maldiventre 
##classe cde
tabsub4<-subset(tabA, subarea=="4")
tabsub4$cd <- tabsub4$c + tabsub4$d

p4<-ggplot(tabsub4, aes(x=sparidi, y=cd)) + 
  geom_point(color='black', size = 1) + 
  geom_smooth(method="lm", se=FALSE, fullrange=TRUE, color='red', span=0.5, lty=2)
p4

cor.test(tabsub4$cd,tabsub4$sparidi,method="pearson")   ###  p-value = 0.02621  



#####GLM SEASCAPE



tabglm$id<-as.factor(tabglm$id)
tabglm$substrato<-as.factor(tabglm$substrato)  
tabglm$settore<-as.factor(tabglm$settore) 
tabglm$sito<-as.factor(tabglm$sito)
tabglm$class<-as.factor(tabglm$class)

summary(tabglm)


library(lme4)

####mcheck 
mcheck<-function (obj, ... ) {   
  rs<-obj$resid
  fv<-obj$fitted
  par(mfrow=c(1,2))
  plot(fv,rs,xlab="Fitted values",ylab="Residuals")
  abline(h=0, lty=2)
  qqnorm(rs,xlab="Normal scores",ylab="Ordered residuals")
  qqline(rs,lty=2)
  par(mfrow=c(1,1))
  invisible(NULL)}

mcheck.glmer<-function (obj, ... ) {   
  rs<-residuals(obj)
  fv<-fitted(obj)
  par(mfrow=c(1,2))
  plot(fv,rs,xlab="Fitted values",ylab="Residuals")
  abline(h=0, lty=2)
  qqnorm(rs,xlab="Normal scores",ylab="Ordered residuals")
  qqline(rs,lty=2)
  par(mfrow=c(1,1))
  invisible(NULL)}

#
####modello  DEFINITIVO CON SOTTOTAGLIE COMMERCIALI ######
tabundersize<-subset(tabglm, class=="undersize")

####glm per undersize in roccia
M0<-glm(den1 ~MPS+PD+IJI+LPI+Partio,family=poisson, data=tabundersize)
summary(M0)
mcheck(M0)### AIC 295.57

M0<-glm(den1 ~MPS+PD+LPI+Partio,family=poisson, data=tabundersize)
summary(M0)
mcheck(M0)### AIC 293.62

M1<-glm(den1 ~MPS+PD,family=poisson, data=tabundersize)
summary(M1)
mcheck(M1)###AIC 290.82 scede, 
####MPS *** positivo
####PD*** positivo

anova(M1,M0, test="Chisq") ###ok diff non significativa tra i due modelli

with(summary(M1), 1 - deviance/null.deviance) ####R-squared: 0.468 Best Model
with(summary(M0), 1 - deviance/null.deviance)####R-squared: 0.476 Full model



##############################################