0

I have data frame named "finalend", I conducted multiple logistic regression model named "model" with all predictors as categorical or binomial variables and "TestAnxiety" as outcome variable, I was asked to calculate power, but I do not know how to do it, after reading some posts here and using copilot:

1- I think I should use simulation

2- The intuition I got is doing simulation including logistic regression model with the same data, but changing only outcome variable "TestAnxiety" for each observation randomly with each iteration based on the probability calculated from the regression model "model" I have done before, I then calculated power for each coefficient

3- I want to know if that is the correct approach and if there is anything wrong in my code

here is my code in r:

> # Set parameters
> n <- nrow(finalend)  # Sample size for each simulation
> B <- 1000  # Number of simulations
> 
> # Get the estimated coefficients from your fitted model
> betas <- coef(model)
> 
> # Initialize power
> power <- rep(0, length(betas)-1)
> 
> # Convert factor variables in finalend into dummy variables
> finalend_dummy <- model.matrix(~. - 1, data = finalend)
> 
> # Loop over simulations
> for (b in 1:B) {
+     
+     # Subset the dataframe using the names of the coefficients
+     predictors <- finalend_dummy[, names(betas)[-1]]
+     
+     # Calculate the log-odds and probability for each observation
+     logit <- betas[1] + predictors %*% betas[-1]
+     prob <- exp(logit) / (1 + exp(logit))
+     
+     # Generate random outcomes based on these probabilities
+     TestAnxiety <- rbinom(n, 1, prob)
+     
+     # Combine TestAnxiety and predictors into a data frame
+     data <- as.data.frame(cbind(TestAnxiety, predictors))
+     
+     # Fit the logistic regression model to the generated data
+     sim_model <- glm(TestAnxiety ~ ., family = binomial(), data = data)
+     
+     # Check if coefficients are significant
+     for (i in 2:length(betas)) {
+         if (summary(sim_model)$coefficients[i, 4] < 0.05) {
+             power[i-1] <- power[i-1] + 1
+         }
+     }
+ }
> 
> # Calculate power
> power <- power / B
> 
> # Print power
> print(power)
 [1] 0.083 0.057 0.044 0.204 0.073 0.244 0.101 0.101 0.064 0.267 0.114 0.109 0.144 0.601 0.859 0.100
[17] 0.286 0.147 0.106 0.082 0.998 1.000 0.259 0.329 0.059 0.701 0.479 0.295 0.995 0.094 0.770 0.370
[33] 0.653 0.259 0.134 0.600 0.082 0.307 0.877 0.351 0.460 0.855 0.855

Edit: I think calculating post hoc power is not a good idea, but instead I would calculate sample size for a follow up study, how I edit this code to calculate it?

NEA
  • 33
  • 8
    have you read about all the reasons that you should not compute post hoc/retrospective power? – Ben Bolker Mar 01 '24 at 00:20
  • @BenBolker I read about that ,but I also read it actually may be beneficial in some situations, so for sake of knowledge I want to know how to calculate it using simulation – NEA Mar 01 '24 at 00:33
  • 3
    @NEA Where have you read that this is useful? – Demetri Pananos Mar 01 '24 at 00:49
  • @DemetriPananos I read it may be beneficial if there is a follow up study (I cannot find the exact post that mentioned that), but also after some research post hoc power seems to be correlated to the p values and that's why using confidence intervals or calculated p value will give the same information we want from post hoc power, do I understand that correctly ? – NEA Mar 01 '24 at 09:23
  • 4
    How would it benefit a follow-up study? If you're planning a follow-up study, you can use the effect size from your current study as the critical effect size for your a priori power calculations for your follow-up study, but the post-hoc power doesn't tell you anything about anything outside your current sample. – Sointu Mar 01 '24 at 09:28
  • @Sointu Ok, then if I want to do a priori power calculation what should I do ? I think I would try simulation with different sample sizes until I get the desired power, also how I will use my current effect size in this simulation, would it be correct like the above code with changing only sample size ? – NEA Mar 01 '24 at 09:42
  • 1
    There are several ways. For simple regression design such as your I personally like ready packages such as WebPower. I only do simulation for complex models such as multilevel models. However, if you want to do it via simulation, this thread seems helpful. – Sointu Mar 01 '24 at 09:49
  • @Sointu Ok but my model is multiple logistic regression and I encountered WebPower and I think it cannot help me in multiple logistic regression. I also encountered this thread ,but I do not know if if I am using effect sizes correctly in simulation in the code above – NEA Mar 01 '24 at 10:01
  • 1
    Maybe post a new question with the specifics of your model? In the meantime, it might be worth it to do some reading regarding power analysis in multiple logistic regression in general. There seems to be quite a lot of resources online, including ready online calculators. – Sointu Mar 01 '24 at 13:29
  • The GPower program with its GUI allows you to compute power for a logistic-regression predictor of interest while taking into account a variety of inputs, including other predictors' combined RSQ. – rolando2 Mar 01 '24 at 14:03
  • @rolando2 Do you mean R squared by RSQ? I will try it – NEA Mar 01 '24 at 15:32
  • @Sointu ok I will try to read more about it – NEA Mar 01 '24 at 15:41
  • @NEA That's right. – rolando2 Mar 01 '24 at 21:46

0 Answers0