---
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.
