Basic power analysis feature?

[DominiqueMakowski]

Disclaimer:

• I personally don't like power analyses and for multiple reasons I don't think it should be a prerequisite for good science
• I know "best practices" in this context are controversial
• I know there's already a lot of tools to do that (in particular very cool shiny apps)

And yet most papers still use G*Power and basic tools rather than complicated simulation-based approaches.

I wouldn't mind a simple and intuitive function in effectsize to do some basic stuff, like:

# Correlation and differences
minimum_sample(min_r=0.1)
minimum_sample(min_d=0.1, alpha=..., beta=...)

This would return the minimum number of observations to detect that effect. We could add some cool power analysis curve plots in see. What do you think?

[bwiernik]

Sure. With our bootstrap CI features, we can also do simulation based power/precision analyses pretty easily

[strengejacke]

We should also look at @debruine's faux package https://debruine.github.io/faux/) and powerlmm.

[mattansb]

The pwr package is a great replacement for G*power - perhaps we can write a vignette or a blog post (lol we have a blog???) about using pwr with effectsize? I have some examples for this from one of my courses (+ an example in the Fei paper!).

With our bootstrap CI features, we can also do simulation based power/precision analyses pretty easily

Um, we don't have those 😅 _{(yet? #578)}

We should also look at @debruine's faux package https://debruine.github.io/faux/) and powerlmm.

Yup - faux, powelmm, Superpower, and simr are all great for mixed-like models. I do think though that anything more complex than a one-way ANOVA / linear model, is way wayyyyyy out of the scope of anything we could properly support...

[DominiqueMakowski]

I do think though that anything more complex than a one-way ANOVA / linear model, is way wayyyyyy out of the scope of anything we could properly support...

Fully agree, even mixed models are borderline IMO - Starting with basic h-tests would already be a useful feature

[bwiernik]

Um, we don't have those 😅

Guess I was thinking of parameters?

[pdwaggoner]

Hi all - its been a while, but I have been thinking about this a little. Here's an idea. Is this along the lines of what you're thinking @DominiqueMakowski ?

# power analysis idea

power <- function(n = 500,  nsim = 100, 
                  alpha = 0.05, beta = 0.2,
                  seed = 1234) {
  
  set.seed(seed)
  
  result <- c()
  
  for (i in 1:nsim) {
    y_0 <- rnorm(n = n, mean = 0, sd = 1)
    y_1 <- y_0 + beta 
    assignment <- rbinom(n = n, size = 1, prob = .5) 
    outcomes <- (y_0 * (1 - assignment)) + (y_1 * assignment)
    p <- t.test(outcomes ~ assignment)$p.value 
    sig <- dplyr::if_else(p <= alpha, 1, 0)
    result <- c(result, sig)
}

  return(
    tibble::tibble(
      power = mean(result),
      full = list(result)
      )
    )
  
}

# result with default vals
power()

# A tibble: 1 × 2
power   full       
<dbl>   <list>     
0.6   <dbl [100]>

[pdwaggoner]

Hi all - I cleaned up the code a bit. Here's a newer and faster version. Let me know thoughts and I can starting working on it and open a pull shortly, only if it would add value to the ES ecosystem:

power <- function(n = 500,  nsim = 100, 
                  alpha = 0.05, beta = 0.2,
                  seed = 1234) {
  
  set.seed(seed)
  
  res <- tibble::tibble(sim = 1:nsim) %>%
    mutate(
      exp_res = 
        purrr::map_dbl(
          sim, ~{
            y_0 <- rnorm(n = n, mean = 0, sd = 1)
            y_1 <- y_0 + beta
            a <- rbinom(n = N, size = 1, prob = 0.5)
            o <- (y_0 * (1 - a)) + (y_1 * a)
            tt <- t.test(o ~ a)
            p <- tt$p.value
            ci <- tt$conf.int # keep? / expand?
            s <- as.numeric(p <= alpha)

            return(mean(s))
            }
          )
    )
  
  return(
    tibble::tibble(
      power = mean(res$exp_res),
      full = list(res$exp_res)
      )
    )
}

[pdwaggoner]

ping @DominiqueMakowski @bwiernik @mattansb @strengejacke

[mattansb]

To reiterate my previous point, instead of having these functions baked in (which can only cover simple cases, already covered by pwr), I think it would be much more useful to users to have a vignette about this topic. It can cover how to use the pwr package, or how to run simple simulations using some of the common packages and how effectsize might be used alongside them (can it?).

[DominiqueMakowski]

As I am not a user / expert in this type of stuff, I don't have a strong opinion regarding its implementation.

In any case, a vignette would definitely be useful. Maybe it can be a starting point, and during the development of the vignette we will see whether there is a need for an easystatsification (adding features to help, support or report existing tools).

[pdwaggoner]

Many thanks to you both. I agree and am fine with that approach. I am already working on the vignette. I'll share via PR when I have something worth sharing, at which point others can review as they'd like

[strengejacke]

Here's (an old) script I used to teach how power analysis and simulations are similar. It contains German comments (maybe Deepl can help you), I probably can translate into English:

https://gist.github.com/strengejacke/1eaa54ea1fc235b0523373c0245912da

It shows a simple example for t-tests, and later how to calculate the power for model parameters (and compares the results to an Anova table that includes power per parameters).

Based on this, you can go on to more complex scenarios for power calculations.

[pdwaggoner]

Thanks @strengejacke - super helpful! This is very similar to my approach, which I am revising now. Stand by for more.

[pdwaggoner]

Alright, just opened PR #605 (and closed erroneous #604 ; apologies!), and requested @mattansb as reviewer. Let me know what else is needed/next. Thanks!

[DominiqueMakowski]

I have a memory of someone describing "fireplots" or something like that, and in my memory it looked kinda like

But when I google fireplot power analysis nothing shows up (but that fig above). Does that ring a bell to anyone?

[pdwaggoner]

@DominiqueMakowski - I have never heard the term "fireplot" before. but if you mean a heatmap, you could create something like this:

library(tidyverse)

# Create a grid of sample sizes and effect sizes
sample_sizes <- seq(10, 100, by = 10)
effect_sizes <- seq(0.1, 2, by = 0.1)

# Generate data for the heatmap
data <- expand.grid(SampleSize = sample_sizes, EffectSize = effect_sizes)

# Calculate power for each combination of sample size and effect size
data <- data |>
  mutate(Power = pwr.t.test(n = SampleSize, d = EffectSize, sig.level = 0.05, type = "two.sample")$power)

# viz
ggplot(data, aes(x = SampleSize, y = EffectSize, fill = Power)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "blue") +
  labs(x = "Sample Size", y = "Effect Size", fill = "Power") +
  theme_minimal()

which gives:

Is this what you mean?

[DominiqueMakowski]

I... don't know if it is what I meant (i.e., if this is was I saw somewhere on twitter) but in any case it looks somewhat nice.

Haha! I found the tweet I was thinking of: https://x.com/dsquintana/status/1641339624231501824?s=20 but it seems like it's mostly (?) for meta-analyses

That said, I think this type of plot is my fav:

Also, here are some other vignettes for inspiration:

• https://www.andywills.info/rminr/power.html
• https://www.andywills.info/rminr/effsize_from_papers.html
• Bayesian power analysis (we should definitely mention that): https://www.andywills.info/rminr/power-bayesian.html

[pdwaggoner]

I am pretty pumped I created that one before I saw your tweet example 😃 I am happy to pull something together to expand the viz capability for power analysis in the package, though it seems like @mattansb doesn't want to add too much re: power analysis to effectsize though given its inclusion in other packages, which I understand.

[mattansb]

We now have a power vignette, that can be expanded on to more examples of using pwr+effectsize (or something else also?).

Thanks @pdwaggoner !

---

I am happy to pull something together to expand the viz capability for power analysis in the package

I think a good place for this would actually be to contribute to the pwr package, that currently does have a plotting method, but it leaves much to be desired...

out <- pwr::pwr.t.test(n = 100, d = 0.3)

plot(out)

Alternatively, if we want that capability in-house, it could be implemented in {see}!

[DominiqueMakowski]

ontribute to the pwr package

Or we add a report() wrapper and a see::plot() method to it