§ INTRODUCTION

§.§ Statistical power

Researchers need to estimate sample size and statistical power as part of good study planning. Power of a test is the probability that it rejects the null hypothesis at a specific significance level under a prespecified alternative hypothesis. Calculating power requires assumptions about relevant quantities (e.g. variance, effect size, sample size). For standard hypothesis tests, including one or two sample t-tests, one-way balanced ANOVA, and chi-squared tests, closed-form equations exist to calculate the power. However, in general, power calculations are more complex. We can estimate power using Monte Carlo (MC) simulations, an approach that is flexible and applicable to a large class of models [Gelman and Hill, 2007]. Each simulation involves first sampling data from a hypothesized data generative model and then fitting an inference model to the data. Out of a large number of simulations, the proportion of times in which a hypothesis test is rejected can be used to estimate power.

§.§ Methods for environmental mixtures

Humans are exposed to a complex mixture of chemicals as part of our daily lives. Statistical models for mixtures aim to achieve several goals, including but not limited to: (1) identifying individual chemicals critical to the health outcome of interest; (2) estimating the health effects of individual (or combinations of) chemicals in the presence of multicollinearity; and (3) examining synergistic interactions between chemicals [Sun et al., 2013]. Several challenges are associated with these goals, and no single approach is best at addressing all of them. For instance, regularized multiple regression with main effects and interactions can examine synergistic effects but tend to be limited to pairwise interactions [Bien et al., 2013, Lim and Hastie, 2015]. Kernel-based regressions can account for high-order interactions and nonlinear effects but tend to not scale well with the number of chemicals or observations [Hamra and Buckley, 2018]. Many models incorporate factor analysis or Bayesian model selection to deal with multicollinearity [Sun et al., 2013, Ferrari and Dunson, 2020, Ferrari and Dunson, 2020]. The choice of the model depends on the goal of analysis and characteristics of the dataset: the dimensionality of the mixtures, the degree of multicollinearity, the sample size, and the hypothesized associations between mixtures and the outcome of interest [Sun et al., 2013]. Thus, it is beneficial to be able to conduct power analysis for different combinations of hypothesized effects and statistical models, while accounting for the degree of multicollinearity in the predictors.

§.§ Our contribution

We present a package containing building blocks to estimate power flexibly when designing studies involving mixtures of chemicals, or more generally correlated predictors. Many R packages implement closed-formed power equations for simple hypothesis tests (package pwr [Champely, 2020]), joint models of longitudinal and survival data (package JMdesign [Cornea et al., 2022]), factorial designs (package powerAnalysis [Fan, 2017]), group sequential design (package PwrGSD [Izmirlian, 2021]), and genomics experiments (package powerGWASinteraction [Kooperberg and Hsu, 2015]). Some packages exist for conducting power analysis using MC simulations (package simr [Green and MacLeod, 2016], skpr [Morgan-Wall and Khoury, 2021], simglm [LeBeau, 2022]). However, none allows users to simulate correlated predictors with a valid joint distribution or gives users as much flexibility in setting up the simulation as our proposed framework. For example, package simr [Green and MacLeod, 2016] only works with generalized linear mixed models. Package skpr [Morgan-Wall and Khoury, 2021] provides options to simulate data from optimal factorial designs and split-plot designs but not correlated mixed-scaled data often seen in health studies. Package simglm [LeBeau, 2022] provides a general framework for power analysis for multi-leveled, mixed-scaled data but assumes independence between predictors. It would also be cumbersome to manually specify all marginal distributions if the number of predictors is large. [Arnold et al., 2011] provided example R code to do power analysis through MC simulations for complex cluster design without considering correlated covariates.

Our package allows researchers to specify or estimate correlation structures for the mixture exposures, generative models for the outcome, and statistical inference models so that even researchers not experienced with R or programming can set up simulations quickly. The building blocks are also general enough to easily interface with other packages, providing access to a wide class of data generative models and statistical models for more advanced users. Because simulation-based power analysis for complex models can be computationally intensive, we allow parallel computation when users have access to multi-core computer processors. In Section <ref>, we describe algorithms for power analysis using MC simulation, simulation of correlated mixed-scaled predictors, automatic scaling of effect size, and state-of-the-art statistical models for chemical mixtures. In Section <ref>, we provide a tutorial on the package's functionalities with example code. In Section <ref>, we provide a tutorial on end-to-end power analyses with synthetic data and NHANES data using our package. In Section <ref>, we provide an example on how to estimate the power of a more involved study and use power analyses as a tool for inference models comparison.

§ METHODOLOGY

§.§ Power analysis using Monte Carlo simulation

We can use Monte Carlo simulation to estimate power given a data generative model and an inference model. Let $s$ be the number of simulations, $n$ be the sample size of one simulation, $y_i, i=1,...,n$, be the outcome of the $ith$ subject, and $\theta = (\theta_X, \theta_Y)$ be the parameters (i.e., intercepts, main effects, interactions, variance components, correlations among predictors, ...) in a model that realistically represents the data generating process (e.g., based on the study design [Arnold et al., 2011]) and a hypothesis about the effect. This data generative process will be constructed by a generative model for the predictors and a generative model for the outcome, discussed in the next section. The parameters $\theta$ may be chosen based on prior knowledge of the effect size or estimates from previous studies [Arnold et al., 2011, Green and MacLeod, 2016].

Power analysis using Monte Carlo simulation

[1]
Set $\theta$ based on prior knowledge and/or estimates from previous studies.
Simulate $n$ independent samples of the predictors $\boldsymbol{x}_i, i=1,...,n$ given $\theta_X$.
Simulate $y_i$ given parameters $\theta_Y$ and predictors $\boldsymbol{x}_i$.
Perform analysis with an inference model and record if the effects meet some “significance” criterion (e.g. quantile-based or highest posterior density credible interval coverage).
Calculate the power as the proportion of simulations meeting the defined criterion.

We may repeat Algorithm <ref> with different values for $\theta$ and $n$ to create a power curve. Power curves are useful to study the trade-offs between power, sample size, and effect size.

§.§ Simulate correlated predictors

We provide three methods to construct a generative model for correlated predictors. One approach to simulate samples from a realistic joint distribution of the predictors is to resample them from existing data. Another option is to estimate the joint distribution of the existing data and simulate from it. In this package, we use a semiparametric Bayesian Gaussian copula (SBGC) [Hoff, 2007] to estimate the joint distribution of mixed binary, ordinal, and continuous data. This method is based on Sklar's Theorem that any multivariate probability distribution can be written in terms of its univariate marginal probability distributions and the copula. For example, let $x_1, ..., x_p$ be $p$ random variables with marginal univariate CDF's $F_1, ..., F_p$, not necessarily continuous, and an unknown joint distribution we want to estimate. The transformed variables $u_j=F_j(x_j)$'s all have uniform marginals. A copula models the joint distribution of $u_1,...,u_p$, which also describes the dependence amongst the $x_j$'s separately from $F_1,...,F_p$. We use a Gaussian copula which has the following form:

\begin{align}
(z_1,...,z_p)^T|\boldsymbol{C} &\sim N_p(\boldsymbol{0}, \boldsymbol{C}) \\
x_j &= F_j^{-1}[\Phi(z_j)], j=1,...,p
\label{eqn:copula}
\end{align}

where $z_j$'s are latent variables, $\boldsymbol{C}$ is a correlation matrix of the latent variables, $F_j^{-1}$ is the pseudo-inverse of $F_j$, and $\Phi$ the CDF of a standard normal. Since the closed form of the $F_j$'s are often unknown in practice, we can estimate them using the empirical distributions from $n$ observations $\hat{F}_j(t) = \frac{1}{n}\sum_{i=1}^n \boldsymbol{1}_{x_{ij} \leq t}$. Hoff proposed using the extended rank likelihood to estimate $\boldsymbol{C}$ without depending on estimates of $F_j$'s [Hoff, 2007], which works well for mixed-scaled data, and implemented the Bayesian inference algorithm in the R package sbgcop [Hoff, 2018]. The uncertainty in the estimation can be propagated through the whole power calculation straightforwardly by sampling from the copula's parameters' posteriors before simulating new data. Additionally, no specification of prior distributions or parametric form of the univariate marginal distributions is necessary, ensuring ease of use. We summarize how to simulate data from a generative model defined by a Gaussian copula in Algorithm <ref>.

Simulate $p$ correlated predictors from a Gaussian copula

[1]
Given the univariate marginal CDF $F_j$'s defined by users or estimated with the empirical CDF.
Sample a correlation matrix $\boldsymbol{C}$ from the SBGC posterior (functionsbgcop::sbgcop.mcmc()) or from Algorithm <ref>.
Sample the latent variables $\boldsymbol{z}_i \sim N_p(0, \boldsymbol{C})$.
Get the predictors $\boldsymbol{x}_i$ through the transformation $x_j = F_j^{-1}[\Phi(z_j)], j=1,...,p$.
Return an $n \times p$ matrix $X$.

Alternatively, if we have no existing data to estimate correlation matrix $\boldsymbol{C}$ of the latent variables and the marginal CDF's $F_j$ of the observed variables, we can specify them manually. For instance, we may assume a continuous predictor to be normally distributed with a given mean and standard deviation, or an ordinal predictor to be a multinomial distribution with given probabilities. As for the correlation matrix $\boldsymbol{C}$, given reasonable guesses on pairwise correlations between predictors, we can sample positive definite matrices with desirable structures using Algorithm <ref> from [Lewandowski et al., 2009]. This can add additional uncertainty into the samples of hypothetical data. Let $\rho_{ij;kL}$ be the partial correlation of the $i^{th}$ and $j^{th}$ variables while holding the $k^{th}$ variable and variables in set $L$ constant, where $L$ contains indices distinct from $\{i,j,k\}$. We can calculate partial correlation $\rho_{ij;L}$ using the following formula from [Lewandowski et al., 2009]:

\begin{equation}
\rho_{ij;L} = \rho_{ij;kL}
\sqrt{(1-\rho_{ik;L}^2)(1-\rho_{jk;L}^2)} + \rho_{ik;L}\rho_{jk;L}
\label{eqn:partial-r-formula}
\end{equation}

The partial correlation $\rho_{ij}$ where $L$ is the empty set is equivalent to the correlation $r_{ij}$ between the $i^{th}$ and $j^{th}$ variables. Equation <ref> shows that a $p \times p$ random correlation matrix can be parametrized in terms of pairwise correlations $\rho_{12},...,\rho_{1p}$ and partial correlations $\rho_{23;1}$,...,$\rho_{2p;1}$,$\rho_{34;12}$,..., $\rho_{3p;12}$,...,$\rho_{(p-1)p;1...p-2}$. These choices of partial correlations correspond to a C-vine [Bedford and Cooke, 2002]. They can independently take values on the interval $(-1,1)$ following any density with the right support [Lewandowski et al., 2009]. We sample them from independent $Be(\alpha, \beta)$ distributions (scaled to be on $(-1, 1)$) with appropriate hyperparameters to induce a correlation matrix with desirable structures. Applying Equation <ref> iteratively to smaller index sets $L$, we can calculate every element of the correlation matrix from these partial correlations (Algorithm <ref>). For example, $\rho_{23} = \rho_{23;1} \sqrt{(1-\rho_{12}^2)(1-\rho_{13}^2)} + \rho_{12}\rho_{13}$ where $k=1, L = \emptyset$, and $\rho_{23;1}, \rho_{12}, \rho_{13}$ are Beta random variables. More details can be found in [Joe, 2006, Lewandowski et al., 2009]. When choosing $\{\alpha, \beta\}$, take into consideration that the mean of a Beta random variable is $\frac{\alpha}{\alpha + \beta}$, and the larger $\alpha$ or $\beta$ is, the smaller the variability of the samples. To automatically set reasonable parameters $\{\alpha_{ij}, \beta_{ij}\}$ for $\rho_{ij;1...i-1}$, we use the following heuristics given guesses of pairwise correlations, which are easier to specify than partial correlations. First, we calculate a partial correlation $\hat{\rho}_{ij;1...i-1}$ using the conditional multivariate normal formula [Eaton, 2007] from the rough guesses of pairwise correlations. If the guessed correlation matrix is singular, the generalized inverse can be used in the conditional normal formula. We then set $\alpha_{ij} = m \hat{\rho}_{ij;1...i-1}$ and $\beta_{ij} = m - \alpha_{ij}$ where $m$ controls the magnitude of the parameters. We find the heuristics to work well when the guessed pairwise correlations are not too extreme (e.g. 0.99). Simply put, we can turn rough guesses of pairwise correlations into a valid sampling distribution for the predictors.

Sample random correlation matrix with C-vine

[1]
Set hyperparmeters $\alpha, \beta$ appropriately for each partial correlation
Sample partial correlation $\rho_{ij;1,...,i-1} \sim Be(\alpha_{ij}, \beta_{ij})$ on $(-1,1)$.
Calculate the correlation $r_{ij} = r_{ji} = \rho_{i,j}$ recursively using Equation <ref>.
Return the correlation matrix $\textbf{R}$ where $\textbf{R}_{ij} = r_{ij}$.

§.§ A general notion of effect size

It is not straightforward to set an effect size for nonlinear or complex mean functions. For the model $Y = f(X) + \epsilon$ with independent, additive Gaussian noise $\epsilon \sim N(0, \sigma_\epsilon^2)$, we propose using the signal-to-noise ratio as a general proxy to the joint effect size for any mean function. It is defined as:

\begin{equation}
\text{SNR} = \frac{Var(f(X))}{\sigma_\epsilon^2}
\label{eqn:snr-continuous}
\end{equation}

Fixing either the mean function variance or the noise variance, it is straightforward to scale the other so that a SNR $\gamma$ is achieved. We can estimate the variance of any mean function with $\hat{\sigma}_f^2 = \frac{1}{m-1}\sum_{i=1}^m (f(x_i) - \frac{1}{m}\sum f(x_i))^2$, where $x_1,...,x_m$ are i.i.d. realizations of predictors $X$ for a large number $m$. The specific calculations are summarized in Algorithm <ref> and <ref>. We may extend this measure of effect size to binary or count outcomes using the SNR estimator for generalized linear models (GLM) proposed by [Czanner et al., 2008]. The estimator is defined as:

\begin{equation}
\widehat{\text{SNR}} = \frac{Dev(\textbf{w}, \bar{w}1) - Dev(\textbf{w}, \textbf{h}) + 1 - p}{Dev(\textbf{w}, \textbf{h}) + p}
\label{eqn:snr-discrete}
\end{equation}

where $\textbf{w}, \textbf{h}$ are vectors of realizations of outcome $W$ and mean model $h(X)$ respectively, and $\bar{w} = \frac{1}{n}\sum_{i=1}^n w_i$. The residual deviance $Dev( \textbf{w}, \textbf{h})$ is a measure of discrepancy between the observed and predicted data, with the following forms for binary and count outcome [McCullagh and Nelder, 1983]:

\begin{align}
%\text{Gaussian} &: \sum_{i=1}^n (y_i - \mu_i)^2 \\
\text{Bernoulli} &: 2\sum_{i=1}^n w_i log(\frac{w_i}{h_i}) + (1-w_i)log(\frac{1-w_i}{1-h_i}) \\
\text{Poisson} &: 2 \sum_{i=1}^n y_i log(\frac{w_i}{h_i}) - (w_i - h_i)
\label{eqn:deviance-discrete}
\end{align}

For binary or count outcomes, we can adjust parameters in the mean model $h$ until we get a desirable estimated $\widehat{\text{SNR}}$.

Scale the noise variance

[1]
Given a model of the outcome $Y = f(X) + \epsilon$ with independent, additive Gaussian noise $\epsilon \sim N(0, \sigma_\epsilon^2)$, the predictors $X \sim N(0,\Sigma_x)$, and the desired SNR $\gamma$.
Sample $x_i \sim N(0,\Sigma_x), i=1,...,m$ for large $m$.
Estimate $\hat{\sigma}_f^2 = \frac{1}{m-1}\sum_{i=1}^m (f(x_i) - \frac{1}{m}\sum f(x_i))^2$.
Set the noise variance to $\sigma^2_\epsilon = \frac{\hat{\sigma}_f^2}{\gamma}$.

Scale the mean function

[1]
Given a model of the outcome $Y = f(X) + \epsilon$ with independent, additive Gaussian noise $\epsilon \sim N(0, \sigma_\epsilon^2)$, the predictors $X \sim N(0,\Sigma_x)$, a fixed noise variance $\sigma_\epsilon^2$, and the desired SNR $\gamma$.
Sample $x_i \sim N(0,\Sigma_x), i=1,...,m$ for large $m$.
Estimate $\hat{\sigma}_f^2 = \frac{1}{m-1}\sum_{i=1}^m (f(x_i) - \frac{1}{m}\sum f(x_i))^2$.
Calculate a scaling factor $ s = \gamma \frac{ \sigma_\epsilon^2}{\hat{\sigma}_f^2}$.
Set $f_s = \sqrt{s} f$.

§.§ Statistical methods for environmental mixtures

Our package interfaces with several existing and newly developed analysis strategies for assessing associations between correlated exposures and health outcomes. These strategies can be used as inference models for Algorithm <ref>. We provide an overview of methods and their unique strengths. We do not go into the mathematical description, and instead give intuitions that help researchers select the correct model given characteristics of the data. For a more detailed review of recent methodological developments, see [Joubert et al., 2022].

§.§.§ Bayesian Model Averaging

Model selection is a common strategy to deal with multicollinearity or large dimensionality in the predictors. Instead of selecting for one best model that might be sensitive to the given observed data set, Bayesian Model Averaging (BMA) provides a posterior over the space of all possible models [Hoeting et al., 1999]. Thus, BMA allows for direct model selection and robust effect estimation by averaging the posteriors of the parameters under all models weighted by their posterior model probabilities. Unlike other model selection methods (i.e. forward/backward selection), BMA quantifies uncertainty more fully, including model uncertainty. Interactions and nonlinear effects can be included by transformation of the design matrix. Additionally, BMA can be used with many models, such as linear regression, generalized linear models, and survival models. Our package integrates the implementation of BMA in package BMA [Raftery et al., 2021]. As the dimension of the model space grows exponentially with the number of predictors, for high-dimensional data sets, implementations of BMA often approximate the model posterior by doing preliminary stepwise model selection (i.e. in package BMA) or by sampling the model space. In the presence of high correlations among many predictors, BMA is prone to produce low marginal inclusion probabilities for the correlated predictors, making it difficult to identify critical predictors of the outcome.

§.§.§ Bayesian Kernel Machine Regression

Bayesian kernel machine regression (BKMR) can flexibly model nonlinear, high-order interaction effects of multiple predictors jointly using a kernel function [Bobb et al., 2015]. BKMR incorporates component-wise variable selection, producing a posterior inclusion probability for each individual variable, allowing for identification of critical components. Similar to BMA, BKMR is prone to produce low posterior inclusion probabilities for a large number of highly correlated variables even when many or all of them are important. To alleviate this issue, BKMR also implements hierarchical variable selection, which selects on groups of highly correlated variables. Researchers will need to specify the groups based on domain knowledge of the data set. BKMR can be prone to overfitting [Ferrari and Dunson, 2020]. Our package integrates the implementation of BKMR from the package bkmr available on CRAN [Bobb et al., 2018].

§.§.§ MixSelect

MixSelect (MS) decomposes the fitted mean function of the outcome into linear main effects, pairwise interactions, and a nonlinear deviation [Ferrari and Dunson, 2020]. MS also include variable selection on the linear main effects and hierarchical variable selection on the interactions. This means that a pairwise interaction is included only if the main effects of one or both variables are selected. The nonlinear deviation models the residuals of the linear terms, using a kernel function similar to BKMR. Unlike BKMR, MS uses Bayesian model selection to decide whether to include the nonlinear deviation in the model. Thus, it automatically simplifies to a linear model with variable selection if only linear effects are present in the data. While BKMR relies solely on visualizations for interpretation of marginal effects, MS produces posteriors for the linear terms from which the mean, credible intervals, and other statistics can be calculated. However, MS does not offer group-wise variable selection like BKMR. Like BKMR, MS suffers with data of large dimensionality, sample size, and extensive multicollinearity. We use the implementation provided by the authors at <https://github.com/fedfer/MixSelect>.

§.§.§ Factor Analysis for Interactions

Factor analysis is a common choice to deal with multicollinearity in the predictors. It can infer group structures in the predictors, i.e. latent factors, while being informed by their relationships to the outcome, effectively reducing the dimension of the data. [Ferrari and Dunson, 2020] proposed a factor analysis for interaction (FIN) model that accounts for pairwise interactions between the factors and transforms the parameters to produce linear main effects, quadratic main effects and pairwise interactions in the original predictors. Thus, it examines effects of combinations of predictors as well as individual predictor at a reduced model complexity. Reasonable assumptions about group structures in the predictors must hold. FIN does not produce exact variable selection, though it induces shrinkage on the effect estimates. FIN might not be effective at identifying individual critical predictors, especially when other correlated predictors are not important. We use the implementation provided in the package infinitefactor [Poworoznek, 2020] available on CRAN.

§.§.§ Bayesian Weighted Sums

Bayesian Weighted Sums (BWS) aims at estimating a combined effect and individual variable importance by summing multiple predictors into a single exposure variable using quantitative weights [Hamra et al., 2021]. Thus, BWS reduces the dimension of the problem and produces a parsimoneous model. BWS may be thought of as a factor model that assumes only one latent factor affects the outcome. The weights can be interpreted as proportions of the combined effect explained by each component. Unlike Quantile G-Computation (discussed next), BWS provides uncertainty quantification on the weights. BWS is currently limited to modeling linear additive effects and constrains all components of the mixtures to have effects of the same sign. We use the implementation of BWS in the package bws [Nguyen, 2022] available on CRAN.

§.§.§ Quantile G-Computation

Similar to BWS, Quantile G-Computation (QGC) aims to estimate the combined effects of multiple predictors. QGC first transforms the predictors into categorical variables whose categories correspond to quantiles of the original predictors. QGC then fits a generalized linear model to the categorical variables, and the sum of all regression coefficients produces the g-computation estimator. Since the underlying model is a generalized linear model, if these three assumptions hold: 1) the relationship is in fact linear, 2) quantization of predictors is appropriate, and 3) there are no unmeasured confounders, then the g-computation estimator gives a causal dose-response parameter [Keil et al., 2020]. Unlike BWS, QGC allows predictors to have effects of different signs and allows the inclusion of nonlinear and interaction terms. Identification of critical components may be difficult because the weights lack confidence intervals or p-values. We use the R package qgcomp available on CRAN [Keil, 2021].

§ TUTORIAL: FUNCTIONALITIES

We will walk through all functionality of the package. There are three data models that are building blocks to running a power analysis with this package: a generative model for the predictors, a generative model for the outcome given the predictors, and an inference model.

§.§ Specifying a predictors model

If sufficient existing data are available, simulated predictors may be created by resampling. For example, if we were planning a study of metabolic disease and phthalates, we could sample phthalate levels from NHANES data:

R> xmod <- MixtureModel(method = "resampling", data = nhanes,
+    resamp_prob = weights)

A vector of sampling probabilities of the observations, which must sum to 1, may be given to argument resamp_prob. If limited data exist, predictors may be sampled from a semi-parametric Bayesian copula parametrized by the empirical univariate marginals and a correlation matrix [Hoff, 2007]. The marginals and correlation matrix of the latent variables may be estimated directly as discussed in Section 2.2. This method can handle mixed-scaled data (i.e., continuous, binary, ordinal). The generated predictors are on the same scale as the original data. Note that categorical and ordinal data must be formatted as numeric (e.g. using one-hot-encoding). Users may pass a list of named arguments, including the number of MCMC samples (nsamp), how much to thin a chain (odens), and whether to print MCMC progress (verb) to the argument sbg_args. See functionsbgcop.mcmc() from package sbgcop for the full list of arguments.

R> sbg_args <- list(nsamp = 2000, odens = 2, verb = TRUE)
R> xmod <- MixtureModel(method = "estimation", data = existing_data,
+    sbg_args = sbg_args)

Alternatively, if no existing data are available to estimate the correlation matrix of the latent variables and univariate marginals, users can specify them manually using Algorithm <ref>. The argument G takes a $p \times p$ matrix of rough guesses of pairwise correlations between all variables. Larger values for the argument m corresponds to smaller deviation from G in the random samples. In the code below, we generate random matrices with high pairwise correlations around 0.9, as well as a random matrix with group structures:

R> set.seed(1)
R> p <- 40
R> variable_names <- paste0("x", 1:p)
R> G <- diag(1, p)
R> G[upper.tri(G)] <- 0.9
R> G[lower.tri(G)] <- 0.9
R> cvine_marginals <- lapply(1:p, function(i) "qnorm(mean = 0, sd = 1)")
R> colnames(G) <- rownames(G) <- names(cvine_marginals) <- variable_names
R> small_devi <- MixtureModel(method = "cvine", G = G, m = 100,
+    cvine_marginals = cvine_marginals)
R> large_devi <- MixtureModel(method = "cvine", G = G, m = 10,
+    cvine_marginals = cvine_marginals)
R> G_group <- matrix(c(rep(c(rep(0.8, p/2), rep(0.3, p/2)), p/2),
+    rep(c(rep(0.3, p/2), rep(0.7, p/2)), p/2)), p, p)
R> diag(G_group) <- 1
R> group <- MixtureModel(method = "cvine", G = G_group, m = 100,
+    cvine_marginals = cvine_marginals)

Examples of random correlation matrices generated using Algorithm <ref>.

The sampled correlations concentrate more around 0.9 when m=100 than when m=10 (Figure <ref> left and Figure <ref> middle). The right-most subplot of Figure <ref> shows an example matrix created by an input G with group structures. The argument cvine_marginals takes a named list of strings describing the quantile function and named parameters of the marginal distribution of each variable. Valid strings include standard distributions from the R package stats. See our documentation for a complete list of supported distributions. For example, to create independent binary, count, ordinal and categorical variables, we use the Bernoulli, Poisson, and multinomial distributions as follows:

R> xmod <- MixtureModel(method = "cvine", G = diag(1, 4), m = 50,
+    cvine_marginals = list(binary = "qbinom(size = 1, prob = 0.2)",
+                           count = "qpois(lambda = 1)",
+                           ordinal = "qmultinom(probs = c(0.6, 0.2, 0.2))",
+                           categorical = "qmultinom(probs = c(0.6, 0.2, 0.2))"),
+    cvine_dtypes = list(categorical = "factor"))

All variables are assumed numeric unless specified in cvine_dtypes otherwise. In our example above, we want the categorical variable to be factor.

§.§ Specifying an outcome model

Users need to specify a hypothesized model for the outcome given the predictors. This can be thought of as the “true” relationship between the predictors and outcome. This model may be informed from prior knowledge or pilot studies. When little information is available, investigators may wish to test a variety of outcome models to get a sense of robustness to assumptions about predictor-outcome relationships. The model for the mean is given as a string of R code. For example, if the continuous outcomes are modeled as random draws from a Gaussian $N(2(x_1 + x_2), 1)$, we can use the following code:

R> ymod <- OutcomeModel(f = "2*(x1 + x2)", family = "gaussian", sigma = 1)

The function can generate binary and count data using an appropriate link function (i.e., logit) when family is set to "binomial" or "poisson" respectively. A fitted model of the outcome to existing data may also be used but needs to be wrapped inside a function that takes a predictor matrix  and returns a vector of mean values. Below is an example of wrapping an lm model to be used with our package:

R> pilot_mod <- lm(y   x1, data = pilot_study)
R> pilot_f <- function(x)
+    predict(pilot_mod, as.data.frame(x))$fit
+  }
R> pilot_sigma <- sigma(pilot_f)
R> ymod <- OutcomeModel(f = pilot_f, family = "gaussian", sigma = pilot_sigma)
\end{Code}

\subsection{Specifying an inference model}
While the outcome and the inference models are ideally the same, in practice our model will at best approximate the ``true'' predictor-response relationships. In the code below, we show how to use BKMR as the inference model. The \verb|model| argument also accepts one of the following values: \verb|"glm"|, \verb|"mixselect"|, \verb|"qgcomp"|, \verb|"bws"|, \verb|"fin"|, and \verb|"bma"|. Advanced users may also define custom inference models (see first example in Section \ref{sec:synthetic}). Additional arguments specific to the inference model (e.g. \code{iter}, \code{varsel}) can be passed as done below:

\begin{Code}
R> imod <- InferenceModel(model = "bkmr", iter = 10000, varsel = TRUE)
\end{Code}

The built-in inference functions return values, such as p-value or posterior inclusion probability, as ``significance'' criteria. Table~\ref{tab:sig-criteria} list the criteria available for each built-in inference model.

\begin{table}[ht]
\centering
\begin{tabular}{| p{0.1\linewidth} | p{0.55\linewidth}| p{0.25\linewidth}|}
\hline
Inference method & Individual effect & Value to \verb|crit| \\
\hline
BMA & posterior inclusion probability (PIP) of a main effect or an interaction term $\geq \alpha$. & \verb|"pip"|\\
\hline
GLM & the smallest p-value of the regression coefficients (i.e., main effect, interactions) involving a predictor $\leq \alpha$. & \verb|"pval"|, \verb|"main_pval"| (main effect only), \verb|"int_pval"| (interactions only)\\
\hline
BKMR & PIP of a non-zero length-scale parameter $\geq \alpha$. & \verb|"pip"| (component-wise), \verb|"group_pip"| (group-wise)\\
\hline
MS & PIP of a main effect or a length-scale parameter $\geq \alpha$. & \verb|"pip"|, \verb|"linear_pip"| (main effect alone), \verb|"gp_pip"| (lenghth-scale alone)\\
\hline
FIN & $(1-\alpha)$ \% CI of a main effect or an interaction term doesn't include zero. Equivalent to $min(Pr(\beta >0|.), Pr(\beta < 0|.)) \leq \alpha$ where $Pr(\beta < 0|.)$ is the posterior probability of an effect being less than 0. & \verb|"beta"|, \verb|"linear_beta"| (main effect alone)\\
\hline
BWS & $(1-\alpha)$ \% CI of the joint effect doesn't include zero. Equivalent to $min(Pr(\beta >0|.), Pr(\beta < 0|.)) \leq \alpha$ where $Pr(\beta < 0|.)$ is the posterior probability of the joint effect being less than 0. & \verb|"beta"|\\
\hline
QGC & p-value of the joint effect $\leq \alpha$. & \verb|"pval"|\\
\hline
\end{tabular}
\caption{Definitions of effect criteria for built-in inference models given a threshold $\alpha$.}
\label{tab:sig-criteria}
\end{table}

\subsection{Running a power analysis}
Function\code{sim_power()} executes steps 2 to 6 of Algorithm \ref{alg:power_mc} given the predictors generative model, outcome generative model(s), inference model, sample size(s), and the number of MC iterations. In the code below, we simulate 2000 data points specified by the predictors and outcome model, fit an inference model and repeat the whole process 100 times:

\begin{Code}
R> res <- sim_power(xmod, ymod, imod, s = 100, n = 2000,
+    snr_iter = 100000, cores = 2, errorhandling = "stop",
+    cluster_export = c("pilot_mod"))
\end{Code}

Each simulation is independent and thus may be parallelized for faster computation. When multiple cores are available on the computer, increase \verb|cores| to parallelize the simulation. Internally, we use the \pkg{doSNOW} [Corporation and Weston, 2022] and \pkg{foreach} [Microsoft} and Weston, 2020] packages for this feature. While using parallelism, if an error occurs in any iteration, \verb|errorhandling| specifies whether to remove that iteration (\verb|"remove"|), return the error message verbatim in the output (\verb|"pass"|), or terminate the loop (\verb|"stop"|). The \verb|"pass"| and \verb|"stop"| options are useful for debugging. Larger values for \verb|snr_iter| results in more precise estimates of the SNR. Finally, any global variable not explicitly passed in as an argument (e.g. a custom outcome model) needs to be exported to the parallel cluster using \code{cluster_export}.

\subsection{Creating a power curve}

It is often desirable to study the trade-offs between power, Type I error rate, sample size, and effect size. Function\code{sim_curve()} generates a power curve by running the power analysis described previously on a combination of sample sizes and effect sizes. Different effect sizes may be baked into the outcome model and measured by the signal to noise ratio. Function\code{sim_curve()} takes a list of outcome model and a vector of sample sizes, but only one predictors model and one inference model at a time. Below we conduct power analyses for 9 combinations of sample sizes 50, 100, and 200 and 3 outcome generative models:

\begin{Code}
R> res_curve <- sim_curve(xmod = xmod, imod = imod,
+    ymod = list(ymod1, ymod2, ymod3), s = 1000, n = c(50, 100, 200))
\end{Code}

\subsection{Estimating the signal-to-noise ratio}

As we discuss in Section \ref{sec:method}, the SNR is a good proxy for a generalization of an effect size to non-linear outcome model. To estimate the SNR of a data generating process and its bootstrap standard error (s.e.) , we can pass a predictors model and an outcome model to the following function:

\begin{Code}
R> estimate_snr(xmod, ymod, m = 100000, R = 1000)
\end{Code}

Argument \verb|m| define the number of samples drawn from the predictors and outcome models to estimate SNR and \verb|R| the number of bootstrap replicates. For illustration, we will estimate the SNR of the following data-generating process using different values for \code{m}:

\begin{Schunk}
\begin{Sinput}
R> set.seed(1)
R> xmod <- MixtureModel(method = "resampling",
+    data = data.frame(x1 = rnorm(200000, mean = 0, sd = 1),
+                      x2 = rnorm(200000, mean = 0, sd = 1)))
R> ymod <- OutcomeModel(f = "0.3 * x1 + 0.3 * x2", sigma = 1, family = "gaussian")
R> for (m in c(500, 5000, 50000, 100000)){
+    estimate_snr(ymod, xmod, m = m, R = 1000)
+  }
\end{Sinput}
\end{Schunk}

Since the predictors are independent standard normal distributions, and the noise variance is 1, we can calculate the true SNR as $[0.3^2(1) + 0.3^2(1)]/1 = 0.18$. Figure~\ref{fig:snr-compare} shows the estimated SNR and 1000-bootstrap s.e. for \code{m} $∈{500, 5000, 50000, 10000, 200000}$. A larger \verb|m| results in a more precise estimate. When the mixture model is defined based on resampling, it may not be possible to choose a large \verb|m| without duplicating observations and underestimating the signal.

\begin{figure}[!ht]
\centering
\includegraphics{article-vizSNR}
\caption{The estimated SNR for the linear model example is unbiased but the standard error might be large with a small sample of simulated data. The red horizontal line is the ground truth SNR.}
\label{fig:snr-compare}
\end{figure}

If the outcome generative model has Gaussian white noise, we can automatically scale the mean function (use\code{scale_f()}) or the noise variance (use\code{scale_sigma()}) to get a desired SNR as discussed in Algorithms \ref{alg:scale_noise} and \ref{alg:scale_mean}. The code below modifies a given outcome model and returns a new one with the specified SNR:

\begin{Code}
R> desired_snr <- 0.5
R> new_ymod <- scale_f(desired_snr, cur_ymod, xmod)
R> new_ymod <- scale_sigma(desired_snr, cur_ymod, xmod)
\end{Code}

\subsection{Getting the summary}

The\code{summary()} function produces summaries statistics and plots for the power analysis. It takes a `\code{Sim}' object or a `\code{SimCurve}' object returned by\code{sim_power()} or\code{sim_curve()}. The \verb|crit| argument specifies the significance criterion; the \verb|thres| argument specifies the significance threshold; and \verb|how| describes how the criterion should be compared to the threshold (\verb|"lesser"| or \verb|"greater"|) . For example, for BMA, if PIP $≥0.5$ is considered ``significant'', then the criterion is PIP and the threshold is 0.5. For Bayesian parametric models, credible interval can also be used as a criterion, and p-values for frequentist models (e.g. QGC). Table~\ref{tab:sig-criteria} shows criteria available for each built-in inference model and corresponding values to give to \verb|crit|. The code below shows how to extract summaries for a BMA with PIP greater than 0.5, a GLM with p-values less than 0.01, and a BWS model with zero-coverage of the 90\% or larger credible interval as ``significance'' thresholds:

\begin{Code}
R> summary(bma_res, crit = "pip", thres = 0.5, how = "greater")
R> summary(glm_res, crit = "pval", thres = 0.01, how = "lesser")
R> summary(bws_res, crit = "beta", thres = 0.1, how = "lesser")
\end{Code}

Different thresholds will affect the power of the test. We may plot how the power changes for different thresholds:

\begin{Code}
R> plot_summary(bma_res, crit = "pip", thres = c(0.5, 0.6, 0.7, 0.8), how = "greater")
\end{Code}

\section[Synthetic Examples]{Tutorial: Synthetic Examples} \label{sec:synthetic}

In this section, we provide three examples of doing end-to-end power analysis with the functionality presented in the last section. The first example shows how to use C-vine for simulating predictors, to define a simple custom inference function, and to generate a power curve. The second example shows how to estimate realistic dependence structure from existing data. The third example shows a more complex example of a linear regression and resampling of existing data.

\subsection{Power curve of F-test for linear model}

In the first example, we generate predictors using the C-vine method and calculate power in an F-test for the general linear model. A closed-form power formula exists for the F-test and can be used to validate power estimated with MC simulations. The code below defines a generative model for two continuous, one binary, and one ordinal predictor with 3 levels, all independent of one another:

\begin{Schunk}
\begin{Sinput}
R> xmod <- MixtureModel(method = "cvine", G = diag(1, 4), m = 1000,
+    cvine_marginals = list(x1 = "qnorm(mean=0, sd=1)",
+                           x2 = "qnorm(mean=0, sd=1)",
+                           x3 = "qbinom(size=1, prob=0.7)",
+                           x4 = "qmultinom(probs=c(0.5, 0.3, 0.2))"))
\end{Sinput}
\end{Schunk}

\begin{figure}[!ht]
\centering
\includegraphics{article-exOneXFig}
\caption{Output of plotting the C-vine mixture model with \code{mplot(xmod)}. Note that the C-vine correlation matrix is for latent variables in Algorithm~\ref{alg:copula}.}
\label{fig:cvine-plot}
\end{figure}

Figure~\ref{fig:cvine-plot} shows the marginals and correlation matrix of the mixture model. Note that the C-vine correlation matrix is for latent variables in Algorithm~\ref{alg:copula}. We hypothesize the continuous outcome $y$ to be a linear function of $x_1$, $x_2$ and Gaussian white noise with variance 1. Mathematically, $y_i ∼N(β_1 x_i1 + β_2 x_i2, 1)$ where $i$ indexes the observation, and $β_1, β_2$ are the regression coefficients. The code below defines three hypotheses on the outcome model with decreasing SNRs:

\begin{Schunk}
\begin{Sinput}
R> ymod_list <- list(
+    OutcomeModel(f = "0.3 * x1 + 0.3 * x2", sigma = 1, family = "gaussian"),
+    OutcomeModel(f = "0.2 * x1 + 0.2 * x2", sigma = 1, family = "gaussian"),
+    OutcomeModel(f = "0.1 * x1 + 0.1 * x2", sigma = 1, family = "gaussian")
+  )
\end{Sinput}
\end{Schunk}

Since $x_1 ⊥x_2$ and $x_1, x_2, y$ are Gaussian, the SNR can be calculated exactly analytically, e.g., $Var(0.3x_1 + 0.3x_2)/σ^2 = 0.3^2(1) + 0.3^2(1)/1 = 0.18$ for the first, 0.08 for the second, and 0.02 for the third outcome model. Nevertheless, note that the empirical estimates of SNRs might deviate slightly from the truth, as discussed in the last section. Note also that in the ``true'' generative model, $x_3$ and $x_4$ are unrelated to $y$, but this is unknown to the analyst. Next, we define an inference model including all four predictors. Since the F-test is not a built-in inference model, we will need to define it in a way that is compatible with the \class{InferenceModel} object. Specifically, it needs  to take two arguments \code{x}, \code{y} for a predictors matrix and a vector of outcomes respectively, and to return a named list of significance criteria and values to be compared to a threshold $α$ later. For the F-test, one significance criterion is the p-value from the overall F-test:

\begin{Schunk}
\begin{Sinput}
R> ftest <- function(y, x) {
+    dat <- as.data.frame(cbind(y, x))
+    lm_mod <- lm(y ~ ., data = dat)
+    fstat <- summary.lm(lm_mod)$fstat
+    fpval <- pf(fstat[1], fstat[2], fstat[3], lower.tail = F)
+    names(fpval) <- "F-test"
+    return(list(pval = fpval))
+  }
R> imod <- InferenceModel(model = ftest, name = "F-test")
\end{Sinput}
\end{Schunk}

We now have all the components to run power analysis for sample sizes 50, 100, 150, and 200. Since we consider 4 sample sizes and 3 hypothesized outcome models, we will use function\code{sim_curve()} to estimate 12 power calculations using 500 MC simulations each:

\begin{Schunk}
\begin{Sinput}
R> set.seed(1)
R> curve <- sim_curve(xmod, ymod_list, imod, s = 1000,
+    n = c(50, 100, 150, 200), cores = 1)
\end{Sinput}
\end{Schunk}

We can extract the estimated power of the ovarll F-test based on the significance criterion ``p-value $\leq 0.05$'' as a data frame with\code{summary()}:

\begin{Schunk}
\begin{Sinput}
R> curve_df <- summary(curve, crit = "pval", thres = 0.05, how = "lesser")
\end{Sinput}
\begin{Soutput}
*** POWER CURVE ANALYSIS SUMMARY ***
Number of Monte Carlo simulations: 1000
Number of observations in each simulation: 50 100 150 200
Data generating process estimated SNR (for each outcome model): 0.18 0.08 0.02
Inference model: F-test
Significance criterion: pval
Significance threshold: 0.05

|       | power|   n|  snr|
|F-test | 0.562|  50| 0.18|
|F-test | 0.305|  50| 0.08|
|F-test | 0.095|  50| 0.02|
|F-test | 0.914| 100| 0.18|
|F-test | 0.538| 100| 0.08|
|F-test | 0.152| 100| 0.02|
|F-test | 0.985| 150| 0.18|
|F-test | 0.785| 150| 0.08|
|F-test | 0.201| 150| 0.02|
|F-test | 0.999| 200| 0.18|
|F-test | 0.897| 200| 0.08|
|F-test | 0.311| 200| 0.02|
\end{Soutput}
\end{Schunk}

Figure~\ref{fig:syn-1} shows that the power curve estimated with MC simulations closely resembles the analytical power calculations:

\begin{figure}[!ht]
\centering
\includegraphics{article-exOneFig}
\caption{Comparison of power curves for F-test using closed-form formula (theoretical) and MC simulations in C-vine example.}
\label{fig:syn-1}
\end{figure}

\subsection{Estimation of NHANES mixtures with F-test for linear model}

In the second example, we generate predictors by resampling the National Health and Nutrition Examination Survey (NHANES) data. We download publicly available demographic, examination and laboratory data from https://www.cdc.gov/nchs/nhanes.htm. We combine data from the 2015-2016 and 2017-2018 cycles and adjust the survey weights as instructed in the official documentation https://wwwn.cdc.gov/nchs/nhanes/tutorials/module3.aspx. The data is included in our package as \code{"nhanes1518"}. We use a data subset that includes 5483 observations with no missing data and the following variables: age (\code{RIDAGEYR}), gender (\code{RIAGENDR}), household income (\code{INDHHIN2}), BMI (\code{BMXBMI}), Creatinine (\code{URXUCR}), and 16 phthalate metabolites. We choose this data subset because it has continuous, binary, and ordinal predictors and has extensive correlations between many of them. We log-transform the phthalate metabolites and Creatinine. We use the Bayesian Gaussian copula to simulate synthetic data that preserves the dependence structure. Figure~\ref{fig:resamp-cor} shows similar Spearman's correlation matrices between the original data (left) and a simulated data set of 500 observations (right):

\begin{Schunk}
\begin{Sinput}
R> set.seed(1)
R> data("nhanes1518", package = "mpower")
R> nhanes <- nhanes1518 %>%
+    select(starts_with("URX"), BMXBMI, INDHHIN2, RIDAGEYR, RIAGENDR) %>%
+    filter(complete.cases(.)) %>%
+    mutate(across(starts_with("URX"), log))
R> xmod <- MixtureModel(method = "estimation" , data = nhanes,
+    sbg_args = list(nsamp = 2000))
\end{Sinput}
\end{Schunk}

\begin{figure}[!ht]
\centering
\includegraphics{article-exTwoXFig}
\caption{Correlation matrices of mixture data from the NHANES 2015-2018 cycles (left) and of 500 observations simulated from an estimation mixture model (right).}
\label{fig:resamp-cor}
\end{figure}

We then define the mean function for a continuous outcome as a linear combination between age, gender (male), metabolites concentration (using Creatinine), and three correlated phthalate metabolites:

\begin{Schunk}
\begin{Sinput}
R> ymod <- OutcomeModel(sigma = 1, family = "gaussian",
+    f = "0.02*RIDAGEYR + 0.1*I(RIAGENDR==1) + 0.1*URXUCR +
+    0.15*URXMHH + 0.07*URXMOH - 0.1*URXMHP")
\end{Sinput}
\end{Schunk}

We can use the\code{I()} function to define effects of different levels in a categorical variable. We consider the F-test defined in the first example again for the power analysis using 1000 MC simulations for a sample size of 100:

\begin{Schunk}
\begin{Sinput}
R> set.seed(1)
R> power <- sim_power(xmod, ymod, imod, s = 1000, n = 100,
+    cores = 1, errorhandling = "stop")
\end{Sinput}
\end{Schunk}

The estimated SNR is between 0.27 and 0.28. The closed-form formula for the F-test for a linear model at the 0.05 significance level with 21 predictors, 100 observations and a SNR between 0.27 and 0.28 returns a power between 0.81 to 0.85. At the same significance level, the MC estimate is:

\begin{Schunk}
\begin{Sinput}
R> summary_df <- summary(power, crit = "pval", thres = 0.05, how = "lesser")
\end{Sinput}
\begin{Soutput}
*** POWER ANALYSIS SUMMARY ***
Number of Monte Carlo simulations: 1000
Number of observations in each simulation: 100
Data generating process estimated SNR: 0.28
Inference model: F-test
Significance criterion: pval

Significance threshold:  0.05

|       | power|
|F-test | 0.822|
\end{Soutput}
\end{Schunk}

We can also see how the power estimate varies with different significance thresholds (Figure~\ref{fig:thres-syn2}):

\begin{figure}[H]
\centering
\includegraphics{article-exTwoSummaryFig}
\caption{Power of the F-test for different significant thresholds.}
\label{fig:thres-syn2}
\end{figure}

We see from the two examples above that MC simulation is a reliable way to estimate power of a test.

\subsection{Resampling of NHANES with t-tests in linear model}

In this last example, we resample data from NHANES and calculate the power to detect the associations through conditional t-tests in a linear regression. We use a subset of the NHANES data available in \proglang{R} package \pkg{NHANES} [Pruim, 2015] for easy replication of the example. We recode ordinal variables \code{Education} and \code{HHIncome} as numeric variables in this example:

\begin{Schunk}
\begin{Sinput}
R> data("NHANES", package = "NHANES")
R> nhanes_demo <- NHANES %>%
+    select(BMI, Education, HHIncome, HomeOwn, Poverty, Age) %>%
+    filter(complete.cases(.)) %>%
+    mutate(Education = case_when(Education == "8th Grade" ~ 1,
+                                 Education == "9 - 11th Grade" ~ 2,
+                                 Education == "High School" ~ 3,
+                                 Education == "Some College" ~ 4,
+                                 Education == "College Grad" ~ 5),
+           HHIncome = case_when(HHIncome == "more 99999" ~ 12,
+                                HHIncome == "75000-99999" ~ 11,
+                                HHIncome == "65000-74999" ~ 10,
+                                HHIncome == "55000-64999" ~ 9,
+                                HHIncome == "45000-54999" ~ 8,
+                                HHIncome == "35000-44999" ~ 7,
+                                HHIncome == "25000-34999" ~ 6,
+                                HHIncome == "20000-24999" ~ 5,
+                                HHIncome == "15000-19999" ~ 4,
+                                HHIncome == "10000-14999" ~ 3,
+                                HHIncome == " 5000-9999" ~ 2,
+                                HHIncome == " 0-4999" ~ 1),
+           HomeOwn = ifelse(HomeOwn == "Own", 1, 0)
+           )
R> xmod <- MixtureModel(method = "resampling",
+    data = nhanes_demo %>% select(- BMI))
\end{Sinput}
\end{Schunk}

\begin{figure}[!ht]
\centering
\includegraphics{article-exThreeXFig}
\caption{Correlation matrices of data from the NHANES demographic data (left) and of 200 observations resampled from the same data (right).}
\label{fig:corr3}
\end{figure}

This example will do power analysis of the t-test for linear regression to identify an association between BMI and poverty status, conditional on age, education, household income, and home-owner status. NHANES data have survey sampling weights that maybe given to the model through the \verb|resamp_prob| argument. Though the resampling \class{MixtureModel} accepts non-numeric data, we re-code the ordinal variables above as integers to calculate all pair-wise correlations. Figure~\ref{fig:corr3} shows the Spearman correlation matrix of the original data compared to the resampled data. Note that several of the predictors have moderate to high correlations. Instead of using a string as we did in examples before, we will use a linear model fitted to the data as the input to the\code{OutcomeModel()}. Mathematically, we want the outcome model to be:

\begin{align}
y_{i,bmi}|. &\sim N(\mu_i, \sigma^2)\\
\mu_i &= \beta_0 + \beta_1 x_{i,pov} + \beta_2 x_{i, age} + \beta_3 x_{i,age}^2 + \beta_4 x_{i,inc} + \beta_5 x_{i, hom} + \\
& \beta_6 x_{i, edu} + \beta_{12} x_{i, poverty}x_{i,age} + \beta_{13} x_{i, poverty}x_{i,age}^2 + \beta_{14} x_{i, pov}x_{i,inc} +\\
& \beta_{15} x_{i, pov}x_{i,hom} + \beta_{16} x_{i, pov}x_{i,edu}
\label{eqn:bmi-lm}
\end{align}

where $\beta_j$'s and $\sigma$ use the estimates from the linear model. We then need to wrap the fitted linear model inside a function that takes a data frame \verb|x| and returns a vector of expected outcome so that it is a compatible input to\code{OutcomeModel()}:

\begin{Schunk}
\begin{Sinput}
R> lm_demo <- lm(BMI ~ Poverty*(poly(Age, 2) + HHIncome + HomeOwn + Education),
+    data = nhanes_demo)
R> f_demo <- function(x) {
+    predict.lm(lm_demo, newdata = x)
+  }
R> ymod <- OutcomeModel(f = f_demo, family = "gaussian", sigma = sigma(lm_demo))
\end{Sinput}
\end{Schunk}

We pass the ``true'' linear model to a GLM \class{InferenceModel}. Note that we must always refer to the outcome as \verb|y| in any formula we give to the function:

\begin{Schunk}
\begin{Sinput}
R> imod <- InferenceModel(model = "glm", family = "gaussian",
+    formula = y ~ Poverty*(poly(Age, 2) + HHIncome + HomeOwn + Education))
\end{Sinput}
\end{Schunk}

We run power curve for sample sizes 100, 1100, ..., 5100 (incremented by 1000 each):

\begin{Schunk}
\begin{Sinput}
R> set.seed(1)
R> curve <- sim_curve(xmod, ymod, imod, s = 500, n = seq(100, 6000, 1000),
+    cores = 1)
\end{Sinput}
\end{Schunk}

We plot the power curve:
\begin{figure}[H]
\centering
\begin{Schunk}
\begin{Sinput}
R> plot_summary(curve, crit = "pval", thres = 0.05, how = "lesser")
\end{Sinput}
\end{Schunk}
\includegraphics{article-exThreeFig}
\caption{Power curve of t-test for main effects and interactions (i.e., a predictor is ``significant'' if either the main effect or the interaction has p-value below the 0.05 significance level) in linear models of BMI accounting for age, education, homeownership, household income, and poverty status.}
\label{fig:syn-3}
\end{figure}

In the\code{summary()} function above, the significance criterion \verb|pval| tests for significance of any of the main effect of poverty or interactions between poverty and other covariates at the \verb|thres| = 0.05 significance level. The power curve is summarized in Figure~\ref{fig:syn-3}. We see that over 2100 observations are enough to detect the association between poverty and BMI given other covariates with 80\% power.

\section[An example in reproducing a study]{An example in reproducing a study} \label{sec:wu20}

A recent study [Wu \emph{et~al.}, 2020] considers the relationships between mixed exposures to various chemicals and obesity in children and adolescents using NHANES 2005-2010 data. Obesity in the study is measured by BMI z-scores and an indicator of BMI z-scores being over the $95^{th}$ percentile. All exposures are log transformed for the analysis. The results from various statistical models suggest that Dichlorophenol25 and MEP may be important risk factors for obesity. Suppose that we want to design a study to replicate these results, possibly in a specific population of interest. An important first step is to calculate a sample size to ensure that the study is well-powered and that we avoid dedicating resources to collecting more data than needed. We may specify the target power to be $.80$ and the target Type I error rate to be $0.05$, which are common practices. We also need to specify effect sizes appropriately. One may simply use the estimated effects of 2,5-DCP and MEP from the original study using the NHANES 2005-2010 data provided by the authors. Specifically, the original study finds one unit increase in Dichlorophenol25 (on log scale) is associated with a $e^{0.55} = 1.73$ times increase and MEP (on log scale) with a $e^{0.31} = 1.36$ times increase in the odds of obesity in an multivariate logistic regression adjusted for all covariates and chemicals \citep[see Additional file 1 Table S1]{wu2020nhanes}. We additionally hypothesize that there may be a modest interaction between MEP and Methylparaben. Thus, we will use the following conditional mean function for the binary response obesity:

\begin{equation}
logit(y_i \mid .) = 0.55x_{i,dichlorophenol25} + 0.31x_{i, MEP} + 0.2x_{i,MEP}x_{i,Methylparaben}
\label{eqn:obesity-glm}
\end{equation}

We conduct power analysis for three sample sizes $200, 500$, and $1000$ using both a logistic regression as well as BWS as inference models. We can summarize the power at each sample size by evaluating the p-values at the 0.05 significance level for each predictor separately in the logistic regression. Likewise, we look at the $95\%$ credible interval of the joint effect in the BWS model. The following code demonstrate how to do that with our package:

\begin{Schunk}
\begin{Sinput}
R> set.seed(1)
R> data_url <- paste0("https://static-content.springer.com/esm/",
+                     "art%3A10.1186%2Fs12940-020-00642-6/MediaObjects/",
+                     "12940_2020_642_MOESM2_ESM.xlsx")
R> nhanes <- openxlsx::read.xlsx(data_url)
R> chems <- c("UrinaryBisphenolA", "UrinaryBenzophenone3",
+             "Methylparaben", "Propylparaben",
+             "dichlorophenol25", "dichlorophenol24",
+             "MBzP", "MEP", "MiBP")
R> xmod <- MixtureModel(data = nhanes[, chems], method = "resampling")
R> ymod <- OutcomeModel(
+    f = "0.55*dichlorophenol25 + 0.31*MEP + 0.2*MEP*Methylparaben",
+    family = "binomial")
R> bws_logit <- InferenceModel(model = "bws", iter = 5000, chains = 1,
+    refresh = 0, family = "binomial")
R> glm_logit <- InferenceModel(model = "glm", family = "binomial")
R> n_cores <- 1
R> s <- 100
R> n <- c(200, 500, 1000)
R> bws_out <- sim_curve(xmod = xmod, ymod = ymod, imod = bws_logit,
+    s = s, n = n, cores = n_cores)
R> glm_out <- sim_curve(xmod = xmod, ymod = ymod, imod = glm_logit,
+    s = s, n = n, cores = n_cores)
\end{Sinput}
\end{Schunk}

Figure~\ref{fig:nhanes-example} shows that the power to detect the joint effect in BWS is much higher than the power to detect individual effects in a logistic regression for small samples. The weights from BWS, however, could be misleading since they erroneously show Dichlorophenol24, a compound closely related to Dichlorophenol25, having a moderate contribution to the joint effect. Thus, if the researchers are mainly interested in the joint effect of a chemical mixtures, they may save resources by collecting fewer observations and using an inference model such as BWS or GQC. However, if individual effects are of interest, a large sample is needed, and logistic regression may be preferred.

\begin{figure}[!ht]
\centering
\includegraphics{article-exNHANESFig}
\caption{Top: Power curves for each predictor from a logistic regression at a 0.05 significance level. Bottom left: Power curve of the joint effect from a BWS using a $95\%$ credible interval. The dotted red line represents the popular 0.8 threshold for having adequate power. Bottom right: The posterior median BWS weights of the predictors.}
\label{fig:nhanes-example}
\end{figure}

\section{Summary}

We provide an \proglang{R} package that allows researchers to flexibly conduct power analysis for a large class of models using MC simulations. It allows users to simulate mixed-scaled predictors with realistic dependence structures, which is important to account for when testing for individual effects. Another unique strength of the package is that it interfaces with several state-of-the-art inference models for exposure mixtures, so that researchers not experienced with \proglang{R} or programming can set up simulations in just a few lines of code.

%% -- Bibliography -------------------------------------------------------------
%% - References need to be provided in a .bib BibTeX database.
%% - All references should be made with \cite, \citet, \citep, \citealp etc.
%%   (and never hard-coded). See the FAQ for details.
%% - JSS-specific markup (\proglang, \pkg, \code) should be used in the .bib.
%% - Titles in the .bib should be in title case.
%% - DOIs should be included where available.

\clearpage
\begin{thebibliography}{38}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

[Arnold \emph{et~al.}, 2011]
Arnold BF, Hogan DR, Jr JMC, Hubbard AE (2011).
\newblock \enquote{Simulation Methods to Estimate Design Power: An Overview for
Applied Research.}
\newblock \emph{BMC Medical Research Methodology}, \textbf{11}(94).
\newblock \doi{https://doi.org/10.1186/1471-2288-11-94}.

[Bedford and Cooke, 2002]
Bedford T, Cooke RM (2002).
\newblock \enquote{Vines: A New Graphical Model for Dependent Random
\newblock \emph{The Annals of Statistics}, \textbf{30}(4), 1031--1068.

[Bien \emph{et~al.}, 2013]
Bien J, Taylor J, Tibshirani R (2013).
\newblock \enquote{A LASSO for Hierarchical Interactions.}
\newblock \emph{The Annals of Statistics}, \textbf{41}(3), 1111.
\newblock \doi{10.1214/13-AOS1096}.

[Bobb \emph{et~al.}, 2018]
Bobb JF, Henn BC, Valeri L, Coull BA (2018).
\newblock \enquote{Statistical Software for Analyzing the Health Effects of
Multiple Concurrent Exposures via Bayesian Kernel Machine Regression.}
\newblock \emph{Environmental Health}, \textbf{17}(67).
\newblock \doi{10.1186/s12940-018-0413-y}.

[Bobb \emph{et~al.}, 2015]
Bobb JF, Valeri L, Henn BC, Christiani DC, Wright RO, Mazumdar M, Godleski JJ,
Coull BA (2015).
\newblock \enquote{Bayesian Kernel Machine Regression for Estimating the Health
Effects of Multi-Pollutant Mixtures.}
\newblock \emph{Biostatistics}, \textbf{16}(3), 493--508.
\newblock \doi{10.1093/biostatistics/kxu058}.

[Champely, 2020]
Champely S (2020).
\newblock \emph{\pkg{pwr}: Basic Functions for Power Analysis}.
\newblock \proglang{R} package version 1.3-0,
\urlprefix\url{https://CRAN.R-project.org/package=pwr}.

[Cornea \emph{et~al.}, 2022]
Cornea EA, Chen LM, Qaqish BF, Chu H, Ibrahim JG (2022).
\newblock \emph{\pkg{JMdesign}: Joint Modeling of Longitudinal and Survival
Data - Power Calculation}.
\newblock \proglang{R} package version 1.2,
\urlprefix\url{https://CRAN.R-project.org/package=JMdesign}.

[Corporation and Weston, 2022]
Corporation M, Weston S (2022).
\newblock \emph{\pkg{doSNOW}: Foreach Parallel Adaptor for the \pkg{snow}
\newblock \proglang{R} package version 1.0.20,
\urlprefix\url{https://CRAN.R-project.org/package=doSNOW}.

[Czanner \emph{et~al.}, 2008]
Czanner G, Sarma SV, Eden UT, Brown EN (2008).
\newblock \enquote{A Signal-to-Noise Ratio Estimator for Generalized Linear
Model Systems.}
\newblock \emph{Proceedings of the World Congress on Engineering}, \textbf{2}.

[Eaton, 2007]
Eaton ML (2007).
\newblock \enquote{Multivariate Statistics: A Vector Space Approach.}
\newblock \emph{Institute of Mathematical Statistics Lecture Notes - Monograph
Series}, \textbf{53}, 512.
\newblock \doi{10.1214/lnms/1196285102}.

[Fan, 2017]
Fan FY (2017).
\newblock \emph{\pkg{powerAnalysis}: Power Analysis in Experimental Design}.
\newblock \proglang{R} package version 0.2.1,
\urlprefix\url{https://CRAN.R-project.org/package=powerAnalysis}.

[Ferrari and Dunson, 2020]
Ferrari F, Dunson DB (2020{\natexlab{a}}).
\newblock \enquote{Bayesian Factor Analysis for Inference on Interactions.}
\newblock \emph{Journal of the American Statistical Association},
\textbf{0}(0), 1--12.
\newblock \doi{10.1080/01621459.2020.1745813}.

[Ferrari and Dunson, 2020]
Ferrari F, Dunson DB (2020{\natexlab{b}}).
\newblock \enquote{Identifying Main Effects and Interactions among Exposures
Using Gaussian Processes.}
\newblock \emph{The Annals of Applied Statistics}, \textbf{14}(4), 1743--1758.
\newblock \doi{https://doi.org/10.1214/20-AOAS1363}.

[Gelman and Hill, 2007]
Gelman A, Hill J (2007).
\newblock \emph{Data Analysis Using Regression and Multilevel/Hierarchical
\newblock Cambridge University Press.

[Green and MacLeod, 2016]
Green P, MacLeod CJ (2016).
\newblock \enquote{\pkg{SIMR}: An \proglang{R} Package for Power Analysis of
Generalized Linear Mixed Models by Simulation.}
\newblock \emph{Methods in Ecology and Evolution}, \textbf{7}(4), 493--498.
\newblock \doi{https://doi.org/10.1111/2041-210X.12504}.

[Hamra and Buckley, 2018]
Hamra GB, Buckley JP (2018).
\newblock \enquote{Environmental Exposure Mixtures: Questions and Methods to
Address Them.}
\newblock \emph{Current Epidemiology Reports}, \textbf{5}(2), 160--165.
\newblock \urlprefix\url{https://doi.org/10.1007/s40471-018-0145-0}.

[Hamra \emph{et~al.}, 2021]
Hamra GB, MacLehose RF, Croen L, Kauffman EM, Newschaffer C (2021).
\newblock \enquote{Bayesian Weighted Sums: a Flexible Approach to Estimate
Summed Mixture Effects.}
\newblock \emph{International Journal of Environmental Research and Public
Health}, \textbf{18}(4), 1373.
\newblock \doi{10.3390/ijerph18041373}.

[Hoeting \emph{et~al.}, 1999]
Hoeting JA, Madigan D, Raftery AE, Volinsky CT (1999).
\newblock \enquote{Bayesian Model Averaging: A Tutorial.}
\newblock \emph{Statistical Science}, \textbf{14}(4), 382--401.
\newblock ISSN 08834237.
\newblock \urlprefix\url{http://www.jstor.org/stable/2676803}.

[Hoff, 2018]
Hoff P (2018).
\newblock \emph{\pkg{sbgcop}: Semiparametric Bayesian Gaussian Copula
Estimation and Imputation}.
\newblock \proglang{R} package version 0.980,
\urlprefix\url{https://CRAN.R-project.org/package=sbgcop}.

[Hoff, 2007]
Hoff PD (2007).
\newblock \enquote{Extending the Rank Likelihood for Semiparametric Copula
\newblock \emph{The Annals of Applied Statistics}, \textbf{1}(1), 265 -- 283.
\newblock \doi{10.1214/07-AOAS107}.
\newblock \urlprefix\url{https://doi.org/10.1214/07-AOAS107}.

[Izmirlian, 2021]
Izmirlian G (2021).
\newblock \emph{\pkg{PwrGSD}: Power in a Group Sequential Design}.
\newblock \proglang{R} package version 2.3.6,
\urlprefix\url{https://CRAN.R-project.org/package=PwrGSD}.

[Joe, 2006]
Joe H (2006).
\newblock \enquote{Generting Random Correlation Matrices Based on Partial
\newblock \emph{Journal of Multivariate Analysis}, \textbf{97}, 2177--2189.

[Joubert \emph{et~al.}, 2022]
Joubert BR, Kioumourtzoglou MA, Chamberlain T, Chen HY, Gennings C, Turyk ME,
Miranda ML, Webster TF, Ensor KB, Dunson DB, Coull BA (2022).
\newblock \enquote{Powering Research through Innovative Methods for Mixtures in
Epidemiology (PRIME) Program: Novel and Expanded Statistical Methods.}
\newblock \emph{International Journal of Environmental Research and Public
Health}, \textbf{19}(3).
\newblock ISSN 1660-4601.
\newblock \doi{10.3390/ijerph19031378}.
\newblock \urlprefix\url{https://www.mdpi.com/1660-4601/19/3/1378}.

[Keil, 2021]
Keil A (2021).
\newblock \emph{\pkg{qgcomp}: Quantile G-Computation}.
\newblock \proglang{R} package version 2.7.0,
\urlprefix\url{https://CRAN.R-project.org/package=qgcomp}.

[Keil \emph{et~al.}, 2020]
Keil AP, Buckley JP, O'Brien KM, Ferguson KK, Zhao S, White AJ (2020).
\newblock \enquote{A Quantile-Based G-Computation Approach to Addressing the
Effects of Exposure Mixtures.}
\newblock \emph{Environmental Health Perspectives}, \textbf{128}(4).
\newblock \doi{10.1289/EHP5838}.

[Kooperberg and Hsu, 2015]
Kooperberg C, Hsu L (2015).
\newblock \emph{\pkg{powerGWASinteraction}: Power Calculations for GxE and GxG
Interactions for GWAS}.
\newblock \proglang{R} package version 1.1.3,
\urlprefix\url{https://CRAN.R-project.org/package=powerGWASinteraction}.

[LeBeau, 2022]
LeBeau B (2022).
\newblock \emph{\pkg{simglm}: Simulate Models Based on the Generalized Linear
\newblock \proglang{R} package version 0.8.9,
\urlprefix\url{https://CRAN.R-project.org/package=simglm}.

[Lewandowski \emph{et~al.}, 2009]
Lewandowski D, Kurowicka D, Joe H (2009).
\newblock \enquote{Generating Random Correlation Matrices Based on Vines and
Extended Onion Method.}
\newblock \emph{Journal of Multivariate Analysis}, \textbf{100}(9), 1989--2001.
\newblock \doi{https://doi.org/10.1016/j.jmva.2009.04.008}.

[Lim and Hastie, 2015]
Lim M, Hastie T (2015).
\newblock \enquote{Learning Interactions via Hierarchical Group-LASSO
\newblock \emph{Journal of Computational and Graphical Statistics},
\textbf{24}(3), 627--654.
\newblock \doi{doi: 10.1080/10618600.2014.938812}.

[McCullagh and Nelder, 1983]
McCullagh P, Nelder JA (1983).
\newblock \emph{Generalized Linear Models}.
\newblock Chapman and Hall.

[Microsoft} and Weston, 2020]
{Microsoft}, Weston S (2020).
\newblock \emph{\pkg{foreach}: Provides Foreach Looping Construct}.
\newblock \proglang{R} package version 1.5.1,
\urlprefix\url{https://CRAN.R-project.org/package=foreach}.

[Morgan-Wall and Khoury, 2021]
Morgan-Wall T, Khoury G (2021).
\newblock \enquote{Optimal Design Generation and Power Evaluation in
\proglang{R}: The \pkg{skpr} Package.}
\newblock \emph{Journal of Statistical Software}, \textbf{99}(1), 1–36.
\newblock \doi{10.18637/jss.v099.i01}.
\newblock
\urlprefix\url{https://www.jstatsoft.org/index.php/jss/article/view/v099i01}.

[Nguyen, 2022]
Nguyen PH (2022).
\newblock \emph{bws: Bayesian Weighted Sums}.
\newblock \proglang{R} package version 0.1.0,
\urlprefix\url{https://CRAN.R-project.org/package=bws}.

[Poworoznek, 2020]
Poworoznek E (2020).
\newblock \emph{\pkg{infinitefactor}: Bayesian Infinite Factor Models}.
\newblock \proglang{R} package version 1.0,
\urlprefix\url{https://CRAN.R-project.org/package=infinitefactor}.

[Pruim, 2015]
Pruim R (2015).
\newblock \emph{\pkg{NHANES}: Data from the US National Health and Nutrition
Examination Study}.
\newblock \proglang{R} package version 2.1.0,
\urlprefix\url{https://CRAN.R-project.org/package=NHANES}.

[Raftery \emph{et~al.}, 2021]
Raftery A, Hoeting J, Volinsky C, Painter I, Yeung KY (2021).
\newblock \emph{\pkg{BMA}: Bayesian Model Averaging}.
\newblock \proglang{R} package version 3.18.15,
\urlprefix\url{https://CRAN.R-project.org/package=BMA}.

[Sun \emph{et~al.}, 2013]
Sun Z, Tao Y, Li S, Ferguson KK, Meeker JD, Park SK, Batterman SA, Mukherjee B
\newblock \enquote{Statistical Strategies for Constructing Health Risk Models
with Multiple Pollutants and Their Interactions: Possible Choices and
\newblock \emph{Environmental Health}, \textbf{12}(1).
\newblock \doi{10.1186/1476-069X-12-85}.

[Wu \emph{et~al.}, 2020]
Wu B, Jiang Y, Jin X, He L (2020).
\newblock \enquote{Using Three Statistical Methods to Analyze the Association
Between Exposure to 9 Compounds and Obesity in Children and Adolescents:
NHANES 2005-2010.}
\newblock \emph{Environmental Health}, \textbf{19}(94).
\newblock \doi{https://doi.org/10.1186/s12940-020-00642-6}.

\end{thebibliography}

%% -----------------------------------------------------------------------------

\end{document}