Code for simulating figures 2 and 3
source("Functions For Simulations.R")
set.seed(12345)
n <- 500 #total number of vertices in one graph
kshuff <- c(10, 50, 100, 150, 200) #number shuffled in null
lshuff <- c(1, 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
block1 <- c(1:250) #list of vertices in block 1
block2 <- c(251:500) #list of vertices in block 2
nMC<-200 #total number of repetitions
null_reg <- list()
alt_reg <- list()
null_adj <- list()
alt_adj <- list()
error<-c(0.01, 0.02, 0.03, -0.01, -0.02, -0.03) #error levels
kk<-length(kshuff)
for(i in 1:length(error))
{
null_reg[[i]]<-matrix(0,nMC,kk)
null_adj[[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)
alt_adj[[i]]<-list()
alt_adj[[i]][[1]]<-matrix(0,nMC,3)
alt_adj[[i]][[2]]<-matrix(0,nMC,5)
alt_adj[[i]][[3]]<-matrix(0,nMC,6)
alt_adj[[i]][[4]]<-matrix(0,nMC,7)
alt_adj[[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
null_adj[[i]][h,j] <- norm(as.matrix(A - B[ord_null, ord_null]), "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
alt_adj[[i]][[j]][h,q] <- norm(as.matrix(A - C[ord_alt, ord_alt]), "F")^2
}
}
}
}
#====================================================================================================================
#-----------------------------------------Calculating critical values for null---------------------------------------
#====================================================================================================================
null_crit_reg<-list()
null_crit_adj<-list()
for(i in 1:length(error))
{
null_crit_reg[[i]]<-rep(0,5)
null_crit_adj[[i]]<-rep(0,5)
for(j in 1:length(kshuff))
{
null_crit_reg[[i]][j]<-quantile(null_reg[[i]][,j], 0.95)
null_crit_adj[[i]][j]<-quantile(null_adj[[i]][,j], 0.95)
}
}
#=====================================================================================================================
#---------------------------Calculating power by counting percentage of entries in alt -------------------------------
#-------------------------that are greater than their respective critical values in null------------------------------
#=====================================================================================================================
Pow_reg <-list()
Pow_adj <-list()
# Compute least favorable shuffling under H0
for(i in 1:length(error)){
for(j in 1:length(kshuff)){
null_crit_reg[[i]][j]<-max(null_crit_reg[[i]][1:j])
null_crit_adj[[i]][j]<-max(null_crit_adj[[i]][1:j])
}
}
for(i in 1:length(error))
{
Pow_reg[[i]]<-matrix(NA,5,8)
Pow_adj[[i]]<-matrix(NA,5,8)
for(j in 1:length(kshuff))
{
kk<-kshuff[j]
ll<-sum(lshuff<=kk)
for(q in 1:ll)
{
Pow_reg[[i]][j,q]<-sum(alt_reg[[i]][[j]][,q] > null_crit_reg[[i]][j])/nMC
Pow_adj[[i]][j,q]<-sum(alt_adj[[i]][[j]][,q] > null_crit_adj[[i]][j])/nMC
}
}
}
Pow_reg_sd<-list()
Pow_adj_sd<-list()
for(i in 1:length(error)){
Pow_reg_sd[[i]] <- matrix(0,length(kshuff), length(lshuff))
Pow_adj_sd[[i]] <- matrix(0,length(kshuff), length(lshuff))
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_sd[[i]][j,q]<-sqrt(m1*(1-m1)/nMC)
m2<-sum(alt_adj[[i]][[j]][,q] > null_crit_adj[[i]][j])/nMC
Pow_adj_sd[[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"))
t <- "Error ="
t_adj <- "Error = "
plot <- list()
plotadj <- list()
for (i in 1:length(Pow_reg))
{
p <- Pow_reg[[i]]
pp <- Pow_reg_sd[[i]]
e1 <- as.data.frame(t(rbind(2*lshuff,p)))
e2 <- as.data.frame(t(rbind(2*lshuff,pp)))
colnames(e1) <- c("lshuff","20", "100", "200", "300", "400")
colnames(e2) <- c("lshuff","20", "100", "200", "300", "400")
e11 <- melt(e1, id.vars = "lshuff")
e22 <- melt(e2, id.vars = "lshuff")
colnames(e22) <- c("lshuff","variable", "se")
e11$se<-e22$se
e11$variable <- factor(e11$variable, levels = unique(e11$variable))
e <- error[i]
tfinal <- paste(t, as.character(e), sep=" ")
plot[[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")
p_adj <- Pow_adj[[i]]
pp_adj <- Pow_adj_sd[[i]]
e1_adj <- as.data.frame(t(rbind(2*lshuff,p_adj)))
e2_adj <- as.data.frame(t(rbind(2*lshuff,pp_adj)))
colnames(e1_adj) <- c("lshuff", "20", "100", "200", "300", "400")
colnames(e2_adj) <- c("lshuff", "20", "100", "200", "300", "400")
e11_adj <- melt(e1_adj, id.vars = "lshuff")
e22_adj <- melt(e2_adj, id.vars = "lshuff")
colnames(e22_adj) <- c("lshuff","variable", "se")
e11_adj$se<-e22_adj$se
e11_adj$variable <- factor(e11_adj$variable, levels = unique(e11_adj$variable))
t_adjfinal <- paste0(t_adj, as.character(error[i]))
plotadj[[i]] <- ggplot(e11_adj, 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=t_adjfinal, 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_adj$variable), name="Number shuffled in H0")+
scale_shape_manual(values=c("20"= 18,"100"=16, "200"= 8, "300"= 17, "400"= 15),breaks = levels(e11_adj$variable), name="Number shuffled in H0") +
scale_linetype_manual(values=c("20"= 1,"100"=2, "200"= 3, "300"= 4,
"400"=5),breaks = levels(e11_adj$variable),
name="Number shuffled in H0")
}
library(patchwork)
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
##
## area
plot[[1]] + plot[[2]] + plot[[3]] + plot[[4]] + plot[[5]] + plot[[6]] + plot_layout(ncol = 3, guides =
"collect")
plotadj[[1]] + plotadj[[2]] + plotadj[[3]] + plotadj[[4]] + plotadj[[5]] + plotadj[[6]] + plot_layout(ncol = 3, guides =
"collect")
