R: How to calculate power to detect a change in mean beta reg.?

Question

I'm interested in calculating the power to detect a change in my mean response variable by a certain effect size, but I'm not sure how to start. My Y variable is mean plant cover 0-100% cover (mean of 10 random replicate surveys per site, as a decimal, including 0 and 1), and I'd like to calculate the power to detect an effect size for a future hypothetical change in cover (one/paired-sample test). Plant cover is repeatedly sampled every year at the same locations (random intercept).

I typed this question into chatgpt, but I'd like to verify this is the correct method. It responded with suggesting I make two separate models (one in which my covariate has no effect, and another where it has a specified effect), then calculate the difference between the two mean responses from each model. It then says to calculate power using a t-test.

Does this answer seem correct? How would I accurately simulate data that looks like my raw data (how do I choose alpha and beta in rbeta)? How would I make my year covariate have a specific effect on cover (chatgpt assumed a continuous covariate and set it to 0.5)? Also, I know my response isn't a probability, but shouldn't I use a power calculation for a proportion as this is 0->1 data (i.e., use pwr::pwr.p.test)?

I tried simulating the data to make it easier to post this question, but it didn't have the intended effect in the summary output.

set.seed(123)

n <- 100
x1 <- rep(seq(1, 10), each = 10)
x2 <- seq(1, 10)

y <- rbeta(n, 2, 5)  # Simulate response variable (beta distributed)

# Create data frame
data <- data.frame(cover = y, year = x1, site = x2)
data$site <- as.factor(data$site)

I thought this would be better than betareg as it allows for random intercepts

library(glmmTMB)
# Fit beta regression model
model <- glmmTMB(y ~ year + (1|site), ziformula = ~1, data = data, family=beta_family(link="logit"))

Chatgpt:

# Simulate Data for Power Calculation:
# To calculate power, you typically need to simulate data under the null hypothesis and alternative hypothesis.

# Null Hypothesis (H0): Assume the predictor variable x1 has no effect.
# Alternative Hypothesis (H1): Assume the predictor variable x1 has an effect of a certain size.
# Let's assume you want to calculate power for detecting a difference in the mean response variable y between two values of x1: x1 = 0 and x1 = 0.5.

# Simulate data under H0 (no effect of x1)
data_null <- data.frame(y = rbeta(n, 2, 5), x1 = rnorm(n))

# Simulate data under H1 (effect of x1 = 0.5)
data_alt <- data.frame(y = rbeta(n, 2, 5), x1 = rnorm(n, mean = 0.5))



model_null <- betareg(y ~ x1, data = data_null)
model_alt <- betareg(y ~ x1, data = data_alt)

# Calculate difference in mean response variable (y) between x1 = 0.5 and x1 = 0 under H0
mean_y_alt <- predict(model_alt, newdata = data_null)
mean_y_null <- predict(model_null, newdata = data_null)
effect_size <- mean(mean_y_alt - mean_y_null)

# Standard deviation of y
sd_y <- sd(predict(model, type = "response"))

# Sample size
n <- nrow(data)

# Significance level
alpha <- 0.05

# Calculate power
power <- pwr.t.test(n = n, delta = effect_size, sd = sd_y, sig.level = alpha, type = "two.sample")$power
power

This is a better test for proportion data, no?

power <- pwr.p.test(h = effect_size, n = n, sig.level = alpha, alternative = "two.sided")$power
power

Answer 1

I think that ChatGPT has misled you somewhat (not surprisingly).

I tried simulating the data to make it easier to post this question, but it didn't have the intended effect in the summary output.

That is a great idea. I am very much a fan of doing power calculations via simulation. It is not always easy, but by doing it that way you can control whatever you want to control. In other words it is very flexible. So, the main problem with your simulation lies here:

y <- rbeta(n, 2, 5)  # Simulate response variable (beta distributed)

Here you have simulated y but it does not depend on year or site. In order for this to work we need to incorporate these dependencies. Luckily it is not too hard to do, so I have adapted your code. It computes the power to detect an effect size of 0.08 for year. Hopefully it is self-explanatory, but if not, feel free to ask (or point out any deficiencies).

library(glmmTMB)

set.seed(15)
n_sites <- 10                   # Number of sites 
n_years <- 5                    # Number of years
n_rep <- 5                      # Number of replicate surveys per site per year

n <- n_sites * n_years * n_rep  # Total number of observations

year <- rep(rep(1:n_years, each = n_rep), n_sites) 
site <- rep(1:n_sites, each = n_years * n_rep)     

n_sim <- 50                   # Number of simulations
vec_pvalues <- numeric(n_sim) # Vector to store the p-values
eff_size <- 0.08

for (i in 1:n_sim) {
  
  # Simulate random intercepts for `site`
  random_effects <- rnorm(n_sites, 0, 0.5)  
  
  # Linear predictor (year effect and random effect for site)
  eta <- eff_size * year + random_effects[site]
  
  # Inverse logit transformation to bound the mean between 0 and 1
  mu <- plogis(eta)
 
  # Precision parameter for the beta distribution
  phi <- 10  

  # Simulate beta-distributed response variable y
  y <- rbeta(n, mu * phi, (1 - mu) * phi)
  
  data <- data.frame(cover = y, year = year, site = site)
  data$site <- as.factor(data$site)
  
  # Fit beta regression model
  model <- glmmTMB(cover ~ year + (1|site), ziformula = ~1, data = data, family=beta_family(link="logit"))
  
  # Extract p-value for the effect of `year`
  vec_pvalues[i] <- coef(summary(model))$cond[2, 4] 
}
# Calculate power
mean(na.omit(vec_pvalues) < 0.05)

That results in 0.82 in this instance.

Ideally there should be some code to check for convergence issues. In this instance there were none, so I didn't bother with that.