Package {mvboxcox}

NHANES Depression and Blood Metals Data

Description

A sampling-weighted analytic dataset derived from the 2005-2006 and 2007-2008 cycles of the US National Health and Nutrition Examination Survey (NHANES). The data are used in the bivariate logistic Box-Cox application described by Xu and Zhang.

Usage

depress

Format

A data frame with 8,893 rows and 6 variables:

depression

Integer indicator: 1 for a Patient Health Questionnaire-9 (PHQ-9) score of at least 10, and 0 otherwise.

mercury

Total blood mercury concentration in micrograms per liter (mcg/L).

blood_lead

Blood lead concentration in micrograms per deciliter (mcg/dL).

age

Age in years.

gender

Integer indicator: 1 for male and 0 for female.

weight

Four-year Day 1 dietary sampling weight, constructed for the combined NHANES cycles.

Details

The analytic sample contains adults aged 20 years or older with a valid Day 1 dietary sampling weight, complete PHQ-9 information, and nonmissing positive blood mercury and blood lead measurements. The package's survey-weighted implementation uses weight but does not incorporate NHANES strata or primary sampling-unit identifiers, so analyses should be interpreted as sampling-weighted rather than as complete design-based survey analyses.

Source

National Center for Health Statistics, National Health and Nutrition Examination Survey, 2005-2006 and 2007-2008 cycles.

References

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

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.

Examples

data(depress)
str(depress)
summary(depress)

Boxcox transformation

Description

if (li == 0) log(xi) else (xi**li - 1)/li

Usage

ld(xi, li)

Arguments

Value

numeric vector of Box-Cox transformed values

Build the full Box-Cox model formula

Description

Apply a set of lambda exponents to the positive predictor terms of formulaA (via ld()) and combine with the untransformed terms of formulaB into a single formula suitable for glm().

Usage

mvbc.fullformula(formulaA, formulaB, lambdas)

Arguments

Value

a formula object with the positive predictor terms of formulaA replaced by ld(term, lambda) and the terms of formulaB appended

Single train/test glm fit at a fixed lambda

Description

Fits glm() to train at a given set of lambda exponents (weighted, if survey = TRUE and weights is supplied) and, if test is supplied, scores the fit on test and attaches sspr/ssdr (see mvbc.trainer.ssr()) to the returned model object.

Usage

mvbc.glm.step(
  formulaA,
  formulaB = NULL,
  train,
  test = NULL,
  lambdas,
  weights = NULL,
  test_weights = NULL,
  survey = FALSE
)

Arguments

Value

the fitted glm object, with sspr and ssdr attached when test is supplied

Median effects from a fitted BLBC model

Description

Computes coefficient-like median effects for the positive predictors in a model fitted by mvbc.train(). The calculation uses the empirical median, or the sampling-weighted empirical median when weights is supplied, and therefore does not require a log-normal distribution.

Usage

mvbc.median.effect(object, data, q = 1, weights = NULL, na.rm = TRUE)

Arguments

Details

For predictor X_j and reference transformation q, the median effect is

\delta_j(q) = \beta_j m_j^{\lambda_j-q},

where m_j is the marginal median of X_j. The plug-in estimator replaces \beta_j, \lambda_j, and m_j with their fitted or empirical counterparts. If X_j is log-normal with location parameter \mu_j, then m_j = \exp(\mu_j) and the expression reduces to \beta_j\exp\{(\lambda_j-q)\mu_j\}.

With weights, medians are calculated as weighted empirical 0.5 quantiles: values are sorted ascending and the median is the value at which the weighted empirical CDF first reaches 0.5 (equivalent to a one-stage survey design with ids = ~1). The returned effects are on the chosen reference-transformation scale. Comparisons between predictors remain dependent on their measurement units.

Value

A data frame with one row per transformed predictor and columns predictor, beta, lambda, q, median, and median_effect.

Examples

fit <- structure(
  list(
    formulaA = y ~ x1 + x2,
    lambda.fits = data.frame(x1 = 0.5, x2 = 1),
    beta.fits = data.frame(
      foldid = 0L, intercept = -1, x1 = 2, x2 = -0.5
    )
  ),
  class = "mvbc.train.model"
)
dat <- data.frame(
  y = c(0, 0, 1, 1, 1),
  x1 = c(1, 2, 3, 7, 12),
  x2 = c(2, 3, 4, 5, 8)
)
mvbc.median.effect(fit, dat, q = 1)

Optimizer

Description

The Multivariate Boxcox Optimizer accepts an instance of an mvbc.trainer.model and uses Thin Plate Splines (TPS) to search for more optimal lambda estimates.

Usage

mvbc.optimizer(
  model,
  data,
  tpsargs = list(scale.type = "unscaled"),
  optimargs = list(method = "L-BFGS-B", lower = 0, upper = 2),
  weights = NULL,
  survey = FALSE
)

Arguments

Details

The mvbc.optimizer() function performs the following:

1. Starting with an instance of an mvbc.trainer.model repeat the following steps for each minima value in the instance.

2. Generate a Thin Plate Spline (TPS) surface using the lambda grid as the independent x variables in the TPS, and the some-of-square residual minima as the dependent y variable.

3. Use the general purpose optim() function to choose updated estimates of lambda to minimize each of the some-of-square residuals.

4. The updated lambda estimates are used to fit another set of beta values and compute some-of-square residuals as is done in the original mvbc.trainer().

5. Append the updated lambda and beta estimates and residuals to the existing training results. The updated estimates may or may not have improved residuals.

Value

mvbc.optimizer() returns an updated instance of the class mvbc.trainer.model as defined in the original mvbc.trainer().

Examples

vBeta = c(-2.2,-0.4,-0.2,-0.005)
vLambda = c(0.5,1.5)
vMean =  c(-0.08, -0.01, 50)
vSd =    c(0.932543,    0.8,       18.11569)
vNames = c("mercury",   "lead", "age")
sim = mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, 500, 1)
summary(sim)

# ignore "glm.fit: fitted probabilities numerically 0 or 1 occurred"
suppressWarnings({
  trained = mvbc.trainer(
    Ybin ~ mercury + lead, ~ age, sim$data,
    griddomain = seq(0, 2, length.out = 4), K = 3
  )
})
trained = mvbc.optimizer(trained, sim$data)

summary(trained)

Predict probability from a trained mvbc.train.model

Description

Predict \hat{p} on test data using a model fitted by mvbc.train(), i.e. one whose lambda.fits and beta.fits have already been resolved to a single chosen lambda tuple (unlike mvbc.trainer.predict.p_hat(), which predicts from a specific tupleid inside an mvbc.trainer.model that may hold many candidate tuples).

Usage

mvbc.predict(trainModel, data)

Arguments

Details

Given the definition of \bf{Z}(\vec{\lambda},\vec{\beta},\bf{X}) and \hat{p}(\bf{Z}) = 1 / (1 + e^{-\bf{Z}}) as shown mvbc.simulator(), mvbc.predict():

1. Builds the full formula using the chosen \hat{\lambda} in trainModel$lambda.fits.

2. Builds a model matrix from data using that formula.

3. Uses the full-data estimate \hat{\beta} in trainModel$beta.fits to compute \bf{Z}(\hat{\lambda},\hat{\beta},\bf{X}) and \hat{p}(\bf{X}).

4. Returns the estimated probabilities.

Value

mvbc.predict() returns a vector of estimated probabilities, one per row of data.

Examples

vBeta = c(-2.2,-0.4,-0.2,-0.005)
vLambda = c(0.5,1.5)
vMean =  c(-0.08, -0.01, 50)
vSd =    c(0.932543,    0.8,       18.11569)
vNames = c("mercury",   "lead", "age")
trainset = mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, 500, 1)
testset = mvbc.simulator(simModel = trainset)

suppressWarnings({
  trained = mvbc.train(
    Ybin ~ mercury + lead, ~ age, trainset$data,
    griddomain = seq(0, 2, length.out = 4), K = 3
  )
})
p_hat = mvbc.predict(trained, testset$data)
head(p_hat)

Simulator

Description

The Multivariate Boxcox Simulator mvbc.simulator generates binomial samples from a logistic model built from both Boxcox transformed log normal random variables and normal random variables. The simulator is given a set of model parameters, and returns a class object with the parameters, the random variables, the logistic output \bf{Z} and the binomial samples \bf{Y_{bin}}. The simulator includes controlled seeding to allow for generation of the same sequence of simulated data each time.

Also read the introduction.

Usage

mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, n, seed, simModel = NULL)

Arguments

Details

The basic Multivariate Boxcox logistic model is expressed as:

\bf{Z}(\vec{\lambda},\vec{\beta},\bf{X}) = \beta_0 + \beta_1 \bf{X_1}^{(\lambda_1)} + ... + \beta_c \bf{X_c}^{(\lambda_c)} + \beta_{c+1} \bf{X_{c+1}} + ... + \beta_{c+k} \bf{X_{c+k}}

where the probability of a positive event p = E(y = 1|\bf{X}) is estimated as:

\bf{\hat{p}}(\bf{Z}) = 1 / (1 + e^{-\bf{Z}})

and the Boxcox transformation x^{(\lambda)} is:

• (x^{\lambda} - 1)/{\lambda} \text{, if} \ \lambda \neq 0

• ln(x) \text{, if} \ \lambda = 0

The contributing variables are:

• \bf{X_a} are the log normal random variables where 1 < a < c

• \bf{X_b} are the normal random variables where c+1 < b < c+k

• c is the number of log normal random variables

• k is the number of normal random variables

• \vec{\beta} is a vector of the logistic model coefficients

• \vec{\mu} is a vector of the means of the random variables \bf{X_a}, \bf{X_b}

• \vec{\sigma} is a vector of the standard deviations of the random variables \bf{X_a}, \bf{X_b}

• \vec{\lambda} is a vector of Boxcox exponents for the log normal random variables \bf{X_a}

The simulator will generate n samples for \bf{X_a} and \bf{X_b} and then compute the vector \bf{Z}. The samples are generated using the R functions rlnorm(n, \mu, \sigma) and rnorm(n, \mu, \sigma) respectively. The vector \bf{Z} is then used to generate n binomial samples \bf{Y_{bin}} using the R function:

\bf{Y_{bin}} = rbinom(n, 1, 1 / (1 + e^{-\bf{Z}}))

An initial first usage call to mvbc.simulator() will set the given seed using set.seed(). Subsequent calls using the object will increment the seed first, and then call set.seed().

Value

mvbc.simulator() returns an object of class mvbc.simulator.model which contains the following components:

Examples

vBeta = c(-2.2,-0.4,-0.2,-0.005)
vLambda = c(0.5,1.5)
vMean =  c(-0.08, -0.01, 50)
vSd =    c(0.932543,    0.8,       18.11569)
vNames = c("mercury",   "lead", "age")
sim = mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, 1000, 1)
summary(sim)

Tester

Description

The Multivariate Boxcox Tester provides a means to predict binary outcomes for a test data set, given a trained model from the Multivariate Boxcox Trainer Core. The Tester can also run multiple train and test cycles using the Multivariate Boxcox Simulator to generate matching training data and test data. A variety of evaluation metrics are available including the Pearson and Deviance residuals and the F1-score and related metrics.

Also read the introduction.

Usage

mvbc.tester(
  formulaA,
  formulaB = NULL,
  simModel,
  nTrain,
  nTest,
  griddomain = seq(0, 2, l = 5),
  K = 5,
  depth = 2,
  cores = 1,
  minima = c("ssdr"),
  tpsargs = list(scale.type = "unscaled"),
  optimargs = list(method = "L-BFGS-B", lower = 0, upper = 2)
)

Arguments

Details

The mvbc.tester() function performs the following:

1. Using simModel as a basis, call mvbc.simulator() to generate nTrain training data sets and nTest test data sets. Relies on the embedded seed to generate repeatable results.

2. For each of the nTrain training data sets, fit an mvbc.trainer.model mvbc.trainer() and then optimize the fit using mvbc.optimizer().

3. For the given minima:

   1. For each nTest test data set, estimate \hat{p} using mvbc.trainer.predict.p_hat() then estimate \hat{\bf{Y}}_{bin}. Compute sum of square Pearson and Deviance residuals, and a confusion matrix report cm. Save results in an instance of a mvbc.tester.result

4. Return an instance of mvbc.tester.result.

When estimating \hat{\bf{Y}}_{bin} for use in the confusion matrix, the trainModel$prior will be used to define the cut-off for \hat{p}.

Value

mvbc.tester() returns an object of class mvbc.tester.result which contains the following components:

Examples

vBeta = c(-2.2,-0.4,-0.2,-0.005)
vLambda = c(0.45,1.50)
vMean =  c(-0.08, -0.01, 50)
vSd =    c(0.932543,    0.8,       18.11569)
vNames = c("mercury",   "lead", "age")
sim = mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, 300, 1)

# one training sample and one test sample
# ignore "glm.fit: fitted probabilities numerically 0 or 1 occurred"
suppressWarnings({
    test.result = mvbc.tester(
      Ybin ~ mercury + lead, ~ age, sim,
      nTrain = 1, nTest = 1,
      griddomain = seq(0, 2, length.out = 3),
      K = 2, depth = 1
    )
})
summary(test.result)

maxF1 = which.max(test.result$lambda.fits$F1)
test.result$cms[maxF1]

Train

Description

Full pipeline: runs mvbc.trainer() to grid-search candidate lambdas, then mvbc.optimizer() to refine the selected minima with a Thin Plate Spline, refits the regression coefficients on the full data set at the selected lambda, and returns a fitted mvbc.train.model with lambda.fits and beta.fits populated at the chosen minima.

Usage

mvbc.train(
  formulaA,
  formulaB = NULL,
  data,
  griddomain = seq(0, 2, l = 10),
  K = 5,
  depth = 1,
  cores = 1,
  minima = c("ssdr"),
  tpsargs = list(scale.type = "unscaled"),
  optimargs = list(method = "L-BFGS-B", lower = 0, upper = 2),
  seed = 111,
  weights = NULL,
  survey = FALSE
)

Arguments

Value

mvbc.train() returns an object of class mvbc.train.model, which extends mvbc.trainer.model (see mvbc.trainer()) with lambda.fits and beta.fits set to the estimates at the chosen minima. beta.fits contains the single full-data refit, identified by foldid = 0.

Examples

vBeta = c(-2.2,-0.4,-0.2,-0.005)
vLambda = c(0.5,1.5)
vMean =  c(-0.08, -0.01, 50)
vSd =    c(0.932543,    0.8,       18.11569)
vNames = c("mercury",   "lead", "age")
sim = mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, 500, 1)

# ignore "glm.fit: fitted probabilities numerically 0 or 1 occurred"
suppressWarnings({
  trained = mvbc.train(
    Ybin ~ mercury + lead, ~ age, sim$data,
    griddomain = seq(0, 2, length.out = 4), K = 3, depth = 1
  )
})
summary(trained)

Trainer

Description

The Multivariate Box-Cox Trainer Core mvbc.trainer accepts positive continuous predictors \bf{X_a} to be Box-Cox transformed, optional untransformed covariates \bf{X_b}, and a binary outcome \bf{Y_{bin}}. Using K-fold cross-validation, the trainer estimates logistic models for \bf{Z}(\vec{\lambda},\vec{\beta},\bf{X}) consisting of the estimated Box-Cox exponents \hat{\lambda} and K sets of estimated coefficients \hat{\beta}.

Also read the introduction.

Usage

mvbc.trainer(
  formulaA,
  formulaB = NULL,
  data,
  griddomain = seq(0, 2, l = 10),
  K = 5,
  depth = 1,
  cores = 1,
  minima = c("ssdr"),
  seed = 111,
  weights = NULL,
  survey = FALSE
)

Arguments

Details

Recall from mvbc.simulator() that the basic Multivariate Boxcox logistic model is expressed as:

\bf{Z}(\vec{\lambda},\vec{\beta},\bf{X}) = \beta_0 + \beta_1 \bf{X_1}^{(\lambda_1)} + ... + \beta_c \bf{X_c}^{(\lambda_c)} + \beta_{c+1} \bf{X_{c+1}} + ... + \beta_{c+k} \bf{X_{c+k}}

where:

• c is the number of positive predictors to be transformed

• k is the number of additional untransformed covariates

• x^{(\lambda)} is the Box-Cox transformation

• \hat{p}(\bf{Z}) = 1 / (1 + e^{-\bf{Z}}) is the estimated probability of a positive event

The mvbc.trainer() algorithm is summarized here:

1. Creates a Lambda grid matrix from the griddomain. The grid matrix contains all the possible candidate tuples of Lambda, with additional columns sspr and ssdr to contain sum-of-squares of Pearson and Deviance residuals. For example, if each lambda domain contains 10 values, then there would be G = 10^c tuples.

2. Split the data into K-folds (reuse folds if available in object$folds).

3. Create a gridfits matrix where each row will store tupleid, foldid, and the \hat{\beta} estimates from a glm() model.

4. For each fold, create training-set and test-set...

   1. For each lambda tupleid in grid...

      1. Fit \bf{Z}(\vec{\lambda},\vec{\beta},\bf{X}) to the training-set using glm().

      2. Use the glm model to estimate probabilities \bf{\hat{p}} from the test-set.

      3. Use \bf{Y}_{bin} and \bf{\hat{p}} from the test-set to compute Pearson and Deviance residuals and accumulate in grid.

      4. Save the glm model coefficients as the \hat{\beta} estimates in the gridfits matrix.

5. Find lambda tupleids for the minimum sspr and ssdr and save in minimas list.

6. While depth > 1 then:

   1. Decrement depth.

   2. Compute a finer griddomain centered on chosen minima. Repeat steps 3 to 5 and merge the the grid, gridfits, and minimas results.

Value

mvbc.trainer() returns an object of class mvbc.trainer.model which contains the following components:

Examples

vBeta = c(-2.2,-0.4,-0.2,-0.005)
vLambda = c(0.5,1.5)
vMean =  c(-0.08, -0.01, 50)
vSd =    c(0.932543,    0.8,       18.11569)
vNames = c("mercury",   "lead", "age")
sim = mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, 500, 1)
summary(sim)

# ignore "glm.fit: fitted probabilities numerically 0 or 1 occurred"
suppressWarnings({
  trained = mvbc.trainer(
    Ybin ~ mercury + lead, ~ age, sim$data,
    griddomain = seq(0, 2, length.out = 4), K = 3
  )
})
summary(trained)

Workhorse to fit instance of mvbc.trainer.model

Description

Workhorse to fit instance of mvbc.trainer.model

Usage

mvbc.trainer.fit(
  trainModel,
  data,
  gridid.start = 1,
  weights = NULL,
  survey = FALSE
)

Arguments

Value

An updated object of class mvbc.trainer.model. Its grid component contains the accumulated cross-validation residual criteria for each lambda tuple, and gridfits contains the corresponding fold-specific regression coefficients. All other training metadata in trainModel is retained.

Predict outcomes

Description

Predict \hat{Y_{bin}} outcomes from training model

Usage

mvbc.trainer.predict.Ybin_hat(prior, p_hat)

Arguments

Value

integer vector of predicted binary outcomes (0/1)

Predict probability

Description

Predict \hat{p} on test data using training model

Usage

mvbc.trainer.predict.p_hat(trainModel, data, tupleid)

Arguments

Details

Given the definition of \bf{Z}(\vec{\lambda},\vec{\beta},\bf{X}) and \hat{p}(\bf{Z}) = 1 / (1 + e^{-\bf{Z}}) as shown mvbc.simulator(), the mvbc.trainer.predict.p_hat() function performs the following:

1. Look up the estimated \hat{\lambda} values associated with the given tupleid.

2. Look up the K sets of estimated \hat{\beta} values associated with the given tupleid.

3. For each of the K sets of \hat{\beta} values, compute an estimated value of \bf{Z}(\hat{\lambda},\hat{\beta},\bf{X}) and \hat{p}(\bf{X}) from the test data

4. Compute the mean of the K estimates of \hat{p}(\bf{Z}) as \hat{p}_{mean}(\bf{Z}) and return this vector of mean \hat{p} values.

Value

mvbc.trainer.predict.p_hat() returns vector of \hat{p_{mean}}(\bf{Z}) estimates.

Examples

vBeta = c(-2.2,-0.4,-0.2,-0.005)
vLambda = c(0.5,1.5)
vMean =  c(-0.08, -0.01, 50)
vSd =    c(0.932543,    0.8,       18.11569)
vNames = c("mercury",   "lead", "age")
trainset = mvbc.simulator(vLambda, vBeta, vMean, vSd, vNames, 500, 1)
testset = mvbc.simulator(simModel = trainset)

# train and predict one pair
suppressWarnings({
  trained = mvbc.trainer(
    Ybin ~ mercury + lead, ~ age, trainset$data,
    griddomain = seq(0, 2, length.out = 4), K = 3
  )
})
p_hat = mvbc.trainer.predict.p_hat(trained, testset$data, trained$minimas[["ssdr"]])
head(p_hat)

Sum of Squares residuals

Description

Return sum-of-squares of Pearson and Deviance Residuals

Usage

mvbc.trainer.ssr(Ybin, p_hat, weights = NULL)

Arguments

Value

named numeric vector c(sspr, ssdr) with the sum-of-squares of the Pearson and Deviance residuals

Plot instance of mvbc.sim.result

Description

Plot instance of mvbc.sim.result

Usage

## S3 method for class 'mvbc.sim.result'
plot(x, ...)

Arguments

Value

No return value, called for side effects

Plot instance of mvbc.tester.result

Description

Plot instance of mvbc.tester.result

Usage

## S3 method for class 'mvbc.tester.result'
plot(x, ...)

Arguments

Value

No return value, called for side effects

Summarize instance of mvbc.sim.result

Description

Summarize instance of mvbc.sim.result

Usage

## S3 method for class 'mvbc.sim.result'
summary(object, ...)

Arguments

Value

No return value, called for side effects

Summarize instance of mvbc.simulator.model

Description

Summarize instance of mvbc.simulator.model

Usage

## S3 method for class 'mvbc.simulator.model'
summary(object, ...)

Arguments

Value

No return value, called for side effects

Summarize instance of mvbc.tester.result

Description

Summarize instance of mvbc.tester.result

Usage

## S3 method for class 'mvbc.tester.result'
summary(object, ...)

Arguments

Value

No return value, called for side effects

Summarize instance of mvbc.trainer.model

Description

Summarize instance of mvbc.trainer.model

Usage

## S3 method for class 'mvbc.trainer.model'
summary(object, ...)

Arguments

Value

No return value, called for side effects