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.

Two families of fitting functions are available:

plsi.lr.auto(), plsi.logistic.auto(), plsi.log.auto() — recommended for most users. These use mgcv::gam() with REML-based automatic smoothness selection for the link function, so you only need to supply an upper bound k on basis complexity rather than tune exact degrees of freedom.
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() / plsi.lr.v2()); the binary- and count-outcome analogs (plsi.logistic.auto(), plsi.log.auto()) follow the same interface and are introduced briefly 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()`

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)

# 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_auto$si.fun.edf
#> [1] 3.34882

# single-index coefficients: which exposures drive the index, and how much
model_auto$si.coefficient
#>                          Estimate   Std.Error    t value Pr(>|t|)  Lower.95CI
#> X4_a.tocopherol        0.75566630 0.131081926  5.7648398   <.0001  0.49875045
#> X3_g.tocopherol        0.46466249 0.086207601  5.3900408   <.0001  0.29569869
#> X2_retinol             0.26633906 0.059855598  4.4496934   <.0001  0.14902425
#> X8_3.3.4.4.5.pncb      0.13655925 0.052184072  2.6168761   0.0089  0.03428035
#> X5_PCB99               0.01160050 0.056452956  0.2054897   0.8372 -0.09904526
#> X1_trans.b.carotene   -0.01785651 0.004643782 -3.8452514   0.0001 -0.02695815
#> X9_1.2.3.4.7.8.hxcdf  -0.09325569 0.050205380 -1.8574840   0.0632 -0.19165643
#> X6_PCB156             -0.15558488 0.078007253 -1.9944925   0.0461 -0.30847629
#> X7_PCB206             -0.17435195 0.082331927 -2.1176712   0.0342 -0.33571957
#> X10_2.3.4.6.7.8.hxcdf -0.24437058 0.051827916 -4.7150376   <.0001 -0.34595143
#>                         Upper.95CI Contribution proportion
#> X4_a.tocopherol        1.012582158                   0.571
#> X3_g.tocopherol        0.633626279                   0.216
#> X2_retinol             0.383653877                   0.071
#> X8_3.3.4.4.5.pncb      0.238838153                   0.019
#> X5_PCB99               0.122246265                   0.000
#> X1_trans.b.carotene   -0.008754863                   0.000
#> X9_1.2.3.4.7.8.hxcdf   0.005145045                   0.009
#> X6_PCB156             -0.002693477                   0.024
#> X7_PCB206             -0.012984342                   0.030
#> X10_2.3.4.6.7.8.hxcdf -0.142789726                   0.060

# confounder coefficients
model_auto$confounder.coefficient
#>                         Estimate Std. Error   t value Pr(>|t|)  Lower.95CI
#> AGE.c                 0.16891442 0.02479677  6.811954   <.0001  0.12031366
#> SEX.Female           -0.15617521 0.04665633 -3.347353   0.0009 -0.24761994
#> RACE.NH.Black        -0.14604013 0.06145787 -2.376264   0.0177 -0.26649534
#> RACE.MexicanAmerican  0.07995826 0.06015162  1.329279   0.1841 -0.03793675
#> RACE.OtherRace        0.31675619 0.13400282  2.363802   0.0183  0.05411549
#> RACE.Hispanic         0.51328748 0.14005064  3.665013   0.0003  0.23879327
#>                       Upper.95CI
#> AGE.c                 0.21751519
#> SEX.Female           -0.06473049
#> RACE.NH.Black        -0.02558491
#> RACE.MexicanAmerican  0.19785328
#> RACE.OtherRace        0.57939688
#> RACE.Hispanic         0.78778170

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.

Manual control: `plsi.lr.v2()`

If you need a fixed-df spline instead — for example, running many fits in a simulation where mgcv::gam()’s per-call overhead adds up — plsi.lr.v2() provides the same model with a spline.num/spline.degree natural cubic spline in place of the penalized smooth:

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 = 2023)

model_v2$confounder.coefficient
#>                          Estimate Std. Error   t value Pr(>|t|)  Lower.95CI
#> zAGE.c                 0.17020052 0.02484554  6.850344   <.0001  0.12150415
#> zSEX.Female           -0.15666216 0.04673903 -3.351849   0.0008 -0.24826897
#> zRACE.NH.Black        -0.14619532 0.06155057 -2.375206   0.0178 -0.26683222
#> zRACE.MexicanAmerican  0.08055743 0.06025605  1.336919   0.1816 -0.03754226
#> zRACE.OtherRace        0.31728977 0.13414158  2.365335   0.0183  0.05437710
#> zRACE.Hispanic         0.51409085 0.14045182  3.660265   0.0003  0.23881034
#>                        Upper.95CI
#> zAGE.c                 0.21889689
#> zSEX.Female           -0.06505535
#> zRACE.NH.Black        -0.02555841
#> zRACE.MexicanAmerican  0.19865711
#> zRACE.OtherRace        0.58020244
#> zRACE.Hispanic         0.78937136

The remaining sections use model_auto, but every plotting function below works identically with model_v2.

Visualizing results

Single-index coefficients

si.coef.plot(model_auto$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_auto$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_auto, dat, exp_name = "X4_a.tocopherol", type = "linear")

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

e.main.plot(model_auto, 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_auto, 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_auto, 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_auto, 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_auto, dat, type = "linear")

Binary and count outcomes

plsi.logistic.auto() and plsi.log.auto() follow the same interface as plsi.lr.auto() for binary and count outcomes respectively, and every plotting function above accepts a matching type = "logistic" or type = "log" argument. For a binary outcome:

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

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

and similarly for plsi.log.auto() with type = "log" and family = "nb" (negative binomial, the default) or family = "poisson". See ?plsi.logistic.auto and ?plsi.log.auto for full examples, including how si.coefficient and confounder.coefficient report odds ratios and rate ratios respectively for these two outcome types.

An Introduction to EPLSIM Package

Yuyan Wang, Mengling Liu

2026-07-16

Introduction

Preparing confounders

Fitting the model

Recommended: `plsi.lr.auto()`

Manual control: `plsi.lr.v2()`

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

An Introduction to EPLSIM Package

Yuyan Wang, Mengling Liu

2026-07-16

Introduction

Preparing confounders

Fitting the model

Recommended: plsi.lr.auto()

Manual control: plsi.lr.v2()

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

Recommended: `plsi.lr.auto()`

Manual control: `plsi.lr.v2()`