mvboxcox implements an interpretable nonlinear extension of logistic
regression for binary outcomes and positive continuous predictors, described
in Bivariate logistic Box-Cox regression for interpretable nonlinear
exposure-response modeling by Shiyu Xu and Xuekui Zhang. It is particularly
useful when predictor-outcome relationships may be nonlinear, including
applications involving right-skewed exposures. The package extends the
univariate logistic Box-Cox model of Xing et al. (2021) [2] to two
positive predictors, each with its own data-driven Box-Cox transformation
parameter, estimated by K-fold cross-validation with adaptive grid refinement
and thin-plate spline (TPS) refinement.
Install the package with:
install.packages("mvboxcox")
For a positive predictor (x > 0), the Box-Cox transformation is (x^{(\lambda)} = (x^\lambda - 1)/\lambda) for (\lambda > 0), and (\log(x)) for (\lambda = 0). The bivariate model is
[ \operatorname{logit}{\Pr(Y_i = 1)} = \beta_0 + \beta_1 X_{i1}^{(\lambda_1)} + \beta_2 X_{i2}^{(\lambda_2)}
where (X_{i1}, X_{i2}) are positive exposures (formulaA) and
(\mathbf Z_i) holds covariates left untransformed (formulaB).
(\lambda_j) describes the shape of the exposure-response relationship;
(\beta_j) describes the association on the transformed scale. Because the
magnitude of (\beta_j) depends on that scale, coefficients from different
transformation shapes should not be compared directly.
Let (m_j) denote the marginal median of a positive exposure (X_j). The median effect relative to a reference transformation (q) is
[ \delta_j(q) = \beta_j m_j^{\lambda_j-q}. ]
At (q=1), this is the effective linear slope at the median exposure and provides a common scale for comparing effects across transformation shapes. No parametric distribution is required. 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}). Numerical comparisons between different predictors remain dependent on their measurement units.
The shape parameter for an exposure is identified through its association with the outcome. When the corresponding coefficient is close to zero, the shape parameter can be weakly identified. A flat cross-validation surface, a boundary estimate, or variability in the selected lambda may therefore reflect limited exposure-outcome information rather than numerical failure.
library(mvboxcox)
sim <- 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 = 2000, seed = 1
)
sim_test <- mvbc.simulator(simModel = sim)
fit <- suppressWarnings(
mvbc.train(Ybin ~ mercury + lead, ~ age, sim$data, depth = 2)
)
p_hat <- mvbc.predict(fit, sim_test$data)
mvbc.trainer.ssr(sim_test$data$Ybin, p_hat)
mvbc.median.effect(fit, sim$data, q = 1)
mvbc.train(), mvbc.trainer(), and mvbc.optimizer() also accept an
observation weights vector; with survey = TRUE they fit each candidate
model with a weighted glm() instead of an unweighted glm(), and weight the
cross-validation residual criteria accordingly. See
vignette("introduction", package = "mvboxcox") for a full walkthrough,
including a weighted example.
The survey-weighted implementation incorporates observation weights but does not accept survey strata or primary sampling-unit identifiers. Results should therefore be interpreted as sampling-weighted rather than as complete design-based survey estimates.
The bundled depress data frame (8,893 rows, 6 columns) is the analytic sample
from the paper's NHANES 2005-2006 and 2007-2008 application. It contains
depression status, blood mercury and blood_lead, age, gender (1 =
male, 0 = female), and the combined-cycle Day 1 dietary sampling weight.
data(depress, package = "mvboxcox")
summary(depress)
The paper's sampling-weighted BLBC model can be fitted with:
fit_nhanes <- mvbc.train(
depression ~ mercury + blood_lead,
~ age + factor(gender),
data = depress,
weights = depress$weight,
survey = TRUE
)
mvbc.median.effect(
fit_nhanes,
depress,
q = 1,
weights = depress$weight
)
| Function | Purpose |
|---|---|
mvbc.simulator() |
Generate data from a logistic Box-Cox model |
ld() |
Apply the Box-Cox transformation |
mvbc.trainer() |
Evaluate lambda tuples by K-fold cross-validation |
mvbc.optimizer() |
Refine a trained lambda surface using TPS and L-BFGS-B |
mvbc.train() |
Run the trainer and optimizer as an end-to-end pipeline |
mvbc.predict() |
Predict probabilities from an mvbc.train.model
|
mvbc.median.effect() |
Compute empirical or sampling-weighted median effects |
mvbc.tester() |
Evaluate repeated simulated train/test experiments |
mvbc.trainer.ssr() |
Calculate Pearson and deviance residual criteria |
Use help(package = "mvboxcox") for the complete reference.
[1] 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.
[2] 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.
[3] Xu, S., & Zhang, X. Bivariate logistic Box-Cox regression for interpretable nonlinear exposure-response modeling. Manuscript.