plsi.log.auto.constraint: Partial linear single index log-linear (count) regression...

In EPLSIM: Partial Linear Single Index Models for Environmental Mixture Analysis

Partial linear single index log-linear (count) regression with automatic smoothness selection and built-in identifiability constraint

Description

Count-outcome analog of plsi.lr.auto.constraint(): the same model as plsi.log.auto() (penalized regression spline link via mgcv::gam with a log link, negative binomial or Poisson family, REML smoothing-parameter selection), but with the single-index direction searched over a reparameterization that satisfies the identifiability constraint by construction instead of normalizing after the fact.

Usage

plsi.log.auto.constraint(
  data,
  Y.name,
  X.name,
  Z.name,
  family = c("nb", "poisson"),
  k = 10,
  bs = "cr",
  initial.random.num = 5,
  seed = 2026,
  maxit = 100
)

Arguments

data	A data set including all needed variables
Y.name	Variable name for the count outcome. Must be non-negative integers (0, 1, 2, ...).
X.name	Variable name vector for exposures (at least 2).
Z.name	Variable name vector for confounders
family	Count distribution to use. "nb" (default) fits a negative binomial GAM (mgcv::nb()), which estimates its own dispersion parameter and is the safer default for real, usually overdispersed, count data. "poisson" fits a Poisson GAM, valid only if the mean-equals-variance assumption holds.
k	Upper bound on basis dimension for the smooth link function. Default 10.
bs	Smooth basis type passed to mgcv::s(); default "cr" (cubic regression spline).
initial.random.num	Number of random initials for the single-index direction search
seed	A single integer used to seed the RNG so results (random initials and downstream optim() convergence) are reproducible. Set to NULL to skip seeding and use whatever RNG state is already active in the calling environment. Default is 2026.
maxit	Maximum number of Nelder-Mead iterations per starting value for the \sigma-space search. Default 100, matching plsi.lr.auto.constraint(). With many exposures this is often too few for the simplex to converge, which can leave the numerical information matrix indefinite and the single-index SEs inflated; check optim.convergence (and the warning issued when the selected start did not converge) and increase maxit if needed. Ignored when there are exactly two exposures (Brent search).

Details

The single-index coefficient \beta is constrained to \|\beta\| = 1 and \beta_1 > 0 for identifiability. The search is reparameterized in terms of \sigma = (\sigma_1, \dots, \sigma_{p-1}) via \beta(\sigma) = \{(1 - \|\sigma\|^2)^{1/2}, \sigma_1, \dots, \sigma_{p-1}\}^T. Both constraints hold automatically for any \sigma with \|\sigma\| < 1, so optim() runs unconstrained over \sigma. Standard errors for \beta are obtained from the Hessian in \sigma-space via the delta method.

The anchor exposure, the one assigned \beta_1, the coordinate constrained > 0, is chosen as the candidate exposure with the largest absolute correlation with the count outcome, so an exposure strongly but negatively associated with the outcome can still be selected as the anchor. Which exposure was selected is reported in attr(si.coefficient, "note") and in anchor.note of the returned list.

The confounder-free single-index model is built exactly as in plsi.log.auto(): the confounder contribution from the full model is fixed as a known offset() and a second count gam (same family) is refit directly against the observed y.

Standard errors are computed from a numerical Hessian in a local chart of the unit sphere anchored on the exposure with the largest |\beta_j| (not necessarily the \beta_1 > 0 anchor), then mapped to \beta by the delta method. This keeps the finite-difference steps away from the \|\sigma\| = 1 boundary when the estimated \beta_1 is small; the point estimates are unaffected.

With exactly two exposures, \sigma is one-dimensional; the search then uses optim(method = "Brent") on (-1, 1) instead of Nelder-Mead, which optim() flags as unreliable in one dimension.

The returned object works with all plotting functions using type = "log".

Value

A list of model estimation and prediction results, structured analogously to plsi.log.auto()'s output, plus:

• anchor.note: Which exposure was used as the anchor (\beta_1 > 0) and why; also attached as attr(si.coefficient, "note").

• initial.sigma.table: The starting values actually used by the \sigma-space search.

• optim.convergence: For every starting value, the optim() convergence code (0 = converged, 1 = maxit reached), the number of log-likelihood evaluations, and the final log-likelihood.

The remaining elements are as in plsi.log.auto():

• si.coefficient: Single-index direction estimates and inference.

• confounder.coefficient: Confounder log-rate coefficients, plus rate ratios and their 95% CIs.

• si.fun: The estimated single-index link function, on both the log scale (fit/se/lwr/upr) and the count scale (count.fit/count.lwr/count.upr).

• si.fun.model: A confounder-free count gam on the single index alone (confounder contribution held as a fixed offset); predict() from this needs only single_index_estimated (and confounder_offset = 0).

• full.model: The full fitted gam (smooth + confounders).

Author(s)

Yuyan Wang

Examples

## Not run: 
data(nhanes.new)
data <- nhanes.new

X.name <- c("X1_trans.b.carotene", "X2_retinol", "X3_g.tocopherol", "X4_a.tocopherol",
            "X5_PCB99", "X6_PCB156", "X7_PCB206",
            "X8_3.3.4.4.5.pncb", "X9_1.2.3.4.7.8.hxcdf", "X10_2.3.4.6.7.8.hxcdf")
Z.name <- c("AGE.c", "SEX.Female", "RACE.NH.Black",
           "RACE.MexicanAmerican", "RACE.OtherRace", "RACE.Hispanic" )

# demo count outcome (illustrative only as nhanes.new has no native count
# variable), simulated from a true single-index combination of the exposures
set.seed(2026)
beta_true <- c(0.30, -0.20, 0.10, 0.40, -0.30, 0.20, -0.10, 0.25, -0.15, 0.35)
beta_true <- beta_true / sqrt(sum(beta_true^2))
x_std <- scale(data[, X.name])
single_index_true <- as.vector(x_std %*% beta_true)
log_rate <- 0.3 + 0.4 * sin(single_index_true) + 0.05 * data$AGE.c
data$n.events <- stats::rpois(nrow(data), lambda = exp(log_rate))

model_log_con <- plsi.log.auto.constraint(data = data, Y.name = "n.events",
                   X.name = X.name, Z.name = Z.name, family = "nb",
                   k = 10, bs = "cr", initial.random.num = 1, seed = 2026)
model_log_con$si.coefficient
si.fun.plot(model_log_con$si.fun, type = "log")

## End(Not run)

EPLSIM documentation built on Nov. 5, 2026, 1:08 a.m.