#############################################################################
#####   Chi-square vs. logistic regression
#############################################################################
##### Do this simulation with them


rm(list=ls())

n<-100
x<-rbinom(n,1,0.5)
py<-ifelse(x==1,0.7,0.5)
y<-rbinom(n,1,py)
table(x,y)

chisq.test(x,y)

mod<-glm(y~x,family="binomial")
summary(mod)
summary(mod)$coeff[2,4]

pval.chisq<-pval.chisq2<-pval.logis<-NA
for (i in 1:1000){
  x<-rbinom(n,1,0.5)
  py<-ifelse(x==1,0.7,0.5)
  y<-rbinom(n,1,py)
  table(x,y)
  
  pval.chisq[i]<-chisq.test(x,y)$p.value
  pval.chisq2[i]<-chisq.test(x,y,correct=FALSE)$p.value
  
  mod<-glm(y~x, family="binomial")
  summary(mod)
  pval.logis[i]<-summary(mod)$coeff[2,4]
}

cor(pval.chisq,pval.logis)
plot(pval.chisq,pval.logis)
abline(c(0,1))

cor(pval.chisq2,pval.logis)
plot(pval.chisq2,pval.logis)
abline(c(0,1))



##########################################################################
#####  Simulate Power and Sample Size
##########################################################################
##### Have them do this themselves
#####  1. What is power to detect a difference (2-sided p-value <0.05) 
#####       between p1=0.7 and p0=0.5, with approximately equal sized groups
#####       and n=100?
#####  2. What is power if there are twice as many with x=1 compared to x=0?
#####  3. What sample size is needed to have 80% power with equal sized groups,
#####       2-sided alpha=0.05?
#####  4. With a sample size of n=100, equal sized groups, p0=0.5, what value
#####       of p1 will you have 80% power to detect with a 2-sided alpha=0.05?
#####       What is the detectable odds ratio in this setting?
##############################################################################


rm(list=ls())

n<-100
p<-0.5
p1<-0.7
p0<-0.5

pval.logis<-NA
for (i in 1:10000){
  x<-rbinom(n,1,p)
  py<-ifelse(x==1,p1,p0)
  y<-rbinom(n,1,py)
  
  mod<-glm(y~x, family="binomial")
  summary(mod)
  pval.logis[i]<-summary(mod)$coeff[2,4]
}
mean(pval.logis<0.05)




############################################################################
#####      Power and sample size calculations using functions
############################################################################

## install.packages("pwr")
library(pwr)

p0<-0.5
p1<-0.7
n1<-n0<-50
pwr.2p2n.test(n1=50, n2=50, ES.h(p0, p1), sig.level=0.05, alternative="two.sided")

pwr.2p2n.test(n1=33, n2=67, ES.h(p0, p1), sig.level=0.05, alternative="two.sided")

pwr.2p.test(power=0.8, ES.h(p0, p1), sig.level=0.05, alternative="two.sided")
##### Notice that the function changed a little

pwr.2p2n.test(n1=50, n2=50, power=0.8, sig.level=0.05, alternative="two.sided")

?ES.h
# 2*asin(sqrt(0.5))-2*asin(sqrt(p1)) = h
h<-0.5603 # (plus or minus)
(sin((h+2*asin(sqrt(0.5)))/2))^2




####  Power for continuous outcome

library(pwr)
pwr.t2n.test(n1=21, n2=24, d=NULL, sig.level=.05, power=.8, alternative="two.sided")
pwr.t2n.test(n1=21, n2=31, d=NULL, sig.level=.05, power=.8, alternative="two.sided")
pwr.t2n.test(n1=24, n2=31, d=NULL, sig.level=.05, power=.8, alternative="two.sided")