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) #######