plsi.log.auto: Partial linear single index log-linear regression with...

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

Partial linear single index log-linear regression with automatic smoothness selection

Description

Count-outcome analog of plsi.lr.auto() / plsi.logistic.auto(): same partial linear single index model, but for non-negative count data via a log link. The link function is estimated with a penalized regression spline (mgcv::gam, REML smoothing-parameter selection), so there is no spline.num to hand-tune, only an upper bound k on basis complexity, with REML shrinking unneeded wiggliness toward a straight line on the log scale.

Usage

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

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
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 count data, which is usually overdispersed (variance > mean) relative to Poisson. "poisson" fits a Poisson GAM, valid only if the mean-equals- variance assumption actually holds; if unsure, leave the default.
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.

Details

As with plsi.logistic.auto(), the "confounder-free" single-index model cannot be built by subtracting confounder effects from y on the response scale (counts aren't additive that way either), and it cannot be built by residualizing the full model's own fitted linear predictor, that would leave a value that is already an exact smooth function of the single index, so refitting a smooth to it would just retrace the same curve with near-zero (and spuriously tiny) standard errors. Instead, the confounder contribution from the full model is fixed as a known offset() and a second, confounder-free count GAM (same family as the full model) is refit directly against the actual observed y, preserving real count sampling variability. Predictions from this model are made at offset = 0, consistent with how the offset was netted out of every training point.

Value

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

• 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), obtained by exponentiating the log-scale CI so it stays non-negative.

• si.fun.model: A confounder-free gam of the count outcome on the single index alone (via a fixed offset for confounder contribution), predict() from this needs only single_index_estimated, not the original confounders. See Details.

• full.model: The full fitted gam (smooth + confounders), kept for reference/diagnostics.

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, run through a
nonlinear log-rate link, plus a confounder effect, so the example actually has an
exposure-outcome relationship for plsi.log.auto() to recover, rather than
depending on the confounder alone.

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)) # unit norm, matching the model's
# own identifiability constraint
x_std <- scale(data[, X.name]) # standardize so the demo index is on a sane, comparable scale
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))

Y.name <- "n.events"

k <- 10
bs <- "cr"
initial.random.num <- 1
seed <- 2026

model_log_auto <- plsi.log.auto(data = data, Y.name = Y.name, X.name = X.name,
                      Z.name = Z.name, family = "nb", k = k, bs = bs,
                      initial.random.num = initial.random.num, seed = seed)

## End(Not run)

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