# The "Many Coauthors" project reports that these data were collected by Francesca Gino.

# Code by Nick Brown (nicholasjlbrown@gmail.com), November 2023.
# Version history:
# 2023-11-09 21:30 UTC  First public release

# Convert data from SPSS if needed
# rio::convert("190404 - Study X Data.sav", "Paper029-Study3.csv")

df <- read.csv("Paper029-Study3.csv")

# "In 11 sessions, an odd number of individuals showed up"
df <- df[!is.na(df$Durationinseconds_B),]     #11 NAs

# Remove one record (line 84 of the CSV file) where the dependent variable is NA
df <- df[!is.na(df$EmoAcc),]

# Remove one record (line 193 of the CSV file) where the condition is NA
df <- df[!is.na(df$Condition),]

df.full <- df

# Exclude the same participants as the article ("RAs detected inattention").
df.repl <- df[(df$Exclude_LabNotes == 0),]         #37 further exclusions

# Make a dataset to analyse *only* the excluded participants.
df.exc <- df[(df$Exclude_LabNotes == 1),]

cat("Replicating the principal result of the study...\n")
df.repl.dis <- df.repl[df.repl$Condition == 1,]
df.repl.hon <- df.repl[df.repl$Condition == 2,]
tt.repl <- t.test(df.repl.dis$EmoAcc, df.repl.hon$EmoAcc, var.equal=TRUE)
print(tt.repl)

cat("\n")
cat("Calculating the effect for all participants pairs, including those excluded by the RAs...\n")
df.full.dis <- df.full[df.full$Condition == 1,]
df.full.hon <- df.full[df.full$Condition == 2,]
tt.full <- t.test(df.full.dis$EmoAcc, df.full.hon$EmoAcc, var.equal=TRUE)
print(tt.full)

cat("\n")
cat("Calculating the effect for the excluded participant pairs only...\n")
df.exc.dis <- df.exc[df.exc$Condition == 1,]
df.exc.hon <- df.exc[df.exc$Condition == 2,]
tt.exc <- t.test(df.exc.dis$EmoAcc, df.exc.hon$EmoAcc, var.equal=TRUE)
print(tt.exc)

cD <- lsr::cohensD(df.exc.dis$EmoAcc, df.exc.hon$EmoAcc)
cat("Cohen's d for excluded participants = ", sprintf("%.3f", cD), "\n", sep="")

loops <- 0
if (loops <= 0) {
  cat("\n")
  cat("Set the value of \"loops\" to something above 0 to run the simulations.\n")
} else {
  n.full <- nrow(df)
  n.repl <- nrow(df.repl)

  p.actual <- tt.repl$p.value
  p.target <- 0.05

  df.min <- df[c("Condition", "EmoAcc")]    # reduce data to bare minimum for speed

  hits.actual <- 0
  hits.target <- 0
  p.min <- 1

  set.seed(1)
  for (i in 1:loops) {
    if ((i %% 100) == 0) {
      cat(".")    # count progress
    }

    inc <- sample(1:n.full, n.repl)     # random sample of 211 participant pairs out of 248

    df.sim.inc <- df.min[inc,]
    p <- t.test(df.sim.inc[df.sim.inc$Condition == 1,]$EmoAcc, df.sim.inc[df.sim.inc$Condition == 2,]$EmoAcc, var.equal=TRUE)$p.value
    if (p < p.actual) {     # original p value from the article
      hits.actual <- hits.actual + 1
    }
    if (p < p.target) {     # p value that would give a publishable result
      hits.target <- hits.target + 1
    }
    if (p < p.min) {
      p.min <- p
    }
  }
  
  cat("\n")
  cat("Simulated exclusions of random participant pairs:\n")
  cat("  Hits (p < ", sprintf("%.8f", p.actual), "): ", hits.actual, "\n", sep="")
  cat("  Hits (p < ", sprintf("%.8f", p.target), "): ", hits.target, "\n", sep="")
  cat("  Smallest p = ", sprintf("%.8f", p.min), "\n", sep="")
  cat("\n")
}