rm(list=ls())
set.seed(100)

# STEP 1: Create a population
# Population size = 1000
# Half exposed, half with EM
# Qualitative effect modificaiton
# Group with no EM: True OR=2.25, Prevalence ratio = 2
# Group with EM: True OR=0.44, Prevalence ratio = 0.5

# Define population
number.unexposed.without.outcome.no.em <- 225
number.unexposed.with.outcome.no.em <- 25
number.exposed.without.outcome.no.em <- 200
number.exposed.with.outcome.no.em <- 50

number.unexposed.without.outcome.with.em <- 200
number.unexposed.with.outcome.with.em <- 50
number.exposed.without.outcome.with.em <- 225
number.exposed.with.outcome.with.em <- 25

# STEP 2: Calculate OR
true.odds.ratio.no.em <- 50/200/(25/225) 
true.odds.ratio.with.em <- 25/225/(50/200)
cat(c("True Odds Ratio without Effect Modifier = ", true.odds.ratio.no.em))
cat("\n")
cat(c("True Odds Ratio with Effect Modifier = ", true.odds.ratio.with.em))
cat("\n")

# STEP 3: Sample without replacement 10x and calculate ORs and CIs
# Create population
pop <- c(rep("T000",number.unexposed.without.outcome.no.em),rep("T010",number.unexposed.with.outcome.no.em),rep("T100",number.exposed.without.outcome.no.em),rep("T110",number.exposed.with.outcome.no.em),rep("T001",number.unexposed.without.outcome.with.em),rep("T011",number.unexposed.with.outcome.with.em),rep("T101",number.exposed.without.outcome.with.em),rep("T111",number.exposed.with.outcome.with.em))

# Calculate ORs and CIs
ors.no.em <- c()
cis.no.em <- data.frame()
ors.with.em <- c()
cis.with.em <- data.frame()
for(ii in 1:10){
	# Take a sample of 100
	sample.pop <- sample(pop,100,replace=F)
	# Estimate in those with no EM
	no.em.tab <- table(sample.pop)[c("T000","T010","T100","T110")]
	no.em.tab[which(is.na(no.em.tab))] <- 0
	ft <- fisher.test(matrix(no.em.tab,nrow=2,byrow=T))
	ors.no.em <- c(ors.no.em,ft$estimate)
	cis.no.em <- rbind(cis.no.em,ft$conf.int)
	# Estimate in those with with EM
	with.em.tab <- table(sample.pop)[c("T001","T011","T101","T111")]
	with.em.tab[which(is.na(with.em.tab))] <- 0
	ft <- fisher.test(matrix(with.em.tab,nrow=2,byrow=T))
	ors.with.em <- c(ors.with.em,ft$estimate)
	cis.with.em <- rbind(cis.with.em,ft$conf.int)
}

# Cleaning up the results format a bit: 
names(ors.no.em) <- NULL
names(ors.with.em) <- NULL
colnames(cis.no.em) <- c("Lower Confidence Limit","Upper Conifidence Limit")
colnames(cis.with.em) <- c("Lower Confidence Limit","Upper Conifidence Limit")


# STEP 4: display results
cat("\n")
cat("\n")
cat("NO EFFECT MODIFIER")
cat("\n")
cat("Sample ORs\n")
cat(round(ors.no.em,2))
cat("\n")
cat("Mean of Sample ORs\n")
cat(round(mean(ors.no.em,na.rm=T),2))
cat("\n")
cat("95% CIs\n")
print(round(cis.no.em,2))
cat("\n")
cat("\n")
cat("WITH EFFECT MODIFIER")
cat("\n")
cat("Sample ORs\n")
cat(round(ors.with.em,2))
cat("\n")
cat("Mean of Sample ORs\n")
cat(round(mean(ors.with.em,na.rm=T),2))
cat("\n")
cat("95% CIs\n")
print(round(cis.with.em,2))
cat("\n")
cat("\n")
cat("ADDITIONAL SUMMARIES")
cat("\n")
cat("Number of times that the OR for those without EM is greater than the or for those with EM: ")
cat(sum(ors.no.em>ors.with.em))
cat("\n")
cat("Number of times that the OR for those without EM is greater than 1: ")
cat(sum(ors.no.em>1))
cat("\n")
cat("Number of times that the OR for those with EM is less than 1: ")
cat(sum(ors.with.em<1))
cat("\n")