Hi all,

I did phylogenetic comparative analysis of genome size and environmental factors using more than 120 species. The function I used is from the package phylolm:

mod <- phylolm(log_GS ~ factor_scaled, data = df, phy = target_tree, model = "BM", boot = 10000)

The result is non-significant. In order to better understand the required effect size to detect with my current sample size, I've been exploring power analysis using simulation. I found a code snippet for power analysis, which is quite straightforward. In the following code, I assume that Y is genome size, and X is environmental factor in my case.

my_power_function <- function(effect, sample_size) {
  sig_results <- c()
  for (i in 1:500) {
    X <- runif(sample_size, 0, 1)
    Y <- effect * X + rnorm(sample_size, 0, 3)
    tib <- data.frame(X = X, Y = Y)
    model <- feols(Y ~ X, data = tib, se = 'hetero')
    sig_results[i] <- tidy(model)$p.value[2] <= 0.05
  }
  sig_results %>% mean() %>% return()
}

power_levels <- c()
effects_to_try <- c(.4, .8, 1.2, 1.6, 2)

for (i in 1:5) {
  power_levels[i] <- my_power_function(effects_to_try[i], 120)
}

power_results <- data.frame(effect = effects_to_try,
                                               power = power_levels)

However, I am facing some challenges in generating the simulated dataset. The environmental factor in my study are not uniformly or normally distributed, and I am also uncertain about how to incorporate the residual variance with the phylogenetic signal.

If anyone has any ideas or clues on how I can modify the code or any other approaches to achieve my goal, I would greatly appreciate your input. Thanks in advance!