Low power = inaccurate effect size estimates

effect_sizes.R

#------------------------------------------------------------------------------
rm(list = ls())

set.seed(50)

# function for generating random draws from multivariate distribution
# n = number of draws; p = number of variables
# u = mean of each variable; s = SD of each variable
# corMat = correlation matrix
mvrnorm <- function(n, p, u, s, corMat) {
  Z <- matrix(rnorm(n * p), p, n)
  t(u + s * t(chol(corMat)) %*% Z)
}

#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
### simulation setup

# declare simulation parameters
means <- c(100, 600)
sds <- c(20, 80)
n_sims <- 1000

# declare correlation coefficient & generate cor. matrix
cor <- 0.3
cor_mat <- matrix(c(1, cor, 
                    cor, 1), nrow = 2, ncol = 2, byrow = TRUE)

# sample sizes to simulate
sample_sizes <- c(10, 20, 30, 50, 85, 170, 500, 1000)

# number of simulations

# create variable to store data in
final_data <- matrix(nrow = n_sims, ncol = length(sample_sizes))
colnames(final_data) <- sample_sizes
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
### simulation execution

for(i in 1:length(sample_sizes)){
  for(j in 1:n_sims){

    # get the experiment data
    sim_data <- mvrnorm(sample_sizes[i], p = 2, u = means, s = sds, 
                        corMat = cor_mat)
    
    # perform the correlation
    sim_cor <- cor.test(sim_data[, 1], sim_data[, 2], method = "pearson")
    
    # if the correlation is significant, store the effect size
    if(sim_cor$p.value < 0.05){
      final_data[j, i] <- as.numeric(round(abs(sim_cor$estimate), 3))
    }
  }
}
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
## draw the plot
boxplot(final_data, na.action=na.omit, xlab = "Simulated Sample Size", 
        ylab = "Measured Effect Size")
abline(h = cor, col = "red", lty = 2, lwd = 3)

legend("topright", "True Effect Size", cex = 1.5,  lty = 2, bty="n", 
       col = "red", lwd = 3)
#------------------------------------------------------------------------------

Just a head's up, n_sims is not included in the code. I just entered n_sims <- 1000 and the code works fine, but you might want to add it in so it works right away.

Thanks for this nice coding piece. Just echo what is said by BartlettJE.

Sorry for delay on this - Now done :)