#!/usr/bin/env Rscript

# Perform power analysis to assess whether we have enough data to study bots

# Copyright (C) 2018  Nathan TeBlunthuis

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.

source("lib-00-utils.R")
library(effects)
library(texreg)
if(!exists("r")){
    source("lib-00-utils.R")
    source("01_build_datasets.R")
}

p.outcome <- r$newcomer.summary.stats$p.survives
p.dv <- r$newcomer.summary.stats$p.bot.reverted
n <- r$halfak.model@gof[5]

sample.ds <- function(n,p.outcome,p.dv,eff = -0.01){
    dv <- rbinom(n=n,size=1,prob=p.dv)
    iv <- rbinom(n,size=1,prob=p.outcome)
    m1 <- glm(iv ~ 1, family=binomial(link='logit'))
    eta <- eff*dv + coef(m1)[1]
    p <- exp(eta)/(1+exp(eta))
    tmp <- runif(n)
    y <- (tmp < p)
    fit <- glm(y ~ dv,family=binomial(link='logit'))
    summary(fit)$coefficients[2,4]
}

eff <- -0.68
remember(exp(-eff),"power.analysis.effect")
pwr.test.sig.level <- 0.05
remember(pwr.test.sig.level)
n.power.sim <- 1000
remember(n.power.sim)

out <- replicate(n.power.sim,sample.ds(n,p.outcome,p.dv,eff=eff))
remember(mean(out<pwr.test.sig.level),"pwr.test")