I'm running power simulations for ANCOVA using the guidance provided here: https://egap.org/content/power-analysis-simulations-r

Using Rstudio Version 1.1.456, I'm experiencing the following error ("Error in summary(fit.sim)$coefficients[2, 4] : subscript out of bounds") when the number of sims > 100. I'm hoping someone can help me with this error as I'd like to use sims <- 1000 (or greater) for the power simulations. Thanks in advance.

rm(list=ls())

possible.ns <- seq(from=10, to=250, by=5)

powers <- rep(NA, length(possible.ns))

powers.cov <- rep(NA, length(possible.ns)) # Need a second empty vector

alpha <- 0.05

sims <- 1000

### Outer loop to vary the number of subjects

for (j in 1:length(possible.ns)){

N <- possible.ns[j] # Pick the jth value for N

significant.experiments <- rep(NA, sims) # Empty object to count significant experiments

significant.experiments.cov <- rep(NA, sims) # Need a second empty vector here too

### Inner loop to conduct experiments "sims" times over for each N

for (i in 1:sims){

age <- sample(x=18:40, size=N, replace=TRUE) # Generate "age" covariate

baseline <- sample(x=1:15, size=N, replace=TRUE) # Generate "baseline" covariate

effectofage <- 0.5 # Hypothesize the "effect" of age on dv

effectofbaseline <- 0.8 # Hypothesize the "effect" of baseline on dv

```
## Hypothesize Control Outcome as a function of age, baseline, and error
Y0 <- effectofage*age + effectofbaseline*baseline + rnorm(n=N, mean=6, sd=1.5)
## This is all the same ##
tau <- 2 # Hypothesize treatment effect
Y1 <- Y0 + tau # treatment potential outcome
Z.sim <- rbinom(n=N, size=1, prob=.5) # Do a random assignment
Y.sim <- Y1*Z.sim + Y0*(1-Z.sim) # Reveal outcomes according to assignment
fit.sim <- lm(Y.sim ~ Z.sim) # Do analysis (Simple regression)
## This is the novel analysis -- including three covariates to increase precision ##
fit.sim.cov <- lm(Y.sim ~ Z.sim + age + baseline)
## extract p-values and calculate significance ##
p.value <- summary(fit.sim)$coefficients[2,4]
p.value.cov <- summary(fit.sim.cov)$coefficients[2,4]
significant.experiments[i] <- (p.value <= alpha)
significant.experiments.cov[i] <- (p.value.cov <= alpha)
```

}

powers[j] <- mean(significant.experiments)

powers.cov[j] <- mean(significant.experiments.cov)

}

plot(possible.ns, powers, ylab="Power", xlab="Sample Size", col="black", type="l", ylim=c(0,1))

abline(h=.90, col="black", lty = "dashed")

points(possible.ns, powers.cov, type="l", ylab="Power", xlab="Sample Size", col="grey")

Daniel