I am trying to replicate the following graph that I found in a paper: Power Graph

After emailing the author, he described the process to generate this graph by saying:

"...we did simulations with these different SD values superimposed on the slopes and then computed power"

I've tried to replicate it but I'm having trouble getting similar results. So far what I've done is generated many linear models at each slope steepness (4xSD, 2xSD etc.) and for each number of sample size (3, 4 etc.) and then calculated the power of these regressions using the pwr.f2.test function in the pwr package in R.

For example, if I wanted to calculate the first square point on the graph (around (3,0.18)) I generated many sets of 3 points normally distributed around a line of a particular slope (determined by the SD chosen, 4xSD etc.). I then calculated the R^2 value of the fit of the points on the line and calculated the mean R^2 of all of the generated fits.

I then used the mean R^2 to calculate mean effect size (f^2 = R^2/(1-R^2)) and added it to the power equation: pwr.2f.test(u, v, f2, sig) where u is the numerator degree of freedom (p - 1 = 2 - 1), v is the denominator df (n - p = 3 - 2 = 1), f2 is mean effect size and sig is level of significance (0.05) in this example.

Is this correct? The power values I get, however tend to 1 much quicker than this graph. I've also tried calculating the power for each model generated individually and getting the mean of all the powers but these results seem even further away from the graph. Sorry if this is incoherent, I'm relatively new to statistics.


Turns out I was doing the method right and I was messing up which axis I was doing the regression on! But I'll leave this up in case other people want to search for this method.

Below is my results plotted in base R: Replicate graph

Codesnippet of what i did:

CalcPower <- function(num, sd, slope){ #input sample size, standard deviation and slope you want to calculate power for
  vec <- vector(mode = "numeric", length = 5000) #initialize vector to hold simulated R^2 results
  y <- seq(0, length.out = num, by = 1) #setting axis
  for(x in 1:5000){
    a <- rnorm(n=num, sd=sd) + seq(0, length.out = num, by = -0.5)*slope #simulating variance around line of particular slope
    vec[x] <- summary(lm(a ~ y))$adj.r.squared #storing results to vector
  f2 <- mean(vec)/(1-mean(vec)) #Calculating mean effect size from R^2
  return(pwr.f2.test(u = 1, v = num-2, f2 = f2, sig.level = 0.05)) #Calculating power

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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