knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "man/figures/README-", dev = "png", dpi = 200, fig.align = "center", knitr::opts_chunk$set(comment = NA) ) library(ggplot2) library(BGGM)
The R
package BGGM provides tools for making Bayesian inference in
Gaussian graphical models (GGM). The methods are organized around two general
approaches for Bayesian inference: (1) estimation and (2) hypothesis
testing. The key distinction is that the former focuses on either
the posterior or posterior predictive distribution [@Gelman1996a; see section 5 in @rubin1984bayesianly], whereas the latter focuses on model comparison with the Bayes factor [@Jeffreys1961; @Kass1995].
A Gaussian graphical model captures conditional (in)dependencies among a set of variables. These are pairwise relations (partial correlations) controlling for the effects of all other variables in the model.
The Gaussian graphical model is used across the sciences, including (but not limited to) economics [@millington2020partial], climate science [@zerenner2014gaussian], genetics [@chu2009graphical], and psychology [@rodriguez2020formalizing].
To install the latest release version (2.1.1
) from CRAN use
install.packages("BGGM")
The current developmental version can be installed with
if (!requireNamespace("remotes")) { install.packages("remotes") } remotes::install_github("donaldRwilliams/BGGM")
The methods in BGGM build upon existing algorithms that are well-known in the literature. The central contribution of BGGM is to extend those approaches:
Bayesian estimation with the novel matrix-F prior distribution [@Mulder2018]
Bayesian hypothesis testing with the matrix-F prior distribution [@Williams2019_bf]
Comparing Gaussian graphical models [@Williams2019; @williams2020comparing]
Extending inference beyond the conditional (in)dependence structure [@Williams2019]
Posterior uncertainty intervals for the partial correlations
The computationally intensive tasks are written in c++
via the R
package Rcpp [@eddelbuettel2011rcpp] and the c++
library Armadillo [@sanderson2016armadillo]. The Bayes factors are computed with the R
package BFpack [@mulder2019bfpack]. Furthermore, there are plotting functions
for each method, control variables can be included in the model (e.g., ~ gender
),
and there is support for missing values (see bggm_missing
).
Continuous: The continuous method was described in @Williams2019. Note that this is based on the customary Wishart distribution.
Binary: The binary method builds directly upon @talhouk2012efficient that, in turn, built upon the approaches of @lawrence2008bayesian and @webb2008bayesian (to name a few).
Ordinal: The ordinal methods require sampling thresholds. There are two approach included in BGGM. The customary approach described in @albert1993bayesian (the default) and the 'Cowles' algorithm described in @cowles1996accelerating.
Mixed: The mixed data (a combination of discrete and continuous) method was introduced in @hoff2007extending. This is a semi-parametric copula model (i.e., a copula GGM) based on the ranked likelihood. Note that this can be used for only ordinal data (not restricted to "mixed" data).
The following includes brief examples for some of the methods in BGGM.
An ordinal GGM is estimated with
library(BGGM) library(ggplot2) # data Y <- ptsd[,1:5] + 1 # ordinal fit <- estimate(Y, type = "ordinal", analytic = FALSE)
Notice the + 1
. This is required, because the first category must be 1
when type = "ordinal"
. The partial correlations can the be summarized with
summary(fit) #> BGGM: Bayesian Gaussian Graphical Models #> --- #> Type: ordinal #> Analytic: FALSE #> Formula: #> Posterior Samples: 250 #> Observations (n): #> Nodes (p): 5 #> Relations: 10 #> --- #> Call: #> estimate(Y = Y, type = "ordinal", analytic = FALSE, iter = 250) #> --- #> Estimates: #> Relation Post.mean Post.sd Cred.lb Cred.ub #> B1--B2 0.258 0.079 0.105 0.418 #> B1--B3 0.028 0.086 -0.127 0.189 #> B2--B3 0.517 0.058 0.406 0.616 #> B1--B4 0.356 0.070 0.210 0.486 #> B2--B4 -0.076 0.078 -0.219 0.063 #> B3--B4 0.246 0.077 0.107 0.385 #> B1--B5 0.131 0.080 -0.020 0.279 #> B2--B5 0.127 0.083 -0.040 0.284 #> B3--B5 0.202 0.079 0.063 0.366 #> B4--B5 0.349 0.070 0.209 0.474 #> ---
The returned object can also be plotted, which allows for visualizing the posterior uncertainty interval for each relation. An example is provided below in Posterior uncertainty intervals. The partial correlation matrix is accessed with
pcor_mat(fit)
load(file = "readme_models/pcor_ex.rda") knitr::kable(pcor_ex, row.names = TRUE)
The graph is selected with
E <- select(fit)
and then plotted
# "communities" comm <- substring(colnames(Y), 1, 1) plot(select(fit), groups = comm, edge_magnify = 5, palette = "Pastel1", node_size = 12)
This basic "workflow" can be used with all methods and data types. A more involved network plot is provided below.
There is also an analytic solution that is based on the Wishart distribution. This
simple solution provides competitive performance with "state-of-the-art" methods,
assuming that n (observations) > p (variables). The one
caveat is that it works only for type = "continuous"
(the default).
# analytic fit <- estimate(Y, analytic = TRUE) # network plot plot(select(fit))
This is quite handy when (1) only the conditional dependence structure is of interest and (2) an immediate solution is desirable. An example of (2) is provided in Posterior Predictive Check.
The Bayes factor based methods allow for determining the conditional independence structure (evidence for the null hypothesis).
# now 10 nodes Y <- ptsd[,1:10] # exploratory hypothesis testing fit<- explore(Y) # select E <- select(fit, alternative = "exhaustive")
The option alternative = "exhaustive"
compares three hypotheses: (1) a null relation; (2) a positive
relation; and (3) a negative relation.
summary(E) #> BGGM: Bayesian Gaussian Graphical Models #> --- #> Type: ordinal #> Alternative: exhaustive #> --- #> Call: #> select.explore(object = fit, alternative = "exhaustive") #> --- #> Hypotheses: #> H0: rho = 0 #> H1: rho > 0 #> H2: rho < 0 #> --- #> #> Relation Post.mean Post.sd Pr.H0 Pr.H1 Pr.H2 #> B1--B2 0.263 0.080 0.000 0.999 0.001 #> B1--B3 0.020 0.081 0.710 0.173 0.116 #> B2--B3 0.523 0.073 0.000 1.000 0.000 #> B1--B4 0.362 0.070 0.000 1.000 0.000 #> B2--B4 -0.082 0.068 0.459 0.061 0.480 #> B3--B4 0.252 0.073 0.000 1.000 0.000 #> B1--B5 0.129 0.072 0.120 0.847 0.033 #> B2--B5 0.118 0.078 0.223 0.726 0.051 #> B3--B5 0.213 0.077 0.001 0.996 0.003 #> B4--B5 0.348 0.072 0.000 1.000 0.000
The posterior hypothesis probabilities are provided in the last three columns.
When using plot(E)
, there is a network plot for each hypothesis.
A central contribution of BGGM is confirmatory hypothesis testing of (in)equality constraints [@Hoijtink2011]. By this we are referring to testing expectations, as opposed to feeding the data to, say, estimate
, and seeing what happens to emerge.
In this example, the focus is on suicidal thoughts (PHQ9
) in a comorbidity network. Here is an example set of hypotheses
# data (+ 1) Y <- depression_anxiety_t1 + 1 # example hypotheses hyp <- c("PHQ2--PHQ9 > PHQ1--PHQ9 > 0; PHQ2--PHQ9 = PHQ1--PHQ9 = 0")
There are two hypotheses separated by (;
). The first expresses that the relation PHQ2--PHQ9
("feeling down, depressed, or hopeless" and "suicidal thoughts") is larger than PHQ1--PHQ9
("little interest or pleasure in doing things" and "suicidal thoughts"). In other words, that the partial correlation is larger for PHQ2--PHQ9
. There is an additional constraint to positive values (> 0
) for both relations. The second hypothesis is then a "null" model.
load(file = "readme_models/fit_hyp1.rda") fit <- fit_hyp1
# (try to) confirm fit <- confirm(Y = Y, hypothesis = hyp, type = "ordinal")
The object fit
is then printed
fit #> BGGM: Bayesian Gaussian Graphical Models #> Type: ordinal #> --- #> Posterior Samples: 250 #> Observations (n): 403 #> Variables (p): 16 #> Delta: 15 #> --- #> Call: #> confirm(Y = Y + 1, hypothesis = hyp, type = "ordinal", #> iter = 250) #> --- #> Hypotheses: #> #> H1: PHQ2--PHQ9>PHQ1--PHQ9>0 #> H2: PHQ2--PHQ9=PHQ1--PHQ9=0 #> H3: complement #> --- #> Posterior prob: #> #> p(H1|data) = 0.895 #> p(H2|data) = 0.002 #> p(H3|data) = 0.103 #> --- #> Bayes factor matrix: #> H1 H2 H3 #> H1 1.000 529.910 8.666 #> H2 0.002 1.000 0.016 #> H3 0.115 61.147 1.000 #> --- #> note: equal hypothesis prior probabilities
The posterior hypothesis probability is 0.895
which provides some evidence for the order
constraint. The Bayes factor matrix then divides the posterior probabilities. This provide a measure
of relative support for which hypothesis the data were more likely under.
Finally, the results can be plotted
plot(fit) + scale_fill_brewer(palette = "Set2", name = "Posterior Prob") + ggtitle("Confirmatory: Comorbidity Network")
This demonstrates that all the plot()
functions in BGGM return ggplot
objects that can be further customized. Note that BGGM is not focused on making publication ready plots. Typically the bare minimum is provided that can then be honed in.
This method compares groups by computing the difference for each relation in the model. In other words, there are pairwise contrasts for each partial correlation, resulting in a posterior distribution for each difference.
In all examples in this section, personality networks are compared for males and females.
# data Y <- bfi # males Ymales <- subset(Y, gender == 1, select = -c(gender, education)) # females Yfemales <- subset(Y, gender == 2, select = -c(gender, education))
Fit the model
fit <- ggm_compare_estimate(Ymales, Yfemales)
Then plot the results, in this case the posterior distribution for each difference
# plot summary plot(summary(fit))
Note that it is also possible to use select
for the object fit
and then plot the results. This
produces a network plot including the selected differences.
Furthermore, it is also possible to plot the partial correlations (not the differences). This
is accomplished by using plot
with the summary computed from an estimate
object
(see above).
The predictive check method uses Jensen-Shannon divergence (i.e., symmetric Kullback-Leibler divergence Wikipedia) and the sum of squared error (for the partial correlation matrices) to compare groups [@williams2020comparing].
The following compares the groups
fit <- ggm_compare_ppc(Ymales, Yfemales)
Then print the summary output with
fit #> BGGM: Bayesian Gaussian Graphical Models #> --- #> Test: Global Predictive Check #> Posterior Samples: 500 #> Group 1: 805 #> Group 2: 1631 #> Nodes: 25 #> Relations: 300 #> --- #> Call: #> ggm_compare_ppc(Ymales, Yfemales, iter = 500) #> --- #> Symmetric KL divergence (JSD): #> #> contrast JSD.obs p_value #> Yg1 vs Yg2 0.442 0 #> --- #> #> Sum of Squared Error: #> #> contrast SSE.obs p.value #> Yg1 vs Yg2 0.759 0 #> --- #> note: #> JSD is Jensen-Shannon divergence
In this case, there seems to be decisive evidence that the networks are different (as indicated by the posterior predictive p-value). The predictive distribution can also be plotted
plot(fit, critical = 0.05)$plot_jsd
where the red region is the "critical" area and the black point is the observed KL divergence for the networks. This again shows that the "distance" between the networks is much more than expected, assuming that the groups were actually the same.
This next example is a new feature in BGGM (2.0.0
), that allows for comparing GGMs any way the user wants. All that is required is to (1) decide on a test-statistic and (2) write a custom function.
Here is an example using Hamming distance (Wikipedia), which is essentially the squared error between adjacency matrices (a test for different structures).
First define the custom function
f <- function(Yg1, Yg2){ # remove NA x <- na.omit(Yg1) y <- na.omit(Yg2) # nodes p <- ncol(x) # identity matrix I_p <- diag(p) # estimate graphs fit1 <- estimate(x, analytic = TRUE) fit2 <- estimate(y, analytic = TRUE) # select graphs sel1 <- select(fit1) sel2 <- select(fit2) # Hamming distance sum((sel1$adj[upper.tri(I_p)] - sel2$adj[upper.tri(I_p)])^2) }
Note that (1) analytic = TRUE
is being used, which is needed in this case because two graphs are
estimated for each iteration (or draw from the posterior predictive distribution) and (2) f
requires two datasets as the input and returns a single number (the chosen test-statistic). The next step is to compute the observed Hamming distance
load(file = "readme_models/obs.rda")
# observed difference obs <- f(Ymales, Yfemales)
then compare the groups
fit <- ggm_compare_ppc(Ymales, Yfemales, FUN = f, custom_obs = obs) # print fit #> BGGM: Bayesian Gaussian Graphical Models #> --- #> Test: Global Predictive Check #> Posterior Samples: 250 #> Group 1: 805 #> Group 2: 1631 #> Nodes: 25 #> Relations: 300 #> --- #> Call: #> ggm_compare_ppc(Ymales, Yfemales, iter = 250, FUN = f, custom_obs = obs) #> --- #> Custom: #> #> contrast custom.obs p.value #> Yg1 vs Yg2 75 0.576 #> ---
In this case, the p-value does not indicate that the groups are different for this test-statistic. This may seem contradictory to the previous results, but it is important to note that Hamming distance asks a much different question related to the adjacency matrices (no other information, such as edge weights, is considered).
The Bayes factor based methods allow for determining the conditional independence structure (evidence for the null hypothesis), in this case for group equality.
Fit the model
fit <- ggm_compare_explore(Ymales, Yfemales)
Then plot the results
plot(summary(fit)) + theme_bw() + theme(panel.grid = element_blank(), axis.text.y = element_blank())
Here the posterior probability for a difference is visualized for each relation in the GGM. Note
that it is also possible to use select
for the object fit
and then plot the results. This
produces a network plot including the selected differences, in addition to a plot depicting the
relations for which there was evidence for the null hypothesis.
A central contribution of BGGM is confirmatory hypothesis testing of (in)equality constraints [@Hoijtink2011],
in this case for comparing groups. By this we are referring to testing expectations, as opposed to feeding the data to, say, estimate
, and seeing what happens to emerge.
In this example, the focus is on agreeableness in a personality network. Here is a set of hypotheses
hyp <- c("g1_A2--A4 > g2_A2--A4 > 0 & g1_A4--A5 > g2_A4--A5 > 0; g1_A4--A5 = g2_A4--A5 = 0 & g1_A2--A4 = g2_A2--A4 = 0")
where the variables are A2
("inquire about others' well being"), A4
("love children"),
and A5
("make people feel at ease"). The first hypothesis states that the conditionally
dependent effects are larger for female than males (note the &
), with the additional
constraint to positive values, whereas the second hypothesis is a "null" model.
The hypothesis is tested with the following
fit <- ggm_compare_confirm(Yfemales, Ymales, hypothesis = hyp) # print fit #> BGGM: Bayesian Gaussian Graphical Models #> Type: continuous #> --- #> Posterior Samples: 500 #> Group 1: 1631 #> Group 2: 805 #> Variables (p): 25 #> Relations: 300 #> Delta: 15 #> --- #> Call: #> ggm_compare_confirm(Yfemales, Ymales, hypothesis = hyp, iter = 500) #> --- #> Hypotheses: #> #> H1: g1_A2--A4>g2_A2--A4>0&g1_A4--A5>g2_A4--A5>0 #> H2: g1_A4--A5=g2_A4--A5=0&g1_A2--A4=g2_A2--A4=0 #> H3: complement #> --- #> Posterior prob: #> #> p(H1|data) = 0.989 #> p(H2|data) = 0 #> p(H3|data) = 0.011 #> --- #> Bayes factor matrix: #> H1 H2 H3 #> H1 1.000 1.180798e+14 92.115 #> H2 0.000 1.000000e+00 0.000 #> H3 0.011 1.281873e+12 1.000 #> --- #> note: equal hypothesis prior probabilities
The posterior hypothesis probability is 0.989 which provides strong evidence for the hypothesis that predicted these "agreeableness" relations would be larger in females than in males. This can also be plotted, as in Confirmatory (one group). See @rodriguez2020formalizing for a full treatment of confirmatory testing in substantive applications.
In this example, predictability is computed for each node in the network [see here for rationale @haslbeck2018well]. Currently BGGM computes Bayesian variance explained for all data types [@gelman_r2_2019].
The following computes predictability for binary data
# binary Y <- women_math # fit model fit <- estimate(Y, type = "binary") # compute r2 r2 <- predictability(fit, iter = 500) # plot plot(r2, type = "ridgeline")
See Partial Correlation Differences
A new feature to BGGM allows for computing user defined network statistics, given a partial correlation or weighted adjacency matrix.
Here is an example for bridge centrality [@jones2019bridge]. The first step is to define the function
# need this package library(networktools) # custom function f <- function(x, ...){ bridge(x, ...)$`Bridge Strength` }
Note that x
takes the matrix and f
can return either a single number or a number for each node. The next step is to fit the model and compute the network statistic
# data Y <- ptsd # clusters communities <- substring(colnames(Y), 1, 1) # estimate the model fit <- estimate(Y) # bridge strength net_stat <- roll_your_own(fit, FUN = f, select = TRUE, communities = communities)
The function f
is provided to FUN
and communities
is passed to brigde
(inside of f
) via ...
. The results can be printed
# print net_stat #> BGGM: Bayesian Gaussian Graphical Models #> --- #> Network Stats: Roll Your Own #> Posterior Samples: 100 #> --- #> Estimates: #> #> Node Post.mean Post.sd Cred.lb Cred.ub #> 1 0.340 0.097 0.166 0.546 #> 2 0.319 0.100 0.176 0.513 #> 3 0.000 0.000 0.000 0.000 #> 4 0.337 0.086 0.189 0.489 #> 5 0.559 0.133 0.332 0.791 #> 6 0.188 0.073 0.029 0.320 #> 7 0.505 0.138 0.241 0.781 #> 8 0.153 0.070 0.022 0.286 #> 9 0.175 0.063 0.041 0.281 #> 10 0.000 0.000 0.000 0.000 #> 11 0.365 0.107 0.178 0.627 #> 12 0.479 0.093 0.280 0.637 #> 13 0.155 0.074 0.022 0.301 #> 14 0.000 0.000 0.000 0.000 #> 15 0.374 0.097 0.175 0.550 #> 16 0.174 0.065 0.034 0.295 #> 17 0.000 0.000 0.000 0.000 #> 18 0.491 0.132 0.238 0.745 #> 19 0.613 0.113 0.408 0.825 #> 20 0.144 0.066 0.038 0.289 #> ---
And then plotted
plot(net_stat)
There are additional examples in the documentation.
Here is an example of a more involved network plot. In this case, the graph is estimated with a semi-parametric copula (type = "mixed"
), where two control variables are included in the model.
load(file = "readme_models/fit_plt_ex.rda") fit <- fit_plt_ex # select graph E <- select(fit)
# personality (includes gender and education) Y <- bfi # fit copula GGM fit <- estimate(Y, type = "mixed") # select graph E <- select(fit)
The graph is then plotted
# extract communities comm <- substring(colnames(Y), 1, 1) # plot plot(E, # enlarge edges edge_magnify = 5, # cluster nodes groups = comm, # change layout layout = "circle")$plt + # plot title ggtitle("Semi-Parametric Copula") + # add custom labels scale_color_brewer(breaks = c("A", "C", "E", "N", "O", "e", "g"), labels = c("A", "C", "E", "N", "O", "Education", "Gender"), palette = "Set2")
Note that layout
can be changed to any option provided in the R
package sna [@sna].
The primary focus of BGGM is Gaussian graphical modeling (the inverse covariance matrix). The residue is a suite of useful methods not explicitly for GGMs. For example,
Bivariate correlations for binary
(tetrachoric), ordinal
(polychoric), mixed
(rank based),
and continuous
(Pearson's) data.
Here is an example for computing tetrachoric correlations:
load(file = "readme_models/binary_cors.rda")
# binary data Y <- women_math[1:500,] cors <- zero_order_cors(Y, type = "binary", iter = 250) cors$R
row.names(cors$R_mean) <- 1:6 knitr::kable(round(cors$R_mean, 3), col.names = c("1", "2", "3","4","5", "6"), row.names = TRUE)
The object cors
also includes the sampled correlation matrices (in this case 250) in an array.
Multivariate regression for binary (probit), ordinal (probit), mixed (rank likelihood), and continuous data.
Here is an example for a multivariate probit model with an ordinal outcome, where
E5
("take charge") and N5
("panic easily") are predicted by gender
and education
:
load("readme_models/mv_probit.rda")
# personality data Y <- bfi # variables Y <- subset(Y, select = c("E5", "N5", "gender", "education")) mv_probit <- estimate(Y, formula = ~ gender + as.factor(education), type = "ordinal")
Note that BGGM does not use the customary model.matrix
formulation. This is for good reason, as
each variable in the GGM does not need to be written out. Here we effectively "tricked" BGGM to
fit a multivariate probit model (each variable included in formula
is removed from Y
).
regression_summary(mv_probit)
This basic idea can also be used to fit regression models with a single outcome.
All of the data types (besides continuous) model latent data. That is, unobserved data that is assumed to be Gaussian distributed. For example, a tetrachoric correlation (binary data) is a special case of a polychoric correlation (ordinal data). Both relations are between "theorized normally distributed continuous latent variables Wikepedia. In both instances, the corresponding partial correlation between observed variables is conditioned on the remaining variables in the latent space. This implies that interpretation is similar to continuous data, but with respect to latent variables. We refer interested users to [see page 2364, section 2.2, in @webb2008bayesian].
BGGM was built specifically for social-behavioral scientists. Of course, the methods can be used by all researchers. However, there is currently not support for high-dimensional data (i.e., more variables than observations) that are common place in, say, the genetics literature. These data are rare in the social-behavioral sciences. In the future, support for high-dimensional data may be added to BGGM.
Bug reports and feature requests can be made by opening an issue on Github. To contribute towards the development of BGGM, you can start a branch with a pull request and we can discuss the proposed changes there.
BGGM is the only R
package to implement all of these algorithms and methods. The mixed
data approach
is also implemented in the package sbgcop [base R
, @hoff2007extending]. The R
package BDgraph implements a Gaussian copula graphical model in c++
[@mohammadi2015bdgraph], but not the binary or ordinal approaches. Furthermore, BGGM is the only package for confirmatory testing and comparing graphical models with the methods described in @williams2020comparing.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.