The goal of GGMnonreg is to estimate non-regularized graphical models. Note that the title is a bit of a misnomer, in that Ising and mixed graphical models are also supported.
Graphical modeling is quite common in fields with wide data, that is, when there are more variables than observations. Accordingly, many regularization-based approaches have been developed for those kinds of data. There are key drawbacks of regularization when the goal is inference, including, but not limited to, the fact that obtaining a valid measure of parameter uncertainty is very (very) difficult.
More recently, graphical modeling has emerged in psychology (Epskamp et al. 2018), where the data is typically long or low-dimensional [p < n; Donald R. Williams et al. (2019); Donald R. Williams and Rast (2019)]. The primary purpose of GGMnonreg is to provide methods specifically for low-dimensional data (e.g., those common to psychopathology networks).
The following are also included
To install the latest release version (1.1.0) from CRAN use
install.packages("GGMnonreg")
You can install the development version from GitHub with:
# install.packages("devtools")
devtools::install_github("donaldRwilliams/GGMnonreg")
An Ising model is fitted with the following
library(GGMnonreg)
# make binary
Y <- ifelse(ptsd[,1:5] == 0, 0, 1)
# fit model
fit <- ising_search(Y, IC = "BIC",
progress = FALSE)
fit
#> 1 2 3 4 5
#> 1 0.000000 1.439583 0.000000 1.273379 0.000000
#> 2 1.439583 0.000000 1.616511 0.000000 1.182281
#> 3 0.000000 1.616511 0.000000 1.716747 1.077322
#> 4 1.273379 0.000000 1.716747 0.000000 1.662550
#> 5 0.000000 1.182281 1.077322 1.662550 0.000000
Note the same code, more or less, is also used for GGMs and mixed graphical models.
It is common to compute predictability, or variance explained, for each node in the network. An advantage of GGMnonreg is that a measure of uncertainty is also provided.
# data
Y <- na.omit(bfi[,1:5])
# fit model
fit <- ggm_inference(Y, boot = FALSE)
# predictability
predictability(fit)
#> Estimate Est.Error Ci.lb Ci.ub
#> 1 0.13 0.01 0.11 0.15
#> 2 0.33 0.01 0.30 0.36
#> 3 0.38 0.01 0.35 0.41
#> 4 0.18 0.01 0.15 0.20
#> 5 0.29 0.01 0.26 0.32
Confidence intervals for each relation are obtained with
# data
Y <- na.omit(bfi[,1:5])
# fit model
fit <- ggm_inference(Y, boot = TRUE,
method = "spearman",
B = 100, progress = FALSE)
#> | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |========= | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============== | 19% | |============== | 20% | |=============== | 21% | |================ | 22% | |================ | 23% | |================= | 24% | |================== | 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% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
confint(fit)
#> 2.5% 97.5%
#> [1,] -0.29205233 -0.212951642
#> [2,] -0.14926375 -0.075919715
#> [3,] 0.25508466 0.326055944
#> [4,] -0.03778214 0.032287534
#> [5,] 0.12219960 0.204592760
#> [6,] 0.13387070 0.206412243
#> [7,] -0.07641808 0.006226351
#> [8,] 0.10997005 0.195192621
#> [9,] 0.34388225 0.418927425
#> [10,] 0.08130533 0.153621543
These can then be plotted with, say, ggplot2 (left to the user).
When mining data, or performing an automatic search, it is difficult to make inference on the network parameters (e.g., confidence are not easily computed). To summarize data mining, GGMnonreg provides edge inclusion “probabilities” (proportion bootstrap samples for which each relation was detected).
# data
Y <- na.omit(bfi[,1:5])
# fit model
fit <- eip(Y, method = "spearman",
B = 100, progress = FALSE)
fit
#> eip
#> 1 1.00
#> 2 1.00
#> 3 1.00
#> 4 0.06
#> 5 1.00
#> 6 1.00
#> 7 0.31
#> 8 1.00
#> 9 1.00
#> 10 1.00
Note in all cases, the provided estimates correspond to the upper-triangular elements of the network.
GGMnonreg allows for computing expected network replicability (ENR), i.e., the number of effects that will be detected in any number of replications. This is an analytic solution.
The first step is defining a true network
# first make the true network
cors <- cor(GGMnonreg::ptsd)
# inverse
inv <- solve(cors)
# partials
pcors <- -cov2cor(inv)
# set values to zero
pcors <- ifelse(abs(pcors) < 0.05, 0, pcors)
Then obtain ENR
fit_enr <- enr(net = pcors, n = 500, replications = 2)
fit_enr
#> Average Replicability: 0.45
#> Average Number of Edges: 53 (SD = 3.7)
#>
#> ----
#>
#> Cumulative Probability:
#>
#> prop.edges edges Pr(R > prop.edges)
#> 0.0 0 1.00
#> 0.1 12 1.00
#> 0.2 23 1.00
#> 0.3 35 1.00
#> 0.4 47 0.91
#> 0.5 58 0.05
#> 0.6 70 0.00
#> 0.7 82 0.00
#> 0.8 94 0.00
#> 0.9 105 0.00
#> ----
#> Pr(R > prop.edges):
#> probability of replicating more than the
#> correpsonding proportion (and number) of edges
Note this is inherently frequentist. As such, over the long run, 45 % of the edges will be replicated on average. Then we can further infer that, in hypothetical replication attempts, more than half of the edges will be replicated only 5 % of the time.
ENR can also be plotted
plot_enr(fit_enr)
Here is the basic idea of ENR
# 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 1000 to, say, 5,000
for(i in 1:1000){
# 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)
}
Notice that replication of two networks is being assessed over the long run. In other words, if we draw two random samples, what is the expected replicability.
Compare analytic to simulation
# 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))
)
#> simulation analytic
#> 1 1.000 1.000
#> 2 1.000 1.000
#> 3 1.000 1.000
#> 4 1.000 1.000
#> 5 0.897 0.912
#> 6 0.055 0.052
#> 7 0.000 0.000
#> 8 0.000 0.000
#> 9 0.000 0.000
#> 10 0.000 0.000
# average replicability (simulation)
mean(R / length(index))
#> [1] 0.4482051
# average replicability (analytic)
fit_enr$ave_pwr
#> [1] 0.4485122
ENR works with any correlation, assuming there is an estimate of the standard error.
# data
Y <- ptsd
# estimate graph
fit <- ggm_inference(Y, boot = FALSE)
# get info for plotting
plot(fit, edge_magnify = 5)
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
Bug reports and feature requests can be made by opening an issue on Github. To contribute towards the development of GGMnonreg, you can start a branch with a pull request and we can discuss the proposed changes there.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.