11
$\begingroup$

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?

$\endgroup$
2

5 Answers 5

12
$\begingroup$

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.

$\endgroup$
10
  • 1
    $\begingroup$ 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. $\endgroup$ Mar 27 at 4:10
  • 2
    $\begingroup$ +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. $\endgroup$
    – Alexis
    Mar 27 at 14:52
  • $\begingroup$ 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$. $\endgroup$ Mar 27 at 15:13
  • 1
    $\begingroup$ @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. $\endgroup$ Mar 27 at 17:18
  • 2
    $\begingroup$ @ShawnHemelstrand I've added one $\endgroup$ Mar 27 at 17:30
10
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Mar 27 at 4:07
  • 1
    $\begingroup$ 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). $\endgroup$
    – Sointu
    Mar 27 at 8:25
  • 1
    $\begingroup$ Ah I just read the documentation and noticed it was based on Cohen's $f^2$ so that makes more sense now. $\endgroup$ Mar 27 at 9:08
7
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.