1
$\begingroup$

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
$\endgroup$
6
  • $\begingroup$ If you have 10 years (as in your example), are you interested in treating year as a continuous covariate (i.e. testing for linear [on the log-odds scale] trends over time), or in estimating whether there is an overall significant difference among years (treating years as either a categorical fixed effect or as a random-effects grouping variable, i.e. year-as-cluster)? (I can't say I'm super-enthusiastic about cleaning up after ChatGPT, although I appreciate that you made the effort ...) $\endgroup$
    – Ben Bolker
    Commented Jul 10 at 21:43
  • $\begingroup$ The general idea of simulating data with a specific effect size and seeing how often you can reject the null hypothesis is a good idea, though. The combination of simulate_new() + glmmTMB() + a for-loop (or replicate or map or whatever) should do it ... the t-test and 'power of a proportion test' stuff seem to be red herrings. $\endgroup$
    – Ben Bolker
    Commented Jul 10 at 21:45
  • $\begingroup$ @BenBolker, sorry, I wasn't keen on using ChatGPT, but had to start somewhere. I wasn't sure what might be easier (detecting a change in slope between year 1 and 10, or comparing means of two groups/factors of 5 years worth of data; i.e. a sort of "before" and "after" comparison). Either works more me. Thank you for looking it over. $\endgroup$
    – Nate
    Commented Jul 10 at 22:36
  • $\begingroup$ I started to write an answer but realized I was making it very complicated. I think you should figure out first what kind of test you want to run (year 1 vs year 10? average of year 1-5 vs year 6-10? slope of trend from year 1 to year 10?) and what kind of random effects you think are reasonable/would put in the model. Then it will be relatively easy to show you how to generate simulations/fit models/test effects/summarize effects to compute the power for a specified set of assumptions. But I think "how should I test this?" is a separate question from "how do I calc power of this test?" $\endgroup$
    – Ben Bolker
    Commented Jul 11 at 0:22
  • 2
    $\begingroup$ The answers to your questions about "how do I decide what parameters to use?" are that you have to (1) understand what the parameters mean and (2) play around with options until you get results that look like they would be realistic for your context. $\endgroup$
    – Ben Bolker
    Commented Jul 11 at 0:23

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.