beta_binomial_power.R

pwr <- function(S, a, mu = 0.086, phi=240){
  D <- S * (S - 1)
  2 * (D + 1) * (D/(phi + 1) + 1) * mu * ( 1 - mu) + ( a * mu + mu)  * (1 - a * mu - mu) * (2.8)^2 / ((D + 1)*a*mu)^2
}

pwr(6, 0.10)
curve(pwr(x, 0.4), xlim = c(3,30))

curve(pwr(x, 0.2), xlim = c(3,30), log = "y")


curve(pwr(x, 0.2), xlim = c(3,30), log = "y")

library(tidyverse)

tidyr::expand_grid(S = 3:30, p=c(.05, .1, .2, .4, .5, .7)) %>% 
  mutate(N = pwr(S, p)) %>% 
  ggplot(aes(x = S, y = N, colour = p, group = p)) + 
  geom_line(lwd = 2) + 
  coord_trans(y = "log") +
  scale_color_distiller(palette = "Greens") + 
  labs(x = "Species richness",
       y = "Sample size",
       colour = "effect size %") + 
  theme_dark() +
  scale_y_log10(
  breaks = scales::trans_breaks("log10", function(x) 10^x),
  labels = scales::trans_format("log10", scales::math_format(10^.x))
)+ 
  geom_hline(yintercept = 30, lty = 2)