Introduction

mvboxcox fits bivariate logistic Box-Cox (BLBC) regression models for a binary outcome and two positive continuous predictors. It is intended for settings in which the predictor-outcome relationships may be nonlinear but a compact, interpretable parametric model is preferred to a fully nonparametric fit.

The package estimates a separate Box-Cox shape parameter for each positive predictor. Candidate shape parameters are evaluated by K-fold cross-validation, with optional adaptive grid refinement and thin-plate spline (TPS) refinement. The package also provides simulation, prediction, testing, and sampling-weighted fitting tools.

Installation

Install the package with:

install.packages("mvboxcox")

Model

For a positive value \(x\), the Box-Cox transformation is

\[ x^{(\lambda)} = \begin{cases} (x^\lambda - 1)/\lambda, & \lambda \ne 0, \\ \log(x), & \lambda = 0. \end{cases} \]

For two positive predictors \(X_{i1}\) and \(X_{i2}\) and a covariate vector \(\mathbf Z_i\), the BLBC model is

\[ \operatorname{logit}\{\Pr(Y_i=1)\} = \beta_0 + \beta_1 X_{i1}^{(\lambda_1)} + \beta_2 X_{i2}^{(\lambda_2)} + \boldsymbol{\gamma}^{\mathsf T}\mathbf Z_i. \]

In the package interface:

  • formulaA includes the binary outcome and the positive predictors that receive Box-Cox transformations.
  • formulaB contains covariates that remain on their original model-matrix scale and does not include the outcome.

For example, Ybin ~ mercury + lead specifies the transformed predictors, whereas ~ age adds age without a Box-Cox transformation.

Interpreting the parameters

Each \(\lambda_j\) describes the shape of a predictor-response relationship, whereas \(\beta_j\) describes its direction and strength on the corresponding transformed scale. The magnitude of \(\beta_j\) depends on \(\lambda_j\), so raw coefficients associated with different transformation shapes should not be compared directly.

Let \(m_j\) denote the marginal median of a positive predictor \(X_j\). A common-scale summary is the median effect relative to a reference transformation \(q\):

\[ \delta_j(q) = \beta_j m_j^{\lambda_j-q}. \]

At \(q=1\), \(\delta_j(1)\) is the effective linear slope at the median predictor value. This definition requires no parametric distribution for \(X_j\). If \(X_j\) is log-normal with location parameter \(\mu_j\), then \(m_j=\exp(\mu_j)\) and the expression reduces to \(\delta_j(q)=\beta_j\exp\{(\lambda_j-q)\mu_j\}\). Fixing \(q\) provides a common reference scale across transformation shapes, although numerical comparisons between different predictors remain dependent on their measurement units.

Weak-signal interpretation

The shape parameter is identified through the predictor’s association with the outcome. When the corresponding coefficient is close to zero, the data contain little information about the shape, and the estimated \(\lambda_j\) may be unstable. A flat cross-validation surface, a boundary estimate, or variation in the selected shape can therefore indicate a weak predictor-outcome signal rather than an optimization failure.

Simulate, fit, and predict

The following example generates independent training and test datasets from the same bivariate logistic Box-Cox model.

library(mvboxcox)

sim_train <- mvbc.simulator(
  vLambda = c(0.5, 1.5),
  vBeta = c(-2.2, -0.4, -0.2, -0.005),
  vMean = c(-0.08, -0.01, 50),
  vSd = c(0.93, 0.8, 18.12),
  vNames = c("mercury", "lead", "age"),
  n = 1000,
  seed = 1
)

sim_test <- mvbc.simulator(simModel = sim_train)

mvbc.train() runs the complete fitting pipeline: it searches the candidate lambda grid by K-fold cross-validation and then refines the selected minimum using a TPS surface and bounded optimization. A small grid is used here so the vignette builds quickly.

fit <- mvbc.train(
  Ybin ~ mercury + lead,
  ~ age,
  data = sim_train$data,
  griddomain = seq(0, 2, length.out = 5),
  K = 3,
  depth = 1,
  seed = 1
)
p_hat <- mvbc.predict(fit, sim_test$data)
head(p_hat)
#>         [,1]
#> 1 0.07567006
#> 2 0.06729733
#> 3 0.03173852
#> 4 0.11982941
#> 5 0.12117104
#> 6 0.05727863
mvbc.trainer.ssr(sim_test$data$Ybin, p_hat)
#>      sspr      ssdr 
#> 1149.6138  570.6824

The fitted object stores one selected shape estimate for each transformed predictor in lambda.fits. Its beta.fits component contains the coefficient estimate obtained by refitting the model on the full training data at the selected lambda tuple. The value foldid = 0 identifies this full-data refit. mvbc.predict() calculates probabilities using this final coefficient estimate.

fit$lambda.fits
#>     mercury lead
#> 1 0.8235748    0
fit$beta.fits
#>   foldid intercept    mercury      lead         age
#> 1      0 -2.755889 -0.4583821 -0.385433 0.006281219

fit$grid[which.min(fit$grid$ssdr), ]
#>      mercury lead     sspr    ssdr
#> 26 0.8235748    0 996.7963 547.778

For the simulated predictors, the empirical-median plug-in estimates at \(q=1\) can be obtained directly from the fitted object.

median_effect_q1 <- mvbc.median.effect(
  fit,
  sim_train$data,
  q = 1
)
median_effect_q1
#>   predictor       beta    lambda q    median median_effect
#> 1   mercury -0.4583821 0.8235748 1 0.8932834    -0.4675999
#> 2      lead -0.3854330 0.0000000 1 0.9631115    -0.4001956

Built-in NHANES data

The package includes depress, the analytic dataset used in the paper’s NHANES application. It contains 8,893 adults from the 2005-2006 and 2007-2008 cycles and six variables:

  • depression: indicator equal to 1 for a PHQ-9 score of at least 10;
  • mercury: total blood mercury concentration in micrograms per liter;
  • blood_lead: blood lead concentration in micrograms per deciliter;
  • age: age in years;
  • gender: 1 for male and 0 for female; and
  • weight: combined-cycle Day 1 dietary sampling weight.
data(depress, package = "mvboxcox")
dim(depress)
#> [1] 8893    6
head(depress)
#>   depression mercury blood_lead age gender    weight
#> 1          0    7.13       1.01  37      1 18979.728
#> 2          0    1.40       2.30  60      0 19877.753
#> 3          0    2.46       2.17  60      0 16671.291
#> 4          0    6.24       2.10  55      0  3749.384
#> 5          0   11.20       1.52  64      1 66567.646
#> 6          0    1.48       3.69  51      1  5494.848
table(depress$depression)
#> 
#>    0    1 
#> 8175  718

See ?depress for the variable definitions and data-source details.

Sampling-weighted fitting

mvbc.train(), mvbc.trainer(), and mvbc.optimizer() accept an optional observation-weight vector. When weights is supplied with survey = TRUE, candidate models are fitted with a weighted glm() instead of an unweighted glm(), and the held-out residual criteria are weighted accordingly.

The sampling-weighted BLBC model used for the NHANES application is fitted below using the same \(10\times10\) initial grid, five cross-validation folds, and refinement depth \(D=2\) as in the manuscript.

fit_nhanes <- mvbc.train(
  depression ~ mercury + blood_lead,
  ~ age + factor(gender),
  data = depress,
  weights = depress$weight,
  survey = TRUE,
  griddomain = seq(0, 2, length.out = 10),
  K = 5,
  depth = 2,
  seed = 1
)

The selected transformation parameters, full-data coefficient refit, and sampling-weighted median effects are:

fit_nhanes$lambda.fits
#>     mercury blood_lead
#> 1 0.2543137          0
fit_nhanes$beta.fits
#>   foldid intercept    mercury blood_lead         age factor(gender)1
#> 1      0 -1.937421 -0.3758578  0.4173978 -0.01068292      -0.6647706

mvbc.median.effect(
  fit_nhanes,
  depress,
  q = 1,
  weights = depress$weight
)
#>    predictor       beta    lambda q median median_effect
#> 1    mercury -0.3758578 0.2543137 1   0.97    -0.3844924
#> 2 blood_lead  0.4173978 0.0000000 1   1.36     0.3069101

The current survey-weighted implementation uses observation weights but does not accept survey strata or primary sampling-unit identifiers. Results should therefore be interpreted as sampling-weighted model estimates rather than as a complete design-based NHANES survey analysis. As with other observational analyses, the fitted associations are not by themselves causal effects.

Main functions

  • mvbc.simulator() generates data from a logistic Box-Cox model.
  • mvbc.trainer() evaluates candidate lambda tuples by cross-validation.
  • mvbc.optimizer() refines a trained lambda surface using TPS and L-BFGS-B.
  • mvbc.train() runs the trainer and optimizer as one pipeline.
  • mvbc.predict() predicts probabilities from a fitted model.
  • mvbc.median.effect() computes empirical or sampling-weighted median effects.
  • mvbc.tester() performs repeated simulation-based train/test evaluations.
  • mvbc.trainer.ssr() calculates Pearson and deviance residual criteria.

Use help(package = "mvboxcox") for the complete function reference.

References

Box, G. E. P., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society: Series B (Methodological), 26(2), 211-243.

Xing, L., Zhang, X., Burstyn, I., & Gustafson, P. (2021). On logistic Box-Cox regression for flexibly estimating the shape and strength of exposure-disease relationships. Canadian Journal of Statistics, 49(3), 808-825.

Xu, S., & Zhang, X. Bivariate logistic Box-Cox regression for interpretable nonlinear exposure-response modeling. Manuscript.