library(cluster)
#library(ape)
#library(sna)
library(lme4)
library(igraph)

########
### You just change these lines ####
########
netsims<-diffsimcont20.2
currentstack<-stack20.2

netsims<-diffsimcont20.5
currentstack<-stack20.5

netsims<-diffsimcont20.10
currentstack<-stack20.10

netsims<-diffsimcont50.2
currentstack<-stack50.2

netsims<-diffsimcont50.5
currentstack<-stack50.5

netsims<-diffsimcont50.10
currentstack<-stack50.10

netsims<-diffsimcont100.2
currentstack<-stack100.2

netsims<-diffsimcont100.5
currentstack<-stack100.5

netsims<-diffsimcont100.10
currentstack<-stack100.10


########

Nmat<-dim(netsims)[3]

tot.res<-matrix(nrow=Nmat,ncol=6)

colnames(tot.res)<-c("interceptBeta","networkBeta","treeBeta","interceptP","networkP","treeP")

int<-rep(1,dim(netsims)[1])



for(i in 1:Nmat) {

#library(sna)
#detach(package:sna)

sn.adj1<-get.adjacency(graph.edgelist(currentstack[,,i],directed=F),sparse=F)

if(i < Nmat) {
	sn.adj2<-get.adjacency(graph.edgelist(currentstack[,,i+1],directed=F),sparse=F)
		} else {
	sn.adj2<-get.adjacency(graph.edgelist(currentstack[,,i-1],directed=F),sparse=F)
			}
	
#sn.rnorm1<-sn.adj1/rowSums(sn.adj1)
#sn.rnorm1[is.na(sn.rnorm1)]<-0

#sn.rnorm2<-sn.adj2/rowSums(sn.adj2)
#sn.rnorm2[is.na(sn.rnorm2)]<-0

sn.adj1<-1-sn.adj1
sn.adj2<-1-sn.adj2


names.vector<-1:nrow(sn.adj1)

rows<-matrix(rep(names.vector,ncol(sn.adj1)),ncol=ncol(sn.adj1))
cols<-matrix(rep(names.vector,ncol(sn.adj1)),ncol=ncol(sn.adj1),byrow=T)

outcome.vector<-as.vector(daisy(as.data.frame(netsims[,,i]),metric="euclidean"))

temp.data<-data.frame(outcome.vector,sn.adj1[lower.tri(sn.adj1)],sn.adj2[lower.tri(sn.adj2)],rows[lower.tri(rows)],cols[lower.tri(cols)])

colnames(temp.data)<-c("outcome","sn.adj1","sn.adj2","rows","cols")

net.mod<-lmer(outcome~sn.adj1+sn.adj2+(1|rows)+(1|cols),data=temp.data)

### This line will test the performance of simple linear regression
#net.mod<-lm(outcome~sn.adj1+sn.adj2,data=temp.data)
#####

test<-summary(net.mod)$coefficients
tot.res[i,1]<-test["(Intercept)","Estimate"]
tot.res[i,2]<-test["sn.adj1","Estimate"]
tot.res[i,3]<-test["sn.adj2","Estimate"]

tot.res[i,4]<-2*(1-pnorm(abs(test["(Intercept)","Estimate"]
/test["(Intercept)","Std. Error"])))
 
tot.res[i,5]<-2*(1-pnorm(abs(test["sn.adj1","Estimate"]
/test["sn.adj1","Std. Error"])))

tot.res[i,6]<-2*(1-pnorm(abs(test["sn.adj2","Estimate"]
/test["sn.adj2","Std. Error"]))) 

			}


tot.res<-as.data.frame(tot.res)

tot.res<-na.omit(tot.res)
nrow(tot.res)

power<-nrow(tot.res[tot.res$networkP<0.05 & tot.res$networkBeta>0,])/nrow(tot.res)
t1e.tree<-nrow(tot.res[tot.res$treeP<0.05,])/nrow(tot.res)
t1e.network<-nrow(tot.res[tot.res$networkP<0.05 & tot.res$networkBeta<0,])/nrow(tot.res)

power
t1e.network
t1e.tree

########
### You just change these lines ####
########
tot.res.cont20.2<-na.omit(tot.res)
tot.res.cont20.5<-na.omit(tot.res)
tot.res.cont20.10<-na.omit(tot.res)
tot.res.cont50.2<-na.omit(tot.res)
tot.res.cont50.5<-na.omit(tot.res)
tot.res.cont50.10<-na.omit(tot.res)
tot.res.cont100.2<-na.omit(tot.res)
tot.res.cont100.5<-na.omit(tot.res)
tot.res.cont100.10<-na.omit(tot.res)
#######