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;
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)## [1] 1038unique(clean_data_1$PROLIFIC_PID[
clean_data_1$group == "ASD" &
!str_detect(clean_data_1$comorbidities, fixed("ASD (Autism Spectrum Disorder)"))
])## [1] "63d7dc6c0210225d447bf7f3" "6679fce572cc2e7049c02f63"
## [3] "614ce6a7a33bac955957f7a6" "60d0de75282711f66dce1d3c"
## [5] "666af4d6bf41af7b36d832be" "66c61757724d29c2c00376a6"
## [7] "56c893c9dfea6c000db24a52" "5d6ed96c3fe1ac001add3be8"
## [9] "663b3cc669a891b14e3838ed" "66e738df937b233c0ad3427e"
## [11] "66dbf8dd02d9c5ef0a509925" "5fff2877a99bcb65a14a223c"
## [13] "651b43625c3bc6e2217a7e24" "67007c60de5487026b2a9fd2"
## [15] "5f62c1eea3901f34d53e1c18" "66c22291803cf74e9e9df57c"
## [17] "6548218c0c9527cd69c29f55" "63225b65b2d6e889f7289909"
## [19] "66e2ef616d6c87918a06cac8" "5edf975930f30b270a6d90c7"
## [21] "669d37d5e19a4d6c692e49e4" "5dcf0f73af16d3067b87757f"
## [23] "66a3da0f282774c7648ac2e7" "647869e6f69318bacc341d53"
## [25] "66375fb8f0a305603fd38b51" "6091ab55ea2e97910519f8e7"
## [27] "630a8feaa8eef385f72dc1a5" "66697669e9a10915e94dcf9b"
## [29] "6633bd3bc5ec5e52e25592a4" "66c0ff51cdc16a2baa0a44f6"
## [31] "60393d92ed55293a04c77920" "65fd86d8e2090876b355bd0c"
## [33] "6228f5506248f4f80a39008a" "66a66d138294ce022ab541fe"
## [35] "6272eccadd98507271a651c7" "66a68c7589bc09c861687bed"
## [37] "6686a08d1b81f0919733a674" "5e3e964a5ed1320e1757da5a"
## [39] "5c484b16c41b2a000183c837" "654ce755d6963ac49dc27300"
## [41] "6697c424a4294edb51a3b960" "66116329fca98e8641c2483e"
## [43] "65cf6a9d504d43cd1bf99d77" "654c3a40b7939b14ad8ec864"
## [45] "665c89bfc048d14185cac307" "66f59be8dca82c836105816f"
## [47] "5aa990fb4eecca0001de185c" "6153144b6c389cb97e3169c3"
## [49] "5f0151a9e5751c2a5313c45e" "66fafe2f993dfb430b808f37"
## [51] "663945510818036ee5453ff7" "668dacb2b85145da2a254290"
## [53] NAtable(clean_data_1$PROLIFIC_PID[
clean_data_1$group == "ASD" &
!str_detect(clean_data_1$comorbidities, fixed("ASD (Autism Spectrum Disorder)"))
])##
## 56c893c9dfea6c000db24a52 5aa990fb4eecca0001de185c 5c484b16c41b2a000183c837
## 20 20 20
## 5d6ed96c3fe1ac001add3be8 5dcf0f73af16d3067b87757f 5e3e964a5ed1320e1757da5a
## 20 20 20
## 5edf975930f30b270a6d90c7 5f0151a9e5751c2a5313c45e 5f62c1eea3901f34d53e1c18
## 20 18 20
## 5fff2877a99bcb65a14a223c 60393d92ed55293a04c77920 6091ab55ea2e97910519f8e7
## 20 20 20
## 60d0de75282711f66dce1d3c 614ce6a7a33bac955957f7a6 6153144b6c389cb97e3169c3
## 20 20 20
## 6228f5506248f4f80a39008a 6272eccadd98507271a651c7 630a8feaa8eef385f72dc1a5
## 20 20 20
## 63225b65b2d6e889f7289909 63d7dc6c0210225d447bf7f3 647869e6f69318bacc341d53
## 20 20 20
## 651b43625c3bc6e2217a7e24 6548218c0c9527cd69c29f55 654c3a40b7939b14ad8ec864
## 20 20 20
## 654ce755d6963ac49dc27300 65cf6a9d504d43cd1bf99d77 65fd86d8e2090876b355bd0c
## 20 20 20
## 66116329fca98e8641c2483e 6633bd3bc5ec5e52e25592a4 66375fb8f0a305603fd38b51
## 20 20 20
## 663945510818036ee5453ff7 663b3cc669a891b14e3838ed 665c89bfc048d14185cac307
## 20 20 20
## 66697669e9a10915e94dcf9b 666af4d6bf41af7b36d832be 6679fce572cc2e7049c02f63
## 20 20 20
## 6686a08d1b81f0919733a674 668dacb2b85145da2a254290 6697c424a4294edb51a3b960
## 20 20 20
## 669d37d5e19a4d6c692e49e4 66a3da0f282774c7648ac2e7 66a66d138294ce022ab541fe
## 20 20 20
## 66a68c7589bc09c861687bed 66c0ff51cdc16a2baa0a44f6 66c22291803cf74e9e9df57c
## 20 20 20
## 66c61757724d29c2c00376a6 66dbf8dd02d9c5ef0a509925 66e2ef616d6c87918a06cac8
## 20 20 20
## 66e738df937b233c0ad3427e 66f59be8dca82c836105816f 66fafe2f993dfb430b808f37
## 20 20 20
## 67007c60de5487026b2a9fd2
## 20# 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## [1] 345table(clean_data_1$group)/20 # 200 NT and 145 ASD participants##
## ASD NT
## 145 200# Make categorical variables factors
str(clean_data_1$gender)## chr [1:6900] "Female" "Female" "Female" "Female" "Female" "Female" ...clean_data_1$gender <- factor(clean_data_1$gender)
str(clean_data_1$gender)## Factor w/ 4 levels "Female","Male",..: 1 1 1 1 1 1 1 1 1 1 ...str(clean_data_1$group)## chr [1:6900] "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:6900] 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:6900] 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:6900] 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:6900] 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:6900] 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]]
}
# Compute Cronbach’s alpha
alpha_result <- psych::alpha(clean_data_1[, c("aq1", "aq2_inv", "aq3_inv", "aq4_inv", "aq5_inv", "aq6_inv",
"aq7", "aq8", "aq9_inv", "aq10")])
alpha_result$total$raw_alpha## [1] 0.8166325##### 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;
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)## [1] 3712sum(clean_data_2$group == "ASD" & clean_data_2$autism == 0, na.rm = TRUE) # I confirm 3712 rows (116 participants)## [1] 3712clean_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## [1] 268# Make categorical variables factors
str(clean_data_2$text)## num [1:8576] 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:8576] 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:8576] 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:8576] 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:8576] 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]]
}
# Compute Cronbach’s alpha
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.7873506# Compute total alpha of 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_1[!duplicated(clean_data_1$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.8069646### STUDY 1
# Create dataset with 1 row per participant only ASD data (group and aq_score)
aq_data_1 <- clean_data_1 %>%
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
## 143 192 1 9# 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 145 6.13 2.62 0.218 5.70 6.56### 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
## 2 117 147 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 74 5.61 2.85 0.331 4.95 6.27
## 2 NT 194 2.79 1.78 0.128 2.54 3.04### 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: 613 × 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
## # ℹ 603 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 394 2.84 1.86 0.0939 2.66 3.03
## 2 ASD 219 5.95 2.70 0.183 5.59 6.31### STUDY 1
model_1 <- lm(aq_score ~ group, data = aq_data_1)
Anova(model_1, type = 3) # Significant## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 1676.20 1 329.80 < 2.2e-16 ***
## group 880.25 1 173.19 < 2.2e-16 ***
## Residuals 1743.31 343
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 2
model_1b <- lm(aq_score ~ group, data = aq_data_2)
Anova(model_1b, type = 3) # Significant## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 2327.36 1 515.054 < 2.2e-16 ***
## group 425.82 1 94.236 < 2.2e-16 ***
## Residuals 1201.97 266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### AGGREGATE (BOTH STUDIES TOGETHER)
model_1c <- lm(aq_score ~ group, data = aq_data_total)
Anova(model_1c, type = 3) # Significant## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 3183.8 1 657.23 < 2.2e-16 ***
## group 1362.9 1 281.36 < 2.2e-16 ***
## Residuals 2959.8 611
## ---
## 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 learning disorder, and marginal effect of anxiety and sleep disorder## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 813.36 1 162.1983 < 2e-16 ***
## autism 566.06 1 112.8821 < 2e-16 ***
## adhd 0.06 1 0.0115 0.91473
## anxiety 15.60 1 3.1102 0.07873 .
## depression 0.55 1 0.1106 0.73972
## ocd 8.85 1 1.7654 0.18487
## epilepsy 5.30 1 1.0575 0.30454
## gastrointestinal 0.68 1 0.1362 0.71233
## sleep_dis 16.64 1 3.3193 0.06938 .
## sensory_process_dis 0.14 1 0.0276 0.86817
## learn_dis 19.97 1 3.9820 0.04681 *
## gender 14.18 3 0.9428 0.42013
## Residuals 1659.84 331
## ---
## 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) # No other effects## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 4.50 1 0.9937 0.3198
## autism 306.72 1 67.7279 9.866e-15 ***
## adhd 4.17 1 0.9218 0.3379
## anxiety 0.03 1 0.0060 0.9383
## depression 1.33 1 0.2936 0.5884
## ocd 0.00 1 0.0000 0.9995
## epilepsy 0.13 1 0.0283 0.8665
## gastrointestinal 0.20 1 0.0448 0.8325
## sleep_dis 8.78 1 1.9391 0.1650
## sensory_process_dis 0.07 1 0.0163 0.8985
## learn_dis 0.04 1 0.0096 0.9219
## gender 25.57 3 1.8822 0.1330
## Residuals 1150.29 254
## ---
## 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) # Significant effect of sleep disorder and marginal effect of learning disorder## Anova Table (Type III tests)
##
## Response: aq_score
## Sum Sq Df F value Pr(>F)
## (Intercept) 1453.62 1 304.5782 < 2e-16 ***
## autism 911.03 1 190.8897 < 2e-16 ***
## adhd 1.86 1 0.3900 0.53251
## anxiety 8.56 1 1.7927 0.18111
## depression 0.13 1 0.0275 0.86828
## ocd 5.47 1 1.1464 0.28474
## epilepsy 6.35 1 1.3300 0.24926
## gastrointestinal 0.06 1 0.0130 0.90932
## sleep_dis 30.27 1 6.3417 0.01205 *
## sensory_process_dis 0.35 1 0.0726 0.78763
## learn_dis 14.13 1 2.9603 0.08585 .
## gender 33.36 4 1.7474 0.13798
## Residuals 2854.00 598
## ---
## 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_1, 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) 302.3087 1 <2e-16 ***
## text 624.8260 1 <2e-16 ***
## purpose 515.9398 1 <2e-16 ***
## text:purpose 2.3694 1 0.1237
## ---
## 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_1, family = 'binomial')
Anova(model_4, type = 3) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 213.3216 1 <2e-16 ***
## group 2.3341 1 0.1266
## text 333.1003 1 <2e-16 ***
## purpose 276.7290 1 <2e-16 ***
## text:purpose 1.0005 1 0.3172
## group:text 1.9249 1 0.1653
## group:purpose 1.8824 1 0.1701
## group:text:purpose 0.0118 1 0.9136
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### 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## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 468.2058 1 <2e-16 ***
## text 1187.8163 1 <2e-16 ***
## purpose 407.0032 1 <2e-16 ***
## text:purpose 0.1621 1 0.6872
## ---
## 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) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 187.0741 1 <2e-16 ***
## group 0.6320 1 0.4266
## text 369.3217 1 <2e-16 ***
## purpose 123.3891 1 <2e-16 ***
## text:purpose 1.9261 1 0.1652
## group:text 0.1715 1 0.6788
## group:purpose 0.1003 1 0.7515
## group:text:purpose 1.9096 1 0.1670
## ---
## 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_3c <- lmer(cw_resp ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1)
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) 497.502 1 < 2.2e-16 ***
## text 2804.793 1 < 2.2e-16 ***
## purpose 2150.041 1 < 2.2e-16 ***
## text:purpose 29.185 1 6.578e-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_1)
Anova(model_4c, type = 3) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 471.7002 1 < 2.2e-16 ***
## group 1.3594 1 0.2436499
## text 1636.7627 1 < 2.2e-16 ***
## purpose 1269.4167 1 < 2.2e-16 ***
## text:purpose 14.4395 1 0.0001447 ***
## group:text 0.0439 1 0.8339657
## group:purpose 0.2536 1 0.6145808
## group:text:purpose 0.2333 1 0.6291235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### 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) 1208.04 1 < 2.2e-16 ***
## text 6091.91 1 < 2.2e-16 ***
## purpose 1138.45 1 < 2.2e-16 ***
## text:purpose 23.68 1 1.138e-06 ***
## ---
## 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) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 687.5457 1 < 2.2e-16 ***
## group 0.2278 1 0.6331917
## text 1702.4099 1 < 2.2e-16 ***
## purpose 338.0830 1 < 2.2e-16 ***
## text:purpose 14.6226 1 0.0001313 ***
## group:text 0.0833 1 0.7728555
## group:purpose 0.5955 1 0.4402859
## group:text:purpose 2.2171 1 0.1364909
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 1
model_5 <- glmer(response ~ moral_response * letter_response + (1 | scene) + (1 | subject_nr),
data = clean_data_1, family = 'binomial')
Anova(model_5, type = 3) # Both main effects, no interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0431 1 0.8355
## moral_response 276.0415 5 <2e-16 ***
## letter_response 666.9376 5 <2e-16 ***
## moral_response:letter_response 26.9181 25 0.3600
## ---
## 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_1, 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.9929
## group 0.0003 1 0.9857
## moral_response 197.1834 5 <2e-16 ***
## letter_response 408.9429 5 <2e-16 ***
## moral_response:letter_response 26.2777 25 0.3929
## group:moral_response 1.9023 5 0.8625
## group:letter_response 2.3708 5 0.7958
## group:moral_response:letter_response 14.6738 25 0.9487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### 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) 660.85 1 < 2.2e-16 ***
## response_literal_meaning 3076.47 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) 279.0761 1 < 2e-16 ***
## group 4.2906 1 0.03832 *
## response_literal_meaning 967.6552 1 < 2e-16 ***
## group:response_literal_meaning 6.6058 1 0.01016 *
## ---
## 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_1)
Anova(model_5c, type = 3) # Both main effects, no interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0375 1 0.8465
## moral_response 849.1193 5 < 2.2e-16 ***
## letter_response 1954.1410 5 < 2.2e-16 ***
## moral_response:letter_response 104.2957 25 1.174e-11 ***
## ---
## 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_1)
Anova(model_6c, type = 3) # Marginal group * letter_response interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cw_resp
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0098 1 0.9210
## group 0.3121 1 0.5764
## moral_response 544.0282 5 < 2.2e-16 ***
## letter_response 1125.2261 5 < 2.2e-16 ***
## moral_response:letter_response 70.6920 25 3.035e-06 ***
## group:moral_response 7.2345 5 0.2038
## group:letter_response 3.2997 5 0.6539
## group:moral_response:letter_response 21.4951 25 0.6647
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### 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) 2410.7 1 < 2.2e-16 ***
## response_literal_meaning 17711.3 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) 1135.9693 1 < 2e-16 ***
## group 3.1003 1 0.07828 .
## response_literal_meaning 4614.2034 1 < 2e-16 ***
## group:response_literal_meaning 5.7630 1 0.01637 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 2
model_5c <- glmer(response ~ response_literal_meaning * response_upset_rating + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_5c, type = 3) # Both main effects and interaction## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 531.322 1 < 2.2e-16 ***
## response_literal_meaning 1841.824 1 < 2.2e-16 ***
## response_upset_rating 254.757 3 < 2.2e-16 ***
## response_literal_meaning:response_upset_rating 30.258 3 1.218e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6c <- glmer(response ~ group * (response_literal_meaning * response_upset_rating) + (1 | scene) + (1 | subject_nr),
data = clean_data_2, family = 'binomial')
Anova(model_6c, type = 3) # No moderation of group## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: response
## Chisq Df Pr(>Chisq)
## (Intercept) 188.1589 1 < 2.2e-16 ***
## group 3.1941 1 0.073902 .
## response_literal_meaning 548.2422 1 < 2.2e-16 ***
## response_upset_rating 80.2023 3 < 2.2e-16 ***
## response_literal_meaning:response_upset_rating 10.5757 3 0.014256 *
## group:response_literal_meaning 10.4502 1 0.001226 **
## group:response_upset_rating 1.8991 3 0.593600
## group:response_literal_meaning:response_upset_rating 1.8194 3 0.610724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### STUDY 1
model_7 <- clmm(letter_response ~ text * purpose + (1 | scene) + (1 | subject_nr),
data = clean_data_1)
Anova(model_7, type = 3) # Both main effects## Analysis of Deviance Table (Type III tests)
##
## Response: letter_response
## Df Chisq Pr(>Chisq)
## text 1 2038.8168 <2e-16 ***
## purpose 1 382.6831 <2e-16 ***
## text:purpose 1 1.6402 0.2003
## ---
## 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) 373.246 1 <2e-16 ***
## text 1667.245 1 <2e-16 ***
## purpose 210.730 1 <2e-16 ***
## text:purpose 1.537 1 0.2151
## ---
## 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_1)
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 652.66 < 2.2e-16 ***
## purpose 1 2390.42 < 2.2e-16 ***
## text:purpose 1 131.31 < 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 1378.51 < 2.2e-16 ***
## purpose 1 2427.78 < 2.2e-16 ***
## text:purpose 1 107.84 < 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 0.3464 0.5562LOOK 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) 148.8489 1 < 2e-16 ***
## group 0.3223 1 0.57023
## purpose_display 0.5877 1 0.44332
## text 321.5305 1 < 2e-16 ***
## purpose 97.0124 1 < 2e-16 ***
## group:purpose_display 1.6065 1 0.20498
## group:text 0.0000 1 0.99762
## group:purpose 0.1557 1 0.69311
## purpose_display:text 1.9648 1 0.16100
## purpose_display:purpose 0.4539 1 0.50049
## group:purpose_display:text 3.1144 1 0.07761 .
## group:purpose_display:purpose 1.1043 1 0.29333
## ---
## 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) 118.9733 1 <2e-16 ***
## group 0.1707 1 0.6795
## purpose_display 0.1310 1 0.7174
## text 459.3778 1 <2e-16 ***
## group:purpose_display 0.4547 1 0.5001
## group:text 0.1717 1 0.6786
## purpose_display:text 1.5687 1 0.2104
## group:purpose_display:text 1.8907 1 0.1691
## ---
## 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) 14.9433 1 0.0001108 ***
## group 0.5182 1 0.4716189
## purpose_display 0.2698 1 0.6034724
## purpose 66.2821 1 3.908e-16 ***
## group:purpose_display 0.1095 1 0.7407144
## group:purpose 0.0898 1 0.7643987
## purpose_display:purpose 0.0574 1 0.8107158
## group:purpose_display:purpose 0.0372 1 0.8470814
## ---
## 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) 498.7629 1 <2e-16 ***
## group 0.6894 1 0.4064
## purpose_display 0.0163 1 0.8985
## text 1532.5463 1 <2e-16 ***
## purpose 237.4549 1 <2e-16 ***
## group:purpose_display 0.0435 1 0.8347
## group:text 0.0095 1 0.9224
## group:purpose 0.0504 1 0.8224
## purpose_display:text 0.8599 1 0.3538
## purpose_display:purpose 0.1224 1 0.7265
## group:purpose_display:text 0.8859 1 0.3466
## group:purpose_display:purpose 0.0063 1 0.9365
## ---
## 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 8576 -7065.31 14170.62 3245(9738) 3.01e-03 3.0e+03
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject_nr (Intercept) 0.4556 0.6750
## scene (Intercept) 0.1043 0.3229
## Number of groups: subject_nr 268, scene 8
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## groupNT -0.51641 0.26468 -1.951 0.051042
## purpose_display1 -0.15874 0.28960 -0.548 0.583589
## text1 3.52139 0.22840 15.418 < 2e-16
## purpose1 4.87535 0.23544 20.708 < 2e-16
## groupNT:purpose_display1 0.18238 0.35490 0.514 0.607326
## groupNT:text1 0.16229 0.27757 0.585 0.558769
## purpose_display1:text1 -0.06288 0.32511 -0.193 0.846643
## groupNT:purpose1 0.19891 0.28254 0.704 0.481426
## purpose_display1:purpose1 0.50965 0.33402 1.526 0.127057
## text1:purpose1 -1.12131 0.32329 -3.468 0.000524
## groupNT:purpose_display1:text1 -0.01284 0.39517 -0.032 0.974084
## groupNT:purpose_display1:purpose1 -0.33204 0.40366 -0.823 0.410742
## groupNT:text1:purpose1 -0.16047 0.38247 -0.420 0.674803
## purpose_display1:text1:purpose1 -0.04697 0.47063 -0.100 0.920500
## groupNT:purpose_display1:text1:purpose1 -0.06998 0.55223 -0.127 0.899162
##
## 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 2.2808 0.2436 9.362
## 1|2 3.5056 0.2462 14.241
## 2|3 4.8316 0.2484 19.450exp(coef(model_12))## 0|1 1|2
## 9.7848783 33.3024973
## 2|3 groupNT
## 125.4134135 0.5966560
## purpose_display1 text1
## 0.8532146 33.8315565
## purpose1 groupNT:purpose_display1
## 131.0204128 1.2000701
## groupNT:text1 purpose_display1:text1
## 1.1761965 0.9390583
## groupNT:purpose1 purpose_display1:purpose1
## 1.2200732 1.6647046
## text1:purpose1 groupNT:purpose_display1:text1
## 0.3258525 0.9872444
## groupNT:purpose_display1:purpose1 groupNT:text1:purpose1
## 0.7174564 0.8517413
## purpose_display1:text1:purpose1 groupNT:purpose_display1:text1:purpose1
## 0.9541159 0.9324142#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_GrTePu model_12_DisTePu model_12_TePu model_12_full
## 14166.99 14168.38 14170.37 14170.62
## model_12_GrPu model_12_GrDisTe model_12_DisPu model_12_Pu
## 17260.16 17261.78 17268.14 17269.04
## model_12_GrTe model_12_Te model_12_GrDisTe model_12_DisTe
## 20736.43 20740.26 20741.28 20741.59
## model_12_Gr model_12_GrDis model_12_Dis
## 21903.00 21906.53 21911.77# 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 8576 -7071.50 14166.99 1910(5733) 1.35e-03 6.9e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject_nr (Intercept) 0.4543 0.6740
## scene (Intercept) 0.1039 0.3223
## Number of groups: subject_nr 268, scene 8
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## groupNT -0.42695 0.20067 -2.128 0.0334 *
## text1 3.48313 0.16591 20.994 < 2e-16 ***
## purpose1 5.12080 0.17448 29.349 < 2e-16 ***
## groupNT:text1 0.16048 0.19870 0.808 0.4193
## groupNT:purpose1 0.03924 0.20365 0.193 0.8472
## text1:purpose1 -1.14224 0.23596 -4.841 1.29e-06 ***
## groupNT:text1:purpose1 -0.19789 0.27611 -0.717 0.4736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 2.3580 0.2010 11.73
## 1|2 3.5812 0.2041 17.55
## 2|3 4.9045 0.2068 23.71exp(coef(model_12_GrTePu))## 0|1 1|2 2|3
## 10.5695353 35.9161683 134.8926452
## groupNT text1 purpose1
## 0.6524952 32.5613370 167.4686875
## groupNT:text1 groupNT:purpose1 text1:purpose1
## 1.1740775 1.0400213 0.3191043
## groupNT:text1:purpose1
## 0.8204608# 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 8576 -7072.19 14168.38 1336(4011) 1.25e-02 4.8e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject_nr (Intercept) 0.4794 0.6924
## scene (Intercept) 0.1044 0.3231
## Number of groups: subject_nr 268, scene 8
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## purpose_display1 -0.03714 0.16732 -0.222 0.824
## text1 3.62937 0.13345 27.197 < 2e-16 ***
## purpose1 5.01026 0.13906 36.028 < 2e-16 ***
## purpose_display1:text1 -0.06164 0.18451 -0.334 0.738
## purpose_display1:purpose1 0.27519 0.18696 1.472 0.141
## text1:purpose1 -1.22669 0.17289 -7.095 1.29e-12 ***
## purpose_display1:text1:purpose1 -0.11881 0.24418 -0.487 0.627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 2.6449 0.1702 15.54
## 1|2 3.8706 0.1740 22.24
## 2|3 5.1962 0.1774 29.29exp(coef(model_12_DisTePu))## 0|1 1|2
## 14.0824291 47.9707002
## 2|3 purpose_display1
## 180.5842609 0.9635378
## text1 purpose1
## 37.6889817 149.9442153
## purpose_display1:text1 purpose_display1:purpose1
## 0.9402234 1.3167772
## text1:purpose1 purpose_display1:text1:purpose1
## 0.2932624 0.8879800# 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)