enr: Expected Network Replicability

Description Usage Arguments Value Note References Examples

View source: R/enr.R

Description

Compute expected network replicability in any number of replication attempts. This works for any kind of correlation, assuming it is possible to obtain the standard error (analytically or with a bootstrap).

Usage

1
enr(net, n, alpha = 0.05, replications = 2, type = "pearson")

Arguments

net

True network of dimensions p by p.

n

Integer. The samples size, assumed equal in the replication attempts.

alpha

The desired significance level (defaults to 0.05). Note that 1 - alpha corresponds to specificity.

replications

Integer. The desired number of replications.

type

Character string. Which type of correlation coefficients to be computed. Options include "pearson" (default) and "spearman".

Value

An list of class enr, including a variety of information used by other functions (e.g., to plot the results).

Note

This method was introduced in \insertCitewilliams2020learning;textualGGMnonreg.

References

\insertAllCited

Examples

 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
# correlations
cors <- cor(GGMnonreg::ptsd)

# inverse
inv <- solve(cors)

# partials
pcors <-  -cov2cor(inv)

# set values to zero
pcors <- ifelse(abs(pcors) < 0.05, 0, pcors)

fit_enr <- enr(net = pcors, n = 500, replications = 2)


# intuition for the method

# location of edges
index <- which(pcors[upper.tri(diag(20))] != 0)

# convert network into correlation matrix
diag(pcors) <- 1
cors_new <- corpcor::pcor2cor(pcors)

# replicated edges
R <- NA

# increase 100 to, say, 5,000
for(i in 1:100){

  # two replications
  Y1 <- MASS::mvrnorm(500, rep(0, 20), cors_new)
  Y2 <- MASS::mvrnorm(500, rep(0, 20), cors_new)

  # estimate network 1
  fit1 <- ggm_inference(Y1, boot = FALSE)

  # estimate network 2
  fit2 <- ggm_inference(Y2, boot = FALSE)

  # number of replicated edges (detected in both networks)
  R[i] <- sum(
    rowSums(
      cbind(fit1$adj[upper.tri(diag(20))][index],
            fit2$adj[upper.tri(diag(20))][index])
    ) == 2)
}


# combine simulation and analytic
cbind.data.frame(
  data.frame(simulation = sapply(seq(0, 0.9, 0.1), function(x) {
    mean(R > round(length(index) * x) )
  })),
  data.frame(analytic = round(fit_enr$cdf, 3))
)

# average replicability (simulation)
mean(R / length(index))

# average replicability (analytic)
fit_enr$ave_pwr

donaldRwilliams/GGMnonreg documentation built on May 13, 2021, 11:57 a.m.