Introduction to Logistic Box-Cox Regression with lboxcox

Introduction

The lboxcox package fits logistic Box-Cox (LBC) regression models for a binary outcome and a strictly positive continuous predictor. The model is useful when the predictor-outcome relationship may be nonlinear but a compact, interpretable parametric form is preferred to a fully nonparametric fit.

Ordinary logistic regression assumes that a continuous predictor has a linear effect on the log-odds scale. LBC regression relaxes this assumption by applying a Box-Cox transformation to the primary predictor and estimating its shape parameter from the data. The original model and its median-effect interpretation were developed by Xing et al. (2021).

Estimating the shape parameter requires nonlinear optimization, which can be sensitive to starting values. The current package therefore also implements the multi-start and bootstrap-aggregation procedures developed by Xu, Wang, and Xing. These additions provide four related fitting strategies: single-fit maximum likelihood (LBC-ML), multi-start fitting (LBC-MS), bootstrap aggregation of LBC-ML fits (LBC-EL), and bootstrap aggregation with multi-start fitting within each resample (LBC-CM).

Model

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

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

For a binary outcome \(Y_i\), a positive primary predictor \(X_i\), and adjustment covariates \(\mathbf Z_i\), the LBC model is

\[ \operatorname{logit}\{\Pr(Y_i=1)\} = \beta_0 + \beta_1 X_i^{(\lambda)} + \boldsymbol{\gamma}^{\mathsf T}\mathbf Z_i. \]

In a package formula such as y ~ x + z1 + z2, the first term on the right-hand side is treated as the primary predictor and receives the Box-Cox transformation. The remaining terms are adjustment covariates.

Parameter interpretation

The shape parameter \(\lambda\) controls the form of the predictor-outcome relationship. Values near 0 correspond to a logarithmic transformation, \(\lambda=1\) gives a linear term, and larger values permit increasingly convex relationships. The coefficient \(\beta_1\) gives the direction and strength of association on the transformed scale and should be interpreted together with the fitted value of \(\lambda\).

Median effect

For an approximately log-normal predictor with log-scale location \(\mu\), the median effect on the original predictor scale is

\[ \Delta^* = \beta_1 \exp\{(\lambda-1)\mu\}. \]

median_effect() evaluates this summary at the weighted mean of the log predictor and returns a Wald-type 95% confidence interval. Direct maximum-likelihood fits use the joint likelihood Hessian. Profile-grid, cross-validation, and ensemble refits return an interval conditional on the selected lambda.

Data requirements

The response must be binary, and the primary continuous predictor must be strictly positive. weight_column_name may identify a column containing non-negative observation weights or be a numeric vector with one value per row. Use NULL or 1 for an unweighted analysis:

fit <- lbc_maxlik(
  y ~ x + z1 + z2,
  weight_column_name = NULL,
  data = mydata
)

The package incorporates observation weights but does not accept survey strata or primary sampling-unit identifiers. The resulting fits are therefore sampling-weighted model estimates rather than complete design-based survey estimates.

Fitting the models

library(lboxcox)
#> Loading required package: survey
#> Loading required package: grid
#> Loading required package: Matrix
#> Loading required package: survival
#> 
#> Attaching package: 'survey'
#> The following object is masked from 'package:graphics':
#> 
#>     dotchart
data(depress)

formula_lbc <- depression ~ mercury + age + factor(gender)

LBC-ML: single-fit maximum likelihood

lbc_maxlik() fits one LBC model by maximum likelihood. By default, survey-weighted logistic fits over a lambda grid are used to construct a starting vector before direct likelihood maximization.

fit_ml <- lbc_maxlik(
  formula_lbc,
  weight_column_name = "weight",
  data = depress,
  seed = 1
)

fit_ml$estimate
#>         Beta_0         Beta_1         Lambda            age factor(gender) 
#>   -2.208855660   -0.365758960    0.203435933   -0.003480318   -0.491747914

LBC-MS: multi-start fitting

lbc_train_ms() evaluates the likelihood from multiple starting lambda values, selects the lambda associated with the highest achieved log-likelihood, and returns a full-data weighted logistic refit conditional on that lambda.

fit_ms <- lbc_train_ms(
  formula_lbc,
  weight_column_name = "weight",
  data = depress,
  svy_lambda_vector = seq(0, 2, length.out = 10)
)

fit_ms$estimate
#>          Beta_0          Beta_1          Lambda             age factor(gender)1 
#>    -2.206312885    -0.365895449     0.208246689    -0.003513522    -0.492101233

LBC-EL and LBC-CM: bootstrap aggregation

lbc_train_bagging() fits LBC-ML models to 100 bootstrap samples. This is the LBC-EL procedure. lbc_train_all() applies multi-start fitting within each bootstrap sample and implements LBC-CM; it is consequently more computationally intensive.

Both functions return a list containing a full-data refit at the median of the bootstrap-specific lambda estimates, the individual bootstrap fits, and the number of bootstrap calls that did not return a usable fit. The median lambda is a descriptive summary. Ensemble predictions are obtained by averaging predicted probabilities across the available bootstrap fits.

set.seed(1)
fit_el <- lbc_train_bagging(
  formula_lbc,
  weight_column_name = "weight",
  data = depress,
  cores = 2
)

fit_cm <- lbc_train_all(
  formula_lbc,
  weight_column_name = "weight",
  data = depress,
  cores = 2
)

The bootstrap examples are not evaluated when this vignette is built because each procedure fits 100 resampled datasets.

Prediction and model evaluation

Use lboxcox_maxLik.predict() with an LBC-ML or LBC-MS fit. Use lboxcox_maxLik_el.predict() with an LBC-EL or LBC-CM result.

p_ml <- lboxcox_maxLik.predict(fit_ml, depress, formula_lbc)
p_ms <- lboxcox_maxLik.predict(fit_ms, depress, formula_lbc)

head(p_ml)
#>            [,1]
#> [1,] 0.02383239
#> [2,] 0.07275958
#> [3,] 0.05847010
#> [4,] 0.03872629
#> [5,] 0.01688090
#> [6,] 0.04620477

p_el <- lboxcox_maxLik_el.predict(fit_el, depress, formula_lbc)

devr() computes the sum of absolute deviance residuals (SADR). Lower values indicate better predictive performance when models are evaluated on the same observations.

devr(depress$depression, p_ml)
#> [1] 4718.065
devr(depress$depression, p_ms)
#> [1] 4718.08

As an alternative to likelihood-based estimation, lboxcox_cv.fit() selects lambda from a user-supplied grid by minimizing cross-validated SADR.

fit_cv <- lboxcox_cv.fit(
  mydata = depress,
  ixx = depress$mercury,
  iyy = depress$depression,
  formula = formula_lbc,
  weight_column_name = "weight",
  lambda_vector = seq(0, 2, length.out = 10),
  k = 5
)

p_cv <- lboxcox_cv.predict(fit_cv, depress, formula_lbc)

For a direct LBC-ML fit, the median-effect summary is obtained with:

median_effect(
  formula_lbc,
  weight_column_name = "weight",
  data = depress,
  trained_model = fit_ml
)
#> median effect  lower 95% ci  upper 95% ci 
#>    -0.3666334    -0.4616603    -0.2716064

Built-in NHANES data

The bundled depress data frame contains the 8,893 adults aged 20 years or older used in the NHANES application. The analytic sample combines the 2005–2006 and 2007–2008 survey cycles and contains:

depression: indicator equal to 1 for a Patient Health Questionnaire-9 (PHQ-9) score of at least 10;
mercury: total blood mercury concentration in micrograms per litre;
age: age in years;
gender: 1 for male and 0 for female; and
weight: four-year Day 1 dietary sampling weight, formed as WTDRD1 / 2 for the two combined cycles.

summary(depress)
#>    depression         mercury            age            gender      
#>  Min.   :0.00000   Min.   : 0.140   Min.   :20.00   Min.   :0.0000  
#>  1st Qu.:0.00000   1st Qu.: 0.490   1st Qu.:34.00   1st Qu.:0.0000  
#>  Median :0.00000   Median : 0.890   Median :49.00   Median :0.0000  
#>  Mean   :0.08074   Mean   : 1.478   Mean   :49.47   Mean   :0.4871  
#>  3rd Qu.:0.00000   3rd Qu.: 1.670   3rd Qu.:64.00   3rd Qu.:1.0000  
#>  Max.   :1.00000   Max.   :38.700   Max.   :85.00   Max.   :1.0000  
#>      weight        
#>  Min.   :   293.9  
#>  1st Qu.:  6893.1  
#>  Median : 15051.3  
#>  Mean   : 21598.5  
#>  3rd Qu.: 27847.6  
#>  Max.   :169230.1

See ?depress for the variable definitions and source details.

Main functions

Function	Purpose
`lbc_maxlik()`	Fit one LBC model by maximum likelihood
`lbc_train_ms()`	Fit an LBC model from multiple starting values
`lbc_train_bagging()`	Fit the LBC-EL bootstrap procedure
`lbc_train_all()`	Fit the LBC-CM combined procedure
`lboxcox_maxLik.predict()`	Predict from an LBC-ML or LBC-MS fit
`lboxcox_maxLik_el.predict()`	Average predictions across bootstrap fits
`lboxcox_cv.fit()`	Select lambda by cross-validated SADR
`lboxcox_cv.predict()`	Predict from a cross-validated fit
`devr()`	Calculate SADR
`median_effect()`	Calculate the median-effect summary and confidence interval

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. https://doi.org/10.1002/cjs.11587

Xu, S., Wang, J., & Xing, L. Ensemble Logistic Box-Cox Model for Improved Prediction and Estimation. Manuscript in preparation.

Lumley, T. (2011). Complex Surveys: A Guide to Analysis Using R. John Wiley & Sons.