# FBMS-Flexible Bayesian Model Selection and Model Averaging

This vignette introduces the FBMS package and shows how to perform Flexible Bayesian Model Selection (BMS) and Bayesian Model Averaging (BMA) for linear, generalized, nonlinear, fractional–polynomial, mixed–effect, and logic–regression models. More details are provided in the paper: “FBMS: An R Package for Flexible Bayesian Model Selection and Model Averaging” available on arxiv: more explicit examples accompanying the paper can be found on github https://github.com/jonlachmann/GMJMCMC/tree/data-inputs/tests_current

# Setup

Load technical markdown settings and a custom function for printing only what is needed through `printn()`.

# Bayesian Model Selection and Averaging

1. Consider a class of models: \(\Omega: m_1(Y|X,\theta_1), \dots, m_k(Y|X,\theta_k)\).
2. Assign priors to models \(P(m_1), \dots, P(m_k)\) and parameters \(P(\theta_j|m_j)\).
3. Obtain joint posterior \(P(m_j,\theta_j|D)\).
4. Make inference on a quantity of interest \(\Delta\).

Show equations

\[ P(\Delta|D) = \sum_{m\in\Omega} P(m|D)\, \int_{\Theta_m} P(\Delta|m,\theta,D)\, P(\theta|m,D)\, d\theta. \]

## Bayesian Generalized Nonlinear Model (BGNLM)

Reference: [@hubin2020flexible]

Show equations

\[ \begin{aligned} Y_i \mid \mu_i, \phi &\sim \mathfrak{f}(y \mid \mu_i,\phi), \quad i = 1,\dots,n,\\ \mathsf{h}(\mu_i) &= \beta_0 + \sum_{j=1}^{q} \gamma_j \beta_j F_j(\mathbf{x}_i, \boldsymbol{\alpha}_j) + \sum_{k=1}^{r} \delta_k. \end{aligned} \]

- Predictors (features): \(F_j(\mathbf{x}_i, \boldsymbol{\alpha}_j)\), \(j=1,\dots,q\)
- Random effects: \(\delta_k\)
- Total models: \(2^{q+r}\)

### Feature space

Depending on allowed nonlinear functions \(\mathbb{G}\): neural nets (sigmoid), decision trees (thresholds), MARS (hinges), fractional polynomials, logic regression, etc.

### Transformations available in `FBMS`

Name | Function | Name | Function
--- | --- | --- | ---
sigmoid | \(1 / (1 + \exp(-x))\) | sqroot | \(|x|^{1/2}\)
relu | \(\max(x, 0)\) | troot | \(|x|^{1/3}\)
nrelu | \(\max(-x, 0)\) | sin_deg | \(\sin(x/180 \cdot \pi)\)
hs | \(x > 0\) | cos_deg | \(\cos(x/180 \cdot \pi)\)
nhs | \(x < 0\) | exp_dbl | \(\exp(-|x|)\)
gelu | \(x \Phi(x)\) | gauss | \(\exp(-x^2)\)
ngelu | \(-x \Phi(-x)\) | erf | \(2 \Phi(\sqrt{2}x) - 1\)
pm2 | \(x^{-2}\) | p0pm2 | \(p0(x) \cdot x^{-2}\)
pm1 | \(\text{sign}(x) |x|^{-1}\) | p0pm05 | \(p0(x) \cdot |x|^{-0.5}\)
pm05 | \(|x|^{-0.5}\) | p0p0 | \(p0(x)^2\)
p05 | \(|x|^{0.5}\) | p0p05 | \(p0(x) \cdot |x|^{0.5}\)
p2 | \(x^2\) | p0p1 | \(p0(x) \cdot x\)
p3 | \(x^3\) | p0p2 | \(p0(x) \cdot x^2\)
|  | p0p3 | \(p0(x) \cdot x^3\)

Custom functions can be added to \(\mathbb{G}\).

## Priors

### Model priors

Let \(M = (\gamma_1, \dots, \gamma_q)\) (for linear models \(q=p\)).

Penalizing complexity priors

Show equation

\[ P(M) \propto \mathbb{I}(|\boldsymbol\gamma_{1:q}| \le Q) \prod_{j=1}^q r^{\gamma_j c(F_j(\mathbf{x}, \boldsymbol{\alpha}_j))}, \quad 0 < r < 1, \]

- \(c(F_j(\cdot))\): complexity measure (linear models: \(c(x)=1\); BGNLM counts algebraic operators).

The following parameter will be used to change the prior penalization.

    # model prior complexity penalty
    model_prior = list(r = 0.005) #default is 1/n

### Parameter priors

Mixtures of g-priors

Show equations

\[ \begin{aligned} P(\beta_0, \phi \mid M) &\propto \phi^{-1}, \\ P(\boldsymbol{\beta} \mid g) &\sim \mathbb{N}_{|M|}\!\left(\mathbf{0},\, g \cdot \phi\, J_n(\boldsymbol{\beta})^{-1}\right), \\ \frac{1}{1+g} &\sim \text{tCCH}\!\left(\frac{a}{2}, \frac{b}{2}, \rho, \frac{s}{2}, v, \kappa\right). \end{aligned} \]

- Robust-g (convenient default): \(a=1, b=2, \rho=1.5, s=0, v=\frac{n+1}{|M|+1}, \kappa=1\).

Jeffreys prior

Show equations

\[ P(\phi \mid M) = \phi^{-1}, \quad P(\beta_0, \boldsymbol{\beta} \mid M) = |J_n(\beta_0, \boldsymbol{\beta})|^{1/2}. \]

All priors from the table below (default is the g-prior)

### Parameter priors available (and where to tune)

Prior (Alias) | \(a\) | \(b\) | \(\rho\) | \(s\) | \(v\) | \(k\) | Families
--- | --- | --- | --- | --- | --- | --- | ---
Default: |  |  |  |  |  |  |
`g-prior` | \(g\) (default: \(\max(n, p^2)\)) |  |  |  |  |  | GLM
tCCH-Related Priors: |  |  |  |  |  |  |
`CH` | \(a\) | \(b\) | 0 | \(s\) | 1 | 1 | GLM
`uniform` | 2 | 2 | 0 | 0 | 1 | 1 | GLM
`Jeffreys` | 0 | 2 | 0 | 0 | 1 | 1 | GLM
`beta.prime` | \(1/2\) | \(n - p_M - 1.5\) | 0 | 0 | 1 | 1 | GLM
`benchmark` | 0.02 | \(0.02 \max(n, p^2)\) | 0 | 0 | 1 | 1 | GLM
`TG` | \(2a\) | 2 | 0 | \(2s\) | 1 | 1 | GLM
`ZS-adapted` | 1 | 2 | 0 | \(n + 3\) | 1 | 1 | GLM
`robust` | 1 | 2 | 1.5 | 0 | \((n+1)/(p_M+1)\) | 1 | GLM
`hyper-g-n` | 1 | 2 | 1.5 | 0 | 1 | \(1/n\) | GLM
`intrinsic` | 1 | 1 | 1 | 0 | \((n + p_M + 1)/(p_M + 1)\) | \((n + p_M + 1)/n\) | GLM
`tCCH` | \(a\) | \(b\) | \(\rho\) | \(s\) | \(v\) | \(k\) | GLM
Other Priors: |  |  |  |  |  |  |
`EB-local` | \(a\) |  |  |  |  |  | GLM
`EB-global` | \(a\) |  |  |  |  |  | G
`JZS` | \(a\) |  |  |  |  |  | G
`ZS-null` | \(a\) |  |  |  |  |  | G
`ZS-full` | \(a\) |  |  |  |  |  | G
`hyper-g` | \(a\) |  |  |  |  |  | GLM
`hyper-g-laplace` | \(a\) |  |  |  |  |  | G
`AIC` | None |  |  |  |  |  | GLM
`BIC` | None |  |  |  |  |  | GLM
`Jeffreys-BIC` | Var |  |  |  |  |  | GLM

Here \(p_M\) is the number of predictors excluding the intercept. “G” denotes Gaussian-only; “GLM” additionally includes binomial, Poisson, and gamma.

How to switch priors in code (be explicit):

    # g-prior with g = 1000
    beta_prior = list(type = "g-prior", alpha = 1000)

    # Robust prior (tune by Table parameters)
    beta_prior = list(type = "robust")

    # Jeffreys-BIC
    beta_prior = list(type = "Jeffreys-BIC")

    # Generic tCCH (provide all hyperparameters explicitly)
    beta_prior = list(type = "tCCH", a = 2, b = 2, rho = 0, s = 0, v = 1, k = 1)

# Inference algorithms

### Model posterior

Show equations

Marginal likelihood \[ P(D|M) = \int_{\Theta_M} P(D|\theta_M, M) \, P(\theta_M|M) \, d\theta_M. \]

Posterior \[ P(M|D) = \frac{P(D|M) P(M)}{\sum_{M' \in \Omega} P(D|M') P(M')}. \]

Approximation over discovered models \(\Omega^*\) \[ P(M|D) \approx \frac{P(D|M) P(M)}{\sum_{M' \in \Omega^*} P(D|M') P(M')}, \quad M \in \Omega^*. \]

Marginal inclusion probability \[ P(\gamma_j=1|D) \approx \sum_{M \in \Omega^*: \gamma_j=1} P(M|D). \]

### MCMC, MJMCMC, GMJMCMC

- Variable selection spans exponential number of models; naive MCMC can get trapped.
- MJMCMC: random mode jumps + local improvements; valid MH acceptance [@hubin2018mode].
- GMJMCMC: embeds MJMCMC in a genetic algorithm to traverse huge spaces [@hubin2020logic; @hubin2020flexible].
- RGMJMCMC: reversible version [@hubin2021reversible].
- Subsampling: efficient for tall data [@lachmann2022subsampling].

### Parallelization

Run multiple chains and aggregate unique models \(\Omega^*\):

Show equation

\[ \widehat{P}(\Delta|D) = \sum_{M \in \Omega^*} P(\Delta|M,D) \, \widehat{P}(M|D). \]

    # Parameters for parallel runs are set to a single thread and single core to comply with CRAN requirenments (please tune for your machine if you have more capacity)
    runs  <- 1  # 1 set for simplicity; use rather 16 or more
    cores <- 1  # 1 set for simplicity; use rather 8 or more

# Example 1 — BGNLM: recovering Kepler’s third law

Data: \(n=939\) exoplanets; variables include `semimajoraxis`, `period`, `hoststar_mass`, etc.
Target relationship: \({a \approx K_2 \left(P^2 M_h\right)^{1/3}}\).

We shall run a single chain GMJMCMC

    # Load example
    data <- FBMS::exoplanet

    # Choose a small but expressive transform set for a quick demo
    transforms <- c("sigmoid", "sin_deg", "exp_dbl", "p0", "troot", "p3")

    # ---- fbms() call (simple GMJMCMC) ----
    # Key parameters (explicit):
    # - formula      : semimajoraxis ~ 1 + .       # response and all predictors
    # - data         : data                        # dataset
    # - beta_prior   : list(type = "g-prior")      # parameter prior
    # - model_prior  : list(r = 1/dim(data)[1])    # model prior
    # - method       : "gmjmcmc"                   # exploration strategy
    # - transforms   : transforms                  # nonlinear feature dictionary
    # - P            : population size per generation (search breadth)
    result_single <- fbms(
      formula     = semimajoraxis ~ 1 + .,
      data        = data,
      beta_prior  = list(type = "g-prior", alpha = dim(data)[1]),
      model_prior = list(r = 1/dim(data)[1]),
      method      = "gmjmcmc",
      transforms  = transforms,
      P           = 10
    )

    # Summarize
    printn(summary(result_single))

    ##
    ## Best   population: 7  log marginal posterior: 3001.64
    ##
    ##                            feats.strings marg.probs
    ## 1          troot((period*sigmoid(mass))) 1.00000000
    ## 2        ((period*hoststar_mass)*period) 1.00000000
    ## 3             p3((period*hoststar_mass)) 1.00000000
    ## 4          (troot(period)*troot(period)) 1.00000000
    ## 5                 (period*hoststar_mass) 1.00000000
    ## 6 (sigmoid(mass)*(period*hoststar_mass)) 1.00000000
    ## 7                          hoststar_mass 0.91282456
    ## 8                   hoststar_temperature 0.01243301

and a parallel GMJMCMC

    # ---- fbms() call (parallel GMJMCMC) ----
    # Key parameters (explicit):
    # - formula      : semimajoraxis ~ 1 + .       # response and all predictors
    # - data         : data                        # dataset
    # - beta_prior   : list(type = "g-prior")      # parameter prior
    # - model_prior  : list(r = 1/dim(data)[1])    # model prior
    # - method       : "gmjmcmc"                   # exploration strategy
    # - transforms   : transforms                  # nonlinear feature dictionary
    # - runs, cores  : parallelization controls
    # - P            : population size per generation (search breadth)
    result_parallel <- fbms(
      formula     = semimajoraxis ~ 1 + .,
      data        = data,
      beta_prior  = list(type = "g-prior", alpha = dim(data)[1]),
      model_prior = list(r = 1/dim(data)[1]),
      method      = "gmjmcmc.parallel",
      transforms  = transforms,
      runs        = runs,     # by default the rmd has runs =  1; increase for convergence
      cores       = cores,    # by default the rmd has cores = 1; increase for convergence
      P           = 10
    )

    # Summarize
    printn(summary(result_parallel))

    ##
    ## Best   population: 10  thread: 1  log marginal posterior: 2144.064
    ##
    ##                                              feats.strings marg.probs
    ## 1                                                   period  1.0000000
    ## 2                                            (period*mass)  1.0000000
    ## 3                            ((period*mass)*(period*mass))  1.0000000
    ## 4                                          sigmoid(period)  1.0000000
    ## 5                                               p0(period)  1.0000000
    ## 6                 ((period*mass)*(radius*(period*period)))  0.9999902
    ## 7                                 (radius*(period*period))  0.9999844
    ## 8 (hoststar_temperature*p0(((period*mass)*(period*mass))))  0.9945983
    ## 9                                             eccentricity  0.9554446

Plot output

    plot(result_parallel)

Convergence plots

    diagn_plot(result_parallel)

# Example 2 — Bayesian linear models

Reference: [@hubin2018mode]

Model

\[ \begin{aligned} Y_i \mid \mu_i, \phi &\sim \mathfrak{f}(y\mid \mu_i, \phi), \quad i=1,\dots,n, \\ \mathsf{h}(\mu_i) &= \beta_0 + \sum_{j=1}^{p} \gamma_j \beta_j x_{ij}. \end{aligned} \]

We simulate data with a known sparse truth and run MJMCMC with an explicit g-prior.

    library(mvtnorm)

    n <- 100               # sample size
    p <- 20                # number of covariates
    k <- 5                 # size of true submodel

    correct.model <- 1:k
    beta.k <- (1:5)/5

    beta <- rep(0, p)
    beta[correct.model] <- beta.k

    set.seed(123)
    x <- rmvnorm(n, rep(0, p))
    y <- x %*% beta + rnorm(n)

    # Standardize
    y <- scale(y)
    X <- scale(x) / sqrt(n)

    df <- as.data.frame(cbind(y, X))
    colnames(df) <- c("Y", paste0("X", seq_len(ncol(df) - 1)))

    printn(correct.model)

    ## [1] 1 2 3 4 5

    printn(beta.k)

    ## [1] 0.2 0.4 0.6 0.8 1.0

Run MJMCMC with a g-prior (g = 100)

    # ---- fbms() call (MJMCMC) ----
    # Explicit prior choice:
    #   beta_prior = list(type = "g-prior", alpha = 100)
    # To switch to another prior, e.g. robust:
    #   beta_prior = list(type = "robust")
    result.lin <- fbms(
      formula    = Y ~ 1 + .,
      data       = df,
      method     = "mjmcmc",
      N          = 5000,                         # number of iterations
      beta_prior = list(type = "g-prior", alpha = 100)
    )

Plot results

    plot(result.lin)

Summarize with posterior effects

    # 'effects' specifies quantiles for posterior modes of effects across models
    printn(summary(result.lin, effects = c(0.5, 0.025, 0.975)))

    ##
    ## Best log marginal posterior:  56.2259
    ##
    ## $PIP
    ##    feats.strings marg.probs
    ## 1             X4 1.00000000
    ## 2             X3 1.00000000
    ## 3             X5 1.00000000
    ## 4             X2 0.95262397
    ## 5             X1 0.16880319
    ## 6            X18 0.10464730
    ## 7            X20 0.10105529
    ## 8            X15 0.08767336
    ## 9            X10 0.08257391
    ## 10           X16 0.08200567
    ## 11           X13 0.07741378
    ## 12           X12 0.07474969
    ## 13           X11 0.07090813
    ## 14           X19 0.06903813
    ## 15            X7 0.06899870
    ## 16           X14 0.06586870
    ## 17            X9 0.05690408
    ## 18            X8 0.05592665
    ## 19            X6 0.05154572
    ## 20           X17 0.04739590
    ##
    ## $EFF
    ##    Covariate quant_0.5 quant_0.025 quant_0.975
    ## 1  intercept         0           0           0
    ## 2         X1         0           0      0.6606
    ## 3         X2    1.7024           0      1.7383
    ## 4         X3    4.8118      4.7487      4.8867
    ## 5         X4    4.9803      4.9086      5.0868
    ## 6         X5    4.7744      4.3659      4.8573
    ## 7         X6         0     -0.0694           0
    ## 8         X7         0     -0.1382           0
    ## 9         X8         0           0       0.072
    ## 10        X9         0     -0.1875           0
    ## 11       X10         0           0      0.5609
    ## 12       X11         0           0       0.444
    ## 13       X12         0     -0.2502           0
    ## 14       X13         0     -0.3655           0
    ## 15       X14         0           0      0.2716
    ## 16       X15         0           0      0.5654
    ## 17       X16         0           0      0.5074
    ## 18       X17         0           0      0.2266
    ## 19       X18         0           0      0.4858
    ## 20       X19         0     -0.2733           0
    ## 21       X20         0     -0.5012           0

Run parallel MJMCMC

    # ---- fbms() call (parallel MJMCMC) ----
    # Explicit prior choice:
    #   beta_prior = list(type = "g-prior", alpha = 100)
    # To switch to another prior, e.g. robust:
    #   beta_prior = list(type = "robust")
    #   method = mjmcmc.parallel
    #   runs, cores  : parallelization controls
    result.lin.par <- fbms(
      formula    = Y ~ 1 + .,
      data       = df,
      method     = "mjmcmc.parallel",
      N          = 5000,                         # number of iterations
      beta_prior = list(type = "g-prior", alpha = 100),
      runs = runs,
      cores = cores
    )
    printn(summary(result.lin.par, effects = c(0.5, 0.025, 0.975)))

    ##
    ## Best log marginal posterior:  56.2259
    ##
    ## $PIP
    ##    feats.strings marg.probs
    ## 1             X4 1.00000000
    ## 2             X3 1.00000000
    ## 3             X5 1.00000000
    ## 4             X2 0.95518016
    ## 5             X1 0.17518018
    ## 6            X10 0.10751407
    ## 7            X20 0.10044924
    ## 8            X18 0.09363710
    ## 9            X16 0.09089584
    ## 10           X15 0.08729800
    ## 11           X13 0.08418363
    ## 12           X19 0.07302402
    ## 13           X14 0.07262783
    ## 14           X11 0.06881624
    ## 15           X17 0.06104115
    ## 16            X9 0.05570830
    ## 17           X12 0.05215211
    ## 18            X6 0.04938128
    ## 19            X7 0.04711947
    ## 20            X8 0.04447923
    ##
    ## $EFF
    ##    Covariate quant_0.5 quant_0.025 quant_0.975
    ## 1  intercept         0           0           0
    ## 2         X1         0           0      0.6606
    ## 3         X2    1.7027           0      1.7394
    ## 4         X3    4.8118      4.7487      4.8885
    ## 5         X4    4.9803      4.9043      5.0888
    ## 6         X5    4.7744      4.3659      4.8576
    ## 7         X6         0     -0.0694           0
    ## 8         X7         0     -0.1346           0
    ## 9         X8         0           0       0.072
    ## 10        X9         0     -0.1875           0
    ## 11       X10         0           0      0.5625
    ## 12       X11         0           0      0.4382
    ## 13       X12         0     -0.2501           0
    ## 14       X13         0     -0.3742           0
    ## 15       X14         0           0      0.2739
    ## 16       X15         0           0      0.5567
    ## 17       X16         0           0      0.5074
    ## 18       X17         0           0      0.2266
    ## 19       X18         0           0      0.4718
    ## 20       X19         0     -0.2733           0
    ## 21       X20         0      -0.511           0

# Example 3 — Bayesian Fractional Polynomials (FP)

Reference: [@hubin2023fractional]

We augment the linear example to follow an FP truth and fit with GMJMCMC.

FP model class

\[ \begin{aligned} Y_i \mid \mu_i, \phi &\sim \mathfrak{f}(y|\mu_i,\phi),\\ \mathsf{h}(\mu_i) &= \beta_0 + \sum_{j=1}^{p} \sum_{k \in K} \gamma_{jk}\, \beta_{jk}\, \rho_k(x_{ij}), \quad \text{with } K = \mathbf{F}_0 \cup \mathbf{F}_1 \cup \mathbf{F}_2, \end{aligned} \]

    # Create FP-style response with known structure, covariates are from previous example
    df$Y <- p05(df$X1) + df$X1 + pm05(df$X3) + p0pm05(df$X3) + df$X4 +
             pm1(df$X5) + p0(df$X6) + df$X8 + df$X10 + rnorm(nrow(df))

    # Allow common FP transforms
    transforms <- c(
      "p0", "p2", "p3", "p05", "pm05", "pm1", "pm2", "p0p0",
      "p0p05", "p0p1", "p0p2", "p0p3", "p0p05", "p0pm05", "p0pm1", "p0pm2"
    )

    # Generation probabilities — here only modifications and mutations
    probs <- gen.probs.gmjmcmc(transforms)
    probs$gen <- c(0, 1, 0, 1)

    # Feature-generation parameters
    params <- gen.params.gmjmcmc(ncol(df) - 1)
    params$feat$D <- 1   # max depth 1 features

Run GMJMCMC (single-thread)

    result <- fbms(
      formula    = Y ~ 1 + .,
      data       = df,
      method     = "gmjmcmc",
      transforms = transforms,
      beta_prior = list(type = "Jeffreys-BIC"),
      probs      = probs,
      params     = params,
      P          = 10
    )

    printn(summary(result))

    ##
    ## Best   population: 10  log marginal posterior: -206.4266
    ##
    ##    feats.strings   marg.probs
    ## 1       pm05(X3) 1.0000000000
    ## 2       p0p0(X3) 1.0000000000
    ## 3        pm1(X5) 1.0000000000
    ## 4     p0pm05(X6) 0.9999959549
    ## 5      p0pm1(X3) 0.9359835862
    ## 6            X17 0.8981659238
    ## 7            X15 0.8795347468
    ## 8            X10 0.8584612830
    ## 9             X9 0.8507841966
    ## 10    p0pm2(X13) 0.7855348631
    ## 11           X20 0.1825368449
    ## 12            X4 0.0362927730
    ## 13            X7 0.0166240434
    ## 14           X13 0.0079441065
    ## 15            X2 0.0002276410
    ## 16     p0p3(X20) 0.0001784105

Parallel GMJMCMC

    result_parallel <- fbms(
      formula    = Y ~ 1 + .,
      data       = df,
      method     = "gmjmcmc.parallel",
      transforms = transforms,
      beta_prior = list(type = "Jeffreys-BIC"),
      probs      = probs,
      params     = params,
      P          = 10,
      runs       = runs,
      cores      = cores
    )

    printn(summary(result_parallel))

    ##
    ## Best   population: 10  thread: 1  log marginal posterior: -506.0695
    ##
    ##   feats.strings  marg.probs
    ## 1    p0pm05(X3) 1.000000000
    ## 2            X5 1.000000000
    ## 3       pm2(X5) 1.000000000
    ## 4            X3 0.398001125
    ## 5      p0p0(X3) 0.001258547

# Example 4 — Mixed-effects FP with interactions

Dataset: \(n=659\) kids; response \(y\): standardized height; covariates: `c.bf, c.age, m.ht, m.bmi, reg`. Random intercept for `dr.`
We specify a custom estimator that uses a mixed model (via `lme4`), and plug it into `fbms()` with `family = "custom"`. We pass extra parameters of the estimator through `mlpost_params = list(dr = dr,r = r)`

    # Custom approximate log marginal likelihood for mixed model using Laplace approximation
    mixed.model.loglik.lme4 <- function (y, x, model, complex, mlpost_params) {
      if (sum(model) > 1) {
        x.model <- x[, model]
        data <- data.frame(y, x = x.model[, -1], dr = mlpost_params$dr)
        mm <- lmer(as.formula(paste0("y ~ 1 +",
                              paste0(names(data)[2:(ncol(data)-1)], collapse = "+"),
                              " + (1 | dr)")), data = data, REML = FALSE)
      } else {
        data <- data.frame(y, dr = mlpost_params$dr)
        mm <- lmer(y ~ 1 + (1 | dr), data = data, REML = FALSE)
      }
      # log marginal likelihood (Laplace approx) + log model prior
      mloglik <- as.numeric(logLik(mm)) - 0.5 * log(length(y)) * (ncol(data) - 2)
      if (length(mlpost_params$r) == 0) mlpost_params$r <- 1 / nrow(x)
      lp <- log_prior(mlpost_params, complex)
      list(crit = mloglik + lp, coefs = fixef(mm))
    }

    library(lme4)

    ## Loading required package: Matrix

    data(Zambia, package = "cAIC4")

    df <- as.data.frame(sapply(Zambia[1:5], scale))

    transforms <- c("p0","p2","p3","p05","pm05","pm1","pm2",
                    "p0p0","p0p05","p0p1","p0p2","p0p3",
                    "p0p05","p0pm05","p0pm1","p0pm2")

    probs <- gen.probs.gmjmcmc(transforms)
    probs$gen <- c(1, 1, 0, 1) # include modifications and interactions

    params <- gen.params.gmjmcmc(ncol(df) - 1)
    params$feat$D <- 1
    params$feat$pop.max <- 10

    result2a <- fbms(
      formula      = z ~ 1 + .,
      data         = df,
      method       = "gmjmcmc.parallel",
      transforms   = transforms,
      probs        = probs,
      params       = params,
      P            = 10,
      N            = 100,
      runs         = runs,
      cores        = cores,
      family       = "custom",
      loglik.pi    = mixed.model.loglik.lme4,
      model_prior  = list(r = 1 / nrow(df)), # model_prior is passed to mlpost_params
      extra_params = list(dr = droplevels(Zambia$dr)) # extra_params are passed to mlpost_params
    )

    printn(summary(result2a, tol = 0.05, labels = names(df)[-1]))

    ##
    ## Best   population: 6  thread: 1  log marginal posterior: -876.5476
    ##
    ##   feats.strings marg.probs
    ## 1          c.bf  1.0000000
    ## 2         c.age  1.0000000
    ## 3         m.bmi  1.0000000
    ## 4          m.ht  1.0000000
    ## 5 (c.age*c.age)  0.9999999

# Example 5 — Bayesian Logic Regression

Reference: [@hubin2020logic]

Model

\[ \mathsf{h}(\mu_i) = \beta_0 + \sum_{j=1}^{q} \gamma_j \beta_j L_{ji}, \quad L_{ji} \text{ are logic trees (e.g., } (x_{i1}\wedge x_{i2}) \vee x_{i3}^c ). \]

We generate Boolean covariates and a known logic signal, define a custom estimator with a logic prior, and fit via GMJMCMC.

    n <- 2000
    p <- 50

    set.seed(1)
    X2 <- as.data.frame(matrix(rbinom(n * p, size = 1, prob = runif(n * p, 0, 1)), n, p))
    y2.Mean <- 1 + 7*(X2$V4*X2$V17*X2$V30*X2$V10) + 9*(X2$V7*X2$V20*X2$V12) +
               3.5*(X2$V9*X2$V2) + 1.5*(X2$V37)

    Y2 <- rnorm(n, mean = y2.Mean, sd = 1)
    df <- data.frame(Y2, X2)

    # Train/test split
    df.training <- df[1:(n/2), ]
    df.test     <- df[(n/2 + 1):n, ]
    df.test$Mean <- y2.Mean[(n/2 + 1):n]

Custom estimator with logic regression priors

    estimate.logic.lm <- function(y, x, model, complex, mlpost_params) {
      suppressWarnings({
        mod <- fastglm(as.matrix(x[, model]), y, family = gaussian())
      })
      mloglik <- -(mod$aic + (log(length(y)) - 2) * (mod$rank)) / 2
      wj <- complex$width
      lp <- sum(log(factorial(wj))) - sum(wj * log(4 * mlpost_params$p) - log(4))
      logpost <- mloglik + lp
      if (logpost == -Inf) logpost <- -10000
      list(crit = logpost, coefs = mod$coefficients)
    }

Run GMJMCMC

    set.seed(5001)

    # Only "not" operator; "or" is implied by De Morgan via "and" + "not"
    transforms <- c("not")
    probs <- gen.probs.gmjmcmc(transforms)
    probs$gen <- c(1, 1, 0, 1) # no projections

    params <- gen.params.gmjmcmc(p)
    params$feat$pop.max <- 50
    params$feat$L <- 15

    result <- fbms(
      formula     = Y2 ~ 1 + .,
      data        = df.training,
      method      = "gmjmcmc",
      transforms  = transforms,
      N           = 500,
      P           = 10,
      family      = "custom",
      loglik.pi   = estimate.logic.lm,
      probs       = probs,
      params      = params,
      model_prior = list(p = p)
    )

    printn(summary(result))

    ##
    ## Best   population: 9  log marginal posterior: -1912.786
    ##
    ##           feats.strings   marg.probs
    ## 1        ((V4*V10)*V17) 1.0000000000
    ## 2               (V9*V2) 1.0000000000
    ## 3                   V37 1.0000000000
    ## 4        ((V20*V7)*V12) 1.0000000000
    ## 5  (V33*(V30*(V4*V17))) 0.9999999179
    ## 6                   V30 0.9999969055
    ## 7                   V17 0.0137256351
    ## 8                   V32 0.0025183771
    ## 9                   V35 0.0020386609
    ## 10                  V12 0.0013621295
    ## 11                  V45 0.0012730261
    ## 12                  V18 0.0012697929
    ## 13                  V16 0.0010825818
    ## 14                  V41 0.0010734783
    ## 15                  V44 0.0009967761
    ## 16                  V31 0.0009461623
    ## 17                   V8 0.0008938150
    ## 18                  V13 0.0007713390
    ## 19                   V3 0.0006926010
    ## 20                  V24 0.0006547320
    ## 21                  V21 0.0006439233
    ## 22                  V49 0.0006393027
    ## 23                   V9 0.0006283463
    ## 24                  V36 0.0006209956
    ## 25                  V15 0.0006131099

    # Extract models
    mpm   <- get.mpm.model(result, y = df.training$Y2, x = df.training[,-1],
                           family = "custom", loglik.pi = estimate.logic.lm, params = list(p = 50))
    printn(mpm$coefs)

    ##          (Intercept)       ((V4*V10)*V17)              (V9*V2)
    ##            0.6578016            2.7049544            3.5767605
    ##                  V30                  V37 (V33*(V30*(V4*V17)))
    ##            0.6051133            1.3582049            2.1279042
    ##       ((V20*V7)*V12)
    ##            9.1213744

    mpm2  <- get.mpm.model(result, y = df.training$Y2, x = df.training[,-1])
    printn(mpm2$coefs)

    ##          (Intercept)       ((V4*V10)*V17)              (V9*V2)
    ##            0.6578016            2.7049544            3.5767605
    ##                  V30                  V37 (V33*(V30*(V4*V17)))
    ##            0.6051133            1.3582049            2.1279042
    ##       ((V20*V7)*V12)
    ##            9.1213744

    mbest <- get.best.model(result)
    printn(mbest$coefs)

    ##            Intercept       ((V4*V10)*V17)              (V9*V2)
    ##            0.6578016            2.7049544            3.5767605
    ##                  V30                  V37 (V33*(V30*(V4*V17)))
    ##            0.6051133            1.3582049            2.1279042
    ##       ((V20*V7)*V12)
    ##            9.1213744

Prediction

    # Correct link is identity for Gaussian
    pred       <- predict(result, x = df.test[,-1], link = function(x) x)
    pred_mpm   <- predict(mpm,   x = df.test[,-1], link = function(x) x)
    pred_best  <- predict(mbest, x = df.test[,-1], link = function(x) x)

    # RMSEs
    printn(sqrt(mean((pred$aggr$mean - df.test$Y2)^2)))

    ## [1] 1.528321

    printn(sqrt(mean((pred_mpm     - df.test$Y2)^2)))

    ## [1] 1.528355

    printn(sqrt(mean((pred_best    - df.test$Y2)^2)))

    ## [1] 1.528355

    printn(sqrt(mean((df.test$Mean - df.test$Y2)^2)))

    ## [1] 1.028222

    # Errors to the true mean (oracle)
    printn(sqrt(mean((pred$aggr$mean - df.test$Mean)^2)))

    ## [1] 1.076062

    printn(sqrt(mean((pred_best      - df.test$Mean)^2)))

    ## [1] 1.076113

    printn(sqrt(mean((pred_mpm       - df.test$Mean)^2)))

    ## [1] 1.076113

    # Quick diagnostic plot
    plot(pred$aggr$mean, df.test$Y2,
         xlab = "Predicted (BMA)", ylab = "Observed")
    points(pred$aggr$mean, df.test$Mean, col = 2)
    points(pred_best,      df.test$Mean, col = 3)
    points(pred_mpm,       df.test$Mean, col = 4)

# Example 6 — Full BGNLM classification (Bernoulli): `spam` data

We fit a binomial BGNLM and compare BMA, best-model, and MPM predictions.
Important: specify the correct link in `predict()` (here logistic).

    library(kernlab)
    data("spam")

    df <- spam[, c(58, 1:57)]
    n  <- nrow(df)
    p  <- ncol(df) - 1

    colnames(df) <- c("y", paste0("x", 1:p))
    df$y <- as.numeric(df$y == "spam")

    to3 <- function(x) x^3
    transforms <- c("sigmoid","sin_deg","exp_dbl","p0","troot","to3")
    probs <- gen.probs.gmjmcmc(transforms)
    probs$gen <- c(1, 1, 1, 1)

    params <- gen.params.gmjmcmc(p)
    params$feat$check.col <- FALSE

    set.seed(6001)
    result <- fbms(
      formula    = y ~ 1 + .,
      data       = df,
      method     = "gmjmcmc",
      family     = "binomial",
      beta_prior = list(type = "Jeffreys-BIC"),
      transforms = transforms,
      probs      = probs,
      params     = params
    )

    printn(summary(result))

    ##
    ## Best   population: 2  log marginal posterior: -1048.305
    ##
    ##             feats.strings  marg.probs
    ## 1                      x5 1.000000000
    ## 2                      x8 1.000000000
    ## 3                     x55 1.000000000
    ## 4            sigmoid(x21) 1.000000000
    ## 5                     x16 1.000000000
    ## 6                     x20 1.000000000
    ## 7                     x21 1.000000000
    ## 8                     x23 1.000000000
    ## 9                     x25 1.000000000
    ## 10                    x27 1.000000000
    ## 11            sigmoid(x7) 1.000000000
    ## 12                p0(x17) 1.000000000
    ## 13               to3(x25) 1.000000000
    ## 14                    x41 1.000000000
    ## 15                    x42 1.000000000
    ## 16                    x44 1.000000000
    ## 17                    x12 1.000000000
    ## 18                     x6 1.000000000
    ## 19           sigmoid(x52) 1.000000000
    ## 20                    x53 1.000000000
    ## 21                    x52 1.000000000
    ## 22                    x45 1.000000000
    ## 23              (x50*x25) 1.000000000
    ## 24                p0(x19) 1.000000000
    ## 25                p0(x22) 0.999999632
    ## 26 sigmoid(1+1*x16+1*x55) 0.999999632
    ## 27               (x24*x3) 0.998106918
    ## 28                     x2 0.998091170
    ## 29                    x28 0.998090532
    ## 30             troot(x27) 0.998090188
    ## 31                    x39 0.998090175
    ## 32                    x32 0.998090164
    ## 33                    x46 0.994872577
    ## 34                    x19 0.994843934
    ## 35              troot(x8) 0.807086602
    ## 36                    x24 0.807017530
    ## 37                    x49 0.806845676
    ## 38                    x29 0.036497149
    ## 39                    x40 0.017892855
    ## 40                    x10 0.014132910
    ## 41           sigmoid(x46) 0.003217954

Prediction accuracy

    # Logistic link
    invlogit <- function(x) 1/(1 + exp(-x))

    # Model averaging
    pred <- predict(result, x = df[,-1], link = invlogit)
    printn(mean(round(pred$aggr$mean) == df$y))

    ## [1] 0.9397957

    # Best model
    bm <- get.best.model(result)
    preds <- predict(bm, df[,-1], link = invlogit)
    printn(mean(round(preds) == df$y))

    ## [1] 0.9402304

    # Median Probability Model
    mpm <- get.mpm.model(result, family = "binomial", y = df$y, x = df[,-1])
    preds_mpm <- predict(mpm, df[,-1], link = invlogit)
    printn(mean(round(preds_mpm) == df$y))

    ## [1] 0.9402304

# References

- Hubin, A., Storvik, G. (2018). Mode jumping MCMC for Bayesian variable selection in GLMs.
- Hubin, A., Frommlet, F., & Storvik, G. (2020). Flexible Bayesian Model Averaging for Generalized Nonlinear Models.
- Hubin, A., et al. (2020). Bayesian Logic Regression via GMJMCMC.
- Hubin, A., et al. (2021). Reversible GMJMCMC.
- Lachmann, J., et al. (2022). Subsampling for tall data in GMJMCMC.