##############################################################################
#
# Purpose:
# Use pilot project data to calculate power of a full study through simulation
#
# Parts:
# (0) - Setup
# (1) - Get the pilot data and clean it
# (2) - Run the model on the pilot data and extract effects
# (3) - Set up and run the simulation
# ====> Set variables at the arrows <====
#
##############################################################################
rm(list=ls())
set.seed(424242) 

library(readr)

# (1) - Get the pilot data and clean it
#source('~/Research/tor_wikipedia_edits/handcoded_edits/inter_coder_reliability_ns0.R')
#source ('/data/users/mgaughan/kkex_data_110823_3')
data1 <- read_csv('../power_data_110923_mmt.csv',show_col_types = FALSE)
data2 <- read_csv('/data_ext/users/kcz1100/kaylea_dissertation/collabnetXS/raw_data/inst_all_packages_full_results.csv')
#d$nd <- to_logical(d$not.damaging, custom_true=c("Y")) 
#levels(d$source) <- c("IP-based Editors", "New Editors", "Registered Editors", "Tor-based Editors")

data1$up.fac.mean <- as.numeric(data2$up.fac.mean[match(data1$pkg, data2$pkg)])
# (2) - Run the model on the pilot data
pilotM <- glm(up.fac.mean ~ ((mmt) / (milestones/age)),   # give the anticipated regression a try
                  family=gaussian(link='identity'), data=data1)
summary(pilotM) #we expect effect sizes on this order

pilot.b0 <- coef(summary(pilotM))[1,1] 
pilot.b1 <- coef(summary(pilotM))[2,1]
pilot.b2 <- coef(summary(pilotM))[3,1]
pilot.b3 <- coef(summary(pilotM))[4,1] 


# (3) - Set up and run the simulation

source('powerAnalysis.R') #my little "lib"

#====>
nSims <- 5000 #how many simulations to run
n <- 100 #a guess for necessary sample size (per group)
#makeData(10) #DEBUGGING CODE -- you can uncomment this if you want to see it work
#<====

#print("Levels are:")
#print(levels(d$source))
powerCheck(n, nSims)

#Sample values
powerCheck(50, 100) 
powerCheck(80, 1000) 
powerCheck(200, 5000)