#START OF CODE####

#load required packages
require(sna)
require(tnet)

#Network regression####
#for question 1: relationships between the intensity of pre- and post-copulatory competition within pairs of competing males

#2006
#load data files
mtime06<-as.matrix(read.table("mtime06.txt", header=T))#temporal overlap
mspace06<-as.matrix(read.table("mspace06.txt", header=T))#spatial overlap
mfnet06sym<-as.matrix(read.table("mfnet06sym.txt", header=T)) #male fighting network, symmetrical
mmnet06sym<-as.matrix(read.table("mmnet06sym.txt", header=T)) #male sperm competition network, symmetrical

#Mean centre each matrix and divide by its own standard deviation (see Methods)
s.mfnet06sym<- (mfnet06sym-mean(mfnet06sym))/sd(mfnet06sym)
s.mtime06 <- (mtime06 -mean(mtime06))/sd(mtime06)
s.mspace06 <- (mspace06 -mean(mspace06))/sd(mspace06)

#combine into a single list
sperm.pred06<- list(s.mfnet06sym,s.mtime06,s.mspace06)

#enter into network regression, with the network of sperm competition as a response
sperm.model06<-netlm(mmnet06sym,sperm.pred06,intercept=T,mode="graph",nullhyp="qapspp")
summary(sperm.model06)

#2013

mtime13<-as.matrix(read.table("mtime13.txt", header=T)) #temporal overlap
mspace13<-as.matrix(read.table("mspace13.txt", header=T)) #spatial overlap
mfnet13sym<-as.matrix(read.table("mfnet13sym.txt", header=T)) #male fighting network, symmetrical
mmnet13sym<-as.matrix(read.table("mmnet13sym.txt", header=T)) #male sperm competition network, symmetrical

#following steps are the same as for 2006

s.mfnet13sym<- (mfnet13sym-mean(mfnet13sym))/sd(mfnet13sym)
s.mtime13 <- (mtime13 -mean(mtime13))/sd(mtime13)
s.mspace13 <- (mspace13 -mean(mspace13))/sd(mspace13)

sperm.pred13<- list(s.mfnet13sym,s.mtime13,s.mspace13)

sperm.model13<-netlm(mmnet13sym,sperm.pred13,intercept=T,mode="graph",nullhyp="qapspp")
summary(sperm.model13)

#Correlations of degree between networks####
#For question 2: the correlation within-individuals in engagement in pre- and post-copulatory competition

#2006
mfnet06asym<-as.matrix(read.table("mfnet06asym.txt", header=T)) #male fighting network, asymmetrical
mmnet06asym<-as.matrix(read.table("mmnet06asym.txt", header=T)) #male sperm competition network, asymmetrical

#calculate the unweighted degree
mmnet06asym.deg<-sna::degree(mmnet06asym,gmode="digraph",cmode="outdegree",ignore.eval=T) 
mfnet06asym.deg<-sna::degree(mfnet06asym,gmode="digraph",cmode="outdegree",ignore.eval=T)

#test the correlation
cor.test(mmnet06asym.deg,mfnet06asym.deg,method="spearman")

#calculate the weighted degree
mmnet06asym.wdeg<-sna::degree(mmnet06asym,gmode="digraph",cmode="outdegree",ignore.eval=F) 
mfnet06asym.wdeg<-sna::degree(mfnet06asym,gmode="digraph",cmode="outdegree",ignore.eval=F)

#test the correlation
cor.test(mmnet06asym.wdeg,mfnet06asym.wdeg,method="spearman") 

#2013
mfnet13asym<-as.matrix(read.table("mfnet13asym.txt", header=T)) #male fighting network, asymmetrical
mmnet13asym<-as.matrix(read.table("mmnet13asym.txt", header=T)) #male sperm competition network, asymmetrical

#following steps are the same as for 2006
#unweighted degree
mmnet13asym.deg<-sna::degree(mmnet13asym,gmode="digraph",cmode="outdegree",ignore.eval=T) 
mfnet13asym.deg<-sna::degree(mfnet13asym,gmode="digraph",cmode="outdegree",ignore.eval=T)

cor.test(mmnet13asym.deg,mfnet13asym.deg,method="spearman")

#weighted degree
mmnet13asym.wdeg<-sna::degree(mmnet13asym,gmode="digraph",cmode="outdegree",ignore.eval=F) 
mfnet13asym.wdeg<-sna::degree(mfnet13asym,gmode="digraph",cmode="outdegree",ignore.eval=F)

cor.test(mmnet13asym.wdeg,mfnet13asym.wdeg,method="spearman") 

#Mating network degree correlation####
#For question 3: if promiscuous males mate with promiscuous females

#WARNING these can take a long time to run especially for 2013.

#2006

mfmnet06<-as.matrix(read.table("mfmnet06.txt", header=T)) #male-female mating network
matlo06<-as.matrix(read.table("matlo06.txt", header=T)) #male-female temporal overlap
matso06<-as.matrix(read.table("matso06.txt", header=T)) #male-female spatial overlap

matlso06<-matlo06 * matso06 #create matrix where cells are non-zero only if the temporal and spatial overlap values are non-zero

length(matlso06[matlso06!=0]) - length(mfmnet06[mfmnet06!=0]) #this tells us how many zeroes
length(mfmnet06[mfmnet06!=0]) #this tells us how many ones

#create 1000 new networks, based on the matrix of spatial*temporal overlap
#but each time some of the links are changed to zero, to give the same density as the original network

n.mats06<-array(0,c(nrow(mfmnet06),ncol(mfmnet06),1000))

for (k in 1:1000) {
  
  for (i in 1:ncol(n.mats06)) {
    for (j in 1:nrow(n.mats06)) {
      if (matlso06[j,i]>0) { n.mats06[j,i,k]<-sample(c(rep(0,347),#number of zeroes
							rep(1,102)#number of ones
								),1)
      }
    }
  }
  
}

#for each of the new networks, loop through and calculate the degree correlation
cors.mf06<-vector(mode="numeric",1000)
for (i in 1:1000){
  m06id.tm<-as.tnet(n.mats06[,,i])
  m06id.f.tm<-as.tnet(cbind(m06id.tm[,2],m06id.tm[,1])) #create a two-mode network of mating, with females as the focus
  
  m06id.f.dj<-degree_tm(m06id.tm,measure="degree")[,2][match(m06id.f.tm[,2],degree_tm(m06id.tm,measure="degree")[,1])]
  #create vector of the male side of each edge, with the degree of the male replacing his ID. Using female focused edgelist hence males are column 2 and the j nodes
  m06id.f.di<-degree_tm(m06id.f.tm,measure="degree")[,2][match(m06id.f.tm[,1],degree_tm(m06id.f.tm,measure="degree")[,1])]
  #create vector of the female side of each edge, with the degree of the female replacing her ID. Using female focused edgelist hence females are column 1 and the i nodes
  cors.mf06[i]<-cor(m06id.f.di,m06id.f.dj,method="spearman")
}
cors.mf06

#calculate the degree correlation in the original network
mf06id.tm<-as.tnet(mfmnet06)
mf06id.f.tm<-as.tnet(cbind(mf06id.tm[,2],mf06id.tm[,1],mf06id.tm[,3]),type="weighted two-mode tnet") 

mf06id.f.dj<-degree_tm(mf06id.tm,measure="degree")[,2][match(mf06id.f.tm[,2],degree_tm(mf06id.tm,measure="degree")[,1])]
mf06id.f.di<-degree_tm(mf06id.f.tm,measure="degree")[,2][match(mf06id.f.tm[,1],degree_tm(mf06id.f.tm,measure="degree")[,1])]
real.cor.mf06<-cor(mf06id.f.di,mf06id.f.dj,method="spearman")

Dbig06 <- sum(cors.mf06 >= real.cor.mf06)
Dbig06 / length(cors.mf06) #permutation p value
#note this may not be exactly the same as quoted in the paper or each time the analysis is run
#This is because the random permutations will produce a slightly different result each time
#the overall conclusions should not be altered.


#2013

mfmnet13<-as.matrix(read.table("mfmnet13.txt", header=T))
matlo13<-as.matrix(read.table("matlo13.txt", header=T))
matso13<-as.matrix(read.table("matso13.txt", header=T))


#This follows the same method as for 2006, expect the number of zeroes and ones is changed to achieve the correct density

matlso13<-matlo13 * matso13

length(matlso13[matlso13!=0]) - length(mfmnet13[mfmnet13!=0]) #this tells us how many zeroes
length(mfmnet13[mfmnet13!=0]) #this tells us how many ones

n.mats13<-array(0,c(nrow(mfmnet13),ncol(mfmnet13),1000))

for (k in 1:1000) {
  
  for (i in 1:ncol(n.mats13)) {
    for (j in 1:nrow(n.mats13)) {
      if (matlso13[j,i]>0) { n.mats13[j,i,k]<-sample(c(rep(0,1021),rep(1,249)),1)
      }
    }
  }
  
}

cors.mf13<-vector(mode="numeric",1000)
for (i in 1:1000){
  m13id.tm<-as.tnet(n.mats13[,,i])
  m13id.f.tm<-as.tnet(cbind(m13id.tm[,2],m13id.tm[,1]))
  m13id.f.dj<-degree_tm(m13id.tm,measure="degree")[,2][match(m13id.f.tm[,2],degree_tm(m13id.tm,measure="degree")[,1])]
  m13id.f.di<-degree_tm(m13id.f.tm,measure="degree")[,2][match(m13id.f.tm[,1],degree_tm(m13id.f.tm,measure="degree")[,1])]
  cors.mf13[i]<-cor(m13id.f.di,m13id.f.dj,method="spearman")
}
cors.mf13

mf13id.tm<-as.tnet(mfmnet13)
mf13id.f.tm<-as.tnet(cbind(mf13id.tm[,2],mf13id.tm[,1],mf13id.tm[,3]),type="weighted two-mode tnet") 
mf13id.f.dj<-degree_tm(mf13id.tm,measure="degree")[,2][match(mf13id.f.tm[,2],degree_tm(mf13id.tm,measure="degree")[,1])]
mf13id.f.di<-degree_tm(mf13id.f.tm,measure="degree")[,2][match(mf13id.f.tm[,1],degree_tm(mf13id.f.tm,measure="degree")[,1])]
real.cor.mf13<-cor(mf13id.f.di,mf13id.f.dj,method="spearman")

Dbig13 <- sum(cors.mf13 >= real.cor.mf13)
Dbig13 / length(cors.mf13) #permutation p value
#This may not be exactly the same as the p-value quoted in the paper, for reasons mentioned above