knitr::opts_chunk$set(echo = TRUE)
This vignette describes how to run power analyses for 2-way interactions that include covariates. It assumes you are already familiar with power analyses for 2-way interactions. If you aren't, check out our tutorial paper or the main vignette.
A two-way interaction analysis with two covariates (for example, but any number is allowable) take the form:
$$ Y \sim \beta_0 + X_1\beta_1 + X_2\beta_2 + X_1X_2\beta_3+ C_1\beta_5 + C_2\beta_6 + C_1X_1\beta_7 + C_1X_2\beta_8 + C_2X_1\beta_9 + C_2X_2\beta_{10} + \epsilon $$
Where $C_i$ are the covariates, and $C_iX_i$ are interactions between the covariates and the main variables of interest. The inclusion of the $C_iX_i$ terms may surprise some users. However, it is well established that failure to include them can lead to false-positive results, because of omitted variable bias. This can occur not only when $C_iX_i$ is independent predictor of $Y$, but also simply when $cor(X_i,C_i)$ is non-zero (i.e., when a covariate is correlated with either of the main interaction terms). A few relevant citations: Hull et al., 1992, Yzerbyt et al., 2004, Keller, 2014, and blogs: Baranger, and Simonsohn.
First, the function generate.interaction.cov.input()
is used to setup the input correlations for the power analysis. The only input is the number of covariates. This generates a named list. All correlations are by default 0
and all reliabilites are 1
.
library(InteractionPoweR) power.input = generate.interaction.cov.input(c.num = 2) # number of covariates head(power.input$correlations) head(power.input$reliability)
Modify the list as-needed for your own power analysis:
power.input$correlations$r.y.x1 = .3 power.input$correlations$r.y.x2 = .2 power.input$correlations$r.y.c1 = .1 power.input$correlations$r.y.c2 = .2 power.input$correlations$r.y.x1x2 = seq(0.01,.1,.01) # inputs can be single values or vectors power.input$correlations$r.x1.x2 = .2 power.input$correlations$r.c1.c2 = .3 power.input$correlations$r.x1.c1 = .1 power.input$correlations$r.x2.c1 = .1
Now we run the power analysis with the power_interaction_r2_covs()
function:
pwr.analysis = power_interaction_r2_covs(cov.input =power.input, # our parameters N = c(1000) ) pwr.analysis
As always, we can plot the power curve:
plot_power_curve(power_data = pwr.analysis,x = "r.y.x1x2")
And we can estimate where the power curve intersects some desired level of power:
power_estimate(power_data = pwr.analysis,x = "r.y.x1x2",power_target = 0.8)
Even though analytic power analyses are fast, with so many parameters we can easily reach a point where the analysis will take a long time to run. Say for instance we are interested in testing out 10 different values of 4 different parameters - that's 10,000 power analyses. At that point, it can be helpful to run things in parallel. Parallel analyses can be run using the cl
flag. cl=4
is probably a reasonable value for most personal computers.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.