05_power_simulation.R

   1 #!/usr/bin/env Rscript
   2
   3 # Perform power analysis to assess whether we have enough data to study bots
   4
   5 # Copyright (C) 2018  Nathan TeBlunthuis
   6
   7 # This program is free software: you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation, either version 3 of the License, or
  10 # (at your option) any later version.
  11
  12 # This program is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16
  17 # You should have received a copy of the GNU General Public License
  18 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
  19
  20 source("lib-00-utils.R")
  21 library(effects)
  22 library(texreg)
  23 if(!exists("r")){
  24     source("lib-00-utils.R")
  25     source("01_build_datasets.R")
  26 }
  27
  28 p.outcome <- r$newcomer.summary.stats$p.survives
  29 p.dv <- r$newcomer.summary.stats$p.bot.reverted
  30 n <- r$halfak.model@gof[5]
  31
  32 sample.ds <- function(n,p.outcome,p.dv,eff = -0.01){
  33     dv <- rbinom(n=n,size=1,prob=p.dv)
  34     iv <- rbinom(n,size=1,prob=p.outcome)
  35     m1 <- glm(iv ~ 1, family=binomial(link='logit'))
  36     eta <- eff*dv + coef(m1)[1]
  37     p <- exp(eta)/(1+exp(eta))
  38     tmp <- runif(n)
  39     y <- (tmp < p)
  40     fit <- glm(y ~ dv,family=binomial(link='logit'))
  41     summary(fit)$coefficients[2,4]
  42 }
  43
  44 eff <- -0.68
  45 remember(exp(-eff),"power.analysis.effect")
  46 pwr.test.sig.level <- 0.05
  47 remember(pwr.test.sig.level)
  48 n.power.sim <- 1000
  49 remember(n.power.sim)
  50
  51 out <- replicate(n.power.sim,sample.ds(n,p.outcome,p.dv,eff=eff))
  52 remember(mean(out<pwr.test.sig.level),"pwr.test")
  53