library(rms)
#Power simulation for Yi = alpha + beta*Ai + Ei


fnPower <- function(intercept, beta, sd,
		nTotalSubjects, nTreated,
		nIterations, alpha=0.05){
  start.time <- Sys.time()
  pVec <- betaVec <- seVec <- rep(NA,nIterations)
  tx <- c(rep(1,nTreated),rep(0,nTotalSubjects - nTreated))
  for(i in 1:nIterations){
    resid <- rep(rnorm(nTotalSubjects,0,sd))
    y <- intercept + beta*tx + resid
    o <- ols(y~tx,x=TRUE,y=TRUE)
    v <- robcov(fit=o, cluster=1:length(tx))
    pVec[i] <- 2*pnorm(-abs(v$coeff[2]/sqrt(diag(v$var))[2]))
    betaVec[i] <- v$coef[2]
    seVec[i] <- sqrt(diag(v$var))[2]
  }
  power <- length(pVec[pVec<alpha])/length(pVec)
  return(list(power=power,p=pVec,beta=betaVec,se=seVec))
}

# Should return about alpha
outCalibrate <-
fnPower(intercept=-0.25,beta=0,sd=1,
	nTotalSubjects=100,nTreated=50,
	nIterations=10000, alpha=0.1)
outCalibrate$power
hist(outCalibrate$p)



library(pwr)
pwr.t.test(500, 0.15)

outCalibrate <-
fnPower(intercept=-0.25,beta=0.15,sd=1,
  nTotalSubjects=1000,nTreated=500,
  nIterations=10000, alpha=0.05)
outCalibrate$power


# Power curve for sample size
sizes <- seq(100, 3000, by=50)
powerVals <- vector(length=length(sizes))
for (i in 1:length(sizes)) {
  powerVals[i] <- fnPower(intercept=1,beta=0.15,sd=1,
	nTotalSubjects=sizes[i],nTreated=sizes[i]/2,
	nIterations=100)$power
}
powerCalc <- data.frame(sizes, powerVals)
ggplot(powerCalc) + aes(x=sizes, y=powerVals) + geom_smooth()