library(emmeans)
library(tidyverse)
library(dplyr)
library(car)
library(stringr)
library(lmerTest)
library(MASS)
library(ordinal)
library(MuMIn)
library(psych)##### Import data
raw_data_1 = read.csv('asd1_merged.csv')
##### Clean data
##### Exclude according to 3 criteria:
##### 1) Fail to complete study (or took it multiple times);
length(unique(raw_data_1$PROLIFIC_PID)) # 403 different people took the survey## [1] 403max(table(raw_data_1$PROLIFIC_PID)) # The max number of trials a participant had was 32## [1] 32sum(table(raw_data_1$PROLIFIC_PID) == 20) # 398 participants had 20 trials## [1] 398table(raw_data_1$PROLIFIC_PID)[table(raw_data_1$PROLIFIC_PID) != 20] # 5 had a different number of trials##
## 5589cdf8fdf99b18bd86d031 5e55463217d334029c5d9a35 62029c754dc28df6b6246e96
## 1 4 4
## 667ad3a189817cf129ae61e6 66e06758c868de3476e86e02
## 32 26# I keep only those with 20 trials (to avoid incomplete runs and participants taking the survey multiple times)
clean_data_1 <- raw_data_1[raw_data_1$PROLIFIC_PID %in% names(table(raw_data_1$PROLIFIC_PID)[table(raw_data_1$PROLIFIC_PID) == 20]), ]
any(table(clean_data_1$PROLIFIC_PID) != 20)## [1] FALSElength(unique(clean_data_1$PROLIFIC_PID)) # 398 participants are left## [1] 398##### 2) Complete whole study in < 120 seconds (minimum completion time)
# I cannot find this info in the data, so I do not perform exclusions
##### 3) Provide four or more responses in under 500 milliseconds (minimum response time, for each DV)
# Violation judgment
sum(clean_data_1$rt < 500, na.rm = TRUE) # There are 6 trials with rt < 500 ms for violation judgments## [1] 6table(clean_data_1$PROLIFIC_PID[clean_data_1$rt < 500]) # 5 of those trials come from 1 participant alone##
## 66904bbd1eb9bb790b9b0f6f 66bf5fde0a42438253e0d94e
## 5 1# Self-reported confidence
sum(clean_data_1$rt_conf < 500, na.rm = TRUE) # There are 0 trials with rt < 500 ms for self-reported confidence judgments## [1] 0# Literal and Moral judgments
sum(clean_data_1$lm_rt < 500, na.rm = TRUE) # There are 0 trials with rt < 500 ms for literal & moral judgments## [1] 0# I exclude that 1 participant with 5 trials < 500 ms (all 20 trials, not just those 5)
clean_data_1 <- clean_data_1[!(clean_data_1$PROLIFIC_PID %in%
names(table(clean_data_1$PROLIFIC_PID[clean_data_1$rt < 500])[
table(clean_data_1$PROLIFIC_PID[clean_data_1$rt < 500]) >= 4
])), ]
length(unique(clean_data_1$PROLIFIC_PID)) # 397 participants are left## [1] 397##### Additionally exclude those whose Prolific says they are ASD, but they claim they are not
#sum(clean_data_1$group == "ASD" &
# !str_detect(clean_data_1$comorbidities, fixed("ASD (Autism Spectrum Disorder)")),
# na.rm = TRUE) # 1038 rows (52 participants, with 1 showing 18 instead of 20 trials because they seemingly did not answer the comorbidity question at all)
#
#unique(clean_data_1$PROLIFIC_PID[
# clean_data_1$group == "ASD" &
# !str_detect(clean_data_1$comorbidities, fixed("ASD (Autism Spectrum Disorder)"))
#])
#
#table(clean_data_1$PROLIFIC_PID[
# clean_data_1$group == "ASD" &
# !str_detect(clean_data_1$comorbidities, fixed("ASD (Autism Spectrum Disorder)"))
#])
# I exclude those 52 participants
#clean_data_1 <- clean_data_1[!(
# clean_data_1$PROLIFIC_PID %in% unique(clean_data_1$PROLIFIC_PID[
# clean_data_1$group == "ASD" &
# !str_detect(clean_data_1$comorbidities, fixed("ASD (Autism Spectrum Disorder)"))
# ])
#), ]
#length(unique(clean_data_1$PROLIFIC_PID)) # 345 participants are left
# Make categorical variables factors
str(clean_data_1$gender)## chr [1:7940] "Female" "Female" "Female" "Female" "Female" "Female" ...clean_data_1$gender <- factor(clean_data_1$gender)
str(clean_data_1$gender)## Factor w/ 5 levels "","Female","Male",..: 2 2 2 2 2 2 2 2 2 2 ...str(clean_data_1$group)## chr [1:7940] "NT" "NT" "NT" "NT" "NT" "NT" "NT" "NT" "NT" "NT" "NT" "NT" ...clean_data_1$group <- factor(clean_data_1$group, levels = c("NT", "ASD"))
str(clean_data_1$group)## Factor w/ 2 levels "NT","ASD": 1 1 1 1 1 1 1 1 1 1 ...str(clean_data_1$text)## int [1:7940] 1 0 0 1 0 1 0 1 0 0 ...clean_data_1$text <- factor(clean_data_1$text)
str(clean_data_1$text)## Factor w/ 2 levels "0","1": 2 1 1 2 1 2 1 2 1 1 ...str(clean_data_1$purpose)## int [1:7940] 0 1 0 1 1 0 0 1 0 1 ...clean_data_1$purpose <- factor(clean_data_1$purpose)
str(clean_data_1$purpose)## Factor w/ 2 levels "0","1": 1 2 1 2 2 1 1 2 1 2 ...str(clean_data_1$letter_response)## int [1:7940] 5 0 0 5 0 5 0 5 0 0 ...clean_data_1$letter_response <- factor(clean_data_1$letter_response)
str(clean_data_1$letter_response)## Factor w/ 6 levels "0","1","2","3",..: 6 1 1 6 1 6 1 6 1 1 ...class(clean_data_1$letter_response)## [1] "factor"clean_data_1$letter_response <- as.ordered(clean_data_1$letter_response)
class(clean_data_1$letter_response)## [1] "ordered" "factor"str(clean_data_1$moral_response)## int [1:7940] 4 3 0 4 4 0 0 5 0 4 ...clean_data_1$moral_response <- factor(clean_data_1$moral_response)
str(clean_data_1$moral_response)## Factor w/ 6 levels "0","1","2","3",..: 5 4 1 5 5 1 1 6 1 5 ...class(clean_data_1$moral_response)## [1] "factor"clean_data_1$moral_response <- as.ordered(clean_data_1$moral_response)
class(clean_data_1$moral_response)## [1] "ordered" "factor"# Recompute cw_resp with the new formula (response_confidence * (response - .5) * 2)
clean_data_1 <- clean_data_1 %>%
mutate(cw_resp = confidence * (response - 0.5) * 2)
# Anonymize participants’ ID
clean_data_1$subject_nr <- match(clean_data_1$PROLIFIC_PID,
unique(clean_data_1$PROLIFIC_PID))
clean_data_1$PROLIFIC_PID <- NULL
# Correct columns’ names
clean_data_1 <- clean_data_1 %>%
rename(scene = rule)
# Recompute AQ-10 scores
# 0 = Definitely Agree
# 1 = Slightly Agree
# 2 = Slightly Disagree
# 3 = Definitely Disagree
# SCORING: Only 1 point can be scored for each question. Score 1 point for Definitely or
# Slightly Agree on each of items 1, 7, 8, and 10. Score 1 point for Definitely or Slightly
# Disagree on each of items 2, 3, 4, 5, 6, and 9. If the individual scores 6 or above, consider
# referring them for a specialist diagnostic assessment.
clean_data_1 <- clean_data_1 %>%
rename(aq_score_old = aq_score)
clean_data_1$aq_score <-
rowSums(sapply(clean_data_1[c("aq1","aq7","aq8","aq10")], function(x) x %in% c(0,1))) +
rowSums(sapply(clean_data_1[c("aq2","aq3","aq4","aq5","aq6","aq9")], function(x) x %in% c(2,3)))
str(clean_data_1$aq_score)## num [1:7940] 2 2 2 2 2 2 2 2 2 2 ...# Reverse code items 2, 3, 4, 5, 6, and 9.
items_to_invert <- c("aq2", "aq3", "aq4", "aq5", "aq6", "aq9")
for (item in items_to_invert) {
new_col <- paste0(item, "_inv")
clean_data_1[[new_col]] <- 3 - clean_data_1[[item]]
}
# Fix row names/numbers
tail(rownames(clean_data_1))## [1] "8022" "8023" "8024" "8025" "8026" "8027"rownames(clean_data_1) <- NULL
table(clean_data_1$group)/20 # There is something wrong with a participant in the ASD group: it is due to NAs##
## NT ASD
## 200.0 196.9##### Deal with participants with NAs (hence those who skipped questions or took only 1/2 surveys)
# Remove the column remarks, since it obviously has many NAs
clean_data_1$remarks <- NULL
# Identify NA columns per participant
na_summary <- clean_data_1 %>%
group_by(subject_nr) %>%
summarise(
na_columns = list(names(.)[sapply(across(everything()), function(x) any(is.na(x)))]),
.groups = "drop"
) %>%
# Keep only participants with at least one NA
dplyr::filter(lengths(na_columns) > 0) %>%
# Convert the list of NA columns into a single string per participant
dplyr::mutate(na_columns_str = sapply(na_columns, paste, collapse = ", ")) %>%
dplyr::select(subject_nr, na_columns_str)
glimpse(na_summary)## Rows: 80
## Columns: 2
## $ subject_nr <int> 9, 12, 14, 17, 27, 29, 30, 31, 37, 46, 54, 57, 64, 67, …
## $ na_columns_str <chr> "moral_response, letter_response, theory2, lm_rt", "mor…# Most of them simply did not take Part 2 of Study 1 (78/80), but some have other missing values as well (2/80)
na_summary %>%
filter(
sapply(str_split(na_columns_str, ",\\s*"), function(cols) any(!cols %in% c("moral_response", "letter_response", "theory2", "lm_rt")))
) %>%
nrow() # I confirm there are 2 who somehow did not complete the survey## [1] 2# Who are these 2?
na_summary$subject_nr[
sapply(str_split(na_summary$na_columns_str, ",\\s*"),
function(cols) any(!cols %in% c("moral_response", "letter_response", "theory2", "lm_rt")))
]## [1] 391 397# Which rows in clean_data_1 do they correspond to?
which(clean_data_1$subject_nr == "391") # The rows are not continuous, so they quit and rejoined the survey, but did not finish answering## [1] 7801 7802 7803 7804 7805 7806 7807 7808 7809 7810 7811 7812 7813 7814 7815
## [16] 7816 7817 7818 7939 7940which(clean_data_1$subject_nr == "397") # The rows are continuous, so they just skipped questions## [1] 7919 7920 7921 7922 7923 7924 7925 7926 7927 7928 7929 7930 7931 7932 7933
## [16] 7934 7935 7936 7937 7938# I exclude data from those two participants
clean_data_1 <- clean_data_1 %>%
filter(!subject_nr %in% c("391", "397"))
# Rename clean_data_1 to clean_data_1a (containing all people who took Part 1 of Study 1)
clean_data_1a <- clean_data_1
rm(clean_data_1)
length(unique(clean_data_1a$subject_nr)) # 395 participants are left## [1] 395table(clean_data_1a$group)/20 # 200 NT and 195 ASD took Part 1 of Study 1##
## NT ASD
## 200 195# Create a second dataset, named clean_data_1b, containing only participants who took both Part 1 and Part 2
clean_data_1b <- clean_data_1a %>%
filter(!subject_nr %in% na_summary$subject_nr)
length(unique(clean_data_1b$subject_nr)) # 317 participants are left## [1] 317table(clean_data_1b$group)/20 # 164 NT and 153 ASD took Part 1 of Study 1##
## NT ASD
## 164 153##### Compute Cronbach’s alpha of AQ-10 for Study 1
alpha_result <- psych::alpha(clean_data_1a[, c("aq1", "aq2_inv", "aq3_inv", "aq4_inv", "aq5_inv", "aq6_inv",
"aq7", "aq8", "aq9_inv", "aq10")])
alpha_result$total$raw_alpha## [1] 0.7986231##### Import data
raw_data_2 = bind_rows(read.csv('study2a-batch1-asd.csv') %>%
mutate(group = 'ASD'),
read.csv('study2a-batch2-nt2.csv') %>%
mutate(group = 'NT'))
##### Rename variables, keep only relevant columns, and create comorbidity flags
clean_data_2 = raw_data_2 %>%
mutate(relabel_literal = str_detect(stimulus, "<div>\n <p>How"),
relabel_upset = str_detect(stimulus, "<div>\n <p style=\"color: grey;\">How")) %>%
mutate(task = if_else(relabel_literal, 'upset_rating',
if_else(relabel_upset, 'literal_meaning', task)))
clean_data_2 <- clean_data_2 %>%
filter(task %in% c('literal_meaning', 'upset_rating', 'violation') | str_detect(stimulus, 'confident')) %>%
mutate(task = if_else(str_detect(stimulus, 'confident'), 'confidence', task)) %>%
dplyr::select(group, gender, comorbidities, purpose_present, lateralization,
trial_index, rule, condition, run_id, PROLIFIC_PID, task, rt, response,
aq1, aq2, aq3, aq4, aq5, aq6, aq7, aq8, aq9, aq10) %>%
mutate(condition = na_if(condition, "1"),
rule = na_if(rule, ""),
rt = as.numeric(rt)) %>% # Convert "" to NA
fill(condition, .direction = "down") %>%
fill(rule, .direction = "down") %>%
pivot_wider(names_from = 'task', values_from = c('trial_index', 'response', 'rt')) %>%
mutate(response = case_when(lateralization == 'Left0' ~ 1 - as.numeric(response_violation),
lateralization == 'Right1' ~ as.numeric(response_violation)),
response_confidence = as.numeric(response_confidence),
response_literal_meaning = 1 - as.numeric(response_literal_meaning),
response_upset_rating = as.numeric(response_upset_rating),
text = case_when(condition %in% c('purpose_and_text_violate',
'purpose_comply_text_violate') ~ 1,
condition %in% c('purpose_and_text_comply',
'purpose_violate_text_comply') ~ 0),
purpose = case_when(condition %in% c('purpose_and_text_comply',
'purpose_comply_text_violate') ~ 0,
condition %in% c('purpose_and_text_violate',
'purpose_violate_text_comply') ~ 1),
cw_resp = response_confidence * (response - .5) * 2) %>%
group_by(run_id, PROLIFIC_PID) %>%
mutate(trial = dense_rank(trial_index_violation)) %>%
ungroup() %>%
mutate(purpose_display = case_when(purpose_present == 'Block1' & trial <= 16 ~ 1,
purpose_present == 'Block1' & trial > 16 ~ 0,
purpose_present == 'Block2' & trial <= 16 ~ 0,
purpose_present == 'Block2' & trial > 16 ~ 1))
clean_data_2 <- clean_data_2 %>%
rowwise() %>%
mutate(aq_score = sum(c_across(aq1:aq10), na.rm = TRUE)) %>%
ungroup() %>%
mutate(
# Binary comorbidity flags as 0/1 factors
autism = factor(as.integer(str_detect(comorbidities, fixed("ASD (Autism Spectrum Disorder)")))),
adhd = factor(as.integer(str_detect(comorbidities, fixed("Attention Deficit Hyperactivity Disorder")))),
anxiety = factor(as.integer(str_detect(comorbidities, fixed("Anxiety Disorder")))),
depression = factor(as.integer(str_detect(comorbidities, fixed("Depression")))),
ocd = factor(as.integer(str_detect(comorbidities, fixed("Obsessive-Compulsive Disorder")))),
epilepsy = factor(as.integer(str_detect(comorbidities, fixed("Epilepsy or Seizure Disorder")))),
gastrointestinal = factor(as.integer(str_detect(comorbidities, fixed("Gastrointestinal Issues")))),
sleep_dis = factor(as.integer(str_detect(comorbidities, fixed("Sleep Disorders")))),
sensory_process_dis = factor(as.integer(str_detect(comorbidities, fixed("Sensory Processing Disorder")))),
learn_dis = factor(as.integer(str_detect(comorbidities, fixed("Learning Disability (e.g., Dyslexia)")))),
group = factor(group), # NT vs ASD
gender = factor(gender) # ensure categorical
)
clean_data_2 <- clean_data_2 %>%
left_join(
raw_data_2 %>%
group_by(PROLIFIC_PID) %>%
summarize(time_elapsed = max(time_elapsed, na.rm = TRUE)),
by = "PROLIFIC_PID"
)
##### Exclude according to 3 criteria:
##### 1) Fail to complete study (or took it multiple times);
length(unique(clean_data_2$PROLIFIC_PID)) # 417 different people took the survey## [1] 417table(table(clean_data_2$PROLIFIC_PID)) # 384 participants have 32 rows (trials)##
## 1 3 4 8 9 15 16 18 19 27 28 32 33 34 35 36 39 42 44 49
## 5 1 4 1 2 1 1 2 1 1 2 384 1 1 2 1 1 2 1 2
## 64
## 1# I keep only those with 32 trials (to avoid incomplete runs and participants taking the survey multiple times)
clean_data_2 <- clean_data_2[clean_data_2$PROLIFIC_PID %in% names(table(clean_data_2$PROLIFIC_PID)[table(clean_data_2$PROLIFIC_PID) == 32]), ]
any(table(clean_data_2$PROLIFIC_PID) != 32)## [1] FALSElength(unique(clean_data_2$PROLIFIC_PID)) # 384 participants are left## [1] 384##### 2) Complete whole study in < 180 seconds (minimum completion time)
sum(clean_data_2$time_elapsed < 180000) # No participant completed the study in < 180 seconds## [1] 0##### 3) Provide four or more responses in under 500 milliseconds (minimum response time)
# Violation judgment
sum(clean_data_2$rt_violation < 500, na.rm = TRUE) # There are 6 trials with rt < 500 ms for violation judgments## [1] 6# Self-reported confidence
sum(clean_data_2$rt_confidence < 500, na.rm = TRUE) # 0## [1] 0# Literal judgments
sum(clean_data_2$rt_literal_meaning < 500, na.rm = TRUE) # 1## [1] 1# Upsetness judgments
sum(clean_data_2$rt_upset_rating < 500, na.rm = TRUE) # 4## [1] 4# Who are these people and how many trials with rt < 500 ms are they associated with?
clean_data_2 %>%
mutate(fast_count = (rt_violation < 500) +
(rt_confidence < 500) +
(rt_literal_meaning < 500) +
(rt_upset_rating < 500)) %>%
group_by(PROLIFIC_PID) %>%
summarise(n_occurrences = sum(fast_count)) %>%
filter(n_occurrences > 0)## # A tibble: 8 × 2
## PROLIFIC_PID n_occurrences
## <chr> <int>
## 1 5be44162fa676700011d80d7 2
## 2 6128188a420f8d3f27f8d8b1 1
## 3 63d3fa5e6c990bc96e11beec 1
## 4 66ba4608bc582e7db54007e5 3
## 5 671bbaf1bec15d847b871b9b 1
## 6 672e5f7de203449f90bf93a0 1
## 7 6735b2a6112cf17b7b8df8d1 1
## 8 677b2d64aaa09f45f5a53d3e 1# No participant had more than 4 trials with rt < 500 ms, so no exclusions are required
##### Additionally exclude those whose Prolific says they are ASD, but they claim they are not
#sum(clean_data_2$group == "ASD" &
# !str_detect(clean_data_2$comorbidities, fixed("ASD (Autism Spectrum Disorder)")),
# na.rm = TRUE) # 3712 rows (116 participants)
#sum(clean_data_2$group == "ASD" & clean_data_2$autism == 0, na.rm = TRUE) # I confirm 3712 rows (116 participants)
#clean_data_2 <- clean_data_2[!(clean_data_2$group == "ASD" & clean_data_2$autism == 0), ]
#length(unique(clean_data_2$PROLIFIC_PID)) # 268 participants are left
# Make categorical variables factors
str(clean_data_2$text)## num [1:12288] 1 0 0 1 1 0 1 0 1 1 ...clean_data_2$text <- factor(clean_data_2$text)
str(clean_data_2$text)## Factor w/ 2 levels "0","1": 2 1 1 2 2 1 2 1 2 2 ...str(clean_data_2$purpose)## num [1:12288] 1 0 1 0 1 1 0 0 1 0 ...clean_data_2$purpose <- factor(clean_data_2$purpose)
str(clean_data_2$purpose)## Factor w/ 2 levels "0","1": 2 1 2 1 2 2 1 1 2 1 ...str(clean_data_2$purpose_display)## num [1:12288] 1 1 1 1 1 1 1 1 1 1 ...clean_data_2$purpose_display <- factor(clean_data_2$purpose_display)
str(clean_data_2$purpose_display)## Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...str(clean_data_2$response_upset_rating)## num [1:12288] 3 3 3 0 3 3 1 0 3 2 ...clean_data_2$response_upset_rating <- factor(clean_data_2$response_upset_rating)
str(clean_data_2$response_upset_rating)## Factor w/ 4 levels "0","1","2","3": 4 4 4 1 4 4 2 1 4 3 ...class(clean_data_2$response_upset_rating)## [1] "factor"clean_data_2$response_upset_rating <- as.ordered(clean_data_2$response_upset_rating)
class(clean_data_2$response_upset_rating)## [1] "ordered" "factor"# Anonymize participants’ ID
clean_data_2$subject_nr <- match(clean_data_2$PROLIFIC_PID,
unique(clean_data_2$PROLIFIC_PID))
clean_data_2$PROLIFIC_PID <- NULL
# Correct columns’ names
clean_data_2 <- clean_data_2 %>%
rename(scene = rule)
# Create a new column with Upsetness as a factor
clean_data_2$response_upset_factor <- as.factor(clean_data_2$response_upset_rating)
str(clean_data_2$response_upset_factor)## Ord.factor w/ 4 levels "0"<"1"<"2"<"3": 4 4 4 1 4 4 2 1 4 3 ...# Recompute AQ-10 scores
# 0 = Definitely Agree
# 1 = Slightly Agree
# 2 = Slightly Disagree
# 3 = Definitely Disagree
# SCORING: Only 1 point can be scored for each question. Score 1 point for Definitely or
# Slightly Agree on each of items 1, 7, 8, and 10. Score 1 point for Definitely or Slightly
# Disagree on each of items 2, 3, 4, 5, 6, and 9. If the individual scores 6 or above, consider
# referring them for a specialist diagnostic assessment.
clean_data_2 <- clean_data_2 %>%
rename(aq_score_old = aq_score)
clean_data_2$aq_score <-
rowSums(sapply(clean_data_2[c("aq1","aq7","aq8","aq10")], function(x) x %in% c(0,1))) +
rowSums(sapply(clean_data_2[c("aq2","aq3","aq4","aq5","aq6","aq9")], function(x) x %in% c(2,3)))
str(clean_data_2$aq_score)## num [1:12288] 2 2 2 2 2 2 2 2 2 2 ...# Reverse code items 2, 3, 4, 5, 6, and 9.
items_to_invert <- c("aq2", "aq3", "aq4", "aq5", "aq6", "aq9")
for (item in items_to_invert) {
new_col <- paste0(item, "_inv")
clean_data_2[[new_col]] <- 3 - clean_data_2[[item]]
}
table(clean_data_2$group)/20 # There is something wrong with a participant in the ASD group: it is due to NAs##
## ASD NT
## 304.0 310.4##### Deal with participants with NAs (hence those who skipped questions or took only 1/2 surveys)
# Identify NA columns per participant
na_summary <- clean_data_2 %>%
group_by(subject_nr) %>%
summarise(
na_columns = list(names(.)[sapply(across(everything()), function(x) any(is.na(x)))]),
.groups = "drop"
) %>%
# Keep only participants with at least one NA
dplyr::filter(lengths(na_columns) > 0) %>%
# Convert the list of NA columns into a single string per participant
dplyr::mutate(na_columns_str = sapply(na_columns, paste, collapse = ", ")) %>%
dplyr::select(subject_nr, na_columns_str)
glimpse(na_summary)## Rows: 2
## Columns: 2
## $ subject_nr <int> 190, 382
## $ na_columns_str <chr> "aq1, aq2, aq3, aq4, aq5, aq6, aq7, aq8, aq9, aq10, aq_…# 2 people did not complete the survey.
# I exclude data from those two participants
clean_data_2 <- clean_data_2 %>%
filter(!subject_nr %in% c("190", "382"))
##### Compute Cronbach’s alpha of AQ-10 for Study 2
alpha_result <- psych::alpha(clean_data_2[, c("aq1", "aq2_inv", "aq3_inv", "aq4_inv", "aq5_inv", "aq6_inv",
"aq7", "aq8", "aq9_inv", "aq10")])
alpha_result$total$raw_alpha## [1] 0.7331627##### Compute total alpha of AQ-10 for Study 1 and Study 2 together
alpha_columns <- c("aq1", "aq2_inv", "aq3_inv", "aq4_inv", "aq5_inv", "aq6_inv",
"aq7", "aq8", "aq9_inv", "aq10")
cd1_alpha <- clean_data_1a[!duplicated(clean_data_1a$subject_nr), ]
cd1_alpha <- cd1_alpha[, alpha_columns]
cd2_alpha <- clean_data_2[!duplicated(clean_data_2$subject_nr), ]
cd2_alpha <- cd2_alpha[, alpha_columns]
AQ_tot_alpha <- rbind(cd1_alpha, cd2_alpha)
alpha_result <- psych::alpha(AQ_tot_alpha)
alpha_result$total$raw_alpha## [1] 0.7728161### STUDY 1
# Create dataset with 1 row per participant only ASD data (group and aq_score)
aq_data_1 <- clean_data_1a %>%
group_by(group, subject_nr) %>%
summarise(
aq_score = first(aq_score),
gender = first(gender), # keep gender
comorbid = first(comorbidities), # for flag creation
.groups = "drop"
) %>%
mutate(
# Binary comorbidity flags as 0/1 factors
autism = factor(as.integer(str_detect(comorbid, fixed("ASD (Autism Spectrum Disorder)")))),
adhd = factor(as.integer(str_detect(comorbid, fixed("Attention Deficit Hyperactivity Disorder")))),
anxiety = factor(as.integer(str_detect(comorbid, fixed("Anxiety Disorder")))),
depression = factor(as.integer(str_detect(comorbid, fixed("Depression")))),
ocd = factor(as.integer(str_detect(comorbid, fixed("Obsessive-Compulsive Disorder")))),
epilepsy = factor(as.integer(str_detect(comorbid, fixed("Epilepsy or Seizure Disorder")))),
gastrointestinal = factor(as.integer(str_detect(comorbid, fixed("Gastrointestinal Issues")))),
sleep_dis = factor(as.integer(str_detect(comorbid, fixed("Sleep Disorders")))),
sensory_process_dis = factor(as.integer(str_detect(comorbid, fixed("Sensory Processing Disorder")))),
learn_dis = factor(as.integer(str_detect(comorbid, fixed("Learning Disability (e.g., Dyslexia)"))))
) %>%
dplyr::select(-comorbid)
table(aq_data_1$gender)##
## Female Male No answer Non-binary
## 0 156 228 1 10# Compute statistics per group levels (NT vs. ASD)
aq_by_group_1 <- aq_data_1 %>%
group_by(group) %>%
summarise(
n = n(),
mean = mean(aq_score, na.rm = TRUE),
sd = sd(aq_score, na.rm = TRUE),
se = sd / sqrt(n),
ci95_low = mean - qt(.975, df = n - 1) * se,
ci95_high = mean + qt(.975, df = n - 1) * se,
.groups = "drop"
)
aq_by_group_1## # A tibble: 2 × 7
## group n mean sd se ci95_low ci95_high
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NT 200 2.90 1.95 0.138 2.62 3.17
## 2 ASD 195 5.65 2.63 0.188 5.28 6.02### STUDY 2
# Create dataset with 1 row per participant only ASD data (group and aq_score)
aq_data_2 <- clean_data_2 %>%
dplyr::select(group, subject_nr, aq_score, gender, autism, adhd, anxiety, depression,
ocd, epilepsy, gastrointestinal, sleep_dis, sensory_process_dis, learn_dis) %>%
dplyr::group_by(subject_nr) %>%
dplyr::slice(1) %>%
dplyr::ungroup()
table(aq_data_2$gender)##
## Female Male Non-binary
## 1 184 195 2# Compute statistics per group levels (NT vs. ASD)
aq_by_group_2 <- aq_data_2 %>%
group_by(group) %>%
summarise(
n = n(),
mean = mean(aq_score, na.rm = TRUE),
sd = sd(aq_score, na.rm = TRUE),
se = sd / sqrt(n),
ci95_low = mean - qt(.975, df = n - 1) * se,
ci95_high = mean + qt(.975, df = n - 1) * se,
.groups = "drop"
)
aq_by_group_2## # A tibble: 2 × 7
## group n mean sd se ci95_low ci95_high
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ASD 189 4.46 2.38 0.173 4.12 4.80
## 2 NT 193 2.80 1.77 0.128 2.55 3.05### CREATE AGGREGATE DATASET
# Create aggregate aq_data_total
# Find the max subject_nr in aq_data_1
max_id <- max(aq_data_1$subject_nr, na.rm = TRUE)
# Offset subject_nr in aq_data_2
aq_data_2_offset <- aq_data_2 %>%
dplyr::mutate(subject_nr = subject_nr + max_id)
# Combine the datasets
aq_data_total <- dplyr::bind_rows(aq_data_1, aq_data_2_offset)
aq_data_total## # A tibble: 777 × 14
## group subject_nr aq_score gender autism adhd anxiety depression ocd
## <fct> <int> <dbl> <fct> <fct> <fct> <fct> <fct> <fct>
## 1 NT 1 2 Female 0 0 0 0 0
## 2 NT 2 6 Female 0 0 0 0 0
## 3 NT 3 2 Female 0 0 0 0 0
## 4 NT 4 2 Male 0 0 0 0 0
## 5 NT 5 2 Female 0 0 0 0 0
## 6 NT 6 2 Female 0 0 0 0 0
## 7 NT 7 2 Male 0 0 0 0 0
## 8 NT 8 3 Female 0 0 0 0 0
## 9 NT 9 1 Female 0 0 0 1 1
## 10 NT 10 4 Male 0 0 0 0 0
## # ℹ 767 more rows
## # ℹ 5 more variables: epilepsy <fct>, gastrointestinal <fct>, sleep_dis <fct>,
## # sensory_process_dis <fct>, learn_dis <fct># # Create aggregate aq_by_group_total
aq_by_group_total <- aq_data_total %>%
group_by(group) %>%
summarise(
n = n(),
mean = mean(aq_score, na.rm = TRUE),
sd = sd(aq_score, na.rm = TRUE),
se = sd / sqrt(n),
ci95_low = mean - qt(.975, df = n - 1) * se,
ci95_high = mean + qt(.975, df = n - 1) * se,
.groups = "drop"
)
aq_by_group_total## # A tibble: 2 × 7
## group n mean sd se ci95_low ci95_high
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NT 393 2.85 1.86 0.0939 2.67 3.03
## 2 ASD 384 5.07 2.58 0.131 4.81 5.32### STUDY 1
model_1 <- lm(aq_score ~ group + autism, data = aq_data_1)
Anova(model_1, type = 3) # Both main effects## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 1676.21 1 334.741 < 2.2e-16 ***
## group 74.53 1 14.884 0.0001337 ***
## autism 130.16 1 25.993 5.345e-07 ***
## Residuals 1962.93 392
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 2
model_1b <- lm(aq_score ~ group + autism, data = aq_data_2)
Anova(model_1b, type = 3) # Both main effects## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 1592.90 1 400.007 < 2.2e-16 ***
## group 60.81 1 15.271 0.0001103 ***
## autism 160.22 1 40.235 6.421e-10 ***
## Residuals 1509.25 379
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### AGGREGATE (BOTH STUDIES TOGETHER)
model_1c <- lm(aq_score ~ group + autism, data = aq_data_total)
Anova(model_1c, type = 3) # Both main effects## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 3191.9 1 706.563 < 2.2e-16 ***
## group 124.5 1 27.556 1.973e-07 ***
## autism 403.0 1 89.214 < 2.2e-16 ***
## Residuals 3496.5 774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 1
model_2 <- lm(aq_score ~ autism + adhd + anxiety + depression +
ocd + epilepsy + gastrointestinal + sleep_dis +
sensory_process_dis + learn_dis + gender, data = aq_data_1)
Anova(model_2, type = 3) # Effect of self-report autism, and marginal effect of adhd and learning disorder## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 984.46 1 191.7956 < 2e-16 ***
## autism 518.44 1 101.0041 < 2e-16 ***
## adhd 14.92 1 2.9069 0.08902 .
## anxiety 13.61 1 2.6516 0.10427
## depression 0.02 1 0.0044 0.94729
## ocd 5.93 1 1.1557 0.28304
## epilepsy 3.02 1 0.5881 0.44363
## gastrointestinal 0.55 1 0.1069 0.74385
## sleep_dis 0.72 1 0.1399 0.70858
## sensory_process_dis 1.01 1 0.1961 0.65812
## learn_dis 14.77 1 2.8771 0.09066 .
## gender 10.91 3 0.7086 0.54733
## Residuals 1955.62 381
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 2
model_2b <- lm(aq_score ~ autism + adhd + anxiety + depression +
ocd + epilepsy + gastrointestinal + sleep_dis +
sensory_process_dis + learn_dis + gender, data = aq_data_2)
Anova(model_2b, type = 3) # Effect of self-report autism, and marginal effect of anxiety## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 9.00 1 2.1804 0.14063
## autism 259.20 1 62.7969 2.755e-14 ***
## adhd 1.52 1 0.3678 0.54456
## anxiety 11.82 1 2.8633 0.09147 .
## depression 2.75 1 0.6674 0.41448
## ocd 0.50 1 0.1206 0.72857
## epilepsy 1.56 1 0.3780 0.53904
## gastrointestinal 1.73 1 0.4185 0.51809
## sleep_dis 3.11 1 0.7537 0.38587
## sensory_process_dis 1.76 1 0.4258 0.51447
## learn_dis 4.49 1 1.0885 0.29749
## gender 22.80 3 1.8415 0.13918
## Residuals 1518.96 368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### AGGREGATE (BOTH STUDIES TOGETHER)
model_2c <- lm(aq_score ~ autism + adhd + anxiety + depression +
ocd + epilepsy + gastrointestinal + sleep_dis +
sensory_process_dis + learn_dis + gender, data = aq_data_total)
Anova(model_2c, type = 3) # Effect of self-report autism, and significant effect of anxiety## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 9.0 1 1.9415 0.16392
## autism 870.7 1 187.8171 < 2e-16 ***
## adhd 11.0 1 2.3826 0.12310
## anxiety 22.9 1 4.9390 0.02655 *
## depression 1.1 1 0.2364 0.62694
## ocd 6.4 1 1.3806 0.24037
## epilepsy 0.3 1 0.0540 0.81635
## gastrointestinal 2.0 1 0.4273 0.51353
## sleep_dis 5.7 1 1.2267 0.26839
## sensory_process_dis 0.9 1 0.1951 0.65886
## learn_dis 2.3 1 0.4921 0.48320
## gender 23.9 4 1.2914 0.27178
## Residuals 3532.4 762
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 1
model_3 <- glmer(response ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1a, family = 'binomial')
Anova(model_3, type = 3) # Both main effects## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 318.021 1 < 2e-16 ***
## text 769.625 1 < 2e-16 ***
## purpose 623.305 1 < 2e-16 ***
## text:purpose 2.852 1 0.09126 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4 <- glmer(response ~ group * (text * purpose) + (1 | scene) + (1 | subject_nr),
data = clean_data_1a, family = 'binomial')
Anova(model_4, type = 3) # Group moderates both the main effect of Text and of Purpose## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 214.1513 1 < 2e-16 ***
## group 7.1067 1 0.00768 **
## text 334.4301 1 < 2e-16 ***
## purpose 277.7153 1 < 2e-16 ***
## text:purpose 1.0059 1 0.31589
## group:text 5.8448 1 0.01562 *
## group:purpose 6.4138 1 0.01132 *
## group:text:purpose 0.0005 1 0.98218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4e <- glmer(response ~ group * text + (1 | scene) + (1 | subject_nr),
data = clean_data_1a, family = 'binomial')
Anova(model_4e, type = 3) # The moderation disappears## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 79.9361 1 <2e-16 ***
## group 0.5896 1 0.4426
## text 999.4425 1 <2e-16 ***
## group:text 0.7445 1 0.3882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4f <- glmer(response ~ group * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1a, family = 'binomial')
Anova(model_4f, type = 3) # The moderation disappears## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 60.4937 1 7.381e-15 ***
## group 0.9726 1 0.3240
## purpose 704.7650 1 < 2.2e-16 ***
## group:purpose 1.4110 1 0.2349
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_4, model_4e, model_4f) # model_4 with group * (text * purpose) is the best fitting one## df AIC
## model_4 10 6383.559
## model_4e 6 8639.867
## model_4f 6 9439.210### STUDY 2
model_3b <- glmer(response ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_3b, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 524.5024 1 < 2e-16 ***
## text 1916.5411 1 < 2e-16 ***
## purpose 692.0763 1 < 2e-16 ***
## text:purpose 4.1278 1 0.04218 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4b <- glmer(response ~ group * (text * purpose) + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_4b, type = 3) # Moderation on text and on the text * purpose interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 339.2819 1 < 2.2e-16 ***
## group 17.7643 1 2.5e-05 ***
## text 1059.5782 1 < 2.2e-16 ***
## purpose 409.5808 1 < 2.2e-16 ***
## text:purpose 12.6458 1 0.0003764 ***
## group:text 10.5682 1 0.0011505 **
## group:purpose 2.0273 1 0.1544928
## group:text:purpose 6.0188 1 0.0141542 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4g <- glmer(response ~ group * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_4g, type = 3) # Both main effects, no moderation## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 27.7882 1 1.353e-07 ***
## group 5.5152 1 0.01885 *
## purpose 408.3558 1 < 2.2e-16 ***
## group:purpose 0.7557 1 0.38469
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4h <- glmer(response ~ group * text + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_4h, type = 3) # Both main effects, moderation on text## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 142.597 1 < 2.2e-16 ***
## group 37.077 1 1.135e-09 ***
## text 1920.958 1 < 2.2e-16 ***
## group:text 23.138 1 1.508e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_4b, model_4g, model_4h) # model_4b is the best fitting one## df AIC
## model_4b 10 9122.289
## model_4g 6 16043.284
## model_4h 6 10796.794# Replicate the effect in a linear model (with confidence weighting)
### STUDY 1
model_3c <- lmer(cw_resp ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1a)
Anova(model_3c, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 480.200 1 < 2.2e-16 ***
## text 3074.859 1 < 2.2e-16 ***
## purpose 2287.960 1 < 2.2e-16 ***
## text:purpose 30.926 1 2.68e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4c <- lmer(cw_resp ~ group * (text * purpose) + (1 | scene) + (1 | subject_nr),
data = clean_data_1a)
Anova(model_4c, type = 3) # Marginal effect of group, no moderation## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 458.3546 1 < 2.2e-16 ***
## group 3.3359 1 0.0677843 .
## text 1586.9076 1 < 2.2e-16 ***
## purpose 1230.7508 1 < 2.2e-16 ***
## text:purpose 13.9997 1 0.0001828 ***
## group:text 0.2746 1 0.6002440
## group:purpose 2.1782 1 0.1399825
## group:text:purpose 0.0950 1 0.7579329
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4i <- lmer(cw_resp ~ group * text + (1 | scene) + (1 | subject_nr),
data = clean_data_1a)
Anova(model_4i, type = 3) # No moderation## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 119.1045 1 <2e-16 ***
## group 0.8364 1 0.3604
## text 1850.9461 1 <2e-16 ***
## group:text 0.7367 1 0.3907
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4j <- lmer(cw_resp ~ group * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1a)
Anova(model_4j, type = 3) # Marginal moderation on purpose## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 87.0914 1 < 2e-16 ***
## group 2.5528 1 0.11010
## purpose 1253.8331 1 < 2e-16 ***
## group:purpose 3.4189 1 0.06445 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_4c, model_4i, model_4j) # model_4c is the best fitting one## df AIC
## model_4c 11 87188.81
## model_4i 7 90357.98
## model_4j 7 91273.10### STUDY 2
model_3d <- lmer(cw_resp ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_3d, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 1081.914 1 < 2.2e-16 ***
## text 7942.452 1 < 2.2e-16 ***
## purpose 1785.186 1 < 2.2e-16 ***
## text:purpose 92.699 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4d <- lmer(cw_resp ~ group * (text * purpose) + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_4d, type = 3) # Moderation on purpose and on the text * purpose interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 828.6614 1 < 2.2e-16 ***
## group 9.8715 1 0.001679 **
## text 3854.4362 1 < 2.2e-16 ***
## purpose 1062.0086 1 < 2.2e-16 ***
## text:purpose 112.7508 1 < 2.2e-16 ***
## group:text 1.2669 1 0.260350
## group:purpose 15.2468 1 9.433e-05 ***
## group:text:purpose 28.9552 1 7.407e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4k <- lmer(cw_resp ~ group * text + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_4k, type = 3) # Moderation on text## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 324.423 1 < 2.2e-16 ***
## group 49.712 1 1.780e-12 ***
## text 4868.202 1 < 2.2e-16 ***
## group:text 39.733 1 2.912e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_4l <- lmer(cw_resp ~ group * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_4l, type = 3) # No moderation## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 39.1442 1 3.936e-10 ***
## group 6.6840 1 0.009728 **
## purpose 582.8051 1 < 2.2e-16 ***
## group:purpose 0.0092 1 0.923486
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_4d, model_4k, model_4l) # model_4c is slightly better## df AIC
## model_4d 11 134724.2
## model_4k 7 137122.1
## model_4l 7 143884.6### STUDY 1
model_5 <- glmer(response ~ moral_response * letter_response + (1 | scene) + (1 | subject_nr),
data = clean_data_1b, family = 'binomial')
Anova(model_5, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 0.3254 1 0.56836
## moral_response 288.2253 5 < 2e-16 ***
## letter_response 739.2384 5 < 2e-16 ***
## moral_response:letter_response 41.5428 25 0.02011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6 <- glmer(response ~ group * (moral_response * letter_response) + (1 | scene) + (1 | subject_nr),
data = clean_data_1b, family = 'binomial')
Anova(model_6, type = 3) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0001 1 0.9909
## group 0.1581 1 0.6909
## moral_response 194.9494 5 <2e-16 ***
## letter_response 407.7273 5 <2e-16 ***
## moral_response:letter_response 26.0167 25 0.4067
## group:moral_response 6.0413 5 0.3022
## group:letter_response 2.2839 5 0.8086
## group:moral_response:letter_response 17.6383 25 0.8574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6f <- glmer(response ~ group * moral_response + (1 | scene) + (1 | subject_nr),
data = clean_data_1b, family = 'binomial')
Anova(model_6f, type = 3) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 4.9408 1 0.02623 *
## group 0.4035 1 0.52529
## moral_response 849.7932 5 < 2e-16 ***
## group:moral_response 7.7656 5 0.16963
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6g <- glmer(response ~ group * letter_response + (1 | scene) + (1 | subject_nr),
data = clean_data_1b, family = 'binomial')
Anova(model_6g, type = 3) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 0.1039 1 0.7472
## group 0.2199 1 0.6391
## letter_response 1051.7721 5 <2e-16 ***
## group:letter_response 3.9423 5 0.5578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_6, model_6f, model_6g) # model_6 is the best fitting one## df AIC
## model_6 74 4875.164
## model_6f 14 6665.433
## model_6g 14 5905.108### STUDY 2
model_5b <- glmer(response ~ response_literal_meaning + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_5b, type = 3) # Significant main effect## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 963.6 1 < 2.2e-16 ***
## response_literal_meaning 4448.2 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6b <- glmer(response ~ group * response_literal_meaning + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_6b, type = 3) # Significant moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 564.332 1 < 2.2e-16 ***
## group 18.477 1 1.719e-05 ***
## response_literal_meaning 2362.786 1 < 2.2e-16 ***
## group:response_literal_meaning 15.713 1 7.373e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Replicate the effect in a linear model (with confidence weighting)
### STUDY 1
model_5c <- lmer(cw_resp ~ moral_response * letter_response + (1 | scene) + (1 | subject_nr),
data = clean_data_1b)
Anova(model_5c, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 0.071 1 0.7898
## moral_response 866.405 5 < 2.2e-16 ***
## letter_response 2139.504 5 < 2.2e-16 ***
## moral_response:letter_response 126.649 25 1.464e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6c <- lmer(cw_resp ~ group * (moral_response * letter_response) + (1 | scene) + (1 | subject_nr),
data = clean_data_1b)
Anova(model_6c, type = 3) # Moderation on moral_response interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0135 1 0.90742
## group 0.0232 1 0.87897
## moral_response 520.0263 5 < 2.2e-16 ***
## letter_response 1077.9032 5 < 2.2e-16 ***
## moral_response:letter_response 66.9953 25 1.065e-05 ***
## group:moral_response 12.0498 5 0.03411 *
## group:letter_response 5.2123 5 0.39052
## group:moral_response:letter_response 16.8778 25 0.88618
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6h <- lmer(cw_resp ~ group * moral_response + (1 | scene) + (1 | subject_nr),
data = clean_data_1b)
Anova(model_6h, type = 3) # Moderation on moral_response## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 11.1450 1 0.0008426 ***
## group 1.0843 1 0.2977380
## moral_response 2099.9882 5 < 2.2e-16 ***
## group:moral_response 12.9029 5 0.0243058 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6i <- lmer(cw_resp ~ group * letter_response + (1 | scene) + (1 | subject_nr),
data = clean_data_1b)
Anova(model_6i, type = 3) # No moderation## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 0.3148 1 0.5747
## group 0.1195 1 0.7296
## letter_response 2896.0486 5 <2e-16 ***
## group:letter_response 7.3514 5 0.1958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_6c, model_6h, model_6i) # model_6c is the best fitting one## df AIC
## model_6c 75 68866.19
## model_6h 15 71981.76
## model_6i 15 70950.81### STUDY 2
model_5d <- lmer(cw_resp ~ response_literal_meaning + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_5d, type = 3) # Main effect## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 3413.9 1 < 2.2e-16 ***
## response_literal_meaning 23377.9 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6d <- lmer(cw_resp ~ group * response_literal_meaning + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_6d, type = 3) # Significant moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 2128.389 1 < 2.2e-16 ***
## group 20.518 1 5.906e-06 ***
## response_literal_meaning 10757.256 1 < 2.2e-16 ***
## group:response_literal_meaning 25.781 1 3.824e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 2
model_5e <- glmer(response ~ response_literal_meaning * response_upset_rating + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_5e, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 654.049 1 < 2.2e-16 ***
## response_literal_meaning 2701.328 1 < 2.2e-16 ***
## response_upset_rating 408.497 3 < 2.2e-16 ***
## response_literal_meaning:response_upset_rating 52.639 3 2.189e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6e <- glmer(response ~ group * (response_literal_meaning * response_upset_rating) + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_6e, type = 3) # Moderation on response_literal_meaning## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 342.0923 1 < 2.2e-16
## group 25.6991 1 3.990e-07
## response_literal_meaning 1331.8493 1 < 2.2e-16
## response_upset_rating 247.5337 3 < 2.2e-16
## response_literal_meaning:response_upset_rating 29.4850 3 1.771e-06
## group:response_literal_meaning 47.2214 1 6.340e-12
## group:response_upset_rating 1.2485 3 0.7414
## group:response_literal_meaning:response_upset_rating 3.4765 3 0.3238
##
## (Intercept) ***
## group ***
## response_literal_meaning ***
## response_upset_rating ***
## response_literal_meaning:response_upset_rating ***
## group:response_literal_meaning ***
## group:response_upset_rating
## group:response_literal_meaning:response_upset_rating
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6j <- glmer(response ~ group * response_literal_meaning + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_6j, type = 3) # Moderation on response_literal_meaning## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 564.332 1 < 2.2e-16 ***
## group 18.477 1 1.719e-05 ***
## response_literal_meaning 2362.786 1 < 2.2e-16 ***
## group:response_literal_meaning 15.713 1 7.373e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6k <- glmer(response ~ group * response_upset_rating + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_6k, type = 3) # Moderation on response_upset_rating## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 0.7455 1 0.3879
## group 0.4812 1 0.4879
## response_upset_rating 1666.2515 3 < 2.2e-16 ***
## group:response_upset_rating 31.9375 3 5.395e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_6e, model_6j, model_6k) # model_6e is the best fitting one## df AIC
## model_6e 18 6295.133
## model_6j 6 8015.614
## model_6k 10 11419.005### STUDY 1
model_7 <- clmm(letter_response ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1b)
Anova(model_7, type = 3) # Both main effects## Analysis of Deviance Table (Type III tests)
##
## Response: letter_response
## Df Chisq Pr(>Chisq)
## text 1 2283.5546 <2e-16 ***
## purpose 1 427.7998 <2e-16 ***
## text:purpose 1 0.0998 0.752
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 2
model_7b <- glmer(response_literal_meaning ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_7b, type = 3) # Both main effects## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response_literal_meaning
## Chisq Df Pr(>Chisq)
## (Intercept) 374.6022 1 <2e-16 ***
## text 2480.1146 1 <2e-16 ***
## purpose 360.8950 1 <2e-16 ***
## text:purpose 0.0653 1 0.7983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_8 <- clmm(moral_response ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1b)
Anova(model_8, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III tests)
##
## Response: moral_response
## Df Chisq Pr(>Chisq)
## text 1 698.65 < 2.2e-16 ***
## purpose 1 2638.08 < 2.2e-16 ***
## text:purpose 1 127.50 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_8b <- clmm(response_upset_rating ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_8b, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III tests)
##
## Response: response_upset_rating
## Df Chisq Pr(>Chisq)
## text 1 2104.83 < 2.2e-16 ***
## purpose 1 3349.53 < 2.2e-16 ***
## text:purpose 1 132.78 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_8c <- clmm(response_upset_rating ~ purpose_display + (1 | scene) + (1 | subject_nr),
data = clean_data_2)
Anova(model_8c, type = 3) # No effect## Analysis of Deviance Table (Type III tests)
##
## Response: response_upset_rating
## Df Chisq Pr(>Chisq)
## purpose_display 1 1.0652 0.302LOOK AT THE GOOGLE DOC (PART 4 OF ANALYSES)
TO DO
model_9 <- glmer(response ~ group * purpose_display * (text + purpose) + (1 | scene) + (1 | subject_nr),
clean_data_2, family = 'binomial')
Anova(model_9, type = 3) # Main effects of text and purpose, and marginal group:purpose_display:text interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 270.4493 1 < 2.2e-16 ***
## group 15.9362 1 6.551e-05 ***
## purpose_display 0.9929 1 0.3190342
## text 852.6919 1 < 2.2e-16 ***
## purpose 297.7614 1 < 2.2e-16 ***
## group:purpose_display 2.0709 1 0.1501309
## group:text 12.6507 1 0.0003754 ***
## group:purpose 3.3894 1 0.0656185 .
## purpose_display:text 0.9808 1 0.3220092
## purpose_display:purpose 0.9481 1 0.3302027
## group:purpose_display:text 2.0940 1 0.1478806
## group:purpose_display:purpose 1.6444 1 0.1997172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_9b <- glmer(response ~ group * purpose_display * text + (1 | scene) + (1 | subject_nr),
clean_data_2, family = 'binomial')
Anova(model_9b, type = 3) # Main effect of text## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 123.5592 1 < 2.2e-16 ***
## group 18.0237 1 2.182e-05 ***
## purpose_display 0.1087 1 0.741686
## text 1004.1622 1 < 2.2e-16 ***
## group:purpose_display 0.4256 1 0.514148
## group:text 8.2829 1 0.004002 **
## purpose_display:text 0.2611 1 0.609381
## group:purpose_display:text 0.5823 1 0.445423
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_9c <- glmer(response ~ group * purpose_display * purpose + (1 | scene) + (1 | subject_nr),
clean_data_2, family = 'binomial')
Anova(model_9c, type = 3) # Main effect of purpose## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 21.5378 1 3.469e-06 ***
## group 2.1114 1 0.1462
## purpose_display 0.0605 1 0.8057
## purpose 211.5466 1 < 2.2e-16 ***
## group:purpose_display 0.0863 1 0.7689
## group:purpose 0.7029 1 0.4018
## purpose_display:purpose 0.1573 1 0.6917
## group:purpose_display:purpose 0.1005 1 0.7513
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Replicate with self-reported confidence in judgments
model_10 <- lmer(cw_resp ~ group * purpose_display * (text + purpose) + (1 | scene) + (1 | subject_nr),
clean_data_2)
Anova(model_10, type = 3) # Main effects of text and purpose## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 608.8391 1 < 2.2e-16 ***
## group 17.0105 1 3.717e-05 ***
## purpose_display 0.4455 1 0.5045
## text 2991.2386 1 < 2.2e-16 ***
## purpose 636.6484 1 < 2.2e-16 ***
## group:purpose_display 0.7002 1 0.4027
## group:text 20.2563 1 6.773e-06 ***
## group:purpose 0.4510 1 0.5019
## purpose_display:text 0.3245 1 0.5689
## purpose_display:purpose 0.1606 1 0.6886
## group:purpose_display:text 0.3263 1 0.5678
## group:purpose_display:purpose 0.6551 1 0.4183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Upsetness (response_upset_factor) as DV (ordinal regression models)
# model_11 <- polr(response_upset_factor ~ group * purpose_display * text * purpose, clean_data_2)
# model_11_null <- polr(response_upset_factor ~ 1, data = clean_data_2, Hess = TRUE)
# anova(model_null, model_11)
class(clean_data_2$response_upset_factor)## [1] "ordered" "factor"clean_data_2$response_upset_factor <- as.ordered(clean_data_2$response_upset_factor)
class(clean_data_2$response_upset_factor)## [1] "ordered" "factor"model_12 <- clmm(response_upset_factor ~ group * purpose_display * text * purpose + (1 | scene) + (1 | subject_nr),
clean_data_2)
summary(model_12)## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: response_upset_factor ~ group * purpose_display * text * purpose +
## (1 | scene) + (1 | subject_nr)
## data: clean_data_2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 12224 -10797.12 21634.23 5424(16296) 2.83e-02 1.9e+03
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject_nr (Intercept) 0.3846 0.6202
## scene (Intercept) 0.0768 0.2771
## Number of groups: subject_nr 382, scene 8
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## groupNT -0.986889 0.190859 -5.171 2.33e-07
## purpose_display1 -0.004453 0.148974 -0.030 0.9762
## text1 3.218224 0.125938 25.554 < 2e-16
## purpose1 3.755197 0.128507 29.222 < 2e-16
## groupNT:purpose_display1 0.017561 0.253727 0.069 0.9448
## groupNT:text1 0.382711 0.202246 1.892 0.0585
## purpose_display1:text1 -0.129271 0.175216 -0.738 0.4606
## groupNT:purpose1 1.185523 0.204510 5.797 6.76e-09
## purpose_display1:purpose1 0.336264 0.178303 1.886 0.0593
## text1:purpose1 -0.889670 0.175219 -5.077 3.82e-07
## groupNT:purpose_display1:text1 0.075917 0.284982 0.266 0.7899
## groupNT:purpose_display1:purpose1 -0.131995 0.288383 -0.458 0.6472
## groupNT:text1:purpose1 -0.349807 0.268965 -1.301 0.1934
## purpose_display1:text1:purpose1 -0.207997 0.247575 -0.840 0.4008
## groupNT:purpose_display1:text1:purpose1 0.060764 0.380219 0.160 0.8730
##
## groupNT ***
## purpose_display1
## text1 ***
## purpose1 ***
## groupNT:purpose_display1
## groupNT:text1 .
## purpose_display1:text1
## groupNT:purpose1 ***
## purpose_display1:purpose1 .
## text1:purpose1 ***
## groupNT:purpose_display1:text1
## groupNT:purpose_display1:purpose1
## groupNT:text1:purpose1
## purpose_display1:text1:purpose1
## groupNT:purpose_display1:text1:purpose1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 1.7816 0.1507 11.82
## 1|2 2.9301 0.1529 19.16
## 2|3 4.2204 0.1551 27.21exp(coef(model_12))## 0|1 1|2
## 5.9391102 18.7300696
## 2|3 groupNT
## 68.0599140 0.3727344
## purpose_display1 text1
## 0.9955573 24.9836975
## purpose1 groupNT:purpose_display1
## 42.7426338 1.0177165
## groupNT:text1 purpose_display1:text1
## 1.4662542 0.8787356
## groupNT:purpose1 purpose_display1:purpose1
## 3.2723970 1.3997086
## text1:purpose1 groupNT:purpose_display1:text1
## 0.4107913 1.0788727
## groupNT:purpose_display1:purpose1 groupNT:text1:purpose1
## 0.8763452 0.7048239
## purpose_display1:text1:purpose1 groupNT:purpose_display1:text1:purpose1
## 0.8122091 1.0626485#model_12_null <- clmm(response_upset_factor ~ 1 + (1 | scene) + (1 | subject_nr),
# clean_data_2)
#anova(model_12_null, model_12)
#drop1(model_12, test = "Chisq") # 4-way interaction is not significant. Removing it reduces AIC but just negligibly.
#drop1(update(model_12, . ~ . - group:purpose_display:text), test = "Chisq")
# Create nested models for model_12
# 1 - Create formulas for all nested models
formulas <- list(
model_12_full = response_upset_factor ~ group * purpose_display * text * purpose +
(1 | scene) + (1 | subject_nr),
model_12_GrDisTe = response_upset_factor ~ group * purpose_display * text +
(1 | scene) + (1 | subject_nr),
model_12_GrDisTe = response_upset_factor ~ group * purpose_display * purpose +
(1 | scene) + (1 | subject_nr),
model_12_GrTePu = response_upset_factor ~ group * text * purpose +
(1 | scene) + (1 | subject_nr),
model_12_DisTePu = response_upset_factor ~ purpose_display * text * purpose +
(1 | scene) + (1 | subject_nr),
model_12_GrDis = response_upset_factor ~ group * purpose_display +
(1 | scene) + (1 | subject_nr),
model_12_GrTe = response_upset_factor ~ group * text +
(1 | scene) + (1 | subject_nr),
model_12_GrPu = response_upset_factor ~ group * purpose +
(1 | scene) + (1 | subject_nr),
model_12_DisTe = response_upset_factor ~ purpose_display * text +
(1 | scene) + (1 | subject_nr),
model_12_DisPu = response_upset_factor ~ purpose_display * purpose +
(1 | scene) + (1 | subject_nr),
model_12_TePu = response_upset_factor ~ text * purpose +
(1 | scene) + (1 | subject_nr),
model_12_Gr = response_upset_factor ~ group +
(1 | scene) + (1 | subject_nr),
model_12_Dis = response_upset_factor ~ purpose_display +
(1 | scene) + (1 | subject_nr),
model_12_Te = response_upset_factor ~ text +
(1 | scene) + (1 | subject_nr),
model_12_Pu = response_upset_factor ~ purpose +
(1 | scene) + (1 | subject_nr)
)
# 2 - Fit all models and store them in a list
model_12_list <- lapply(formulas, function(f) clmm(f, data = clean_data_2))
# 3 - Assign proper names to the list elements
names(model_12_list) <- names(formulas)
# 4 - Compare models by AIC
aic_table <- sapply(model_12_list, AIC)
aic_table <- sort(aic_table)
aic_table## model_12_full model_12_GrTePu model_12_DisTePu model_12_TePu
## 21634.23 21637.15 21739.05 21748.92
## model_12_GrDisTe model_12_GrPu model_12_DisPu model_12_Pu
## 25955.51 25957.29 26017.51 26021.66
## model_12_GrTe model_12_GrDisTe model_12_DisTe model_12_Te
## 29599.71 29601.03 29605.42 29607.63
## model_12_Gr model_12_GrDis model_12_Dis
## 31513.12 31516.04 31521.46# In sum: Text and Purpose are very important, Group interacts with them,
# while Goal (purpose_display) does not improve the model much.
# 1st best fitting model
model_12_GrTePu <- clmm(response_upset_factor ~ group * text * purpose +
(1 | scene) + (1 | subject_nr), clean_data_2)
summary(model_12_GrTePu)## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: response_upset_factor ~ group * text * purpose + (1 | scene) +
## (1 | subject_nr)
## data: clean_data_2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 12224 -10806.57 21637.15 1790(5386) 3.84e-03 4.1e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject_nr (Intercept) 0.3835 0.6193
## scene (Intercept) 0.0761 0.2759
## Number of groups: subject_nr 382, scene 8
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## groupNT -0.97774 0.14242 -6.865 6.64e-12 ***
## text1 3.15081 0.09075 34.721 < 2e-16 ***
## purpose1 3.91795 0.09417 41.607 < 2e-16 ***
## groupNT:text1 0.42127 0.14322 2.941 0.00327 **
## groupNT:purpose1 1.12187 0.14536 7.718 1.18e-14 ***
## text1:purpose1 -0.99026 0.12480 -7.935 2.10e-15 ***
## groupNT:text1:purpose1 -0.32274 0.19040 -1.695 0.09006 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 1.7838 0.1310 13.61
## 1|2 2.9309 0.1336 21.93
## 2|3 4.2185 0.1361 31.01exp(coef(model_12_GrTePu))## 0|1 1|2 2|3
## 5.9522638 18.7439359 67.9321593
## groupNT text1 purpose1
## 0.3761594 23.3550531 50.2969899
## groupNT:text1 groupNT:purpose1 text1:purpose1
## 1.5238903 3.0705897 0.3714785
## groupNT:text1:purpose1
## 0.7241626# 2nd best fitting model
model_12_DisTePu <- clmm(response_upset_factor ~ purpose_display * text * purpose +
(1 | scene) + (1 | subject_nr), clean_data_2)
summary(model_12_DisTePu)## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: response_upset_factor ~ purpose_display * text * purpose + (1 |
## scene) + (1 | subject_nr)
## data: clean_data_2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 12224 -10857.53 21739.05 2049(6163) 8.89e-03 4.3e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject_nr (Intercept) 0.38873 0.6235
## scene (Intercept) 0.07515 0.2741
## Number of groups: subject_nr 382, scene 8
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## purpose_display1 -0.002304 0.119599 -0.019 0.9846
## text1 3.325290 0.098632 33.714 < 2e-16 ***
## purpose1 4.265558 0.101972 41.831 < 2e-16 ***
## purpose_display1:text1 -0.083814 0.135990 -0.616 0.5377
## purpose_display1:purpose1 0.277257 0.137845 2.011 0.0443 *
## text1:purpose1 -0.989267 0.131517 -7.522 5.4e-14 ***
## purpose_display1:text1:purpose1 -0.185556 0.185160 -1.002 0.3163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 2.2038 0.1328 16.60
## 1|2 3.3406 0.1355 24.66
## 2|3 4.6234 0.1380 33.50exp(coef(model_12_DisTePu))## 0|1 1|2
## 9.0589644 28.2370619
## 2|3 purpose_display1
## 101.8426773 0.9976984
## text1 purpose1
## 27.8070702 71.2046266
## purpose_display1:text1 purpose_display1:purpose1
## 0.9196024 1.3195061
## text1:purpose1 purpose_display1:text1:purpose1
## 0.3718491 0.8306425# Required for dredge
#options(na.action = "na.fail")
# 1. Fit the full model
#model_12_full <- clmm(
# response_upset_factor ~ group * purpose_display * text * purpose +
# (1 | scene) + (1 | subject_nr),
# data = clean_data_2
#)
# 2. Generate all hierarchical submodels
#model_12_set <- dredge(model_12_full)
# 3. Show all models ranked by AIC
#model_12_set
# 4. Save results to object
#model_12_aic_table <- as.data.frame(model_12_set)
# View top models
#head(model_12_aic_table)
# 5. Extract the best model
#model_12_best <- get.models(model_12_set, 1)[[1]]
#summary(model_12_best)
# 6. Extract models within ΔAIC < 2 (often considered equivalent)
#model_12_top_models <- get.models(model_12_set, subset = delta < 2)
#model_12_top_models
# 7. Plot AIC ranking
#plot(model_12_set)
# Literal binary violation judgment (response_literal_meaning) as DV (logistic regression models)
#model_13 <- glmer(response_literal_meaning ~ group * purpose_display * (text + purpose) + (1 | scene) + (1 | subject_nr),
# clean_data_2, family = 'binomial')
#Anova(model_13, type = 3)
#model_13b <- glmer(response_literal_meaning ~ group * (text + purpose) + (1 | scene) + (1 | subject_nr),
# clean_data_2, family = 'binomial')
#Anova(model_13b, type = 3)
#model_13c <- glmer(response_literal_meaning ~ group * text + (1 | scene) + (1 | subject_nr),
# clean_data_2, family = 'binomial')
#Anova(model_13c, type = 3)
#model_13d <- glmer(response_literal_meaning ~ group * purpose + (1 | scene) + (1 | subject_nr),
# clean_data_2, family = 'binomial')
#Anova(model_13d, type = 3)