library(tidyverse)
library(nlme)
library(lme4)
library(data.table)


generate.data <- function(data, num_subjects, num_trials) {
  # Pick random subjects
  samplesubjects <- sample(unique(data$subject), num_subjects, replace=TRUE)
  
  # Pull rows amd add to list
  sampledata <- as.list(num_subjects * num_trials)
  simsubject <- 1
  for (samplesubject in samplesubjects) {
    sampletrials <- sample((data[data$subject == samplesubject,]$trial), num_trials, replace=TRUE)
    simtrial <- 1
    for (sampletrial in sampletrials) {
      row <- data[(data$subject == samplesubject) & (data$trial == sampletrial),]
      row$subject <- simsubject
      sampledata[[(simsubject - 1) * num_trials + simtrial]] <- row
      simtrial <- simtrial + 1
    }
    simsubject <- simsubject + 1
  }
  
  # Return list as datatable
  return(data.table::rbindlist(sampledata))
}

# Run Monte Carlo simulation to calculate power and coefficient confidence intervals
run.simulation <- function(data, num_subjects, num_trials) {
  # Generate a simulated dataset
  sim_data <- generate.data(data, num_subjects, num_trials)
  
  # Exp 1 & 2
  # Fit a full and null model and test for significance
  #M1.full <- lme(rt ~ hand * response, random = ~1 | subject, method="ML", data=sim_data)
  #M1.null <- lme(rt ~ hand + response, random = ~1 | subject, method="ML", data=sim_data)
  #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, random = ~1 | subject, method="REML", data=sim_data)
  #coefs <- model$coefficients
  
  # Exp 3
  # Fit a full and null model and test for significance
  M1.full <- lme(rt ~ hand * response + space * hand, random = ~1 | subject, method="ML", data = sim_data)
  M1.null <- lme(rt ~ response + space * hand, random = ~1 | subject, method = "ML", data = sim_data)
  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 + space * hand, random = ~1 | subject, method = "REML", data = sim_data)
  coefs <- model$coefficients
  
  # Return the model and p-value of the anova
  return(list("coefs"=coefs, "p"=p.val))
}

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

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

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

# Print the mean difference
m <- rep(0, num_iterations)
for (i in 1:length(sim_results)) {
  sim_data <- sim_results[[i]]$model$data
  m1 <- mean(sim_data$rt[sim_data$hand == "RIGHT" & sim_data$response == "RIGHT"])
  m2 <- mean(sim_data$rt[sim_data$hand == "LEFT" & sim_data$response == "RIGHT"])
  m[i] <- m1 - m2
}
m <- sort(m)
ci_lower <- m[max(1, round(num_iterations * 0.05))]
ci_upper <- m[round(num_iterations * 0.95)]
ci_lower
ci_upper

# Print the coefs
coef <- rep(0, num_iterations)
for (i in 1:length(sim_results)) {
  coef[i] <- sim_results[[i]]$coefs$fixed[2]
}
coef <- sort(coef)
ci_lower <- coef[max(1, round(num_iterations * 0.05))]
ci_upper <- coef[round(num_iterations * 0.95)]
ci_lower
ci_upper


# ------------------------------------------------------------------------------------------------------------
# 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), 27, 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 + (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[4]
  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), 27, 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 + (1|subject) + (1|trial), data = sampledata)
  fixed <- getME(lmerfit1, "beta")
  diff <- fixed[4]
  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]