Is it possible to calculate the power for a specific beta coefficient in a multiple linear regression?

I want to conduct a multiple linear regression where I am only interested in interpreting one beta coefficient but I will adjust for two additional covariates (i.e. three variables total). Is it possible to calculate the power I would have given a sample size of $$n = 150$$ to detect a standardised beta coefficient of (let's say) $$0.3$$ in my explanatory variable of interest given these two other covariates? If so, what additional information do I require and is there an R package or code that will run the calculations for me?

Yes, this is possible.

Eric Vittinghoff and coauthors include sample size equations at the end of the chapters in their book Regression Methods in Biostatistics.

The power $$\gamma$$ attained given a sample size and smallest detectable effect is

$$\gamma=1-\Phi\left[z_{1-\alpha / 2}-\left|\beta_j^a\right| \sigma_{x_j} \sqrt{n\left(1-\rho_j^2\right)} / \sigma_{y \mid \mathbf{x}}\right]$$

Here:

• $$\Phi$$ is the standard normal CDF
• $$z_{1-\alpha/2}$$ is the $$1-\alpha/2$$ quantile of a standard normal. If you use $$\alpha=0.05$$ then this quantity is 1.96
• $$\beta^\alpha_j$$ is the adjusted direct effect of a unit increase in the covariate of interest
• $$\sigma_x$$ is the standard deviation of the predictor for which you would like to perform the hypothesis test. This can always be rescaled to be 1 in your data by standardizing the variable if it is continuous.
• $$n$$ is your sample size
• $$(1-\rho^2)$$ is the variance inflation factor, and $$\rho$$ is the multiple correlation between your covariate of interest and the remaining 2 covariates.
• $$\sigma_{y\mid x}$$ is the residual standard deviation.

Here is an example in R. I assume that the variable under investigation is a binary exposure with 50% probability of success. This variable is randomly assigned and so $$\rho=0$$. This is equivalent to a two sample t-test.

Here is the code to produce the power calculation and a simulation to verify the power calculation empirically.

library(tidyverse)

a <- 0.05
n <- 150
sigma <- 0.5
s <- 1
beta <- 0.5

pwr <- 1 - pnorm(qnorm(1-a/2) - abs(beta)*sigma*sqrt(n)/s)
pwr
#> [1] 0.864747

sims <- map_dfr(1:1000, ~{
x <- rbinom(n, 1, 0.5)
y <- beta*x + rnorm(n, 0, s)
fit <- lm(y~x)

broom::tidy(fit) %>%
filter(term=='x')

})

mean(sims$p.value<0.05) #> [1] 0.854  Created on 2024-03-27 with reprex v2.0.2 The difference between the simulation and power calculation is small enough to be negligible in my opinion. • While the equation is informative, I think OP is looking for an applied solution. Perhaps there is a way to demonstrate this in R as they suggest. Mar 27 at 4:10 • +1 Demetri Pananos, can you amplify a little what "on the natural scale" means? I suspect there may be variations in language between disciplines that may make that phrase less understandable to some. Mar 27 at 14:52 • By "natural" scale, I mean that$\beta$is the associated effect of the covariate$x$in which$x$is unchanged due to scaling or transformation. It would be the expected coefficient were you to use "raw"$x$in your regression. The interpretation of$\beta$is changed were you to standardize$x$, for example, as then it would be the associated change in$y$given a 1 standard deviation change in$x$. Mar 27 at 15:13 • @LukasLohse Yes, in the limit the test statistic for coefficients is a z statistic. The quantity in the square root is the variance inflation factor. If you have a better name as opposed to "natural scale" I'm happy to hear it. Additionally, the$\beta$should not be the assumed coefficient, it should be the smallest effect of interest. Further discussion on this point can get lengthy, and so if you like I can maybe find some resources on this. Mar 27 at 17:18 • @ShawnHemelstrand I've added one Mar 27 at 17:30 Yes, it's possible. You can use for instance the WebPower package's wp.regression function. It requires the effect size in f2 form though. Standardized beta of 0.3 can be seen as medium effect but on the smaller side of medium, so a corresponding f2 would be something like ~ 0.10. The test is run as the difference in f2 between a model with covariates only and a model with your predictor of interest and the covariates. p1 is the number of predictors in full model and p2 is the number of predictors in model with just the covariates. So, with 2 covariates library(WebPower) fpower<-wp.regression(n=150, p1=3, p2=2, f2=0.10) fpower Power for multiple regression n p1 p2 f2 alpha power 150 3 2 0.1 0.05 0.9704481  Edit. I didn't know about the pwrss packaged linked by Mathemagician777, it seems better for this purpose as you can directly put in the standardized coefficient. • Wouldn't this not look at the individual power of the coefficients in question though? This seems to be a more model-based version of power if my understanding is right. Mar 27 at 4:07 • It does test the power for the individual coefficient, it just conceptualizes its effect size as the r-squared increase when adding the target predictor (vs. without it). Mar 27 at 8:25 • Ah I just read the documentation and noticed it was based on Cohen's$f^2$so that makes more sense now. Mar 27 at 9:08 Quick Answer One can just use a simulated power analysis to specify every piece of the model that they are weighting if they are interested in a single coefficient. For example, if we want to specify the exact beta coefficients of each predictor, their standard error, the sample size, or any other metric, we could just simulate it in R directly. This video gives an excellent example for both simple regression and multiple linear regression. The R code for their simulations can be found here. Worked Example A shortened version I just came up with can be found below. First you setup the function which randomly generates the data necessary for the simulation. Here the arguments are as follows: • n: sample size • beta0: the intercept • beta1: the slope of $$x_1$$ • beta2: the slope of $$x_2$$ • error: the residual standard error of the model Which you can create with this code (and can be modified to be more specific depending on your use case). For example, here I assume the SE for each predictor is the same, but this too can be modified: #### Data Generator Function #### gen.data <- function(n, beta0, beta1, beta2, error){ x1 <- rnorm(n) x2 <- rnorm(n) e <- rnorm(n, error) y <- (beta0) + (beta1 * x1) + (beta2 * x2) + e df <- data.frame(x1, x2, y) return(df) }  Then you setup the simulation parameters by running $$B = 1000$$ runs, an alpha level of $$a = .05$$, and three numeric vectors to store the $$p$$ values for each coefficient. #### Setup Simulation #### boot <- 1000 alpha <- 0.05 # significance level p0 <- p1 <- p2 <- numeric(boot) n <- 100 # sample size beta0 <- 0 # intercept beta1 <- .2 # slope 1 beta2 <- .2 # slope 2 error <- 10 # RSE  Then you just run a for-loop which first generates the data based off our parameters, fits the data to a regression, and saves the $$p$$ values from each fit. #### Run Sim #### for(i in 1:boot){ df <- gen.data( n = n, beta0 = beta0, beta1 = beta1, beta2 = beta2, error = error ) fit <- lm(y ~ x1 + x2, df) p_values <- summary(fit)$coefficients[,4]
p0[i] <- p_values[1]
p1[i] <- p_values[2]
p2[i] <- p_values[3]
}


You can then calculate the power by taking the mean of all the values which are below your cutoff.

#### Calculate Power ####
power0 <- mean(p0 < alpha)
power1 <- mean(p1 < alpha)
power2 <- mean(p2 < alpha)

#### Check Power ####
power1
power2


This should print off the slopes of $$\beta_1$$ and $$\beta_2$$ accordingly. For more info, check out the linked video for a full explanation.

I will summarize two different approaches utilizing the pwrss R package. A comprehensive tutorial was provided by @Mathemagician777 in one of the comments. The challenging aspect involves manually calculating the R-squared value, which I wish pwrss had automated.

If the main predictor is continuous:

library(pwrss)

k <- 3 # number of predictors
m <- 1 # number of predictors under scrutiny
beta1 <- 0.30 # the standardized regression coefficient
r2 <- beta1^2 # the R-squared is at least beta1^2
n <- 150 # sample size

pwrss.t.reg(beta1 = beta1, k = k, r2 = r2, n = n)

# Linear Regression Coefficient (t Test)
# H0: beta1 = beta0
# HA: beta1 != beta0
# ------------------------------
#  Statistical power = 0.969
#  n = 150
# ------------------------------
# Alternative = “not equal”
# Degrees of freedom = 146
# Non-centrality parameter = 3.852
# Type I error rate = 0.05
# Type II error rate = 0.031

pwrss.f.reg(f2 = r2 / (1 - r2), k = k, m = m, n = n)

# Hierarchical Linear Regression (F test)
# R-squared Change
# H0: r2 = 0
# HA: r2 > 0
# ------------------------------
#  Statistical power = 0.969
#  n = 150
# ------------------------------
# Numerator degrees of freedom = 1
# Denominator degrees of freedom = 146
# Non-centrality parameter = 14.835
# Type I error rate = 0.05
# Type II error rate = 0.031


If the main predictor is binary:

p <- 0.50 # proportion of assignment
varx <- p*(1-p) # variance of the binary predictor
r2 <- beta1^2 * varx # the R-squared is at least beta^2 * varx

pwrss.t.reg(beta1 = beta1, sdx = sqrt(varx), k = k, r2 = r2, n = n)

# Linear Regression Coefficient (t Test)
# H0: beta1 = beta0
# HA: beta1 != beta0
# ------------------------------
#  Statistical power = 0.455
#  n = 150
# ------------------------------
# Alternative = “not equal”
# Degrees of freedom = 146
# Non-centrality parameter = 1.858
# Type I error rate = 0.05
# Type II error rate = 0.545

pwrss.f.reg(f2 =  r2 / (1 - r2), k = k, m = m, n = n)

# Hierarchical Linear Regression (F test)
# R-squared Change
# H0: r2 = 0
# HA: r2 > 0
# ------------------------------
#  Statistical power = 0.455
#  n = 150
# ------------------------------
# Numerator degrees of freedom = 1
# Denominator degrees of freedom = 146
# Non-centrality parameter = 3.453
# Type I error rate = 0.05
# Type II error rate = 0.545


I hope it helps.

OP is asking specifically about performing a power analysis when the hypothesis is about a standardized beta in a multiple regression. Some of the other answers assume that the variables in the model are all uncorrelated, which may or may not be the case. The tricky part is that instead of specifying the standardized beta, most packages expect correlation coefficients ($$r$$) or some other effect size (and figuring out that effect size - $$R^2$$ or $$f^2$$ - in many cases requires first having all the correlation coefficients in hand). One approach then is to solve for what $$r$$ would yield the desired beta. Other than beta, we also need to know the correlation ($$r$$) between all other variables in the model ($$y$$, $$x1$$, $$x2$$, and $$x3$$).

# correlations between all the variables
yx2=.2
yx3=.3
x1x2=.05
x1x3=.06
x2x3=.1

beta=.1 # this is the beta we are interested in

# finds the difference between the desired beta and the
#   resulting beta for a given correlation between y and x1

find.beta = function(yx1,yx2,yx3,x1x2,x1x3,x2x3,beta){

cor.mat = matrix(data = c(1   , yx1,yx2,  yx3,
yx1, 1   ,x1x2, x1x3,
yx2,  x1x2, 1  , x2x3,
yx3,  x1x3,  x2x3, 1) ,nrow = 4,ncol = 4,byrow = T)
betas= solve(cor.mat[-1,-1],cor.mat[1,-1])
return(betas[1] - beta)
}

# solve for the correlation that will yield the desired beta
new.r = uniroot(f = find.beta,lower = -1,upper = 1,yx2=yx2,yx3=yx3,
x1x2=x1x2,x1x3=x1x3,x2x3=x2x3,beta=beta)

new.r$root [1] 0.125  In this example, yx1 ($$r$$) of 0.125 yields a b1 of 0.1. Now we do a standard power analysis using the correlation we found: install.packages('pwr2ppl') pwr2ppl::MRC(ry1 = new.r$root,ry2 = yx2,ry3 = yx3,r12 = x1x2,
r13 = x1x3,r23 =x2x3,n = 100,alpha = 0.05 )

Sample size is 100
Power R2 = 0.888
Power b1 = 0.179
Power b2 = 0.408
Power b3 = 0.817


With a sample size of N=100, we have a power of 0.179 to detect our beta of 0.1.

Edit

If you instead have the standardized effect size for all three variables, you can also solve for all three regression coefficients at once, and then compute power the same way.

x1x2=.05
x1x3=.06
x2x3=.1

betas = c(0.1,0.1672727,0.2772727)

cor.mat = matrix(data = c(1   ,x1x2, x1x3,
x1x2, 1  , x2x3,
x1x3,  x2x3, 1) ,nrow = 3,ncol = 3,byrow = T)

xy.cor = solve(solve(cor.mat),betas)
xy.cor
[1] 0.125 0.200 0.300

pwr2ppl::MRC(ry1 = xy.cor[1],ry2 = xy.cor[2],ry3 = xy.cor[3],r12 = x1x2,
r13 = x1x3,r23 =x2x3,n = 100,alpha = 0.05 )
Sample size is 100
Power R2 = 0.888
Power b1 = 0.179
Power b2 = 0.408
Power b3 = 0.817


Yielding the same power as previously