library(tidyverse)
library(nlme)
library(lme4)

num_subjects <- 43
num_iterations <- 1000
num_trials <- 100
exp.pre <- read.csv("experiment3.csv")

# Run Monte Carlo simulation to calculate power and coefficient confidence intervals
run_simulation <- function(data) {
  # Pick random subjects
  subjects <- sample(unique(data$subject), num_subjects, replace=TRUE)
  
  # Pull rows amd add to new data set
  sampledata <- data[0:0,]
  i <- 1
  for (samplesubject in subjects) {
    sampletrials <- sample((data[data$subject == samplesubject,]$trial), num_trials, replace=TRUE)
    for (eachnumber in sampletrials) {
      row <- data[(data$subject == samplesubject) & (data$trial == eachnumber),]
      row$subject <- i
      sampledata <- rbind(sampledata, row)
    }
    i <- i + 1
  }
  
  # Fit a full and null model and test for significance
  M1.full <- lme(rt ~ hand * response + response * space, random = ~1 | subject, method = "ML", data = sampledata)
  M1.null <- lme(rt ~ response * space + hand, random = ~1 | subject, method = "ML", data = sampledata)
  comparison <- anova(M1.full, M1.null)
  p.val <- comparison$`p-value`[2]
  
  # Fit a full model with REML for the values of coefficients
  model <- lme(rt ~ hand * response + response * space, random = ~1 | subject, method = "REML", data = sampledata)
  
  # Return the model and p-value of the anova
  return(list("model"=model, "p"=p.val))
}

# Run the simulations
sim_results <- list()
for (i in 1:num_iterations) {
  sim_results[[i]] <- run_simulation(exp.pre)
  print(i)
}

p <- rep(0, num_iterations)
for (i in 1:length(sim_results)) {
  p[i] <- sim_results[[i]]$p
}
sum(p < 0.05) / num_iterations

coef <- rep(0, num_iterations)
for (i in 1:length(sim_results)) {
  coef[i] <- sim_results[[i]]$model$coefficients$fixed[4]
}
coef <- sort(coef)
ci_lower <- coef[max(1, round(num_iterations * 0.05))]
ci_upper <- coef[round(num_iterations * 0.95)]


# ------------------------------------------------------------------------------------------------------------
# Run Monte Carlo simulation for power 
simulation_results <- rep(0, 1000)
run_simulation <- function() {
  #pick random subject
  subjects <- sample(unique(data$subject), 43, replace=TRUE)
  
  #pull 100 rows amd add to new data set
  sampledata <- data[0:0,]
  i <- 1
  for (samplesubject in subjects) {
    sampletrials <- sample((data[data$subject == samplesubject,]$trial), 100, replace=TRUE)
    for (eachnumber in sampletrials) {
      row <- data[(data$subject == samplesubject) & (data$trial == eachnumber),]
      row$subject <- i
      sampledata <- rbind(sampledata, row)
    }
    i <- i + 1
  }
  
  #create full and null models and compute chi-square comparison 
  
  
  #if significant count this trial as 1; 0 otherwise
  return(comparison$`p-value`[2] < .05)
}

for (i in 1:1000) {
  simulation_results[i] <- run_simulation()
  print(i)
}
sum(simulation_results)/length(simulation_results)


# Run Monte Carlo simulation for effect size CI
effect_size_results <- rep(0, 1000)

get_effect_size <- function() {
  #pick random subject
  subjects <- sample(unique(data$subject), 43, replace=TRUE)
  
  #pull 100 rows amd add to new data set
  sampledata <- data[0:0,]
  i <- 1
  for (samplesubject in subjects) {
    sampletrials <- sample((data[data$subject == samplesubject,]$trial), 100, replace=TRUE)
    for (eachnumber in sampletrials) {
      row <- data[(data$subject == samplesubject) & (data$trial == eachnumber),]
      row$subject <- i
      sampledata <- rbind(sampledata, row)
    }
    i <- i + 1
  }
  
  #run model and extract relevant statistics for computing effect size
  lmerfit1 <- lmer(rt ~ hand * response + response * space + (1|subject) + (1|trial), data = sampledata)
  randomvar <- VarCorr(lmerfit1)
  print(randomvar, comp=c("Variance"))
  randomvar <- as.data.frame(randomvar)
  trialvar <- randomvar$vcov[1]
  subjectvar <- randomvar$vcov[2]
  resid <- randomvar$vcov[3]
  fixed <- getME(lmerfit1, "beta")
  diff <- fixed[6]
  effectsize <- diff/sqrt(trialvar + subjectvar + resid)
  return(effectsize)
}

for (i in 1:1000) {
  effect_size_results[i] <- get_effect_size()
  print(i)
}
mean(effect_size_results)
effect_size_results <- data.frame(effect_size_results)
effect_size_ordered <- sort(effect_size_results$effect_size_results)
CIUpper <- effect_size_ordered[950]
CILower <- effect_size_ordered[50]


diff_results <- rep(0, 1000)

get_diff <- function() {
  #pick random subject
  subjects <- sample(unique(data$subject), 43, replace=TRUE)
  
  #pull 100 rows amd add to new data set
  sampledata <- data[0:0,]
  i <- 1
  for (samplesubject in subjects) {
    sampletrials <- sample((data[data$subject == samplesubject,]$trial), 100, replace=TRUE)
    for (eachnumber in sampletrials) {
      row <- data[(data$subject == samplesubject) & (data$trial == eachnumber),]
      row$subject <- i
      sampledata <- rbind(sampledata, row)
    }
    i <- i + 1
  }
  
  lmerfit1 <- lmer(rt ~ hand * response + response * space + (1|subject) + (1|trial), data = sampledata)
  fixed <- getME(lmerfit1, "beta")
  diff <- fixed[6]
  return(diff)
}

for (i in 1:1000) {
  diff_results[i] <- get_diff()
  print(i)
}
mean(diff_results)
diff_results <- data.frame(diff_results)
diff_ordered <- sort(diff_results$diff_results)
CIUpper <- diff_ordered[950]
CILower <- diff_ordered[50]