An Introduction to EPLSIM Package

Introduction

EPLSIM fits partial linear single index models (PLSIM) for evaluating the health effects of environmental mixtures. The model combines a set of exposures into a single index via a linear combination, relates that index to the outcome through a flexible (possibly nonlinear) link function, and adjusts for confounders through a separate linear term. This vignette walks through a typical workflow: preparing confounders, fitting a model, and visualizing the results.

Three families of fitting functions are available:

plsi.lr.auto.constraint(), plsi.logistic.auto.constraint(), plsi.log.auto.constraint(): recommended for most users, for continuous, binary, and count outcomes respectively. The link function is estimated with mgcv::gam() using REML-based automatic smoothness selection, so you only supply an upper bound k on basis complexity. The single-index direction is searched over a reparameterization that satisfies the identifiability constraint (||beta|| = 1, beta_1 > 0) by construction, rather than normalizing and sign-flipping the estimate after the search.
plsi.lr.auto(), plsi.logistic.auto(), plsi.log.auto(): the same models, with the direction searched over unconstrained beta and normalized afterwards.
plsi.lr.v2(): retained for cases where you want direct control over the link function’s degrees of freedom (e.g. simulation studies calling the fitting function many times, or when you already know the right basis complexity from prior work). It uses a fixed-df natural cubic spline (splines::ns()) instead of a penalized smooth.

This vignette focuses on the continuous-outcome case (plsi.lr.auto.constraint()); the binary- and count-outcome analogs follow the same interface and are introduced at the end.

Preparing confounders

confounder.trans() standardizes continuous confounders and creates dummy variables for categorical ones, returning both the updated data and the resulting variable names to use as Z.name.

dat.cov <- data.frame(
  age = c(1.5, 2.3, 3.1, 4.8, 5.2),
  sex = c(1, 2, 1, 2, 2),
  race = c(1, 2, 3, 4, 5)
)

# specify the confounder vector
Z.name <- c("age", "sex", "race")

# set levels and make the reference level first for categorical confounders
dat.cov$sex <- factor(dat.cov$sex, 1:2, c('Male', 'Female'))
dat.cov$race <- factor(dat.cov$race, 1:5, c("NH-White", "NH-Black",
                                             "MexicanAmerican", "OtherRace", "Hispanic"))

# transform the confounder vector and check
cov_m <- confounder.trans(Z_continuous = c("age"), Z_discrete = c("sex", "race"), data = dat.cov)
Z.name <- cov_m$New.Name
dat.cov <- cov_m$Updated.data
print(Z.name)
#> [1] "age.c"                "sex.Female"           "race.NH-Black"       
#> [4] "race.MexicanAmerican" "race.OtherRace"       "race.Hispanic"

The package’s built-in example data, nhanes.new, already has its confounders prepared this way, so the rest of this vignette uses it directly.

Fitting the model

data(nhanes.new)
dat <- nhanes.new

# outcome
Y.name <- "log.triglyceride"

# exposures making up the single index
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")

# confounders, already prepared as in the section above
Z.name <- c("AGE.c", "SEX.Female", "RACE.NH.Black",
            "RACE.MexicanAmerican", "RACE.OtherRace", "RACE.Hispanic")

Recommended: `plsi.lr.auto.constraint()`

model_con <- plsi.lr.auto.constraint(data = dat, Y.name = Y.name,
                                     X.name = X.name, Z.name = Z.name,
                                     k = 10, bs = "cr",
                                     initial.random.num = 1, seed = 2026)

# which exposure anchors the index (its coefficient is constrained > 0)
model_con$anchor.note
#> [1] "First/anchor exposure (single-index coefficient constrained > 0 for identifiability, ||beta|| = 1): 'X4_a.tocopherol', selected as the candidate exposure with the largest absolute correlation with log.triglyceride (r = 0.564)."

# effective degrees of freedom REML actually used for the link function
# (will be <= k; substantially less than k indicates the relationship is
# close to linear)
model_con$si.fun.edf
#> [1] 3.17364

# single-index coefficients: which exposures drive the index, and how much
model_con$si.coefficient
#>                          Estimate   Std.Error    t value Pr(>|t|)  Lower.95CI
#> X4_a.tocopherol        0.75545431 0.029802790 25.3484427   <.0001  0.69704191
#> X3_g.tocopherol        0.47331321 0.032249099 14.6767887   <.0001  0.41010614
#> X2_retinol             0.29191585 0.044625336  6.5414824   <.0001  0.20445180
#> X8_3.3.4.4.5.pncb      0.14910275 0.048192599  3.0938930   0.0020  0.05464699
#> X5_PCB99               0.01226986 0.073151473  0.1677323   0.8668 -0.13110439
#> X1_trans.b.carotene   -0.01722678 0.003618891 -4.7602375   <.0001 -0.02431968
#> X9_1.2.3.4.7.8.hxcdf  -0.08432514 0.045943063 -1.8354270   0.0664 -0.17437189
#> X6_PCB156             -0.09419057 0.075006809 -1.2557602   0.2092 -0.24120121
#> X7_PCB206             -0.16927863 0.072098731 -2.3478725   0.0189 -0.31058955
#> X10_2.3.4.6.7.8.hxcdf -0.22963393 0.040029701 -5.7365888   <.0001 -0.30809070
#>                        Upper.95CI Contribution proportion
#> X4_a.tocopherol        0.81386670                   0.571
#> X3_g.tocopherol        0.53652029                   0.224
#> X2_retinol             0.37937990                   0.085
#> X8_3.3.4.4.5.pncb      0.24355850                   0.022
#> X5_PCB99               0.15564412                   0.000
#> X1_trans.b.carotene   -0.01013388                   0.000
#> X9_1.2.3.4.7.8.hxcdf   0.00572161                   0.007
#> X6_PCB156              0.05282008                   0.009
#> X7_PCB206             -0.02796771                   0.029
#> X10_2.3.4.6.7.8.hxcdf -0.15117716                   0.053

# confounder coefficients
model_con$confounder.coefficient
#>                        Estimate Std. Error   t value Pr(>|t|)  Lower.95CI
#> AGE.c                 0.1163091 0.02532260  4.593095   <.0001  0.06667772
#> SEX.Female           -0.1456704 0.04663300 -3.123762   0.0019 -0.23706938
#> RACE.NH.Black        -0.1377246 0.06158133 -2.236466   0.0256 -0.25842174
#> RACE.MexicanAmerican  0.1002847 0.06023325  1.664939   0.0963 -0.01777031
#> RACE.OtherRace        0.3213728 0.13420315  2.394674   0.0169  0.05833945
#> RACE.Hispanic         0.5185808 0.14021225  3.698541   0.0002  0.24376985
#>                       Upper.95CI
#> AGE.c                 0.16594050
#> SEX.Female           -0.05427137
#> RACE.NH.Black        -0.01702736
#> RACE.MexicanAmerican  0.21833968
#> RACE.OtherRace        0.58440614
#> RACE.Hispanic         0.79339179

k is a ceiling, not a target: REML shrinks the fitted smooth toward a straight line wherever the data don’t support extra wiggliness, so setting k generously (e.g. 10) and letting REML choose is usually preferable to guessing a lower value up front.

Identifiability and the anchor exposure

A single-index model is only identified up to the scale and sign of beta: multiplying beta by -1 and reflecting the link function gives the same fit. The constrained functions fix this by writing

beta(sigma) = (sqrt(1 - ||sigma||^2), sigma_1, ..., sigma_{p-1})

and optimizing over sigma. Any sigma with ||sigma|| < 1 gives a beta with unit length and a positive first coordinate, so the constraint never has to be imposed after the fact. Standard errors for beta are obtained from the numerical Hessian and mapped back to beta by the delta method.

The first coordinate, the anchor, is the exposure with the largest absolute correlation with the outcome. An exposure strongly but negatively associated with the outcome can therefore be the anchor; its sign is absorbed by the beta_1 > 0 convention, so the signs of the other coefficients should be read relative to it.

Other fitting functions

plsi.lr.auto() fits the same model without the built-in constraint, and plsi.lr.v2() uses a fixed-df natural cubic spline in place of the penalized smooth, which is useful when running many fits in a simulation where mgcv::gam()’s per-call overhead adds up:

model_auto <- plsi.lr.auto(data = dat, Y.name = Y.name, X.name = X.name, Z.name = Z.name,
                           k = 10, bs = "cr", initial.random.num = 1, seed = 2026)

model_v2 <- plsi.lr.v2(data = dat, Y.name = Y.name, X.name = X.name, Z.name = Z.name,
                       spline.num = 5, spline.degree = 3, initial.random.num = 1, seed = 2026)

Every plotting function below works identically with the output of any of these fitting functions.

Visualizing results

Single-index coefficients

si.coef.plot(model_con$si.coefficient)

Bars show each exposure’s contribution to the single index, with 95% confidence intervals; the color indicates the sign of the coefficient.

Single-index function

The single-index function shows the estimated (possibly nonlinear) relationship between the combined exposure index and the outcome.

si.fun.plot(model_con$si.fun, type = "linear")

One exposure’s main effect

e.main.plot() shows the predicted outcome as a single exposure varies, holding the rest of the index at its confounder-adjusted reference level.

e.main.plot(model_con, dat, exp_name = "X4_a.tocopherol", type = "linear")

e.main.plot(model_con, dat, exp_name = "X5_PCB99", type = "linear")

e.main.plot(model_con, dat, exp_name = "X10_2.3.4.6.7.8.hxcdf", type = "linear")

Interaction between two exposures

e.interaction.plot() shows one exposure’s effect at the 25th/50th/75th percentile of a second exposure, and vice versa, side by side.

e.interaction.plot(model_con, dat, "X4_a.tocopherol", "X3_g.tocopherol", type = "linear")

The two panels are not interchangeable: the left panel varies the first named exposure while holding the second at its quantiles, and the right panel does the reverse. Swapping the argument order swaps which exposure is varied in which panel:

e.interaction.plot(model_con, dat, "X3_g.tocopherol", "X4_a.tocopherol", type = "linear")

Interquartile effects across all exposures

interquartile.quartile.plot() summarizes, for every exposure at once, the difference in predicted outcome between its 75th and 25th percentile, at the 25th/50th/75th percentile of every other exposure.

interquartile.quartile.plot(model_con, dat, type = "linear")

Mixture overall effect

mixture.overall.plot() shows the predicted outcome as all exposures move together through their joint quantiles (e.g. the 10th percentile of every exposure simultaneously, the 15th percentile of every exposure simultaneously, and so on).

mixture.overall.plot(model_con, dat, type = "linear")

Binary and count outcomes

plsi.logistic.auto.constraint() and plsi.log.auto.constraint() follow the same interface as plsi.lr.auto.constraint(), and every plotting function above accepts a matching type = "logistic" or type = "log" argument. nhanes.new has no native binary or count outcome, so the examples below simulate one from a known 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))
si_true <- as.vector(scale(dat[, X.name]) %*% beta_true)

Binary outcome

dat$high.tg <- rbinom(nrow(dat), size = 1,
                      prob = plogis(-0.2 + 0.5 * sin(si_true) + 0.05 * dat$AGE.c))

model_logistic <- plsi.logistic.auto.constraint(data = dat, Y.name = "high.tg",
                                                X.name = X.name, Z.name = Z.name,
                                                k = 10, bs = "cr",
                                                initial.random.num = 1, seed = 2026,
                                                maxit = 300)

# did the selected starting value converge?
model_logistic$optim.convergence

model_logistic$si.coefficient
model_logistic$confounder.coefficient   # includes odds ratios

si.fun.plot(model_logistic$si.fun, type = "logistic")
e.main.plot(model_logistic, dat, exp_name = "X4_a.tocopherol", type = "logistic")

Count outcome

dat$n.events <- rpois(nrow(dat), lambda = exp(0.3 + 0.4 * sin(si_true) + 0.05 * dat$AGE.c))

model_log <- plsi.log.auto.constraint(data = dat, Y.name = "n.events",
                                      X.name = X.name, Z.name = Z.name,
                                      family = "nb", k = 10, bs = "cr",
                                      initial.random.num = 1, seed = 2026,
                                      maxit = 300)

model_log$optim.convergence
model_log$si.coefficient
model_log$confounder.coefficient        # includes rate ratios

si.fun.plot(model_log$si.fun, type = "log")
e.main.plot(model_log, dat, exp_name = "X4_a.tocopherol", type = "log")

family = "nb" (negative binomial, the default) estimates its own dispersion parameter and is the safer choice for real, usually overdispersed, count data; use family = "poisson" only if the mean-equals-variance assumption holds.

For both outcome types, maxit sets the number of Nelder-Mead iterations per starting value. With many exposures the default of 100 may be too few for the search to converge; check optim.convergence and increase maxit if the selected start reports a convergence code of 1.

An Introduction to EPLSIM Package

Yuyan Wang

2026-11-04

Introduction

Preparing confounders

Fitting the model

Recommended: `plsi.lr.auto.constraint()`

Identifiability and the anchor exposure

Other fitting functions

Visualizing results

Single-index coefficients

Single-index function

One exposure’s main effect

Interaction between two exposures

Interquartile effects across all exposures

Mixture overall effect

Binary and count outcomes

Binary outcome

Count outcome

An Introduction to EPLSIM Package

Yuyan Wang

2026-11-04

Introduction

Preparing confounders

Fitting the model

Recommended: plsi.lr.auto.constraint()

Identifiability and the anchor exposure

Other fitting functions

Visualizing results

Single-index coefficients

Single-index function

One exposure’s main effect

Interaction between two exposures

Interquartile effects across all exposures

Mixture overall effect

Binary and count outcomes

Binary outcome

Count outcome

Recommended: `plsi.lr.auto.constraint()`