#Comparing Modularity Models: 8 group power test library(geomorph) library(EMMLi) source('compare.CR.r') #simulate data with known covariance sim.data<-function(n,s){ a<-matrix(rnorm(n*nrow(s)),ncol=nrow(s)) b<-a%*%chol(s) return(b) } #builds with Sw and Sb with random noise in all components s.build<-function(p,sb,sw,mod.gp){ s<-matrix(sb,nrow=p, ncol=p) for (i in 1:nlevels(mod.gp)){ s[which(mod.gp==levels(mod.gp)[i]),which(mod.gp==levels(mod.gp)[i])]<-sw } low<-s[lower.tri(s)] for (i in 1:length(low)){ low[i]<-rnorm(1,mean=low[i],sd=0.01) } s[lower.tri(s)]<-low s[upper.tri(s)]<-t(s)[upper.tri(s)] diag(s)<-rnorm(p,mean=1,sd=0.01) s<-matrix(Matrix::nearPD(s)$mat,nrow=p,ncol=p) return(s) } #FOR random covariance matrix: 8 module iter<-999 N <- c(50, 100, 250) p <- c(16,32,64) nsimul <- 500 ###### Power Simulation s.b <- 0.6 s.w <- c(0,.025,.05,.075,.1) #NOTE, this is ADDED to SB CR.Power <- CR.Misspec <- EMMLi.Power <- EMMLi.Misspec <- array(NA, dim = c(length(p), length(N),length(s.w))) row.names(CR.Power) <- row.names(CR.Misspec) <- row.names(EMMLi.Power) <- row.names(EMMLi.Misspec) <- paste("p_",p,sep="") colnames(CR.Power) <- row.names(CR.Misspec) <-colnames(EMMLi.Power) <- colnames(EMMLi.Misspec) <- paste("N_",N,sep="") dimnames(CR.Power)[[3]] <- dimnames(CR.Misspec)[[3]] <- dimnames(EMMLi.Power)[[3]] <- dimnames(EMMLi.Misspec)[[3]] <- paste("SW_",s.w,sep="") EMMLi.best <- NULL for (i in 1:length(p)){ for (j in 1:length(N)){ #hypotheses m2 <- gl(2,p[i]/2) m4 <- gl(4,p[i]/4) m8 <- gl(8,p[i]/8) #data for (ii in 1:length(s.w)){ S <- s.build(p[i],sb=s.b,sw=(s.b+s.w[ii]),mod.gp=m8) Y <- lapply(1:nsimul, function(k) sim.data(n=N[j],S)) CR.m2 <- lapply(1:nsimul, function(k) modularity.test(Y[[k]], partition.gp = m2, iter = iter, print.progress = FALSE)) CR.m4 <- lapply(1:nsimul, function(k) modularity.test(Y[[k]], partition.gp = m4, iter = iter, print.progress = FALSE)) CR.m8 <- lapply(1:nsimul, function(k) modularity.test(Y[[k]], partition.gp = m8, iter = iter, print.progress = FALSE)) CR.comp <- lapply(1:nsimul, function(k) compare.CR(CR.null = TRUE, CR.m2[[k]],CR.m4[[k]],CR.m8[[k]])) #Power: true (2-gp) vs. null CR.Power[i,j,ii] <- length(which(lapply(1:nsimul, function(k) CR.comp[[k]]$pairwise.P[1,4]<=0.05 )=="TRUE")) / nsimul #how often is true model better than null? #Misspecification: best vs. true (8-gp) CR.min <- unlist(lapply(1:nsimul, function(k) which.min(CR.comp[[k]]$sample.z) )) CR.p <- unlist(lapply(1:nsimul, function(k) CR.comp[[k]]$pairwise.P[4,CR.min[k]])) CR.Misspec[i,j,ii] <- length(which(CR.p <=0.05)) / nsimul #significant that are not input model ## EMMLi. best model X.1 <-as.character(seq(1:p[i])) models <- data.frame(X.1,m2, m4, m8) R <- lapply(1:nsimul, function(k) cor(Y[[k]])) EMMLi.out <- lapply(1:nsimul, function(k) EMMLi(R[[k]],N[j], models)$results ) best <- unlist(lapply(1:nsimul, function(k) which.min(EMMLi.out[[k]][,5]))) #NOTE: for power, ANY model with 8 modules counted as best (a favorable view for EMMLi. I could be more strict) EMMLi.Power[i,j,ii] <- ifelse(ii==1,(1-(length(which(best== 1)) / nsimul)) ,length(which(best==8 | best==9 | best==10 | best==11)) / nsimul ) EMMLi.Misspec[i,j,ii] <- 1 - EMMLi.Power[i,j,ii] #Misspecification is the converse best.out <- rep(0,11); best.out[sort(unique(best))] <- table(best) EMMLi.best <- rbind(EMMLi.best, c(p[i],N[j],s.w[ii], best.out)) #to ensure that all options are saved print(paste("i =",i)) print(paste("j =",j)) print(paste("ii =",ii)) } } } colnames(EMMLi.best) <- c("p","N","s.w","1","2","3","4","5","6","7","8","9","10","11") write.csv(EMMLi.best,"emmli.best-8gp-.6sb.csv") write.csv(EMMLi.Power,"emmli.power-8gp-.6sb.csv") write.csv(EMMLi.Misspec,"emmli.misspec-8gp-.6sb.csv") write.csv(CR.Power, "CR.power-8gp-.6sb.csv") write.csv(CR.Misspec,"CR.misspec-8gp-.6sb.csv")