# Calculating the (mean) power of simulated linear regressions using pwr package in R

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

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:

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
}