diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index 9b069bd..33e399a 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -1,5 +1,3 @@
-# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag.
-# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions
on:
push:
branches:
@@ -24,64 +22,84 @@ jobs:
fail-fast: false
matrix:
config:
+ # - {os: macOS-latest, r: 'release'}
- {os: windows-latest, r: 'release'}
- - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
- - {os: ubuntu-20.04, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest", http-user-agent: "R/4.1.0 (ubuntu-20.04) R (4.1.0 x86_64-pc-linux-gnu x86_64 linux-gnu) on GitHub Actions" }
- - {os: macOS-latest, r: 'release'}
+ - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'}
+ - {os: ubuntu-latest, r: 'release'}
+ # - {os: ubuntu-latest, r: 'oldrel-1'}
env:
- R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
RSPM: ${{ matrix.config.rspm }}
GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+ R_REMOTES_NO_ERRORS_FROM_WARNINGS: false
+ NOT_CRAN: true
steps:
- uses: actions/checkout@v2
- - uses: r-lib/actions/setup-r@v1
+ - uses: r-lib/actions/setup-r@master
with:
r-version: ${{ matrix.config.r }}
- - uses: r-lib/actions/setup-pandoc@v1
+ - uses: r-lib/actions/setup-pandoc@master
- name: Query dependencies
run: |
install.packages('remotes')
saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
- writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
shell: Rscript {0}
- - name: Restore R package cache
- uses: actions/cache@v2
+ - name: Cache R packages
+ if: runner.os != 'Windows'
+ uses: actions/cache@v1
with:
path: ${{ env.R_LIBS_USER }}
- key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
- restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
-
- - name: Install system dependencies
- if: runner.os == 'Linux'
- run: |
- while read -r cmd
- do
- eval sudo $cmd
- done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))')
+ key: ${{ runner.os }}-r-${{ matrix.config.r }}-3-${{ hashFiles('.github/depends.Rds') }}
+ restore-keys: ${{ runner.os }}-r-${{ matrix.config.r }}-3-
- - name: Install JAGS (windows-latest)
+ - name: Set path for RTools 4.0
+ if: runner.os == 'Windows'
+ run: echo "C:/rtools40/usr/bin;C:/rtools40/mingw64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
+
+ - name: Install most Windows system dependencies
if: runner.os == 'Windows'
run: |
- curl.exe -o wjags.exe --url https://deac-fra.dl.sourceforge.net/project/mcmc-jags/JAGS/4.x/Windows/JAGS-4.2.0-Rtools33.exe
- wjags.exe /S
- del wjags.exe
+ pacman -Syu mingw-w64-x86_64-make --noconfirm
+ g++ --version
+ Get-Command g++ | Select-Object -ExpandProperty Definition
+ mingw32-make --version
+ Get-Command mingw32-make | Select-Object -ExpandProperty Definition
+ (New-Object System.Net.WebClient).DownloadFile('https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/JAGS-4.3.1.exe', 'C:\JAGS-4.3.1.exe')
+ shell: powershell
+
+ - name: Install JAGS on Windows
+ if: runner.os == 'Windows'
+ run: C:\JAGS-4.3.1.exe /S
shell: cmd
- - name: Install JAGS (macOS-latest)
+ - name: Install Linux system dependencies
+ if: runner.os == 'Linux'
+ env:
+ RHUB_PLATFORM: linux-x86_64-ubuntu-gcc
+ run: |
+ sudo apt-get install -y libglpk-dev
+ Rscript -e "remotes::install_github('r-hub/sysreqs')"
+ sysreqs=$(Rscript -e "cat(sysreqs::sysreq_commands('DESCRIPTION'))")
+ sudo -s eval "$sysreqs"
+ sudo -s apt-get install jags
+ sudo -s apt-get install libv8-dev
+ sudo -s apt-get install libnode-dev
+ sudo -s apt-get install libcurl4-openssl-dev
+
+ - name: Install Mac system dependencies
if: runner.os == 'macOS'
- run : |
- rm '/usr/local/bin/gfortran'
- brew install automake jags
-
+ run: |
+ brew install jags
+ brew install libgit2
+
- name: Install dependencies
run: |
- remotes::install_deps(dependencies = TRUE)
+ remotes::install_deps(dependencies = TRUE, type = ifelse(tolower(Sys.info()["sysname"]) == "linux", "source", "binary"))
remotes::install_cran("rcmdcheck")
shell: Rscript {0}
@@ -90,7 +108,6 @@ jobs:
env:
_R_CHECK_CRAN_INCOMING_REMOTE_: false
run: |
- options(crayon.enabled = TRUE)
rcmdcheck::rcmdcheck(
args = c("--no-manual", "--as-cran", "--ignore-vignettes"),
build_args = c("--no-build-vignettes"),
@@ -102,17 +119,15 @@ jobs:
if: runner.os == 'Windows'
env:
_R_CHECK_CRAN_INCOMING_REMOTE_: false
- JAGS_ROOT: "/c/progra~1/JAGS/JAGS-4.2.0"
+ JAGS_ROOT: "/c/progra~1/JAGS/JAGS-4.3.1"
JAGS_MAJOR_VERSION: 4
run: |
- options(crayon.enabled = TRUE)
rcmdcheck::rcmdcheck(
args = c("--no-manual", "--as-cran", "--ignore-vignettes"),
build_args = c("--no-build-vignettes"),
error_on = "warning",
check_dir = "check")
- shell: Rscript {0}
-
+ shell: Rscript {0}"
- name: Upload check results
if: failure()
diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index 6406aee..b9ee825 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -20,49 +20,85 @@ jobs:
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:
+ RSPM: ${{ matrix.config.rspm }}
GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+ R_REMOTES_NO_ERRORS_FROM_WARNINGS: false
+ NOT_CRAN: true
+
steps:
- uses: actions/checkout@v2
- - uses: r-lib/actions/setup-r@v1
-
- - uses: r-lib/actions/setup-pandoc@v1
-
- - name: Install JAGS (windows-latest)
- if: runner.os == 'Windows'
- run: |
- curl.exe -o wjags.exe --url https://deac-fra.dl.sourceforge.net/project/mcmc-jags/JAGS/4.x/Windows/JAGS-4.2.0-Rtools33.exe
- wjags.exe /S
- del wjags.exe
- shell: cmd
+ - uses: r-lib/actions/setup-r@master
+ with:
+ r-version: ${{ matrix.config.r }}
- - name: Install JAGS (macOS-latest)
- if: runner.os == 'macOS'
- run : |
- rm '/usr/local/bin/gfortran'
- brew install jags
+ - uses: r-lib/actions/setup-pandoc@master
- name: Query dependencies
run: |
install.packages('remotes')
saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
- writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
shell: Rscript {0}
- - name: Restore R package cache
- uses: actions/cache@v2
+ - name: Cache R packages
+ if: runner.os != 'Windows'
+ uses: actions/cache@v1
with:
path: ${{ env.R_LIBS_USER }}
- key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
- restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
+ key: ${{ runner.os }}-r-${{ matrix.config.r }}-3-${{ hashFiles('.github/depends.Rds') }}
+ restore-keys: ${{ runner.os }}-r-${{ matrix.config.r }}-3-
+ - name: Set path for RTools 4.0
+ if: runner.os == 'Windows'
+ run: echo "C:/rtools40/usr/bin;C:/rtools40/mingw64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
+
+ - name: Install most Windows system dependencies
+ if: runner.os == 'Windows'
+ run: |
+ pacman -Syu mingw-w64-x86_64-make --noconfirm
+ g++ --version
+ Get-Command g++ | Select-Object -ExpandProperty Definition
+ mingw32-make --version
+ Get-Command mingw32-make | Select-Object -ExpandProperty Definition
+ (New-Object System.Net.WebClient).DownloadFile('https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/JAGS-4.3.1.exe', 'C:\JAGS-4.3.1.exe')
+ shell: powershell
+
+ - name: Install JAGS on Windows
+ if: runner.os == 'Windows'
+ run: C:\JAGS-4.3.1.exe /S
+ shell: cmd
+
+ - name: Install Linux system dependencies
+ if: runner.os == 'Linux'
+ env:
+ RHUB_PLATFORM: linux-x86_64-ubuntu-gcc
+ run: |
+ sudo apt-get install -y libglpk-dev
+ Rscript -e "remotes::install_github('r-hub/sysreqs')"
+ sysreqs=$(Rscript -e "cat(sysreqs::sysreq_commands('DESCRIPTION'))")
+ sudo -s eval "$sysreqs"
+ sudo -s apt-get install jags
+ sudo -s apt-get install libv8-dev
+ sudo -s apt-get install libnode-dev
+ sudo -s apt-get install libcurl4-openssl-dev
+
+ - name: Install Mac system dependencies
+ if: runner.os == 'macOS'
+ run: |
+ brew install jags
+ brew install libgit2
+
- name: Install dependencies
run: |
- install.packages(c("remotes"))
- remotes::install_deps(dependencies = TRUE)
+ remotes::install_deps(dependencies = TRUE, type = ifelse(tolower(Sys.info()["sysname"]) == "linux", "source", "binary"))
+ remotes::install_cran("rcmdcheck")
remotes::install_cran("covr")
shell: Rscript {0}
diff --git a/DESCRIPTION b/DESCRIPTION
index 6172e89..982074d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: BayesTools
Title: Tools for Bayesian Analyses
-Version: 0.2.10
+Version: 0.2.11
Description: Provides tools for conducting Bayesian analyses. The package contains
functions for creating a wide range of prior distribution objects, mixing posterior
samples from 'JAGS' and 'Stan' models, plotting posterior distributions, and etc...
diff --git a/NAMESPACE b/NAMESPACE
index b732e72..f09f264 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -33,15 +33,18 @@ export(JAGS_diagnostics)
export(JAGS_diagnostics_autocorrelation)
export(JAGS_diagnostics_density)
export(JAGS_diagnostics_trace)
+export(JAGS_estimates_table)
export(JAGS_evaluate_formula)
export(JAGS_fit)
export(JAGS_formula)
export(JAGS_get_inits)
+export(JAGS_inference_table)
export(JAGS_marglik_parameters)
export(JAGS_marglik_parameters_formula)
export(JAGS_marglik_priors)
export(JAGS_marglik_priors_formula)
export(JAGS_parameter_names)
+export(JAGS_summary_table)
export(JAGS_to_monitor)
export(add_column)
export(ccdf)
@@ -68,12 +71,14 @@ export(interpret)
export(is.prior)
export(is.prior.PEESE)
export(is.prior.PET)
+export(is.prior.discrete)
export(is.prior.dummy)
export(is.prior.factor)
export(is.prior.none)
export(is.prior.orthonormal)
export(is.prior.point)
export(is.prior.simple)
+export(is.prior.spike_and_slab)
export(is.prior.vector)
export(is.prior.weightfunction)
export(lines_prior_list)
@@ -110,6 +115,7 @@ export(prior_factor)
export(prior_informed)
export(prior_informed_medicine_names)
export(prior_none)
+export(prior_spike_and_slab)
export(prior_weightfunction)
export(qpoint)
export(quant)
@@ -121,6 +127,8 @@ export(rtwo.sided)
export(rtwo.sided_fixed)
export(runjags_estimates_empty_table)
export(runjags_estimates_table)
+export(runjags_inference_empty_table)
+export(runjags_inference_table)
export(sd)
export(transform_orthonormal_samples)
export(var)
diff --git a/NEWS.md b/NEWS.md
index 4bc613f..ea8b6a2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,15 @@
+## version 0.2.11
+### Features
+- adding Bernoulli prior distribution
+- adding spike and slab type of prior distributions (without marginal likelihood computations/model-averaging capabilities)
+- new vignette comparing Bayes factor computation via marginal likelihood and spike and slab priors
+
+### Fixes
+- when a transformation is applied, JAGS summary tables now produce the mean of the transformed variable (previous versions incorrectly returned transformation of the mean)
+
+### Changes
+- runjags_XXX_table functions are now also exported as JAGS_XXX_functions for consistency with the rest of the code
+
## version 0.2.10
### Features
- trace, density, and autocorrelation diagnostic plots for JAGS models
diff --git a/R/JAGS-diagnostics.R b/R/JAGS-diagnostics.R
index aeeba16..bde5ded 100644
--- a/R/JAGS-diagnostics.R
+++ b/R/JAGS-diagnostics.R
@@ -221,6 +221,20 @@ JAGS_diagnostics_autocorrelation <- function(fit, parameter, plot_type = "base",
# extract the relevant parameters
+ if(is.prior.spike_and_slab(prior_list[[parameter]])){
+
+ if(sum(model_samples[,colnames(model_samples) == paste0(parameter, "_indicator")] != 0) < 10)
+ stop("The parameter with a spike and slab prior did not result in enough samples under the slab for producing a diagnostic figure.")
+
+ # change the samples between conditional
+ model_samples[
+ model_samples[,colnames(model_samples) == paste0(parameter, "_indicator")] == 0,
+ colnames(model_samples) == parameter] <- NA
+
+ # modify the parameter list
+ prior_list[[parameter]] <- prior_list[[parameter]]$variable
+ }
+
if(is.prior.factor(prior_list[[parameter]])){
if(attr(prior_list[[parameter]], "levels") > 2){
model_samples <- model_samples[,paste0(parameter, "[", 1:(attr(prior_list[[parameter]], "levels")-1), "]"),drop = FALSE]
@@ -339,14 +353,14 @@ JAGS_diagnostics_autocorrelation <- function(fit, parameter, plot_type = "base",
for(i in 1:ncol(plot_data)){
if(is.null(xlim)){
- x_range <- range(plot_data[,i])
+ x_range <- range(plot_data[,i], na.rm = TRUE)
}else{
x_range <- xlim
}
for(j in seq_along(unique(chain))){
- temp_args <- list(x = plot_data[chain == j,i], n = n_points, from = x_range[1], to = x_range[2])
+ temp_args <- list(x = plot_data[chain == j,i], n = n_points, from = x_range[1], to = x_range[2], na.rm = TRUE)
temp_density <- do.call(stats::density, temp_args)
x_den <- temp_density$x
@@ -404,7 +418,7 @@ JAGS_diagnostics_autocorrelation <- function(fit, parameter, plot_type = "base",
for(i in 1:ncol(plot_data)){
if(is.null(ylim)){
- y_range <- range(plot_data[,i])
+ y_range <- range(plot_data[,i], na.rm = TRUE)
}else{
y_range <- ylim
}
@@ -458,7 +472,7 @@ JAGS_diagnostics_autocorrelation <- function(fit, parameter, plot_type = "base",
for(j in seq_along(unique(chain))){
temp_x <- 0:lags
- temp_y <- stats::acf(plot_data[chain == j,i], lag.max = lags, plot = FALSE)$acf[, , 1L]
+ temp_y <- stats::acf(plot_data[chain == j,i], lag.max = lags, plot = FALSE, na.action = stats::na.pass)$acf[, , 1L]
temp_autocor <- list(
diff --git a/R/JAGS-fit.R b/R/JAGS-fit.R
index 837ebfd..e08bdf6 100644
--- a/R/JAGS-fit.R
+++ b/R/JAGS-fit.R
@@ -359,7 +359,18 @@ JAGS_add_priors <- function(syntax, prior_list){
syntax_end <- substr(syntax, opening_bracket + 1, nchar(syntax))
# create the priors relevant syntax
+ syntax_priors <- .JAGS_add_priors.fun(prior_list)
+
+ # merge everything back together
+ syntax <- paste0(syntax_start, "\n", syntax_priors, "\n", syntax_end)
+
+ return(syntax)
+}
+
+.JAGS_add_priors.fun <- function(prior_list){
+
syntax_priors <- ""
+
for(i in seq_along(prior_list)){
if(is.prior.weightfunction(prior_list[[i]])){
@@ -370,6 +381,10 @@ JAGS_add_priors <- function(syntax, prior_list){
syntax_priors <- paste(syntax_priors, .JAGS_prior.PP(prior_list[[i]]))
+ }else if(is.prior.spike_and_slab(prior_list[[i]])){
+
+ syntax_priors <- paste(syntax_priors, .JAGS_prior.spike_and_slab(prior_list[[i]], names(prior_list)[i]))
+
}else if(is.prior.factor(prior_list[[i]])){
syntax_priors <- paste(syntax_priors, .JAGS_prior.factor(prior_list[[i]], names(prior_list)[i]))
@@ -385,13 +400,8 @@ JAGS_add_priors <- function(syntax, prior_list){
}
}
- # merge everything back together
- syntax <- paste0(syntax_start, "\n", syntax_priors, "\n", syntax_end)
-
- return(syntax)
+ return(syntax_priors)
}
-
-
.JAGS_prior.simple <- function(prior, parameter_name){
.check_prior(prior)
@@ -410,6 +420,7 @@ JAGS_add_priors <- function(syntax, prior_list){
"invgamma" = paste0("inv_",parameter_name," ~ dgamma(",prior$parameter[["shape"]],",",prior$parameter[["scale"]],")"),
"exp" = paste0(parameter_name," ~ dexp(",prior$parameter[["rate"]],")"),
"beta" = paste0(parameter_name," ~ dbeta(",prior$parameter[["alpha"]],",",prior$parameter[["beta"]],")"),
+ "bernoulli" = paste0(parameter_name," ~ dbern(",prior$parameter[["probability"]],")"),
"uniform" = paste0(parameter_name," ~ dunif(",prior$parameter[["a"]],",",prior$parameter[["b"]],")")
)
@@ -583,7 +594,33 @@ JAGS_add_priors <- function(syntax, prior_list){
return(syntax)
}
+.JAGS_prior.spike_and_slab <- function(prior, parameter_name){
+ .check_prior(prior)
+ if(!is.prior.spike_and_slab(prior))
+ stop("improper prior provided")
+ check_char(parameter_name, "parameter_name")
+
+ if(is.prior.PET(prior[["variable"]]) | is.prior.PEESE(prior[["variable"]]) | is.prior.weightfunction(prior[["variable"]]))
+ stop("Spike and slab functionality is not implemented for publication bias prior distributions.")
+ if(is.prior.spike_and_slab(prior[["variable"]]))
+ stop("Spike and slab prior distribution cannot be nested inside of a spike and slab prior distribution.")
+
+
+ prior_variable_list <- prior["variable"]
+ prior_inclusion_list <- prior["inclusion"]
+ names(prior_variable_list) <- paste0(parameter_name, "_variable")
+ names(prior_inclusion_list) <- paste0(parameter_name, "_inclusion")
+
+ syntax <- paste0(
+ .JAGS_add_priors.fun(prior_variable_list),
+ .JAGS_add_priors.fun(prior_inclusion_list),
+ parameter_name, "_indicator ~ dbern(", paste0(parameter_name, "_inclusion"), ")\n",
+ parameter_name, " = ", parameter_name, "_variable * ", parameter_name, "_indicator\n"
+ )
+
+ return(syntax)
+}
.check_JAGS_syntax <- function(syntax){
@@ -638,52 +675,62 @@ JAGS_get_inits <- function(prior_list, chains, seed){
inits <- vector("list", chains)
for(j in 1:chains){
- temp_inits <- list()
+ temp_inits <- .JAGS_get_inits.fun(prior_list)
- for(i in seq_along(prior_list)){
+ temp_inits[[".RNG.seed"]] <- seed + j
+ temp_inits[[".RNG.name"]] <- "base::Super-Duper"
- if(is.prior.point(prior_list[[i]])){
+ inits[[j]] <- temp_inits
+ }
- next
+ return(inits)
+}
- }else if(is.prior.weightfunction(prior_list[[i]])){
+.JAGS_get_inits.fun <- function(prior_list){
- temp_inits <- c(temp_inits, .JAGS_init.weightfunction(prior_list[[i]]))
+ temp_inits <- list()
- }else if(is.prior.PET(prior_list[[i]]) | is.prior.PEESE(prior_list[[i]])){
+ for(i in seq_along(prior_list)){
- temp_inits <- c(temp_inits, .JAGS_init.PP(prior_list[[i]]))
+ if(is.prior.point(prior_list[[i]])){
- }else if(is.prior.factor(prior_list[[i]])){
+ next
- temp_inits <- c(temp_inits, .JAGS_init.factor(prior_list[[i]], names(prior_list)[i]))
+ }else if(is.prior.weightfunction(prior_list[[i]])){
- }else if(is.prior.vector(prior_list[[i]])){
+ temp_inits <- c(temp_inits, .JAGS_init.weightfunction(prior_list[[i]]))
- temp_inits <- c(temp_inits, .JAGS_init.vector(prior_list[[i]], names(prior_list)[i]))
+ }else if(is.prior.PET(prior_list[[i]]) | is.prior.PEESE(prior_list[[i]])){
- }else if(is.prior.simple(prior_list[[i]])){
+ temp_inits <- c(temp_inits, .JAGS_init.PP(prior_list[[i]]))
- temp_inits <- c(temp_inits, .JAGS_init.simple(prior_list[[i]], names(prior_list)[i]))
+ }else if(is.prior.spike_and_slab(prior_list[[i]])){
- }
- }
+ temp_inits <- c(temp_inits, .JAGS_init.spike_and_slab(prior_list[[i]], names(prior_list)[i]))
- temp_inits[[".RNG.seed"]] <- seed + j
- temp_inits[[".RNG.name"]] <- "base::Super-Duper"
+ }else if(is.prior.factor(prior_list[[i]])){
- inits[[j]] <- temp_inits
- }
+ temp_inits <- c(temp_inits, .JAGS_init.factor(prior_list[[i]], names(prior_list)[i]))
- return(inits)
-}
+ }else if(is.prior.vector(prior_list[[i]])){
+
+ temp_inits <- c(temp_inits, .JAGS_init.vector(prior_list[[i]], names(prior_list)[i]))
+
+ }else if(is.prior.simple(prior_list[[i]])){
+ temp_inits <- c(temp_inits, .JAGS_init.simple(prior_list[[i]], names(prior_list)[i]))
+ }
+ }
+
+ return(temp_inits)
+}
.JAGS_init.simple <- function(prior, parameter_name){
.check_prior(prior)
if(!is.prior.simple(prior))
stop("improper prior provided")
+ check_char(parameter_name, "parameter_name")
if(prior[["distribution"]] == "point"){
@@ -714,6 +761,7 @@ JAGS_get_inits <- function(prior_list, chains, seed){
.check_prior(prior)
if(!is.prior.vector(prior))
stop("improper prior provided")
+ check_char(parameter_name, "parameter_name")
if(prior[["distribution"]] == "point"){
@@ -799,6 +847,23 @@ JAGS_get_inits <- function(prior_list, chains, seed){
return(init)
}
+.JAGS_init.spike_and_slab <- function(prior, parameter_name){
+
+ .check_prior(prior)
+ if(!is.prior.spike_and_slab(prior))
+ stop("improper prior provided")
+
+ prior_variable <- prior["variable"]
+ names(prior_variable) <- paste0(parameter_name, "_variable")
+ init <- .JAGS_get_inits.fun(prior_variable)
+
+ if(!is.prior.point(prior[["inclusion"]])){
+ init[[paste0(parameter_name, "_inclusion")]] <- rng(prior[["inclusion"]], 1)
+ }
+
+
+ return(init)
+}
#' @title Create list of monitored parameters for 'JAGS' model
@@ -836,6 +901,10 @@ JAGS_to_monitor <- function(prior_list){
monitor <- c(monitor, .JAGS_monitor.PP(prior_list[[i]]))
+ }else if(is.prior.spike_and_slab(prior_list[[i]])){
+
+ monitor <- c(monitor, .JAGS_monitor.spike_and_slab(prior_list[[i]], names(prior_list)[i]))
+
}else if(is.prior.factor(prior_list[[i]])){
monitor <- c(monitor, .JAGS_monitor.factor(prior_list[[i]], names(prior_list)[i]))
@@ -860,8 +929,7 @@ JAGS_to_monitor <- function(prior_list){
.check_prior(prior)
if(!(is.prior.simple(prior) | is.prior.vector(prior) | is.prior.factor(prior)))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(prior[["distribution"]] == "invgamma"){
monitor <- c(parameter_name, paste0("inv_", parameter_name))
@@ -912,7 +980,27 @@ JAGS_to_monitor <- function(prior_list){
return(monitor)
}
+.JAGS_monitor.spike_and_slab <- function(prior, parameter_name){
+ .check_prior(prior)
+ if(!is.prior.spike_and_slab(prior))
+ stop("improper prior provided")
+ check_char(parameter_name, "parameter_name")
+
+ prior_variable <- prior["variable"]
+ prior_inclusion <- prior["inclusion"]
+ names(prior_variable) <- paste0(parameter_name, "_variable")
+ names(prior_inclusion) <- paste0(parameter_name, "_inclusion")
+
+ monitor <- c(
+ parameter_name,
+ JAGS_to_monitor(prior_variable),
+ JAGS_to_monitor(prior_inclusion),
+ paste0(parameter_name, "_indicator")
+ )
+
+ return(monitor)
+}
#' @title Check and list 'JAGS' fitting settings
#'
diff --git a/R/JAGS-formula.R b/R/JAGS-formula.R
index f5ec7b1..9d1c001 100644
--- a/R/JAGS-formula.R
+++ b/R/JAGS-formula.R
@@ -87,21 +87,45 @@ JAGS_formula <- function(formula, parameter, data, prior_list){
# assign factor contrasts to the data based on prior distributions
if(any(predictors_type == "factor")){
+
for(factor in names(predictors_type[predictors_type == "factor"])){
- if(is.prior.dummy(prior_list[[factor]])){
+
+ # select the corresponding prior for the variable
+ if(is.prior.spike_and_slab(prior_list[[factor]])){
+ this_prior <- prior_list[[factor]][["variable"]]
+ }else{
+ this_prior <- prior_list[[factor]]
+ }
+
+ if(is.prior.dummy(this_prior)){
stats::contrasts(data[,factor]) <- "contr.treatment"
- }else if(is.prior.orthonormal(prior_list[[factor]])){
+ }else if(is.prior.orthonormal(this_prior)){
stats::contrasts(data[,factor]) <- "contr.orthonormal"
- }else if(is.prior.point(prior_list[[factor]])){
- class(prior_list[[factor]]) <- c(class(prior_list[[factor]]), "prior.factor")
+ }else if(is.prior.point(this_prior)){
+ # change the prior class to behave as factor
+ class(this_prior) <- c(class(this_prior), "prior.factor")
+ if(is.prior.spike_and_slab(prior_list[[factor]])){
+ this_prior -> prior_list[[factor]][["variable"]]
+ }else{
+ this_prior -> prior_list[[factor]]
+ }
}else{
stop(paste0("Unsupported prior distribution defined for '", factor, "' factor variable. See '?prior_factor' for details."))
}
}
}
if(any(predictors_type == "continuous")){
+
for(continuous in names(predictors_type[predictors_type == "continuous"])){
- if(!is.prior.simple(prior_list[[continuous]]) || is.prior.factor(prior_list[[continuous]])){
+
+ # select the corresponding prior for the variable
+ if(is.prior.spike_and_slab(prior_list[[continuous]])){
+ this_prior <- prior_list[[continuous]][["variable"]]
+ }else{
+ this_prior <- prior_list[[continuous]]
+ }
+
+ if(is.prior.factor(this_prior)){
stop(paste0("Unsupported prior distribution defined for '", continuous, "' continuous variable. See '?prior' for details."))
}
}
@@ -143,39 +167,49 @@ JAGS_formula <- function(formula, parameter, data, prior_list){
# add remaining terms (omitting the intercept indexed as NA)
for(i in unique(terms_indexes[terms_indexes > 0])){
- attr(prior_list[[model_terms[i]]], "interaction") <- grepl("__xXx__", model_terms[i])
+ # extract the corresponding prior distribution for a given coefficient
+ if(is.prior.spike_and_slab(prior_list[[model_terms[i]]])){
+ this_prior <- prior_list[[model_terms[i]]][["variable"]]
+ }else{
+ this_prior <- prior_list[[model_terms[i]]]
+ }
+
+ # check whether the term is an interaction or not
+ attr(this_prior, "interaction") <- grepl("__xXx__", model_terms[i])
if(model_terms_type[i] == "continuous"){
- # continuous variables or interactions of continuous variables are simple predictors
+ # continuous variables or interactions of continuous variables are simple predictors
JAGS_data[[paste0(parameter, "_data_", model_terms[i])]] <- model_matrix[,terms_indexes == i]
+
formula_syntax <- c(formula_syntax, paste0(
- if(!is.null(attr(prior_list[[model_terms[i]]], "multiply_by"))) paste0(attr(prior_list[[model_terms[i]]], "multiply_by"), " * "),
+ if(!is.null(attr(this_prior, "multiply_by"))) paste0(attr(this_prior, "multiply_by"), " * "),
parameter, "_", model_terms[i],
" * ",
parameter, "_data_", model_terms[i], "[i]"
))
}else if(model_terms_type[i] == "factor"){
+
# factor variables or interactions with a factor requires factor style prior
# add levels information attributes to factors
- attr(prior_list[[model_terms[i]]], "levels") <- sum(terms_indexes == i) + 1
- if(attr(prior_list[[model_terms[i]]], "interaction")){
+ attr(this_prior, "levels") <- sum(terms_indexes == i) + 1
+ if(attr(this_prior, "interaction")){
level_names <- list()
for(sub_term in strsplit(model_terms[i], "__xXx__")[[1]]){
if(model_terms_type[sub_term] == "factor"){
level_names[[sub_term]] <- levels(data[,sub_term])
}
}
- attr(prior_list[[model_terms[i]]], "level_names") <- level_names
+ attr(this_prior, "level_names") <- level_names
}else{
- attr(prior_list[[model_terms[i]]], "level_names") <- levels(data[,model_terms[i]])
+ attr(this_prior, "level_names") <- levels(data[,model_terms[i]])
}
JAGS_data[[paste0(parameter, "_data_", model_terms[i])]] <- model_matrix[,terms_indexes == i, drop = FALSE]
formula_syntax <- c(formula_syntax, paste0(
- if(!is.null(attr(prior_list[[model_terms[i]]], "multiply_by"))) paste0(attr(prior_list[[model_terms[i]]], "multiply_by"), " * "),
+ if(!is.null(attr(this_prior, "multiply_by"))) paste0(attr(this_prior, "multiply_by"), " * "),
"inprod(",
parameter, "_", model_terms[i],
", ",
@@ -186,6 +220,12 @@ JAGS_formula <- function(formula, parameter, data, prior_list){
stop("Unrecognized model term.")
}
+ # update the corresponding prior distribution back into the prior list
+ if(is.prior.spike_and_slab(prior_list[[model_terms[i]]])){
+ this_prior -> prior_list[[model_terms[i]]][["variable"]]
+ }else{
+ this_prior -> prior_list[[model_terms[i]]]
+ }
}
@@ -288,35 +328,55 @@ JAGS_evaluate_formula <- function(fit, formula, parameter, data, prior_list){
# check that passed data correspond to the specified priors (factor levels etc...) and set the proper contrasts
if(any(predictors_type == "factor")){
- # check the proper data input
+
+ # check the proper data input for each factor prior
for(factor in names(predictors_type[predictors_type == "factor"])){
+
+ # select the corresponding prior in the variable
+ if(is.prior.spike_and_slab(prior_list_formula[[factor]])){
+ this_prior <- prior_list_formula[[factor]][["variable"]]
+ }else{
+ this_prior <- prior_list_formula[[factor]]
+ }
+
if(is.factor(data[,factor])){
- if(all(levels(data[,factor]) %in% attr(prior_list_formula[[factor]], "level_names"))){
+ if(all(levels(data[,factor]) %in% attr(this_prior, "level_names"))){
# either the formatting is correct, or the supplied levels are a subset of the original levels
# reformat to check ordering and etc...
- data[,factor] <- factor(data[,factor], levels = attr(prior_list_formula[[factor]], "level_names"))
+ data[,factor] <- factor(data[,factor], levels = attr(this_prior, "level_names"))
}else{
# there are some additional levels
stop(paste0("Levels specified in the '", factor, "' factor variable do not match the levels used for model specification."))
}
- }else if(all(unique(data[,factor]) %in% attr(prior_list_formula[[factor]], "level_names"))){
+ }else if(all(unique(data[,factor]) %in% attr(this_prior, "level_names"))){
# the variable was not passed as a factor but the values matches the factor levels
- data[,factor] <- factor(data[,factor], levels = attr(prior_list_formula[[factor]], "level_names"))
+ data[,factor] <- factor(data[,factor], levels = attr(this_prior, "level_names"))
}else{
# there are some additional mismatching values
stop(paste0("Levels specified in the '", factor, "' factor variable do not match the levels used for model specification."))
}
}
+
# set the contrast
- if(is.prior.orthonormal(prior_list_formula[[factor]])){
+ if(is.prior.orthonormal(this_prior)){
stats::contrasts(data[,factor]) <- "contr.orthonormal"
- }else if(is.prior.dummy(prior_list_formula[[factor]])){
+ }else if(is.prior.dummy(this_prior)){
stats::contrasts(data[,factor]) <- "contr.treatment"
}
}
if(any(predictors_type == "continuous")){
+
+ # check the proper data input for each continuous prior
for(continuous in names(predictors_type[predictors_type == "continuous"])){
- if(!is.prior.simple(prior_list_formula[[continuous]]) || is.prior.factor(prior_list_formula[[continuous]])){
+
+ # select the corresponding prior in the variable
+ if(is.prior.spike_and_slab(prior_list_formula[[continuous]])){
+ this_prior <- prior_list_formula[[continuous]][["variable"]]
+ }else{
+ this_prior <- prior_list_formula[[continuous]]
+ }
+
+ if(is.prior.factor(this_prior)){
stop(paste0("Unsupported prior distribution defined for '", continuous, "' continuous variable. See '?prior' for details."))
}
}
diff --git a/R/JAGS-marglik.R b/R/JAGS-marglik.R
index c4a1b98..79d312b 100644
--- a/R/JAGS-marglik.R
+++ b/R/JAGS-marglik.R
@@ -109,6 +109,9 @@ JAGS_bridgesampling <- function(fit, log_posterior, data = NULL, prior_list = NU
all_prior_list <- prior_list
}
+ if(any(sapply(all_prior_list, is.prior.discrete)))
+ stop("Discrete or spike and slab priors are not supported with bridgesampling.")
+
### extract relevant variables and upper and lower bound
bridgesampling_posterior <- JAGS_bridgesampling_posterior(posterior = posterior, prior_list = all_prior_list, add_parameters = add_parameters, add_bounds = add_bounds)
if(ncol(bridgesampling_posterior) == 0)
@@ -248,6 +251,8 @@ JAGS_bridgesampling_posterior <- function(posterior, prior_list, add_parameters
stop("lb' and 'ub' must be numeric vectors.")
}
+ if(any(sapply(prior_list, is.prior.spike_and_slab)))
+ stop("Marginal likelihood computation for spike and slab priors is not implemented.")
# get information about the specified parameters
parameters_names <- .JAGS_bridgesampling_posterior_info(prior_list)
@@ -302,6 +307,10 @@ JAGS_bridgesampling_posterior <- function(posterior, prior_list, add_parameters
add_parameter <- .JAGS_bridgesampling_posterior_info.PP(prior_list[[i]])
+ }else if(is.prior.spike_and_slab(prior_list[[i]])){
+
+ add_parameter <- .JAGS_bridgesampling_posterior_info.spike_and_slab(prior_list[[i]], names(prior_list)[i])
+
}else if(is.prior.factor(prior_list[[i]])){
add_parameter <- .JAGS_bridgesampling_posterior_info.factor(prior_list[[i]], names(prior_list)[i])
@@ -333,8 +342,7 @@ JAGS_bridgesampling_posterior <- function(posterior, prior_list, add_parameters
.check_prior(prior)
if(!is.prior.simple(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(prior[["distribution"]] == "invgamma"){
@@ -359,8 +367,7 @@ JAGS_bridgesampling_posterior <- function(posterior, prior_list, add_parameters
.check_prior(prior)
if(!is.prior.vector(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(prior$parameters[["K"]] == 1){
parameter <- parameter_name
@@ -381,8 +388,7 @@ JAGS_bridgesampling_posterior <- function(posterior, prior_list, add_parameters
.check_prior(prior)
if(!is.prior.factor(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(is.prior.point(prior) | is.prior.dummy(prior)){
@@ -462,6 +468,33 @@ JAGS_bridgesampling_posterior <- function(posterior, prior_list, add_parameters
names(attr(parameter, "lb")) <- parameter
names(attr(parameter, "ub")) <- parameter
+ return(parameter)
+}
+.JAGS_bridgesampling_posterior_info.spike_and_slab <- function(prior, parameter_name){
+
+ .check_prior(prior)
+ if(!is.prior.spike_and_slab(prior))
+ stop("improper prior provided")
+ check_char(parameter_name, "parameter_name")
+
+ if(!is.prior.point(prior[["inclusion"]])){
+
+ parameter_variable <- .JAGS_bridgesampling_posterior_info.simple(prior[["variable"]], paste0(parameter_name, "_variable"))
+ parameter_inclusion <- .JAGS_bridgesampling_posterior_info.simple(prior[["inclusion"]], paste0(parameter_name, "_inclusion"))
+
+ parameter <- c(parameter_variable, parameter_inclusion)
+
+ attr(parameter, "lb") <- c(attr(parameter_variable, "lb"), attr(parameter_inclusion, "lb"))
+ attr(parameter, "ub") <- c(attr(parameter_variable, "ub"), attr(parameter_inclusion, "ub"))
+
+ names(attr(parameter, "lb")) <- c(names(attr(parameter_variable, "lb")), names(attr(parameter_inclusion, "lb")))
+ names(attr(parameter, "ub")) <- c(names(attr(parameter_variable, "ub")), names(attr(parameter_inclusion, "ub")))
+
+ }else{
+ parameter <- .JAGS_bridgesampling_posterior_info.simple(prior[["variable"]], paste0(parameter_name, "_variable"))
+ }
+
+
return(parameter)
}
@@ -510,6 +543,10 @@ JAGS_marglik_priors <- function(samples, prior_list){
marglik <- marglik + .JAGS_marglik_priors.PP(samples, prior_list[[i]])
+ }else if(is.prior.spike_and_slab(prior_list[[i]])){
+
+ marglik <- marglik + .JAGS_marglik_priors.spike_and_slab(samples, prior_list[[i]], names(prior_list)[i])
+
}else if(is.prior.factor(prior_list[[i]])){
marglik <- marglik + .JAGS_marglik_priors.factor(samples, prior_list[[i]], names(prior_list)[i])
@@ -534,8 +571,7 @@ JAGS_marglik_priors <- function(samples, prior_list){
.check_prior(prior)
if(!is.prior.simple(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(prior[["distribution"]] == "invgamma"){
@@ -562,8 +598,7 @@ JAGS_marglik_priors <- function(samples, prior_list){
.check_prior(prior)
if(!is.prior.vector(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(prior$parameters[["K"]] == 1){
marglik <- lpdf(prior, samples[[ parameter_name ]])
@@ -578,8 +613,7 @@ JAGS_marglik_priors <- function(samples, prior_list){
.check_prior(prior)
if(!is.prior.factor(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(is.prior.point(prior) | is.prior.dummy(prior)){
@@ -641,6 +675,24 @@ JAGS_marglik_priors <- function(samples, prior_list){
return(marglik)
}
+.JAGS_marglik_priors.spike_and_slab <- function(samples, prior, parameter_name){
+
+ .check_prior(prior)
+ if(!is.prior.spike_and_slab(prior))
+ stop("improper prior provided")
+ check_char(parameter_name, "parameter_name")
+
+ marglik <- 0
+ if(!is.prior.point(prior[["inclusion"]])){
+ marglik <- marglik + .JAGS_marglik_priors.simple(samples, prior[["inclusion"]], paste0(parameter_name, "_inclusion"))
+ }
+
+ if(samples[[ paste0(if(prior[["variable"]][["distribution"]] == "invgamma") "inv_" else "", parameter_name, "_variable") ]] != 0){
+ marglik <- marglik + .JAGS_marglik_priors.simple(samples, prior[["variable"]], paste0(parameter_name, "_variable"))
+ }
+
+ return(marglik)
+}
#' @rdname JAGS_marglik_priors
JAGS_marglik_priors_formula <- function(samples, formula_prior_list){
@@ -700,6 +752,10 @@ JAGS_marglik_parameters <- function(samples, prior_list){
parameters <- c(parameters, .JAGS_marglik_parameters.PP(samples, prior_list[[i]]))
+ }else if(is.prior.spike_and_slab(prior_list[[i]])){
+
+ parameters <- c(parameters, .JAGS_marglik_parameters.spike_and_slab(samples, prior_list[[i]], names(prior_list)[i]))
+
}else if(is.prior.factor(prior_list[[i]])){
parameters <- c(parameters, .JAGS_marglik_parameters.factor(samples, prior_list[[i]], names(prior_list)[i]))
@@ -724,8 +780,7 @@ JAGS_marglik_parameters <- function(samples, prior_list){
.check_prior(prior)
if(!is.prior.simple(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
parameter <- list()
@@ -744,8 +799,7 @@ JAGS_marglik_parameters <- function(samples, prior_list){
.check_prior(prior)
if(!is.prior.vector(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
parameter <- list()
@@ -762,8 +816,7 @@ JAGS_marglik_parameters <- function(samples, prior_list){
.check_prior(prior)
if(!is.prior.factor(prior))
stop("improper prior provided")
- if(!is.character(parameter_name) | length(parameter_name) != 1)
- stop("'parameter_name' must be a character vector of length 1.")
+ check_char(parameter_name, "parameter_name")
if(is.prior.point(prior) | is.prior.dummy(prior)){
@@ -838,6 +891,21 @@ JAGS_marglik_parameters <- function(samples, prior_list){
return(parameter)
}
+.JAGS_marglik_parameters.spike_and_slab <- function(samples, prior, parameter_name){
+
+ .check_prior(prior)
+ if(!is.prior.spike_and_slab(prior))
+ stop("improper prior provided")
+ check_char(parameter_name, "parameter_name")
+
+ parameter <- list()
+ parameter[paste0(parameter_name, "_variable")] <- .JAGS_marglik_parameters.simple(samples, prior[["variable"]], paste0(parameter_name, "_variable"))
+ if(!is.prior.point(prior[[parameter_name]][["inclusion"]])){
+ parameter[paste0(parameter_name, "_inclusion")] <- .JAGS_marglik_parameters.simple(samples, prior[["inclusion"]], paste0(parameter_name, "_inclusion"))
+ }
+
+ return(parameter)
+}
#' @rdname JAGS_marglik_parameters
JAGS_marglik_parameters_formula <- function(samples, formula_data_list, formula_prior_list, prior_list_parameters){
diff --git a/R/priors-density.R b/R/priors-density.R
index b86a6df..ff7e68c 100644
--- a/R/priors-density.R
+++ b/R/priors-density.R
@@ -90,6 +90,10 @@ density.prior <- function(x,
x_range <- c(0, 1)
}else if(!individual & is.prior.weightfunction(x)){
x_range <- c(0, 1)
+ }else if(is.prior.spike_and_slab(x)){
+ x_range <- range(x[["variable"]][["truncation"]]["lower"], x[["variable"]][["truncation"]]["upper"], 0)
+ }else if(is.prior.discrete(x)){
+ x_range <- c(x[["truncation"]]["lower"], x[["truncation"]]["upper"])
}else{
x_range <- range(x, if(is.null(x_range_quant)) .range.prior_quantile_default(x) else x_range_quant)
}
@@ -98,7 +102,11 @@ density.prior <- function(x,
# get the x_seq for plotting
if(is.null(x_seq)){
- x_seq <- seq(x_range[1], x_range[2], length.out = n_points)
+ if(is.prior.discrete(x)){
+ x_seq <- seq(x_range[1], x_range[2], by = 1)
+ }else{
+ x_seq <- seq(x_range[1], x_range[2], length.out = n_points)
+ }
}
# specify it on the transformed range if requested
@@ -113,6 +121,8 @@ density.prior <- function(x,
out <- .density.prior.weightfunction(x, x_seq, x_range, n_points, n_samples, force_samples, individual)
}else if(is.prior.PET(x) | is.prior.PEESE(x)){
out <- .density.prior.PETPEESE(x, x_seq, x_range, n_points, n_samples, force_samples, individual, transformation, transformation_arguments, truncate_end)
+ }else if(is.prior.spike_and_slab(x)){
+ out <- .density.prior.spike_and_slab(x, x_seq, x_range, n_points, n_samples, force_samples, transformation, transformation_arguments, truncate_end)
}else if(is.prior.orthonormal(x)){
out <- .density.prior.orthonormal(x, x_seq, x_range, n_points, n_samples, force_samples, transformation, transformation_arguments, truncate_end)
}else if(is.prior.point(x)){
@@ -131,6 +141,11 @@ density.prior <- function(x,
x_sam <- rng(x, n_samples)
x_den <- stats::density(x_sam, n = n_points, from = x_range[1], to = x_range[2])$y
}else{
+
+ if(is.prior.discrete(x)){
+ x_seq <- unique(round(x_seq))
+ }
+
x_den <- mpdf(x, x_seq)
x_sam <- NULL
}
@@ -443,7 +458,32 @@ density.prior <- function(x,
return(out)
}
+.density.prior.spike_and_slab <- function(x, x_seq, x_range, n_points, n_samples, force_samples, transformation, transformation_arguments, truncate_end){
+ density_variable <- .density.prior.simple(x[["variable"]], x_seq, x_range, n_points, n_samples, force_samples, transformation, transformation_arguments, truncate_end)
+ density_inclusion <- .density.prior.point(prior(distribution = "spike", parameters = list(location = 0)), x_seq, x_range, n_points, n_samples, force_samples, transformation, transformation_arguments)
+
+ density_variable$y <- density_variable[["y"]] * mean(x[["inclusion"]])
+ density_inclusion$y <- density_inclusion[["y"]] * (1 - mean(x[["inclusion"]]))
+
+ attr(density_variable, "y_range") <- attr(density_variable, "y_range") * mean(x[["inclusion"]])
+ attr(density_inclusion, "y_range") <- attr(density_inclusion, "y_range") * (1 - mean(x[["inclusion"]]))
+
+ # create the output object
+ out <- list(
+ call = call("density", print(x, silent = TRUE)),
+ variable = density_variable,
+ inclusion = density_inclusion
+ )
+
+
+ class(out) <- c("density", "density.prior.spike_and_slab")
+ attr(out, "x_range") <- x_range
+ attr(out, "y_range_variable") <- attr(density_variable, "y_range")
+ attr(out, "y_range_inclusion") <- attr(density_inclusion, "y_range")
+
+ return(out)
+}
#' @title Prior range
diff --git a/R/priors-plot.R b/R/priors-plot.R
index 657c01e..ab4617f 100644
--- a/R/priors-plot.R
+++ b/R/priors-plot.R
@@ -60,6 +60,9 @@ plot.prior <- function(x, plot_type = "base",
if(is.null(xlim) & is.null(x_seq)){
if((is.prior.PET(x) | is.prior.PEESE(x) | is.prior.weightfunction(x)) & !individual){
xlim <- c(0, 1)
+ }else if(is.prior.spike_and_slab(x)){
+ xlim <- range(c(0,range(x[["variable"]], quantiles = x_range_quant)))
+ xlim <- range(pretty(xlim))
}else{
xlim <- range(x, quantiles = x_range_quant)
xlim <- range(pretty(xlim))
@@ -140,6 +143,25 @@ plot.prior <- function(x, plot_type = "base",
}
}
+ # spike and slab prior plots
+ if(is.prior.spike_and_slab(x)){
+ plots <- .plot.prior.spike_and_slab(x = x, plot_type = plot_type, plot_data = plot_data, par_name = par_name, ...)
+ if(plot_type == "ggplot"){
+ return(plots)
+ }else{
+ return(invisible())
+ }
+ }
+
+ # discrete prior plots
+ if(is.prior.discrete(x)){
+ plots <- .plot.prior.discrete(x = x, plot_type = plot_type, plot_data = plot_data, par_name = par_name, ...)
+ if(plot_type == "ggplot"){
+ return(plots)
+ }else{
+ return(invisible())
+ }
+ }
# default prior plots
if(is.prior.simple(x)){
@@ -244,6 +266,49 @@ plot.prior <- function(x, plot_type = "base",
return(plot)
}
}
+.plot.prior.discrete <- function(x, plot_type, plot_data, par_name = NULL, ...){
+
+ # get default plot settings
+ dots <- list(...)
+
+ xlim <- attr(plot_data, "x_range")
+ ylim <- attr(plot_data, "y_range")
+
+ short_name <- if(is.null(dots[["short_name"]])) FALSE else dots[["short_name"]]
+ parameter_names <- if(is.null(dots[["parameter_names"]])) FALSE else dots[["parameter_names"]]
+
+ main <- if(!is.null(attr(plot_data, "steps"))) print(x, plot = TRUE, short_name = short_name, parameter_names = parameter_names) else ""
+ xlab <- if(!is.null(attr(plot_data, "steps"))) bquote(omega["["*.(attr(plot_data, "steps")[1])*","~.(attr(plot_data, "steps")[2])*"]"]) else bquote(.(if(!is.null(par_name)){bquote(.(par_name)~"~")})~.(print(x, plot = TRUE, short_name = short_name, parameter_names = parameter_names)))
+ ylab <- "Probability"
+
+ # add it to the user input if desired
+ if(is.null(dots[["main"]])) dots$main <- main
+ if(is.null(dots[["xlab"]])) dots$xlab <- xlab
+ if(is.null(dots[["ylab"]])) dots$ylab <- ylab
+ if(is.null(dots[["xlim"]])) dots$xlim <- xlim
+ if(is.null(dots[["ylim"]])) dots$ylim <- ylim
+
+
+ if(plot_type == "base"){
+
+ .plot.prior_empty("simple", dots)
+ .lines.prior.discrete(plot_data, ...)
+ plot <- NULL
+
+ }else if(plot_type == "ggplot"){
+
+ plot <- .ggplot.prior_empty("simple", dots)
+ plot <- plot + .geom_prior.discrete(plot_data, ...)
+
+ }
+
+ # return the plots
+ if(plot_type == "base"){
+ return(invisible())
+ }else if(plot_type == "ggplot"){
+ return(plot)
+ }
+}
.plot.prior.weightfunction <- function(x, plot_type, plot_data, rescale_x, par_name = NULL, ...){
# get default plot settings
@@ -388,6 +453,41 @@ plot.prior <- function(x, plot_type = "base",
return(plot)
}
}
+.plot.prior.spike_and_slab <- function(x, plot_type, plot_data, par_name = NULL, ...){
+
+ # get default plot settings
+ dots <- list(...)
+
+ xlim <- attr(plot_data, "x_range")
+ ylim <- attr(plot_data, "y_range")
+
+ short_name <- if(is.null(dots[["short_name"]])) FALSE else dots[["short_name"]]
+ parameter_names <- if(is.null(dots[["parameter_names"]])) FALSE else dots[["parameter_names"]]
+
+ main <- if(!is.null(attr(plot_data, "steps"))) print(x, plot = TRUE, short_name = short_name, parameter_names = parameter_names) else ""
+ xlab <- if(!is.null(attr(plot_data, "steps"))) bquote(omega["["*.(attr(plot_data, "steps")[1])*","~.(attr(plot_data, "steps")[2])*"]"]) else bquote(.(if(!is.null(par_name)){bquote(.(par_name)~"~")})~.(print(x, plot = TRUE, short_name = short_name, parameter_names = parameter_names)))
+ ylab <- "Density"
+
+ # add it to the user input if desired
+ if(is.null(dots[["main"]])) dots$main <- main
+ if(is.null(dots[["xlab"]])) dots$xlab <- xlab
+ if(is.null(dots[["ylab"]])) dots$ylab <- ylab
+ if(is.null(dots[["xlim"]])) dots$xlim <- xlim
+ if(is.null(dots[["ylim"]])) dots$ylim <- ylim
+
+ args_prior <- dots
+ args_prior$plot_data <- plot_data
+ args_prior$plot_type <- plot_type
+
+ plot <- do.call(.plot_prior_list.both, args_prior)
+
+ # return the plots
+ if(plot_type == "base"){
+ return(invisible())
+ }else if(plot_type == "ggplot"){
+ return(plot)
+ }
+}
.plot.prior_empty <- function(type, dots = list(), ...){
@@ -552,7 +652,8 @@ plot.prior <- function(x, plot_type = "base",
col.fill = "grey80",
lwd = 1,
lty = 1,
- scale_y2 = 1.10
+ scale_y2 = 1.10,
+ width = 0.20
))
}
.get_scale_y2 <- function(plot_data, ...){
@@ -633,6 +734,9 @@ lines.prior <- function(x, xlim = NULL, x_seq = NULL, x_range_quant = NULL, n_po
if(is.null(xlim) & is.null(x_seq)){
if((is.prior.PET(x) | is.prior.PEESE(x) | is.prior.weightfunction(x)) & !individual){
xlim <- c(0, 1)
+ }else if(is.prior.spike_and_slab(x)){
+ xlim <- range(c(0,range(x[["variable"]], quantiles = x_range_quant)))
+ xlim <- range(pretty(xlim))
}else{
xlim <- range(x, quantiles = x_range_quant)
xlim <- range(pretty(xlim))
@@ -676,6 +780,18 @@ lines.prior <- function(x, xlim = NULL, x_seq = NULL, x_range_quant = NULL, n_po
return(invisible())
}
+ # plot discrete prior
+ if(is.prior.discrete(x)){
+ .lines.prior.discrete(plot_data, ...)
+ return(invisible())
+ }
+
+ # plot spike and slab prior
+ if(is.prior.spike_and_slab(x)){
+ .lines.prior.spike_and_slab(plot_data, ...)
+ return(invisible())
+ }
+
# default prior plots
if(is.prior.simple(x)){
if(inherits(plot_data, "density.prior.simple")){
@@ -760,6 +876,20 @@ geom_prior <- function(x, xlim = NULL, x_seq = NULL, x_range_quant = NULL, n_po
}
+ # plot discrete prior
+ if(is.prior.discrete(x)){
+ geom <- .geom_prior.discrete(plot_data, ...)
+ return(geom)
+ }
+
+
+ # plot spike and slab prior
+ if(is.prior.spike_and_slab(x)){
+ geom <- .geom_prior.spike_and_slab(plot_data, ...)
+ return(geom)
+ }
+
+
# default prior plots
if(is.prior.simple(x)){
if(inherits(plot_data, "density.prior.simple")){
@@ -787,6 +917,23 @@ geom_prior <- function(x, xlim = NULL, x_seq = NULL, x_range_quant = NULL, n_po
return(invisible())
}
+.lines.prior.discrete <- function(plot_data, ...){
+
+ dots <- list(...)
+ col <- if(!is.null(dots[["col"]])) dots[["col"]] else .plot.prior_settings()[["col"]]
+ width <- if(!is.null(dots[["width"]])) dots[["width"]] else .plot.prior_settings()[["width"]]
+ lwd <- if(!is.null(dots[["lwd"]])) dots[["lwd"]] else .plot.prior_settings()[["lwd"]]
+ lty <- if(!is.null(dots[["lty"]])) dots[["lty"]] else .plot.prior_settings()[["lty"]]
+
+ graphics::rect(
+ xleft = plot_data$x - width/2,
+ ybottom = 0,
+ xright = plot_data$x + width/2,
+ ytop = plot_data$y,
+ lwd = lwd, lty = lty, col = col)
+
+ return(invisible())
+}
.lines.prior.point <- function(plot_data, scale_y2 = 1, ...){
dots <- list(...)
@@ -878,6 +1025,13 @@ geom_prior <- function(x, xlim = NULL, x_seq = NULL, x_range_quant = NULL, n_po
return(invisible())
}
+.lines.prior.spike_and_slab <- function(plot_data, ...){
+
+ .lines.prior.simple(plot_data[["variable"]], ...)
+ .lines.prior.point(plot_data[["inclusion"]], ...)
+
+ return(invisible())
+}
# ggplot prior plot elements
.geom_prior.simple <- function(plot_data, ...){
@@ -898,6 +1052,25 @@ geom_prior <- function(x, xlim = NULL, x_seq = NULL, x_range_quant = NULL, n_po
return(geom)
}
+.geom_prior.discrete <- function(plot_data, ...){
+
+ dots <- list(...)
+ col <- if(!is.null(dots[["col"]])) dots[["col"]] else .plot.prior_settings()[["col"]]
+ width <- if(!is.null(dots[["width"]])) dots[["width"]] else .plot.prior_settings()[["width"]]
+ lwd <- if(!is.null(dots[["lwd"]])) dots[["lwd"]] else .plot.prior_settings()[["lwd"]]
+ lty <- if(!is.null(dots[["lty"]])) dots[["lty"]] else .plot.prior_settings()[["lty"]]
+
+ geom <- geom <- ggplot2::geom_bar(
+ data = data.frame(
+ x = plot_data$x,
+ y = plot_data$y),
+ mapping = ggplot2::aes_string(
+ x = "x",
+ weight = "y"),
+ size = lwd, linetype = lty, color = col, fill = col, width = width)
+
+ return(geom)
+}
.geom_prior.point <- function(plot_data, scale_y2 = 1, ...){
dots <- list(...)
@@ -1048,3 +1221,12 @@ geom_prior <- function(x, xlim = NULL, x_seq = NULL, x_range_quant = NULL, n_po
return(geom)
}
+.geom_prior.spike_and_slab <- function(plot_data, ...){
+
+ geom <- list(
+ .geom_prior.simple(plot_data[["variable"]], ...),
+ .geom_prior.point(plot_data[["inclusion"]], ...)
+ )
+
+ return(geom)
+}
diff --git a/R/priors-print.R b/R/priors-print.R
index 4e4a429..18e1554 100644
--- a/R/priors-print.R
+++ b/R/priors-print.R
@@ -56,6 +56,8 @@ print.prior <- function(x, short_name = FALSE, parameter_names = FALSE, plot = F
output <- .print.prior.simple(x, short_name, parameter_names, plot, digits_estimates, silent)
}else if(is.prior.weightfunction(x)){
output <- .print.prior.weightfunction(x, short_name, parameter_names, plot, digits_estimates, silent)
+ }else if(is.prior.spike_and_slab(x)){
+ output <- .print.prior.spike_and_slab(x, short_name, parameter_names, plot, digits_estimates, silent)
}
@@ -99,6 +101,7 @@ print.prior <- function(x, short_name = FALSE, parameter_names = FALSE, plot = F
"invgamma" = "Ig",
"point" = "S",
"beta" = "B",
+ "bernoulli" = "Br",
"exp" = "E",
"uniform" = "U",
"mnormal" = "mN",
@@ -116,6 +119,7 @@ print.prior <- function(x, short_name = FALSE, parameter_names = FALSE, plot = F
"invgamma" = "InvGamma",
"point" = "Spike",
"beta" = "Beta",
+ "bernoulli" = "Bernoulli",
"exp" = "Exponential",
"uniform" = "Uniform",
"mnormal" = "mNormal",
@@ -297,3 +301,17 @@ print.prior <- function(x, short_name = FALSE, parameter_names = FALSE, plot = F
return(output)
}
+.print.prior.spike_and_slab <- function(x, short_name, parameter_names, plot, digits_estimates, silent){
+
+ variable <- print(x[["variable"]], short_name, parameter_names, plot, digits_estimates, silent = TRUE)
+ inclusion <- print(x[["inclusion"]], short_name, parameter_names, plot, digits_estimates, silent = TRUE)
+
+ ### combine the results together
+ if(!plot){
+ output <- paste0(variable, " * ", inclusion)
+ }else{
+ output <- bquote(.(variable) ~ "*" ~ .(inclusion))
+ }
+
+ return(output)
+}
diff --git a/R/priors-tools.R b/R/priors-tools.R
index d9a44ba..2073b5e 100644
--- a/R/priors-tools.R
+++ b/R/priors-tools.R
@@ -161,6 +161,15 @@
}
}
+.check_parameter_range <- function(parameter, name, lower, upper, include_bounds = FALSE){
+ if(include_bounds){
+ if(any(parameter < lower) | any(parameter > upper))
+ stop(paste0("The '", name, "' must be higher than ", lower, " and lower than ", upper, "."), call. = FALSE)
+ }else{
+ if(any(parameter <= lower) | any(parameter >= upper))
+ stop(paste0("The '", name, "' must be higher or equal to than ", lower, " and lower or equal to than ", upper, "."), call. = FALSE)
+ }
+}
#' @title Reports whether x is a a prior object
#'
@@ -189,6 +198,7 @@
#'
#' @export is.prior
#' @export is.prior.simple
+#' @export is.prior.discrete
#' @export is.prior.vector
#' @export is.prior.point
#' @export is.prior.none
@@ -198,6 +208,7 @@
#' @export is.prior.factor
#' @export is.prior.orthonormal
#' @export is.prior.dummy
+#' @export is.prior.spike_and_slab
#' @name is.prior
NULL
@@ -218,6 +229,10 @@ is.prior.simple <- function(x){
inherits(x, "prior.simple")
}
#' @rdname is.prior
+is.prior.discrete <- function(x){
+ inherits(x, "prior.discrete")
+}
+#' @rdname is.prior
is.prior.vector <- function(x){
inherits(x, "prior.vector")
}
@@ -245,6 +260,11 @@ is.prior.orthonormal <- function(x){
is.prior.dummy <- function(x){
inherits(x, "prior.dummy")
}
+#' @rdname is.prior
+is.prior.spike_and_slab <- function(x){
+ inherits(x, "prior.spike_and_slab")
+}
+
.check_prior <- function(prior, name = "prior"){
if(!is.prior(prior))
diff --git a/R/priors.R b/R/priors.R
index 5a41898..e3b5e8e 100644
--- a/R/priors.R
+++ b/R/priors.R
@@ -95,6 +95,8 @@ prior <- function(distribution, parameters, truncation = list(lower = -Inf, uppe
distribution <- "gamma"
}else if(distribution %in% c("beta")){
distribution <- "beta"
+ }else if(distribution %in% c("bernoulli", "bernouli", "bern")){
+ distribution <- "bernoulli"
}else if(distribution %in% c("exp", "exponential")){
distribution <- "exp"
}else if(distribution %in% c("uniform", "unif")){
@@ -322,6 +324,59 @@ prior_factor <- function(distribution, parameters, truncation = list(lower = -In
return(output)
}
+
+#' @title Creates a spike and slab prior distribution
+#'
+#' @description \code{prior_spike_and_slab} creates a spike and slab prior
+#' distribution corresponding to the specification in
+#' \insertCite{kuo1998variable;textual}{BayesTools} (see
+#' \insertCite{ohara2009review;textual}{BayesTools} for further details). I.e.,
+#' a prior distribution is multiplied by an independent indicator with values
+#' either zero or one.
+#'
+#' @param prior_parameter a prior distribution for the parameter
+#' @param prior_inclusion a prior distribution for the inclusion probability. The
+#' inclusion probability must be bounded within 0 and 1 range. Defaults to
+#' \code{prior("spike", parameters = list(location = 0.5))} which corresponds to 1/2
+#' prior probability of including the slab prior distribution (but other prior
+#' distributions, like beta etc can be also specified).
+#'
+#'
+#' @examples
+#' # create a spike and slab prior distribution
+#' p1 <- prior_spike_and_slab(
+#' prior(distribution = "normal", parameters = list(mean = 0, sd = 1)),
+#' prior_inclusion = prior(distribution = "beta", parameters = list(alpha = 1, beta = 1))
+#' )
+#'
+#' @return return an object of class 'prior'.
+#'
+#' @inheritParams prior
+#' @seealso [prior()]
+#' @export
+prior_spike_and_slab <- function(prior_parameter,
+ prior_inclusion = prior(distribution = "spike", parameters = list(location = 0.5)),
+ prior_weights = 1){
+ if(!is.prior(prior_parameter))
+ stop("'prior_parameter' must be a prior distribution")
+ if(!is.prior(prior_inclusion))
+ stop("'prior_inclusion' must be a prior distribution")
+ if(is.prior.point(prior_inclusion) && (prior_inclusion$parameters[["location"]] < 0 | prior_inclusion$parameters[["location"]] > 1))
+ stop("The probability parameter of 'prior_inclusion' must be within 0 and 1.")
+ if(!is.prior.point(prior_inclusion) && (prior_inclusion$truncation[["lower"]] < 0 | prior_inclusion$truncation[["upper"]] > 1))
+ stop("The range of the probability parameter (set via the 'truncation' argument) of 'prior_inclusion' must be within 0 and 1.")
+
+ output <- list(
+ variable = prior_parameter,
+ inclusion = prior_inclusion
+ )
+
+ class(output) <- c("prior", "prior.spike_and_slab")
+
+ return(output)
+}
+
+
#### functions for constructing prior distributions ####
.prior_normal <- function(parameters, truncation){
@@ -523,6 +578,27 @@ prior_factor <- function(distribution, parameters, truncation = list(lower = -In
return(output)
}
+.prior_bernoulli <- function(parameters, truncation){
+
+ output <- list()
+
+ # check overall settings
+ parameters <- .check_and_name_parameters(parameters, "probability", "bernoulli")
+ truncation <- .check_and_set_truncation(truncation, lower = 0, upper = 1)
+
+ # check individual parameters
+ .check_parameter(parameters$probability, "probability")
+ .check_parameter_range(parameters$probability, "probability", lower = 0, upper = 1, include_bounds = TRUE)
+
+ # add the values to the output
+ output$distribution <- "bernoulli"
+ output$parameters <- parameters
+ output$truncation <- truncation
+
+ class(output) <- c("prior", "prior.simple", "prior.discrete")
+
+ return(output)
+}
.prior_uniform <- function(parameters, truncation){
output <- list()
@@ -842,6 +918,7 @@ rng.prior <- function(x, n, ...){
"gamma" = stats::rgamma(nn, shape = prior$parameters[["shape"]], rate = prior$parameters[["rate"]]),
"invgamma" = extraDistr::rinvgamma(nn, alpha = prior$parameters[["shape"]], beta = prior$parameters[["scale"]]),
"beta" = stats::rbeta(nn, shape1 = prior$parameters[["alpha"]], shape2 = prior$parameters[["beta"]]),
+ "bernoulli" = stats::rbinom(nn, size = 1, prob = prior$parameters[["probability"]]),
"exp" = stats::rexp(nn, rate = prior$parameters[["rate"]]),
"uniform" = stats::runif(nn, min = prior$parameters[["a"]], max = prior$parameters[["b"]]),
"point" = rpoint(n, location = prior$parameters[["location"]])
@@ -849,6 +926,11 @@ rng.prior <- function(x, n, ...){
x <- c(x, temp_x[temp_x >= prior$truncation[["lower"]] & temp_x <= prior$truncation[["upper"]]])
}
+ # make sure the enough samples were generated
+ if(length(x) < n){
+ x <- rng(prior, n)
+ }
+
x <- x[1:n]
}else if(is.prior.orthonormal(prior)){
@@ -912,6 +994,10 @@ rng.prior <- function(x, n, ...){
"two.sided.fixed" = rtwo.sided_fixed(n, omega = prior$parameters[["omega"]])
)
+ }else if(is.prior.spike_and_slab(prior)){
+
+ x <- rng(prior[["variable"]], n) * stats::rbinom(n, size = 1, prob = rng(prior[["inclusion"]], n))
+
}
return(x)
@@ -941,6 +1027,7 @@ cdf.prior <- function(x, q, ...){
"gamma" = stats::pgamma(q[!q_lower], shape = prior$parameters[["shape"]], rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"invgamma" = extraDistr::pinvgamma(q[!q_lower], alpha = prior$parameters[["shape"]], beta = prior$parameters[["scale"]], lower.tail = TRUE, log.p = FALSE),
"beta" = stats::pbeta(q[!q_lower], shape1 = prior$parameters[["alpha"]], shape2 = prior$parameters[["beta"]], lower.tail = TRUE, log.p = FALSE),
+ "bernoulli" = stats::pbinom(q[!q_lower], size = 1, prob = prior$parameters[["probability"]], lower.tail = TRUE, log.p = FALSE),
"exp" = stats::pexp(q[!q_lower], rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"uniform" = stats::punif(q[!q_lower], min = prior$parameters[["a"]], max = prior$parameters[["b"]], lower.tail = TRUE, log.p = FALSE),
"point" = ppoint(q[!q_lower], location = prior$parameters[["location"]], lower.tail = TRUE, log.p = FALSE)
@@ -955,6 +1042,10 @@ cdf.prior <- function(x, q, ...){
stop("Only marginal cdfs are implemented for prior weightfunctions.")
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No cdfs are implemented for spike and slab priors.")
+
}
return(p)
@@ -984,6 +1075,7 @@ ccdf.prior <- function(x, q, ...){
"gamma" = stats::pgamma(q[!q_higher], shape = prior$parameters[["shape"]], rate = prior$parameters[["rate"]], lower.tail = FALSE, log.p = FALSE),
"invgamma" = extraDistr::pinvgamma(q[!q_higher], alpha = prior$parameters[["shape"]], beta = prior$parameters[["scale"]], lower.tail = FALSE, log.p = FALSE),
"beta" = stats::pbeta(q[!q_higher], shape1 = prior$parameters[["alpha"]], shape2 = prior$parameters[["beta"]], lower.tail = FALSE, log.p = FALSE),
+ "bernoulli" = stats::pbinom(q[!q_higher], size = 1, prob = prior$parameters[["probability"]], lower.tail = FALSE, log.p = FALSE),
"exp" = stats::pexp(q[!q_higher], rate = prior$parameters[["rate"]], lower.tail = FALSE, log.p = FALSE),
"uniform" = stats::punif(q[!q_higher], min = prior$parameters[["a"]], max = prior$parameters[["b"]], lower.tail = FALSE, log.p = FALSE),
"point" = ppoint(q[!q_higher], location = prior$parameters[["location"]], lower.tail = FALSE, log.p = FALSE)
@@ -998,6 +1090,10 @@ ccdf.prior <- function(x, q, ...){
stop("Only marginal ccdf functions are implemented for prior weightfunctions.")
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No ccdf are implemented for spike and slab priors.")
+
}
return(p)
@@ -1021,6 +1117,7 @@ lpdf.prior <- function(x, y, ...){
"gamma" = stats::dgamma(x, shape = prior$parameters[["shape"]], rate = prior$parameters[["rate"]], log = TRUE),
"invgamma" = extraDistr::dinvgamma(x, alpha = prior$parameters[["shape"]], beta = prior$parameters[["scale"]], log = TRUE),
"beta" = stats::dbeta(x, shape1 = prior$parameters[["alpha"]], shape2 = prior$parameters[["beta"]], log = TRUE),
+ "bernoulli" = stats::dbinom(x, size = 1, prob = prior$parameters[["probability"]], log = TRUE),
"exp" = stats::dexp(x, rate = prior$parameters[["rate"]], log = TRUE),
"uniform" = stats::dunif(x, min = prior$parameters[["a"]], max = prior$parameters[["b"]], log = TRUE),
"point" = dpoint(x, location = prior$parameters[["location"]], log = TRUE)
@@ -1055,6 +1152,10 @@ lpdf.prior <- function(x, y, ...){
stop("Only marginal lpdf are implemented for prior weightfunctions.")
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No lpdf are implemented for spike and slab priors.")
+
}
return(log_lik)
@@ -1093,6 +1194,7 @@ quant.prior <- function(x, p, ...){
"gamma" = stats::qgamma(p, shape = prior$parameters[["shape"]], rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"invgamma" = extraDistr::qinvgamma(p, alpha = prior$parameters[["shape"]], beta = prior$parameters[["scale"]], lower.tail = TRUE, log.p = FALSE),
"beta" = stats::qbeta(p, shape1 = prior$parameters[["alpha"]], shape2 = prior$parameters[["beta"]], lower.tail = TRUE, log.p = FALSE),
+ "bernoulli" = stats::qbinom(p, size = 1, prob = prior$parameters[["probability"]], lower.tail = TRUE, log.p = FALSE),
"exp" = stats::qexp(p, rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"uniform" = stats::qunif(p, min = prior$parameters[["a"]], max = prior$parameters[["b"]], lower.tail = TRUE, log.p = FALSE),
"point" = qpoint(p, location = prior$parameters[["location"]], lower.tail = TRUE, log.p = FALSE)
@@ -1132,6 +1234,10 @@ quant.prior <- function(x, p, ...){
stop("Only marginal quantile functions are implemented for prior weightfunctions.")
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No quantile functions are implemented for spike and slab priors.")
+
}
return(q)
@@ -1148,6 +1254,7 @@ quant.prior <- function(x, p, ...){
"gamma" = stats::pgamma(prior$truncation[["lower"]], shape = prior$parameters[["shape"]], rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"invgamma" = extraDistr::pinvgamma(prior$truncation[["lower"]], alpha = prior$parameters[["shape"]], beta = prior$parameters[["scale"]], lower.tail = TRUE, log.p = FALSE),
"beta" = stats::pbeta(prior$truncation[["lower"]], shape1 = prior$parameters[["alpha"]], shape2 = prior$parameters[["beta"]], lower.tail = TRUE, log.p = FALSE),
+ "bernoulli" = stats::pbinom(prior$truncation[["lower"]] - 1, size = 1, prob = prior$parameters[["probability"]], lower.tail = TRUE, log.p = FALSE),
"exp" = stats::pexp(prior$truncation[["lower"]], rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"uniform" = stats::punif(prior$truncation[["lower"]], min = prior$parameters[["a"]], max = prior$parameters[["b"]], lower.tail = TRUE, log.p = FALSE),
"point" = ppoint(prior$truncation[["lower"]], location = prior$parameters[["location"]], lower.tail = TRUE, log.p = FALSE)
@@ -1169,6 +1276,7 @@ quant.prior <- function(x, p, ...){
"gamma" = stats::pgamma(prior$truncation[["upper"]], shape = prior$parameters[["shape"]], rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"invgamma" = extraDistr::pinvgamma(prior$truncation[["upper"]], alpha = prior$parameters[["shape"]], beta = prior$parameters[["scale"]], lower.tail = TRUE, log.p = FALSE),
"beta" = stats::pbeta(prior$truncation[["upper"]], shape1 = prior$parameters[["alpha"]], shape2 = prior$parameters[["beta"]], lower.tail = TRUE, log.p = FALSE),
+ "bernoulli" = stats::pbinom(prior$truncation[["upper"]] + 1, size = 1, prob = prior$parameters[["probability"]], lower.tail = TRUE, log.p = FALSE),
"exp" = stats::pexp(prior$truncation[["upper"]], rate = prior$parameters[["rate"]], lower.tail = TRUE, log.p = FALSE),
"uniform" = stats::punif(prior$truncation[["upper"]], min = prior$parameters[["a"]], max = prior$parameters[["b"]], lower.tail = TRUE, log.p = FALSE),
"point" = ppoint(prior$truncation[["upper"]], location = prior$parameters[["location"]], lower.tail = TRUE, log.p = FALSE)
@@ -1199,6 +1307,7 @@ quant.prior <- function(x, p, ...){
"gamma" = isTRUE(all.equal(prior$truncation[["lower"]], 0)) & is.infinite(prior$truncation[["upper"]]),
"invgamma" = isTRUE(all.equal(prior$truncation[["lower"]], 0)) & is.infinite(prior$truncation[["upper"]]),
"beta" = isTRUE(all.equal(prior$truncation[["lower"]], 0)) & isTRUE(all.equal(prior$truncation[["upper"]], 1)),
+ "bernoulli" = isTRUE(all.equal(prior$truncation[["lower"]], 0)) & isTRUE(all.equal(prior$truncation[["upper"]], 1)),
"exp" = isTRUE(all.equal(prior$truncation[["lower"]], 0)) & is.infinite(prior$truncation[["upper"]]),
"uniform" = TRUE,
"point" = TRUE,
@@ -1262,6 +1371,10 @@ mcdf.prior <- function(x, q, ...){
"mt" = extraDistr::plst(q, df = prior$parameters[["df"]], mu = 0, sigma = par2)
)
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No cdf are implemented for spike and slab priors.")
+
}
return(p)
@@ -1314,6 +1427,10 @@ mccdf.prior <- function(x, q, ...){
"mnormal" = stats::pnorm(q, mean = 0, sd = par2, lower.tail = FALSE),
"mt" = extraDistr::plst(q, df = prior$parameters[["df"]], mu = 0, sigma = par2, lower.tail = FALSE)
)
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No mccdf are implemented for spike and slab priors.")
+
}
return(p)
@@ -1367,6 +1484,10 @@ mlpdf.prior <- function(x, y, ...){
"mnormal" = stats::dnorm(x, mean = 0, sd = par2, log = TRUE),
"mt" = extraDistr::dlst(x, df = prior$parameters[["df"]], mu = 0, sigma = par2, log = TRUE)
)
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No lpdf are implemented for spike and slab priors.")
+
}
return(log_lik)
@@ -1434,6 +1555,10 @@ mquant.prior <- function(x, p, ...){
"mt" = extraDistr::qlst(p, df = prior$parameters[["df"]], mu = 0, sigma = par2)
)
+ }else if(is.prior.spike_and_slab(prior)){
+
+ stop("No quantile functions are implemented for spike and slab priors.")
+
}
return(q)
@@ -1565,6 +1690,7 @@ mean.prior <- function(x, ...){
"gamma" = x$parameters[["shape"]] / x$parameters[["rate"]],
"invgamma" = ifelse(x$parameters[["shape"]] > 1, x$parameters[["scale"]]/(x$parameters[["shape"]] - 1), NaN),
"beta" = x$parameters[["alpha"]] / (x$parameters[["alpha"]] + x$parameters[["beta"]]),
+ "bernoulli" = x$parameters[["probability"]],
"exp" = 1 / x$parameters[["rate"]],
"uniform" = (x$parameters[["a"]] + x$parameters[["b"]]) / 2,
"point" = x$parameters[["location"]]
@@ -1614,6 +1740,10 @@ mean.prior <- function(x, ...){
m <- 0
}
+ }else if(is.prior.spike_and_slab(x)){
+
+ m <- mean(x[["variable"]]) * mean(x[["inclusion"]])
+
}
return(m)
@@ -1698,6 +1828,7 @@ var.prior <- function(x, ...){
"gamma" = x$parameters[["shape"]] / x$parameters[["rate"]]^2,
"invgamma" = ifelse(x$parameters[["shape"]] > 2, x$parameters[["scale"]]^2 / (x$parameters[["shape"]] - 1)^2 * (x$parameters[["shape"]] - 2), NaN),
"beta" = (x$parameters[["alpha"]] * x$parameters[["beta"]]) / ((x$parameters[["alpha"]] + x$parameters[["beta"]])^2 * (x$parameters[["alpha"]] + x$parameters[["beta"]] + 1)),
+ "bernoulli" = (x$parameters[["probability"]] * (1 - x$parameters[["probability"]]) ),
"exp" = 1 / x$parameters[["rate"]]^2,
"uniform" = (x$parameters[["b"]] - x$parameters[["a"]])^2 / 12,
"point" = 0
@@ -1765,6 +1896,17 @@ var.prior <- function(x, ...){
"mt" = var.prior(x("t", parameters = c(location = 0, scale = par2, df = x$parameters[["df"]]))))
}
+ }else if(is.prior.spike_and_slab(x)){
+
+ # the inclusion is always beta -> indicators are betabinom
+ var_inclusion <- with(x[["inclusion"]][["parameters"]], (alpha * beta * (alpha + beta + 1) ) / ( (alpha + beta)^2 * (alpha + beta + 1) ) )
+
+
+ var <-
+ (var(x[["variable"]]) + mean(x[["variable"]])^2) *
+ (var_inclusion + mean(x[["inclusion"]])^2) -
+ (mean(x[["variable"]])^2 * mean(x[["inclusion"]])^2)
+
}
return(var)
diff --git a/R/summary-tables.R b/R/summary-tables.R
index 77eb2f6..6ac9ab9 100644
--- a/R/summary-tables.R
+++ b/R/summary-tables.R
@@ -410,10 +410,14 @@ ensemble_diagnostics_table <- function(models, parameters, title = NULL, footnot
#' posterior probability \code{prior_prob} and \code{post_prob},
#' and model inclusion Bayes factor \code{inclusion_BF}
#' @param fit runjags model fit
+#' @param conditional summarizes estimates conditional on being included
+#' in the model for spike and slab priors. Defaults to \code{FALSE}.
#' @param transformations named list of transformations to be applied
#' to specific parameters
#' @param model_description named list with additional description
#' to be added to the table
+#' @param remove_inclusion whether estimates of the inclusion probabilities
+#' should be excluded from the summary table. Defaults to \code{FALSE}.
#' @inheritParams BayesTools_ensemble_tables
#'
#'
@@ -423,9 +427,14 @@ ensemble_diagnostics_table <- function(models, parameters, title = NULL, footnot
#' returns an empty estimates table. All of the tables are objects of
#' class 'BayesTools_table'.
#'
+#' @export JAGS_summary_table
+#' @export JAGS_estimates_table
+#' @export JAGS_inference_table
#' @export model_summary_table
#' @export runjags_estimates_table
+#' @export runjags_inference_table
#' @export runjags_estimates_empty_table
+#' @export runjags_inference_empty_table
#' @name BayesTools_model_tables
#'
#' @seealso [BayesTools_ensemble_tables]
@@ -477,7 +486,7 @@ model_summary_table <- function(model, model_description = NULL, title = NULL, f
next
}else if(is.prior.weightfunction(prior_list[[i]]) | is.prior.PET(prior_list[[i]]) | is.prior.PEESE(prior_list[[i]])){
temp_prior <- print(prior_list[[i]], silent = TRUE, short_name = short_name)
- }else if(is.prior.simple(prior_list[[i]]) | is.prior.vector(prior_list[[i]]) | is.prior.factor(prior_list[[i]])){
+ }else if(is.prior.simple(prior_list[[i]]) | is.prior.vector(prior_list[[i]]) | is.prior.factor(prior_list[[i]]) | is.prior.spike_and_slab(prior_list[[i]])){
temp_prior <- paste0(names(prior_list)[i], " ~ " , print(prior_list[[i]], silent = TRUE, short_name = short_name))
}else if(is.prior.point(prior_list[[i]])){
temp_prior <- paste0(names(prior_list)[i], " = " , print(prior_list[[i]], silent = TRUE, short_name = short_name))
@@ -522,7 +531,7 @@ model_summary_table <- function(model, model_description = NULL, title = NULL, f
}
#' @rdname BayesTools_model_tables
-runjags_estimates_table <- function(fit, transformations = NULL, title = NULL, footnotes = NULL, warnings = NULL, remove_spike_0 = TRUE, transform_orthonormal = FALSE, formula_prefix = TRUE){
+runjags_estimates_table <- function(fit, transformations = NULL, title = NULL, footnotes = NULL, warnings = NULL, conditional = FALSE, remove_spike_0 = TRUE, transform_orthonormal = FALSE, formula_prefix = TRUE, remove_inclusion = FALSE){
# most of the code is shared with .diagnostics_plot_data function (keep them in sync on update)
@@ -542,6 +551,7 @@ runjags_estimates_table <- function(fit, transformations = NULL, title = NULL,
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(remove_spike_0, "remove_spike_0")
+ check_bool(conditional, "conditional")
check_bool(transform_orthonormal, "transform_orthonormal")
check_bool(formula_prefix, "formula_prefix")
@@ -555,12 +565,65 @@ runjags_estimates_table <- function(fit, transformations = NULL, title = NULL,
runjags_summary[,"Median"] <- NA
}
+ # simplify spike and slab priors to simple priors -- the samples and summary can be dealth with as any other prior
+ for(par in names(prior_list)){
+ if(is.prior.spike_and_slab(prior_list[[par]])){
+
+ # prepare parameter names
+ if(is.prior.factor(prior_list[[par]][["variable"]])){
+ if((attr(prior_list[[par]][["variable"]], "levels") - 1) == 1){
+ par_names <- par
+ }else{
+ par_names <- paste0(par, "[", 1:(attr(prior_list[[par]][["variable"]], "levels") - 1), "]")
+ }
+ }else{
+ par_names <- par
+ }
+
+ # change the samples between conditional/averaged based on the preferences
+ if(conditional){
+
+ # set the spike samples to NA
+ model_samples[
+ model_samples[,colnames(model_samples) == paste0(par, "_indicator")] == 0,
+ colnames(model_samples) %in% par_names] <- NA
+
+ # recompute the summaries
+ runjags_summary[par_names, "Mean"] <- mean(model_samples[,par_names], na.rm = TRUE)
+ runjags_summary[par_names, "Median"] <- stats::median(model_samples[,par_names], na.rm = TRUE)
+ runjags_summary[par_names, "SD"] <- sd(model_samples[,par_names], na.rm = TRUE)
+ runjags_summary[par_names, "Lower95"] <- stats::quantile(model_samples[,par_names], .025, na.rm = TRUE)
+ runjags_summary[par_names, "Upper95"] <- stats::quantile(model_samples[,par_names], .975, na.rm = TRUE)
+ }
+
+ # remove the indicator
+ runjags_summary <- runjags_summary[rownames(runjags_summary) != paste0(par, "_indicator"),,drop=FALSE]
+ model_samples <- model_samples[colnames(runjags_summary) != paste0(par, "_indicator"),,drop=FALSE]
+
+ # remove the latent variable
+ runjags_summary <- runjags_summary[!rownames(runjags_summary) %in% gsub(par, paste0(par, "_variable"), par_names),,drop=FALSE]
+ model_samples <- model_samples[!colnames(runjags_summary) %in% gsub(par, paste0(par, "_variable"), par_names),,drop=FALSE]
+
+ # remove/rename the inclusions probabilities
+ if(remove_inclusion){
+ runjags_summary <- runjags_summary[rownames(runjags_summary) != paste0(par, "_inclusion"),,drop=FALSE]
+ model_samples <- model_samples[colnames(runjags_summary) != paste0(par, "_inclusion"),,drop=FALSE]
+ }else{
+ rownames(runjags_summary)[rownames(runjags_summary) == paste0(par, "_inclusion")] <- paste0(par, " (inclusion)")
+ colnames(model_samples)[colnames(model_samples) == paste0(par, "_inclusion")] <- paste0(par, " (inclusion)")
+ }
+
+ # modify the parameter list
+ prior_list[[par]] <- prior_list[[par]]$variable
+ }
+ }
+
# apply transformations
if(!is.null(transformations)){
for(par in names(transformations)){
model_samples[,par] <- do.call(transformations[[par]][["fun"]], c(list(model_samples[,par]), transformations[[par]][["arg"]]))
- runjags_summary[par, "Mean"] <- do.call(transformations[[par]][["fun"]], c(list(runjags_summary[par, "Mean"]), transformations[[par]][["arg"]]))
- runjags_summary[par, "SD"] <- sd(model_samples[,par])
+ runjags_summary[par, "Mean"] <- mean(model_samples[,par], na.rm = TRUE)
+ runjags_summary[par, "SD"] <- sd(model_samples[,par], na.rm = TRUE)
runjags_summary[par, "Median"] <- do.call(transformations[[par]][["fun"]], c(list(runjags_summary[par, "Median"]), transformations[[par]][["arg"]]))
runjags_summary[par, "MCerr"] <- do.call(transformations[[par]][["fun"]], c(list(runjags_summary[par, "MCerr"]), transformations[[par]][["arg"]]))
runjags_summary[par, "MC.ofSD"] <- runjags_summary[par, "MCerr"] / runjags_summary[par, "SD"]
@@ -596,21 +659,41 @@ runjags_estimates_table <- function(fit, transformations = NULL, title = NULL,
model_samples <- cbind(model_samples, transformed_samples)
# update summary
- transformed_chains <- lapply(split(data.frame(transformed_samples), sort(rep(1:length(fit[["mcmc"]]), fit[["sample"]]))), coda::mcmc)
- transformed_summary <- summary(runjags::combine.mcmc(transformed_chains, collapse.chains = FALSE))
- transformed_summary <- cbind(
- Lower95 = transformed_summary$quantiles[,"2.5%"],
- Median = transformed_summary$quantiles[,"50%"],
- Upper95 = transformed_summary$quantiles[,"97.5%"],
- Mean = transformed_summary$statistics[,"Mean"],
- SD = transformed_summary$statistics[,"SD"],
- Mode = NA,
- MCerr = transformed_summary$statistics[,"Naive SE"],
- MC.ofSD = transformed_summary$statistics[,"Naive SE"] / transformed_summary$statistics[,"SD"],
- SSeff = unname(coda::effectiveSize(coda::as.mcmc(transformed_samples))),
- AC.10 = coda::autocorr.diag(coda::as.mcmc(transformed_samples), lags = 10)[1,],
- psrf = if(length(fit$mcmc)) unname(coda::gelman.diag(transformed_chains, multivariate = FALSE)$psrf[,"Point est."]) else NA
- )
+ if(anyNA(transformed_samples)){
+ # remove NA's introduced by conditional models and spike & slab priors -- also removes the
+ transformed_chains <- lapply(split(data.frame(transformed_samples), sort(rep(1:length(fit[["mcmc"]]), fit[["sample"]]))), function(x) coda::mcmc(stats::na.omit(x)))
+ transformed_summary <- summary(runjags::combine.mcmc(transformed_chains, collapse.chains = FALSE))
+ transformed_summary <- cbind(
+ Lower95 = transformed_summary$quantiles[,"2.5%"],
+ Median = transformed_summary$quantiles[,"50%"],
+ Upper95 = transformed_summary$quantiles[,"97.5%"],
+ Mean = transformed_summary$statistics[,"Mean"],
+ SD = transformed_summary$statistics[,"SD"],
+ Mode = NA,
+ MCerr = NA,
+ MC.ofSD = NA,
+ SSeff = NA,
+ AC.10 = NA,
+ psrf = NA
+ )
+ }else{
+ transformed_chains <- lapply(split(data.frame(transformed_samples), sort(rep(1:length(fit[["mcmc"]]), fit[["sample"]]))), coda::mcmc)
+ transformed_summary <- summary(runjags::combine.mcmc(transformed_chains, collapse.chains = FALSE))
+ transformed_summary <- cbind(
+ Lower95 = transformed_summary$quantiles[,"2.5%"],
+ Median = transformed_summary$quantiles[,"50%"],
+ Upper95 = transformed_summary$quantiles[,"97.5%"],
+ Mean = transformed_summary$statistics[,"Mean"],
+ SD = transformed_summary$statistics[,"SD"],
+ Mode = NA,
+ MCerr = transformed_summary$statistics[,"Naive SE"],
+ MC.ofSD = transformed_summary$statistics[,"Naive SE"] / transformed_summary$statistics[,"SD"],
+ SSeff = unname(coda::effectiveSize(coda::as.mcmc(transformed_samples))),
+ AC.10 = coda::autocorr.diag(coda::as.mcmc(transformed_samples), lags = 10)[1,],
+ psrf = if(length(fit$mcmc)) unname(coda::gelman.diag(transformed_chains, multivariate = FALSE)$psrf[,"Point est."]) else NA
+ )
+ }
+
rownames(transformed_summary) <- transformed_names
par_index <- which.max(rownames(runjags_summary) %in% par_names)
@@ -625,14 +708,14 @@ runjags_estimates_table <- function(fit, transformations = NULL, title = NULL,
# change HPD to quantile intervals
for(par in rownames(runjags_summary)){
- runjags_summary[par, "Lower95"] <- stats::quantile(model_samples[,par], .025)
- runjags_summary[par, "Upper95"] <- stats::quantile(model_samples[,par], .975)
+ runjags_summary[par, "Lower95"] <- stats::quantile(model_samples[,par], .025, na.rm = TRUE)
+ runjags_summary[par, "Upper95"] <- stats::quantile(model_samples[,par], .975, na.rm = TRUE)
}
# remove un-wanted columns
runjags_summary <- runjags_summary[,!colnames(runjags_summary) %in% c("Mode", "AC.10"),drop = FALSE]
- # remove un-wanted estimates (or support values)
+ # remove un-wanted estimates (or support values) - spike and slab priors already dealt with
for(i in seq_along(prior_list)){
if(is.prior.weightfunction(prior_list[[i]])){
# remove etas
@@ -710,6 +793,72 @@ runjags_estimates_table <- function(fit, transformations = NULL, title = NULL,
return(runjags_summary)
}
+#' @rdname BayesTools_model_tables
+runjags_inference_table <- function(fit, title = NULL, footnotes = NULL, warnings = NULL, formula_prefix = TRUE){
+
+ # check fits
+ if(!inherits(fit, "runjags"))
+ stop("'fit' must be a runjags fit")
+ if(!inherits(fit, "BayesTools_fit"))
+ stop("'fit' must be a BayesTools fit")
+ prior_list <- attr(fit, "prior_list")
+ check_list(prior_list, "prior_list")
+ if(!all(sapply(prior_list, is.prior)))
+ stop("'prior_list' must be a list of priors.")
+ check_char(title, "title", allow_NULL = TRUE)
+ check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
+ check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
+ check_bool(formula_prefix, "formula_prefix")
+
+ # return empty table if none of the priors is spike and slab
+ if(!any(sapply(prior_list, is.prior.spike_and_slab))){
+ runjags_summary <- runjags_inference_empty_table(title = title, footnotes = footnotes, warnings = warnings)
+ return(runjags_summary)
+ }
+
+ # extract samples
+ model_samples <- suppressWarnings(coda::as.mcmc(fit))
+ runjags_summary <- data.frame(matrix(nrow = 0, ncol = 4))
+
+ for(par in names(prior_list)){
+ if(is.prior.spike_and_slab(prior_list[[par]])){
+
+ temp_prior_prob <- mean(prior_list[[par]][["inclusion"]])
+ temp_post_prob <- mean(model_samples[,paste0(par, "_indicator")])
+
+ runjags_summary <- rbind(runjags_summary, data.frame(
+ Parameter = par,
+ prior_prob = temp_prior_prob,
+ post_prob = temp_post_prob,
+ inclusion_BF = (temp_post_prob / (1-temp_post_prob)) / (temp_prior_prob / (1-temp_prior_prob))
+ ))
+ }
+ }
+
+ runjags_summary$Parameter <- format_parameter_names(
+ parameters = runjags_summary$Parameter,
+ formula_parameters = unique(unlist(lapply(prior_list, attr, which = "parameter"))),
+ formula_prefix = formula_prefix)
+
+ class(runjags_summary) <- c("BayesTools_table", "BayesTools_runjags_summary", class(runjags_summary))
+ attr(runjags_summary, "type") <- c("string", "prior_prob", "post_prob", "inclusion_BF")
+ attr(runjags_summary, "rownames") <- FALSE
+ attr(runjags_summary, "title") <- title
+ attr(runjags_summary, "footnotes") <- footnotes
+ attr(runjags_summary, "warnings") <- warnings
+
+ return(runjags_summary)
+}
+
+#' @rdname BayesTools_model_tables
+JAGS_estimates_table <- runjags_estimates_table
+
+#' @rdname BayesTools_model_tables
+JAGS_inference_table <- runjags_inference_table
+
+#' @rdname BayesTools_model_tables
+JAGS_summary_table <- model_summary_table
+
#' @rdname BayesTools_model_tables
runjags_estimates_empty_table <- function(title = NULL, footnotes = NULL, warnings = NULL){
@@ -726,6 +875,22 @@ runjags_estimates_empty_table <- function(title = NULL, footnotes = NULL, warnin
return(empty_table)
}
+#' @rdname BayesTools_model_tables
+runjags_inference_empty_table <- function(title = NULL, footnotes = NULL, warnings = NULL){
+
+ empty_table <- data.frame(matrix(nrow = 0, ncol = 4))
+ colnames(empty_table) <- c("Parameter", "prior_prob", "post_prob", "BF")
+
+ class(empty_table) <- c("BayesTools_table", "BayesTools_runjags_summary", class(empty_table))
+ attr(empty_table, "type") <- c("string", "prior_prob", "post_prob", "BF")
+ attr(empty_table, "rownames") <- FALSE
+ attr(empty_table, "title") <- title
+ attr(empty_table, "footnotes") <- footnotes
+ attr(empty_table, "warnings") <- warnings
+
+ return(empty_table)
+}
+
#' @title Print a BayesTools table
#'
@@ -782,12 +947,13 @@ format_BF <- function(BF, logBF = FALSE, BF01 = FALSE, inclusion = FALSE){
check_bool(logBF, "logBF")
check_bool(BF01, "BF01")
- name <- ifelse(inclusion, "Inclusion BF", "BF")
-
if(BF01){
BF <- 1/BF
name <- ifelse(inclusion, "Exclusion BF", "1/BF")
+ }else{
+ name <- ifelse(inclusion, "Inclusion BF", "BF")
}
+
if(logBF){
BF <- log(BF)
name <- paste0("log(", name, ")")
@@ -918,8 +1084,8 @@ add_column <- function(table, column_title, column_values, column_position = NUL
"prior_prob" = "Prior prob.",
"post_prob" = "Post. prob.",
"marglik" = "log(marglik)",
- "BF" = attr(values, "name"),
- "inclusion_BF" = attr(values, "name"),
+ "BF" = if(is.null(attr(values, "name"))) "BF" else attr(values, "name"),
+ "inclusion_BF" = if(is.null(attr(values, "name"))) "Inclusion BF" else attr(values, "name"),
"n_models" = "Models",
"ESS" = "ESS",
"R_hat" = "R-hat",
diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib
index e0b570c..a53db73 100644
--- a/inst/REFERENCES.bib
+++ b/inst/REFERENCES.bib
@@ -45,28 +45,28 @@ @article{erp2017estimates
}
@article{maier2020robust,
- title = {Robust {B}ayesian Meta-Analysis: Addressing Publication Bias with Model-Averaging},
+ title = {Robust {B}ayesian meta-analysis: {A}ddressing publication bias with model-averaging},
author = {Maier, Maximilian and Barto{\v{s}}, Franti{\v{s}}ek and Wagenmakers, Eric-Jan},
- year = {in press},
+ year = {2022},
journal = {Psychological Methods},
- doi = {10.31234/osf.io/u4cns}
+ url = {10.1037/met0000405}
}
@unpublished{bartos2021no,
- title = {Need to choose: {R}obust {B}ayesian meta-analysis with competing publication bias adjustment methods},
+ title = {Robust Bayesian meta-analysis: {M}odel-averaging across complementary publication bias adjustment methods},
author = {Barto{\v{s}}, Franti{\v{s}}ek and Maier, Maximilian and Wagenmakers, Eric-Jan and Doucouliagos, Hristos and Stanley, Tom D.},
- year = {2021},
- publisher = {PsyArXiv},
+ year = {in press},
+ publisher = {Research Synthesis Methods},
note = {preprint at \url{https://doi.org/10.31234/osf.io/kvsp7}},
doi = {10.31234/osf.io/kvsp7}
}
@unpublished{bartos2020adjusting,
- title = {Adjusting for publication bias in {JASP} & {R} -- selection models, {PET-PEESE}, and robust {B}ayesian meta-analysis},
+ title = {Adjusting for publication bias in {JASP} & {R} -- {S}election models, {PET-PEESE}, and robust {B}ayesian meta-analysis},
author = {Barto{\v{s}}, Franti{\v{s}}ek and Maier, Maximilian and Quintana, Daniel S and Wagenmakers, Eric-Jan},
- year = {2021},
- journal = {PsyArXiv},
+ year = {in press},
+ journal = {Advances in Methods and Practices in Psychological Science},
note = {preprint at \url{https://doi.org/10.31234/osf.io/75bqn}},
doi = {10.31234/osf.io/75bqn}
}
@@ -269,3 +269,84 @@ @misc{BayesFactor
year = {2015},
url = {https://CRAN.R-project.org/package=BayesFactor}
}
+
+@article{kuo1998variable,
+ title = {Variable selection for regression models},
+ author = {Kuo, Lynn and Mallick, Bani},
+ journal = {Sankhy{\=a}: The Indian Journal of Statistics, Series B},
+ pages = {65--81},
+ year = {1998},
+ publisher = {JSTOR}
+}
+
+@article{ohara2009review,
+ title = {A review of Bayesian variable selection methods: what, how and which},
+ author = {O'Hara, Robert B and Sillanp{\"a}{\"a}, Mikko J},
+ journal = {Bayesian Analysis},
+ volume = {4},
+ number = {1},
+ pages = {85--117},
+ year = {2009},
+ publisher = {International Society for Bayesian Analysis},
+ doi = {10.1214/09-BA403}
+}
+
+@inproceedings{JAGS,
+ title = {{JAGS}: A program for analysis of {B}ayesian graphical models using {G}ibbs sampling},
+ author = {Plummer, Martyn},
+ year = 2003,
+ booktitle = {Proceedings of the 3rd international workshop on distributed statistical computing},
+ volume = 124,
+ number = {125.10},
+ pages = {1--10},
+ organization = {Vienna, Austria}
+}
+
+@article{rouder2019teaching,
+ title = {Teaching {B}ayes’ theorem: Strength of evidence as predictive accuracy},
+ author = {Rouder, Jeffrey N and Morey, Richard D},
+ year = 2019,
+ journal = {The American Statistician},
+ publisher = {Taylor \& Francis},
+ volume = 73,
+ number = 2,
+ pages = {186--190},
+ url = {https://doi.org/10.1080/00031305.2017.1341334}
+}
+
+@article{kass1995bayes,
+ title = {{B}ayes factors},
+ author = {Kass, R. E. and Raftery, A. E.},
+ year = 1995,
+ journal = {Journal of the American Statistical Association},
+ volume = 90,
+ number = 430,
+ pages = {773--795},
+ url = {https://doi.org/10.1080/01621459.1995.10476572}
+}
+
+@article{wrinch1921on,
+ title = {On Certain Fundamental Principles of Scientific Inquiry},
+ author = {Wrinch, D. and Jeffreys, H.},
+ year = 1921,
+ journal = {Philosophical Magazine},
+ volume = 42,
+ pages = {369--390},
+ url = {https://doi.org/10.1080/14786442108633773}
+}
+
+@misc{RoBMA,
+ title = {{RoBMA}: {A}n {R} package for robust {B}ayesian meta-analyses},
+ author = {Barto{\v{s}}, Franti{\v{s}}ek and Maier, Maximilian},
+ year = {2020},
+ note = {R package version 2.1.1},
+ url = {https://CRAN.R-project.org/package=RoBMA},
+}
+
+@misc{RoBSA,
+ title = {{RoBSA}: {A}n {R} package for robust {B}ayesian survival-analyses},
+ author = {Barto{\v{s}}, Franti{\v{s}}ek},
+ year = {2022},
+ note = {R package version 1.0.0},
+ url = {https://CRAN.R-project.org/package=RoBSA},
+}
diff --git a/man/BayesTools_model_tables.Rd b/man/BayesTools_model_tables.Rd
index d891475..2c67440 100644
--- a/man/BayesTools_model_tables.Rd
+++ b/man/BayesTools_model_tables.Rd
@@ -4,7 +4,12 @@
\alias{BayesTools_model_tables}
\alias{model_summary_table}
\alias{runjags_estimates_table}
+\alias{runjags_inference_table}
+\alias{JAGS_estimates_table}
+\alias{JAGS_inference_table}
+\alias{JAGS_summary_table}
\alias{runjags_estimates_empty_table}
+\alias{runjags_inference_empty_table}
\title{Create BayesTools model tables}
\usage{
model_summary_table(
@@ -24,12 +29,56 @@ runjags_estimates_table(
title = NULL,
footnotes = NULL,
warnings = NULL,
+ conditional = FALSE,
remove_spike_0 = TRUE,
transform_orthonormal = FALSE,
+ formula_prefix = TRUE,
+ remove_inclusion = FALSE
+)
+
+runjags_inference_table(
+ fit,
+ title = NULL,
+ footnotes = NULL,
+ warnings = NULL,
+ formula_prefix = TRUE
+)
+
+JAGS_estimates_table(
+ fit,
+ transformations = NULL,
+ title = NULL,
+ footnotes = NULL,
+ warnings = NULL,
+ conditional = FALSE,
+ remove_spike_0 = TRUE,
+ transform_orthonormal = FALSE,
+ formula_prefix = TRUE,
+ remove_inclusion = FALSE
+)
+
+JAGS_inference_table(
+ fit,
+ title = NULL,
+ footnotes = NULL,
+ warnings = NULL,
+ formula_prefix = TRUE
+)
+
+JAGS_summary_table(
+ model,
+ model_description = NULL,
+ title = NULL,
+ footnotes = NULL,
+ warnings = NULL,
+ remove_spike_0 = TRUE,
+ short_name = FALSE,
formula_prefix = TRUE
)
runjags_estimates_empty_table(title = NULL, footnotes = NULL, warnings = NULL)
+
+runjags_inference_empty_table(title = NULL, footnotes = NULL, warnings = NULL)
}
\arguments{
\item{model}{model object containing a list of \code{priors}
@@ -62,9 +111,15 @@ be printed. Defaults to \code{TRUE}.}
\item{transformations}{named list of transformations to be applied
to specific parameters}
+\item{conditional}{summarizes estimates conditional on being included
+in the model for spike and slab priors. Defaults to \code{FALSE}.}
+
\item{transform_orthonormal}{whether factors with orthonormal prior
distributions should be transformed to differences from the grand
mean}
+
+\item{remove_inclusion}{whether estimates of the inclusion probabilities
+should be excluded from the summary table. Defaults to \code{FALSE}.}
}
\value{
\code{model_summary_table} returns a table with
diff --git a/man/is.prior.Rd b/man/is.prior.Rd
index 767baa1..51a1475 100644
--- a/man/is.prior.Rd
+++ b/man/is.prior.Rd
@@ -5,6 +5,7 @@
\alias{is.prior.point}
\alias{is.prior.none}
\alias{is.prior.simple}
+\alias{is.prior.discrete}
\alias{is.prior.vector}
\alias{is.prior.PET}
\alias{is.prior.PEESE}
@@ -12,6 +13,7 @@
\alias{is.prior.factor}
\alias{is.prior.orthonormal}
\alias{is.prior.dummy}
+\alias{is.prior.spike_and_slab}
\title{Reports whether x is a a prior object}
\usage{
is.prior(x)
@@ -22,6 +24,8 @@ is.prior.none(x)
is.prior.simple(x)
+is.prior.discrete(x)
+
is.prior.vector(x)
is.prior.PET(x)
@@ -35,6 +39,8 @@ is.prior.factor(x)
is.prior.orthonormal(x)
is.prior.dummy(x)
+
+is.prior.spike_and_slab(x)
}
\arguments{
\item{x}{an object of test}
diff --git a/man/prior_spike_and_slab.Rd b/man/prior_spike_and_slab.Rd
new file mode 100644
index 0000000..bd78f70
--- /dev/null
+++ b/man/prior_spike_and_slab.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/priors.R
+\name{prior_spike_and_slab}
+\alias{prior_spike_and_slab}
+\title{Creates a spike and slab prior distribution}
+\usage{
+prior_spike_and_slab(
+ prior_parameter,
+ prior_inclusion = prior(distribution = "spike", parameters = list(location = 0.5)),
+ prior_weights = 1
+)
+}
+\arguments{
+\item{prior_parameter}{a prior distribution for the parameter}
+
+\item{prior_inclusion}{a prior distribution for the inclusion probability. The
+inclusion probability must be bounded within 0 and 1 range. Defaults to
+\code{prior("spike", parameters = list(location = 0.5))} which corresponds to 1/2
+prior probability of including the slab prior distribution (but other prior
+distributions, like beta etc can be also specified).}
+
+\item{prior_weights}{prior odds associated with a given distribution.
+The value is passed into the model fitting function, which creates models
+corresponding to all combinations of prior distributions for each of
+the model parameters and sets the model priors odds to the product
+of its prior distributions.}
+}
+\value{
+return an object of class 'prior'.
+}
+\description{
+\code{prior_spike_and_slab} creates a spike and slab prior
+distribution corresponding to the specification in
+\insertCite{kuo1998variable;textual}{BayesTools} (see
+\insertCite{ohara2009review;textual}{BayesTools} for further details). I.e.,
+a prior distribution is multiplied by an independent indicator with values
+either zero or one.
+}
+\examples{
+# create a spike and slab prior distribution
+p1 <- prior_spike_and_slab(
+ prior(distribution = "normal", parameters = list(mean = 0, sd = 1)),
+ prior_inclusion = prior(distribution = "beta", parameters = list(alpha = 1, beta = 1))
+)
+
+}
+\seealso{
+\code{\link[=prior]{prior()}}
+}
diff --git a/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-1.svg b/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-1.svg
new file mode 100644
index 0000000..733e2a5
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-1.svg
@@ -0,0 +1,64 @@
+
+
diff --git a/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-2.svg b/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-2.svg
new file mode 100644
index 0000000..0409b8e
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-2.svg
@@ -0,0 +1,190 @@
+
+
diff --git a/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-3.svg b/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-3.svg
new file mode 100644
index 0000000..436898a
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-diagnostics/diagnostics-plot-spike-and-slab-3.svg
@@ -0,0 +1,92 @@
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-fit-formula-spike-and-slab-1.svg b/tests/testthat/_snaps/JAGS-fit/jags-fit-formula-spike-and-slab-1.svg
new file mode 100644
index 0000000..3284ec0
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-fit/jags-fit-formula-spike-and-slab-1.svg
@@ -0,0 +1,196 @@
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-fit-formula-spike-and-slab-2.svg b/tests/testthat/_snaps/JAGS-fit/jags-fit-formula-spike-and-slab-2.svg
new file mode 100644
index 0000000..bf82515
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-fit/jags-fit-formula-spike-and-slab-2.svg
@@ -0,0 +1,378 @@
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-1.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-1.svg
index 49ec17d..323b00a 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-1.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-1.svg
@@ -22,26 +22,28 @@
(0, 1)
samples[, names(priors)[i]]
Density
-
-
-
-
-
--2
-0
-2
-4
-
+
+
+
+
+
+
+-4
+-2
+0
+2
+4
+
-
-
-
-
+
+
+
+
0.0
-0.1
-0.2
-0.3
-0.4
+0.1
+0.2
+0.3
+0.4
@@ -49,46 +51,54 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-10.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-10.svg
index 76f1043..663ce55 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-10.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-10.svg
@@ -33,19 +33,19 @@
3
4
5
-
+
-
-
-
-
-
+
+
+
+
+
0.00
-0.05
-0.10
-0.15
-0.20
-0.25
+0.05
+0.10
+0.15
+0.20
+0.25
@@ -53,46 +53,46 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-12.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-12.svg
index ef64fd5..78f12c6 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-12.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-12.svg
@@ -57,43 +57,43 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-13.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-13.svg
index 190ce3b..69552f1 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-13.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-13.svg
@@ -23,28 +23,28 @@
(1, 1)
samples[, names(priors)[i]]
Density
-
+
-
-
-
-
+
+
+
+
0
-2
-4
-6
-8
-
+2
+4
+6
+8
+
-
-
-
-
+
+
+
+
0.0
-0.2
-0.4
-0.6
-0.8
+0.2
+0.4
+0.6
+0.8
@@ -52,49 +52,51 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-14.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-14.svg
new file mode 100644
index 0000000..608a2e4
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-14.svg
@@ -0,0 +1,56 @@
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-2.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-2.svg
index 89950a7..3d8ef7d 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-2.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-2.svg
@@ -27,32 +27,32 @@
]
samples[, names(priors)[i]]
Density
-
+
-
-
-
-
-
-
-
-1.0
-1.5
-2.0
-2.5
-3.0
-3.5
-4.0
-4.5
-
+
+
+
+1
+2
+3
+4
+
-
-
-
+
+
+
+
+
+
+
0.0
-0.5
-1.0
-1.5
+0.2
+0.4
+0.6
+0.8
+1.0
+1.2
+1.4
@@ -60,75 +60,43 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-3.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-3.svg
index ea77232..28f6b49 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-3.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-3.svg
@@ -22,28 +22,32 @@
(0, 0.5)
samples[, names(priors)[i]]
Density
-
-
-
-
-
-
-0
-2
-4
-6
-8
-
+
+
+
+
+
+
+
+0
+1
+2
+3
+4
+5
+
-
-
-
-
+
+
+
+
+
0.0
-0.2
-0.4
-0.6
-0.8
+0.2
+0.4
+0.6
+0.8
+1.0
@@ -51,45 +55,62 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-4.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-4.svg
index a2f13cb..09e1233 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-4.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-4.svg
@@ -22,30 +22,30 @@
(0, 0.5, 5)
samples[, names(priors)[i]]
Density
-
-
-
-
-
-
-
-
--4
--2
-0
-2
-4
-6
-8
-
+
+
+
+
+
+-5
+0
+5
+10
+
-
-
-
+
+
+
+
+
+
0.0
-0.2
-0.4
-0.6
+0.1
+0.2
+0.3
+0.4
+0.5
+0.6
@@ -53,76 +53,38 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-5.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-5.svg
index 8781522..e6704fa 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-5.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-5.svg
@@ -36,17 +36,17 @@
-4
-2
0
-
+
-
-
-
-
+
+
+
+
0.0
-0.2
-0.4
-0.6
-0.8
+0.2
+0.4
+0.6
+0.8
@@ -54,56 +54,56 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-6.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-6.svg
index bd923fc..e594b95 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-6.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-6.svg
@@ -22,32 +22,32 @@
(2, 1)
samples[, names(priors)[i]]
Density
-
+
-
-
-
-
-
-
+
+
+
+
+
+
0
-2
-4
-6
-8
-10
-12
-
+2
+4
+6
+8
+10
+12
+
-
-
-
-
+
+
+
+
0.0
-0.1
-0.2
-0.3
-0.4
+0.1
+0.2
+0.3
+0.4
@@ -55,68 +55,67 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-7.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-7.svg
index 0b016b7..53a68b0 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-7.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-7.svg
@@ -34,15 +34,15 @@
2.0
2.5
3.0
-
+
-
-
-
+
+
+
0.0
-0.5
-1.0
-1.5
+0.5
+1.0
+1.5
@@ -51,45 +51,45 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-8.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-8.svg
index 94d9522..482d9e4 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-8.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-8.svg
@@ -24,30 +24,36 @@
Density
-
+
+
-
+
+
0
-2
-4
-6
-8
-
+1
+2
+3
+4
+5
+6
+
-
-
-
-
-
-
+
+
+
+
+
+
+
0.0
-0.2
-0.4
-0.6
-0.8
-1.0
-1.2
+0.2
+0.4
+0.6
+0.8
+1.0
+1.2
+1.4
@@ -55,48 +61,69 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-9.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-9.svg
index 15954c3..27d98d2 100644
--- a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-9.svg
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-9.svg
@@ -35,15 +35,15 @@
0.6
0.8
1.0
-
+
-
-
-
+
+
+
0.0
-0.5
-1.0
-1.5
+0.5
+1.0
+1.5
@@ -51,55 +51,55 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-spike-and-slab-1.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-spike-and-slab-1.svg
new file mode 100644
index 0000000..1963214
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-spike-and-slab-1.svg
@@ -0,0 +1,107 @@
+
+
diff --git a/tests/testthat/_snaps/JAGS-fit/jags-model-prior-spike-and-slab-2.svg b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-spike-and-slab-2.svg
new file mode 100644
index 0000000..1524108
--- /dev/null
+++ b/tests/testthat/_snaps/JAGS-fit/jags-model-prior-spike-and-slab-2.svg
@@ -0,0 +1,382 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-1.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-1.svg
new file mode 100644
index 0000000..81fccef
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-1.svg
@@ -0,0 +1,65 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-10.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-10.svg
new file mode 100644
index 0000000..daf462d
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-10.svg
@@ -0,0 +1,93 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-11.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-11.svg
new file mode 100644
index 0000000..8e9f1fb
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-11.svg
@@ -0,0 +1,96 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-2.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-2.svg
new file mode 100644
index 0000000..8a048a3
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-2.svg
@@ -0,0 +1,65 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-3.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-3.svg
new file mode 100644
index 0000000..b4fb12d
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-3.svg
@@ -0,0 +1,65 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-4.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-4.svg
new file mode 100644
index 0000000..c06549e
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-4.svg
@@ -0,0 +1,65 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-5.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-5.svg
new file mode 100644
index 0000000..f08988b
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-5.svg
@@ -0,0 +1,68 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-6.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-6.svg
new file mode 100644
index 0000000..665c07d
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-6.svg
@@ -0,0 +1,93 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-7.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-7.svg
new file mode 100644
index 0000000..83e3084
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-7.svg
@@ -0,0 +1,93 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-8.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-8.svg
new file mode 100644
index 0000000..2243aa7
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-8.svg
@@ -0,0 +1,93 @@
+
+
diff --git a/tests/testthat/_snaps/priors-plot/priors-plot-15-9.svg b/tests/testthat/_snaps/priors-plot/priors-plot-15-9.svg
new file mode 100644
index 0000000..373b58c
--- /dev/null
+++ b/tests/testthat/_snaps/priors-plot/priors-plot-15-9.svg
@@ -0,0 +1,92 @@
+
+
diff --git a/tests/testthat/_snaps/priors-print/priors-print-4.svg b/tests/testthat/_snaps/priors-print/priors-print-4.svg
new file mode 100644
index 0000000..5547696
--- /dev/null
+++ b/tests/testthat/_snaps/priors-print/priors-print-4.svg
@@ -0,0 +1,36 @@
+
+
diff --git a/tests/testthat/_snaps/priors/prior-bernoulli-1.svg b/tests/testthat/_snaps/priors/prior-bernoulli-1.svg
new file mode 100644
index 0000000..e6a0557
--- /dev/null
+++ b/tests/testthat/_snaps/priors/prior-bernoulli-1.svg
@@ -0,0 +1,53 @@
+
+
diff --git a/tests/testthat/_snaps/priors/prior-bernoulli-2.svg b/tests/testthat/_snaps/priors/prior-bernoulli-2.svg
new file mode 100644
index 0000000..616f0db
--- /dev/null
+++ b/tests/testthat/_snaps/priors/prior-bernoulli-2.svg
@@ -0,0 +1,55 @@
+
+
diff --git a/tests/testthat/_snaps/priors/prior-spike-and-slab-1.svg b/tests/testthat/_snaps/priors/prior-spike-and-slab-1.svg
new file mode 100644
index 0000000..a6152a5
--- /dev/null
+++ b/tests/testthat/_snaps/priors/prior-spike-and-slab-1.svg
@@ -0,0 +1,104 @@
+
+
diff --git a/tests/testthat/test-JAGS-diagnostics.R b/tests/testthat/test-JAGS-diagnostics.R
index 2ecc4ef..7a7fb0c 100644
--- a/tests/testthat/test-JAGS-diagnostics.R
+++ b/tests/testthat/test-JAGS-diagnostics.R
@@ -150,3 +150,55 @@ test_that("JAGS diagnostics work", {
expect_doppelganger("diagnostics-ggplot-autocorrelation-3.1",temp_plot[[1]])
expect_doppelganger("diagnostics-ggplot-autocorrelation-3.2",temp_plot[[2]])
})
+
+test_that("JAGS diagnostics work (spike and slab)", {
+
+ skip_on_os(c("mac", "linux", "solaris")) # multivariate sampling does not exactly match across OSes
+ skip_on_cran()
+
+ set.seed(1)
+
+ data_formula <- data.frame(
+ x_cont = rnorm(300),
+ x_fac2t = factor(rep(c("A", "B"), 150), levels = c("A", "B")),
+ x_fac3t = factor(rep(c("A", "B", "C"), 100), levels = c("A", "B", "C"))
+ )
+ data <- list(
+ y = rnorm(300, .20 * data_formula$x_cont + ifelse(data_formula$x_fac3t == "A", 0.0, ifelse(data_formula$x_fac3t == "B", -0.2, 0.4)), ifelse(data_formula$x_fac2t == "A", 0.5, 1)),
+ N = 300
+ )
+
+
+ # create model with mix of a formula and free parameters ---
+ formula_list <- list(
+ mu = ~ x_cont
+ )
+ formula_data_list <- list(
+ mu = data_formula
+ )
+ formula_prior_list <- list(
+ mu = list(
+ "intercept" = prior("normal", list(0, 5)),
+ "x_cont" = prior_spike_and_slab(prior("normal", list(0, 1)), prior_inclusion = prior("beta", list(1,1)))
+ )
+ )
+ prior_list <- list(
+ sigma = prior("lognormal", list(0, 1))
+ )
+ model_syntax <- paste0(
+ "model{\n",
+ "for(i in 1:N){\n",
+ " y[i] ~ dnorm(mu[i], 1/pow(sigma, 2))\n",
+ "}\n",
+ "}"
+ )
+
+ fit <- JAGS_fit(
+ model_syntax = model_syntax, data = data, prior_list = prior_list,
+ formula_list = formula_list, formula_data_list = formula_data_list, formula_prior_list = formula_prior_list)
+
+ ### density plots
+ expect_doppelganger("diagnostics-plot-spike-and-slab-1", function() JAGS_diagnostics_density(fit, parameter = "mu_x_cont"))
+ expect_doppelganger("diagnostics-plot-spike-and-slab-2", function() JAGS_diagnostics_autocorrelation(fit, parameter = "mu_x_cont"))
+ expect_doppelganger("diagnostics-plot-spike-and-slab-3", function() JAGS_diagnostics_trace(fit, parameter = "mu_x_cont"))
+})
diff --git a/tests/testthat/test-JAGS-fit.R b/tests/testthat/test-JAGS-fit.R
index a66b09d..131f5bb 100644
--- a/tests/testthat/test-JAGS-fit.R
+++ b/tests/testthat/test-JAGS-fit.R
@@ -17,7 +17,8 @@ test_that("JAGS model functions work (simple)", {
p10 = prior("uniform", list(1, 5)),
p11 = prior("point", list(1)),
PET = prior_PET("normal", list(0, 1)),
- PEESE = prior_PEESE("gamma", list(1, 1))
+ PEESE = prior_PEESE("gamma", list(1, 1)),
+ p12 = prior("bernoulli", list(0.75))
)
model_syntax <- JAGS_add_priors(model_syntax, priors)
@@ -32,7 +33,11 @@ test_that("JAGS model functions work (simple)", {
for(i in seq_along(priors)){
expect_doppelganger(paste0("JAGS-model-prior-",i), function(){
- hist(samples[,names(priors)[i]], breaks = 50, main = print(priors[[i]], plot = TRUE), freq = FALSE)
+ if(is.prior.discrete(priors[[i]])){
+ barplot(table(samples[,names(priors)[i]])/length(samples[,names(priors)[i]]), main = print(priors[[i]], plot = T), width = 1/(max(samples[,names(priors)[i]])+1), space = 0, xlim = c(-0.25, max(samples[,names(priors)[i]])+0.25))
+ }else{
+ hist(samples[,names(priors)[i]], breaks = 50, main = print(priors[[i]], plot = TRUE), freq = FALSE)
+ }
lines(priors[[i]], individual = TRUE)
})
}
@@ -199,6 +204,65 @@ test_that("JAGS model functions work (weightfunctions)", {
}
})
+test_that("JAGS model functions work (spike and slab)", {
+
+ skip_if_not_installed("rjags")
+ priors <- list(
+ "mu" = prior_spike_and_slab(prior("normal", list(0, 1)), prior_inclusion = prior("beta", list(1,1))),
+ "beta" = prior_spike_and_slab(prior_factor("mnormal", list(0, 1), contrast = "orthonormal"),
+ prior_inclusion = prior("beta", list(1,1)))
+ )
+
+ attr(priors$beta$variable, "levels") <- 3
+
+ for(i in 1:length(priors)){
+ model_syntax <- "model{}"
+ model_syntax <- JAGS_add_priors(model_syntax, priors[i])
+ monitor <- JAGS_to_monitor(priors[i])
+ inits <- JAGS_get_inits(priors[i], chains = 2, seed = 1)
+
+ set.seed(1)
+ model <- rjags::jags.model(file = textConnection(model_syntax), inits = inits, n.chains = 2, quiet = TRUE)
+ samples <- rjags::coda.samples(model = model, variable.names = monitor, n.iter = 5000, quiet = TRUE, progress.bar = "none")
+ samples <- do.call(rbind, samples)
+
+
+ if(i == 1){
+ expect_doppelganger(paste0("JAGS-model-prior_spike-and-slab-",i), function(){
+ temp_samples <- samples[,names(priors)[i]]
+ hs <- hist(temp_samples[temp_samples != 0], breaks = 50, plot = FALSE)
+ hs$density <- hs$density * mean(temp_samples != 0)
+ plot(hs, main = print(priors[[i]], plot = TRUE), freq = FALSE, ylim = range(c(0, max(hs$density), mean(temp_samples == 0))))
+ lines(priors[[i]], individual = TRUE)
+ })
+ }else{
+ expect_doppelganger(paste0("JAGS-model-prior_spike-and-slab-",i), function(){
+
+ oldpar <- graphics::par(no.readonly = TRUE)
+ on.exit(graphics::par(mfrow = oldpar[["mfrow"]]))
+ par(mfrow = c(1, 3))
+
+ temp_samples <- samples[,paste0(names(priors)[i], "[", 1:2, "]")] %*% t(contr.orthonormal(1:3))
+
+ hs1 <- hist(temp_samples[temp_samples[,1] != 0, 1], breaks = 50, plot = FALSE)
+ hs1$density <- hs1$density * mean(temp_samples[,1] != 0)
+ plot(hs1, main = print(priors[[i]], plot = TRUE), freq = FALSE, ylim = range(c(0, max(hs1$density), mean(temp_samples == 0))))
+ lines(priors[[i]], individual = TRUE)
+
+ hs2 <- hist(temp_samples[temp_samples[,2] != 0, 2], breaks = 50, plot = FALSE)
+ hs2$density <- hs2$density * mean(temp_samples[,1] != 0)
+ plot(hs2, main = print(priors[[i]], plot = TRUE), freq = FALSE, ylim = range(c(0, max(hs2$density), mean(temp_samples == 0))))
+ lines(priors[[i]], individual = TRUE)
+
+ hs3 <- hist(temp_samples[temp_samples[,3] != 0, 3], breaks = 50, plot = FALSE)
+ hs3$density <- hs3$density * mean(temp_samples[,1] != 0)
+ plot(hs3, main = print(priors[[i]], plot = TRUE), freq = FALSE, ylim = range(c(0, max(hs3$density), mean(temp_samples == 0))))
+ lines(priors[[i]], individual = TRUE)
+ })
+ }
+ }
+})
+
test_that("JAGS fit function works" , {
set.seed(1)
@@ -449,6 +513,90 @@ test_that("JAGS fit function integration with formula works" , {
})
+test_that("JAGS fit function integration with formula and spike and slab works" , {
+
+ skip_on_os(c("mac", "linux", "solaris")) # multivariate sampling does not exactly match across OSes
+ skip_on_cran()
+
+ set.seed(1)
+
+ data_formula <- data.frame(
+ x_cont1 = rnorm(300),
+ x_fac2t = factor(rep(c("A", "B"), 150), levels = c("A", "B")),
+ x_fac3t = factor(rep(c("A", "B", "C"), 100), levels = c("A", "B", "C"))
+ )
+ data <- list(
+ y = rnorm(300, .20 * data_formula$x_cont1 + ifelse(data_formula$x_fac3t == "A", 0.0, ifelse(data_formula$x_fac3t == "B", -0.2, 0.2)), ifelse(data_formula$x_fac2t == "A", 0.5, 1)),
+ N = 300
+ )
+
+
+ # create model with mix of a formula and free parameters ---
+ formula_list1 <- list(
+ mu = ~ x_cont1 + x_fac3t
+ )
+ formula_data_list1 <- list(
+ mu = data_formula
+ )
+ formula_prior_list1 <- list(
+ mu = list(
+ "intercept" = prior("normal", list(0, 5)),
+ "x_cont1" = prior_spike_and_slab(prior("normal", list(0, 1)), prior_inclusion = prior("beta", list(1,1))),
+ "x_fac3t" = prior_spike_and_slab(prior_factor("mnormal", list(0, 1), contrast = "orthonormal"),
+ prior_inclusion = prior("spike", list(0.5)))
+ )
+ )
+ attr(formula_prior_list1$mu$x_fac3t$variable, "multiply_by") <- "sigma"
+ prior_list1 <- list(
+ sigma = prior("lognormal", list(0, 1))
+ )
+ model_syntax1 <- paste0(
+ "model{\n",
+ "for(i in 1:N){\n",
+ " y[i] ~ dnorm(mu[i], 1/pow(sigma, 2))\n",
+ "}\n",
+ "}"
+ )
+
+ fit1 <- JAGS_fit(
+ model_syntax = model_syntax1, data = data, prior_list = prior_list1,
+ formula_list = formula_list1, formula_data_list = formula_data_list1, formula_prior_list = formula_prior_list1)
+
+ posterior1 <- suppressWarnings(coda::as.mcmc(fit1))
+
+ expect_doppelganger("JAGS-fit-formula-spike-and-slab-1", function(){
+
+ oldpar <- graphics::par(no.readonly = TRUE)
+ on.exit(graphics::par(mfcol = oldpar[["mfcol"]]))
+ par(mfrow = c(1, 3))
+
+ hist(posterior1[,"mu_x_cont1"], freq = FALSE, main = "x_cont1")
+
+ hist(posterior1[,"mu_x_cont1_variable"], freq = FALSE, main = "x_cont1_variable")
+
+ hist(posterior1[,"mu_x_cont1_inclusion"], freq = FALSE, main = "x_cont1_inclusion")
+ })
+
+ expect_doppelganger("JAGS-fit-formula-spike-and-slab-2", function(){
+
+ oldpar <- graphics::par(no.readonly = TRUE)
+ on.exit(graphics::par(mfcol = oldpar[["mfcol"]]))
+ par(mfrow = c(2, 3))
+
+ temp_samples <- posterior1[,paste0("mu_x_fac3t[", 1:2, "]")] %*% t(contr.orthonormal(1:3))
+ temp_samples_variable <- posterior1[,paste0("mu_x_fac3t_variable[", 1:2, "]")] %*% t(contr.orthonormal(1:3))
+
+ hist(temp_samples[,1], freq = FALSE, main = "x_fac3t[A]")
+ hist(temp_samples[,2], freq = FALSE, main = "x_fac3t[B]")
+ hist(temp_samples[,3], freq = FALSE, main = "x_fac3t[C]")
+
+ hist(temp_samples_variable[,1], freq = FALSE, main = "x_fac3t_variable[A]")
+ hist(temp_samples_variable[,2], freq = FALSE, main = "x_fac3t_variable[B]")
+ hist(temp_samples_variable[,3], freq = FALSE, main = "x_fac3t_variable[C]")
+ })
+
+})
+
test_that("JAGS parallel fit function works" , {
skip("requires parallel processing")
diff --git a/tests/testthat/test-JAGS-marglik.R b/tests/testthat/test-JAGS-marglik.R
index 819c2c3..6a88ebd 100644
--- a/tests/testthat/test-JAGS-marglik.R
+++ b/tests/testthat/test-JAGS-marglik.R
@@ -16,6 +16,7 @@ test_that("JAGS model functions work (simple)", {
p10 = prior("uniform", list(1, 5)),
PET = prior_PET("normal", list(0, 1)),
PEESE = prior_PEESE("gamma", list(1, 1))
+ #p12 = prior("bernoulli", list(0.75)) discrete priors are not supported with bridgesampling
)
log_posterior <- function(parameters, data){
return(0)
@@ -98,6 +99,35 @@ test_that("JAGS model functions work (factor)", {
})
+test_that("JAGS model functions work (spike and slab)", {
+ skip("Marginal likelihood computation for spike and slab priors is not implemented.")
+ skip_if_not_installed("rjags")
+ all_priors <- list(
+ p1 = prior_spike_and_slab(prior("normal", list(0, 1)), prior_inclusion = prior("beta", list(1, 1))),
+ p2 = prior_spike_and_slab(prior("gamma", list(3, 4)), prior_inclusion = prior("beta", list(5, 1))),
+ p3 = prior_spike_and_slab(prior("invgamma", list(4, 5)), prior_inclusion = prior("point", list(.3)))
+ )
+
+ log_posterior <- function(parameters, data){
+ return(0)
+ }
+
+
+ for(i in seq_along(all_priors)){
+ prior_list <- all_priors[i]
+ model_syntax <- JAGS_add_priors("model{}", prior_list)
+ monitor <- JAGS_to_monitor(prior_list)
+ inits <- JAGS_get_inits(prior_list, chains = 2, seed = 1)
+
+ set.seed(1)
+ model <- rjags::jags.model(file = textConnection(model_syntax), inits = inits, n.chains = 2, quiet = TRUE)
+ samples <- rjags::coda.samples(model = model, variable.names = monitor, n.iter = 10000, quiet = TRUE, progress.bar = "none")
+ marglik <- JAGS_bridgesampling(samples, prior_list = prior_list, data = list(), log_posterior = log_posterior)
+ expect_equal(marglik$logml, 0, tolerance = 1e-3)
+ }
+
+})
+
test_that("JAGS model functions work (weightfunctions)", {
skip_if_not_installed("rjags")
diff --git a/tests/testthat/test-priors-plot.R b/tests/testthat/test-priors-plot.R
index c758c79..96b2fd0 100644
--- a/tests/testthat/test-priors-plot.R
+++ b/tests/testthat/test-priors-plot.R
@@ -45,6 +45,26 @@ test_that("Prior plot (simple) function works", {
})
+test_that("Prior plot (discrete) function works", {
+
+ # check the default options
+ p1 <- prior("bernoulli", list(.33))
+
+ expect_doppelganger("priors-plot-15-1", function()plot(p1))
+ expect_doppelganger("priors-plot-15-2", function()plot(p1, short_name = TRUE))
+ expect_doppelganger("priors-plot-15-3", function()plot(p1, parameter_names = TRUE))
+ expect_doppelganger("priors-plot-15-4", function()plot(p1, lwd = 3, lty = 3, col = "blue"))
+ expect_doppelganger("priors-plot-15-5", function()plot(p1, par_name = "name"))
+
+ expect_doppelganger("priors-plot-15-6", plot(p1, plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-7", plot(p1, short_name = TRUE, plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-8", plot(p1, parameter_names = TRUE, plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-9", plot(p1, xlab = "xlab", ylab = "ylab", main = "main", plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-10", plot(p1, lwd = 3, lty = 3, col = "blue", plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-11", plot(p1, par_name = bquote(mu), plot_type = "ggplot"))
+
+})
+
test_that("Prior plot (point) function works", {
# check point distributions
p5 <- prior("point", list(1.5))
@@ -57,6 +77,26 @@ test_that("Prior plot (point) function works", {
expect_doppelganger("priors-plot-6-7", plot(p5, plot_type = "ggplot"))
})
+test_that("Prior plot (spike and slab) function works", {
+
+ # check the default options
+ p1 <- prior("bernoulli", list(.33))
+
+ expect_doppelganger("priors-plot-15-1", function()plot(p1))
+ expect_doppelganger("priors-plot-15-2", function()plot(p1, short_name = TRUE))
+ expect_doppelganger("priors-plot-15-3", function()plot(p1, parameter_names = TRUE))
+ expect_doppelganger("priors-plot-15-4", function()plot(p1, lwd = 3, lty = 3, col = "blue"))
+ expect_doppelganger("priors-plot-15-5", function()plot(p1, par_name = "name"))
+
+ expect_doppelganger("priors-plot-15-6", plot(p1, plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-7", plot(p1, short_name = TRUE, plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-8", plot(p1, parameter_names = TRUE, plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-9", plot(p1, xlab = "xlab", ylab = "ylab", main = "main", plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-10", plot(p1, lwd = 3, lty = 3, col = "blue", plot_type = "ggplot"))
+ expect_doppelganger("priors-plot-15-11", plot(p1, par_name = bquote(mu), plot_type = "ggplot"))
+
+})
+
test_that("Prior plot (weightfunction) function works", {
set.seed(1)
diff --git a/tests/testthat/test-priors-print.R b/tests/testthat/test-priors-print.R
index c28d48e..b583b97 100644
--- a/tests/testthat/test-priors-print.R
+++ b/tests/testthat/test-priors-print.R
@@ -109,4 +109,13 @@ test_that("Prior print function works", {
text(0.5, 0.6, print(p15, plot = TRUE))
text(0.5, 0.5, print(p16, plot = TRUE))
})
+
+ p17 <- prior_spike_and_slab(prior("gamma", list(1, 2), list(0, Inf)),
+ prior_inclusion = prior("beta", list(3, 2)))
+ expect_equal(utils::capture.output(print(p17)), "Gamma(1, 2) * Beta(3, 2)")
+ expect_equal(utils::capture.output(print(p17, short_name = TRUE)), "G(1, 2) * B(3, 2)")
+ expect_doppelganger("priors-print-4", function(){
+ empty_plot()
+ text(0.5, 1, print(p17, plot = TRUE))
+ })
})
diff --git a/tests/testthat/test-priors-tools.R b/tests/testthat/test-priors-tools.R
index 8241b5c..601d364 100644
--- a/tests/testthat/test-priors-tools.R
+++ b/tests/testthat/test-priors-tools.R
@@ -11,6 +11,10 @@ test_that("Prior handling works", {
expect_null(BayesTools:::.check_parameter_negative(0, "par", TRUE))
expect_error(BayesTools:::.check_parameter_negative(.01, "par", TRUE), "The 'par' must be non-positive.")
+ expect_error(BayesTools:::.check_parameter_range(-.01, "par", lower = 0, upper = 1, include_bounds = TRUE), "The 'par' must be higher than 0 and lower than 1.")
+ expect_error(BayesTools:::.check_parameter_range(3, "par", lower = 1, upper = 3, include_bounds = FALSE), "The 'par' must be higher or equal to than 1 and lower or equal to than 3.")
+
+
# check different ordering of names withing lists
expect_equal(prior("normal", list(0, 1)), prior("normal", list(0, 1), list(-Inf, Inf)))
expect_equal(prior("normal", list(0, 1), list(0, Inf)), prior("normal", list(0, 1), list("lower" = 0)))
diff --git a/tests/testthat/test-priors.R b/tests/testthat/test-priors.R
index 03f8f53..974c6f2 100644
--- a/tests/testthat/test-priors.R
+++ b/tests/testthat/test-priors.R
@@ -11,13 +11,23 @@ test_prior <- function(prior, skip_moments = FALSE){
set.seed(1)
# tests rng and print function (for plot)
samples <- rng(prior, 100000)
- hist(samples, main = print(prior, plot = T), breaks = 50, freq = FALSE)
+ if(is.prior.discrete(prior)){
+ barplot(table(samples)/length(samples), main = print(prior, plot = T), width = 1/(max(samples)+1), space = 0, xlim = c(-0.25, max(samples)+0.25))
+ }else if(is.prior.spike_and_slab(prior)){
+ xh <- hist(samples[samples != 0], breaks = 50, plot = FALSE)
+ xh$density <- xh$density * mean(samples != 0)
+ plot(xh, main = print(prior, plot = T), freq = FALSE)
+ }else{
+ hist(samples, main = print(prior, plot = T), breaks = 50, freq = FALSE)
+ }
# tests density function
lines(prior, individual = TRUE)
# tests quantile function
- abline(v = quant(prior, 0.5), col = "blue", lwd = 2)
+ if(!is.prior.spike_and_slab(prior)){
+ abline(v = quant(prior, 0.5), col = "blue", lwd = 2)
+ }
# tests that pdf(q(x)) == x
- if(prior$distribution != "point"){
+ if(!is.prior.point(prior) & !is.prior.discrete(prior) & !is.prior.spike_and_slab(prior)){
expect_equal(.25, cdf(prior, quant(prior, 0.25)), tolerance = 1e-5)
}
# test mean and sd functions
@@ -130,6 +140,13 @@ test_that("Beta prior distribution works", {
})
+test_that("Beta prior distribution works", {
+
+ expect_doppelganger("prior-bernoulli-1", function()test_prior(prior("bernoulli", list(.50))))
+ expect_doppelganger("prior-bernoulli-2", function()test_prior(prior("bernoulli", list(.66))))
+
+})
+
test_that("Uniform prior distribution works", {
expect_doppelganger("prior-uniform-1", function()test_prior(prior("uniform", list(0, 1))))
@@ -143,6 +160,16 @@ test_that("Spike prior distribution works", {
})
+test_that("spike and slab prior distribution works", {
+
+ expect_doppelganger("prior-spike-and-slab-1", function()test_prior(
+ prior_spike_and_slab(
+ prior("gamma", list(2, 2), list(0, Inf)),
+ prior_inclusion = prior("beta", list(2, 1)))
+ ))
+
+})
+
test_that("PET & PEESE prior distribution works", {
expect_doppelganger("prior-PET-1", function()test_prior(prior_PET("normal", list(0, 1))))
diff --git a/tests/testthat/test-summary-tables.R b/tests/testthat/test-summary-tables.R
index 3bbd7ce..bbac81f 100644
--- a/tests/testthat/test-summary-tables.R
+++ b/tests/testthat/test-summary-tables.R
@@ -99,7 +99,7 @@ test_that("Summary tables functions work",{
### test additional settings
# transformations
runjags_summary2t <- runjags_estimates_table(fit2, transformations = list("m" = list(fun = exp)))
- expect_equal(exp(models[[3]]$fit_summary[1,c("Mean","lCI","Median","uCI","MCMC_error")]), runjags_summary2t[1,c("Mean","lCI","Median","uCI","MCMC_error")], tolerance = 1e-5)
+ expect_equal(exp(models[[3]]$fit_summary[1,c("lCI","Median","uCI","MCMC_error")]), runjags_summary2t[1,c("lCI","Median","uCI","MCMC_error")], tolerance = 1e-5)
expect_equal(colnames(models[[3]]$fit_summary), colnames(runjags_summary2t))
expect_equal(rownames(models[[3]]$fit_summary), rownames(runjags_summary2t))
@@ -527,3 +527,101 @@ test_that("Summary tables functions work (formulas + factors)",{
})
+test_that("Summary tables functions work (spike and slab priors)",{
+
+ skip_on_os(c("mac", "linux", "solaris")) # multivariate sampling does not exactly match across OSes
+ skip_on_cran()
+
+ set.seed(1)
+
+ data_formula <- data.frame(
+ x_cont1 = rnorm(60),
+ x_fac2t = factor(rep(c("A", "B"), 30), levels = c("A", "B")),
+ x_fac3o = factor(rep(c("A", "B", "C"), 20), levels = c("A", "B", "C")),
+ x_fac3t = factor(rep(c("A", "B", "C"), 20), levels = c("A", "B", "C"))
+ )
+ data <- list(
+ y = rnorm(60, .4 * data_formula$x_cont1 + ifelse(data_formula$x_fac3t == "A", 0.0, ifelse(data_formula$x_fac3t == "B", -0.2, 0.4)), 1),
+ N = 60
+ )
+
+ # create model with mix of a formula and free parameters ---
+ formula_list0 <- list(mu = ~ x_cont1 + x_fac2t + x_fac3o)
+
+ formula_prior_list0 <- list(
+ mu = list(
+ "intercept" = prior("normal", list(0, 5)),
+ "x_cont1" = prior_spike_and_slab(
+ prior_parameter = prior("normal", list(0, 0.5)),
+ prior_inclusion = prior("beta", list(1, 1))
+ ),
+ "x_fac2t" = prior_spike_and_slab(
+ prior_parameter = prior_factor("normal", contrast = "treatment", list(0, 1)),
+ prior_inclusion = prior("beta", list(1, 1))
+ ),
+ "x_fac3o" = prior_spike_and_slab(
+ prior_parameter = prior_factor("mnormal", contrast = "orthonormal", list(0, 1)),
+ prior_inclusion = prior("spike", list(.5))
+ )
+ )
+ )
+
+
+ prior_list <- list(sigma = prior("lognormal", list(0, 1)))
+ formula_data_list <- list(mu = data_formula)
+
+ model_syntax <- paste0(
+ "model{\n",
+ "for(i in 1:N){\n",
+ " y[i] ~ dnorm(mu[i], 1/pow(sigma, 2))\n",
+ "}\n",
+ "}"
+ )
+
+ log_posterior <- function(parameters, data){
+ sum(stats::dnorm(data$y, parameters[["mu"]], parameters[["sigma"]], log = TRUE))
+ }
+
+ fit0 <- JAGS_fit(
+ model_syntax = model_syntax, data = data, prior_list = prior_list,
+ formula_list = formula_list0, formula_data_list = formula_data_list, formula_prior_list = formula_prior_list0, seed = 1)
+
+ marglik0 <- list(logml = 0)
+ class(marglik0) <- "bridge"
+ # mix posteriors
+ models <- list(
+ list(fit = fit0, marglik = marglik0, fit_summary = runjags_estimates_table(fit0), prior_weights = 1)
+ )
+ models <- models_inference(models)
+
+ model_summary_table(models[[1]])
+ ### checking summary functions
+ # model summary
+ model_summary <- model_summary_table(models[[1]])
+ expect_equal(model_summary[,1], c("Model ", "Prior prob. ", "log(marglik) ", "Post. prob. ", "Inclusion BF ", " "))
+ expect_equal(model_summary[,4], c("Parameter prior distributions", "(mu) intercept ~ Normal(0, 5)", "(mu) x_cont1 ~ Normal(0, 0.5) * Beta(1, 1)", "(mu) x_fac2t ~ treatment contrast: Normal(0, 1) * Beta(1, 1)", "(mu) x_fac3o ~ orthonormal contrast: mNormal(0, 1) * Spike(0.5)", "sigma ~ Lognormal(0, 1)"))
+
+ model_estimates <- runjags_estimates_table(fit0)
+ expect_equal(colnames(model_estimates), c("Mean", "SD", "lCI", "Median", "uCI", "MCMC_error", "MCMC_SD_error", "ESS", "R_hat"))
+ expect_equal(rownames(model_estimates), c("(mu) intercept", "(mu) x_cont1", "(mu) x_cont1 (inclusion)", "(mu) x_fac2t[B]", "(mu) x_fac2t (inclusion)", "(mu) x_fac3o[1]", "(mu) x_fac3o[2]", "(mu) x_fac3o (inclusion)", "sigma"))
+ expect_equal(unname(unlist(model_estimates[7,])), c(-1.776437e-03, 4.269207e-02, 0.000000e+00, 0.000000e+00, 0.000000e+00, 3.428435e-04, 8.000000e-03, 1.550600e+04, 1.001065e+00), tolerance = 1e-3)
+
+ model_estimates <- runjags_estimates_table(fit0, transform_orthonormal = TRUE, conditional = TRUE)
+ expect_equal(colnames(model_estimates), c("Mean", "SD", "lCI", "Median", "uCI", "MCMC_error", "MCMC_SD_error", "ESS", "R_hat"))
+ expect_equal(rownames(model_estimates), c("(mu) intercept", "(mu) x_cont1", "(mu) x_cont1 (inclusion)", "(mu) x_fac2t[B]", "(mu) x_fac2t (inclusion)", "(mu) x_fac3o [dif: A]", "(mu) x_fac3o [dif: B]", "(mu) x_fac3o [dif: C]", "(mu) x_fac3o (inclusion)", "sigma"))
+ expect_equal(unname(unlist(model_estimates[8,])), c(0.0626582174, 0.1661778973, -0.2621073424, 0.0591205499, 0.3954805352, NA, NA, NA, NA), tolerance = 1e-3)
+
+ model_estimates <- runjags_estimates_table(fit0, transform_orthonormal = TRUE, conditional = TRUE, remove_inclusion = TRUE)
+ expect_equal(colnames(model_estimates), c("Mean", "SD", "lCI", "Median", "uCI", "MCMC_error", "MCMC_SD_error", "ESS", "R_hat"))
+ expect_equal(rownames(model_estimates), c("(mu) intercept", "(mu) x_cont1", "(mu) x_fac2t[B]", "(mu) x_fac3o [dif: A]", "(mu) x_fac3o [dif: B]", "(mu) x_fac3o [dif: C]", "sigma"))
+ expect_equal(unname(unlist(model_estimates[2,])), c(3.040927e-01, 1.355633e-01, 3.256895e-02, 3.058346e-01, 5.678668e-01, 2.298187e-03, 1.300000e-02, 5.720000e+03, 1.000421e+00), tolerance = 1e-3)
+
+ model_inference <- runjags_inference_table(fit0)
+ expect_equal(colnames(model_inference), c("Parameter", "prior_prob", "post_prob", "inclusion_BF"))
+ expect_equal(model_inference[,1], c("(mu) x_cont1", "(mu) x_fac2t", "(mu) x_fac3o"))
+ expect_equal(model_inference[,2], c(0.5, 0.5, 0.5))
+ expect_equal(model_inference[,3], c(0.7798125, 0.1864375, 0.0399375), tolerance = 1e-3)
+ expect_equal(model_inference[,4], c(3.54158388, 0.22916187, 0.04159885), tolerance = 1e-3)
+
+})
+
diff --git a/vignettes/ComparisonR.Rmd b/vignettes/ComparisonR.Rmd
index 006b030..1bdae28 100644
--- a/vignettes/ComparisonR.Rmd
+++ b/vignettes/ComparisonR.Rmd
@@ -14,6 +14,7 @@ vignette: >
%\VignetteEngine{knitr::rmarkdown_notangle}
---
+
```{r setup, include = FALSE}
# is_check <- ("CheckExEnv" %in% search()) ||
# any(c("_R_CHECK_TIMINGS_", "_R_CHECK_LICENSE_") %in% names(Sys.getenv()))
@@ -41,7 +42,7 @@ where $\mu$ corresponds to the grand mean, $\sigma$ to the grand standard deviat
To perform the test, we specify two competing hypotheses:
- the null hypothesis $\mathcal{H}_0: \delta = 0$ assuming that the standardized effect size is zero,
- - and the alternative hypothesis $\mathcal{H}_1: \delta \sim g()$ assuming that the standardized effect size is non-zero and prior distribution $g()$ characterizes our believes about the possible values of the standardized effect size.
+ - and the alternative hypothesis $\mathcal{H}_1: \delta \sim g()$ assuming that the standardized effect size is non-zero and prior distribution $g()$ characterizes our hypothesis about the possible values of the standardized effect size.
### Kitchen Rolls Data Set
@@ -191,16 +192,16 @@ To obtain the resulting Bayes factor in favor of the null hypothesis, I can use
bridgesampling::bf(marglik_model_H0, marglik_model_Hp)
```
-or specify a BayesTools model ensemble object that provides that can be further inspected by a variety of functions.
+or specify a BayesTools model ensemble object that can be further inspected by a variety of functions.
```{r}
models_list <- models_inference(list(
list(model = ttest_model_H0, marglik = marglik_model_H0, prior_weights = 1/2),
list(model = ttest_model_Hp, marglik = marglik_model_Hp, prior_weights = 1/2)
))
-enseble_info <- ensemble_inference(models_list, parameters = "delta", is_null_list = list("delta" = c(TRUE, FALSE)))
+ensemble_info <- ensemble_inference(models_list, parameters = "delta", is_null_list = list("delta" = c(TRUE, FALSE)))
-ensemble_inference_table(enseble_info, "delta", BF01 = TRUE)
+ensemble_inference_table(ensemble_info, "delta", BF01 = TRUE)
```
The corresponding results can be also obtained numerically with the `ttestBF()` function from the BayesFactor R package [@BayesFactor].
diff --git a/vignettes/SpikeAndSlab.Rmd b/vignettes/SpikeAndSlab.Rmd
new file mode 100644
index 0000000..33b3d24
--- /dev/null
+++ b/vignettes/SpikeAndSlab.Rmd
@@ -0,0 +1,271 @@
+---
+title: "Bayes factors via spike and slab prior vs. bridge sampling"
+author: "František Bartoš"
+date: "`r Sys.Date()`"
+output:
+ rmarkdown::html_vignette:
+ self_contained: yes
+bibliography: ../inst/REFERENCES.bib
+csl: ../inst/apa.csl
+vignette: >
+ %\VignetteIndexEntry{Bayes factors via spike and slab prior vs. bridge sampling}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+ %\VignetteEngine{knitr::rmarkdown_notangle}
+---
+
+
+```{r setup, include = FALSE}
+# is_check <- ("CheckExEnv" %in% search()) ||
+# any(c("_R_CHECK_TIMINGS_", "_R_CHECK_LICENSE_") %in% names(Sys.getenv()))
+is_check <- F
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>",
+ eval = !is_check,
+ dev = "png"
+)
+if(.Platform$OS.type == "windows"){
+ knitr::opts_chunk$set(dev.args = list(type = "cairo"))
+}
+```
+
+
+One of the main features of BayesTools is assistance in generating JAGS [@JAGS] code based on formulas and prior distribution objects and subsequent estimation of marginal likelihoods with the bridgesampling R package [@bridgesampling]. Marginal likelihoods, $p(\text{data} \mid \mathcal{M})$, are the key ingredient for computing Bayes factors,
+
+$\text{BF}_{10} = \frac{p(\text{data} \mid \mathcal{M}_{1})}{p(\text{data} \mid \mathcal{M}_{0})}$,
+
+which quantify relative predictive performance of two competing models [@wrinch1921on; @kass1995bayes; @rouder2019teaching]. Convenient model specification then allows users and package developers to easily compute Bayes factors and test a wide range of informed hypotheses. See RoBMA [@RoBMA] and [@RoBSA] R packages for implementation examples.
+
+However, when considering a simple regression, the model space of all possible models increases exponentially with additional predictors. I.e., the possibility of including vs. excluding $k$ predictors leads to $2^k$ possible submodels that need to be computed. Even relatively computationally simple models (e.g. ~ 1 min of computation) with 10 possible predictors would result in more than 17 hours of computation. Therefore, we might require more computationally efficient methods when performing variable selection with more than a few covariates. In this vignette, I showcase how to use BayesTools package to specify spike and slab priors that aim to explore most of the model space and obtain posterior inclusion probabilities for each predictor within a single MCMC run [@kuo1998variable; @ohara2009review].
+
+
+```{r}
+library(BayesTools)
+```
+
+## Simulated Data
+
+To keep it simple, let's consider a linear regression with one continuous predictor $x$. We simulate $N = 100$ observations of a dependent variable $y$ under the presence of a small effect $\beta$ of a continuous predictor $x$.
+
+```{r}
+set.seed(-68) # set seed for reproducibility
+
+N <- 100 # number of observations
+x <- rnorm(N) # continuous predictor
+alpha <- -0.5 # intercept
+beta <- 0.15 # small effect
+
+# compute the mean parameter for each predictor value
+mu <- alpha + beta * x
+
+# generate the response for each observation
+y <- rnorm(N, mean = mu, sd = 1)
+```
+
+We quickly verify that our simulated data correspond to the desired settings (up to a random error) with the `lm` function.
+
+```{r}
+summary(lm(y ~ x))
+```
+
+
+## Model Specification
+
+We consider two following models:
+
+ - $\mathcal{M}_0$: $\beta = 0$
+ - and $\mathcal{M}_1$: $\beta \sim g()$,
+
+where $g()$ characterizes our hypothesis about the degree of the effect. In our example, we specify a simple two-sided hypothesis represented by a normal distribution with mean 0 and standard deviation 0.5, e.g., $\beta \sim \text{Normal}(0, 0.5^2)$.
+
+
+## Maginal Likelihoods
+
+First, we compute the Bayes factor model comparison via marginal likelihoods. To do that, we need to specify the likelihood for the response variable $y$,
+
+```{r}
+model_likelihood <-
+"model{
+ for(i in 1:N){
+ y[i] ~ dnorm(mu[i], pow(sigma, -2))
+ }
+}
+"
+```
+
+where `mu` corresponds to the mean parameter (that we specify via a formula in the next step) and `sigma` to a standard deviation of the response variable (that we treat as a nuisance parameter here).
+
+We specify formulas for the `mu` parameter of each of the considered models,
+
+```{r}
+formula_M0 <- list("mu" = ~ 1)
+formula_M1 <- list("mu" = ~ 1 + x)
+```
+
+where `1` corresponds to the intercept (it is not necessary for the second model as it is included by default).
+
+To finish the model specification, we set the prior distribution corresponding to our hypothesis test of the beta parameter, set a broad prior distributions for the nuisance intercept and sigma parameters, and create a list containing data for the model specified within `model_likelihood` (in the first step) and a data frame for the data contained within the formula for mu within `formula_M0` and `formula_M1` (specified in the second step).
+
+```{r}
+# prior on the test parameter
+prior_beta <- prior(distribution = "normal", parameters = list(mean = 0, sd = 0.5))
+
+# priors on the nuisance parameters
+prior_int <- prior(distribution = "normal", parameters = list(mean = 0, sd = 5))
+prior_sigma <- prior(distribution = "normal", parameters = list(mean = 0, sd = 5), truncation = list(0, Inf))
+
+# the data list
+data_list <- list(
+ y = y,
+ N = N
+)
+data_formula <- data.frame(
+ x = x
+)
+```
+
+We estimate the models with the `JAGS_fit` function. Since we are using the formula interface (which allows us to specify multiple formulas for different parameters), we need to pass the arguments as named lists,
+
+```{r}
+M0 <- JAGS_fit(
+ # specification for the `model_likelihood` part
+ model_syntax = model_likelihood,
+ data = list(y = y, N = N),
+ prior_list = list("sigma" = prior_sigma),
+
+ # specification for the `formula_M0` part
+ formula_list = formula_M0,
+ formula_prior_list = list("mu" = list("intercept" = prior_int)),
+ formula_data_list = list("mu" = data_formula),
+
+ # seed for reproducibility
+ seed = 0
+)
+
+M1 <- JAGS_fit(
+ model_syntax = model_likelihood,
+ data = list(y = y, N = N),
+ prior_list = list("sigma" = prior_sigma),
+ formula_list = formula_M1,
+ formula_prior_list = list("mu" = list("intercept" = prior_int, "x" = prior_beta)),
+ formula_data_list = list("mu" = data_formula),
+ seed = 1
+)
+```
+
+We quickly verify that our parameter estimates (from the full model) are similar to the frequentist results obtained via `lm` function earlier.
+
+```{r}
+JAGS_estimates_table(M1)
+```
+
+To obtain the marginal likelihoods and compute Bayes factors, we only need to write the likelihood function corresponding to the JAGS model. Importantly, BayesTools handles all priors and formula related computation automatically, in other words, we do not need to worry about computing the mean parameter based on the intercept and predictors since we already obtain the computed mu in the `parameters[["mu"]]` object (a vector with a value for each y),
+
+```{r}
+log_posterior <- function(parameters, data){
+ sum(dnorm(
+ x = data[["y"]],
+ mean = parameters[["mu"]],
+ sd = parameters[["sigma"]],
+ log = TRUE
+ ))
+}
+```
+
+where the `parameters` arguments is a list containing the parameters and `data` argument is a list containing data. We use `sum(dnorm(..., log = TRUE))` to sum the logarithmic likelihood of all observations.
+
+Finally, we pass our objects to the `JAGS_bridgesampling` function to compute the marginal likelihoods.
+
+```{r}
+marglik_model_H0 <- JAGS_bridgesampling(
+ # specification for the model part
+ fit = M0,
+ log_posterior = log_posterior,
+ data = list(y = y, N = N),
+ prior_list = list("sigma" = prior_sigma),
+
+ # specification for the formula` part
+ formula_list = formula_M0,
+ formula_prior_list = list("mu" = list("intercept" = prior_int)),
+ formula_data_list = list("mu" = data_formula)
+)
+
+marglik_model_H1 <- JAGS_bridgesampling(
+ fit = M1,
+ log_posterior = log_posterior,
+ data = list(y = y, N = N),
+ prior_list = list("sigma" = prior_sigma),
+ formula_list = formula_M1,
+ formula_prior_list = list("mu" = list("intercept" = prior_int, "x" = prior_beta)),
+ formula_data_list = list("mu" = data_formula),
+)
+```
+
+We specify a BayesTools model ensemble object that we interrogate with the `ensemble_inference_table` function for information about the test for the beta parameter.
+
+```{r}
+models_list <- models_inference(list(
+ list(model = M0, marglik = marglik_model_H0, prior_weights = 1/2),
+ list(model = M1, marglik = marglik_model_H1, prior_weights = 1/2)
+))
+ensemble_info <- ensemble_inference(models_list, parameters = "x", is_null_list = list("x" = c(TRUE, FALSE)))
+
+ensemble_inference_table(ensemble_info, parameters = "x")
+```
+
+We find absence of evidence for either of the hypotheses, $\text{BF}_{10} = 1.181$, with posterior probability of $P(\mathcal{M}_{1} \mid \text{data}) = 0.542$ (asuming equal prior probability specified via `prior_weights` in the `models_inference` function previously).
+
+
+## Spike and Slab Priors
+
+The @kuo1998variable's spike an slab prior distribution is specified as a mixture of two prior distributions. A spike, a parameter value of zero corresponding to no effect, and a slab, a parameter value sampled from a continuous density corresponding to the alternative hypothesis. The parametrization uses two independent prior distributions: one for the parameter value, $\beta \sim g()$, and one for the inclusion indicator, $I_\beta \sim f()$, which assigns the prior model probability $P(\mathcal{M}_{1})$ of inclusion.
+
+The inclusion indicator can attain one of two values: either zero or one. Multiplying the parameter with the inclusion indicator, $\beta I_\beta$, then results in setting the parameter to zero when the indicator is zero, or keeping its original value when the indicator is one. The proportion of times the indicator attains the value of one then corresponds to the posterior inclusion probability of the predictor, $P(\mathcal{M}_{1} \mid \text{data})$. Since Bayes factor can be written as the change from prior to posterior odds,
+
+$\text{BF}_{10} = \frac{p(\mathcal{M}_{1} \mid \text{data})}{p(\mathcal{M}_{0} \mid \text{data})} / \frac{p(\mathcal{M}_{1})}{p(\mathcal{M}_{0})}$,
+
+we can also estimate the Bayes factor via the inclusion indicator.
+
+
+Now, we compare the two models using the spike and slab prior. We have already specified the likelihood, data lists, prior distributions for the nuisance parameters, and even the formulas (now we need only formula for the full model) in the previous sections. Therefore, we proceed directly by specifying the spike and slab prior distribution for the predictor. We use the `prior_spike_and_slab` which follows similar notation as the `prior` function. We need to specify the distribution, the parameters (and we could also set truncation if needed) corresponding to the alternative hypothesis. Furthermore, we need to specify the prior distribution for the inclusion via the `prior_inclusion` argument. Here, we use a $\text{Spike}(0.5)$ prior which sets the prior inclusion probability to 1/2.
+
+```{r}
+prior_beta_spike_and_slab <- prior_spike_and_slab(
+ prior_parameter = prior(distribution = "normal", parameters = list(mean = 0, sd = 0.5)),
+ prior_inclusion = prior(distribution = "spike", parameters = list(location = 0.5))
+)
+
+```
+
+Then we can directly proceed to calling the `JAGS_fit` function with the same specification as we used for the `M1` model, however, changing the prior distribution object for the predictor `x` to the newly created `prior_beta_spike_and_slab` prior distribution.
+
+```{r}
+MS <- JAGS_fit(
+ model_syntax = model_likelihood,
+ data = list(y = y, N = N),
+ prior_list = list("sigma" = prior_sigma),
+ formula_list = formula_M1,
+ formula_prior_list = list("mu" = list("intercept" = prior_int, "x" = prior_beta_spike_and_slab)),
+ formula_data_list = list("mu" = data_formula),
+ seed = 1
+)
+```
+
+We can again verify that our parameter estimates match the previous results. Now, we need to set the `conditional = TRUE` argument in the `JAGS_estimates_table` to obtain samples assuming that the spike and slab parameter values are included. (The function summarized the complete posterior distribution by default, i.e., parameter estimates model-averaged across the spike and the slab.)
+
+```{r}
+JAGS_estimates_table(MS, conditional = TRUE)
+```
+
+The estimates are essentially identical to the estimates from the previous models. Finally, we can also obtain summary of the inclusion probabilities via the `JAGS_inference_table` function
+
+```{r}
+JAGS_inference_table(MS)
+```
+
+As before, we find absence of evidence for either of the hypotheses, $\text{BF}_{10} = 1.157$, with posterior probability of $P(\mathcal{M}_{1} \mid \text{data}) = 0.536$.
+
+
+### References
+