Similarity Ratings Analyses
Libraries and Data Files
#load libraries
library(dplyr)
library(ggplot2)
library(emmeans)
library(lme4)
library(lmerTest)
#load data
data<-read.csv("AllData_SimilarityRatings.csv", header = T)
Data Pre-processing
#rename useful variables
data$sbj<-data$participant
data$resp<-data$slider_exp.response
data$block<-data$blocks.thisN +1
data <- mutate_if(data, is.character, as.factor)
data$sbj<-as.factor(data$sbj)
str(data)
## 'data.frame': 9648 obs. of 38 variables:
## $ trl : int 20 12 17 14 21 3 10 13 6 23 ...
## $ cnd : Factor w/ 2 levels "between","within": 2 2 2 2 1 2 2 2 2 1 ...
## $ dim : Factor w/ 2 levels "irrel","rel": 2 1 2 2 2 1 1 2 1 2 ...
## $ size1 : num 3.4 2.6 2.2 2.6 2.2 2.2 2.6 2.2 3 3 ...
## $ hue1 : num 0.1 0.3 0.1 -0.3 -0.1 0.1 -0.1 -0.3 -0.1 -0.1 ...
## $ size2 : num 3.4 3 2.2 2.6 2.2 2.6 3 2.2 3.4 3 ...
## $ hue2 : num 0.3 0.3 0.3 -0.1 0.1 0.1 -0.1 -0.1 -0.1 0.1 ...
## $ practice_trials.thisRepN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisTrialN: logi NA NA NA NA NA NA ...
## $ practice_trials.thisN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisIndex : logi NA NA NA NA NA NA ...
## $ blocks.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisTrialN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisIndex : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisTrialN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisIndex : int 19 11 16 13 20 2 9 12 5 22 ...
## $ break_loop.thisRepN : logi NA NA NA NA NA NA ...
## $ break_loop.thisTrialN : logi NA NA NA NA NA NA ...
## $ break_loop.thisN : logi NA NA NA NA NA NA ...
## $ break_loop.thisIndex : logi NA NA NA NA NA NA ...
## $ thisRow.t : num 68.9 73.1 75.9 79.7 82.8 ...
## $ notes : logi NA NA NA NA NA NA ...
## $ slider_exp.response : int 5 7 5 4 3 4 3 2 7 5 ...
## $ slider_exp.rt : num 3.89 2.6 3.45 2.87 4.87 ...
## $ participant : int 17202 17202 17202 17202 17202 17202 17202 17202 17202 17202 ...
## $ session : Factor w/ 2 levels "Post","Pre": 1 1 1 1 1 1 1 1 1 1 ...
## $ group : Factor w/ 4 levels "1100","1600",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ date : Factor w/ 134 levels "2025-02-25_11h40.18.454",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ expName : Factor w/ 1 level "SimilarityRatings": 1 1 1 1 1 1 1 1 1 1 ...
## $ psychopyVersion : Factor w/ 1 level "2023.2.3": 1 1 1 1 1 1 1 1 1 1 ...
## $ frameRate : num 60 60 60 60 60 ...
## $ expStart : Factor w/ 134 levels "2025-02-25 11h48.02.736634 +0200",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ sbj : Factor w/ 67 levels "341","1003","1206",..: 10 10 10 10 10 10 10 10 10 10 ...
## $ resp : int 5 7 5 4 3 4 3 2 7 5 ...
## $ block : num 1 1 1 1 1 1 1 1 1 1 ...
#re-order the levels of the session factor
data$session <- factor(data$session, levels = c("Pre", "Post"))
#exclude one participant (20335) for not following instructions
data<-droplevels(data[data$sbj!="20335",])
str(data)
## 'data.frame': 9504 obs. of 38 variables:
## $ trl : int 20 12 17 14 21 3 10 13 6 23 ...
## $ cnd : Factor w/ 2 levels "between","within": 2 2 2 2 1 2 2 2 2 1 ...
## $ dim : Factor w/ 2 levels "irrel","rel": 2 1 2 2 2 1 1 2 1 2 ...
## $ size1 : num 3.4 2.6 2.2 2.6 2.2 2.2 2.6 2.2 3 3 ...
## $ hue1 : num 0.1 0.3 0.1 -0.3 -0.1 0.1 -0.1 -0.3 -0.1 -0.1 ...
## $ size2 : num 3.4 3 2.2 2.6 2.2 2.6 3 2.2 3.4 3 ...
## $ hue2 : num 0.3 0.3 0.3 -0.1 0.1 0.1 -0.1 -0.1 -0.1 0.1 ...
## $ practice_trials.thisRepN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisTrialN: logi NA NA NA NA NA NA ...
## $ practice_trials.thisN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisIndex : logi NA NA NA NA NA NA ...
## $ blocks.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisTrialN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisIndex : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisTrialN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisIndex : int 19 11 16 13 20 2 9 12 5 22 ...
## $ break_loop.thisRepN : logi NA NA NA NA NA NA ...
## $ break_loop.thisTrialN : logi NA NA NA NA NA NA ...
## $ break_loop.thisN : logi NA NA NA NA NA NA ...
## $ break_loop.thisIndex : logi NA NA NA NA NA NA ...
## $ thisRow.t : num 68.9 73.1 75.9 79.7 82.8 ...
## $ notes : logi NA NA NA NA NA NA ...
## $ slider_exp.response : int 5 7 5 4 3 4 3 2 7 5 ...
## $ slider_exp.rt : num 3.89 2.6 3.45 2.87 4.87 ...
## $ participant : int 17202 17202 17202 17202 17202 17202 17202 17202 17202 17202 ...
## $ session : Factor w/ 2 levels "Pre","Post": 2 2 2 2 2 2 2 2 2 2 ...
## $ group : Factor w/ 4 levels "1100","1600",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ date : Factor w/ 132 levels "2025-02-25_11h40.18.454",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ expName : Factor w/ 1 level "SimilarityRatings": 1 1 1 1 1 1 1 1 1 1 ...
## $ psychopyVersion : Factor w/ 1 level "2023.2.3": 1 1 1 1 1 1 1 1 1 1 ...
## $ frameRate : num 60 60 60 60 60 ...
## $ expStart : Factor w/ 132 levels "2025-02-25 11h48.02.736634 +0200",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ sbj : Factor w/ 66 levels "341","1003","1206",..: 10 10 10 10 10 10 10 10 10 10 ...
## $ resp : int 5 7 5 4 3 4 3 2 7 5 ...
## $ block : num 1 1 1 1 1 1 1 1 1 1 ...
#retain the full sample for sensitivity analysis
data_full <- data
#exclude five participants ("43423" "48888" "63018" "77337" "77623") for not meeting the learning criterion
data<-droplevels(data[!(data$sbj %in% c("43423","48888", "63018", "77337", "77623")), ])
str(data)
## 'data.frame': 8784 obs. of 38 variables:
## $ trl : int 20 12 17 14 21 3 10 13 6 23 ...
## $ cnd : Factor w/ 2 levels "between","within": 2 2 2 2 1 2 2 2 2 1 ...
## $ dim : Factor w/ 2 levels "irrel","rel": 2 1 2 2 2 1 1 2 1 2 ...
## $ size1 : num 3.4 2.6 2.2 2.6 2.2 2.2 2.6 2.2 3 3 ...
## $ hue1 : num 0.1 0.3 0.1 -0.3 -0.1 0.1 -0.1 -0.3 -0.1 -0.1 ...
## $ size2 : num 3.4 3 2.2 2.6 2.2 2.6 3 2.2 3.4 3 ...
## $ hue2 : num 0.3 0.3 0.3 -0.1 0.1 0.1 -0.1 -0.1 -0.1 0.1 ...
## $ practice_trials.thisRepN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisTrialN: logi NA NA NA NA NA NA ...
## $ practice_trials.thisN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisIndex : logi NA NA NA NA NA NA ...
## $ blocks.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisTrialN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisIndex : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisTrialN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisIndex : int 19 11 16 13 20 2 9 12 5 22 ...
## $ break_loop.thisRepN : logi NA NA NA NA NA NA ...
## $ break_loop.thisTrialN : logi NA NA NA NA NA NA ...
## $ break_loop.thisN : logi NA NA NA NA NA NA ...
## $ break_loop.thisIndex : logi NA NA NA NA NA NA ...
## $ thisRow.t : num 68.9 73.1 75.9 79.7 82.8 ...
## $ notes : logi NA NA NA NA NA NA ...
## $ slider_exp.response : int 5 7 5 4 3 4 3 2 7 5 ...
## $ slider_exp.rt : num 3.89 2.6 3.45 2.87 4.87 ...
## $ participant : int 17202 17202 17202 17202 17202 17202 17202 17202 17202 17202 ...
## $ session : Factor w/ 2 levels "Pre","Post": 2 2 2 2 2 2 2 2 2 2 ...
## $ group : Factor w/ 4 levels "1100","1600",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ date : Factor w/ 122 levels "2025-02-25_13h08.44.108",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ expName : Factor w/ 1 level "SimilarityRatings": 1 1 1 1 1 1 1 1 1 1 ...
## $ psychopyVersion : Factor w/ 1 level "2023.2.3": 1 1 1 1 1 1 1 1 1 1 ...
## $ frameRate : num 60 60 60 60 60 ...
## $ expStart : Factor w/ 122 levels "2025-02-25 13h09.40.421447 +0200",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ sbj : Factor w/ 61 levels "341","1003","1206",..: 10 10 10 10 10 10 10 10 10 10 ...
## $ resp : int 5 7 5 4 3 4 3 2 7 5 ...
## $ block : num 1 1 1 1 1 1 1 1 1 1 ...
Comparison of Responses from “Exclusion-Instructions”
Participants
#check if we can keep data from 5 first sbjs
data$excl_instr<-"no"
#For the first 5 participants, we need to register that the instructions were different.
#These 5 participants are: "17202","57393", "94372", "40242", "89443"
data[data$sbj=="17202" | data$sbj=="57393" | data$sbj=="94372" | data$sbj=="40242" | data$sbj=="89443",]$excl_instr<-"yes"
data$excl_instr<-as.factor(data$excl_instr)
str(data)
## 'data.frame': 8784 obs. of 39 variables:
## $ trl : int 20 12 17 14 21 3 10 13 6 23 ...
## $ cnd : Factor w/ 2 levels "between","within": 2 2 2 2 1 2 2 2 2 1 ...
## $ dim : Factor w/ 2 levels "irrel","rel": 2 1 2 2 2 1 1 2 1 2 ...
## $ size1 : num 3.4 2.6 2.2 2.6 2.2 2.2 2.6 2.2 3 3 ...
## $ hue1 : num 0.1 0.3 0.1 -0.3 -0.1 0.1 -0.1 -0.3 -0.1 -0.1 ...
## $ size2 : num 3.4 3 2.2 2.6 2.2 2.6 3 2.2 3.4 3 ...
## $ hue2 : num 0.3 0.3 0.3 -0.1 0.1 0.1 -0.1 -0.1 -0.1 0.1 ...
## $ practice_trials.thisRepN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisTrialN: logi NA NA NA NA NA NA ...
## $ practice_trials.thisN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisIndex : logi NA NA NA NA NA NA ...
## $ blocks.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisTrialN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisIndex : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisTrialN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisIndex : int 19 11 16 13 20 2 9 12 5 22 ...
## $ break_loop.thisRepN : logi NA NA NA NA NA NA ...
## $ break_loop.thisTrialN : logi NA NA NA NA NA NA ...
## $ break_loop.thisN : logi NA NA NA NA NA NA ...
## $ break_loop.thisIndex : logi NA NA NA NA NA NA ...
## $ thisRow.t : num 68.9 73.1 75.9 79.7 82.8 ...
## $ notes : logi NA NA NA NA NA NA ...
## $ slider_exp.response : int 5 7 5 4 3 4 3 2 7 5 ...
## $ slider_exp.rt : num 3.89 2.6 3.45 2.87 4.87 ...
## $ participant : int 17202 17202 17202 17202 17202 17202 17202 17202 17202 17202 ...
## $ session : Factor w/ 2 levels "Pre","Post": 2 2 2 2 2 2 2 2 2 2 ...
## $ group : Factor w/ 4 levels "1100","1600",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ date : Factor w/ 122 levels "2025-02-25_13h08.44.108",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ expName : Factor w/ 1 level "SimilarityRatings": 1 1 1 1 1 1 1 1 1 1 ...
## $ psychopyVersion : Factor w/ 1 level "2023.2.3": 1 1 1 1 1 1 1 1 1 1 ...
## $ frameRate : num 60 60 60 60 60 ...
## $ expStart : Factor w/ 122 levels "2025-02-25 13h09.40.421447 +0200",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ sbj : Factor w/ 61 levels "341","1003","1206",..: 10 10 10 10 10 10 10 10 10 10 ...
## $ resp : int 5 7 5 4 3 4 3 2 7 5 ...
## $ block : num 1 1 1 1 1 1 1 1 1 1 ...
## $ excl_instr : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
#test if there's a difference in response distribution among participants (excl_instr= "yes" vs. "no")
table_sim<-table(data$resp, data$excl_instr)
chi_sim<-chisq.test(table_sim)
## Warning in chisq.test(table_sim): Chi-squared approximation may be incorrect
print(chi_sim)
##
## Pearson's Chi-squared test
##
## data: table_sim
## X-squared = 18.417, df = 8, p-value = 0.01831
#function to calculate the percentage of expected counts >=5)
check_expected_frequencies <- function(chisq_test_result) {
expected <- chisq_test_result$expected
total_cells <- length(expected)
cells_over_5 <- sum(expected >= 5)
percent_over_5 <- (cells_over_5 / total_cells) * 100
return(percent_over_5)
}
check_expected_frequencies(chi_sim)
## [1] 94.44444
# it should be >80%
#plot distribution
plot_sim <- data %>%
group_by(resp, excl_instr) %>%
summarise(n = n()) %>%
group_by(excl_instr) %>%
mutate(freq = n / sum(n))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by resp and excl_instr.
## ℹ Output is grouped by resp.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(resp, excl_instr))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
ggplot(plot_sim, aes(x = as.factor(resp), y = freq, fill = excl_instr)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "Response Distribution", x = "Response (1...9)", y = "Relative Frequency", fill = "Exclude\nInstructions") +
theme_minimal()
# exclude five early "exclusion instruction" participants
data <- droplevels( data[!(data$sbj %in% c("17202","57393", "94372", "40242", "89443")), ])
str(data)
## 'data.frame': 8064 obs. of 39 variables:
## $ trl : int 7 5 9 11 19 3 2 13 12 8 ...
## $ cnd : Factor w/ 2 levels "between","within": 2 2 2 2 2 2 2 2 2 2 ...
## $ dim : Factor w/ 2 levels "irrel","rel": 1 1 1 1 2 1 1 2 1 1 ...
## $ size1 : num 3 3 2.6 2.6 3 2.2 2.2 2.2 2.6 3 ...
## $ hue1 : num 0.1 -0.3 -0.3 0.1 0.1 0.1 -0.1 -0.3 0.3 0.3 ...
## $ size2 : num 3.4 3.4 3 3 3 2.6 2.6 2.2 3 3.4 ...
## $ hue2 : num 0.1 -0.3 -0.3 0.1 0.3 0.1 -0.1 -0.1 0.3 0.3 ...
## $ practice_trials.thisRepN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisTrialN: logi NA NA NA NA NA NA ...
## $ practice_trials.thisN : logi NA NA NA NA NA NA ...
## $ practice_trials.thisIndex : logi NA NA NA NA NA NA ...
## $ blocks.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisTrialN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ blocks.thisIndex : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisRepN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ trials.thisTrialN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisN : int 0 1 2 3 4 5 6 7 8 9 ...
## $ trials.thisIndex : int 6 4 8 10 18 2 1 12 11 7 ...
## $ break_loop.thisRepN : logi NA NA NA NA NA NA ...
## $ break_loop.thisTrialN : logi NA NA NA NA NA NA ...
## $ break_loop.thisN : logi NA NA NA NA NA NA ...
## $ break_loop.thisIndex : logi NA NA NA NA NA NA ...
## $ thisRow.t : num 52.6 56.5 61 66.1 67.6 ...
## $ notes : logi NA NA NA NA NA NA ...
## $ slider_exp.response : int 9 8 6 7 8 4 3 4 2 3 ...
## $ slider_exp.rt : num 3.6 4.28 4.81 1.19 1.56 ...
## $ participant : int 81139 81139 81139 81139 81139 81139 81139 81139 81139 81139 ...
## $ session : Factor w/ 2 levels "Pre","Post": 2 2 2 2 2 2 2 2 2 2 ...
## $ group : Factor w/ 4 levels "1100","1600",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ date : Factor w/ 112 levels "2025-03-18_11h13.00.672",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ expName : Factor w/ 1 level "SimilarityRatings": 1 1 1 1 1 1 1 1 1 1 ...
## $ psychopyVersion : Factor w/ 1 level "2023.2.3": 1 1 1 1 1 1 1 1 1 1 ...
## $ frameRate : num 60 60 60 60 60 ...
## $ expStart : Factor w/ 112 levels "2025-03-18 11h13.54.820981 +0200",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ sbj : Factor w/ 56 levels "341","1003","1206",..: 42 42 42 42 42 42 42 42 42 42 ...
## $ resp : int 9 8 6 7 8 4 3 4 2 3 ...
## $ block : num 1 1 1 1 1 1 1 1 1 1 ...
## $ excl_instr : Factor w/ 1 level "no": 1 1 1 1 1 1 1 1 1 1 ...
#Check that there are 14 sbjs per group
data %>%
distinct(sbj, group) %>%
count(group)
## group n
## 1 1100 14
## 2 1600 14
## 3 600 14
## 4 RD 14
#Calculate standardized score for the resp variable.
data <- data %>%
group_by(sbj) %>%
mutate(
z_resp = (resp - mean(resp, na.rm = TRUE)) /
sd(resp, na.rm = TRUE)
) %>%
ungroup()
Similarity Ratings Plot
##########################################################
# Similarity Ratings, by Group, Session, and Pair Type
# Violin plots + within-subject trajectories (ggplot2)
##########################################################
#create new factor, condition (pair type) with three levels
data$condition<-as.factor(ifelse(data$cnd=="between", "between", ifelse(data$dim=="rel","within_rel","within_irrel" )))
ylb <- "Standardized Similarity Ratings"
xlb <- "Session"
group_levels <- c("600", "1100", "1600", "RD")
group_labels <- c("Group: 600 ms", "Group: 1100 ms", "Group: 1600 ms", "Group: Response Defined")
cond_levels <- c("between", "within_rel", "within_irrel")
cond_labels <- c("Between Category Pairs\nRelevant Dimension",
"Within Category Pairs\nRelevant Dimension",
"Within Category Pairs\nIrrelevant Dimension")
# Participant means per cell (sbj × group × condition × session)
plot_df <- data %>%
mutate(
group = factor(group, levels = group_levels, labels = group_labels),
condition = factor(condition, levels = cond_levels, labels = cond_labels)
) %>%
group_by(sbj, group, condition, session) %>%
summarise(z_resp = mean(z_resp, na.rm = TRUE), .groups = "drop") %>%
filter(!is.na(z_resp))
p_sim <- ggplot(plot_df, aes(x = session, y = z_resp)) +
geom_violin(aes(fill = session),
trim = FALSE, width = 0.8, color = NA, alpha = 0.55) +
geom_line(aes(group = sbj), color = "grey40", linewidth = 0.30, alpha = 0.70) +
geom_point(aes(group = sbj), color = "grey40", size = 0.80, alpha = 0.70) +
coord_cartesian(ylim = c(-1.5, 1.5)) +
facet_grid(group ~ condition) +
scale_fill_manual(values = c("grey25", "grey75")) +
labs(
title = "Similarity Ratings by Group, Session, and Pair Type",
x = xlb, y = ylb
) +
theme_classic(base_size = 11) +
theme(
legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(face = "bold"), # facet titles bold
axis.title.x = element_text(face = "bold"), # axis titles bold
axis.title.y = element_text(face = "bold"),
axis.text.x = element_text(face = "bold"), # Pre/Post bold
plot.title = element_text(face = "bold", hjust = 0.5)
)
print(p_sim)
#graph for paper 600 x 750
LMMs - Inferential Analyses
data_lmm <- data %>%
ungroup() %>%
mutate(
# Fixed order of factor levels
group = factor(
group,
levels = c("600", "1100", "1600", "RD")
),
session = factor(
session,
levels = c("Pre", "Post")
),
condition = factor(
condition,
levels = c("between", "within_rel", "within_irrel")
),
# Identifier of each individual stimulus
stim1_id = paste0(
sprintf("%.1f", size1), "_",
sprintf("%.1f", hue1)
),
stim2_id = paste0(
sprintf("%.1f", size2), "_",
sprintf("%.1f", hue2)
),
# Order-independent stimulus-pair identifier
pair_id = factor(
ifelse(
stim1_id <= stim2_id,
paste(stim1_id, stim2_id, sep = "__"),
paste(stim2_id, stim1_id, sep = "__")
)
)
) %>%
select(
sbj, group, session, condition,
pair_id, block, resp, z_resp
) %>%
droplevels()
###################################
# Code within-participant effects
###################################
data_lmm <- data_lmm %>%
mutate(
# Post - Pre
session_c = ifelse(session == "Post", 0.5, -0.5),
# Two contrasts spanning the main effect of Pair Type
condition_c1 = case_when(
condition == "between" ~ 2/3,
condition == "within_rel" ~ -1/3,
condition == "within_irrel" ~ -1/3
),
condition_c2 = case_when(
condition == "between" ~ 0,
condition == "within_rel" ~ 0.5,
condition == "within_irrel" ~ -0.5
),
# Post - Pre change separately for each Pair Type
change_between = ifelse(
condition == "between",
session_c,
0
),
change_within_rel = ifelse(
condition == "within_rel",
session_c,
0
),
change_within_irrel = ifelse(
condition == "within_irrel",
session_c,
0
)
)
contrasts(data_lmm$group) <- contr.sum(4)
contrasts(data_lmm$session) <- contr.sum(2)
contrasts(data_lmm$condition) <- contr.sum(3)
##########################################################
# Model 1 (m_lmm_full): Full linear mixed-effects model
##########################################################
m_lmm_full <- lmer(
z_resp ~ group * session * condition +
(
0 +
condition_c1 +
condition_c2 +
change_between +
change_within_rel +
change_within_irrel
| sbj
) +
# General differences among stimulus pairs
(1 | pair_id)+
# Pair-specific Post-Pre change
(0 + session_c | pair_id) +
# Repeated evaluations of the same pair
# by the same participant
(1 | sbj:pair_id),
data = data_lmm,
REML = TRUE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 200000)
)
)
## boundary (singular) fit: see help('isSingular')
#Check convergence and singularity
m_lmm_full@optinfo$conv$lme4$messages
## [1] "boundary (singular) fit: see help('isSingular')"
isSingular(m_lmm_full, tol = 1e-4)
## [1] TRUE
VarCorr(m_lmm_full)
## Groups Name Std.Dev. Corr
## sbj.pair_id (Intercept) 0.080370
## sbj condition_c1 0.264272
## condition_c2 0.805398 0.995
## change_between 0.635871 0.234 0.196
## change_within_rel 0.622590 0.341 0.289 0.958
## change_within_irrel 0.591068 -0.386 -0.329 0.352 0.127
## pair_id session_c 0.020164
## pair_id.1 (Intercept) 0.112224
## Residual 0.747113
#######################################################################
#Model 2 (m_lmm_ind): Model with independent participant random slopes
#######################################################################
m_lmm_ind <- lmer(
z_resp ~ group * session * condition +
# General differences among Pair Types
(0 + condition_c1 | sbj) +
(0 + condition_c2 | sbj) +
# Post-Pre change separately for each Pair Type
(0 + change_between | sbj) +
(0 + change_within_rel | sbj) +
(0 + change_within_irrel | sbj) +
# Differences among the 24 stimulus pairs
(1 | pair_id) +
# Pair-specific Post-Pre change
(0 + session_c | pair_id) +
# Repeated evaluations of the same pair
# by the same participant
(1 | sbj:pair_id),
data = data_lmm,
REML = TRUE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 200000)
)
)
# Check convergence and singularity
m_lmm_ind@optinfo$conv$lme4$messages
## NULL
isSingular(m_lmm_ind, tol = 1e-4)
## [1] FALSE
VarCorr(m_lmm_ind)
## Groups Name Std.Dev.
## sbj.pair_id (Intercept) 0.083721
## sbj change_within_irrel 0.591058
## sbj.1 change_within_rel 0.622586
## sbj.2 change_between 0.635766
## sbj.3 condition_c2 0.814599
## sbj.4 condition_c1 0.261193
## pair_id session_c 0.020160
## pair_id.1 (Intercept) 0.112180
## Residual 0.747119
##########################################################
#Type III tests of fixed effects
##########################################################
anova_lmm <- anova( m_lmm_ind, type = 3, ddf = "Kenward-Roger")
anova_lmm
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## group 18.868 6.2892 3 1207.02 11.2672 2.716e-07 ***
## session 10.197 10.1967 1 137.37 18.2675 3.566e-05 ***
## condition 34.170 17.0849 2 52.12 30.3728 1.801e-09 ***
## group:session 1.894 0.6313 3 154.16 1.1310 0.3384119
## group:condition 16.072 2.6787 6 84.90 4.7345 0.0003324 ***
## session:condition 3.644 1.8220 2 95.65 3.2434 0.0433637 *
## group:session:condition 2.632 0.4386 6 122.03 0.7783 0.5884901
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###############################################################
#Post-Pre change within each Pair Type averaged across Groups
###############################################################
emm_session_condition <- emmeans( m_lmm_ind, ~ session | condition, lmer.df = "kenward-roger", pbkrtest.limit = 10000)
## NOTE: Results may be misleading due to involvement in interactions
change_by_condition <- contrast(emm_session_condition,method = list("Post - Pre" = c(-1, 1)), by = "condition")
summary(change_by_condition,infer = c(TRUE, FALSE))
## condition = between:
## contrast estimate SE df lower.CL upper.CL
## Post - Pre -0.3542 0.0948 48.5 -0.545 -0.1637
##
## condition = within_rel:
## contrast estimate SE df lower.CL upper.CL
## Post - Pre -0.2566 0.0883 51.0 -0.434 -0.0793
##
## condition = within_irrel:
## contrast estimate SE df lower.CL upper.CL
## Post - Pre -0.0459 0.0826 51.6 -0.212 0.1199
##
## Results are averaged over the levels of: group
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
test(change_by_condition, by = NULL, adjust = "holm")
## contrast condition estimate SE df t.ratio p.value
## Post - Pre between -0.3542 0.0948 48.5 -3.738 0.0015
## Post - Pre within_rel -0.2566 0.0883 51.0 -2.905 0.0108
## Post - Pre within_irrel -0.0459 0.0826 51.6 -0.556 0.5809
##
## Results are averaged over the levels of: group
## Degrees-of-freedom method: kenward-roger
## P value adjustment: holm method for 3 tests
#############################################################################
# Theory-Guided Follow-Up Analysis
# Group Moderation of Diagnosticity-Dependent Perceptual Change
############################################################################
###############################################################
# Post-Pre change within each Pair Type and Group
###############################################################
emm_session_condition_group <- emmeans( m_lmm_ind, ~ session | condition * group, lmer.df = "kenward-roger", pbkrtest.limit = 10000)
change_by_condition_group <- contrast(emm_session_condition_group, method = list("Post - Pre" = c(-1, 1)), by = c("condition", "group"))
summary(change_by_condition_group, infer = c(TRUE, TRUE), adjust = "none")
## condition = between, group = 600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.3594 0.189 52.0 -0.7381 0.0193 -1.904 0.0624
##
## condition = within_rel, group = 600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.2638 0.176 52.1 -0.6175 0.0898 -1.497 0.1404
##
## condition = within_irrel, group = 600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.3455 0.165 52.1 -0.6764 -0.0145 -2.095 0.0411
##
## condition = between, group = 1100:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.3334 0.189 52.0 -0.7121 0.0453 -1.767 0.0831
##
## condition = within_rel, group = 1100:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.1196 0.176 52.1 -0.4732 0.2341 -0.678 0.5005
##
## condition = within_irrel, group = 1100:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.0924 0.165 52.1 -0.4234 0.2385 -0.560 0.5777
##
## condition = between, group = 1600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.2548 0.189 52.0 -0.6335 0.1239 -1.350 0.1828
##
## condition = within_rel, group = 1600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.2332 0.176 52.1 -0.5869 0.1204 -1.323 0.1915
##
## condition = within_irrel, group = 1600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre 0.2383 0.165 52.1 -0.0927 0.5692 1.445 0.1545
##
## condition = between, group = RD:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.4692 0.189 52.0 -0.8479 -0.0905 -2.486 0.0162
##
## condition = within_rel, group = RD:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.4099 0.176 52.1 -0.7635 -0.0562 -2.326 0.0240
##
## condition = within_irrel, group = RD:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre 0.0160 0.165 52.1 -0.3150 0.3469 0.097 0.9233
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
###############################################################
# Theory-guided diagnosticity-dependent change within each group
###############################################################
warping_by_group <- contrast(change_by_condition_group, method = list("Diagnosticity-dependent change" = c(-0.5, -0.5, 1)),by = "group")
summary(warping_by_group, infer = c(TRUE, FALSE))
## group = 600:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change -0.0339 0.209 114 -0.4488 0.381
##
## group = 1100:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change 0.1341 0.209 114 -0.2809 0.549
##
## group = 1600:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change 0.4823 0.209 114 0.0674 0.897
##
## group = RD:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change 0.4555 0.209 114 0.0405 0.870
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
###############################################################
# Theory-guided group contrast
###############################################################
theory_group_contrast <- contrast(warping_by_group, method = list( "Mean of 1100, 1600, and RD - 600" = c(-1, 1/3, 1/3, 1/3)), by = NULL)
summary(theory_group_contrast, infer = c(TRUE, TRUE), adjust = "none")
## contrast estimate SE df lower.CL upper.CL t.ratio
## Mean of 1100, 1600, and RD - 600 0.391 0.242 114 -0.0876 0.87 1.619
## p.value
## 0.1083
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
Exploratory Analyses, Accuracy ~ Perceptual Change
# Mean rating per participant, session, and pair type
sim_type <- aggregate(z_resp ~ sbj + group + session + condition, data = data,FUN = mean)
# Put Pre and Post in separate columns
sim_change <- reshape(sim_type, idvar = c("sbj", "group", "condition"), timevar = "session",direction = "wide")
# Calculate Post minus Pre change
sim_change$change <-
sim_change$z_resp.Post -
sim_change$z_resp.Pre
# Put the three pair types in separate columns
warping <- reshape(
sim_change[, c("sbj", "group", "condition", "change")],
idvar = c("sbj", "group"),
timevar = "condition",
direction = "wide"
)
# Calculate diagnosticity-dependent warping
warping$warping_index <-
warping$change.within_irrel -
(
warping$change.between +
warping$change.within_rel
) / 2
# Inspect the warping index
head(warping)
## sbj group change.between change.within_irrel change.within_rel
## 1 1003 1100 -0.6602845 0.1834124 -0.60526078
## 2 3402 1100 0.3057330 0.1528665 0.21019141
## 3 8227 1100 -1.2734562 0.2214706 -0.49830893
## 4 19052 1100 0.5853529 0.0000000 0.08780294
## 5 21505 1100 -0.7074160 -0.8449692 -0.56003769
## 6 24220 1100 -1.0469296 -1.1134013 -0.39883033
## warping_index
## 1 0.8161850
## 2 -0.1050957
## 3 1.1073532
## 4 -0.3365779
## 5 -0.2112423
## 6 -0.3905214
summary(warping$warping_index)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.4210 -0.3221 0.1755 0.2595 0.6630 2.9623
###############################
# Load categorization accuracy
###############################
d_cat <- read.csv("d_cat.csv", header = TRUE)
d_cat$sbj <- as.factor(d_cat$sbj)
warping$sbj <- as.factor(warping$sbj)
# Merge similarity and accuracy data
d_cor <- merge( d_cat,warping,by = c("sbj", "group")
)
# Keep complete observations
d_cor <- d_cor[complete.cases(d_cor$acc, d_cor$warping_index),]
###############################
# Pearson correlation
###############################
cor_accuracy_warping <- cor.test(d_cor$acc, d_cor$warping_index,method = "pearson")
cor_accuracy_warping
##
## Pearson's product-moment correlation
##
## data: d_cor$acc and d_cor$warping_index
## t = 0.30918, df = 54, p-value = 0.7584
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2233323 0.3016038
## sample estimates:
## cor
## 0.04203636
###############################
# Simple linear regression
###############################
model_accuracy_warping <- lm(warping_index ~ acc, data = d_cor)
summary(model_accuracy_warping)
##
## Call:
## lm(formula = warping_index ~ acc, data = d_cor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.70712 -0.58434 -0.08248 0.36673 2.66998
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6079 2.8077 -0.217 0.829
## acc 1.0003 3.2353 0.309 0.758
##
## Residual standard error: 0.8054 on 54 degrees of freedom
## Multiple R-squared: 0.001767, Adjusted R-squared: -0.01672
## F-statistic: 0.09559 on 1 and 54 DF, p-value: 0.7584
# Confidence intervals for regression coefficients
confint(model_accuracy_warping)
## 2.5 % 97.5 %
## (Intercept) -6.237015 5.021156
## acc -5.486143 7.486711
Sensitivity Analyses - Category-Learning Sample
# Prepare full sample
data_full <- data_full %>% group_by(sbj) %>% mutate(z_resp = (resp - mean(resp, na.rm = TRUE)) / sd(resp, na.rm = TRUE)) %>% ungroup()
data_full$group <- factor(data_full$group, levels = c("600", "1100", "1600", "RD"))
data_full$session <- factor(data_full$session, levels = c("Pre", "Post"))
data_full$condition <- factor(ifelse(data_full$cnd == "between", "between", ifelse(data_full$dim == "rel", "within_rel", "within_irrel")), levels = c("between", "within_rel", "within_irrel"))
# Create stimulus-pair identifier
stim1_id <- paste0(sprintf("%.1f", data_full$size1), "_", sprintf("%.1f", data_full$hue1))
stim2_id <- paste0(sprintf("%.1f", data_full$size2), "_", sprintf("%.1f", data_full$hue2))
data_full$pair_id <- factor(ifelse(stim1_id <= stim2_id, paste(stim1_id, stim2_id, sep = "__"), paste(stim2_id, stim1_id, sep = "__")))
# Code within-participant effects
data_full$session_c <- ifelse(data_full$session == "Post", 0.5, -0.5)
data_full$condition_c1 <- ifelse(data_full$condition == "between", 2/3, -1/3)
data_full$condition_c2 <- ifelse(data_full$condition == "between", 0, ifelse(data_full$condition == "within_rel", 0.5, -0.5))
data_full$change_between <- ifelse(data_full$condition == "between", data_full$session_c, 0)
data_full$change_within_rel <- ifelse(data_full$condition == "within_rel", data_full$session_c, 0)
data_full$change_within_irrel <- ifelse(data_full$condition == "within_irrel", data_full$session_c, 0)
data_lmm_sens <- droplevels(data_full)
contrasts(data_lmm_sens$group) <- contr.sum(4)
contrasts(data_lmm_sens$session) <- contr.sum(2)
contrasts(data_lmm_sens$condition) <- contr.sum(3)
##########################################################
# Final LMM fitted to the full sample
##########################################################
m_lmm_sens <- lmer(
z_resp ~ group * session * condition +
(0 + condition_c1 | sbj) + (0 + condition_c2 | sbj) +
(0 + change_between | sbj) + (0 + change_within_rel | sbj) + (0 + change_within_irrel | sbj) +
(1 | pair_id) + (0 + session_c | pair_id) + (1 | sbj:pair_id),
data = data_lmm_sens, REML = TRUE,
control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 200000))
)
m_lmm_sens@optinfo$conv$lme4$messages
## NULL
isSingular(m_lmm_sens, tol = 1e-4)
## [1] FALSE
VarCorr(m_lmm_sens)
## Groups Name Std.Dev.
## sbj.pair_id (Intercept) 0.0890991
## sbj change_within_irrel 0.5849603
## sbj.1 change_within_rel 0.6406591
## sbj.2 change_between 0.6325847
## sbj.3 condition_c2 0.7912737
## sbj.4 condition_c1 0.2632722
## pair_id session_c 0.0095334
## pair_id.1 (Intercept) 0.1095738
## Residual 0.7483707
##########################################################
# Type III tests
##########################################################
anova_lmm_sens <- anova(m_lmm_sens, type = 3, ddf = "Kenward-Roger")
anova_lmm_sens
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## group 23.117 7.7056 3 1427.02 13.7585 7.556e-09 ***
## session 15.504 15.5035 1 160.35 27.6819 4.522e-07 ***
## condition 37.888 18.9440 2 53.13 33.5987 3.704e-10 ***
## group:session 2.297 0.7656 3 183.61 1.3670 0.254332
## group:condition 19.856 3.3093 6 101.68 5.8421 2.903e-05 ***
## session:condition 8.204 4.1021 2 112.04 7.2847 0.001061 **
## group:session:condition 2.154 0.3590 6 145.86 0.6359 0.701362
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##########################################################
# Post-Pre change within each Pair Type
##########################################################
emm_session_condition_sens <- emmeans(m_lmm_sens, ~ session | condition, lmer.df = "kenward-roger", pbkrtest.limit = 10000)
## NOTE: Results may be misleading due to involvement in interactions
change_by_condition_sens <- contrast(emm_session_condition_sens, method = list("Post - Pre" = c(-1, 1)), by = "condition")
summary(change_by_condition_sens, infer = c(TRUE, FALSE))
## condition = between:
## contrast estimate SE df lower.CL upper.CL
## Post - Pre -0.41165 0.0868 56.3 -0.586 -0.238
##
## condition = within_rel:
## contrast estimate SE df lower.CL upper.CL
## Post - Pre -0.33018 0.0835 60.3 -0.497 -0.163
##
## condition = within_irrel:
## contrast estimate SE df lower.CL upper.CL
## Post - Pre -0.00584 0.0754 60.9 -0.157 0.145
##
## Results are averaged over the levels of: group
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
test(change_by_condition_sens, by = NULL, adjust = "holm")
## contrast condition estimate SE df t.ratio p.value
## Post - Pre between -0.41165 0.0868 56.3 -4.742 <0.0001
## Post - Pre within_rel -0.33018 0.0835 60.3 -3.955 0.0004
## Post - Pre within_irrel -0.00584 0.0754 60.9 -0.077 0.9386
##
## Results are averaged over the levels of: group
## Degrees-of-freedom method: kenward-roger
## P value adjustment: holm method for 3 tests
##########################################################
# Theory-guided sensitivity contrast
##########################################################
emm_session_condition_group_sens <- emmeans(m_lmm_sens, ~ session | condition * group, lmer.df = "kenward-roger", pbkrtest.limit = 10000)
change_by_condition_group_sens <- contrast(emm_session_condition_group_sens, method = list("Post - Pre" = c(-1, 1)), by = c("condition", "group"))
summary(change_by_condition_group_sens, infer = c(TRUE, TRUE), adjust = "none")
## condition = between, group = 600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.4262 0.166 61.6 -0.7574 -0.0951 -2.573 0.0125
##
## condition = within_rel, group = 600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.3698 0.159 61.9 -0.6884 -0.0511 -2.320 0.0237
##
## condition = within_irrel, group = 600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.2660 0.144 62.0 -0.5540 0.0219 -1.847 0.0695
##
## condition = between, group = 1100:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.4262 0.181 61.7 -0.7889 -0.0635 -2.349 0.0220
##
## condition = within_rel, group = 1100:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.1993 0.175 61.9 -0.5484 0.1497 -1.142 0.2580
##
## condition = within_irrel, group = 1100:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.0378 0.158 62.0 -0.3532 0.2776 -0.240 0.8113
##
## condition = between, group = 1600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.2905 0.170 61.6 -0.6312 0.0503 -1.704 0.0934
##
## condition = within_rel, group = 1600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.2454 0.164 61.9 -0.5733 0.0825 -1.496 0.1397
##
## condition = within_irrel, group = 1600:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre 0.2067 0.148 62.0 -0.0896 0.5030 1.395 0.1681
##
## condition = between, group = RD:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.5037 0.176 61.7 -0.8549 -0.1524 -2.867 0.0057
##
## condition = within_rel, group = RD:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre -0.5062 0.169 61.9 -0.8442 -0.1683 -2.994 0.0039
##
## condition = within_irrel, group = RD:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Post - Pre 0.0738 0.153 62.0 -0.2316 0.3792 0.483 0.6307
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
warping_by_group_sens <- contrast(change_by_condition_group_sens, method = list("Diagnosticity-dependent change" = c(-0.5, -0.5, 1)), by = "group")
summary(warping_by_group_sens, infer = c(TRUE, FALSE))
## group = 600:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change 0.132 0.184 137 -0.2324 0.496
##
## group = 1100:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change 0.275 0.202 138 -0.1242 0.674
##
## group = 1600:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change 0.475 0.190 137 0.0997 0.850
##
## group = RD:
## contrast estimate SE df lower.CL upper.CL
## Diagnosticity-dependent change 0.579 0.195 138 0.1923 0.965
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
theory_group_contrast_sens <- contrast(warping_by_group_sens, method = list("Mean of 1100, 1600, and RD - 600" = c(-1, 1/3, 1/3, 1/3)), by = NULL)
summary(theory_group_contrast_sens, infer = c(TRUE, TRUE), adjust = "none")
## contrast estimate SE df lower.CL upper.CL t.ratio
## Mean of 1100, 1600, and RD - 600 0.311 0.216 138 -0.116 0.738 1.438
## p.value
## 0.1526
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
#######################################################
# Sensitivity Analysis
# Accuracy-Warping Relationship in the Full Sample
#######################################################
# Mean standardized rating per participant, session, and pair type
sim_type_full <- aggregate(z_resp ~ sbj + group + session + condition, data = data_full, FUN = mean)
# Calculate Post-Pre change
sim_change_full <- reshape(sim_type_full, idvar = c("sbj", "group", "condition"), timevar = "session", direction = "wide")
sim_change_full$change <- sim_change_full$z_resp.Post - sim_change_full$z_resp.Pre
# Put pair types in separate columns
warping_full <- reshape(sim_change_full[, c("sbj", "group", "condition", "change")], idvar = c("sbj", "group"), timevar = "condition", direction = "wide")
# Calculate diagnosticity-dependent warping
warping_full$warping_index <- warping_full$change.within_irrel - (warping_full$change.between + warping_full$change.within_rel) / 2
# Prepare one accuracy score per participant
accuracy_full <- aggregate(acc ~ sbj + group, data = d_cat, FUN = mean, na.rm = TRUE)
# Use common identifier formats
accuracy_full$sbj <- as.character(accuracy_full$sbj)
accuracy_full$group <- as.character(accuracy_full$group)
warping_full$sbj <- as.character(warping_full$sbj)
warping_full$group <- as.character(warping_full$group)
# Merge similarity and accuracy data
d_cor_full <- merge(accuracy_full, warping_full, by = c("sbj", "group"))
d_cor_full <- d_cor_full[complete.cases(d_cor_full$acc, d_cor_full$warping_index), ]
# Check the analysis sample
cat("Number of participants:", nrow(d_cor_full), "\n")
## Number of participants: 66
print(table(d_cor_full$group))
##
## 1100 1600 600 RD
## 15 17 18 16
# Pearson correlation
cor_accuracy_warping_full <- cor.test(d_cor_full$acc, d_cor_full$warping_index, method = "pearson")
print(cor_accuracy_warping_full)
##
## Pearson's product-moment correlation
##
## data: d_cor_full$acc and d_cor_full$warping_index
## t = 0.57373, df = 64, p-value = 0.5682
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1735038 0.3082292
## sample estimates:
## cor
## 0.07153294
# Simple linear regression
model_accuracy_warping_full <- lm(warping_index ~ acc, data = d_cor_full)
summary(model_accuracy_warping_full)
##
## Call:
## lm(formula = warping_index ~ acc, data = d_cor_full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8197 -0.5734 -0.1241 0.4539 2.5570
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5573 1.6041 -0.347 0.729
## acc 1.0696 1.8643 0.574 0.568
##
## Residual standard error: 0.8531 on 64 degrees of freedom
## Multiple R-squared: 0.005117, Adjusted R-squared: -0.01043
## F-statistic: 0.3292 on 1 and 64 DF, p-value: 0.5682
confint(model_accuracy_warping_full)
## 2.5 % 97.5 %
## (Intercept) -3.761984 2.647307
## acc -2.654728 4.793927
Session Information
sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Europe/Athens
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lmerTest_3.2-1 lme4_2.0-6 Matrix_1.7-5 emmeans_2.0.4 ggplot2_4.0.3
## [6] dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyr_1.3.2 sass_0.4.10 generics_0.1.4
## [4] lattice_0.22-9 digest_0.6.39 magrittr_2.0.5
## [7] evaluate_1.0.5 grid_4.6.1 estimability_2.0.0
## [10] RColorBrewer_1.1-3 mvtnorm_1.4-2 fastmap_1.2.0
## [13] jsonlite_2.0.0 backports_1.5.1 purrr_1.2.2
## [16] scales_1.4.0 numDeriv_2016.8-1.1 jquerylib_0.1.4
## [19] reformulas_0.4.4 Rdpack_2.6.6 cli_3.6.6
## [22] rlang_1.3.0 rbibutils_2.4.1 splines_4.6.1
## [25] withr_3.0.3 cachem_1.1.0 yaml_2.3.12
## [28] parallel_4.6.1 pbkrtest_0.5.5 tools_4.6.1
## [31] nloptr_2.2.1 minqa_1.2.8 boot_1.3-32
## [34] broom_1.0.13 vctrs_0.7.3 R6_2.6.1
## [37] lifecycle_1.0.5 MASS_7.3-65 pkgconfig_2.0.3
## [40] pillar_1.11.1 bslib_0.11.0 gtable_0.3.6
## [43] glue_1.8.1 Rcpp_1.1.2 xfun_0.60
## [46] tibble_3.3.1 tidyselect_1.2.1 rstudioapi_0.19.0
## [49] knitr_1.51 farver_2.1.2 htmltools_0.5.9
## [52] nlme_3.1-169 rmarkdown_2.31 labeling_0.4.3
## [55] compiler_4.6.1 S7_0.2.2