--- title: "Bivariate Logistic Box-Cox Regression with mvboxcox" output: html_document: toc: true toc_depth: 3 toc_float: collapsed: false smooth_scroll: true vignette: > %\VignetteIndexEntry{Bivariate Logistic Box-Cox Regression with mvboxcox} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE ) ``` ## 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: ```r 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. ```{r simulation} 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. ```{r fit, results = "hide"} fit <- mvbc.train( Ybin ~ mercury + lead, ~ age, data = sim_train$data, griddomain = seq(0, 2, length.out = 5), K = 3, depth = 1, seed = 1 ) ``` ```{r prediction} p_hat <- mvbc.predict(fit, sim_test$data) head(p_hat) mvbc.trainer.ssr(sim_test$data$Ybin, p_hat) ``` 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. ```{r inspect-fit} fit$lambda.fits fit$beta.fits fit$grid[which.min(fit$grid$ssdr), ] ``` For the simulated predictors, the empirical-median plug-in estimates at $q=1$ can be obtained directly from the fitted object. ```{r median-effect} median_effect_q1 <- mvbc.median.effect( fit, sim_train$data, q = 1 ) median_effect_q1 ``` ## 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. ```{r depress-data} data(depress, package = "mvboxcox") dim(depress) head(depress) table(depress$depression) ``` 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. ```{r nhanes-fit, results = "hide"} 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: ```{r nhanes-results} fit_nhanes$lambda.fits fit_nhanes$beta.fits mvbc.median.effect( fit_nhanes, depress, q = 1, weights = depress$weight ) ``` 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.