---
title: "An Introduction to EPLSIM Package"
author: 'Yuyan Wang'
date: '`r Sys.Date()`'
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{An Introduction to EPLSIM Package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}{inputenc}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(EPLSIM)
```

## 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`.

```{r}
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)
```

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

```{r}
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()`

```{r}
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

# single-index coefficients: which exposures drive the index, and how much
model_auto$si.coefficient

# confounder coefficients
model_auto$confounder.coefficient
```

`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:

```{r}
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
```

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

## Visualizing results

### Single-index coefficients

```{r, fig.width = 6, fig.height = 5}
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.

```{r, fig.width = 6, fig.height = 5}
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.

```{r, fig.width = 6, fig.height = 5}
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.

```{r, fig.width = 8, fig.height = 5}
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:

```{r, fig.width = 8, fig.height = 5}
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.

```{r, fig.width = 7, fig.height = 6}
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).

```{r, fig.width = 6, fig.height = 5}
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:

```{r, eval = FALSE}
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.