### First read in the data files lang.d<-read.csv("LanguageDistances.csv",row.names=1,header=T,comment="",quote="\"",stringsAsFactors=F,check.names=F) space.d<-read.csv("Geographic_Distance.csv",row.names=1,header=T,comment="",quote="\"",stringsAsFactors=F,check.names=F) space.adj<-read.csv("geo_adjacency_matrix.csv",row.names=1,header=T,comment="",quote="\"",stringsAsFactors=F,check.names=F) ### check names are full intersection and matrices are aligned nrow(lang.d) ncol(lang.d) nrow(space.d) ncol(space.d) nrow(space.adj) ncol(space.adj) length(intersect(row.names(lang.d),colnames(lang.d))) length(intersect(row.names(lang.d),row.names(space.d))) length(intersect(row.names(lang.d),colnames(space.d))) length(intersect(row.names(lang.d),row.names(space.adj))) length(intersect(row.names(lang.d),colnames(space.adj))) space.d<-space.d[row.names(lang.d),colnames(lang.d)] space.adj<-space.adj[row.names(lang.d),colnames(lang.d)] lang.adj<-lang.d lang.adj[lang.d>1]<-0 lang.adj[lang.adj>0]<-1 lang.d.norm<-(max(lang.d)-lang.d)/max(lang.d) lang.d.norm[is.na(lang.d.norm)]<-0 diag(lang.d.norm)<-0 space.d.norm<-(max(space.d)-space.d)/max(space.d) space.d.norm[is.na(space.d.norm)]<-0 diag(space.d.norm)<-0 diag(lang.adj)<-0 lang.adj.rnorm<-lang.adj/rowSums(lang.adj) lang.adj.rnorm[is.na(lang.adj.rnorm)]<-0 diag(space.adj)<-0 space.adj.rnorm<-space.adj/rowSums(space.adj) space.adj.rnorm[is.na(space.adj.rnorm)]<-0 #source("FunctionSimCharOnNetworkContinuousModified.txt") #sim.contin.on.network(file.name="lang.d",network=lang.d.norm,influence=0.2,Nchar=1,Ntaxa=nrow(lang.d.norm),Nmat=1000) #sim.contin.on.network(file.name="space.d",network=space.d.norm,influence=0.2,Nchar=1,Ntaxa=nrow(space.d.norm),Nmat=1000) #sim.contin.on.network(file.name="lang.adj",network=lang.adj.rnorm,influence=0.2,Nchar=1,Ntaxa=nrow(lang.adj.rnorm),Nmat=1000) #sim.contin.on.network(file.name="space.adj",network=space.adj.rnorm,influence=0.2,Nchar=1,Ntaxa=nrow(space.adj.rnorm),Nmat=1000) ### Here you pick which sims to load #load("lang.adj_network.sims.array.continuous") #load("space.adj_network.sims.array.continuous") #load("lang.d_network.sims.array.continuous") #load("space.d_network.sims.array.continuous") netsims<-allsims Nmat<-dim(netsims)[3] netsims<-allsims dim(netsims) pc.comps<-matrix(nrow=dim(netsims)[1],ncol=Nmat) for(i in 1:Nmat) { pc.comps[,i]<-netsims[,,i] } bic.la.ld.sa.sd<-c() bic.la.ld.sa<-c() bic.la.ld.sd<-c() bic.la.sa.sd<-c() bic.ld.sa.sd<-c() bic.la.ld<-c() bic.la.sa<-c() bic.ld.sa<-c() bic.sa.sd<-c() bic.sd.la<-c() bic.sd.ld<-c() bic.la<-c() bic.ld<-c() bic.sa<-c() bic.sd<-c() bic.null<-c() library(sna) int<-rep(1,nrow(pc.comps)) for(i in 1:250) { net.reg.la.ld.sa.sd<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm,lang.d.norm,space.adj.rnorm,space.d.norm)) net.reg.la.ld.sa<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm,lang.d.norm,space.adj.rnorm)) net.reg.la.ld.sd<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm,lang.d.norm,space.d.norm)) net.reg.la.sa.sd<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm,space.adj.rnorm,space.d.norm)) net.reg.ld.sa.sd<-lnam(pc.comps[,i],int,W2=list(lang.d.norm,space.adj.rnorm,space.d.norm)) net.reg.la.ld<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm,lang.d.norm)) net.reg.la.sa<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm,space.adj.rnorm)) net.reg.ld.sa<-lnam(pc.comps[,i],int,W2=list(lang.d.norm,space.adj.rnorm)) net.reg.sa.sd<-lnam(pc.comps[,i],int,W2=list(space.adj.rnorm,space.d.norm)) net.reg.sd.la<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm,space.d.norm)) net.reg.sd.ld<-lnam(pc.comps[,i],int,W2=list(lang.d.norm,space.d.norm)) net.reg.la<-lnam(pc.comps[,i],int,W2=list(lang.adj.rnorm)) net.reg.ld<-lnam(pc.comps[,i],int,W2=list(lang.d.norm)) net.reg.sa<-lnam(pc.comps[,i],int,W2=list(space.adj.rnorm)) net.reg.sd<-lnam(pc.comps[,i],int,W2=list(space.d.norm)) bic.la.ld.sa.sd[i]<-(-2)*net.reg.la.ld.sa.sd$lnlik.model+log(nrow(pc.comps))*net.reg.la.ld.sa.sd$df.model bic.la.ld.sa[i]<-(-2)*net.reg.la.ld.sa$lnlik.model+log(nrow(pc.comps))*net.reg.la.ld.sa$df.model bic.la.ld.sd[i]<-(-2)*net.reg.la.ld.sd$lnlik.model+log(nrow(pc.comps))*net.reg.la.ld.sd$df.model bic.la.sa.sd[i]<-(-2)*net.reg.la.sa.sd$lnlik.model+log(nrow(pc.comps))*net.reg.la.sa.sd$df.model bic.ld.sa.sd[i]<-(-2)*net.reg.ld.sa.sd$lnlik.model+log(nrow(pc.comps))*net.reg.ld.sa.sd$df.model bic.la.ld[i]<-(-2)*net.reg.la.ld$lnlik.model+log(nrow(pc.comps))*net.reg.la.ld$df.model bic.la.sa[i]<-(-2)*net.reg.la.sa$lnlik.model+log(nrow(pc.comps))*net.reg.la.sa$df.model bic.ld.sa[i]<-(-2)*net.reg.ld.sa$lnlik.model+log(nrow(pc.comps))*net.reg.ld.sa$df.model bic.sa.sd[i]<-(-2)*net.reg.sa.sd$lnlik.model+log(nrow(pc.comps))*net.reg.sa.sd$df.model bic.sd.la[i]<-(-2)*net.reg.sd.la$lnlik.model+log(nrow(pc.comps))*net.reg.sd.la$df.model bic.sd.ld[i]<-(-2)*net.reg.sd.ld$lnlik.model+log(nrow(pc.comps))*net.reg.sd.ld$df.model bic.la[i]<-(-2)*net.reg.la$lnlik.model+log(nrow(pc.comps))*net.reg.la$df.model bic.ld[i]<-(-2)*net.reg.ld$lnlik.model+log(nrow(pc.comps))*net.reg.ld$df.model bic.sa[i]<-(-2)*net.reg.sa$lnlik.model+log(nrow(pc.comps))*net.reg.sa$df.model bic.sd[i]<-(-2)*net.reg.sd$lnlik.model+log(nrow(pc.comps))*net.reg.sd$df.model bic.null[i]<-(-2)*net.reg.sd$lnlik.null+log(nrow(pc.comps))*net.reg.sd$df.null print(i) } top.models<-matrix(nrow=length(bic.la.ld.sa.sd),ncol=8) model.name.list<-c("bic.la.ld.sa.sd","bic.la.ld.sa","bic.la.ld.sd","bic.la.sa.sd","bic.ld.sa.sd","bic.la.ld","bic.la.sa","bic.ld.sa","bic.sa.sd","bic.sd.la","bic.sd.ld","bic.la","bic.ld","bic.sa","bic.sd","bic.null") for (i in 1:length(bic.la.ld.sa.sd)) { temp.list<-c(bic.la.ld.sa.sd[i],bic.la.ld.sa[i],bic.la.ld.sd[i],bic.la.sa.sd[i],bic.ld.sa.sd[i],bic.la.ld[i],bic.la.sa[i],bic.ld.sa[i],bic.sa.sd[i],bic.sd.la[i],bic.sd.ld[i],bic.la[i],bic.ld[i],bic.sa[i],bic.sd[i],bic.null[i]) names(temp.list)<-model.name.list top.models[i,1]<-names(temp.list)[which.min(temp.list)] top.models[i,5]<-min(temp.list) temp.list2<-temp.list[temp.list!=min(temp.list)] top.models[i,2]<-names(temp.list2)[which.min(temp.list2)] top.models[i,6]<-min(temp.list2) temp.list3<-temp.list2[temp.list2!=min(temp.list2)] top.models[i,3]<-names(temp.list3)[which.min(temp.list3)] top.models[i,7]<-min(temp.list3) temp.list4<-temp.list3[temp.list3!=min(temp.list3)] top.models[i,4]<-names(temp.list4)[which.min(temp.list4)] top.models[i,4]<-names(temp.list4)[which.min(temp.list4)] top.models[i,8]<-min(temp.list4) } length(table(top.models[,1])) length(table(top.models[,2])) table(top.models[,1]) table(top.models[,2]) median(as.numeric(na.omit(top.models[,6]))-as.numeric(na.omit(top.models[,5]))) median(as.numeric(na.omit(top.models[,7]))-as.numeric(na.omit(top.models[,5]))) #write.table(top.models,file="lang.adj.sims.processed.txt",row.names=F,col.names=F,sep="|",quote=F) #write.table(top.models,file="space.adj.sims.processed.txt",row.names=F,col.names=F,sep="|",quote=F) #write.table(top.models,file="space.d.sims.processed.txt",row.names=F,col.names=F,sep="|",quote=F) #write.table(top.models,file="lang.d.sims.processed.txt",row.names=F,col.names=F,sep="|",quote=F)