#This code was written for the use in the paper Debortoli, N. S., Sayles, J. S., Clark, D. G. & Ford, J. D. A systems network approach for climate change vulnerability assessment. Environ Res Lett 13, 104019 (2018).

#We	 calculated	 Pi and	 Z-scores	 following	 the	 work	 of	 Battiston	 et	 al.	 (2014).	Formal	 mathematical	detentions	are	provided	in	the	main	text of Debertolie et al. 2018. These	were	operationalised	in	the	R	language	environment	using	the	 	 packages	 network	 and	 sna	(Butts	 2015,	 2016).	 The	 functions	 are	 defined	 in	 the	 following	 script.	Authors wishing	to	use	this	function	might	consider	citing	both	this	work	and	Battiston	et	al.	(2014). 

# Citation: Battiston F, Nicosia V and Latora V 2014 Structural measures for multiplex networks Phys. Rev. E 89 032804

# Citation: Debortoli, N. S., Sayles, J. S., Clark, D. G. & Ford, J. D. A systems network approach for climate change vulnerability assessment. Environ Res Lett 13, 104019 (2018).

#	The	function	multiplex_Pi	generates	Pi	and	z	scores	for	multiplex	graph	following	the	work	of	Battiston	 et	 al.	 (2014).	 gL	 is	 a	 list	 of	 graphs	 representing	 multiplex	 edges.	 gmode	 and	 cmode	correspond	to	degree()	in	package	sna.	getki	specifies	if	ki	(or	degree	score	per	node	per	graphlayer)	should	be	returned.	Options	are	"no"	"yes"	and	"only."	Default	is	"no".	No	returns	a data	frame	 of	 Pi,	 z,	 oi	 (total	 degree	 overlap),	and	 oi_2	 (oi2	=	 oi	with	 0<-NA).	Yes	 returns	a	list	 that	includes	data	frames	of	the	Pi/z	results	and	ki.	Only	returns	just	ki.

#Following Battistonet al (2014) who expanded the concept of single layer network partpation to multiplex networks, we consider Pi to be zero when a node has edges of only one kind or in one layer of the multiplex network.Pi equals one when a node has an equal number of edges among all the categories or layers in the network. 

#Following Battiston et al (2014), we also consider Z-scores, which normalize the total overlapping degree (i.e., total degree centrality across all network layer) to allow comparisons of networks of different sizes. Unlike Pi, which is bound by 0 and 1, the Z-score is not bounded and its magnitude and range illustrate the variability of total degree overlap in the system.


library(sna)
library(network)

multiplex_Pi<-function(gL,gmode,cmode,getki="no"){
  ki<-degree(gL,g=c(1:length(gL)),gmode=gmode,cmode=cmode);
  oi<-rowSums(ki);	
  M<-length(gL);
  Pi<-((M/(M-1)))*(1-rowSums((ki/oi)^2));
  oi_2<-oi; oi_2[oi_2==0]<-NA; #recode	0toNA	so	don't	affect	mean	or	sd
  o_avg<-mean(oi_2,na.rm	=	TRUE);	#use	oi_2	&	set	na.rm=TRUE
  o_sd<-sd(oi_2,na.rm	=	TRUE);	#use	oi_2	&	set	na.rm=TRUE
  z<-(oi_2-o_avg)/o_sd;
  name<-get.vertex.attribute(gL[[1]],"vertex.names")	#added	to	include	v.name	in	output		
  
  if(getki=="no"){out<-(data.frame(Pi,z,oi,oi_2,name))}
  if(getki=="yes"){out<-(list(data.frame(Pi,z,oi,oi_2,name),	data.frame(ki)))}
  if(getki=="only"){out<-(data.frame(ki))}
  print(out)	
}

#	In	the	following,	the	input	data	would	be	a	directed	graph	and	analysis	would	consider	both	in	and	out	edges	as	specified	in	the	package	sna. 

# Example sintax: multiplex_Pi(s1Nets,gmode="digraph",cmode="freeman")

###############################################
############# sample data #####################
##############################################


#The	following	are some	small	working	examples.	For	simplicity	the	sample	graphs	are	digraphs	and	calibrated	to	work	with	outdegree	to	produce	Pi values	of	0	and	1.		

g1<-as.network(matrix(c(0,1,1,1,
                        1,0,1,1,
                        1,1,0,1,
                        0,0,0,0),	ncol=4,byrow=TRUE))
set.vertex.attribute(g1,	"vertex.names",	c(LETTERS[1:4]))	


g2<-as.network(matrix(c(0,1,1,1,
                        1,0,1,1,
                        0,0,0,0,
                        0,0,0,0),	ncol=4,byrow=TRUE))
set.vertex.attribute(g2,	"vertex.names",	c(LETTERS[1:4]))


g3<-as.network(matrix(c(0,1,1,1,
                        0,0,0,0,
                        0,0,0,0,
                        0,0,0,0),	ncol=4,byrow=TRUE))
set.vertex.attribute(g3,	"vertex.names",	c(LETTERS[1:4]))


testList<-list(g1,g2,g3)


#Test the data
multiplex_Pi(testList,gmode="digraph",cmode="outdegree")

multiplex_Pi(testList,gmode="digraph",cmode="outdegree",getki="yes")


#References	
#Battiston,	F.,	V.	Nicosia,	and	V.	Latora.	2014.	Structural	measures	for	multiplex	networks.	Physical	Review	E 89(3):32804.
#Butts,	C.	T.	2015.	network:	Classes	for	Relational	Data.	The	Statnet	Project	(<URL:	http://statnet.org>).	R	package	version	1.13.0,	<URL:	http://CRAN.R-project.org/package=network>.
#Butts,	 C.	 T.	 2016.	 sna:	 Tools	 for	 Social	 Network	 Analysis.	 R	 package	 version	 2.4	 https://CRAN.Rproject.org/package=sna.
#Debortoli, N. S., Sayles, J. S., Clark, D. G. & Ford, J. D. A systems network approach for climate change vulnerability assessment. Environ Res Lett 13, 104019 (2018).