24_deb_pkg_gov/R/powerAnalysis.R

# This is semi-generic code for doing a power analysis of a logistic regression with 4
# levels in a factor
# when there's some pilot values already available and defined
#modelled heavily the simulation example explained in:
#http://meeting.spsp.org/2016/sites/default/files/Lane%2C%20Hennes%2C%20West%20SPSP%20Power%20Workshop%202016.pdf

library('batman')
library('reshape')

l2p <- function(b) {
 odds <- exp(b) 
 prob <- odds/(1+odds)
 return(prob)
}


#Matt: 
makeDataNew <- function(n) {
  sDF <- data.frame(
    sample(2793.638:3066.417, 1, replace = FALSE),
    sample(27.6519:154.6866, 1, replace = FALSE),
    sample(50.01884:96.77090, 1, replace = FALSE),
    sample(c(1,2), 1, prob=(c(0.5,0.5)), replace = FALSE),
    sample(1.610638: 1.684438, 1, replace = FALSE),
    sample(-0.1961401:-0.1647757, 1, replace = FALSE)
  )
  colnames(sDF) <- c('age', 'contributors', 'collaborators', 'milestones', 'mmt', 'up.fac.mean')
  return(sDF)
}

powerCheck <- function(n, nSims) { #run a power calculation on the dataset given
  #set up some empty arrays b/c R
  signif0 <- rep(NA, nSims)
  signif1 <- rep(NA, nSims)
  signif2 <- rep(NA, nSims)
  signif3 <- rep(NA, nSims)
  signifM <- rep(NA, nSims)
  for (s in 1:nSims) {           # repeatedly we will....
    simData <- makeDataNew(n)       # make some data
    #have updated for kkex through here, now need to look at the underproduction work
    m1.sim <- glm(up.fac.mean ~ ((mmt)/ (milestones/age)),   # give the anticipated regression a try
                  family=gaussian(link='identity'), data=simData)
    p0 <- coef(summary(m1.sim))[1,4]
    p1 <- coef(summary(m1.sim))[1,4]
    p2 <- coef(summary(m1.sim))[1,4]
    p3 <- coef(summary(m1.sim))[1,4]
    signif0[s] <- p0 <=.05
    signif1[s] <- p1 <=.05
    signif2[s] <- p2 <=.05
    signif3[s] <- p3 <=.05
    signifM[s] <- p0 <=.05 & p1 <=.05 & p2 <=.05 & p3 <=.05
  }  
  power <- c(mean(signif0), mean(signif1), mean(signif2), mean(signif3), mean(signifM))
  return(power)
}