Variable selection for linear models and generalized linear models with BIC-based posterior probability

In VariableSelection: Select Variables for Linear Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(VariableSelection)

Introduction

Variable selection methods can screen and identify associated variables from regression variables. A simpler model with fewer variables is easier to understand in the real world.

The VariableSelection package uses Bayesian Information Criterion (BIC), Model Posterior Probability (MPP), and Posterior Inclusion Probability (PIP) to perform model selection for linear models (LM) and generalized linear models (GLM). The package provides the best model, BIC, and MPP for candidate models, and PIP for each predictor. This vignette contains an example to illustrate how to use VariableSelection.

Model

The linear models used in the VariableSelection package are of the form $$ \pmb{Y}=X\pmb{\beta}+ \pmb{\epsilon},$$ where

• $\pmb{Y}$ is the vector of observations.
• $X$ is the matrix of covariates.
• $\pmb{\beta}$ is the vector of regression coefficients.
• $\pmb{\epsilon}$ is the vector of errors.

The generalized linear models used in the VariableSelection package assume that the observations come from a distribution that belongs to the exponential family of distributions. In addition, these GLMs assume that the mean of the observations is related to the covariates through a link function $g$ such that $$E(\pmb{Y}|X) = g^{-1}(X\pmb{\beta}),$$ where

• $\pmb{Y}$ is the vector of observations.
• $X$ is the matrix of covariates.
• $\pmb{\beta}$ is the vector of regression coefficients.
• $E(\pmb{Y}|X)$ is the expected value of $\pmb{Y}$ conditional on $X$.
• $g()$ is the link function.

BIC

The Bayesian Information Criterion (BIC) is defined as:

$$ \text{BIC} = -2 \log(L) + k \log(n), $$

where

• $L$ is the likelihood function of the model computed at the MLE.
• $k$ is the number of parameters in the model.
• $n$ is the number of observations.

Lower value of BIC indicates a better model.

Example

Data

The modelselect.lm() function can take a data frame which contains both the response and predictors. For example, here are the first 5 rows of a data frame, where X1 through X6 are six candidate predictors. Y is the response variable which is simulated from a linear model that contains the predictors X1, X2, and X3.

data("dat")
head(dat,5)

The data frame dat above is attached in the VariableSelection package.

Formula

In this example, we use modelselect.lm to select the predictors in a linear model. The modelselect.lm function takes a formula and dataset as input. In the example below, we consider the model space that contains all possible linear models generated with the predictors X1, X2, X3, and X4.

example1 <- modelselect.lm(formula = Y~X1+X2+X3+X4, data = dat)

The output of modelselect.lm returns a table of BICs and MPPs of competing models and a table of PIPs of candidate predictors.

head(example1$models)
example1$variables

Here are some additional examples on how to write a formula. In the next example, the formula ~. includes all predictors in the data frame.

example2 <- modelselect.lm(formula = Y~., data = dat)
example2$variables

Interactions

The next example includes an interaction term between the predictors X1 and X2.

example3 <- modelselect.lm(formula = Y~X1*X2+X3+X4+X5+X6, data = dat)
example3$variables

Note that because modelselect.lm uses lm notation for formulas, we could also have used notation X1+X2+X3+X4+X5+X6+X1:X2.

Stochastic search

modelselect.lm function uses GA_var as a threshold to do exhaustive model search or stochastic model search. The default value of GA_var is 16. If the number of variables is smaller than GA_var, then modelselect.lm does exhaustive model search, otherwise it uses genetic algorithm (GA) to perform stochastic model search. maxiterations is the maximum number of iterations to run before the GA search is stopped runs_til_stop is the number of consecutive generations without any improvement in the best fitness value before the GA is stopped. monitor is a logical defaulting to TRUE showing the evolution of the search. If monitor = FALSE, any output is suppressed.

example5 <- modelselect.lm(formula = Y~X1+X2+X3+X4, data = dat, GA_var = 16, maxiterations = 2000, runs_til_stop = 1000, monitor = TRUE, popSize = 100)

Model fit

Use lm.best to obtain the model fit for the best model. lm.best takes the result from modelselect.lm as an input. In this example, we obtain the model fit by using method = "models". The fitted model is the model with the highest MPP. The return of lm.best is the same as that from lm. We may use $ to obtain the output statistics, for example $coefficients for regressor coefficients' estimate.

lm_model <- lm.best(object = example1, method = "models")
lm_model$coefficients

In the next example, we obtain the model fit by using method = "variables". The fitted model has the predictors with PIP larger than threshold.

lm_var <- lm.best(object = example2, method = "variables", threshold = 0.9)

Note that because the output of lm.best function is of class lm, we can apply the function summary to the output of lm.best. Note that the summary table should be interpreted as conditional on the best model being the true model.

summary(lm_model)

GLM

Here is an example on how to perform model selection for generalized linear models. The modelselect.glm() function takes a data frame which contains response and predictors. In this toy example, the data frame glmdat contains six candidate predictors X1 to X6. The response variable Y was actually simulated from a Bernoulli GLM with predictors X1, X2, and X3.

data("glmdat")
head(glmdat,5)

Data glmdat above are attached in the VariableSelection package. In the next example, we use the function modelselect.glm to find the best model and to compute the PIPs of the several candidate predictors.

example.glm <- modelselect.glm(formula = Y~., family = "binomial", data = glmdat)
example.glm$variables

Then, we use glm.best to obtain the model fit for the best model. glm.best takes the result from modelselect.glm as an argument.

glm_model <- glm.best(object = example.glm, family = "binomial", method = "models", threshold = 0.95)

The function summary can be applied to the result of glm.best.

summary(glm_model)

Functions

In the VariableSelection package, there are four functions: modelselect.lm() and lm.best() for LM, and modelselect.glm() and glm.best() for GLM.

• Function modelselect.lm() uses BIC to perform model selection for LMs. The function modelselect.lm() takes a formula for the full model and a data frame containing the response variable and predictor variables as input. It returns a list of two tables: 1. a table for candidate models with BIC and posterior probabilities; 2. a table for candidate variables with posterior inclusion probabilities. In addition, it returns the original data with variables in the formula.

• Function lm.best() takes result from modelselect.lm() as object. There are two methods to select the best model. method="models" uses models' BIC or posterior probabilities to select the best model. method="variables" selects the variables with PIP larger than the threshold.

• Function modelselect.glm() uses BIC to perfrom model selection for GLMs. The function modelselect.glm() takes formula, data containing response variable and predictors, and family function for error distribution as input. It returns a list of two tables: 1. a table for candidate models with BIC and posterior probabilities; 2. a table for candidate variables with posterior inclusion probabilities. In addition,it returns the original data with variables in the formula.

• Function glm.best() takes result from modelselect.glm() as object. There are two methods to select the best model. method="models" uses models' BIC or posterior probabilities to select the best model. method="variables" selects the variables with PIP larger than the threshold.

Reference

Xu, Shuangshuang, Tegge, Allison, and Ferreira, M. A. R. (202X). paper, Journal, .

VariableSelection documentation built on Feb. 17, 2026, 5:07 p.m.