-
Notifications
You must be signed in to change notification settings - Fork 2
/
MASE-hyptest-example.R
96 lines (74 loc) · 2.49 KB
/
MASE-hyptest-example.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
# Load packages and MASE source code
source("R/loadAll.R")
# Generate a pair of mixed-membership SBM graphs with same community structure
m = 2 # sample size
n = 300 # number of vertices
d=3 # common subspace dimension
p = 0.4 # within community edge probability
q = 0.1 # between community edge probability
library(MCMCpack)
# Sample mixed-memberships
Z <- rdirichlet(n, alpha = rep(0.1, d))
b = 500 #bootstrap samples
numsims = 100
# Vector of differences between graphs
epsilons = seq(0, 0.1, 0.01)
p.values.list <- list()
i = 1
# Equal distribution on the two graphs ----------------------------------------------
epsilon = 0
B1 = diag(3)*(p-q) + q
B2 = B1
B2[1, 1] = B2[1, 1] + epsilon
# Edge probabilities of the graphs
P1 = Z %*% B1 %*% t(Z)
P2 = Z %*% B2 %*% t(Z)
Ps.true = list(P1, P2)
# Obtain COSIE parameters from true matrices
mase.true = mase(Ps.true, d)
V = mase.true$V
R1 = mase.true$R[[1]]
R2 = mase.true$R[[2]]
# Sample two graphs
set.seed(1989)
A1 = sample_from_P(P1)
A2 = sample_from_P(P2)
# MASE algorithm
mase.embed = mase(list(A1, A2), d = 3)
mase.embed$R
MASE.test.stat = norm(mase.embed$R[[1]] - mase.embed$R[[2]], "F")
# Calculate null distribution of the test
mase.tstat.null.true1 <- mase.test.statistic.distr(V = V, R = R1, B = b)
mase.tstat.null.true2 <- mase.test.statistic.distr(V = V, R = R2, B = b)
# Obtain conservative p-value
mase.pval.equal.distr = max(p.values(MASE.test.stat, mase.tstat.null.true1),
p.values(MASE.test.stat, mase.tstat.null.true2))
print(mase.pval.equal.distr)
# Different distribution on the two graphs -------------------------------------------
epsilon = 0.05
B1 = diag(3)*(p-q) + q
B2 = B1
B2[1, 1] = B2[1, 1] + epsilon
# Edge probabilities of the graphs
P1 = Z %*% B1 %*% t(Z)
P2 = Z %*% B2 %*% t(Z)
Ps.true = list(P1, P2)
# Obtain COSIE parameters from true matrices
mase.true = mase(Ps.true, d)
V = mase.true$V
R1 = mase.true$R[[1]]
R2 = mase.true$R[[2]]
# Sample two graphs
A1 = sample_from_P(P1)
A2 = sample_from_P(P2)
# MASE algorithm
mase.embed = mase(list(A1, A2), d = 3)
mase.embed$R
MASE.test.stat = norm(mase.embed$R[[1]] - mase.embed$R[[2]], "F")
# Calculate null distribution of the test
mase.tstat.null.true1 <- mase.test.statistic.distr(V = V, R = R1, B = b)
mase.tstat.null.true2 <- mase.test.statistic.distr(V = V, R = R2, B = b)
# Obtain conservative p-value
mase.pval.diff.distr = max(p.values(MASE.test.stat, mase.tstat.null.true1),
p.values(MASE.test.stat, mase.tstat.null.true2))
print(mase.pval.diff.distr)