Code to make Figure 9
source("Functions For Simulations.R")
n <- 500 #total number of vertices in one graph
kshuff <- c(10, 50, 100, 150, 200) #number shuffled in null
lshuff <- c(2, 5, 10, 25, 50, 100, 150, 200) #number shuffled in alt
prob1 <- matrix(.55, 250, 250) #block1 probabilities
prob2 <- matrix(.45, 250, 250) #block2 probabilities
prob3 <- matrix(.4, 250, 250) #block1 and 2 probabilities
P <- cbind(rbind(prob1, prob3), rbind(prob3, prob2)) #block probability matrix
nMC<-100 #number of repetitions - Monte Carlos
block1 <- c(1:250) #list of vertices in block 1
block2 <- c(251:500) #list of vertices in block 2
null_reg <- list()
alt_reg <- list()
null_reg_unshuff <- list()
alt_reg_unshuff <- list()
set.seed(12345)
error<-c(0.01, 0.05) #error levels
kk<-length(kshuff)
for(i in 1:length(error))
{
null_reg[[i]]<-matrix(0,nMC,kk)
alt_reg[[i]]<-list()
alt_reg[[i]][[1]]<-matrix(0,nMC,3)
alt_reg[[i]][[2]]<-matrix(0,nMC,5)
alt_reg[[i]][[3]]<-matrix(0,nMC,6)
alt_reg[[i]][[4]]<-matrix(0,nMC,7)
alt_reg[[i]][[5]]<-matrix(0,nMC,8)
null_reg_unshuff[[i]]<-matrix(0,nMC,kk)
alt_reg_unshuff[[i]]<-list()
alt_reg_unshuff[[i]][[1]]<-matrix(0,nMC,3)
alt_reg_unshuff[[i]][[2]]<-matrix(0,nMC,5)
alt_reg_unshuff[[i]][[3]]<-matrix(0,nMC,6)
alt_reg_unshuff[[i]][[4]]<-matrix(0,nMC,7)
alt_reg_unshuff[[i]][[5]]<-matrix(0,nMC,8)
}
for (i in 1:length(error))
{
for(h in 1:nMC)
{
for(j in 1:length(kshuff))
{
# make the graphs
X <- sample_correlated_ieg_pair(n, P, P)$graph1
A <- X[]
Xase <- ase(A,2)$X
P1hat<-Xase%*%t(Xase)
P1hat[P1hat<0]<-0
P1hat[P1hat>1]<-1
X2 <- sample_correlated_ieg_pair(n, P, P)$graph1
B <- X2[]
X2ase <- ase(B,2)$X
P2hat<-X2ase%*%t(X2ase)
P2hat[P2hat<0]<-0
P2hat[P2hat>1]<-1
k<-kshuff[j]
ord <- 1:n #original order of vertices
b1_null <- sample(block1, k) #picking vert. in block 1 to shuff in H0
b2_null <- sample(block2, k) #picking vert. in block 2 to shuff in H0
ord_null <- shuff_ord(ord, b1_null, b2_null) #switching vertices for null
P2hatshuff <- P2hat[ord_null,ord_null]
null_reg[[i]][h,j] <- norm(as.matrix(P1hat - P2hatshuff), "F")^2
unshuffledvert <- which(ord == ord_null)
unshuffperm <- gm(P1hat, P2hatshuff, seeds = unshuffledvert)
perm <- as.matrix(unshuffperm$soft)
P2hatunshuff <- perm %*% P2hatshuff %*% t(perm)
null_reg_unshuff[[i]][h,j] <- norm(as.matrix(P1hat - P2hatunshuff), "F")^2
s <- sum(lshuff<= k)
for(q in 1:s)
{
E <- matrix(error[i], n, n)
P_alt<- P + E
X_alt <- sample_correlated_ieg_pair(n, P_alt,P_alt)$graph1
C <- X_alt[]
X_alt_ase <- ase(C,2)$X
P3hat<-X_alt_ase%*%t(X_alt_ase)
P3hat[P3hat<0]<-0
P3hat[P3hat>1]<-1
l<-lshuff[q]
b1_alt <- sample(block1, l)
b2_alt <- sample(block2, l)
ord_alt <- shuff_ord(ord, b1_alt, b2_alt)
P3hatshuff <- P3hat[ord_alt, ord_alt]
alt_reg[[i]][[j]][h,q] <- norm(as.matrix(P1hat - P3hatshuff), "F")^2
unshuffledvert_alt <- which(ord == ord_alt)
unshuffperm_alt <- gm(P1hat, P3hatshuff, seeds = unshuffledvert_alt)
perm <- as.matrix(unshuffperm_alt$soft)
P3hatunshuff <- perm %*% P3hatshuff %*% t(perm)
alt_reg_unshuff[[i]][[j]][h,q] <- norm(as.matrix(P1hat - P3hatunshuff), "F")^2
}
}
}
}
#=================================================================================================================================
#-----------------------------------------Calculating critical values for null--------------------------------------------------
#================================================================================================================================
null_crit_reg<-list()
null_crit_reg_unshuff<-list()
for(i in 1:length(error))
{
null_crit_reg[[i]]<-rep(0,5)
null_crit_reg_unshuff[[i]]<-rep(0,5)
for(j in 1:length(kshuff))
{
null_crit_reg[[i]][j]<-quantile(null_reg[[i]][,j], 0.95)
null_crit_reg_unshuff[[i]][j]<-quantile(null_reg_unshuff[[i]][,j], 0.95)
}
}
for(i in 1:length(error))
{
for(j in 1:length(kshuff))
{
null_crit_reg[[i]][j]<-max(null_crit_reg[[i]][1:j])
}
}
#=================================================================================================================================
#----Calculating power by counting percentage of entries in alt that are greater than their respective critical values in null----
#=================================================================================================================================
Pow_reg <-list()
Pow_reg_se<-list()
Pow_reg_unshuff <-list()
Pow_reg_unshuff_se <-list()
for(i in 1:length(error))
{
Pow_reg[[i]]<-matrix(NA,5,8)
Pow_reg_unshuff[[i]]<-matrix(NA,5,8)
Pow_reg_se[[i]]<-matrix(NA,5,8)
Pow_reg_unshuff_se[[i]]<-matrix(NA,5,8)
for(j in 1:length(kshuff))
{
kk<-kshuff[j]
ll<-sum(lshuff<=kk)
for(q in 1:ll)
{
m1<-sum(alt_reg[[i]][[j]][,q] > null_crit_reg[[i]][j])/nMC
Pow_reg[[i]][j,q]<-m1
Pow_reg_se[[i]][j,q]<-sqrt(m1*(1-m1)/nMC)
m2<-sum(alt_reg_unshuff[[i]][[j]][,q] > null_crit_reg_unshuff[[i]][j])/nMC
Pow_reg_unshuff[[i]][j,q]<-m2
Pow_reg_unshuff_se[[i]][j,q]<-sqrt(m2*(1-m2)/nMC)
}
}
}
#================================================================================================================================
#-----------------------------------------------------------Plotting Power-------------------------------------------------------
#================================================================================================================================
library(ggplot2)
library(reshape2)
library(viridis)
library(colorspace)
library(extrafont)
this <- rev(qualitative_hcl(5, "Dark 3"))
t1 <- "SBM Power Testing after Shuffling,\nError ="
t2 <- "SBM Power Testing after Shuffling\n and Matching, Error ="
plot1 <- list()
plot2 <-list()
for (i in 1:length(Pow_reg))
{
p <- Pow_reg[[i]]
e1 <- as.data.frame(t(rbind(2*lshuff,p)))
colnames(e1) <- c("lshuff", "20", "100", "200", "300", "400")
e11 <- melt(e1, id.vars = "lshuff")
p2 <- Pow_reg_se[[i]]
e2 <- as.data.frame(t(rbind(2*lshuff,p2)))
colnames(e2) <- c("lshuff", "20", "100", "200", "300", "400")
e22 <- melt(e2, id.vars = "lshuff")
colnames(e22)[3]<-"se"
e11$variable <- factor(e11$variable, levels = unique(e11$variable))
e11$se<-e22$se
e <- error[i]
tfinal <- paste(t1, as.character(e), sep=" ")
plot1[[i]] <- ggplot(e11, aes(x=lshuff, y=value, color=variable)) +
geom_line(aes(linetype=variable),size=1.05) +
geom_point(aes(shape=variable),size=4)+ ylim(0,1) +
geom_errorbar(aes(ymin=value-2*se,ymax=value+2*se),width=20)+
labs(family="serif", title=tfinal, x="Number of Vertices Shuffled in H1", y="Power") +
theme_classic() +
theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5),
panel.background = element_rect(fill = "#FFFFFF"),
legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
legend.background = element_rect(fill = "#FFFFFF"),
axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
plot.background = element_rect(fill = "#FFFFFF"),
axis.line = element_line(colour = "#000000", size = 1),
legend.text = element_text( family="serif", color = "#000000", size=15)) +
scale_color_manual(values=c("20"="#B675E0","100"="#00A6CA",
"200"="#00AA5A",
"300"="#AA9000", "400"="#E16A86"),breaks =
levels(e11$variable), name="Number shuffled in H0")+
scale_shape_manual(values=c("20"= 18,"100"=16, "200"= 8, "300"= 17, "400"=
15), breaks = levels(e11$variable),
name="Number shuffled in H0")+
scale_linetype_manual(values=c("20"= 1,"100"=2, "200"= 3, "300"= 4,
"400"=5),breaks = levels(e11$variable),
name="Number shuffled in H0")
}
for (i in 1:length(Pow_reg))
{
p <- Pow_reg_unshuff[[i]]
e1 <- as.data.frame(t(rbind(2*lshuff,p)))
colnames(e1) <- c("lshuff", "20", "100", "200", "300", "400")
e11 <- melt(e1, id.vars = "lshuff")
p2 <- Pow_reg_unshuff_se[[i]]
e2 <- as.data.frame(t(rbind(2*lshuff,p2)))
colnames(e2) <- c("lshuff", "20", "100", "200", "300", "400")
e22 <- melt(e2, id.vars = "lshuff")
colnames(e22)[3]<-"se"
e11$variable <- factor(e11$variable, levels = unique(e11$variable))
e11$se<-e22$se
e <- error[i]
tfinal <- paste(t2, as.character(e), sep=" ")
plot2[[i]] <- ggplot(e11, aes(x=lshuff, y=value, color=variable)) +
geom_line(aes(linetype=variable),size=1.05) +
geom_point(aes(shape=variable),size=4)+ ylim(0,1) +
geom_errorbar(aes(ymin=value-2*se,ymax=value+2*se),width=20)+
labs(family="serif", title=tfinal, x="Number of Vertices Shuffled in H1", y="Power") +
theme_classic() +
theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5),
panel.background = element_rect(fill = "#FFFFFF"),
legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
legend.background = element_rect(fill = "#FFFFFF"),
axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
plot.background = element_rect(fill = "#FFFFFF"),
axis.line = element_line(colour = "#000000", size = 1),
legend.text = element_text( family="serif", color = "#000000", size=15)) +
scale_color_manual(values=c("20"="#B675E0","100"="#00A6CA",
"200"="#00AA5A",
"300"="#AA9000", "400"="#E16A86"),breaks =
levels(e11$variable), name="Number shuffled in H0")+
scale_shape_manual(values=c("20"= 18,"100"=16, "200"= 8, "300"= 17, "400"=
15), breaks = levels(e11$variable),
name="Number shuffled in H0")+
scale_linetype_manual(values=c("20"= 1,"100"=2, "200"= 3, "300"= 4,
"400"=5),breaks = levels(e11$variable),
name="Number shuffled in H0")
}
library(patchwork)
plot1[[1]] + plot1[[2]] + plot2[[1]] + plot2[[2]]+ plot_layout(ncol = 2, guides ="collect")
