Pre-Registration: Data Cleaning & Analysis Script for Intuitive Behavioural Study of Moral Judgements of Humans vs. AI

STEP 0: Setup

Set Working Directory

wd <- "~"

The code below determines your computer’s operating system, and uses the specific file path accordingly.

set.platform <- function(subdir = "") {
  base_wd <- if (Sys.info()[["sysname"]] == "Darwin") {
    file.path("~/",wd)
  } else if (Sys.info()[["sysname"]] == "Windows") {
    file.path("C:/Users/~/",wd)
  } else {
    stop("Unknown operating system, set your working directory manually.")
  }
  # allow empty subdir for top level, or paste others
  full_wd <- if (subdir == "") base_wd else file.path(base_wd, subdir)
  setwd(full_wd)
}

Load All Packages

# Package names
packages <- c("MASS", "Matrix", "dplyr","faux", "ggplot2", "ggeffects", "tidyverse", "lme4", "lmerTest", "simr", "MBESS","future","future.apply","binom", "see","optimx", "afex","papaja","devtools","emmeans","stargazer","haven","reshape2","bayestestR", "sjPlot","parameters","parallel","performance","ltm")

# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Packages loading
library(dplyr)
library(faux)
library(ggplot2)
library(lme4)
library(lmerTest)
library(tidyverse)
library(simr)
library(MBESS)
library(future)
library(future.apply)
library(binom)
library(see)
library(afex)
library(Matrix)
library(devtools)
install_local("~/Library/CloudStorage/OneDrive-UniversityofExeter/PhD/Behavioural Study/sampleSize_simulations/faux_1.2.1.tar.gz")
library(faux)
library(MASS)
library(emmeans)
library(stargazer)
library(haven)
library(reshape2)
library(bayestestR)
library(sjPlot)
library(parameters)
library(parallel)
library(performance)
library(ltm)

Set Seed

set.seed(123)

STEP 1: Cleaning

Read in Data

# Human 
# Where key for yes=F
dataF.h <- read.csv("xx.csv", header=TRUE, na.strings="")
# Where key for yes=J
dataJ.h <- read.csv("xx.csv", header=TRUE, na.strings="")
# AI 
# Where key for yes=F
dataF <- read.csv("xx.csv", header=TRUE, na.strings="")
# Where key for yes=J
dataJ <- read.csv("xx.csv", header=TRUE, na.strings="")
# demographics
demo <- read.csv("xx.csv", header=TRUE, na.strings="")

Join Key-Segregated Data

# combine data
df <- rbind(dataF.h,dataJ.h,dataF, dataJ)

# keep only data rows with moral judgements, remove irrelevant rows (e.g., instructions, etc)
df <- df %>% 
  filter(Response.Type == "response")

# select data we need 
df <- df %>% 
  dplyr::select(Participant.Private.ID, Task.Name, Trial.Number, Spreadsheet..Dilemma, Spreadsheet..Context, Spreadsheet..DilemmaByContext, Spreadsheet..Prohibited.act, Spreadsheet..PDE, Spreadsheet..PDE.noBA, Spreadsheet..PDE.GE.BE, Spreadsheet..PDE.SE_notME, Spreadsheet..MG.Permissible, Spreadsheet..Agent, Absolute.Reaction.Time, Response, randomiser.zmvq)

# rename columns to make them easier to work with
df <- df %>%
  rename(PBH = Spreadsheet..Prohibited.act,
         PDE = Spreadsheet..PDE,
         BA = Spreadsheet..PDE.noBA,
         MC = Spreadsheet..PDE.GE.BE,
         intent = Spreadsheet..PDE.SE_notME,
         MG = Spreadsheet..MG.Permissible,
         RT = Absolute.Reaction.Time,
         ID = Participant.Private.ID,
         MJ = Response,
         Agent = Spreadsheet..Agent,
         dilemma = Spreadsheet..DilemmaByContext,
         dilemmaAll = Spreadsheet..Dilemma,
         context = Spreadsheet..Context,
         keys = randomiser.zmvq)

# filter demo to only include rows from the Gorilla task where generic demographics were collected
demo <- demo %>% 
  filter(Task.Name == "Generic Demographics")

# keep only data rows in demographics we need not all the other stuff gorilla gives you
demo <- demo %>% 
  filter(Question.Key == "Gender"| Question.Key == "Age" | Question.Key == "cultural-arab"| Question.Key == "cultural-arab-text"|  Question.Key == "cultural-asian"| Question.Key == "cultural-asian-text"| Question.Key == "cultural-black"| Question.Key == "cultural-black-text"|  Question.Key == "cultural-multiple"| Question.Key == "cultural-multiple-text"| Question.Key == "cultural-white"| Question.Key == "cultural-white-text"| Question.Key == "cultural-prefNo"| Question.Key == "cultural-prefNo-text"|  Question.Key == "cultural-another"| Question.Key == "cultural-another-text")

# separate the age and gender into two separate variables - wide form
demowide<-spread(data=demo, key=Question.Key, value=MJ)

# now add the demographics to the data
df <- merge(df,demowide, 
            by="Participant.Private.ID")

We add in the Nationality from Prolific, like our previous study for any exploratory analyses.

# load in demo from prolific
prolific_demo <- read.csv("~.csv")

# filter prolific_demo by column: Status so removes rows w. "RETURNED" or "TIMED-OUT" 
prolific_demo <- prolific_demo %>% 
  filter(Status != "RETURNED" & Status != "TIMED-OUT")

# next, we create a shared column between the two datasets
# prolific incl. a column called Participant.id which corresponds to the ID column called Participant.Public.ID in demo data
# therefore, we will rename the Participant.Public.ID column in demo to Participant.id
demo <- demo %>%
  rename(Participant.id = "Participant.Public.ID")

# now we can add the nationality data from prolific_demo to demo using the Participant.id column to ensure these are added to the correct Ps
demo <- demo %>%
  left_join(prolific_demo %>% dplyr::select(Participant.id, Nationality), by = "Participant.id")

Now let’s sort the cultural demographic details that we collected directly from Gorilla.

# pivot longer for .text columns (get ID, question_col, and text_response)
text_cols <- grep("\\.text$", names(demo), value = TRUE)

text_long <- demo %>%
  dplyr::select(Participant.Private.ID, all_of(text_cols)) %>%
  pivot_longer(cols = -Participant.Private.ID, names_to = "text_col", values_to = "text_response") %>%
  filter(!is.na(text_response) & text_response != "")

# get cultural columns starting with the prefix but excluding .text or .quantised
cultural_cols <- grep("^questionnaire\\.qemh\\.cultural\\.", names(demo), value = TRUE)
cultural_cols <- cultural_cols[!str_detect(cultural_cols, "\\.text$|\\.quantised$")]

# for ps that put other / multiple, we might need to manually add these into the correct column
demo <- demo %>% 
  mutate(
    questionnaire.qemh.cultural.white = case_when(
      Participant.Private.ID == "14606099" ~ "British or English or Scottish or Welsh or Northern Irish",
      Participant.Private.ID == "14606101" ~ "Dutch",
      Participant.Private.ID == "14606130" ~ "Portuguese",
      Participant.Private.ID == "14606511" ~ "Italian",
      .........
      TRUE ~ questionnaire.qemh.cultural.white
    ),
    questionnaire.qemh.cultural.multiple = case_when(
      Participant.Private.ID == "14606099"  ~ NA_character_,
      Participant.Private.ID == "14606626" ~ "Latino",
      Participant.Private.ID == "14609697" ~ "Mexican-Japanese",
      TRUE ~ questionnaire.qemh.cultural.multiple
    ),
    questionnaire.qemh.cultural.arab = case_when(
      ...
      .......
      TRUE ~ questionnaire.qemh.cultural.arab
    ),
    questionnaire.qemh.cultural.asian = case_when(
      ...
      TRUE ~ questionnaire.qemh.cultural.asian
    ),
    questionnaire.qemh.cultural.another = case_when(
      ...
      TRUE ~ questionnaire.qemh.cultural.another
    ),
  )

Next, we remove all the bits we no longer need.

# remove values not needed after manual check
#convert all columns to character class
demo <- demo %>%
  mutate(across(everything(), as.character))
# Identify columns to modify: all except those ending with .text or .quantised
cols_to_modify <- names(demo)[!str_detect(names(demo), "\\.text$|\\.quantised$")]

demo <- demo %>%
  mutate(across(all_of(cols_to_modify),
                ~ if_else(. == "In another way (specify, if you wish).", NA_character_, .)))

# rename the column headers to something shorter and more intuitive 
demo <- demo %>% 
  rename(
    ID = "Participant.Private.ID",
    Gender = "questionnaire.qemh.Gender",
    Gender.quant = "questionnaire.qemh.Gender.quantised",
    Age = "questionnaire.qemh.Age",
    cultural.arab = "questionnaire.qemh.cultural.arab",
    cultural.asian = "questionnaire.qemh.cultural.asian",
    cultural.black = "questionnaire.qemh.cultural.black",
    cultural.multiple = "questionnaire.qemh.cultural.multiple",
    cultural.white = "questionnaire.qemh.cultural.white",
    cultural.prefNo = "questionnaire.qemh.cultural.prefNo",
    cultural.another = "questionnaire.qemh.cultural.another",
    Nationality = "Nationality"
  )

# select only columns we renamed above
demo <- demo %>% 
  dplyr::select(ID, Gender, Gender.quant, Age, Nationality, cultural.arab, cultural.asian, cultural.black,
                cultural.multiple, cultural.white, cultural.prefNo, cultural.another)

demo <- demo %>%
  mutate(cultural.another = if_else(cultural.another == "In another way (specify, if you wish)", 
                                    NA_character_, cultural.another)) %>% 
  mutate(cultural.multiple = if_else(cultural.multiple == "Any other mixed or multiple ethnic backgrounds (specify, if you wish).",
                                     NA_character_, cultural.multiple))



cultural_cols <- c("cultural.arab", "cultural.asian", "cultural.black", 
                   "cultural.multiple", "cultural.white", "cultural.prefNo", 
                   "cultural.another")

# Assign primary cultural category as first non-empty for each participant
demo <- demo %>%
  mutate(
    primary_cultural = case_when(
      !is.na(cultural.white) & cultural.white != "" ~ "White",
      !is.na(cultural.black) & cultural.black != "" ~ "Black",
      !is.na(cultural.asian) & cultural.asian != "" ~ "Asian",
      !is.na(cultural.arab) & cultural.arab != "" ~ "Arab",
      !is.na(cultural.multiple) & cultural.multiple != "" ~ "Multiple",
      !is.na(cultural.prefNo) & cultural.prefNo != "" ~ "No preference",
      !is.na(cultural.another) & cultural.another != "" ~ "Other (specified)",
      TRUE ~ NA_character_
    ),
    multiple_cultural = rowSums(across(all_of(cultural_cols), ~ !is.na(.) & . != "")) > 1
  )

demo <- demo %>%
  rowwise() %>%
  mutate(
    identified_prefs = list(
      c_across(all_of(cultural_cols)) %>%
        keep(~ !is.na(.) && . != "")  # Keep only non-NA non-empty values
    ),
    prefs_string = ifelse(length(identified_prefs) == 0, 
                             NA_character_, 
                             paste(identified_prefs, collapse = "; "))
  ) %>%
  ungroup()

# remove "DATA_EXPIRED" value from Nationality
demo <- demo %>%
  mutate(
    Nationality = if_else(Nationality == "DATA_EXPIRED", NA_character_, Nationality)
  )


# separate the age and gender into two separate variables - wide form
dfDemowide<-spread(data=demo, key=Question.Key, value=Response)

# finally, we add the demographics to the behavioural data
df <- merge(df,dfDemowide, 
            by="ID")

With the data combined, we can now start calculating our demographics.

# calculating mean age
dfDemowide$Age <- as.numeric(dfDemowide$Age)
mean_age <- mean(dfDemowide$Age)
mean_age

# calculating standard deviation of age
sd(dfDemowide$Age)

# find the age of the youngest participant
youngest_age <- min(dfDemowide$Age)
youngest_age

# find the age of the oldest participant
oldest_age <- max(dfDemowide$Age)
oldest_age

# change Other (please specify) to what was specified in next question
dfDemowide$Gender <- ifelse(dfDemowide$Gender == "Other (please specify)", "Non Binary", dfDemowide$Gender)

# get count of number of participants of each gender
gender_count <- table(dfDemowide$Gender)
gender_count

woman_count <- gender_count["Woman"]
woman_count

man_count <- gender_count["Man"]
man_count

nonbinary_count <- gender_count["Non-binary"]
nonbinary_count
prefNoGender_count <- gender_count["Prefer Not To Say"]
prefNoGender_count
anotherWayGender_count <- gender_count["In another way (specify, if you wish)"] 
anotherWayGender_count

print(gender_count)

# cultural demographics

cultural_cols <- c("cultural.arab", "cultural.asian", "cultural.black", 
                   "cultural.multiple", "cultural.white", "cultural.prefNo", 
                   "cultural.another")

# Assign primary cultural category as first non-empty for each participant
demo <- demo %>%
  mutate(
    primary_cultural = case_when(
      !is.na(cultural.white) & cultural.white != "" ~ "White",
      !is.na(cultural.black) & cultural.black != "" ~ "Black",
      !is.na(cultural.asian) & cultural.asian != "" ~ "Asian",
      !is.na(cultural.arab) & cultural.arab != "" ~ "Arab",
      !is.na(cultural.multiple) & cultural.multiple != "" ~ "Multiple",
      !is.na(cultural.prefNo) & cultural.prefNo != "" ~ "No preference",
      !is.na(cultural.another) & cultural.another != "" ~ "Other (specified)",
      TRUE ~ NA_character_
    ),
    multiple_cultural = rowSums(across(all_of(cultural_cols), ~ !is.na(.) & . != "")) > 1
  )

# Summary table excluding NAs in primary_cultural
summary_table <- demo %>%
  filter(!is.na(primary_cultural)) %>%   # Exclude NAs here
  count(primary_cultural, multiple_cultural) %>%
  mutate(percent = round(100 * n / sum(n), 1),
         label = paste0(n, " (", percent, "%)"))

print("Summary of primary cultural group and multiple identity flag excluding NAs:")
print(summary_table)


cultural_cols <- c("cultural.arab", "cultural.asian", "cultural.black", 
                   "cultural.multiple", "cultural.white", "cultural.prefNo", 
                   "cultural.another")

demo <- demo %>%
  rowwise() %>%
  mutate(
    identified_prefs = list(
      c_across(all_of(cultural_cols)) %>%
        keep(~ !is.na(.) && . != "")  # Keep only non-NA non-empty values
    ),
    prefs_string = ifelse(length(identified_prefs) == 0, 
                             NA_character_, 
                             paste(identified_prefs, collapse = "; "))
  ) %>%
  ungroup()

# View results: 
cultural_count <- demo %>%
  dplyr::select(ID, primary_cultural, prefs_string) %>%
  dplyr::arrange(primary_cultural,prefs_string)
cultural_count

Prep Variables for Analysis

# make sure everything is factor
df$Task.Name <- as.factor(df$Task.Name)
df$dilemma <- as.factor(df$dilemma)
df$agent <- as.factor(df$agent)
df$context <- as.factor(df$context)
df$PDE <- as.factor(df$PDE)
df$BA <- as.factor(df$BA)
df$MC <- as.factor(df$MC)
df$intent <- as.factor(df$intent)
df$PBH <- as.factor(df$PBH)
df$MG <- as.factor(df$MG)
df$RT <- abs(as.numeric(df$RT))

Contrasts Coding

#set deviation contrasts for ease of interpretation -.5 vs .5
c<-contr.treatment(2)
my.coding<-matrix(rep(1/2, 2), ncol=1)
my.simple<-c-my.coding
my.simple

Set the Contrast Coding Per Variable

contrasts(df$PDE)<-my.simple
contrasts(df$PDE)
contrasts(df$BA)<-my.simple
contrasts(df$BA)
contrasts(df$MC)<-my.simple
contrasts(df$MC)
contrasts(df$intent)<-my.simple
contrasts(df$intent)
contrasts(df$PBH)<-my.simple
contrasts(df$PBH)
contrasts(df$MG)<-my.simple
contrasts(df$MG)
contrasts(df$agent)<- -1 * my.simple # -1 * makes human = -0.5, rather than AI
contrasts(df$agent)

Re-code Response Variable

df <- df %>%
  mutate(MJ = recode(MJ, "Yes" = "1",
                     "No" = "0"))

df$MJ <- as.factor(df$MJ)

Clean Data

write.csv(df, "cleanedData.csv", row.names = FALSE)

STEP 2: Data Analysis

Mixed Effects Logistic Regression (Three-Way) w. Moral Judgements as DV

hVai <- glmer(MJ ~ PBH*PDE*agent +  (1+PDE*PBH*agent|ID) + (1+agent|dilemma),  data=df, family = binomial(link = "logit"),
               control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))

summary(hVai)

confint(hVai)

Estimated Marginal Means

# by agent
emmeans_hVai <- emmeans(hVai, pairwise ~ PBH*PDE|agent, cov.reduce = range)
emmeans_hVai

CIhVai <- confint(emmeans_hVai)
CIhVai

# by PBH
emmeans_PBHvN <- emmeans(hVai, pairwise ~ agent*PDE|PBH, cov.reduce = range)
emmeans_PBHvN

CIPBHvN <- confint(emmeans_PBHvN)
CIPBHvN

# by PDE
emmeans_PDEvN <- emmeans(hVai, pairwise ~ agent*PBH|PDE, cov.reduce = range)
emmeans_PDEvN

CIhVai <- confint(emmeans_PDEvN)
CIhVai

Simple Contrasts

# by agent
contrasts_hVai <- contrast(emmeans_hVai, interaction = "pairwise", by=c("agent"))
contrasts_hVai

# by PBH
contrasts_PBHvN <- contrast(emmeans_PBHvN, interaction = "pairwise", by=c("PBH"))
contrasts_PBHvN

# by PDE
contrasts_PDEvN <- contrast(emmeans_PDEvN, interaction = "pairwise", by=c("PDE"))
contrasts_PDEvN

Plots

plot(allEffects(hVai))

set_theme(base = theme_classic())
plot_model(hVai, terms = "PBH", type = "pred")
plot_model(hVai, terms = "PDE", type = "pred")
plot_model(hVai, terms = "agent", type = "pred")
plot_model(hVai, type = "int")

Estimated Marginal Means

emmeans_hVai <- emmeans(hVaiRT, pairwise ~ PBH*PDE|agent, cov.reduce = range)

CIhVaiRT <- confint(emmeans_hVai)
CIhVaiRT

Simple Contrasts

contrasts_hVai <- contrast(emmeans_hVai, interaction = "pairwise", by=c("agent"))
contrasts_hVai

Subset Data

ai.df <- subset(df, agent == "AI") 
human.df <- subset(df, agent == "Human")

Mixed Effects Logistic Regression (Two-Way)

Human

Here we reduce the interaction of PDE and PBH as well for the ID random slope.

h_glm <- glmer(MJ ~ PBH*PDE +  (1+PDE+PBH|ID) + (1|dilemma),  data=human.df, family = binomial(link = "logit"),
               control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))

summary(h_glm)

confint(h_glm)

AI

Here we reduce the interaction of PBH and PDE as well for the ID random slope.

ai_glm <- glmer(MJ ~ PBH*PDE +  (1+PDE+PBH|ID) + (1|dilemma),  data=ai.df, family = binomial(link = "logit"),
               control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))

summary(ai_glm)

confint(ai_glm)

Mixed Effects Linear Regression (Three-Way) w. RT as DV

hVaiRT <- glmer(RT ~ PBH*PDE*agent +  (1+PDE*PBH*agent|ID) + (1+agent|dilemma),  data=df)

summary(hVaiRT)

confint(hVaiRT)

STEP 3: Regression Tables

Convert Log Odds to Odds Ratios so they are easier to interpret

OR_hVai <- exp(coef(summary(hVai))[, "Estimate"])

lower <- exp(confint(hVai)[,1])
upper <- exp(confint(hVai)[,2])

Subsetted Data

Human

OR_h <- exp(coef(summary(h_glm))[, "Estimate"])

lower <- exp(confint(h_glm)[,1])
upper <- exp(confint(h_glm)[,2])

AI

OR_ai <- exp(coef(summary(ai_glm))[, "Estimate"])

lower <- exp(confint(ai_glm)[,1])
upper <- exp(confint(ai_glm)[,2])

Change LmerTest Summary to LmerMod for Stargazer

class(hVai) <- "lmerMod"

Subsetted Data

Human

class(h_glm) <- class(ai_glm) <- "lmerMod"

Create Regression Table for Moral Judgements

stargazer(hVai, type = "html", title="Logistic GLMM Human v AI MJs explained by MG", 
          ci=TRUE, intercept.bottom = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digit.separator = "", out = "hVai_MJ_table.html")

Subsetted Data

Human

stargazer(h_glm, type = "html", title="Logistic GLMM Human-Only MJs explained by MG", 
          ci=TRUE, intercept.bottom = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digit.separator = "", out = "h_MJ_table.html")

AI

stargazer(ai_glm, type = "html", title="Logistic GLMM AI-Only MJs explained by MG", 
          ci=TRUE, intercept.bottom = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digit.separator = "", out = "ai_MJ_table.html")

STEP 4: TOaST Tests

TOaST <- equivalence_test(hVai, rule = "classic")
TOaST
plot(TOaST)

# set bounds to OR of 1.68 (log(1.68)=0.52) 
TOaST <- equivalence_test(hVai, rule = "classic", range = c(-0.52,0.52))
TOaST
plot(TOaST)

Subsetted Data

Human

TOaST_h <- equivalence_test(h_glm, rule = "classic")
TOaST_h
plot(TOaST_h)

# set bounds to OR of 1.68 (log(1.68)=0.52)
TOaST_h <- equivalence_test(h_glm, rule = "classic", range = c(-0.52,0.52))
TOaST_h
plot(TOaST_h)

AI

TOaST_ai <- equivalence_test(ai_glm, rule = "classic")
TOaST_ai
plot(TOaST_ai)

# set bounds to OR of 1.68 (log(1.68)=0.52)
TOaST_ai <- equivalence_test(ai_glm, rule = "classic", range = c(-0.52,0.52))
TOaST_ai
plot(TOaST_ai)

STEP 5: Item-Level Analyses

I’m interested in whether each dilemma matches the prediction of MG (for each agent), and further, whether these varies across participants.

Therefore, I use the full and subsetted data here to test variations of hypotheses re. each agent. Then I add in the prediction of the moral grammar to all df

human.df$MG_pred <- ifelse(human.df$PBH == "Yes" & human.df$PDE == "No", 0, 1)
ai.df$MG_pred <- ifelse(ai.df$PBH == "Yes" & ai.df$PDE == "No", 0, 1)
df$MG_pred <- ifelse(df$PBH == "Yes" & df$PDE == "No", 0, 1)

Now let’s add in a variable that tells us whether or not the Ps MJ matches the MG prediction:

human.df$match_MG <- ifelse(human.df$MJ == human.df$MG, 1, 0)
ai.df$match_MG <- ifelse(ai.df$MJ == ai.df$MG, 1, 0)

df$match_MG <- ifelse(df$MJ == df$MG, 1, 0)

Finally, let’s analyse.

full_item.level <- glmer(MJ ~ MG  + (1|ID) + (1|dilemma),
               data = df, family = binomial)
summary(full_item.level)


human_item.level <- glmer(MJ ~ MG + (1|ID) + (1|dilemma),
               data = human.df, family = binomial)
summary(human_item.level)


ai_item.level <- glmer(MJ ~ MG + (1|ID) + (1|dilemma),
               data = ai.df, family = binomial)
summary(ai_item.level)

# dilemma
full_dilemma.level <- glmer(MJ ~ MG  +  (1|dilemma),
               data = df, family = binomial)
summary(full_dilemma.level)
ranef(full_dilemma.level)$dilemma


ai_dilemma.level <- glmer(MJ ~ MG  +  (1|dilemma),
               data = ai.df, family = binomial)
summary(ai_dilemma.level)
ranef(ai_dilemma.level)$dilemma

human_dilemma.level <- glmer(MJ ~ MG  +  (1|dilemma),
               data = human.df, family = binomial)
summary(human_dilemma.level)
ranef(human_dilemma.level)$dilemma

# id
full_ID.level <- glmer(MJ ~ MG  +  (1|ID),
               data = df, family = binomial)
summary(full_ID.level)

From Original Model

dilemma_effects <- ranef(hVai)$dilemma
dilemma_effects

id_effects <- ranef(hVai)$ID
id_effects

# extract fixed effect estimates
fixef_intercept <- fixef(hVai)[1]
fixef_slope <- fixef(hVai)[2]

# Base fixed effect line
fixed_line <- fixef_intercept + fixef_slope

STEP 6: Model Comparisons

First, let’s fit the reduced models.

# agent only
agent <- glmer(MJ ~ agent +  (1+agent|ID) + (1+agent|dilemma),  data=df, family = binomial(link = "logit"),
              control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(agent)
# PBH only
PBH <- glmer(MJ ~ PBH +  (1+PBH|ID) + (1|dilemma),  data=df, family = binomial(link = "logit"),
               control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(PBH)
# PDE only
PDE <- glmer(MJ ~ PDE +  (1+PDE|ID) + (1|dilemma),  data=df, family = binomial(link = "logit"),
               control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(PDE)

# PBH+agent only
PBH_agent <- glmer(MJ ~ PBH*agent +  (1+PBH*agent|ID) + (1+agent|dilemma),  data=df, family = binomial(link = "logit"),
             control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(PBH_agent)
# PDE+agent only
PDE_agent <- glmer(MJ ~ PDE*agent +  (1+PDE*agent|ID) + (1+agent|dilemma),  data=df, family = binomial(link = "logit"),
             control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(PDE_agent)

# PBH+PDE only
PBH_PDE <- glmer(MJ ~ PBH*PDE +  (1+PBH*PDE|ID) + (1|dilemma),  data=df, family = binomial(link = "logit"),
             control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(PBH_PDE)

compare_performance(agent, PBH, PDE, PBH_agent, PDE_agent, PBH_PDE, hVai)

STEP 7: Intent

intent <- glmer(MJ ~ intent +  (1+intent|ID) + (1|dilemma),  data=df, family = binomial(link = "logit"),
              control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(intent)

hVai.intent <- glmer(MJ ~ intent*agent +  (1+intent+agent|ID) + (1+agent|dilemma),  data=df, family = binomial(link = "logit"),
              control = glmerControl(optCtrl = list(maxfun = 2e5), optimizer = "bobyqa"))
summary(hVai.intent)

TOaST.i <- equivalence_test(hVai.intent, rule = "classic", range = c(-0.52,0.52))
TOaST.i
plot(TOaST.i)