#R code for PeerJ manuscript Sergio Diaz Infante et al#
### Interaction network plot and attributes (e.g. for the 1985-1986 period) ###
library(bipartite)
red<-read.delim("chamela8586.txt", row.names = 1)
red 
networklevel(red) 
plotweb(sortweb(red, sort.order="dec"), method="normal", arrow="up", ybig= 1.2, text.rot=90)

################################################################################
### Comparison of a metric value between the two networks ###
#difference between pondered networks, based on Patefield´s nude model (1981)
rm(list= ls())
dados<-read.delim("chamela8586.txt", head=TRUE)
dados
dados2<-read.delim("chamela1617.txt", head=TRUE)
dados2
dados_neo=as.matrix(dados)
dados_neo
dados2_neo=as.matrix(dados2)
dados2_neo

# Choose the metric to be compared between the two networks (e. g. "connectance")

metrics=c("connectance")
metrics
# calculate metrics difference between original networks
orig=abs(networklevel(dados_neo,index=metrics)-networklevel(dados2_neo,index=metrics))
orig

# state permutations number and calculate metrics difference between randomized networks

Nperm = 1000
i=1
randomized.patef=matrix(nrow=length(orig),ncol=Nperm+1)
row.names(randomized.patef)=names(orig)
randomized.patef[,1]=orig
while(i <=Nperm){ 
  dados_aleat=permatfull(dados_neo,fixedmar="both",mtype="count",times=1)
  dados_aleat=dados_aleat$perm[[1]]
  dados2_aleat=permatfull(dados2_neo,fixedmar="both",mtype="count",times=1)
  dados2_aleat=dados2_aleat$perm[[1]]
  linha<-abs(networklevel(dados_aleat, index=metrics)-networklevel(dados2_aleat, index=metrics))
  randomized.patef[,i+1]=linha
  print(i)
  i=i+1
} 
randomized.patef

#Plot comparison between observed value and the randomized metric distribution 
niveis<-row.names(randomized.patef)
for(k in niveis)
{
  if(any(is.na(randomized.patef[k,]) == TRUE))
  {
    print("k tem NA")
  } else {
    plot(density(randomized.patef[k,]), main="Comparison between observed and
         randomized values",)
    abline(v=orig[k], col="red", lwd=2, xlab="")
  }
}

# Calculate randomized difference proportion bigger than the difference between the original networks
significance.patef=matrix(nrow=nrow(randomized.patef),ncol=3)
row.names(significance.patef)=row.names(randomized.patef)
colnames(significance.patef)=c("p (rand <= orig)", "p (rand >= orig)", "p (rand=orig)")

signif.sup=function(x) sum(x>=x[1])/length(x)
signif.inf=function(x) sum(x<=x[1])/length(x)
signif.two=function(x) ifelse(min(x)*2>1,1,min(x)*2)

significance.patef[,1]=apply(randomized.patef,1,signif.inf)
significance.patef[,2]=apply(randomized.patef,1,signif.sup)
significance.patef[,3]=apply(significance.patef[,-3],1,signif.two)

significance.patef

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

### Hummingbird´s niche width ###

library(spaa)

dados<-read.delim("chamela8586.txt", head=TRUE)
dados
dados2<-read.delim("chamela1617.txt", head=TRUE)
dados2
niche.width(dados, method = "shannon")
niche.width(dados2, method = "shannon")

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

### Core-periphery species (e.g. first network corresponding to the 1985-1986 period) ###
library(bipartite)
rede <- read.delim("chamela8586.txt", row.names=1)
rede

# hummingbirds
nucleohigher<-function(x){
  kh<-numeric()
  kh<-specieslevel(x, index="degree",level="higher")/nrow(x)
  kmean<-mean(kh[,1])
  ksd<-sd(kh[,1])
  nucleo<-((kh[,1])-kmean)/ksd  
  especies<-rownames(kh)
  print(data.frame(especies,nucleo)[rev(order(nucleo)),])
}
nucleohigher(rede)

# plants
nucleolower<-function(x){
  kl<-numeric()
  kl<-specieslevel(x, index="degree",level="lower")/ncol(x)
  kmean<-mean(kl[,1])
  ksd<-sd(kl[,1])
  nucleo<-((kl[,1])-kmean)/ksd  
  especies<-rownames(kl)
  print(data.frame(especies,nucleo)[rev(order(nucleo)),])
}
nucleolower(rede)