FBMS (version 1.1)

get.mpm.model: Retrieve the Median Probability Model (MPM)

Description

This function extracts the Median Probability Model (MPM) from a fitted model object. The MPM includes features with marginal posterior inclusion probabilities greater than 0.5. It constructs the corresponding model matrix and computes the model fit using the specified likelihood.

Usage

get.mpm.model(
  result,
  y,
  x,
  labels = F,
  family = "gaussian",
  loglik.pi = gaussian.loglik,
  params = NULL
)

Value

A bgnlm_model object containing:

prob

The log marginal likelihood of the MPM.

model

A logical vector indicating included features.

crit

Criterion label set to "MPM".

coefs

A named numeric vector of model coefficients, including the intercept.

Arguments

result

A fitted model object (e.g., from mjmcmc, gmjmcmc, or related classes) containing the summary statistics and marginal probabilities.

y

A numeric vector of response values. For family = "binomial", it should contain binary (0/1) responses.

x

A data.frame of predictor variables. Columns must correspond to features considered during model fitting.

labels

If specified, custom labels of covariates can be used. Default is FALSE.

family

Character string specifying the model family. Supported options are:

• "gaussian" (default) - for continuous outcomes.

• "binomial" - for binary outcomes.

• "custom" - for user-defined likelihood functions.

If an unsupported family is provided, a warning is issued and the Gaussian likelihood is used by default.

loglik.pi

A function that computes the log-likelihood. Defaults to gaussian.loglik unless family = "binomial", in which case logistic.loglik is used. for custom family the user must specify the same likelihood that was used in the inference.

params

Parameters of loglik.pi, if not specified NULL will be used by default

Examples

if (FALSE) {
# Simulate data
set.seed(42)
x <- data.frame(
  PlanetaryMassJpt = rnorm(100),
  RadiusJpt = rnorm(100),
  PeriodDays = rnorm(100)
)
y <- 1 + 0.5 * x$PlanetaryMassJpt - 0.3 * x$RadiusJpt + rnorm(100)

# Assume 'result' is a fitted object from gmjmcmc or mjmcmc
result <- mjmcmc(cbind(y,x))  

# Get the MPM
mpm_model <- get.mpm.model(result, y, x, family = "gaussian")

# Access coefficients
mpm_model$coefs
}