enr: Expected Network Replicability

In GGMnonreg: Non-Regularized Gaussian Graphical Models

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 \insertCitewilliams2020learning;textualGGMnonreg.

References

\insertAllCited

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

GGMnonreg documentation built on April 8, 2021, 5:06 p.m.