Code for creating Figure 4
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(1, 5, 10, 25, 50, 100, 150, 200) #number shuffled in alt
P1 <-matrix(c(0.55, 0.4,0.4, 0.45),nrow=2)
P2 <-matrix(c(0.6, 0.35,0.35, 0.5),nrow=2)
block1 <- c(1:250) #list of vertices in block 1
block2 <- c(251:500) #list of vertices in block 2
nMC<-500
#### get the critical values
null_Adj <- matrix(0, nMC, length(kshuff))
null_Phat <- matrix(0, nMC, length(kshuff))
set.seed(123)
for(i in 1:nMC){
# make the graphs
X <- sample_sbm(500,P1,c(250,250))
A <- X[]
xhat<- ase(A,2)$X
p1hat <-xhat%*%t(xhat)
X <- sample_sbm(500,P1,c(250,250))
B <- X[]
yhat <- ase(B,2)$X
p2hat <- yhat%*%t(yhat)
nn <-choose(n,2)*2
aa<-norm(as.matrix(A), "F")^2*(nn-norm(as.matrix(A), "F")^2)/nn
bb<-norm(as.matrix(B), "F")^2*(nn-norm(as.matrix(B), "F")^2)/nn
for(j in 1:length(kshuff))
{
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
Balt <- B[ord_null,ord_null]
null_Adj[i,j] <- norm(as.matrix(A - Balt), "F")^2/(aa+bb)
null_Phat[i,j] <- norm(as.matrix(p1hat - p2hat[ord_null,ord_null]), "F")^2
}
}
c95<-function(x) quantile(x,0.95)
crit_adj <- apply(null_Adj,2 , c95)
crit_Phat <- apply(null_Phat,2 , c95)
for(i in 1:5){
crit_adj[i]<-max(crit_adj[1:i])
crit_Phat[i]<-max(crit_Phat[1:i])
}
# alt for Palt
alt1_Adj <- matrix(0, nMC, length(lshuff))
alt1_Phat <- matrix(0, nMC, length(lshuff))
for(i in 1:nMC){
# make the graphs
X <- sample_sbm(500,P1,c(250,250))
A <- X[]
xhat<- ase(A,2)$X
p1hat <-xhat%*%t(xhat)
X <- sample_sbm(500,P2,c(250,250))
B <- X[]
yhat <- ase(B,2)$X
p2hat <- yhat%*%t(yhat)
nn <-choose(n,2)*2
aa<-norm(as.matrix(A), "F")^2*(nn-norm(as.matrix(A), "F")^2)/nn
bb<-norm(as.matrix(B), "F")^2*(nn-norm(as.matrix(B), "F")^2)/nn
for(j in 1:length(lshuff))
{
k<-lshuff[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
Balt <- B[ord_null,ord_null]
alt1_Adj[i,j] <- norm(as.matrix(A - Balt), "F")^2/(aa+bb)
alt1_Phat[i,j] <- norm(as.matrix(p1hat - p2hat[ord_null,ord_null]), "F")^2
}
}
### Computing power
Pow_alt1_Adj<-matrix(NA,5,8)
Pow_alt1_Adj_se<-matrix(NA,5,8)
Pow_alt1_Phat<-matrix(NA,5,8)
Pow_alt1_Phat_se<-matrix(NA,5,8)
for(i in 1:length(kshuff))
{
kk<-kshuff[i]
ll<-sum(lshuff<=kk)
for(q in 1:ll)
{
m1<-sum(alt1_Adj[,q] >crit_adj[i])/nMC
Pow_alt1_Adj[i,q]<-m1
Pow_alt1_Adj_se[i,q]<-sqrt(m1*(1-m1)/nMC)
m2<-sum(alt1_Phat[,q] > crit_Phat[i])/nMC
Pow_alt1_Phat[i,q]<-m2
Pow_alt1_Phat_se[i,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"))
powers <- list(Pow_alt1_Adj, Pow_alt1_Phat)
st_err <- list(Pow_alt1_Adj_se, Pow_alt1_Phat_se)
ptitle <- list("Power loss in Normalized Adj. Matrix Test", "Power loss in Phat Test")
pplot <-list()
kshuff<-2*kshuff
for(i in 1:2){
dat <- as.data.frame(t(rbind(2*lshuff,powers[[i]])))
colnames(dat) <- c("lshuff", kshuff)
dat <- melt(dat, id.vars = "lshuff")
dat$Null_Shuff <- factor(dat$variable, levels = unique(dat$variable))
dat2<-as.data.frame(t(rbind(2*lshuff,st_err[[i]])))
colnames(dat2) <- c("lshuff", kshuff)
dat2 <- melt(dat2, id.vars = "lshuff")
colnames(dat2)<- c("lshuff", "variable", "se")
dat$se<-dat2$se
pplot[[i]] <- ggplot(dat, 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", x="Number of Vertices Shuffled in H1", y="Power")+
ggtitle(ptitle[[i]])+
theme_classic() +
theme(plot.title = element_text(family="serif", face = "bold", size = 15,
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(dat$variable), name="Number shuffled in H0")+
scale_shape_manual(values=c("20"= 18,"100"=16, "200"= 8, "300"= 17, "400"=
15), breaks = levels(dat$variable),
name="Number shuffled in H0")+
scale_linetype_manual(values=c("20"= 1,"100"=2, "200"= 3, "300"= 4,
"400"=5),breaks = levels(dat$variable),
name="Number shuffled in H0")
}
library(patchwork)
pplot[[1]] + pplot[[2]] + plot_layout(ncol = 2, guides ="collect")
