add pkgdown and unit tests

.Rbuildignore

@@ -2,3 +2,7 @@
 ^\.Rproj\.user$
 
 ^Readme.Rmd$
+^\.github$
+^_pkgdown\.yml$
+^docs$
+^pkgdown$

.github/workflows/R-CMD-check.yaml

@@ -0,0 +1,49 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+
+name: R-CMD-check
+
+jobs:
+  R-CMD-check:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: macos-latest,   r: 'release'}
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
+          - {os: ubuntu-latest,   r: 'release'}
+          - {os: ubuntu-latest,   r: 'oldrel-1'}
+
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+
+    steps:
+      - uses: actions/checkout@v3
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::rcmdcheck
+          needs: check
+
+      - uses: r-lib/actions/check-r-package@v2
+        with:
+          upload-snapshots: false

.github/workflows/pkgdown.yaml

@@ -0,0 +1,46 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, pkgdown]
+  pull_request:
+    branches: [main]
+  release:
+    types: [published]
+  workflow_dispatch:
+
+name: pkgdown
+
+jobs:
+  pkgdown:
+    runs-on: ubuntu-latest
+    # Only restrict concurrency for non-PR jobs
+    concurrency:
+      group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }}
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    steps:
+      - uses: actions/checkout@v3
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::pkgdown, any::qgraph, local::.
+          needs: website
+
+      - name: Build site
+        run: pkgdown::build_site_github_pages(new_process = FALSE, run_dont_run = TRUE, install = FALSE)
+        shell: Rscript {0}
+
+      - name: Deploy to GitHub pages 🚀
+        if: github.event_name != 'pull_request'
+        uses: JamesIves/github-pages-deploy-action@v4.4.1
+        with:
+          clean: false
+          branch: gh-pages
+          folder: docs

DESCRIPTION

@@ -17,5 +17,9 @@ LazyData: true
 Encoding: UTF-8
 Suggests: 
     knitr,
-    rmarkdown
+    qgraph,
+    rmarkdown,
+    testthat (>= 3.0.0)
 VignetteBuilder: knitr
+Config/testthat/edition: 3
+Config/Needs/website: tidyverse/tidytemplate

R/bgm.R

@@ -87,8 +87,8 @@
 #'  abline(h = 0, lty = 2, col = "gray")
 #'  abline(h = 1, lty = 2, col = "gray")
 #'  abline(h = .5, lty = 2, col = "gray")
-#'  mtext("Posterior Inclusion Probability", side = 1, line = 3, cex = 1.7)
-#'  mtext("Posterior Mode Edge Weight", side = 2, line = 3, cex = 1.7)
+#'  mtext("Posterior Mode Edge Weight", side = 1, line = 3, cex = 1.7)
+#'  mtext("Posterior Inclusion Probability", side = 2, line = 3, cex = 1.7)
 #'  axis(1)
 #'  axis(2, las = 1)
 #'

README.md

@@ -1,3 +1,6 @@
+<!-- badges: start -->
+[![R-CMD-check](https://github.com/MaartenMarsman/bgms/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/MaartenMarsman/bgms/actions/workflows/R-CMD-check.yaml)
+<!-- badges: end -->
 
 # bgms: Bayesian Analysis of Graphical Models
 

_pkgdown.yml

@@ -0,0 +1,5 @@
+url: https://maartenMarsman.github.io/bgms/
+template:
+  package: tidytemplate
+  bootstrap: 5
+

tests/testthat.R

@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/tests.html
+# * https://testthat.r-lib.org/reference/test_package.html#special-files
+
+library(testthat)
+library(bgms)
+
+test_check("bgms")

tests/testthat/test-bgm.R

@@ -0,0 +1,10 @@
+test_that("inclusion probabilities correlate with posterior mode", {
+  data("Wenchuan", package = "bgms")
+  fit <- bgm(x = Wenchuan, iter = 1e2, burnin = 10)
+
+  posterior_modes      <- fit$interactions[lower.tri(fit$interactions)]
+  posterior_incl_probs <- fit$gamma[lower.tri(fit$gamma)]
+
+  testthat::expect_gte(cor(abs(posterior_modes), posterior_incl_probs, method = "spearman"), .9)
+
+})

vignettes/introduction.Rmd

@@ -81,13 +81,13 @@ In summary, the list includes matrices for interactions, gamma, and thresholds,
 
 ## Analysis
 
-```{r message=FALSE, warning=FALSE, cache=TRUE, eval=FALSE}
+```{r message=FALSE, warning=FALSE, cache=TRUE}
 fit <- bgm.em(x = Wenchuan)
 ```
 
 
 Now let's show the estimated edge weights and their inclusion probabilities in a plot. 
-```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7, eval=FALSE}
+```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7}
 par(mar = c(6, 5, 1, 1))
 plot(x = fit$interactions[lower.tri(fit$interactions)], 
     y = fit$gamma[lower.tri(fit$gamma)], ylim = c(0, 1), 
@@ -109,7 +109,7 @@ When using EM edge selection, a median probability model is a network model that
 
 Here is some code to extract and plot the local median probability network:
 
-```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7, eval=FALSE}
+```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7}
 library(qgraph) #For plotting the estimated network
 
 posterior.inclusion <- fit$gamma[lower.tri(fit$gamma)]
@@ -196,23 +196,23 @@ Column averages of these matrices provide the model-averaged posterior means.
 
 ## Analysis
 
-```{r message=FALSE, warning=FALSE, cache=TRUE, eval=FALSE}
+```{r message=FALSE, warning=FALSE, cache=TRUE}
 fit <-  bgm(x = Wenchuan)
 ```
 
 To save time, we ran the algorithm using the default number of iterations, which is 10,000. However, this may not be enough to fully explore the posterior distribution of the network structures and parameters. To obtain reliable and accurate estimates, we recommend increasing the number of iterations to 100,000 or more. 
 
 Let's reproduce the plot from the previous example.
-```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7, eval=FALSE}
+```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7}
 par(mar = c(6, 5, 1, 1))
 plot(x = fit$interactions[lower.tri(fit$interactions)], 
      y = fit$gamma[lower.tri(fit$gamma)], ylim = c(0, 1), 
      xlab = "", ylab = "", axes = FALSE, pch = 21, bg = "gray", cex = 1.3)
 abline(h = 0, lty = 2, col = "gray")
 abline(h = 1, lty = 2, col = "gray")
 abline(h = .5, lty = 2, col = "gray")
-mtext("Posterior inclusion probability", side = 1, line = 3, cex = 1.7)
-mtext("Posterior mean edge weight", side = 2, line = 3, cex = 1.7)
+mtext("Posterior mean edge weight", side = 1, line = 3, cex = 1.7)
+mtext("Posterior inclusion probability", side = 2, line = 3, cex = 1.7)
 axis(1)
 axis(2, las = 1)
 ```
@@ -234,7 +234,7 @@ The greater variety in inclusion probabilities in Example 2 results from the `bg
 
 Using the posterior inclusion probabilities, we can also identify the median probability network. In this network, we include all edges with a posterior inclusion probability greater than `0.5`. We can reuse the code from earlier to create the median probability model.
 
-```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7, eval=FALSE}
+```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7}
 library(qgraph) #For plotting the estimated network
 
 posterior.inclusion <- fit$gamma[lower.tri(fit$gamma)]
@@ -265,7 +265,7 @@ One of the benefits of using a fully Bayesian approach is that it allows us to c
 
 In the current version of `bgms`, it is assumed that the prior inclusion probabilities are equal to `0.5`. To calculate the inclusion Bayes factors, we can simply convert the estimated posterior inclusion probabilities. For easier visualization, it is common to use the natural logarithm of the Bayes factor when plotting.
 
-```{r message=FALSE, warning=FALSE, eval=FALSE}
+```{r message=FALSE, warning=FALSE}
 # Calculate the inclusion BFs
 prior.odds = 1
 posterior.inclusion = fit$gamma[lower.tri(fit$gamma)]
@@ -277,7 +277,7 @@ log.bayesfactor[log.bayesfactor > 5] = 5
 
 Lets plot the relation between the estimated edge weights and the inclusion Bayes factor.
 
-```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7, eval=FALSE}
+```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7}
 par(mar = c(5, 5, 1, 1) + 0.1)
 plot(fit$interactions[lower.tri(fit$interactions)], log.bayesfactor, pch = 21, bg = "#bfbfbf", 
     cex = 1.3, axes = FALSE, xlab = "", ylab = "", ylim = c(-5, 5.5),
@@ -302,15 +302,15 @@ In this example, we use a cut-off value of `10` for the inclusion Bayes factors.
 
 For most purposes, the default output from `bgm` is sufficient, providing us with the posterior means of edge indicators and parameters. However, in some cases, we may want to use the raw samples from the joint posterior distribution. This could be to estimate the posterior distribution of a specific parameter, assess how many network structures fit the given data, or create Bayes factors for hypotheses involving multiple edges. We can obtain the raw samples by setting `save = TRUE`.
 
-```{r message=FALSE, warning=FALSE, cache=TRUE, eval=FALSE}
+```{r message=FALSE, warning=FALSE, cache=TRUE}
 fit <-  bgm(x = Wenchuan, save = TRUE)
 ```
 
 ### Posterior density of edge weight
 
 We can employ the following code to use the posterior samples for plotting the posterior density of a single edge:
 
-```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7, eval=FALSE}
+```{r message=FALSE, warning=FALSE, fig.width= 7, fig.height= 7}
 den = density(fit$interactions[,1], bw = "SJ")
 i = which.min(abs(den$x - mean(fit$interactions[,1])))[1]
 x = den$x[i]
@@ -334,15 +334,15 @@ Note that the estimate is not very smooth. This is because we only used 10,000 s
 
 We can also use the raw samples to count the number of unique structures `bgm` encountered in 10,000 iterations. 
 
-```{r message=FALSE, warning=FALSE, eval=FALSE}
+```{r message=FALSE, warning=FALSE}
 G = 2 * fit$gamma - 1
 S = unique(G)
 nrow(S)
 ```
 
 There are clearly many different network structures that could fit the data. Let's estimate their posterior probabilities.
 
-```{r message=FALSE, warning=FALSE, eval=FALSE}
+```{r message=FALSE, warning=FALSE}
 Ps = vector(length = nrow(S))
 for(r in 1:nrow(S)) {
   s = S[r, ]