Power of a Interaction term in R

Question

I have analysed a dataset with a linear regression model, including an interaction term between a binary variable and a continuous variable. The interaction was significant. Afterwards, I have fitted 2 separate models of the continuous variable for each of the 2 groups of the binary variable. The 2 slopes have different signs and one of the two slopes is significant.

I need to calculate the power of the significance of this interaction. I prefer to do this with an R function.

Answer 1

I do not think there is an R function to calculate this power. My usual suggestion in such cases is to simulate (e.g., here or here or here). Specifically, pick a sample size $n$, simulate your covariate, your binary variable, the dependence of your outcome on both and the noise. Fit your model and assess whether $p<\alpha$. Do this many times and see how often you detect a significant effect. If you want to determine the necessary sample size to reach a target power, change the $n$ until you reach the power you want.

Yes, this means you have to think about many assumptions (e.g., the distributions of your IVs, whether they are independent, whether linearity makes sense, whether your errors will really be homoskedastic, ...). I would argue that this is a feature, not a bug. You will definitely gain more understanding of your problem than if you found a ready-made power calculator, and it's better to gain this understanding before you run your study than afterwards.

(You are asking before you run your study, for sample size determination, right? So-called "post-hoc power" is meaningless.)

One additional advantage is that the flexibility of this approach allows you to vary your assumptions and analyze whether the power degrades. Or you could simulate other data problems, like missing data etc.

Here is some very simple R code, where I am assuming $n=50$, group membership that is random 50-50, a uniformly distributed covariate that is independent of group membership, and an outcome

$$y=\text{covariate}+\text{groups}*\text{covariate}+\epsilon$$

with $\epsilon\sim N(0,1)$. In this case, power comes out to be 0.166. If you want the standard $\beta=0.80$, you can increase $n$ or change your model assumptions - and think about whether the changed assumptions make sense.

n_sims <- 1000

nn <- 50
alpha <- 0.05
result <- rep(FALSE,n_sims)
pb <- winProgressBar(max=n_sims)
for ( ii in 1:n_sims ) {
    setWinProgressBar(pb,ii,paste(ii,"of",n_sims))
    set.seed(ii)    # for reproducibility

    # this is where the assumptions enter
    covariate <- runif(nn)
    groups <- runif(nn)<0.5
    outcome <- covariate+groups*covariate+rnorm(nn,1)

    # fit the model, determine whether an effect was found
    model <- lm(outcome~groups*covariate)
    p_value <- anova(model,update(model,.~.-groups:covariate))[2,6]
    result[ii] <- p_value<alpha
}
close(pb)

sum(result)/n_sims

Answer 2

The R package InteractionPoweR can do a power analysis for this sort of regression.

install.packages("devtools")
devtools::install_github("dbaranger/InteractionPoweR")
library(InteractionPoweR)

test_power<-power_interaction(
  n.iter = 1000,            # number of simulations per unique combination of input parameters
  alpha = 0.05,             # alpha, for the power analysis
  N = 350,                  # sample size
  r.x1x2.y = .15,           # interaction effect to test (correlation between x1*x2 and y)
  r.x1.y = .2,              # correlation between x1 and y
  r.x2.y = .1,              # correlation between x2 and y
  r.x1.x2 = .2,             # correlation between x1 and x2
  k.x1 = 2,                 # x1 is binary
  seed = 581827             # seed, for reproducibility
  adjust.correlations = TRUE)     # Adjust correlations