#############################################################################
#
#					Po river delta data
#
##############################################################################
####A case study from a to z: data from the Po river delta nested random effects
##### The  first part concerns the biodiversity partitioning using package  entropart
#### computing indices, plotting diversity profiles etc.
###### Read data from csv files, a faunistic table with densities of individuals as abundance measure

delta=read.table("delta1.csv",header=T,sep=";",row.names=1)
env=read.table("env1.csv",header=T,sep=";",row.names=1)

str(delta)
str(env)

######## Building metacommunities
require(entropart)

### Communities correspond to monitoring stations, 23 monitoring stations organized in 10 lagoons divided in 3 ssystems

delta.mc1=MetaCommunity(delta,apply(delta,2,sum))
plot(delta.mc1,col=rainbow(48))
summary(delta.mc1)
#############computing indexes 
#############Entropy profiles at the meta community level with and without Bias correction. 
qq=seq(0,2,by=0.1)
delta.ts=delta.bc.ts=rep(0,length(qq))
for(i in 1:length(qq)){
delta.ts[i]=Tsallis(Ps=delta.mc1$Ps,q=qq[i])
delta.bc.ts[i]=bcTsallis(Ns=delta.mc1$Ns,q=qq[i])
}
################# plotting entropies 
plot(qq,delta.ts,ylim=range(c(delta.ts,delta.bc.ts)),pch=20,type="b",xlab="order",ylab="Tsallis entropy")
 lines(qq,delta.bc.ts,type="b",col=3,lty=5)
 legend("topright",c("With Bias correction","Without Bias Correction"),lty=1,col=c(3,1))
 ###### Notice that there are almost no differences when Coverage is almost 1, i.e no singleton
 Coverage(Ns=delta.mc1$Ns)
 ####################### 
##### we can move from entropies to Hill numbers using the deformed exponential:
delta.hill.ts=delta.ts
for(i in 1:length(qq)){
	delta.hill.ts[i]=expq(delta.ts[i],qq[i])
}
############################
# Let's have a look at the results
######################
plot(qq,delta.hill.ts,type="b",xlab="order",ylab="Hill number")
delta.hill.ts
####################################
#Biodiversity partitioning
####################################

delta.part=DivPart(q=1,delta.mc1,Bias=F)
summary(delta.part)
plot(delta.part,main="q=1")

##################################
#
#Rememeber that diversity profiles can be obtained using the DivProfile function
#
########################################################

divP.delta=DivProfile(qq,delta.mc1)
plot(divP.delta)


###########################
#
# chosing a diversity order we can analyze in further details the values of each diversity measure
#
#######################################

############## Shannon (q=1) - we compute the confidence intervals around each value, using 10000 simulations

divE.delta=DivEst(q=1,MC=delta.mc1,Simulations=1000,Biased=FALSE,Correction="Best")

################### Results visualization ###################
summary(divE.delta)
plot(divE.delta)
##########################################################
#
#					Biodiversity partitioning with communities at lagoon levels
#
############################################################

# Now we modify the communities definitions, instead of working at stations level
# we aggregate stations by taking the average abundance in each lagoon
#### the number of available lagoons can be obtained from the column of lagoons labels in the env dataframe
#### the function levels isolates the different categories in a factor, the number of categories is what we need and we obtain it using length
nl=length(levels(env$Lagoon))
### define a new object
delta.lag=data.frame(matrix(0,nrow(delta),nl))
row.names(delta.lag)=row.names(delta)
colnames(delta.lag)=levels(env$Lagoon)
#### to aggregate by lagoon we apply to each row of the delta dataset the mean
for( i in 1:nrow(delta)){
	delta.lag[i,]=tapply(t(delta[i,]),env$Lagoon,mean)
}
#### create the metacommunity
delta.lag.mc=MetaCommunity(delta.lag,apply(delta.lag,2,sum))
######
#### one species is largely more abbundant than the others
plot(delta.lag.mc,col=colorRampPalette(c("red","green","blue"))(47))
sort(round(apply(delta.lag,1,sum)/sum(apply(delta.lag,1,sum)),5))
### in this way we can see that Streblospio.shrubsolii is largely more abundant in the entire meta-community. 
#We need to see in which communities there are species largely more abundant than the others. First we plot the unweighted communities and we verify which communities are more "populated"

barplot(as.matrix(delta.lag),col=colorRampPalette(c("red","green","blue"))(47),cex.names=0.5)

sort(apply(delta.lag,2,sum))
################
# we can see that Valle di Giralda and Smarlacca are the two lagoons with the most abundant portion of individuals and with species largely more abundant than the others, then by looking directly at the abundances we see that the Alitta succinea (It is the on the third row) is the most abundant species at Smarlacca while Streblospio shrubsolii (it is on the 44th row) is the most abundant in Valle di Giralda. 
#In this two lagoons the probability of an individual to belong to this two species is obtained by dividing the amount of each species in each lagoon by the total abundance of individuals of the species in the same lagoon:
delta.lag[c(3,44),c(3,6)]/c(sum(delta.lag[,3]),sum(delta.lag[,6]))

####### Now we can partition diversity with communities defined by lagoons

divP.delta.lag=DivPart(q=1,delta.lag.mc)
summary(divP.delta.lag)
plot(divP.delta.lag)
#########################
# profiles
#########

divP.delta.lag=DivProfile(qq,delta.lag.mc)
plot(divP.delta.lag)
############# compare stations and lagoons

plot(qq,divP.delta$TotalAlphaDiversity,type="b")
 lines(qq,divP.delta.lag$TotalAlphaDiversity,type="b",col=2)
##########################################################
#
#					Biodiversity partitioning with communities at Systems level
#
############################################################
################
# we adopt the same procedure seen for the lagoons, this time we use the System column in the env data.frame
ns=length(levels(env$System))
### define a new object
delta.sys=data.frame(matrix(0,nrow(delta),ns))
row.names(delta.sys)=row.names(delta)
colnames(delta.sys)=levels(env$System)
#### to aggregate by lagoon we apply to each row of the delta dataset the mean
for( i in 1:nrow(delta)){
	delta.sys[i,]=tapply(t(delta[i,]),env$System,mean)
}
#### create the metacommunity
delta.sys.mc=MetaCommunity(delta.sys,apply(delta.sys,2,sum))
plot(delta.sys.mc,col=colorRampPalette(c("red","green","blue"))(47))
### differences have been smoothed 

###### Now we can partition diversity with communities defined by systems

divP.delta.sys=DivPart(q=1,delta.sys.mc)
summary(divP.delta.sys)
plot(divP.delta.sys)
#########################
# profiles
#############################

divP.delta.sys=DivProfile(qq,delta.sys.mc)
plot(divP.delta.sys)
############# compare stations,lagoons and systems
#### alpha diversity
plot(qq,divP.delta$TotalAlphaDiversity,type="b", ylab=expression(paste(alpha,"-diversity",sep="")),xlab="diversity order",ylim=range(c(divP.delta$TotalAlphaDiversity,divP.delta.lag$TotalAlphaDiversity,divP.delta.sys$TotalAlphaDiversity)))
 lines(qq,divP.delta.lag$TotalAlphaDiversity,type="b",col=2)
lines(qq,divP.delta.sys$TotalAlphaDiversity,type="b",col=4)
legend("topright",c("station","lagoon","system"),lty=1,col=c(1,2,4))

####gamma diversity
plot(qq,divP.delta$GammaDiversity,type="b", ylab=expression(paste(gamma,"-diversity",sep="")),xlab="diversity order",ylim=range(c(divP.delta$GammaDiversity,divP.delta.lag$GammaDiversity,divP.delta.sys$GammaDiversity)))
 lines(qq,divP.delta.lag$GammaDiversity,type="b",col=2)
lines(qq,divP.delta.sys$GammaDiversity,type="b",col=4)
legend("topright",c("station","lagoon","system"),lty=1,col=c(1,2,4))

######## beta diversity
plot(qq,divP.delta$TotalBetaDiversity,type="b", ylab=expression(paste(beta,"-diversity",sep="")),xlab="diversity order",ylim=range(c(1,divP.delta$TotalBetaDiversity,divP.delta.lag$TotalBetaDiversity,divP.delta.sys$TotalBetaDiversity)))
 lines(qq,divP.delta.lag$TotalBetaDiversity,type="b",col=2)
lines(qq,divP.delta.sys$TotalBetaDiversity,type="b",col=4)
legend("top",c("station","lagoon","system"),lty=1,col=c(1,2,4),cex=0.5)

###################################################################################
#
#		As expected there is a considerable change in the values of alpha and beta diversity when changing 
#		the definition of community (spatial scale)
#
###################################################################################

#########################################################
#
#											 Mixed models 
#
##########################################################

################################################################################
#
# we are interested in understanding the influence of Sediment and Vegetation as fixed effects 			
# and spatial aggregation (system/lagoon/station) on the evaluation of biodiversity. 
# We define and estimate a linear mixed effects model on the Tsallis alpha entropy of order 1 (Shannon)
# computed at station level
#
################################################################################

require(nlme)

#### To ease model writing we define two factors
sys=env$System
lag=env$Lagoon

### these 2 factors are nested and will be used in the modeling of the random effect
#### we extract the alpha entropy at community level from the object we built above
alpha.1=delta.part$CommunityAlphaEntropies

##### now we are ready to estimate the model, to asses coefficients significance we set an error level of 0.1 (10%)
mod0.delta=lme(alpha.1~1,random=~1|sys/lag,data=env,method="ML")
summary(mod0.delta)

mod1.delta=lme(alpha.1~Vegetation+Sediment,random=~1|sys/lag,data=env,method="ML")
summary(mod1.delta)
#corner point is emerged macrophytes and mud, the level macroa (macroalgae presence) differs consistently from the corner point
## some difference is observed for the level none as well while Sediment seems to have no effect
## ---> try different models
###################################
# remove Sediment
#######
mod2.delta=lme(alpha.1~Vegetation,random=~1|sys/lag,data=env,method="ML")
summary(mod2.delta)

################### we build a factor discriminating between presence/absence of Macro Algae so to include again the Sediment in the model
####################################

Macroa=ifelse(env$Vegetation=="macroa",1,0)
Macroa=factor(Macroa)

mod3.delta=lme(alpha.1~Macroa+Sediment,random=~1|sys/lag,data=env,method="ML")
summary(mod3.delta)

mod4.delta=lme(alpha.1~Macroa,random=~1|sys/lag,data=env,method="ML")
comp=data.frame(AIC=c(
AIC(mod0.delta),
AIC(mod1.delta),
AIC(mod2.delta),
AIC(mod3.delta),
AIC(mod4.delta)
), BIC=c(
BIC(mod0.delta),
BIC(mod1.delta),
BIC(mod2.delta),
BIC(mod3.delta),
BIC(mod4.delta)
))
row.names(comp)=c("mod0","mod1","mod2","mod3","mod4")
barplot(as.matrix(comp),beside=T,col=rainbow(5),legend=T,args.legend=list(horiz=T,x="top",cex=0.7),ylim=c(0,45))
### model 2 is "best" according to both AIC and BIC criterion
summary(mod2.delta)
######### verify if residuals of the chosen model are gaussian
qqnorm(residuals(mod2.delta))
qqline(residuals(mod2.delta),col=2)
shapiro.test(residuals(mod2.delta)) ### confirmed residuals are gaussian
#######################################################
#
#		To asses the random effects contribution to the alpha entropy measure
#		observe the standard deviation given in the summary of the model
##################################################
ss=c(3.691927e-06,8.342304e-07)
names(ss)=c("system","lagoon")
barplot(ss,beside=T,col=c(5,3))


#### what happens if we change index order?
######## reproduce the example considering q=2 and q=0