#Introduction to Simulation/Monte Carlo
#Studies Using R

#1) Goal/Research Question
#Can the results of omnibus tests and
#pairwise comparisons produce different
#conclusions regarding whether or not
#there are any differences among the
#populations being studied?

#2) Data Generation

#2a) Constants
#number of simulations
nsim<-100 
#nominal Type I error rate
alpha<-.05 
#number of groups
g<-5
#group sample sizes 
N <- rep(10, g)
#population sds, specified by case 
#Farmus method
sd <- rep(rep(3, g), times = N) 
#population means, specified by case 
mu <- rep(c(0, 0, 0, 0, 2), times = N) #population mean for each case
#specification of the independent variable
IV=factor(rep(LETTERS[1:g], times = N))

#2b Variables - none 


#Empty outcome matrix
results<-matrix(NA,nrow=nsim,ncol=3)

#2c Run the simulations via a "for" loop
for (i in 1:nsim) {
  DV=rnorm(sum(N),mu,sd)
  ifelse ((oneway.test(DV ~ IV)$p.value<=alpha), 
          results[i,1]<-1, results[i,1]<-0)
  pairwh<-pairwise.t.test(DV,IV,p.adjust="none")$p.value
  pairwh<-pairwh[!upper.tri(pairwh)]
  ifelse (any(pairwh<=alpha),
          results[i,2]<-1, results[i,2]<-0)
  ifelse (results[i,1]==0 & results[i,2]==1,
          results[i,3]<-1, results[i,3]<-0)
}

#Summarize the Results
#Mean rate over all simulations for each column
sum_results<-colMeans(results)
sum_results<-data.frame(t(sum_results))
names(sum_results)<-c("omnibus","any_pairw","missed_pairw")
sum_results


#We would probably want to manipulate some variables,
#not have everything as a constant

#1) Goal/Research Question
#Can the results of omnibus tests and
#pairwise comparisons produce different
#conclusions regarding whether or not
#there are any differences among the
#populations being studied?

#2) Data Generation

#2a) Constants
#number of simulations
nsim<-100 
#nominal Type I error rate
alpha<-.05 
#number of groups
g<- 5
#number of conditions
ncond<-2
#group sample sizes
N <- matrix(rep(10,g*2),
            nrow=ncond,ncol=5,byrow=TRUE)
N
 

#population sds, specified by case 
#Farmus method
sd <- matrix(rep(rep(3,g*ncond), times = N),
            nrow=2,byrow=TRUE) 
sd


#2b Variables 
#population means, specified by case 
popmn<-c(0,0,0,0,2,0,0,1,2,2)
mu <- matrix(rep(popmn, times = N),
             nrow=ncond, byrow=TRUE)
mu

#Empty outcome matrix
results<-matrix(NA,nrow=nsim*ncond,ncol=3)

#2c Run the simulations via a "for" loop
k<-0
for (i in 1:ncond) {
  for (j in 1:nsim) {
    k<-k+1
    IV=factor(rep(LETTERS[1:g], times = N[i,]))
    DV=rnorm(sum(N[i,]),mu[i,],sd[i,])
    ifelse ((oneway.test(DV ~ IV)$p.value<=alpha), 
          results[k,1]<-1, results[k,1]<-0)
    pairwh<-pairwise.t.test(DV,IV,p.adjust="none")$p.value
    pairwh<-pairwh[!upper.tri(pairwh)]
    ifelse (any(pairwh<=alpha),
          results[k,2]<-1, results[k,2]<-0)
    ifelse (results[k,1]==0 & results[k,2]==1,
          results[k,3]<-1, results[k,3]<-0)
  }
}

#Summarize the Results

#Empty matrix to put the results, including columns
#for the population means since that is a variable
sumresults<-matrix(NA,nrow=2,ncol=g+ncol(results))
#Put the population means into the first g columns
sumresults[1,1:g]<-popmn[1:g]
sumresults[2,1:g]<-popmn[(g+1):(g*2)]
#Add the column means for the results to the last 3 columns
sumresults[1,(g+1):(g+3)]<-colMeans(results[1:nsim,])
sumresults[2,(g+1):(g+3)]<-colMeans(results[(nsim+1):(nsim*2),])
sumresults<-data.frame(sumresults)
names(sumresults)<-c(rep("mu",g),"omnibus","any_pairw","missed_pairw")
sumresults