# Reanalysis of Eskine et al. 10.1177/0956797611398497
# By Nick Brown, nicholasjlbrown@gmail.com, February 2025.

filebase <- "Eskine et al. (2011) psyc sci"
csv.file <- paste0(filebase, ".csv")
spss.file <- paste0(filebase, ".sav")
#rio::convert(spss.file, csv.file)    # This is how I created my CSV file!

df <- read.csv(csv.file, stringsAsFactors=FALSE)
df$Taste <- as.factor(df$Taste)
levels(df$Taste) <- c("Bitter", "Water", "Sweet")
df$subject <- as.factor(df$subject)
df$politic <- as.factor(df$politic)
levels(df$politic) <- c("Conservative", "Liberal", "Moderate", "None")

for (t1 in c("", "bitter_taste", "sweet_taste", "neutral_taste", "Disgust")) {
  cat(t1)
  if (nchar(t1) < 8) {
    cat("\t")
  }

  for (t2 in c(1, 3, 2)) {
    if (t1 == "") {
      m = c("Bitter drink", "Water\t", "Sweet drink")[t2]
      s = ""
    }
    else {
      tt <- df[df$Taste == levels(df$Taste)[t2], t1]
      m = sprintf("%.2f", mean(tt))
      s = sprintf(" (%.2f)", sd(tt))
    }
    cat("\t", m, s, sep="")
  }
  cat("\n", sep="")
}
library(reshape2)
library(ez)
library(emmeans)

# Formatting functions
p_string <- function (p)
{
  result <- if (p < .001) "<.001" else paste0("=", substring(sprintf("%.3f", p), 2))
  return(result)
}

anova_string <- function (F, DFn, DFd, p, ssq1, ssq2)
{
  pesq <- ssq1 / (ssq1 + ssq2)
  result <- paste0("F(", DFn, ",", DFd, ")=", sprintf("%.3f", F), ", p", p_string(p), ", ηp²=", substring(sprintf("%.3f", pesq), 2))
  return(result)
}

ezANOVA_string <- function (model, var)
{
  t <- model$ANOVA[model$ANOVA$Effect == var,]

  result <- anova_string(t$F, t$DFn, t$DFd, t$p, t$SSn, t$SSd)
  return(result)
}

my.cohens_d <- function (v1, v2) {
  n1 <- length(v1)
  n2 <- length(v2)
  pooled.sd <- sqrt((((n1 - 1) * (sd(v1) ^ 2)) + ((n2 - 1) * (sd(v2) ^ 2))) / (n1 + n2 - 2))
  return((mean(v1) - mean(v2)) / pooled.sd)
}

# Build ANOVA-friendly version of data.
meltId <- c("subject", "politic", "gender", "Taste")
meltMeasured <- c(
  "Moral_judge_avg", "Appalling_judge_avg"
)
meltKeep <- c(meltId, meltMeasured)
df.a <- melt(df[meltKeep], id=meltId, measured=meltMeasured)

#options(contrasts=c("contr.sum", "contr.poly"))

for (judgment in c("Moral", "Appalling")) {
  cat("\n", "========== Judgment: ", judgment, " ==========", "\n", sep="")
  cat("Effect of beverage type\n")
  if (judgment == "Appalling") {
    exp.str <- "n/a"
  }

  if (judgment == "Moral") {
    exp.str <- "F(2,51)=7.368, p=.002, ηp²=.224"
  }

  cat(exp.str, " (expected)", "\n", sep="")

  var <- paste0(judgment, "_judge_avg")
  df.aj <- df.a[(df.a$variable == var),]
  anova.j <- suppressMessages(suppressWarnings(
    ezANOVA(dv=value
          , wid=subject
          , between=Taste
          , data=df.aj
          , type=3
          , detailed=TRUE
          , return_aov=TRUE
    )
  ))
  calc.str <- ezANOVA_string(anova.j, "Taste")
  cat(calc.str, " (calculated)", sep="")
  if (exp.str != "n/a") {
    cat("... ", if (calc.str == exp.str) "OK!" else "*** Mismatch ***", sep="")
  }
  cat("\n")

  emm.j <- emmeans(anova.j$aov, ~Taste)
  esumm.j <- summary(emm.j)
  psumm.j <- summary(pairs(emm.j))

  cat("\n")
  cat("Taste\tEMmean\tSD\tContrast 1\t  t\t  p\t  d\tContrast 2\t  t\t  p\t  d\n")
  for (tn in 1:length(levels(df$Taste))) {
    this.taste <- levels(df$Taste)[tn]
    vec <- df[df$Taste == this.taste, var]
    m1.j <- esumm.j[tn, 2]
    pooled.sd.j <- esumm.j[tn, 3] * sqrt(length(vec))
    cat(this.taste
      , "\t", sprintf("%.2f", m1.j)
      , "\t", sprintf("%.2f", sd(vec))
      , sep=""
    )

    for (other.tn in 1:length(levels(df$Taste))) {
      if (other.tn == tn) {
        next
      }
      other.taste <- levels(df$Taste)[other.tn]
      pair <- factor(c(this.taste, other.taste), levels = levels(df$Taste))
      cstring <- paste(pair[order(as.numeric(pair))], collapse = " - ")  # match order of names
      d <- (esumm.j[other.tn, 2] - m1.j) / pooled.sd.j

      for (pr in 1:nrow(psumm.j)) {
        if (psumm.j$contrast[pr] == cstring) {
          cat("\t", "(vs ", other.taste, ")"
            , "\t", sprintf("%.3f", psumm.j$t.ratio[pr])
            , "\t", substring(sprintf("%.3f", psumm.j$p.value[pr]), 2)
            , "\t", sprintf("%5s", sprintf("%.2f", d))
            , sep=""
          )
        }
      }
    }

    cat("\n")
  }

  cat("\n")
  cat("Test whether judgments could be predicted by feelings of physical disgust\n")
  if (judgment == "Moral") {
    exp.str <- "R²=0.275, t(52)=4.445, p<.001, β=0.525"
  }
  if (judgment == "Appalling") {
    exp.str <- "n/a"
  }
  cat(exp.str, " (expected)", "\n", sep="")

  disgust.model <- lm(as.formula(paste0(var, "~Disgust")), data=df)
  sdm <- summary(disgust.model)
  calc.str <- paste0(
        "R²=", sprintf("%.3f", sdm$r.squared)
      , ", t(", sdm$df[2], ")=", sprintf("%.3f", sdm$coefficients["Disgust", 3])
      , ", p", p_string(sdm$coefficients["Disgust", 4])
      , ", β=", sprintf("%.3f", sdm$coefficients["Disgust", 1] * sd(df$Disgust) / sd(df[[var]]))
  )

  cat(calc.str, " (calculated)... ", sep="")
  if (exp.str != "n/a") {
    cat("... ", if (calc.str == exp.str) "OK!" else "*** Mismatch ***", sep="")
  }
  cat("\n")

  if (judgment == "Moral") {
# A 2 ... × 3 between-subjects ANOVA was conducted on moral judgments to determine whether
#  political orientation influenced judgments within each taste condition.
    cat("\n")
    cat("Significant main effect of taste\n")
    exp.str <- "F(2,38)=9.741, p<.001, ηp²=.339"
  }

  cat(exp.str, " (expected)", "\n", sep="")

  df.ajp <- df.aj[(df.aj$politic %in% c("Liberal", "Conservative")),]
  df.ajp$politic <- factor(df.ajp$politic)
  df.ajp$subject <- factor(df.ajp$subject)
  anova.jp <- suppressMessages(suppressWarnings(
    ezANOVA(dv=value
          , wid=subject
          , between=.(Taste, politic)
          , data=df.ajp
          , type=3
          , detailed=TRUE
          , return_aov=TRUE
    )
  ))
  calc.str <- ezANOVA_string(anova.jp, "Taste")
  cat(calc.str, " (calculated)", sep="")
  if (exp.str != "n/a") {
    cat("... ", if (calc.str == exp.str) "OK!" else "*** Mismatch ***", sep="")
  }
  cat("\n")

  emm.jp <- emmeans(anova.jp$aov, ~politic|Taste)
  esumm.jp <- summary(emm.jp)
  psumm.jp <- summary(pairs(emm.jp))

  cat("\n")
  cat("Taste\tMlib\tSDlib\tMcon\tSDcon\t  F\t  p\t  ηp²\n")
  f.ratio <- psumm.jp$t.ratio ^ 2
  p.eta2 <- f.ratio / (f.ratio + psumm.jp$df)

  for (tn in 1:length(levels(df$Taste))) {
    taste <- levels(df$Taste)[tn]
    df.t <- df[df$Taste == taste,]
    vec.con <- df.t[df.t$politic == "Conservative",][[var]]
    vec.lib <- df.t[df.t$politic == "Liberal",][[var]]
    cat(taste
      , "\t", sprintf("%.2f", mean(vec.con))
      , "\t", sprintf("%.2f", sd(vec.con))
      , "\t", sprintf("%.2f", mean(vec.lib))
      , "\t", sprintf("%.2f", sd(vec.lib))
      , "\t", sprintf("%.3f", f.ratio[tn])
      , "\t", sprintf("%.3f", psumm.jp$p.value[tn])
      , "\t", sprintf("%.3f", p.eta2[tn])
      , "\n", sep=""
    )
  }

  cat("\n")
  cat("Politics\tTaste\tM\tSD\n")
  for (pn in 1:2) {
    pol <- levels(df$politic)[pn]
    df.p <- df[df$politic == pol,]
    for (tn in 1:length(levels(df$Taste))) {
      taste <- levels(df$Taste)[tn]
      df.pt <- df.p[df.p$Taste == taste,]
      vec <- df.pt[[var]]
      cat(pol, if (nchar(pol) < 8) "\t" else ""
        , "\t", taste
        , "\t", sprintf("%.2f", mean(vec))
        , "\t", sprintf("%.2f", sd(vec))
        , "\n", sep=""
      )
    }
  }

  cat("\n")
  cat("Politics\t  t\t  df\t  p\t  d\n")
  for (pn in 1:2) {
    pol <- levels(df$politic)[pn]
    df.p <- df[df$politic == pol,]

    vec.dis <- df.p[df.p$Taste == "Bitter",][[var]]
    vec.nodis <- df.p[df.p$Taste != "Bitter",][[var]]
    tt <- t.test(vec.dis, vec.nodis, var.equal=TRUE)
    d <- my.cohens_d(vec.dis, vec.nodis)
    cat(pol, if (nchar(pol) < 8) "\t" else ""
      , "\t", sprintf("%.3f", tt$statistic)
      , "\t", sprintf("%4s", tt$parameter)
      , "\t", sprintf("%.3f", tt$p.value)
      , "\t", sprintf("%.3f", d)
      , "\n", sep=""
    )
  }
}

vignette.cols <- c(9, 11, 13, 15, 17, 19)

split.pol <- FALSE   # set to TRUE to split alphas by political stance
cat("\n")
cat("Judgment\tTaste\t", if (split.pol) "Politics\t" else "", " Alpha\n", sep="")
for (judgment in c("Moral", "Appalling")) {
  if (judgment == "Moral") {
    vcols <- vignette.cols
  }
  if (judgment == "Appalling") {
    vcols <- vignette.cols + 1
  }

  for (taste in levels(df$Taste)) {
    for (pol in c("Any", "Conservative" ,"Liberal")) {
      if (pol == "Any") {
        pol.subset <- TRUE
      }
      else {
        pol.subset <- (df$politic == pol)
      }
      ratings <- df[((df$Taste == taste) & pol.subset), vcols]
      jtp.alpha <- ltm::cronbach.alpha(ratings)
      cat(judgment, if (nchar(judgment) < 8) "\t" else ""
        , "\t" , taste
        , if (split.pol) paste0("\t", pol, if (nchar(pol) < 8) "\t" else "") else ""
        , "\t", sprintf("%6s", sprintf("%.3f", jtp.alpha$alpha))
        , "\n", sep=""
      )
      if (! split.pol) {
        break
      }
    }
  }
}

df.bitter <- df[df$Taste == "Bitter",]
df.notbitter <- df[df$Taste != "Bitter",]

moral.bitter <- unlist(as.vector(df.bitter[vignette.cols]))
appal.bitter <- unlist(as.vector(df.bitter[vignette.cols + 1]))
moral.notbitter <- unlist(as.vector(df.notbitter[vignette.cols]))
appal.notbitter <- unlist(as.vector(df.notbitter[vignette.cols + 1]))