enr {GGMnonreg}

Expected Network Replicability

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

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 Williams (2020).

References

Williams DR (2020). “Learning to live with sampling variability: Expected replicability in partial correlation networks.” PsyArXiv. doi: 10.31234/osf.io/fb4sa, https://doi.org/10.31234/osf.io/fb4sa.

Examples

# 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