Code to make Figure 10
source("Functions For Simulations.R")
set.seed(123)
n <- 500 #number of vertices
ord <- 1:n #original order of vertices
nMC <- 100 #number of repetitions
null <- list()
alt <- list()
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
error<-c(0.01, 0.05) #error levels
null_ase <- list()
alt_ase <- list()
null_ase_unshuff <- list()
alt_ase_unshuff <- list()
kk<-length(kshuff)
for(i in 1:length(error))
{
null_ase[[i]]<-matrix(rnorm(100, mean=50, sd=10),nMC,kk)
alt_ase[[i]]<-list()
alt_ase[[i]][[1]]<-matrix(0,nMC,3)
alt_ase[[i]][[2]]<-matrix(0,nMC,5)
alt_ase[[i]][[3]]<-matrix(0,nMC,6)
alt_ase[[i]][[4]]<-matrix(0,nMC,7)
alt_ase[[i]][[5]]<-matrix(0,nMC,8)
null_ase_unshuff[[i]]<-matrix(0,nMC,kk)
alt_ase_unshuff[[i]]<-list()
alt_ase_unshuff[[i]][[1]]<-matrix(0,nMC,3)
alt_ase_unshuff[[i]][[2]]<-matrix(0,nMC,5)
alt_ase_unshuff[[i]][[3]]<-matrix(0,nMC,6)
alt_ase_unshuff[[i]][[4]]<-matrix(0,nMC,7)
alt_ase_unshuff[[i]][[5]]<-matrix(0,nMC,8)
}
for(e in 1:length(error))
{
for(j in 1:length(kshuff))
{
for(h in 1:nMC)
{
k<-kshuff[j]
shuffle <- scramble(ord, k) #shuffling order for null
X1 <- t(sample_dirichlet(n, alpha=c(1,1,1))) #creating latent position matrix X1
A1 <- as.matrix(sample_dot_product(t(X1))) #sampling from X1
A2 <- as.matrix(sample_dot_product(t(X1))) #sampling from X1 for X=Y test
#X2 <- t(sample_dirichlet(n, alpha=c(1,1,1))) #creating latent position matrix X2
#A2 <- as.matrix(sample_dot_product(t(X1))) #sampling from X1 for L(X)=L(Y) test
A2_shuffled <- A2[shuffle, shuffle] #shuffling X2mod
#creating probability matrix via ASE of A1 and A2mod_shuffled
X1_hat <- ase(A1, ncol(X1))$X
X2_hat <- ase(A2_shuffled, ncol(X1))$X
P1_hat <- X1_hat %*% t(X1_hat)
P2hatshuff <- X2_hat %*% t(X2_hat)
#correcting for probability
P1_hat[P1_hat < 0] <- 0
P1_hat[P1_hat > 1] <- 1
P2hatshuff[P2hatshuff < 0] <- 0
P2hatshuff[P2hatshuff > 1] <- 1
null_ase[[e]][h,j] <- TS2(P1_hat, P2hatshuff) #null test statistic for difference in phats
unshuffledvert <- which(ord == shuffle)
unshuffperm <- gm(P1_hat, P2hatshuff, seeds = unshuffledvert)
perm <- as.matrix(unshuffperm$soft)
P2hatunshuff <- perm %*% P2hatshuff %*% t(perm)
null_ase_unshuff[[e]][h,j] <- TS2(P1_hat, P2hatunshuff)
s <- sum(lshuff<= k)
for(q in 1:s)
{
l<-lshuff[q]
shuffle_alt <- scramble(ord, l) #shuffling order for alt
X3 <- X1 #creating latent position matrix X3 for X=Y test
# X3 <- t(sample_dirichlet(n, alpha=c(1,1,1))) #creating latent position matrix X3 for L(X)=L(Y) test
E <- matrix(error[e], dim(X3)[1], dim(X3)[2]) #creating error matrix that will be added to X3
X3error <- X3 + E #adding error to latent positions matrix X3
A3mod <- as.matrix(sample_dot_product(t(X3error))) #sampling from X3mod
A3mod_shuffled <- A3mod[shuffle_alt, shuffle_alt] #shuffling X3mod
X3_hat_shuffled <- ase(A3mod_shuffled, ncol(X1))$X #creating probability matrix via ASE of A3mod_shuffled
P3hatshuff <- X3_hat_shuffled %*% t(X3_hat_shuffled)
#correcting for probability
P3hatshuff[P3hatshuff < 0] <- 0
P3hatshuff[P3hatshuff > 1] <- 1
alt_ase[[e]][[j]][h,q] <- TS2(P1_hat, P3hatshuff) #alt test for difference in phats
unshuffledvert_alt <- which(ord == shuffle_alt)
unshuffperm_alt <- gm(P1_hat, P3hatshuff, seeds = unshuffledvert_alt)
perm <- as.matrix(unshuffperm_alt$soft)
P3hatunshuff <- perm %*% P3hatshuff %*% t(perm)
alt_ase_unshuff[[e]][[j]][h,q] <- norm(as.matrix(P1_hat - P3hatunshuff), "F")^2
}
}
}
}
#==============================================================================================================
#---------------------------Calculating critical values for null-----------------------------------------------
#==============================================================================================================
null_crit_ase<-list()
null_crit_ase_unshuff<-list()
null_crit_ase_ayushi<-list()
null_crit_ase_unshuff_ayushi<-list()
for(i in 1:length(error))
{
null_crit_ase[[i]]<-rep(0,5)
null_crit_ase_unshuff[[i]]<-rep(0,5)
null_crit_ase_ayushi[[i]]<-rep(0,5)
null_crit_ase_unshuff_ayushi[[i]]<-rep(0,5)
for(j in 1:length(kshuff))
{
null_crit_ase[[i]][j]<-quantile(null_ase[[i]][,j], 0.95)
null_crit_ase_unshuff[[i]][j]<-quantile(null_ase_unshuff[[i]][,j], 0.95)
null_crit_ase_ayushi[[i]][j]<-quantile(null_ase[[i]][,j], 0.05)
null_crit_ase_unshuff_ayushi[[i]][j]<-quantile(null_ase_unshuff[[i]][,j], 0.05)
}
}
#=========================================================================================================
#------------------------------------------Calculating power ---------------------------------------------
#=========================================================================================================
Pow_ase <-list()
Pow_ase_unshuff <-list()
Pow_ase_ayushi <-list()
Pow_ase_unshuff_ayushi <-list()
for(i in 1:length(error))
{
Pow_ase[[i]]<-matrix(NA,5,8)
Pow_ase_unshuff[[i]]<-matrix(NA,5,8)
Pow_ase_ayushi[[i]]<-matrix(NA,5,8)
Pow_ase_unshuff_ayushi[[i]]<-matrix(NA,5,8)
for(j in 1:length(kshuff))
{
kk<-kshuff[j]
ll<-sum(lshuff<=kk)
for(q in 1:ll)
{
Pow_ase[[i]][j,q]<-sum(alt_ase[[i]][[j]][,q] > null_crit_ase[[i]][j])/nMC
Pow_ase_unshuff[[i]][j,q]<-sum(alt_ase_unshuff[[i]][[j]][,q] > null_crit_ase_unshuff[[i]][j])/nMC
above <- alt_ase[[i]][[j]][,q] >= null_crit_ase[[i]][j]
below <- alt_ase[[i]][[j]][,q] < null_crit_ase_ayushi[[i]][j]
above_unshuff <- alt_ase[[i]][[j]][,q] >= null_crit_ase_unshuff[[i]][j]
below_unshuff <- alt_ase[[i]][[j]][,q] < null_crit_ase_unshuff_ayushi[[i]][j]
Pow_ase_ayushi[[i]][j,q]<-sum(above,below)/nMC
Pow_ase_unshuff_ayushi[[i]][j,q]<-sum(above_unshuff, below_unshuff)/nMC
}
}
}
#=========================================================================================================
#------------------------------------------Plotting Power-------------------------------------------------
#=========================================================================================================
library(ggplot2)
library(reshape2)
library(viridis)
library(colorspace)
library(extrafont)
library(RColorBrewer)
library(ggforce)
library(ggbreak)
library(patchwork)
this <- rev(qualitative_hcl(5, "Dark 3"))
t <- "Power Testing after Shuffling RDPG with Error ="
plot <- list()
for (i in 1:length(Pow_ase))
{
p <- Pow_ase[[i]]
e1 <- as.data.frame(t(rbind(lshuff,p)))
colnames(e1) <- c("lshuff", "10", "50", "100", "150", "200")
e11 <- melt(e1, id.vars = "lshuff")
e11$variable <- factor(e11$variable, levels = unique(e11$variable))
e <- error[i]
tfinal <- paste(t, as.character(e), sep=" ")
if(i == 1 | i == 3){yhigh <- 0.4}else{yhigh <- 1}
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, yhigh) +
geom_errorbar(aes(ymin=value-2*sqrt(value*(1-value)/nMC),
ymax=value+2*sqrt(value*(1-value)/nMC)),width=10)+ ylim(0, yhigh)+
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("10"="#B675E0","50"="#00A6CA", "100"="#00AA5A", "150"="#AA9000", "200"="#E16A86"),breaks =
levels(e11$variable))+
scale_shape_manual(values=c("10"= 18,"50"=16, "100"= 8, "150"= 17, "200"= 15),breaks = levels(e11$variable))
}
plot
## [[1]]
##
## [[2]]
t_unshuff <- "Power Testing after Matching RDPG with Error ="
plot_unshuff <- list()
for (i in 1:length(Pow_ase_unshuff))
{
p <- Pow_ase_unshuff[[i]]
e1 <- as.data.frame(t(rbind(lshuff,p)))
colnames(e1) <- c("lshuff", "10", "50", "100", "150", "200")
e11 <- melt(e1, id.vars = "lshuff")
e11$variable <- factor(e11$variable, levels = unique(e11$variable))
e <- error[i]
tfinal <- paste(t_unshuff, as.character(e), sep=" ")
plot_unshuff[[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*sqrt(value*(1-value)/nMC),
ymax=value+2*sqrt(value*(1-value)/nMC)),width=10)+
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("10"="#B675E0","50"="#00A6CA", "100"="#00AA5A", "150"="#AA9000", "200"="#E16A86"),breaks =
levels(e11$variable))+
scale_shape_manual(values=c("10"= 18,"50"=16, "100"= 8, "150"= 17, "200"= 15),breaks = levels(e11$variable))
}
library(patchwork)
plot[[1]] + plot[[2]] + plot_unshuff[[1]] + plot_unshuff[[2]]+ plot_layout(ncol = 2, guides ="collect")
