Predictors of social outcomes
Lena Matyjek
2026-10-06
load("./data/data_all.RData")
DFmain <- data_all
rm(data_all)
# For the social interaction outcomes:
DFoutcomes <- DFmain[,c(1:54)]
# Take one task only (they are the same for T1 and T2)
DFoutcomes = DFoutcomes[DFoutcomes$task=="T1",]
DFoutcomes$task <- NULLData preparation
Two tasks: T1 and T2
Nonverbal data are available for both tasks T1 and T2, but transcripts are only available for T1. To analyse predictors of rapportand social motivation, we first check if it makes sense to take means of the non-verbal values across tasks.
Non verbal behaviours: correlations between the tasks
choose_cols <- c("INTER_smile_av_sync","INTER_laughter_av_sync","INTER_expressivity_av_sync", "INTER_headmov_av_sync",
"S_INTRA_per_smile","S_INTRA_per_laughter","S_INTRA_m_expressivity", "S_INTRA_headmov")
pretty_names <- c(
INTER_smile_av_sync = "Smile synch. (inter)",
INTER_laughter_av_sync = "Laughter synch. (inter)",
INTER_expressivity_av_sync = "Expressivity synch. (inter)",
INTER_headmov_av_sync = "Head movement synch. (inter)",
S_INTRA_per_smile = "Smile (%) (intra)",
S_INTRA_per_laughter = "Laughter (%) (intra)",
S_INTRA_m_expressivity = "Mean expressivity (intra)",
S_INTRA_headmov = "Head movement (intra)"
)
facet_order <- c(
"S_INTRA_per_smile", "INTER_smile_av_sync",
"S_INTRA_per_laughter", "INTER_laughter_av_sync",
"S_INTRA_m_expressivity", "INTER_expressivity_av_sync",
"S_INTRA_headmov", "INTER_headmov_av_sync",
"S_INTRA_distance_to_camera", "INTER_distance_av_sync"
)
DF_wide <- DFmain %>%
dplyr::select(dyad_indx, task, all_of(choose_cols)) %>%
dplyr::distinct(dyad_indx, task, .keep_all = TRUE) %>% # <-- keep one row per dyad/task
pivot_wider(names_from = task,
values_from = all_of(choose_cols),
names_glue = "{.value}_{task}")
# Now build plot_df
plot_df <- bind_rows(lapply(choose_cols, function(col) {
data.frame(
variable = col,
T1 = DF_wide[[paste0(col, "_T1")]],
T2 = DF_wide[[paste0(col, "_T2")]]
)
})) %>%
dplyr::filter(!is.na(T1), !is.na(T2))
# make sure variable is character/factor with these levels
plot_df <- plot_df %>%
dplyr::mutate(variable = factor(variable, levels = names(pretty_names)))
# Build long dataframe for faceted plotting
plot_df <- bind_rows(lapply(choose_cols, function(col) {
data.frame(
variable = col,
T1 = DF_wide[[paste0(col, "_T1")]],
T2 = DF_wide[[paste0(col, "_T2")]]
)
})) %>%
dplyr::mutate(variable = factor(variable, levels = facet_order)) %>%
filter(!is.na(T1), !is.na(T2))
# Plot
ggplot(plot_df, aes(x = T1, y = T2)) +
geom_point(alpha = 0.4, size = 1.2) +
geom_smooth(method = "lm", se = TRUE, colour = "steelblue") +
ggpubr::stat_cor(method = "pearson", size = 4, label.x.npc = "left") +
facet_wrap(~ variable, scales = "free") +
labs(
title = "T1 vs T2 correlations for nonverbal measures",
x = "Task 1", y = "Task 2"
) +
facet_wrap(
~ variable,
ncol=2,
scales = "free",
labeller = as_labeller(pretty_names)
) +
theme_minimal(base_size = 15) +
theme(
panel.border = element_rect(color = "grey35", linewidth = 0.6, fill = NA),
panel.spacing = grid::unit(0.8, "lines"),
strip.background = element_rect(fill = "grey95", color = "grey35", linewidth = 0.6)
)
All of the non-verbal behaviours are correlated between T1 and T2.
For the rapport predictors, however, we will only use T1 in which all measurements are available.
collapse_across_tasks = 0 # 1 = take means of the two tasks, ignore NAs. 0 = take T1 only# Means for the two tasks?
if (collapse_across_tasks == 1) {
DF_taskcollapse <- DFmain %>%
dplyr::group_by(dyad_indx, ID, partner_ID, dyad, Group) %>%
dplyr::summarise(
across(matches("^S_INTRA_|^INTER_|^P_INTRA_|^R_|^Q_"),
~ mean(., na.rm = TRUE)),
.groups = "drop"
) %>% as.data.frame()
DF_pca_analysis <- DF_taskcollapse
} else {
DF_taskcollapse <- DFmain %>%
dplyr::filter(task == "T1") %>%
dplyr::group_by(dyad_indx, ID, partner_ID, dyad, Group) %>%
dplyr::summarise(
across(matches("^S_INTRA_|^INTER_|^P_INTRA_|^R_|^Q_"),
~ mean(., na.rm = TRUE)),
.groups = "drop"
) %>% as.data.frame()
DF_pca_analysis <- DF_taskcollapse
}# Add partner neurotype information
DF_pca_analysis <- DF_pca_analysis %>%
mutate(partner_neurotype = case_when(
dyad == "NAUT" ~ "NAUT",
dyad == "AUT" ~ "AUT",
dyad == "MIX" & Group == "NAUT" ~ "AUT",
dyad == "MIX" & Group == "AUT" ~ "NAUT"
))
DF_pca_analysis$partner_neurotype <- as.factor(DF_pca_analysis$partner_neurotype)Dimension reduction
Disposition: PCA
Step 1: A priori decision to take the first 3 PCA components for analysis.
n_pc_manual = 3### Take data
df = DF_pca_analysis[,grepl("Q_",colnames(DF_pca_analysis))]
#cat(paste("Number of variables before NA rejections: ",ncol(df),".",sep=""))
names1 = colnames(df)
df_complete <- df[, colSums(is.na(df)) == 0]
names2 = colnames(df_complete)
#cat(paste("Number of variables entering PCA: ",ncol(df_complete),".",sep=""))
#cat("Removed columns:", paste(setdiff(names1, names2), collapse = ", "))Step 2: Z-score all the variables.
df_z <- as.data.frame(scale(df_complete))Step 3: Run PCA and see the summary with the proportion of variance explained by each principal component.
pca_result <- prcomp(df_z)
summary(pca_result)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 2.1747 1.4349 1.19230 1.10309 1.03866 1.00009 0.95097 0.87785 0.75060 0.73723 0.66539 0.6132 0.58461
## Proportion of Variance 0.2956 0.1287 0.08885 0.07605 0.06743 0.06251 0.05652 0.04816 0.03521 0.03397 0.02767 0.0235 0.02136
## Cumulative Proportion 0.2956 0.4243 0.51313 0.58918 0.65661 0.71912 0.77564 0.82380 0.85901 0.89298 0.92066 0.9442 0.96551
## PC14 PC15 PC16
## Standard deviation 0.52212 0.39412 0.35190
## Proportion of Variance 0.01704 0.00971 0.00774
## Cumulative Proportion 0.98255 0.99226 1.00000View the loadings (rotation matrix, i.e., the contribution of each variable to the principal components).
pca_result$rotation## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## Q_state_stress 0.24361071 -0.25084548 -0.41854944 -0.07992763 0.194472545 0.12855693 0.227270238 -0.04013858
## Q_state_tiredness 0.21699779 -0.23564112 -0.32869110 -0.17756303 0.387963857 0.28812346 -0.250782356 -0.19247068
## Q_state_anxiety -0.16836126 0.17073749 0.29149222 0.32949701 -0.076727761 0.61832265 -0.167338008 0.01087989
## Q_AUT_familiarity 0.25576581 0.39707299 0.01687436 -0.21441941 -0.232335865 -0.05351296 0.172124426 -0.02954445
## Q_AQ 0.40887616 0.10034261 -0.04357789 0.00562059 -0.159523505 0.10126235 -0.081934786 0.02892629
## Q_SIAS 0.40023952 -0.07068449 0.03785748 0.12620429 -0.144097834 0.15527448 -0.066282306 -0.06980925
## Q_CATQ 0.34455075 -0.16958716 0.03292930 0.02224279 -0.121174016 -0.20717353 0.222197545 0.14947609
## Q_PHQ 0.31505594 -0.06198786 0.08186266 0.27094944 0.008236217 0.26114886 0.356479402 0.09672154
## Q_BFI_extraversion -0.28749362 0.02606598 -0.08639583 -0.14211122 -0.018999202 0.10804708 0.692391229 0.23064202
## Q_BFI_agreeableness -0.17973247 -0.02184557 -0.14414571 -0.43225402 -0.382109400 0.45505842 0.128773001 -0.42849327
## Q_BFI_Conscienctiousness 0.07940026 -0.02330951 0.48262437 -0.04098208 0.586201298 0.20298581 0.265120356 -0.07455319
## Q_BFI_Neuroticism 0.29098220 -0.11149700 0.28293710 -0.06921485 -0.351879382 0.08381212 -0.002409232 0.07083401
## Q_BFI_Openess -0.05551055 -0.33651483 0.10776403 -0.43329943 -0.056279733 0.20124032 -0.224916312 0.70202817
## Q_MARS 0.07795623 -0.13393274 0.50587302 -0.46296389 0.061912587 -0.21443824 -0.022451626 -0.33955395
## Q_ASKQ 0.17310382 0.48821268 -0.08715945 -0.22106218 0.190089822 -0.01483809 0.021168546 0.06153138
## Q_SATA 0.09696574 0.51322598 -0.06378656 -0.23067234 0.184065524 0.12690931 -0.149831296 0.25226981
## PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16
## Q_state_stress -0.238883411 -0.20417601 -0.40737778 -0.278360136 -0.336281275 0.310523969 -0.15470403 0.113627745
## Q_state_tiredness 0.118417246 0.23835769 0.11214826 0.164323642 -0.060074322 -0.555990662 0.03368785 -0.008061273
## Q_state_anxiety -0.384548669 -0.08489133 -0.11329307 -0.200381630 -0.136863414 -0.271399477 -0.11138147 0.140603458
## Q_AUT_familiarity -0.043128729 -0.25790716 0.18752588 0.162774225 -0.631404709 -0.221325617 0.13811731 -0.215736551
## Q_AQ 0.005249095 -0.10277508 0.24631199 0.237607661 0.116127403 0.165589011 0.05496845 0.778411896
## Q_SIAS 0.021529067 0.01704907 -0.03355617 -0.375540297 0.193484420 0.077913762 0.72594583 -0.223372933
## Q_CATQ -0.234464829 -0.39008654 0.01679136 -0.093808774 0.436616417 -0.468682434 -0.28842981 -0.109155204
## Q_PHQ -0.218031387 0.40890910 0.06720590 0.470561945 0.050300850 0.233551051 -0.13406767 -0.310391775
## Q_BFI_extraversion 0.054276604 0.22884288 -0.08072558 -0.072575420 0.050553161 -0.295235437 0.31134876 0.301023568
## Q_BFI_agreeableness 0.070963656 -0.17142813 0.11764712 0.045427020 0.292635932 0.163002476 -0.13644779 -0.158654528
## Q_BFI_Conscienctiousness 0.338657248 -0.36351190 0.18048119 -0.048974250 0.009410675 0.125749264 -0.01071200 -0.006820292
## Q_BFI_Neuroticism 0.563130414 0.27781129 -0.32459741 -0.197063152 -0.157659387 -0.098065450 -0.32585950 0.055106898
## Q_BFI_Openess -0.119315507 -0.04745012 0.21251077 -0.007064598 -0.087157973 0.128002461 0.05163414 -0.094197768
## Q_MARS -0.446057985 0.23328474 -0.22577221 0.069255564 -0.009837144 -0.019090620 0.10279507 0.141198138
## Q_ASKQ -0.135552869 0.36919601 0.33348368 -0.520578884 0.128208299 0.090038469 -0.26856890 -0.032589768
## Q_SATA 0.069332020 -0.11615468 -0.58414977 0.277014076 0.290773984 0.009352228 0.07104518 -0.088835452#round(pca_result$rotation*100,0)A plot of the eigenvalues (variance explained) by each component:
screeplot(pca_result, type = "lines")
A plot of the cumulative proportion of variance explained:
explained_variance <- pca_result$sdev^2 / sum(pca_result$sdev^2) # Variance explained
cumulative_variance <- cumsum(explained_variance) # Cumulative variance
plot(cumulative_variance, type = "b", xlab = "Number of Components", ylab = "Cumulative Proportion of Variance Explained", main = "Cumulative Variance Explained")
How much variance is explained with the first 3 components we decided on a priori?
### Extract the first N principal components:
#n_sel_comp = num_kaiser_components
n_sel_comp = n_pc_manual
pca_scores <- pca_result$x[, 1:n_sel_comp]
#How much variance do these components explain?
explained_variance <- pca_result$sdev^2 / sum(pca_result$sdev^2)
cat(paste("Variance explained with ",n_sel_comp," components: ",round(sum(explained_variance[1:n_sel_comp]),2)*100,"%",sep=""))## Variance explained with 3 components: 51%See the loadings (weights) of the original variables on each selected PCA component.
loadings <- pca_result$rotation
head(loadings[, 1:n_sel_comp],10)## PC1 PC2 PC3
## Q_state_stress 0.2436107 -0.25084548 -0.41854944
## Q_state_tiredness 0.2169978 -0.23564112 -0.32869110
## Q_state_anxiety -0.1683613 0.17073749 0.29149222
## Q_AUT_familiarity 0.2557658 0.39707299 0.01687436
## Q_AQ 0.4088762 0.10034261 -0.04357789
## Q_SIAS 0.4002395 -0.07068449 0.03785748
## Q_CATQ 0.3445507 -0.16958716 0.03292930
## Q_PHQ 0.3150559 -0.06198786 0.08186266
## Q_BFI_extraversion -0.2874936 0.02606598 -0.08639583
## Q_BFI_agreeableness -0.1797325 -0.02184557 -0.14414571#loadings <- round((loadings[, 1:n_sel_comp]*100),0)
#loadingsBi-plots:
#biplot(pca_result, choices = 1:2, main = "Biplot of PC1 and PC2")
ggbiplot(pca_result, choices = 1:2) +
ggtitle("Biplot of PC1 and PC2") +
theme_minimal()
#biplot(pca_result, choices = 1:2, main = "Biplot of PC1 and PC2")
ggbiplot(pca_result, choices = c(1,3)) +
ggtitle("Biplot of PC1 and PC3") +
theme_minimal()
#biplot(pca_result, choices = 1:2, main = "Biplot of PC1 and PC2")
ggbiplot(pca_result, choices = 2:3) +
ggtitle("Biplot of PC2 and PC3") +
theme_minimal()
### Save for the models:
pca_scores <- as.data.frame(pca_scores)
names(pca_scores) <- paste("Q_",names(pca_scores),sep="")
Q_PCAs = pca_scores
Q_PCAs_loadings = loadingsLoadings plot
# prepare and clean data
df_loadings <- as.data.frame(loadings[, 1:3])
df_loadings$Variable <- rownames(df_loadings)
comp_names <- c(
"PC1" = "Autism-like\ntraits",
"PC2" = "Attitudes\ntowards autism",
"PC3" = "Cognitive\nengagement"
)
df_long <- df_loadings %>%
pivot_longer(cols = c("PC1", "PC2", "PC3"), names_to = "PC", values_to = "Loading") %>%
mutate(
# Clean variable label
Var_Clean = gsub("^Q_", "", Variable),
Var_Clean = gsub("BFI_", "BFI: ", Var_Clean),
Var_Clean = gsub("_", " ", Var_Clean),
# Combine questionnaire name with its loading value
#Cell_Label = sprintf("%s (%.2f)", Var_Clean, Loading),
Cell_Label = sprintf("%s", Var_Clean),
Component = factor(comp_names[PC], levels = comp_names)
)
# rank independently within each component
df_ranked <- df_long %>%
group_by(Component) %>%
arrange(Loading) %>%
mutate(Rank = row_number()) %>%
ungroup()
# create the 3-panel ranked heatmap plot
p_ranked_loadings <- ggplot(df_ranked, aes(x = 1, y = Rank, fill = Loading)) +
geom_tile(color = "white", linewidth = 0.8, width = 0.95, height = 0.95) +
geom_text(aes(label = Cell_Label), size = 3.2, color = "black") +
scale_fill_gradient2(
low = "black", mid = "#F7F7F7", high = "orange",
midpoint = 0,
limits = c(-0.6, 0.6),
breaks = seq(-0.6, 0.6, by = 0.3), # Adds more tick marks
name = "Loading",
guide = guide_colorbar(
reverse = TRUE,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(14, "cm"), # Adjust width in cm (or use unit(0.85, "npc") to stretch relative to plot area)
barheight = unit(0.5, "cm"), # Controls thickness of the horizontal bar
ticks = TRUE,
ticks.colour = "grey20",
frame.colour = "grey60", # Subtle border around the bar
frame.linewidth = 0.5
)
) +
facet_wrap(~ Component, nrow = 1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = NULL, y = NULL) +
theme_minimal(base_size = 12) +
theme(
strip.text = element_text(face = "bold", size = 11, lineheight = 1.1),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
panel.spacing = unit(0.8, "lines"),
legend.position = "bottom",
legend.box = "horizontal",
legend.title = element_text(size = 10, face = "bold", margin = margin(b = 5)),
legend.text = element_text(size = 9)
)
print(p_ranked_loadings)
### Save
# Open a high-resolution graphics device (12x12 inches)
png("./figures/disposition_loadings.png", width = 12, height = 12, units = "in", res = 300)
df_ranked <- df_ranked %>%
mutate(
# Wrap long questionnaire labels so they fit neatly inside the tiles
#Var_Clean_Wrapped = str_wrap(Var_Clean, width = 18),
Var_Clean_Wrapped = Var_Clean, # or not
# Use exact color strings: white on dark tiles (Loading < -0.25), black elsewhere
Text_Color = ifelse(Loading < -0.25, "white", "black")
)
ggplot(df_ranked, aes(x = 1, y = Rank, fill = Loading)) +
geom_tile(color = "white", linewidth = 0.8, width = 0.95, height = 0.95) +
geom_text(
aes(label = Var_Clean_Wrapped, color = Text_Color),
size = 7.5,
#fontface = "bold",
lineheight = 0.9
) +
scale_color_identity() +
scale_fill_gradient2(
low = "black", mid = "#F7F7F7", high = "orange",
midpoint = 0,
limits = c(-0.6, 0.6),
breaks = seq(-0.6, 0.6, by = 0.3),
#name = "Loading",
name = NULL,
guide = guide_colorbar(
reverse = FALSE,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(0.5, "cm"), # Width/thickness of the vertical bar
barheight = unit(21, "cm"),
ticks = TRUE,
ticks.colour = "grey20",
frame.colour = "grey60",
frame.linewidth = 0.5
)
) +
facet_wrap(~ Component, nrow = 1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = NULL, y = NULL) +
theme_minimal(base_size = 20) +
theme(
strip.text = element_text(face = "bold", size = 23, lineheight = 1.1),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
panel.spacing = unit(0.8, "lines"),
legend.position = "right",
legend.box = "horizontal",
legend.title = element_blank(),
#legend.title = element_text(size = 18, face = "bold", margin = margin(b = 5)),
legend.text = element_text(size = 20)
)
dev.off()## quartz_off_screen
## 2Behaviour: families
Specify families
Step 1: Convert the cross-correlations averages to absolute values - to reflect “strength of coordination” instead of direction.
# absolute values of cross-correlations?
cc_abs = 1
DF_pca_analysis$INTER_smile_av_sync.abs = abs(DF_pca_analysis$INTER_smile_av_sync)
DF_pca_analysis$INTER_laughter_av_sync.abs = abs(DF_pca_analysis$INTER_laughter_av_sync)
DF_pca_analysis$INTER_expressivity_av_sync.abs = abs(DF_pca_analysis$INTER_expressivity_av_sync)
DF_pca_analysis$INTER_volume.abs = abs(DF_pca_analysis$INTER_volume)
DF_pca_analysis$INTER_distance_av_sync.abs = abs(DF_pca_analysis$INTER_distance_av_sync)
DF_pca_analysis$INTER_headmov_av_sync.abs = abs(DF_pca_analysis$INTER_headmov_av_sync)Step 2: Create theory-driven families of behaviours:
family_vars <- list(
intonation = c(
"P_INTRA_intonation_wiggliness",
"P_INTRA_intonation_Spaciousness",
"P_INTRA_intonation_SDf0",
"P_INTRA_intonation_maxmin_f0",
#"P_INTRA_intonation_trimmed_f0",
"P_INTRA_volume_norm"
),
turn_taking = c(
"INTER_n_turn_transitions",
"INTER_mean_turn_length_diff",
"INTER_total_turn_length_diff",
"INTER_total_speaking_time_diff",
"INTER_mean_FTO",
"INTER_mean_FTO_diff"
),
backchannel_tpt = c(
"INTER_P_pure_bc_normalised_rate",
"INTER_P_change_bc_normalised_rate",
"INTER_P_success_tpt_normalised_rate",
"INTER_P_fail_tpt_normalised_rate",
"INTER_nods" #
),
expressivity = c(
"P_INTRA_per_smile",
"P_INTRA_per_laughter",
"P_INTRA_m_expressivity",
"P_INTRA_n_nods",
"P_INTRA_headmov"
),
synchrony_expressivity = c(
ifelse(cc_abs==1,"INTER_smile_av_sync.abs","INTER_smile_av_sync"),
ifelse(cc_abs==1,"INTER_laughter_av_sync.abs","INTER_laughter_av_sync"),
ifelse(cc_abs==1,"INTER_expressivity_av_sync.abs","INTER_expressivity_av_sync"),
ifelse(cc_abs==1,"INTER_volume.abs","INTER_volume"),
ifelse(cc_abs==1,"INTER_headmov_av_sync.abs","INTER_headmov_av_sync")
)
# movement_distance = c(
# "P_INTRA_distance_to_camera"
# #"P_INTRA_headmov"
# #"P_INTRA_n_nods"
# ),
# synchrony_distance = c(
# ifelse(cc_abs==1,"INTER_distance_av_sync.abs","INTER_distance_av_sync")
# #ifelse(cc_abs==1,"INTER_headmov_av_sync.abs","INTER_headmov_av_sync")
# #"INTER_nods"
# )
)
all_family_predictors <- unlist(family_vars)print(family_vars$synchrony_expressivity)## [1] "INTER_smile_av_sync.abs" "INTER_laughter_av_sync.abs" "INTER_expressivity_av_sync.abs"
## [4] "INTER_volume.abs" "INTER_headmov_av_sync.abs"print(family_vars$turn_taking)## [1] "INTER_n_turn_transitions" "INTER_mean_turn_length_diff" "INTER_total_turn_length_diff"
## [4] "INTER_total_speaking_time_diff" "INTER_mean_FTO" "INTER_mean_FTO_diff"print(family_vars$backchannel_tpt)## [1] "INTER_P_pure_bc_normalised_rate" "INTER_P_change_bc_normalised_rate" "INTER_P_success_tpt_normalised_rate"
## [4] "INTER_P_fail_tpt_normalised_rate" "INTER_nods"print(family_vars$expressivity)## [1] "P_INTRA_per_smile" "P_INTRA_per_laughter" "P_INTRA_m_expressivity" "P_INTRA_n_nods" "P_INTRA_headmov"print(family_vars$intonation)## [1] "P_INTRA_intonation_wiggliness" "P_INTRA_intonation_Spaciousness" "P_INTRA_intonation_SDf0"
## [4] "P_INTRA_intonation_maxmin_f0" "P_INTRA_volume_norm"Step 3: Merge with structure information.
# Bring the trait PCs
# One clean DF from self-PC in traits
DF_pca <- cbind(
DF_pca_analysis[, c(
"dyad_indx", "ID", "partner_ID", "dyad", "Group",
"R_Int_rapport","R_Int_willingness"
)],
Q_PCAs
)
DF_pca <- DF_pca %>%
dplyr::rename(
S_Q_PC1 = Q_PC1,
S_Q_PC2 = Q_PC2,
S_Q_PC3 = Q_PC3
)
# Add partner PCs from the other row in the same dyad:
partner_pcs <- DF_pca %>%
dplyr::transmute(
dyad_indx,
match_ID = ID,
P_Q_PC1 = S_Q_PC1,
P_Q_PC2 = S_Q_PC2,
P_Q_PC3 = S_Q_PC3
)
DF_pca <- DF_pca %>%
dplyr::left_join(
partner_pcs,
by = c("dyad_indx", "partner_ID" = "match_ID")
)
# build dataset with all needed columns
dat_all <- cbind(
DF_pca,
DF_pca_analysis[, all_family_predictors]
)Step 4: Z-score raw behaviour variables.
dat_all_z <- dat_all
dat_all_z[all_family_predictors] <- scale(dat_all_z[all_family_predictors])Step 5: Create family composite scores.
for (fam in names(family_vars)) {
dat_all_z[[paste0("fam_", fam)]] <- rowMeans(
dat_all_z[, family_vars[[fam]], drop = FALSE],
na.rm = FALSE
)
}
family_score_vars <- paste0("fam_", names(family_vars))
S_trait_pc_vars <- c(
"S_Q_PC1", "S_Q_PC2", "S_Q_PC3"
)
P_trait_pc_vars <- c(
"P_Q_PC1", "P_Q_PC2", "P_Q_PC3"
)
# Merge all
analysis_predictors <- c(family_score_vars, S_trait_pc_vars, P_trait_pc_vars)
# keep complete cases for analysis <- there are no missing cases here!
dat_analysis_cc <- na.omit(dat_all_z[, c(
"dyad_indx", "ID", "partner_ID", "dyad", "Group",
"R_Int_rapport", "R_Int_willingness",
analysis_predictors
)])
# make sure grouping vars are factors
dat_analysis_cc$dyad <- factor(dat_analysis_cc$dyad)
dat_analysis_cc$Group <- factor(dat_analysis_cc$Group)
dat_analysis_cc$ID <- factor(dat_analysis_cc$ID)
dat_analysis_cc$partner_ID <- factor(dat_analysis_cc$partner_ID)
dat_analysis_cc$dyad_indx <- factor(dat_analysis_cc$dyad_indx)Step 6: Standardize final predictors for ridge comparability:
dat_analysis_cc[analysis_predictors] <- scale(dat_analysis_cc[analysis_predictors])Correlation plots
Corrplot with all raw variables:
pca_vars <- as.character(all_family_predictors)
pca_vars <- pca_vars[c(1:5, 17:21, 6:16, 22:26)] # reorder to have INTRA fams first, then INTER
cor_matrix <- cor(DF_pca_analysis[, pca_vars], use = "complete.obs")
friendly_names <- c(
P_INTRA_intonation_wiggliness = "Intonation: wiggliness (intra)",
P_INTRA_intonation_Spaciousness = "Intonation: spaciousness (intra)",
P_INTRA_intonation_SDf0 = "Intonation: SD f0 (intra)",
P_INTRA_intonation_maxmin_f0 = "Intonation: f0 range (max–min) (intra)",
P_INTRA_volume_norm = "Volume (normalised) (intra)",
P_INTRA_per_smile = "Smiling (%) (intra)",
P_INTRA_per_laughter = "Laughter (%) (intra)",
P_INTRA_m_expressivity = "Mean expressivity (intra)",
P_INTRA_distance_to_camera = "Distance to camera (intra)",
P_INTRA_headmov = "Head movement (intra)",
P_INTRA_n_nods = "Nods (count) (intra)",
INTER_volume.abs = "Volume synchrony (inter)",
INTER_mean_FTO = "Mean FTO (inter)",
INTER_n_turn_transitions = "Turn transitions (count) (inter)",
INTER_P_pure_bc_normalised_rate = "Backchannels: pure norm. rate (inter)",
INTER_P_change_bc_normalised_rate = "Backchannels: change norm. rate (inter)",
INTER_P_success_tpt_normalised_rate = "Turn-taking: success norm. rate (inter)",
INTER_P_fail_tpt_normalised_rate = "Turn-taking: fail norm. rate (inter)",
INTER_mean_turn_length_diff = "Turn length difference (mean) (inter)",
INTER_total_turn_length_diff = "Turn length difference (total) (inter)",
INTER_mean_FTO_diff = "FTO difference (mean) (inter)",
INTER_total_speaking_time_diff = "Speaking time difference (total) (inter)",
INTER_smile_av_sync = "Smile synchrony (inter)",
INTER_laughter_av_sync = "Laughter synchrony (inter)",
INTER_expressivity_av_sync = "Expressivity synchrony (inter)",
INTER_distance_av_sync = "Distance synchrony (inter)",
INTER_headmov_av_sync = "Head-movement synchrony (inter)",
INTER_nods = "Nods (count normalised) (inter)",
INTER_smile_av_sync.abs = "Smile synchrony (inter)",
INTER_laughter_av_sync.abs = "Laughter synchrony (inter)",
INTER_expressivity_av_sync.abs = "Expressivity synchrony (inter)",
INTER_headmov_av_sync.abs = "Head-movement synchrony (inter)"
)
cor_matrix_nice <- cor_matrix
colnames(cor_matrix_nice) <- friendly_names[colnames(cor_matrix_nice)]
rownames(cor_matrix_nice) <- friendly_names[rownames(cor_matrix_nice)]
cp <- corrplot(cor_matrix_nice,
method = "color",
order = "original",
tl.cex = 0.6,
tl.col = "black",
col = colorRampPalette(c("blue", "white", "red"))(200),
addgrid.col = NA)
# Add squares marking the families
vars_in_plot <- colnames(cor_matrix)
family_vars_present <- lapply(family_vars, function(x) x[x %in% vars_in_plot])
family_names_nice <- lapply(family_vars_present, function(x) friendly_names[x])
rects <- do.call(
rbind,
lapply(family_names_nice, function(x) {
c(x[1], x[length(x)], x[length(x)], x[1])
})
)
corrRect(corrRes = cp, namesMat = rects, col = "black", lwd = 2)
Correlation matrix of behavioural predictors. Rectangles mark the predefined behavioural families, and colours represent the direction and magnitude of pairwise correlations.
### Save
# Open a high-resolution graphics device (12x12 inches)
png("./figures/correlation_matrix_high_res.png", width = 12, height = 12, units = "in", res = 300)
# Draw the plot
cp <- corrplot(
cor_matrix_nice,
method = "color",
order = "original",
tl.cex = 1.5, # Increased slightly since we are drawing on a larger canvas
tl.col = "black",
tl.offset = 0.5, # Pulls the text labels slightly closer to the colored squares
cl.cex = 1.2,
mar = c(0, 0, 0, 0), # Removes the hardcoded outer white margins
#col = colorRampPalette(c("blue", "white", "red"))(200),
col = colorRampPalette(c("lightblue", "white", "black"))(200),
addgrid.col = NA
)
# Add the rectangles
corrRect(corrRes = cp, namesMat = rects, col = "black", lwd = 3) # Increased line width to 3 for visibility
dev.off()## quartz_off_screen
## 2Correlation plot between behaviour families:
friendly_family_names <- c(
fam_expressivity = "INTRA: Expressivity",
fam_intonation = "INTRA: Intonation",
fam_turn_taking = "INTER: Turn-taking",
fam_backchannel_tpt = "INTER: Backchannelling",
fam_synchrony_expressivity = "INTER: Non-verbal synchrony"
)
#summary(dat_analysis_cc[, family_score_vars])
# reorder variables to match the order in friendly_family_names
family_score_vars_ordered <- names(friendly_family_names)
cor_matrix <- cor(
dat_analysis_cc[, family_score_vars_ordered],
use = "complete.obs"
)
cor_matrix_nice <- cor_matrix
colnames(cor_matrix_nice) <- friendly_family_names[colnames(cor_matrix_nice)]
rownames(cor_matrix_nice) <- friendly_family_names[rownames(cor_matrix_nice)]
cp <- corrplot(
cor_matrix_nice,
method = "color",
order = "original",
tl.cex = 0.6,
tl.col = "black",
#col = colorRampPalette(c("blue", "white", "red"))(200),
col = colorRampPalette(c("blue", "white", "red"))(200),
addgrid.col = NA
)
Correlation matrix of family-level behavioural scores. Colours indicate the direction and strength of pairwise correlations between the predefined behavioural families.
# Save
png("./figures/correlation_families_matrix_high_res.png", width = 12, height = 12, units = "in", res = 300)
cp <- corrplot(
cor_matrix_nice,
method = "color",
order = "original",
tl.cex = 2, # Increased slightly since we are drawing on a larger canvas
tl.col = "black",
tl.offset = 0.5, # Pulls the text labels slightly closer to the colored squares
cl.cex = 1.8,
mar = c(0, 0, 0, 0), # Removes the hardcoded outer white margins
col = colorRampPalette(c("pink", "white", "black"))(200),
#col = colorRampPalette(c("blue", "white", "red"))(200),
addgrid.col = NA
)
dev.off()## quartz_off_screen
## 2DF_singular_predictors = DF_pca_analysis
DF_reduced_predictors = dat_analysis_ccANALYSIS 1: Predictors of RAPPORT
Rapport across dyads
model_rapport <- lmer(R_Int_rapport ~ dyad + (1|ID) + (1|partner_ID) + (1|dyad_indx), data = DF_singular_predictors, REML = F)
anova(model_rapport)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## dyad 2144.1 1072 2 84.887 5.0763 0.008271 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1report(anova(model_rapport))## The ANOVA suggests that:
##
## - The main effect of dyad is statistically significant and medium (F(2) = 5.08, p = 0.008; Eta2 (partial) = 0.11, 95% CI [0.02,
## 1.00])
##
## Effect sizes were labelled following Field's (2013) recommendations.Pair-wise comparisons:
emm <- emmeans(model_rapport, ~ dyad)
contr <- pairs(emm, adjust = "holm")
t <- as.data.frame(summary(contr, infer = c(TRUE, TRUE)))
t <- t %>%
select(contrast, estimate, SE, df, t.ratio, p.value) %>%
rename(
Contrast = contrast,
`Est.` = estimate,
`Std. Error` = SE,
`df` = df,
`t value` = t.ratio,
`p value` = p.value
)
d_vals <- get_effect_sizes(model_rapport) %>% filter(Term == "dyad") %>% select(contrast, Std_Estimate)
t <- left_join(t, d_vals, by = c("Contrast" = "contrast"))
names(t)[ncol(t)] <- "Cohen's d"
t %>%
kbl(caption = 'Pairwise comparisons for the effect of dyad type on rapport ratings (with Holm correction)', digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = F, position = "left") %>%
pack_rows("Mixed vs. autistic dyads", 1, 1) %>%
pack_rows("Neurotypical vs. autistic dyads", 2, 2) %>%
pack_rows("Neurotypical vs mixed dyads", 3, 3)| Contrast | Est. | Std. Error | df | t value | p value | Cohen’s d |
|---|---|---|---|---|---|---|
| Mixed vs. autistic dyads | ||||||
| AUT - MIX | 3.240 | 3.964 | 110.556 | 0.817 | 0.415 | 0.22 |
| Neurotypical vs. autistic dyads | ||||||
| AUT - NAUT | -8.846 | 5.475 | 114.291 | -1.616 | 0.218 | -0.61 |
| Neurotypical vs mixed dyads | ||||||
| MIX - NAUT | -12.086 | 3.964 | 110.556 | -3.049 | 0.009 | -0.83 |
Plot:
Rapport_summary <- DFoutcomes %>%
dplyr::group_by(dyad, Group) %>% # Ensure Group is included
dplyr::summarise(
mean = mean(R_Int_rapport, na.rm = TRUE),
lower_se = mean - sd(R_Int_rapport, na.rm = TRUE) / sqrt(n()),
upper_se = mean + sd(R_Int_rapport, na.rm = TRUE) / sqrt(n())
)
p.rapport_dyad_type <- ggplot(data = DFoutcomes, aes(y = R_Int_rapport, x = dyad, fill = Group)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .5) +
geom_point(aes(y = R_Int_rapport, color = Group), position = position_jitter(width = .15), size = 2, alpha = 0.8) +
geom_boxplot(width = .1, outlier.shape = NA, alpha = 0.5, fill = 'white', colour = 'black') +
geom_errorbar(data = Rapport_summary, aes(ymin = lower_se, ymax = upper_se, y = mean), position = position_nudge(x = 0.3), width = 0.1) +
geom_point(data = Rapport_summary, aes(x = dyad, y = mean, color = Group), position = position_nudge(x = 0.3), size = 1.1) +
#expand_limits(x = 5.25) +
#guides(fill = FALSE) +
#guides(color = FALSE) +
scale_color_manual(values = colours_neurotypes) +
scale_fill_manual(values = colours_neurotypes) +
# scale_color_manual(values = colours_dyads[c(1,3)]) +
# scale_fill_manual(values = colours_dyads[c(1,3)]) +
#scale_x_discrete(labels = c("AUT" = "AUT", "MIX" = "MIX", "NT" = "NAUT")) +
xlab("Dyad type\n") +
ylab("Rapport") +
labs(fill = "Participant", colour = "Participant") +
#facet_wrap(.~group, strip.position = "bottom") +
#coord_flip() +
my_theme +
theme(axis.title.x = element_text(size = 25),
axis.title.y = element_text(size = 25),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20))
# Map factor levels to numeric positions
level_map <- setNames(1:3, levels(DFoutcomes$dyad)) # Creates: AUT=1, MIX=2, NT=3
pvalues_df <- t %>%
dplyr::mutate(
group1 = sub(" - .*", "", Contrast),
group2 = sub(".* - ", "", Contrast),
xmin = as.numeric(level_map[group1]), # Convert to numeric position
xmax = as.numeric(level_map[group2]), # Convert to numeric position
p.signif = case_when(
`p value` <= 0.001 ~ "***",
`p value` <= 0.01 ~ "**",
`p value` <= 0.05 ~ "*",
TRUE ~ "ns"
)
)
# Add staggered y positions
max_y <- max(DFoutcomes$R_Int_rapport, na.rm = TRUE)
pvalues_df <- pvalues_df %>%
dplyr::mutate(y.position = c(max_y * 1.05, max_y * 1.12, max_y * 1.19))
# Now plot with numeric xmin/xmax
p.rapport <- p.rapport_dyad_type +
stat_pvalue_manual(
pvalues_df,
label = "p.signif",
xmin = "xmin", # Now refers to numeric column
xmax = "xmax", # Now refers to numeric column
y.position = "y.position",
tip.length = 0.01,
size = 4,
inherit.aes = FALSE
)
p.rapport
Performance summary:
performance::model_performance(model_rapport)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ----------------------------------------------------------------------------
## 1850.3 | 1850.9 | 1873.7 | 0.578 | 0.049 | 0.556 | 10.902 | 14.532See BF:
model_null <- update(model_rapport, . ~ . - dyad)
compare_models_BF(model_null = model_null,
model_full = model_rapport,
model_null_name = "Intercept-only",
model_full_name = "With dyad")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Intercept-only model BIC: 1885.11
## - With dyad model BIC: 1886.21
## - BIC difference (Full - Null): 1.1
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Intercept-only vs With dyad ): 1.7
## - BF₁₀ (Evidence for With dyad vs Intercept-only ): 0.58
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Intercept-only
## - Evidence strength: Weak
## - Bayes Factor: 1.73 :1 in favor of Intercept-only
##
## CONCLUSION: There is weak evidence ( 1.73 :1) in favor
## of the Intercept-only model over the With dyad model.
## =======================================================================See effect sizes:
get_anova_effect_sizes(model_rapport)## Parameter Eta2_partial CI CI_low CI_high Omega2_partial Omega2_CI_low Omega2_CI_high
## 1 dyad 0.11 0.95 0.02 1 0.08 0.01 1get_effect_sizes(model_rapport)## contrast Term Type Unstd_Estimate Unstd_SE Unstd_CI_low Unstd_CI_high p_value Std_Estimate Std_CI_low Std_CI_high
## 1 AUT - MIX dyad factor_contrast 3.24 3.96 -6.18 12.66 0.69 0.22 -0.43 0.87
## 2 AUT - NAUT dyad factor_contrast -8.85 5.48 -21.85 4.16 0.24 -0.61 -1.50 0.29
## 3 MIX - NAUT dyad factor_contrast -12.09 3.96 -21.50 -2.67 0.01 -0.83 -1.48 -0.18So there is an effect of dyad type on repport, so that NT dyads report more rapport than the ones with at least one autistic partner. BF suggest that it’s so small, there is even negligible evidence for the null model. Now the question is: where is this difference coming from? Since participants had only just met and interacted briefly, the source of this difference has to be in the behaviour. We then explore different behaviours to first find general predictors or rapport and then to see which of those differ between the neurotypes.
SRM Model
We first fit a model with no predictors but with 4 random effects:
(1|ID) = random intercept for subject, i.e., the actor effect (1|partner_ID) = random intercept for the interaction partner of the subject, i.e., the partner effect (1|dyad_indx) = random intercept for the dyad, i.e., the relation effect (1|session) = random intercept for the experimental session (in case the groups in different sessions differed systematically)
srm_model_null <- lmer(R_Int_rapport ~ 1 + (1|ID) + (1|partner_ID) + (1|dyad_indx), data = DF_singular_predictors, REML = T)
srm_model_null@call## lmer(formula = R_Int_rapport ~ 1 + (1 | ID) + (1 | partner_ID) +
## (1 | dyad_indx), data = DF_singular_predictors, REML = T)Variance:
VarCorr(srm_model_null)## Groups Name Std.Dev.
## dyad_indx (Intercept) 11.2194
## partner_ID (Intercept) 8.7118
## ID (Intercept) 9.6240
## Residual 14.5079Let’s look at the explained variance in the model:
MuMIn::r.squaredGLMM(srm_model_null)## R2m R2c
## [1,] 0 0.5831028v = MuMIn::r.squaredGLMM(srm_model_null)
v = round(v[2],2)*100
v## [1] 58The random effects explained 58% of variance in the model. We use this SRM to estimate how much variation in rapport is a function of the individuals (autistic or neurotypical) in the role of the actor (the one rating rapport) and of the partner (the one with whom the actor is rating the rapport), and how much of the unique interactions across dyads. Let’s look at that in detail.
# Look at variance per actor, partner, and dyad
variance_components <- as.data.frame(VarCorr(srm_model_null))
variance_components$variance <- variance_components$sdcor^2
#variance_components
variances <- as.data.frame(variance_components)$vcov
total_variance <- sum(variances[1:3])
variance_proportions <- variances / total_variance
variance_table <- data.frame(
Component = variance_components$grp,
Variance = variances,
Proportion = round(variance_proportions,3)
)
print(variance_table[1:3,])## Component Variance Proportion
## 1 dyad_indx 125.87390 0.428
## 2 partner_ID 75.89522 0.258
## 3 ID 92.62168 0.315actor_v = round(variance_table$Proportion[variance_table$Component=="ID"],2)*100
partner_v = round(variance_table$Proportion[variance_table$Component=="partner_ID"],2)*100
dyad_v = round(variance_table$Proportion[variance_table$Component=="dyad_indx"],2)*100
SRM_variances_p1 = paste("Rapport: ","Actor: ",actor_v,"%. Partner: ",partner_v,"%. Dyad: ",dyad_v,"%",sep="")Out of the 58% variance explained by the random effects:
actor effects account for
ractor_v% of the variance, *partner effects account forr
partner_v% of the variance, dyad effects account for
rdyad_v`% of the variance.
Rapport by neurotype in MIX
model_rapport_MIX <- lmer(R_Int_rapport ~ Group + (1|ID) + (1|partner_ID),
data = DF_singular_predictors[DF_singular_predictors$dyad=="MIX",],
REML = F)
model_rapport_MIX@call## lmer(formula = R_Int_rapport ~ Group + (1 | ID) + (1 | partner_ID),
## data = DF_singular_predictors[DF_singular_predictors$dyad ==
## "MIX", ], REML = F)a <- anova(model_rapport_MIX)
report(a)## The ANOVA suggests that:
##
## - The main effect of Group is statistically not significant and very small (F(1) = 2.23e-03, p = 0.963; Eta2 (partial) = 4.35e-05,
## 95% CI [0.00, 1.00])
##
## Effect sizes were labelled following Field's (2013) recommendations.Performance summary:
performance::model_performance(model_rapport_MIX)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -------------------------------------------------------------------------
## 940.4 | 941.0 | 953.6 | 0.579 | 3.297e-05 | 0.579 | 10.740 | 15.046See BF:
model_null <- update(model_rapport_MIX, . ~ . - Group)
compare_models_BF(model_null = model_null,
model_full = model_rapport_MIX,
model_null_name = "Intercept-only",
model_full_name = "With neurotype")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Intercept-only model BIC: 958.19
## - With neurotype model BIC: 962.83
## - BIC difference (Full - Null): 4.64
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Intercept-only vs With neurotype ): 10
## - BF₁₀ (Evidence for With neurotype vs Intercept-only ): 0.098
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Intercept-only
## - Evidence strength: Strong
## - Bayes Factor: 10.19 :1 in favor of Intercept-only
##
## CONCLUSION: There is strong evidence ( 10.19 :1) in favor
## of the Intercept-only model over the With neurotype model.
## =======================================================================See effect sizes:
get_anova_effect_sizes(model_rapport_MIX)## Parameter Eta2_partial CI CI_low CI_high Omega2_partial Omega2_CI_low Omega2_CI_high
## 1 Group 0 0.95 0 1 0 0 1get_effect_sizes(model_rapport_MIX)## contrast Term Type Unstd_Estimate Unstd_SE Unstd_CI_low Unstd_CI_high p_value Std_Estimate Std_CI_low Std_CI_high
## 1 AUT - NAUT Group factor_contrast 0.26 5.73 -11.21 11.74 0.96 0.02 -0.74 0.78Predictors screening (Ridge)
Step 1: Build the analysis dataset
analysis_predictors <- c(family_score_vars,S_trait_pc_vars)
DFrapp = DF_reduced_predictors[,c(1:7,which(names(DF_reduced_predictors) %in% analysis_predictors))]Step 2: Ridge as a predictor screening procedure
Build ridge design matrix (all predictors + all predictor:neurotype interactions):
predictor_main_vars <- analysis_predictors
ridge_formula <- as.formula(
paste(
"R_Int_rapport ~ (",
paste(predictor_main_vars, collapse = " + "),
") * Group"
)
)
x_ridge <- model.matrix(ridge_formula, data = DFrapp)[, -1] # drop intercept
y_ridge <- DFrapp$R_Int_rapport
cat("Design matrix dimensions:", dim(x_ridge), "\n")## Design matrix dimensions: 199 17cat("Number of columns in design matrix:", ncol(x_ridge), "\n\n")## Number of columns in design matrix: 17First 20 column names:
print(colnames(x_ridge)[1:min(20, ncol(x_ridge))])## [1] "fam_intonation" "fam_turn_taking" "fam_backchannel_tpt"
## [4] "fam_expressivity" "fam_synchrony_expressivity" "S_Q_PC1"
## [7] "S_Q_PC2" "S_Q_PC3" "GroupNAUT"
## [10] "fam_intonation:GroupNAUT" "fam_turn_taking:GroupNAUT" "fam_backchannel_tpt:GroupNAUT"
## [13] "fam_expressivity:GroupNAUT" "fam_synchrony_expressivity:GroupNAUT" "S_Q_PC1:GroupNAUT"
## [16] "S_Q_PC2:GroupNAUT" "S_Q_PC3:GroupNAUT"Columns containing dyad interactions:
print(grep(":", colnames(x_ridge), value = TRUE))## [1] "fam_intonation:GroupNAUT" "fam_turn_taking:GroupNAUT" "fam_backchannel_tpt:GroupNAUT"
## [4] "fam_expressivity:GroupNAUT" "fam_synchrony_expressivity:GroupNAUT" "S_Q_PC1:GroupNAUT"
## [7] "S_Q_PC2:GroupNAUT" "S_Q_PC3:GroupNAUT"Step 3: Grouped folds by ID
Fold sizes:
id_ridge <- as.character(DFrapp$ID)
unique_ids_ridge <- unique(id_ridge)
nfolds <- 10
foldid_by_id <- sample(rep(1:nfolds, length.out = length(unique_ids_ridge)))
names(foldid_by_id) <- unique_ids_ridge
foldid_ridge <- unname(foldid_by_id[id_ridge])
print(table(foldid_ridge))## foldid_ridge
## 1 2 3 4 5 6 7 8 9 10
## 21 20 22 21 22 20 19 15 20 19Step 4: Grouped cross-validation (CV) Ridge
cv_ridge <- glmnet::cv.glmnet(
x = x_ridge,
y = y_ridge,
alpha = 0,
family = "gaussian",
standardize = TRUE,
foldid = foldid_ridge
)
coef_min <- as.matrix(coef(cv_ridge, s = "lambda.min"))
coef_1se <- as.matrix(coef(cv_ridge, s = "lambda.1se"))
coef_min_df <- data.frame(
variable = rownames(coef_min),
coef = as.numeric(coef_min[, 1])
)
coef_1se_df <- data.frame(
variable = rownames(coef_1se),
coef = as.numeric(coef_1se[, 1])
)
coef_min_df <- coef_min_df[order(abs(coef_min_df$coef), decreasing = TRUE), ]
coef_1se_df <- coef_1se_df[order(abs(coef_1se_df$coef), decreasing = TRUE), ]Top coefficients at lambda.min:
print(head(coef_min_df, 20))## variable coef
## 1 (Intercept) 66.84088299
## 7 S_Q_PC1 -1.17525537
## 16 S_Q_PC1:GroupNAUT -0.94345101
## 9 S_Q_PC3 0.90434312
## 6 fam_synchrony_expressivity 0.87621884
## 15 fam_synchrony_expressivity:GroupNAUT 0.72845198
## 8 S_Q_PC2 0.69853838
## 10 GroupNAUT 0.65275424
## 17 S_Q_PC2:GroupNAUT 0.58052358
## 2 fam_intonation -0.42666740
## 4 fam_backchannel_tpt 0.29383719
## 3 fam_turn_taking 0.24526873
## 14 fam_expressivity:GroupNAUT -0.20943382
## 5 fam_expressivity 0.18587480
## 11 fam_intonation:GroupNAUT -0.16811087
## 13 fam_backchannel_tpt:GroupNAUT 0.13169529
## 12 fam_turn_taking:GroupNAUT 0.11122702
## 18 S_Q_PC3:GroupNAUT 0.06813819Top coefficients at lambda.1se:
print(head(coef_1se_df, 20))## variable coef
## 1 (Intercept) 67.4422110552763882651561289094388484954833984
## 16 S_Q_PC1:GroupNAUT -0.0000000000000000000000000000000000051147960
## 7 S_Q_PC1 -0.0000000000000000000000000000000000050017237
## 10 GroupNAUT 0.0000000000000000000000000000000000043495562
## 6 fam_synchrony_expressivity 0.0000000000000000000000000000000000033291527
## 9 S_Q_PC3 0.0000000000000000000000000000000000031548305
## 15 fam_synchrony_expressivity:GroupNAUT 0.0000000000000000000000000000000000030253163
## 8 S_Q_PC2 0.0000000000000000000000000000000000027289114
## 17 S_Q_PC2:GroupNAUT 0.0000000000000000000000000000000000025195578
## 2 fam_intonation -0.0000000000000000000000000000000000015273884
## 4 fam_backchannel_tpt 0.0000000000000000000000000000000000014365799
## 13 fam_backchannel_tpt:GroupNAUT 0.0000000000000000000000000000000000011269396
## 11 fam_intonation:GroupNAUT -0.0000000000000000000000000000000000009009670
## 5 fam_expressivity 0.0000000000000000000000000000000000008496410
## 18 S_Q_PC3:GroupNAUT 0.0000000000000000000000000000000000008051178
## 3 fam_turn_taking 0.0000000000000000000000000000000000006280321
## 14 fam_expressivity:GroupNAUT -0.0000000000000000000000000000000000001547509
## 12 fam_turn_taking:GroupNAUT 0.0000000000000000000000000000000000001521861The top results include interactions, which confirms that the rapport predictors will likely be different between neurotypes. But this is just one ridge fit. Let’s check the stability with bootstrapping.
Step 5: Bootstrapping
Grouped bootstrap + grouped CV inside bootstrap for the full ridge screening model:
B <- 500
nfolds <- 10
unique_ids_boot <- unique(as.character(DFrapp$ID))
boot_coef_mat <- matrix(NA, nrow = B, ncol = ncol(x_ridge))
colnames(boot_coef_mat) <- colnames(x_ridge)
boot_info <- data.frame(
iter = 1:B,
n_rows = NA_integer_,
n_id_boot = NA_integer_,
lambda_min = NA_real_,
lambda_1se = NA_real_,
cv_ok = FALSE,
error_msg = NA_character_
)
for (b in 1:B) {
boot_ids <- sample(unique_ids_boot, size = length(unique_ids_boot), replace = TRUE)
boot_list <- lapply(seq_along(boot_ids), function(i) {
id_now <- boot_ids[i]
tmp <- DFrapp[as.character(DFrapp$ID) == id_now, , drop = FALSE]
tmp$ID_boot <- paste0(id_now, "_rep", i)
tmp
})
boot_dat <- do.call(rbind, boot_list)
boot_info$n_rows[b] <- nrow(boot_dat)
boot_info$n_id_boot[b] <- length(unique(boot_dat$ID_boot))
x_boot <- model.matrix(ridge_formula, data = boot_dat)[, -1, drop = FALSE]
y_boot <- boot_dat$R_Int_rapport
id_boot <- as.character(boot_dat$ID_boot)
unique_id_boot <- unique(id_boot)
foldid_by_id_boot <- sample(rep(1:nfolds, length.out = length(unique_id_boot)))
names(foldid_by_id_boot) <- unique_id_boot
foldid_boot <- unname(foldid_by_id_boot[id_boot])
cv_boot <- tryCatch(
glmnet::cv.glmnet(
x = x_boot,
y = y_boot,
alpha = 0,
family = "gaussian",
standardize = TRUE,
foldid = foldid_boot
),
error = function(e) e
)
if (inherits(cv_boot, "error")) {
boot_info$error_msg[b] <- conditionMessage(cv_boot)
next
}
boot_info$cv_ok[b] <- TRUE
boot_info$lambda_min[b] <- cv_boot$lambda.min
boot_info$lambda_1se[b] <- cv_boot$lambda.1se
coef_boot <- as.matrix(coef(cv_boot, s = "lambda.min"))[-1, 1]
boot_coef_mat[b, names(coef_boot)] <- coef_boot
}Summaries:
boot_summary_full <- data.frame(
variable = colnames(boot_coef_mat),
mean_coef = colMeans(boot_coef_mat, na.rm = TRUE),
sd_coef = apply(boot_coef_mat, 2, sd, na.rm = TRUE),
prop_positive = colMeans(boot_coef_mat > 0, na.rm = TRUE),
prop_negative = colMeans(boot_coef_mat < 0, na.rm = TRUE),
mean_abs_coef = colMeans(abs(boot_coef_mat), na.rm = TRUE)
)
boot_summary_full <- boot_summary_full[order(boot_summary_full$mean_abs_coef, decreasing = TRUE), ]
cat("Successful bootstrap CV runs:", sum(boot_info$cv_ok), "out of", B, "\n")## Successful bootstrap CV runs: 500 out of 500cat("Failed runs:", sum(!boot_info$cv_ok), "\n\n")## Failed runs: 0if (sum(!boot_info$cv_ok) > 0) {
cat("Bootstrap errors:\n")
print(sort(table(na.omit(boot_info$error_msg)), decreasing = TRUE))
}Top 25 bootstrap-stable coefficients:
print(head(boot_summary_full, 25))## variable mean_coef sd_coef prop_positive prop_negative
## S_Q_PC1 S_Q_PC1 -3.32478451 2.7122723 0.002 0.998
## GroupNAUT GroupNAUT -0.35498579 3.9249703 0.610 0.390
## S_Q_PC3 S_Q_PC3 2.51364152 1.6580501 0.988 0.012
## S_Q_PC3:GroupNAUT S_Q_PC3:GroupNAUT -0.59326332 2.9646597 0.470 0.530
## S_Q_PC1:GroupNAUT S_Q_PC1:GroupNAUT -0.28238434 3.0817414 0.262 0.738
## fam_synchrony_expressivity fam_synchrony_expressivity 2.02333049 1.1881443 0.998 0.002
## fam_expressivity:GroupNAUT fam_expressivity:GroupNAUT -0.85045669 2.2734049 0.364 0.636
## fam_synchrony_expressivity:GroupNAUT fam_synchrony_expressivity:GroupNAUT 1.25756729 1.6789354 0.838 0.162
## S_Q_PC2 S_Q_PC2 1.23799043 1.3575944 0.926 0.074
## S_Q_PC2:GroupNAUT S_Q_PC2:GroupNAUT 1.22326266 1.9057239 0.848 0.152
## fam_backchannel_tpt:GroupNAUT fam_backchannel_tpt:GroupNAUT -0.11657210 1.8098598 0.524 0.476
## fam_intonation:GroupNAUT fam_intonation:GroupNAUT -0.02523403 1.7651556 0.442 0.558
## fam_turn_taking:GroupNAUT fam_turn_taking:GroupNAUT 0.23796623 1.5901965 0.548 0.452
## fam_intonation fam_intonation -1.04134413 0.9178973 0.066 0.934
## fam_expressivity fam_expressivity 0.62641078 1.3584565 0.694 0.306
## fam_turn_taking fam_turn_taking 0.60834744 1.2323449 0.750 0.250
## fam_backchannel_tpt fam_backchannel_tpt 0.34704778 1.1693522 0.708 0.292
## mean_abs_coef
## S_Q_PC1 3.3451019
## GroupNAUT 2.5453905
## S_Q_PC3 2.5190333
## S_Q_PC3:GroupNAUT 2.0502200
## S_Q_PC1:GroupNAUT 2.0318264
## fam_synchrony_expressivity 2.0238047
## fam_expressivity:GroupNAUT 1.6916055
## fam_synchrony_expressivity:GroupNAUT 1.5463075
## S_Q_PC2 1.4986429
## S_Q_PC2:GroupNAUT 1.4967393
## fam_backchannel_tpt:GroupNAUT 1.2102224
## fam_intonation:GroupNAUT 1.1466315
## fam_turn_taking:GroupNAUT 1.1308569
## fam_intonation 1.1088689
## fam_expressivity 1.0569423
## fam_turn_taking 0.9639584
## fam_backchannel_tpt 0.8671546The ranking gives us now a basis for which predictors to use in the final mixed models. Again, the strongest and most stable signals are mostly dyad-specific interaction terms, not just pooled main effects: so we can’t screen on the entire sample.
Step 6: Informed feature selection
To select the predictors that will enter the mixed models’ analysis, we will choose those that: - have relatively high mean_abs_coef (we choose 2) - have a strong sign stability (at least 0.8 in any direction) - all main effects of any interaction that fits the rules above
# Cutoff for the stability in either direction
sign_stability_cutoff <- 0.8
# Cutoff for the mean coefficient score
mean_abs_cutoff <- 1.5
boot_screen <- boot_summary_full
boot_screen$sign_stability <- pmax(boot_screen$prop_positive, boot_screen$prop_negative)
boot_screen$is_interaction <- grepl(":", boot_screen$variable)
# keep strong and sign-stable terms
selected_core <- boot_screen[
boot_screen$mean_abs_coef >= mean_abs_cutoff &
boot_screen$sign_stability >= sign_stability_cutoff,
]
# extract main effects corresponding to selected interactions
selected_interactions <- selected_core$variable[selected_core$is_interaction]
interaction_bases <- unique(sub(":.*$", "", selected_interactions))
# keep main effects for any selected interaction, if present in boot table
main_effects_to_add <- boot_screen$variable[
boot_screen$variable %in% interaction_bases
]
# final selected list = core selected + hierarchy-preserving main effects
selected_final_vars <- unique(c(selected_core$variable, main_effects_to_add))
selected_final <- boot_screen[boot_screen$variable %in% selected_final_vars, ]
selected_final <- selected_final[order(selected_final$is_interaction, -selected_final$mean_abs_coef), ]
cat("mean_abs_cutoff =", mean_abs_cutoff, "\n")## mean_abs_cutoff = 1.5cat("sign_stability_cutoff =", sign_stability_cutoff, "\n\n")## sign_stability_cutoff = 0.8Selected core terms:
print(
selected_core[order(selected_core$mean_abs_coef, decreasing = TRUE),
c("variable", "mean_coef", "mean_abs_coef", "prop_positive", "prop_negative", "sign_stability")]
)## variable mean_coef mean_abs_coef prop_positive prop_negative
## S_Q_PC1 S_Q_PC1 -3.324785 3.345102 0.002 0.998
## S_Q_PC3 S_Q_PC3 2.513642 2.519033 0.988 0.012
## fam_synchrony_expressivity fam_synchrony_expressivity 2.023330 2.023805 0.998 0.002
## fam_synchrony_expressivity:GroupNAUT fam_synchrony_expressivity:GroupNAUT 1.257567 1.546308 0.838 0.162
## sign_stability
## S_Q_PC1 0.998
## S_Q_PC3 0.988
## fam_synchrony_expressivity 0.998
## fam_synchrony_expressivity:GroupNAUT 0.838Main effects added to preserve hierarchy:
print(setdiff(main_effects_to_add, selected_core$variable))## character(0)Final selected terms:
print(
selected_final[, c("variable", "mean_coef", "mean_abs_coef", "prop_positive", "prop_negative", "sign_stability", "is_interaction")]
)## variable mean_coef mean_abs_coef prop_positive prop_negative
## S_Q_PC1 S_Q_PC1 -3.324785 3.345102 0.002 0.998
## S_Q_PC3 S_Q_PC3 2.513642 2.519033 0.988 0.012
## fam_synchrony_expressivity fam_synchrony_expressivity 2.023330 2.023805 0.998 0.002
## fam_synchrony_expressivity:GroupNAUT fam_synchrony_expressivity:GroupNAUT 1.257567 1.546308 0.838 0.162
## sign_stability is_interaction
## S_Q_PC1 0.998 FALSE
## S_Q_PC3 0.988 FALSE
## fam_synchrony_expressivity 0.998 FALSE
## fam_synchrony_expressivity:GroupNAUT 0.838 TRUECollapse the ridge terms to variable-level output for the mixed models:
selected_terms <- selected_final$variable
collapse_to_predictor <- function(x) {
# interaction term -> keep only variable before :
if (grepl(":", x)) {
return(sub(":.*$", "", x))
}
# dyad dummy terms -> collapse to dyad
if (grepl("^Group", x)) {
return("Group")
}
# all other main effects stay as they are
return(x)
}
collapsed_predictors <- unique(vapply(selected_terms, collapse_to_predictor, character(1)))
# which predictors were selected as interactions with dyad?
selected_interaction_predictors <- unique(
vapply(
selected_terms[grepl(":", selected_terms)],
collapse_to_predictor,
character(1)
)
)
# final variable-level summary table
selected_predictor_summary <- data.frame(
predictor = collapsed_predictors,
include_main_effect = TRUE,
include_group_interaction = collapsed_predictors %in% selected_interaction_predictors
)
# sort a bit more nicely
selected_predictor_summary$type <- ifelse(
selected_predictor_summary$predictor == "Group", "grouping_variable",
ifelse(grepl("^fam_", selected_predictor_summary$predictor), "behavior_family", "trait_PC")
)
selected_predictor_summary <- selected_predictor_summary[
order(selected_predictor_summary$type, selected_predictor_summary$predictor),
]Final predictor list for mixed models:
print(selected_predictor_summary)## predictor include_main_effect include_group_interaction type
## 3 fam_synchrony_expressivity TRUE TRUE behavior_family
## 1 S_Q_PC1 TRUE FALSE trait_PC
## 2 S_Q_PC3 TRUE FALSE trait_PCBehavior families selected:
print(selected_predictor_summary$predictor[selected_predictor_summary$type == "behavior_family"])## [1] "fam_synchrony_expressivity"Trait PCs selected:
print(selected_predictor_summary$predictor[selected_predictor_summary$type == "trait_PC"])## [1] "S_Q_PC1" "S_Q_PC3"Predictors that should be tested with Group interactions:
print(selected_predictor_summary$predictor[selected_predictor_summary$include_group_interaction])## [1] "fam_synchrony_expressivity"# For plotting later:
boot_summary_rapport <- boot_summary_fullBy NEUROTYPE
Rapport ~ predictors x S_neurotype x P_neurotype (LMM)
Model
Prepare the lmer formula:
DFrapp$partner_neurotype <- ifelse(
DFrapp$dyad == "MIX" & DFrapp$Group == "NAUT", "AUT",
ifelse(
DFrapp$dyad == "MIX" & DFrapp$Group == "AUT", "NAUT",
as.character(DFrapp$Group)
)
)
DFrapp$partner_neurotype <- factor(
DFrapp$partner_neurotype,
levels = c("AUT", "NAUT")
)
main_predictors <- selected_predictor_summary$predictor
interaction_predictors <- main_predictors
main_predictors <- c(main_predictors,"partner_neurotype","Group")
fixed_formula_text <- paste0(
"R_Int_rapport ~ ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":partner_neurotype"), collapse = " + "),
"+",
paste(paste0(interaction_predictors, ":Group"), collapse = " + "),
"+ partner_neurotype:Group"
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_rapport ~ fam_synchrony_expressivity + S_Q_PC1 + S_Q_PC3 +
## partner_neurotype + Group + fam_synchrony_expressivity:partner_neurotype +
## S_Q_PC1:partner_neurotype + S_Q_PC3:partner_neurotype + fam_synchrony_expressivity:Group +
## S_Q_PC1:Group + S_Q_PC3:Group + partner_neurotype:Group +
## (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)Fit the model:
m.rapp.select.ALL <- lmer(
formula = mixed_formula,
data = DFrapp,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Performance
Singular fit?
print(isSingular(m.rapp.select.ALL, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.rapp.select.ALL), residuals(m.rapp.select.ALL), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.rapp.select.ALL), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.rapp.select.ALL), main='')
qqline(residuals(m.rapp.select.ALL))
layout(1)Variance components:
print(VarCorr(m.rapp.select.ALL))## Groups Name Std.Dev.
## dyad_indx (Intercept) 10.100
## partner_ID (Intercept) 7.402
## ID (Intercept) 10.894
## Residual 12.250Performance summary:
performance::model_performance(m.rapp.select.ALL)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ---------------------------------------------------------------------------
## 1717.9 | 1721.3 | 1773.9 | 0.718 | 0.200 | 0.647 | 8.689 | 12.250Interpretation
Anova:
a.m.rapp.select.ALL <- anova(m.rapp.select.ALL)
anova_tab <- as.data.frame(a.m.rapp.select.ALL)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the rapport model with neurotype",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| fam_synchrony_expressivity | 929.348 | 929.348 | 1 | 107.082 | 6.193 | 0.014 |
| S_Q_PC1 | 2004.635 | 2004.635 | 1 | 58.576 | 13.358 | 0.001 |
| S_Q_PC3 | 873.967 | 873.967 | 1 | 90.897 | 5.824 | 0.018 |
| partner_neurotype | 195.011 | 195.011 | 1 | 41.072 | 1.299 | 0.261 |
| Group | 571.085 | 571.085 | 1 | 53.032 | 3.805 | 0.056 |
| fam_synchrony_expressivity:partner_neurotype | 3.327 | 3.327 | 1 | 151.980 | 0.022 | 0.882 |
| S_Q_PC1:partner_neurotype | 62.785 | 62.785 | 1 | 109.290 | 0.418 | 0.519 |
| S_Q_PC3:partner_neurotype | 369.562 | 369.562 | 1 | 96.935 | 2.463 | 0.120 |
| fam_synchrony_expressivity:Group | 44.311 | 44.311 | 1 | 159.290 | 0.295 | 0.588 |
| S_Q_PC1:Group | 369.973 | 369.973 | 1 | 58.043 | 2.465 | 0.122 |
| S_Q_PC3:Group | 287.149 | 287.149 | 1 | 95.287 | 1.913 | 0.170 |
| partner_neurotype:Group | 246.650 | 246.650 | 1 | 113.376 | 1.644 | 0.202 |
Plots:
p_combined_ALL_rapport_neurotype <- plot_lmm_significant(
model = m.rapp.select.ALL,
anova_obj = a.m.rapp.select.ALL,
df = DFrapp,
grouping_var = "partner_neurotype",
n_cols_plotting = 3,
show_data = T,
colour_by = "Group"
)
p_combined_ALL_rapport_neurotype
See effect sizes:
get_anova_effect_sizes(m.rapp.select.ALL)## Parameter Eta2_partial CI CI_low CI_high Omega2_partial Omega2_CI_low Omega2_CI_high
## 1 fam_synchrony_expressivity 0.05 0.95 0.01 1 0.05 0.00 1
## 2 S_Q_PC1 0.19 0.95 0.06 1 0.17 0.05 1
## 3 S_Q_PC3 0.06 0.95 0.01 1 0.05 0.00 1
## 4 partner_neurotype 0.03 0.95 0.00 1 0.01 0.00 1
## 5 Group 0.07 0.95 0.00 1 0.05 0.00 1
## 6 fam_synchrony_expressivity:partner_neurotype 0.00 0.95 0.00 1 0.00 0.00 1
## 7 S_Q_PC1:partner_neurotype 0.00 0.95 0.00 1 0.00 0.00 1
## 8 S_Q_PC3:partner_neurotype 0.02 0.95 0.00 1 0.01 0.00 1
## 9 fam_synchrony_expressivity:Group 0.00 0.95 0.00 1 0.00 0.00 1
## 10 S_Q_PC1:Group 0.04 0.95 0.00 1 0.02 0.00 1
## 11 S_Q_PC3:Group 0.02 0.95 0.00 1 0.01 0.00 1
## 12 partner_neurotype:Group 0.01 0.95 0.00 1 0.01 0.00 1#get_effect_sizes(m.rapp.select.ALL)See BF:
INTER: expressivity
model_null <- update(m.rapp.select.ALL, . ~ . - fam_synchrony_expressivity -
fam_synchrony_expressivity:partner_neurotype -
fam_synchrony_expressivity:Group)
compare_models_BF(model_null = model_null,
model_full = m.rapp.select.ALL,
model_null_name = "Without INTER: expressivity",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without INTER: expressivity model BIC: 1815.97
## - Full model model BIC: 1825.59
## - BIC difference (Full - Null): 9.62
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without INTER: expressivity vs Full model ): 123
## - BF₁₀ (Evidence for Full model vs Without INTER: expressivity ): 0.0082
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without INTER: expressivity
## - Evidence strength: Decisive/Extreme
## - Bayes Factor: 122.66 :1 in favor of Without INTER: expressivity
##
## CONCLUSION: There is decisive/extreme evidence ( 122.66 :1) in favor
## of the Without INTER: expressivity model over the Full model model.
## =======================================================================Q1:
model_null <- update(m.rapp.select.ALL, . ~ . - S_Q_PC1 -
S_Q_PC1:partner_neurotype -
S_Q_PC1:Group)
compare_models_BF(model_null = model_null,
model_full = m.rapp.select.ALL,
model_null_name = "Without autism-ike traits (PC1)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without autism-ike traits (PC1) model BIC: 1823.08
## - Full model model BIC: 1825.59
## - BIC difference (Full - Null): 2.51
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without autism-ike traits (PC1) vs Full model ): 3.5
## - BF₁₀ (Evidence for Full model vs Without autism-ike traits (PC1) ): 0.29
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without autism-ike traits (PC1)
## - Evidence strength: Moderate
## - Bayes Factor: 3.51 :1 in favor of Without autism-ike traits (PC1)
##
## CONCLUSION: There is moderate evidence ( 3.51 :1) in favor
## of the Without autism-ike traits (PC1) model over the Full model model.
## =======================================================================Q3:
model_null <- update(m.rapp.select.ALL, . ~ . - S_Q_PC3 -
S_Q_PC3:partner_neurotype -
S_Q_PC3:Group)
compare_models_BF(model_null = model_null,
model_full = m.rapp.select.ALL,
model_null_name = "Without cognitive engagement (PC3)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without cognitive engagement (PC3) model BIC: 1822.39
## - Full model model BIC: 1825.59
## - BIC difference (Full - Null): 3.2
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without cognitive engagement (PC3) vs Full model ): 5
## - BF₁₀ (Evidence for Full model vs Without cognitive engagement (PC3) ): 0.2
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without cognitive engagement (PC3)
## - Evidence strength: Moderate
## - Bayes Factor: 4.96 :1 in favor of Without cognitive engagement (PC3)
##
## CONCLUSION: There is moderate evidence ( 4.96 :1) in favor
## of the Without cognitive engagement (PC3) model over the Full model model.
## =======================================================================With gender + age + language
We add age difference (difference between the two partners), gender matching (0 vs 1), and languague.
# Add age difference and gender matching
# 1) Collapse DF to one row per person per dyad
person_info <- DFmain %>%
dplyr::select(dyad_indx, ID, age, gender, language) %>%
dplyr::group_by(dyad_indx, ID,language) %>%
dplyr::summarise(
age_self = mean(age, na.rm = TRUE),
gender_self = dplyr::first(gender),
.groups = "drop"
)
# 2) Self-join within dyad to get partner’s info
dyad_pairs <- person_info %>%
dplyr::inner_join(
person_info,
by = c("dyad_indx","language"),
suffix = c("_self", "_partner"),
relationship = "many-to-many"
) %>%
dplyr::filter(ID_self != ID_partner)
# 3) Clean to one row per person per dyad
dyad_info_unique <- dyad_pairs %>%
dplyr::transmute(
dyad_indx,
ID = ID_self,
age_self = age_self_self,
gender_self = gender_self_self,
age_partner = age_self_partner,
gender_partner = gender_self_partner,
language = language
)
# 4) Join into DFrapp
DFrapp <- DFrapp %>%
dplyr::select(
-dplyr::any_of(c(
"age_self", "gender_self",
"age_partner", "gender_partner",
"age_diff", "gender_match"
))
) %>%
dplyr::left_join(
dyad_info_unique,
by = c("dyad_indx", "ID")
) %>%
dplyr::mutate(
age_diff = age_self - age_partner,
gender_match = dplyr::if_else(
!is.na(gender_self) &
!is.na(gender_partner) &
gender_self == gender_partner,
1L,
0L
)
)
# # Quick sanity check
# DFrapp %>%
# dplyr::select(
# dyad_indx, ID, partner_ID,
# age_self, age_partner, age_diff,
# gender_self, gender_partner, gender_match
# ) %>%
# head()
DFrapp$age_diff <- abs(DFrapp$age_diff)Fit the model:
fixed_formula_text <- paste0(
"R_Int_rapport ~ ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":partner_neurotype"), collapse = " + "),"+",
paste(paste0(interaction_predictors, ":Group"), collapse = " + "),"+",
# paste(paste0(interaction_predictors, ":age_diff"), collapse = " + "),"+", # uncomment if with interactions
# paste(paste0(interaction_predictors, ":gender_match"), collapse = " + "),"+",
# paste(paste0(interaction_predictors, ":language"), collapse = " + "),
"+ partner_neurotype:Group"
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + age_diff + gender_match + language + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_rapport ~ fam_synchrony_expressivity + S_Q_PC1 + S_Q_PC3 +
## partner_neurotype + Group + fam_synchrony_expressivity:partner_neurotype +
## S_Q_PC1:partner_neurotype + S_Q_PC3:partner_neurotype + fam_synchrony_expressivity:Group +
## S_Q_PC1:Group + S_Q_PC3:Group + +partner_neurotype:Group +
## age_diff + gender_match + language + (1 | ID) + (1 | partner_ID) +
## (1 | dyad_indx)m.rapp.select.ALL.agegender <- lmer(
formula = mixed_formula,
data = DFrapp,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Singular fit?
print(isSingular(m.rapp.select.ALL.agegender, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.rapp.select.ALL.agegender), residuals(m.rapp.select.ALL.agegender), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.rapp.select.ALL.agegender), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.rapp.select.ALL.agegender), main='')
qqline(residuals(m.rapp.select.ALL.agegender))
layout(1)Variance components:
print(VarCorr(m.rapp.select.ALL.agegender))## Groups Name Std.Dev.
## dyad_indx (Intercept) 9.8830
## partner_ID (Intercept) 7.2403
## ID (Intercept) 10.5370
## Residual 12.3403Performance summary:
performance::model_performance(m.rapp.select.ALL.agegender)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ---------------------------------------------------------------------------
## 1711.3 | 1716.1 | 1777.2 | 0.712 | 0.217 | 0.632 | 8.838 | 12.340Did these additional predictors improve the model?
comparison <- anova(m.rapp.select.ALL,m.rapp.select.ALL.agegender)
anova_tab <- as.data.frame(comparison)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "Comparison of rapport models with neurotypes: with and without age+gender predictors.",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | npar | AIC | BIC | logLik | -2*log(L) | Chisq | Df | Pr(>Chisq) |
|---|---|---|---|---|---|---|---|---|
| m.rapp.select.ALL | 17 | 1769.604 | 1825.590 | -867.802 | 1735.604 | NA | NA | NA |
| m.rapp.select.ALL.agegender | 20 | 1772.041 | 1837.907 | -866.020 | 1732.041 | 3.563 | 3 | 0.313 |
No. AIC, BIC, and X2 all point to the superiority of the simpler model.
Anova:
a.m.rapp.select.ALL.agegender = anova(m.rapp.select.ALL.agegender)
anova_tab <- as.data.frame(a.m.rapp.select.ALL.agegender)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the age+gender+language rapport model",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| fam_synchrony_expressivity | 1017.320 | 1017.320 | 1 | 105.518 | 6.680 | 0.011 |
| S_Q_PC1 | 2009.440 | 2009.440 | 1 | 58.660 | 13.195 | 0.001 |
| S_Q_PC3 | 811.100 | 811.100 | 1 | 89.986 | 5.326 | 0.023 |
| partner_neurotype | 168.042 | 168.042 | 1 | 41.480 | 1.103 | 0.300 |
| Group | 595.999 | 595.999 | 1 | 53.236 | 3.914 | 0.053 |
| age_diff | 218.176 | 218.176 | 1 | 107.702 | 1.433 | 0.234 |
| gender_match | 303.635 | 303.635 | 1 | 100.707 | 1.994 | 0.161 |
| language | 82.329 | 82.329 | 1 | 70.132 | 0.541 | 0.465 |
| fam_synchrony_expressivity:partner_neurotype | 0.267 | 0.267 | 1 | 153.526 | 0.002 | 0.967 |
| S_Q_PC1:partner_neurotype | 79.666 | 79.666 | 1 | 108.421 | 0.523 | 0.471 |
| S_Q_PC3:partner_neurotype | 384.035 | 384.035 | 1 | 97.022 | 2.522 | 0.116 |
| fam_synchrony_expressivity:Group | 70.523 | 70.523 | 1 | 161.254 | 0.463 | 0.497 |
| S_Q_PC1:Group | 361.204 | 361.204 | 1 | 58.301 | 2.372 | 0.129 |
| S_Q_PC3:Group | 324.775 | 324.775 | 1 | 94.285 | 2.133 | 0.148 |
| partner_neurotype:Group | 186.396 | 186.396 | 1 | 119.080 | 1.224 | 0.271 |
Language, age and gender match/difference don’t change the main results.
By DYAD TYPE
Rapport ~ predictors x dyad (LMM)
Model
Prepare the lmer formula:
main_predictors <- selected_predictor_summary$predictor
interaction_predictors <- main_predictors
fixed_formula_text <- paste0(
"R_Int_rapport ~ dyad + ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":dyad"), collapse = " + ")
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_rapport ~ dyad + fam_synchrony_expressivity + S_Q_PC1 +
## S_Q_PC3 + fam_synchrony_expressivity:dyad + S_Q_PC1:dyad +
## S_Q_PC3:dyad + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)Fit the model:
m.rapp.select.dyad <- lmer(
formula = mixed_formula,
data = DFrapp,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Performance
Singular fit?
print(isSingular(m.rapp.select.dyad, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.rapp.select.dyad), residuals(m.rapp.select.dyad), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.rapp.select.dyad), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.rapp.select.dyad), main='')
qqline(residuals(m.rapp.select.dyad))
layout(1)Variance components:
print(VarCorr(m.rapp.select.dyad))## Groups Name Std.Dev.
## dyad_indx (Intercept) 9.6893
## partner_ID (Intercept) 7.6873
## ID (Intercept) 11.1385
## Residual 12.8762Performance summary:
performance::model_performance(m.rapp.select.dyad)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ---------------------------------------------------------------------------
## 1728.9 | 1731.8 | 1781.5 | 0.683 | 0.154 | 0.626 | 9.291 | 12.876Interpretation
Anova:
a.m.rapp.select.dyad = anova(m.rapp.select.dyad)
anova_tab <- as.data.frame(a.m.rapp.select.dyad)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the rapport model with dyad type",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| dyad | 588.023 | 294.011 | 2 | 139.234 | 1.773 | 0.174 |
| fam_synchrony_expressivity | 820.952 | 820.952 | 1 | 105.140 | 4.952 | 0.028 |
| S_Q_PC1 | 723.744 | 723.744 | 1 | 102.869 | 4.365 | 0.039 |
| S_Q_PC3 | 1001.569 | 1001.569 | 1 | 100.273 | 6.041 | 0.016 |
| dyad:fam_synchrony_expressivity | 48.961 | 24.480 | 2 | 91.910 | 0.148 | 0.863 |
| dyad:S_Q_PC1 | 0.143 | 0.072 | 2 | 138.716 | 0.000 | 1.000 |
| dyad:S_Q_PC3 | 265.962 | 132.981 | 2 | 113.568 | 0.802 | 0.451 |
Plots:
p_combined_ALL_rapport_dyad <- plot_lmm_significant(
model = m.rapp.select.dyad,
anova_obj = a.m.rapp.select.dyad,
df = DFrapp,
grouping_var = "dyad",
n_cols_plotting = 3,
show_data = T,
colour_by="dyad" #dyad, Group, none
)
p_combined_ALL_rapport_dyad
See effect sizes:
get_anova_effect_sizes(m.rapp.select.dyad)## Parameter Eta2_partial CI CI_low CI_high Omega2_partial Omega2_CI_low Omega2_CI_high
## 1 dyad 0.02 0.95 0.00 1 0.01 0 1
## 2 fam_synchrony_expressivity 0.04 0.95 0.00 1 0.04 0 1
## 3 S_Q_PC1 0.04 0.95 0.00 1 0.03 0 1
## 4 S_Q_PC3 0.06 0.95 0.01 1 0.05 0 1
## 5 dyad:fam_synchrony_expressivity 0.00 0.95 0.00 1 0.00 0 1
## 6 dyad:S_Q_PC1 0.00 0.95 0.00 1 0.00 0 1
## 7 dyad:S_Q_PC3 0.01 0.95 0.00 1 0.00 0 1#get_effect_sizes(m.rapp.select.dyad)See BF:
INTER: expressivity
model_null <- update(m.rapp.select.dyad, . ~ . - fam_synchrony_expressivity - fam_synchrony_expressivity:dyad)
compare_models_BF(model_null = model_null,
model_full = m.rapp.select.dyad,
model_null_name = "Without INTER: expressivity",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without INTER: expressivity model BIC: 1820.45
## - Full model model BIC: 1830.06
## - BIC difference (Full - Null): 9.61
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without INTER: expressivity vs Full model ): 122
## - BF₁₀ (Evidence for Full model vs Without INTER: expressivity ): 0.0082
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without INTER: expressivity
## - Evidence strength: Decisive/Extreme
## - Bayes Factor: 122.17 :1 in favor of Without INTER: expressivity
##
## CONCLUSION: There is decisive/extreme evidence ( 122.17 :1) in favor
## of the Without INTER: expressivity model over the Full model model.
## =======================================================================Q1:
model_null <- update(m.rapp.select.dyad, . ~ . - S_Q_PC1 - S_Q_PC1:dyad)
compare_models_BF(model_null = model_null,
model_full = m.rapp.select.dyad,
model_null_name = "Without autism-ike traits (PC1)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without autism-ike traits (PC1) model BIC: 1820.28
## - Full model model BIC: 1830.06
## - BIC difference (Full - Null): 9.79
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without autism-ike traits (PC1) vs Full model ): 133
## - BF₁₀ (Evidence for Full model vs Without autism-ike traits (PC1) ): 0.0075
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without autism-ike traits (PC1)
## - Evidence strength: Decisive/Extreme
## - Bayes Factor: 133.32 :1 in favor of Without autism-ike traits (PC1)
##
## CONCLUSION: There is decisive/extreme evidence ( 133.32 :1) in favor
## of the Without autism-ike traits (PC1) model over the Full model model.
## =======================================================================Q3:
model_null <- update(m.rapp.select.dyad, . ~ . - S_Q_PC3 - S_Q_PC3:dyad)
compare_models_BF(model_null = model_null,
model_full = m.rapp.select.dyad,
model_null_name = "Without cognitive engagement (PC3)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without cognitive engagement (PC3) model BIC: 1823.57
## - Full model model BIC: 1830.06
## - BIC difference (Full - Null): 6.49
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without cognitive engagement (PC3) vs Full model ): 26
## - BF₁₀ (Evidence for Full model vs Without cognitive engagement (PC3) ): 0.039
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without cognitive engagement (PC3)
## - Evidence strength: Strong
## - Bayes Factor: 25.72 :1 in favor of Without cognitive engagement (PC3)
##
## CONCLUSION: There is strong evidence ( 25.72 :1) in favor
## of the Without cognitive engagement (PC3) model over the Full model model.
## =======================================================================With gender + age + language
main_predictors <- selected_predictor_summary$predictor
interaction_predictors <- main_predictors
fixed_formula_text <- paste0(
"R_Int_rapport ~ dyad + ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":dyad"), collapse = " + ")
#,"+",
# paste(paste0(interaction_predictors, ":age_diff"), collapse = " + "),"+", # uncomment if with interactions
# paste(paste0(interaction_predictors, ":gender_match"), collapse = " + "),"+",
# paste(paste0(interaction_predictors, ":language"), collapse = " + ")
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + age_diff + gender_match + language + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_rapport ~ dyad + fam_synchrony_expressivity + S_Q_PC1 +
## S_Q_PC3 + fam_synchrony_expressivity:dyad + S_Q_PC1:dyad +
## S_Q_PC3:dyad + age_diff + gender_match + language + (1 |
## ID) + (1 | partner_ID) + (1 | dyad_indx)Fit the model:
m.rapp.select.dyad.agegender <- lmer(
formula = mixed_formula,
data = DFrapp,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Singular fit?
print(isSingular(m.rapp.select.dyad.agegender, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.rapp.select.dyad.agegender), residuals(m.rapp.select.dyad.agegender), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.rapp.select.dyad.agegender), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.rapp.select.dyad.agegender), main='')
qqline(residuals(m.rapp.select.dyad.agegender))
layout(1)Variance components:
print(VarCorr(m.rapp.select.dyad.agegender))## Groups Name Std.Dev.
## dyad_indx (Intercept) 9.4709
## partner_ID (Intercept) 7.4647
## ID (Intercept) 10.7951
## Residual 13.0081Performance summary:
performance::model_performance(m.rapp.select.dyad.agegender)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ---------------------------------------------------------------------------
## 1722.7 | 1727.0 | 1785.3 | 0.675 | 0.171 | 0.608 | 9.490 | 13.008Did these adidtional predictors improve the fit?
comparison <- anova(m.rapp.select.dyad,m.rapp.select.dyad.agegender)
anova_tab <- as.data.frame(comparison)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "Comparison of rapport models with dyad type: with and without age+gender+language predictors.",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | npar | AIC | BIC | logLik | -2*log(L) | Chisq | Df | Pr(>Chisq) |
|---|---|---|---|---|---|---|---|---|
| m.rapp.select.dyad | 16 | 1777.368 | 1830.061 | -872.684 | 1745.368 | NA | NA | NA |
| m.rapp.select.dyad.agegender | 19 | 1780.360 | 1842.932 | -871.180 | 1742.360 | 3.009 | 3 | 0.39 |
No.
Anova:
a.m.rapp.select.dyad.agegender = anova(m.rapp.select.dyad.agegender)
anova_tab <- as.data.frame(a.m.rapp.select.dyad.agegender)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the age+gender+language rapport model with dyad type",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| dyad | 513.744 | 256.872 | 2 | 141.825 | 1.518 | 0.223 |
| fam_synchrony_expressivity | 903.813 | 903.813 | 1 | 102.865 | 5.341 | 0.023 |
| S_Q_PC1 | 690.209 | 690.209 | 1 | 104.696 | 4.079 | 0.046 |
| S_Q_PC3 | 975.730 | 975.730 | 1 | 100.252 | 5.766 | 0.018 |
| age_diff | 122.854 | 122.854 | 1 | 102.569 | 0.726 | 0.396 |
| gender_match | 277.108 | 277.108 | 1 | 99.778 | 1.638 | 0.204 |
| language | 167.497 | 167.497 | 1 | 62.977 | 0.990 | 0.324 |
| dyad:fam_synchrony_expressivity | 56.156 | 28.078 | 2 | 92.170 | 0.166 | 0.847 |
| dyad:S_Q_PC1 | 1.322 | 0.661 | 2 | 138.752 | 0.004 | 0.996 |
| dyad:S_Q_PC3 | 213.895 | 106.948 | 2 | 112.930 | 0.632 | 0.533 |
ANALYSIS 2: predictors of SOCIAL MOTIVATION
Predictors screening (Ridge)
Step 1: Build the analysis dataset
analysis_predictors <- c(family_score_vars,S_trait_pc_vars)
DFwill = DF_reduced_predictors[,c(1:7,which(names(DF_reduced_predictors) %in% analysis_predictors))]Step 2: Ridge as a predictor screening procedure
predictor_main_vars <- analysis_predictors
ridge_formula <- as.formula(
paste(
"R_Int_willingness ~ (",
paste(predictor_main_vars, collapse = " + "),
") * Group"
)
)
x_ridge <- model.matrix(ridge_formula, data = DFwill)[, -1] # drop intercept
y_ridge <- DFwill$R_Int_willingness
cat("Design matrix dimensions:", dim(x_ridge), "\n")## Design matrix dimensions: 199 17cat("Number of columns in design matrix:", ncol(x_ridge), "\n\n")## Number of columns in design matrix: 17First 20 column names:
print(colnames(x_ridge)[1:min(20, ncol(x_ridge))])## [1] "fam_intonation" "fam_turn_taking" "fam_backchannel_tpt"
## [4] "fam_expressivity" "fam_synchrony_expressivity" "S_Q_PC1"
## [7] "S_Q_PC2" "S_Q_PC3" "GroupNAUT"
## [10] "fam_intonation:GroupNAUT" "fam_turn_taking:GroupNAUT" "fam_backchannel_tpt:GroupNAUT"
## [13] "fam_expressivity:GroupNAUT" "fam_synchrony_expressivity:GroupNAUT" "S_Q_PC1:GroupNAUT"
## [16] "S_Q_PC2:GroupNAUT" "S_Q_PC3:GroupNAUT"Columns containing dyad interactions:
print(grep(":", colnames(x_ridge), value = TRUE))## [1] "fam_intonation:GroupNAUT" "fam_turn_taking:GroupNAUT" "fam_backchannel_tpt:GroupNAUT"
## [4] "fam_expressivity:GroupNAUT" "fam_synchrony_expressivity:GroupNAUT" "S_Q_PC1:GroupNAUT"
## [7] "S_Q_PC2:GroupNAUT" "S_Q_PC3:GroupNAUT"Step 3: Grouped folds by ID
Fold sizes:
id_ridge <- as.character(DFwill$ID)
unique_ids_ridge <- unique(id_ridge)
nfolds <- 10
foldid_by_id <- sample(rep(1:nfolds, length.out = length(unique_ids_ridge)))
names(foldid_by_id) <- unique_ids_ridge
foldid_ridge <- unname(foldid_by_id[id_ridge])
print(table(foldid_ridge))## foldid_ridge
## 1 2 3 4 5 6 7 8 9 10
## 18 20 23 22 20 22 20 19 18 17Step 4: Grouped cross validation (CV) Ridge
cv_ridge <- glmnet::cv.glmnet(
x = x_ridge,
y = y_ridge,
alpha = 0,
family = "gaussian",
standardize = TRUE,
foldid = foldid_ridge
)
coef_min <- as.matrix(coef(cv_ridge, s = "lambda.min"))
coef_1se <- as.matrix(coef(cv_ridge, s = "lambda.1se"))
coef_min_df <- data.frame(
variable = rownames(coef_min),
coef = as.numeric(coef_min[, 1])
)
coef_1se_df <- data.frame(
variable = rownames(coef_1se),
coef = as.numeric(coef_1se[, 1])
)
coef_min_df <- coef_min_df[order(abs(coef_min_df$coef), decreasing = TRUE), ]
coef_1se_df <- coef_1se_df[order(abs(coef_1se_df$coef), decreasing = TRUE), ]Top coefficients at lambda.min:
print(head(coef_min_df, 20))## variable coef
## 1 (Intercept) 58.24120603015075658959176507778465747833251953
## 8 S_Q_PC2 0.00000000000000000000000000000000000525435935
## 10 GroupNAUT -0.00000000000000000000000000000000000468778696
## 17 S_Q_PC2:GroupNAUT 0.00000000000000000000000000000000000446936080
## 9 S_Q_PC3 0.00000000000000000000000000000000000344067168
## 16 S_Q_PC1:GroupNAUT -0.00000000000000000000000000000000000245052183
## 7 S_Q_PC1 -0.00000000000000000000000000000000000229017417
## 3 fam_turn_taking 0.00000000000000000000000000000000000177397563
## 13 fam_backchannel_tpt:GroupNAUT 0.00000000000000000000000000000000000153735090
## 4 fam_backchannel_tpt 0.00000000000000000000000000000000000134674367
## 18 S_Q_PC3:GroupNAUT 0.00000000000000000000000000000000000122118325
## 15 fam_synchrony_expressivity:GroupNAUT -0.00000000000000000000000000000000000120594273
## 5 fam_expressivity 0.00000000000000000000000000000000000120111785
## 14 fam_expressivity:GroupNAUT 0.00000000000000000000000000000000000105621138
## 12 fam_turn_taking:GroupNAUT 0.00000000000000000000000000000000000104750386
## 2 fam_intonation -0.00000000000000000000000000000000000094226786
## 6 fam_synchrony_expressivity 0.00000000000000000000000000000000000015529310
## 11 fam_intonation:GroupNAUT 0.00000000000000000000000000000000000006121268Top coefficients at lambda.1se:
print(head(coef_1se_df, 20))## variable coef
## 1 (Intercept) 58.24120603015075658959176507778465747833251953
## 8 S_Q_PC2 0.00000000000000000000000000000000000525435935
## 10 GroupNAUT -0.00000000000000000000000000000000000468778696
## 17 S_Q_PC2:GroupNAUT 0.00000000000000000000000000000000000446936080
## 9 S_Q_PC3 0.00000000000000000000000000000000000344067168
## 16 S_Q_PC1:GroupNAUT -0.00000000000000000000000000000000000245052183
## 7 S_Q_PC1 -0.00000000000000000000000000000000000229017417
## 3 fam_turn_taking 0.00000000000000000000000000000000000177397563
## 13 fam_backchannel_tpt:GroupNAUT 0.00000000000000000000000000000000000153735090
## 4 fam_backchannel_tpt 0.00000000000000000000000000000000000134674367
## 18 S_Q_PC3:GroupNAUT 0.00000000000000000000000000000000000122118325
## 15 fam_synchrony_expressivity:GroupNAUT -0.00000000000000000000000000000000000120594273
## 5 fam_expressivity 0.00000000000000000000000000000000000120111785
## 14 fam_expressivity:GroupNAUT 0.00000000000000000000000000000000000105621138
## 12 fam_turn_taking:GroupNAUT 0.00000000000000000000000000000000000104750386
## 2 fam_intonation -0.00000000000000000000000000000000000094226786
## 6 fam_synchrony_expressivity 0.00000000000000000000000000000000000015529310
## 11 fam_intonation:GroupNAUT 0.00000000000000000000000000000000000006121268The top results include interactions, which confirms that the willingness predictors will likely be different between neurotypes. But this is just one ridge fit. Let’s check the stability with bootstrapping.
Step 5: Bootstrapping
Grouped bootstrap + grouped CV inside bootstrap for the full ridge screening model.
B <- 500
nfolds <- 10
unique_ids_boot <- unique(as.character(DFwill$ID))
boot_coef_mat <- matrix(NA, nrow = B, ncol = ncol(x_ridge))
colnames(boot_coef_mat) <- colnames(x_ridge)
boot_info <- data.frame(
iter = 1:B,
n_rows = NA_integer_,
n_id_boot = NA_integer_,
lambda_min = NA_real_,
lambda_1se = NA_real_,
cv_ok = FALSE,
error_msg = NA_character_
)
for (b in 1:B) {
boot_ids <- sample(unique_ids_boot, size = length(unique_ids_boot), replace = TRUE)
boot_list <- lapply(seq_along(boot_ids), function(i) {
id_now <- boot_ids[i]
tmp <- DFwill[as.character(DFwill$ID) == id_now, , drop = FALSE]
tmp$ID_boot <- paste0(id_now, "_rep", i)
tmp
})
boot_dat <- do.call(rbind, boot_list)
boot_info$n_rows[b] <- nrow(boot_dat)
boot_info$n_id_boot[b] <- length(unique(boot_dat$ID_boot))
x_boot <- model.matrix(ridge_formula, data = boot_dat)[, -1, drop = FALSE]
y_boot <- boot_dat$R_Int_willingness
id_boot <- as.character(boot_dat$ID_boot)
unique_id_boot <- unique(id_boot)
foldid_by_id_boot <- sample(rep(1:nfolds, length.out = length(unique_id_boot)))
names(foldid_by_id_boot) <- unique_id_boot
foldid_boot <- unname(foldid_by_id_boot[id_boot])
cv_boot <- tryCatch(
glmnet::cv.glmnet(
x = x_boot,
y = y_boot,
alpha = 0,
family = "gaussian",
standardize = TRUE,
foldid = foldid_boot
),
error = function(e) e
)
if (inherits(cv_boot, "error")) {
boot_info$error_msg[b] <- conditionMessage(cv_boot)
next
}
boot_info$cv_ok[b] <- TRUE
boot_info$lambda_min[b] <- cv_boot$lambda.min
boot_info$lambda_1se[b] <- cv_boot$lambda.1se
coef_boot <- as.matrix(coef(cv_boot, s = "lambda.min"))[-1, 1]
boot_coef_mat[b, names(coef_boot)] <- coef_boot
}Summaries
boot_summary_full <- data.frame(
variable = colnames(boot_coef_mat),
mean_coef = colMeans(boot_coef_mat, na.rm = TRUE),
sd_coef = apply(boot_coef_mat, 2, sd, na.rm = TRUE),
prop_positive = colMeans(boot_coef_mat > 0, na.rm = TRUE),
prop_negative = colMeans(boot_coef_mat < 0, na.rm = TRUE),
mean_abs_coef = colMeans(abs(boot_coef_mat), na.rm = TRUE)
)
boot_summary_full <- boot_summary_full[order(boot_summary_full$mean_abs_coef, decreasing = TRUE), ]
cat("Successful bootstrap CV runs:", sum(boot_info$cv_ok), "out of", B, "\n")## Successful bootstrap CV runs: 500 out of 500cat("Failed runs:", sum(!boot_info$cv_ok), "\n\n")## Failed runs: 0if (sum(!boot_info$cv_ok) > 0) {
cat("Bootstrap errors:\n")
print(sort(table(na.omit(boot_info$error_msg)), decreasing = TRUE))
}Top 25 bootstrap-stable coefficients:
print(head(boot_summary_full, 25))## variable mean_coef sd_coef prop_positive prop_negative
## GroupNAUT GroupNAUT -3.25067254 4.8712603 0.124 0.876
## S_Q_PC3:GroupNAUT S_Q_PC3:GroupNAUT 0.48739705 3.2807362 0.610 0.390
## S_Q_PC2 S_Q_PC2 1.82024026 1.8154380 0.978 0.022
## S_Q_PC3 S_Q_PC3 1.77340616 1.9445139 0.918 0.082
## S_Q_PC1 S_Q_PC1 -1.75600615 2.6367201 0.052 0.948
## S_Q_PC1:GroupNAUT S_Q_PC1:GroupNAUT -0.29465981 2.5267526 0.294 0.706
## fam_expressivity:GroupNAUT fam_expressivity:GroupNAUT -0.38490001 2.4322455 0.472 0.528
## S_Q_PC2:GroupNAUT S_Q_PC2:GroupNAUT 0.72120767 1.9238513 0.804 0.196
## fam_intonation:GroupNAUT fam_intonation:GroupNAUT 0.04528852 1.7542778 0.542 0.458
## fam_expressivity fam_expressivity 0.54462122 1.5874055 0.668 0.332
## fam_backchannel_tpt:GroupNAUT fam_backchannel_tpt:GroupNAUT 0.44165917 1.4681251 0.680 0.320
## fam_turn_taking:GroupNAUT fam_turn_taking:GroupNAUT 0.18789007 1.4957385 0.552 0.448
## fam_synchrony_expressivity:GroupNAUT fam_synchrony_expressivity:GroupNAUT -0.12318604 1.5341493 0.386 0.614
## fam_turn_taking fam_turn_taking 0.69639535 1.0729720 0.828 0.172
## fam_intonation fam_intonation -0.63432587 1.0619472 0.192 0.808
## fam_synchrony_expressivity fam_synchrony_expressivity 0.44765142 0.8529433 0.736 0.264
## fam_backchannel_tpt fam_backchannel_tpt 0.15715121 0.8729096 0.684 0.316
## mean_abs_coef
## GroupNAUT 3.4928008
## S_Q_PC3:GroupNAUT 2.0383883
## S_Q_PC2 1.8371861
## S_Q_PC3 1.8133134
## S_Q_PC1 1.7719956
## S_Q_PC1:GroupNAUT 1.4288803
## fam_expressivity:GroupNAUT 1.4076847
## S_Q_PC2:GroupNAUT 1.3420491
## fam_intonation:GroupNAUT 1.0448970
## fam_expressivity 1.0158789
## fam_backchannel_tpt:GroupNAUT 0.9397192
## fam_turn_taking:GroupNAUT 0.9105197
## fam_synchrony_expressivity:GroupNAUT 0.8838132
## fam_turn_taking 0.8346911
## fam_intonation 0.7997446
## fam_synchrony_expressivity 0.5957906
## fam_backchannel_tpt 0.5731228The ranking gives us now a basis for which predictors to use in the final mixed models. Again, the strongest and most stable signals are mostly dyad-specific interaction terms, not just pooled main effects: so we can’t screen on the entire sample.
Step 6: Informed feature selection
To select the predictors that will enter the mixed models’ analysis, we will choose those that: - have relatively high mean_abs_coef (we choose 2) - have a strong sign stability (at least 0.8 in any direction) - all main effects of any interaction that fits the rules above
# Cutoff for the stability in either direction
sign_stability_cutoff <- 0.8
# Cutoff for the mean coefficient score
mean_abs_cutoff <- 1.5
boot_screen <- boot_summary_full
boot_screen$sign_stability <- pmax(boot_screen$prop_positive, boot_screen$prop_negative)
boot_screen$is_interaction <- grepl(":", boot_screen$variable)
# keep strong and sign-stable terms
selected_core <- boot_screen[
boot_screen$mean_abs_coef >= mean_abs_cutoff &
boot_screen$sign_stability >= sign_stability_cutoff,
]
# extract main effects corresponding to selected interactions
selected_interactions <- selected_core$variable[selected_core$is_interaction]
interaction_bases <- unique(sub(":.*$", "", selected_interactions))
# keep main effects for any selected interaction, if present in boot table
main_effects_to_add <- boot_screen$variable[
boot_screen$variable %in% interaction_bases
]
# final selected list = core selected + hierarchy-preserving main effects
selected_final_vars <- unique(c(selected_core$variable, main_effects_to_add))
selected_final <- boot_screen[boot_screen$variable %in% selected_final_vars, ]
selected_final <- selected_final[order(selected_final$is_interaction, -selected_final$mean_abs_coef), ]
cat("mean_abs_cutoff =", mean_abs_cutoff, "\n")## mean_abs_cutoff = 1.5cat("sign_stability_cutoff =", sign_stability_cutoff, "\n\n")## sign_stability_cutoff = 0.8Selected core terms:
print(
selected_core[order(selected_core$mean_abs_coef, decreasing = TRUE),
c("variable", "mean_coef", "mean_abs_coef", "prop_positive", "prop_negative", "sign_stability")]
)## variable mean_coef mean_abs_coef prop_positive prop_negative sign_stability
## GroupNAUT GroupNAUT -3.250673 3.492801 0.124 0.876 0.876
## S_Q_PC2 S_Q_PC2 1.820240 1.837186 0.978 0.022 0.978
## S_Q_PC3 S_Q_PC3 1.773406 1.813313 0.918 0.082 0.918
## S_Q_PC1 S_Q_PC1 -1.756006 1.771996 0.052 0.948 0.948Main effects added to preserve hierarchy:
print(setdiff(main_effects_to_add, selected_core$variable))## character(0)Final selected terms:
print(
selected_final[, c("variable", "mean_coef", "mean_abs_coef", "prop_positive", "prop_negative", "sign_stability", "is_interaction")]
)## variable mean_coef mean_abs_coef prop_positive prop_negative sign_stability is_interaction
## GroupNAUT GroupNAUT -3.250673 3.492801 0.124 0.876 0.876 FALSE
## S_Q_PC2 S_Q_PC2 1.820240 1.837186 0.978 0.022 0.978 FALSE
## S_Q_PC3 S_Q_PC3 1.773406 1.813313 0.918 0.082 0.918 FALSE
## S_Q_PC1 S_Q_PC1 -1.756006 1.771996 0.052 0.948 0.948 FALSE#Collapse the ridge terms to variable-level output for the mixed models:
selected_terms <- selected_final$variable
collapse_to_predictor <- function(x) {
# interaction term -> keep only variable before :
if (grepl(":", x)) {
return(sub(":.*$", "", x))
}
# dyad dummy terms -> collapse to dyad
if (grepl("^Group", x)) {
return("Group")
}
# all other main effects stay as they are
return(x)
}
collapsed_predictors <- unique(vapply(selected_terms, collapse_to_predictor, character(1)))
# which predictors were selected as interactions with dyad?
selected_interaction_predictors <- unique(
vapply(
selected_terms[grepl(":", selected_terms)],
collapse_to_predictor,
character(1)
)
)
# final variable-level summary table
selected_predictor_summary <- data.frame(
predictor = collapsed_predictors,
include_main_effect = TRUE,
include_group_interaction = collapsed_predictors %in% selected_interaction_predictors
)
# sort a bit more nicely
selected_predictor_summary$type <- ifelse(
selected_predictor_summary$predictor == "Group", "grouping_variable",
ifelse(grepl("^fam_", selected_predictor_summary$predictor), "behavior_family", "trait_PC")
)
selected_predictor_summary <- selected_predictor_summary[
order(selected_predictor_summary$type, selected_predictor_summary$predictor),
]Final predictor list for mixed models:
print(selected_predictor_summary)## predictor include_main_effect include_group_interaction type
## 1 Group TRUE FALSE grouping_variable
## 4 S_Q_PC1 TRUE FALSE trait_PC
## 2 S_Q_PC2 TRUE FALSE trait_PC
## 3 S_Q_PC3 TRUE FALSE trait_PCBehavior families selected:
print(selected_predictor_summary$predictor[selected_predictor_summary$type == "behavior_family"])## character(0)Trait PCs selected:
print(selected_predictor_summary$predictor[selected_predictor_summary$type == "trait_PC"])## [1] "S_Q_PC1" "S_Q_PC2" "S_Q_PC3"Predictors that should be tested with Group interactions:
print(selected_predictor_summary$predictor[selected_predictor_summary$include_group_interaction])## character(0)# For plotting later:
boot_summary_motivation <- boot_summary_fullBy NEUROTYPE
Rapport ~ predictors x S_neurotype x P_neurotype (LMM)
Model
Prepare the lmer formula:
DFwill$partner_neurotype <- ifelse(
DFwill$dyad == "MIX" & DFwill$Group == "NAUT", "AUT",
ifelse(
DFwill$dyad == "MIX" & DFwill$Group == "AUT", "NAUT",
as.character(DFwill$Group)
)
)
DFwill$partner_neurotype <- factor(
DFwill$partner_neurotype,
levels = c("AUT", "NAUT")
)
main_predictors <- selected_predictor_summary$predictor
interaction_predictors <- main_predictors
main_predictors <- c(main_predictors,"partner_neurotype","Group")
fixed_formula_text <- paste0(
"R_Int_willingness ~ ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":partner_neurotype"), collapse = " + "),
"+",
paste(paste0(interaction_predictors, ":Group"), collapse = " + "),
"+ partner_neurotype:Group"
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_willingness ~ Group + S_Q_PC1 + S_Q_PC2 + S_Q_PC3 + partner_neurotype +
## Group + Group:partner_neurotype + S_Q_PC1:partner_neurotype +
## S_Q_PC2:partner_neurotype + S_Q_PC3:partner_neurotype + Group:Group +
## S_Q_PC1:Group + S_Q_PC2:Group + S_Q_PC3:Group + partner_neurotype:Group +
## (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)Fit the model:
m.will.select.ALL <- lmer(
formula = mixed_formula,
data = DFwill,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Performance
Singular fit?
print(isSingular(m.will.select.ALL, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.will.select.ALL), residuals(m.will.select.ALL), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.will.select.ALL), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.will.select.ALL), main='')
qqline(residuals(m.will.select.ALL))
layout(1)Variance components:
print(VarCorr(m.will.select.ALL))## Groups Name Std.Dev.
## dyad_indx (Intercept) 11.1271
## partner_ID (Intercept) 9.6871
## ID (Intercept) 13.9927
## Residual 13.5796Performance summary:
performance::model_performance(m.will.select.ALL)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ---------------------------------------------------------------------------
## 1770.6 | 1774.0 | 1826.6 | 0.745 | 0.175 | 0.692 | 9.427 | 13.580Interpretation
Anova:
a.m.will.select.ALL = anova(m.will.select.ALL)
anova_tab <- as.data.frame(a.m.will.select.ALL)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the social motivation model with neurotypes",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| Group | 812.654 | 812.654 | 1 | 55.511 | 4.407 | 0.040 |
| S_Q_PC1 | 1151.152 | 1151.152 | 1 | 59.586 | 6.242 | 0.015 |
| S_Q_PC2 | 341.211 | 341.211 | 1 | 65.937 | 1.850 | 0.178 |
| S_Q_PC3 | 572.501 | 572.501 | 1 | 94.513 | 3.105 | 0.081 |
| partner_neurotype | 26.128 | 26.128 | 1 | 43.753 | 0.142 | 0.708 |
| Group:partner_neurotype | 236.136 | 236.136 | 1 | 105.886 | 1.281 | 0.260 |
| S_Q_PC1:partner_neurotype | 95.259 | 95.259 | 1 | 117.392 | 0.517 | 0.474 |
| S_Q_PC2:partner_neurotype | 1753.754 | 1753.754 | 1 | 132.711 | 9.510 | 0.002 |
| S_Q_PC3:partner_neurotype | 287.789 | 287.789 | 1 | 100.105 | 1.561 | 0.214 |
| Group:S_Q_PC1 | 371.034 | 371.034 | 1 | 62.130 | 2.012 | 0.161 |
| Group:S_Q_PC2 | 1.644 | 1.644 | 1 | 67.700 | 0.009 | 0.925 |
| Group:S_Q_PC3 | 130.530 | 130.530 | 1 | 96.466 | 0.708 | 0.402 |
Plots:
p_combined_ALL_motivation_neurotype <- plot_lmm_significant(
model = m.will.select.ALL,
anova_obj = a.m.will.select.ALL,
df = DFwill,
grouping_var = "partner_neurotype",
DV = "R_Int_willingness",
n_cols_plotting = 3,
show_data = T,
colour_by = "Group"
)
p_combined_ALL_motivation_neurotype
Check the partner’s neurotype comparisons for autism attitudes
Against 0:
# Estimated S_Q_PC2 slope within each dyad type
tr <- emtrends(
m.will.select.ALL,
~ partner_neurotype,
var = "S_Q_PC2"
)Against 0:
# Test whether each dyad-specific slope differs from zero
summary(tr, infer = c(TRUE, TRUE), adjust = "holm")## partner_neurotype S_Q_PC2.trend SE df lower.CL upper.CL t.ratio p.value
## AUT 8.145 3.32 98.5 0.584 15.71 2.452 0.0319
## NAUT -0.336 3.49 110.6 -8.265 7.59 -0.096 0.9235
##
## Results are averaged over the levels of: Group
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: holm method for 2 testsSee effect sizes:
get_anova_effect_sizes(m.will.select.ALL)## Parameter Eta2_partial CI CI_low CI_high Omega2_partial Omega2_CI_low Omega2_CI_high
## 1 Group 0.07 0.95 0.00 1 0.06 0.00 1
## 2 S_Q_PC1 0.09 0.95 0.01 1 0.08 0.00 1
## 3 S_Q_PC2 0.03 0.95 0.00 1 0.01 0.00 1
## 4 S_Q_PC3 0.03 0.95 0.00 1 0.02 0.00 1
## 5 partner_neurotype 0.00 0.95 0.00 1 0.00 0.00 1
## 6 Group:partner_neurotype 0.01 0.95 0.00 1 0.00 0.00 1
## 7 S_Q_PC1:partner_neurotype 0.00 0.95 0.00 1 0.00 0.00 1
## 8 S_Q_PC2:partner_neurotype 0.07 0.95 0.01 1 0.06 0.01 1
## 9 S_Q_PC3:partner_neurotype 0.02 0.95 0.00 1 0.01 0.00 1
## 10 Group:S_Q_PC1 0.03 0.95 0.00 1 0.02 0.00 1
## 11 Group:S_Q_PC2 0.00 0.95 0.00 1 0.00 0.00 1
## 12 Group:S_Q_PC3 0.01 0.95 0.00 1 0.00 0.00 1#get_effect_sizes(m.will.select.ALL)See BF:
Rater neurotype:
model_null <- update(m.will.select.ALL, . ~ . - Group -
Group:partner_neurotype -
S_Q_PC1:Group -
S_Q_PC2:Group -
S_Q_PC3:Group)
compare_models_BF(model_null = model_null,
model_full = m.will.select.ALL,
model_null_name = "Without rater neurotype",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without rater neurotype model BIC: 1865.02
## - Full model model BIC: 1884.07
## - BIC difference (Full - Null): 19.05
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without rater neurotype vs Full model ): 13681
## - BF₁₀ (Evidence for Full model vs Without rater neurotype ): 0.000073
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without rater neurotype
## - Evidence strength: Decisive/Extreme
## - Bayes Factor: 13680.85 :1 in favor of Without rater neurotype
##
## CONCLUSION: There is decisive/extreme evidence ( 13680.85 :1) in favor
## of the Without rater neurotype model over the Full model model.
## =======================================================================Q1:
model_null <- update(m.will.select.ALL, . ~ . - S_Q_PC1 -
S_Q_PC1:partner_neurotype -
S_Q_PC1:Group)
compare_models_BF(model_null = model_null,
model_full = m.will.select.ALL,
model_null_name = "Without autism-ike traits (PC1)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without autism-ike traits (PC1) model BIC: 1875.11
## - Full model model BIC: 1884.07
## - BIC difference (Full - Null): 8.96
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without autism-ike traits (PC1) vs Full model ): 88
## - BF₁₀ (Evidence for Full model vs Without autism-ike traits (PC1) ): 0.011
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without autism-ike traits (PC1)
## - Evidence strength: Very Strong
## - Bayes Factor: 88.13 :1 in favor of Without autism-ike traits (PC1)
##
## CONCLUSION: There is very strong evidence ( 88.13 :1) in favor
## of the Without autism-ike traits (PC1) model over the Full model model.
## =======================================================================Q2:
model_null <- update(m.will.select.ALL, . ~ . - S_Q_PC2 -
S_Q_PC2:partner_neurotype -
S_Q_PC2:Group)
compare_models_BF(model_null = model_null,
model_full = m.will.select.ALL,
model_null_name = "Without autism-ike traits (PC2)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without autism-ike traits (PC2) model BIC: 1880.58
## - Full model model BIC: 1884.07
## - BIC difference (Full - Null): 3.49
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without autism-ike traits (PC2) vs Full model ): 5.7
## - BF₁₀ (Evidence for Full model vs Without autism-ike traits (PC2) ): 0.17
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without autism-ike traits (PC2)
## - Evidence strength: Moderate
## - Bayes Factor: 5.73 :1 in favor of Without autism-ike traits (PC2)
##
## CONCLUSION: There is moderate evidence ( 5.73 :1) in favor
## of the Without autism-ike traits (PC2) model over the Full model model.
## =======================================================================With gender + age
We add age difference (difference between the two partners), gender matching (0 vs 1), and language.
# Add age difference and gender matching
# 1) Collapse DF to one row per person per dyad
person_info <- DFmain %>%
dplyr::select(dyad_indx, ID, age, gender, language) %>%
dplyr::group_by(dyad_indx, ID, language) %>%
dplyr::summarise(
age_self = mean(age, na.rm = TRUE),
gender_self = dplyr::first(gender),
.groups = "drop"
)
# 2) Self-join within dyad to get partner’s info
dyad_pairs <- person_info %>%
dplyr::inner_join(
person_info,
by = c("dyad_indx","language"),
suffix = c("_self", "_partner"),
relationship = "many-to-many"
) %>%
dplyr::filter(ID_self != ID_partner)
# 3) Clean to one row per person per dyad
dyad_info_unique <- dyad_pairs %>%
dplyr::transmute(
dyad_indx,
ID = ID_self,
age_self = age_self_self,
gender_self = gender_self_self,
age_partner = age_self_partner,
gender_partner = gender_self_partner,
language = language
)
# 4) Join into DFwill
DFwill <- DFwill %>%
dplyr::select(
-dplyr::any_of(c(
"age_self", "gender_self",
"age_partner", "gender_partner",
"age_diff", "gender_match"
))
) %>%
dplyr::left_join(
dyad_info_unique,
by = c("dyad_indx", "ID")
) %>%
dplyr::mutate(
age_diff = age_self - age_partner,
gender_match = dplyr::if_else(
!is.na(gender_self) &
!is.na(gender_partner) &
gender_self == gender_partner,
1L,
0L
)
)
# # Quick sanity check
# DFwill %>%
# dplyr::select(
# dyad_indx, ID, partner_ID,
# age_self, age_partner, age_diff,
# gender_self, gender_partner, gender_match
# ) %>%
# head()
DFwill$age_diff <- abs(DFwill$age_diff)Fit the model:
fixed_formula_text <- paste0(
"R_Int_willingness ~ ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":partner_neurotype"), collapse = " + "),"+",
paste(paste0(interaction_predictors, ":Group"), collapse = " + "),"+",
#
# paste(paste0(interaction_predictors, ":age_diff"), collapse = " + "),"+", # uncomment if with interactions
# paste(paste0(interaction_predictors, ":gender_match"), collapse = " + "),"+",
# paste(paste0(interaction_predictors, ":language"), collapse = " + "),
#
"+ partner_neurotype:Group"
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + age_diff + gender_match + language + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_willingness ~ Group + S_Q_PC1 + S_Q_PC2 + S_Q_PC3 + partner_neurotype +
## Group + Group:partner_neurotype + S_Q_PC1:partner_neurotype +
## S_Q_PC2:partner_neurotype + S_Q_PC3:partner_neurotype + Group:Group +
## S_Q_PC1:Group + S_Q_PC2:Group + S_Q_PC3:Group + +partner_neurotype:Group +
## age_diff + gender_match + language + (1 | ID) + (1 | partner_ID) +
## (1 | dyad_indx)m.will.select.ALL.agegender <- lmer(
formula = mixed_formula,
data = DFwill,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Singular fit?
print(isSingular(m.will.select.ALL.agegender, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.will.select.ALL.agegender), residuals(m.will.select.ALL.agegender), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.will.select.ALL.agegender), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.will.select.ALL.agegender), main='')
qqline(residuals(m.will.select.ALL.agegender))
layout(1)Variance components:
print(VarCorr(m.will.select.ALL.agegender))## Groups Name Std.Dev.
## dyad_indx (Intercept) 11.0995
## partner_ID (Intercept) 9.3291
## ID (Intercept) 13.7316
## Residual 13.6500Performance summary:
performance::model_performance(m.will.select.ALL.agegender)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ---------------------------------------------------------------------------
## 1764.1 | 1768.8 | 1829.9 | 0.743 | 0.191 | 0.682 | 9.532 | 13.650Did adding these factors improve the model?
comparison <- anova(m.will.select.ALL,m.will.select.ALL.agegender)
anova_tab <- as.data.frame(comparison)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "Comparison of social motivation models with neurotypes: with and without age+gender+language predictors.",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | npar | AIC | BIC | logLik | -2*log(L) | Chisq | Df | Pr(>Chisq) |
|---|---|---|---|---|---|---|---|---|
| m.will.select.ALL | 17 | 1828.084 | 1884.07 | -897.042 | 1794.084 | NA | NA | NA |
| m.will.select.ALL.agegender | 20 | 1831.654 | 1897.52 | -895.827 | 1791.654 | 2.431 | 3 | 0.488 |
No.
Anova:
a.m.will.select.ALL.agegender = anova(m.will.select.ALL.agegender)
anova_tab <- as.data.frame(a.m.will.select.ALL.agegender)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the age+gender+language willingness model with neurotypes",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| Group | 749.872 | 749.872 | 1 | 55.789 | 4.025 | 0.050 |
| S_Q_PC1 | 1055.004 | 1055.004 | 1 | 59.358 | 5.662 | 0.021 |
| S_Q_PC2 | 455.623 | 455.623 | 1 | 64.570 | 2.445 | 0.123 |
| S_Q_PC3 | 503.141 | 503.141 | 1 | 94.749 | 2.700 | 0.104 |
| partner_neurotype | 18.100 | 18.100 | 1 | 44.311 | 0.097 | 0.757 |
| age_diff | 60.831 | 60.831 | 1 | 106.751 | 0.326 | 0.569 |
| gender_match | 318.932 | 318.932 | 1 | 94.529 | 1.712 | 0.194 |
| language | 133.439 | 133.439 | 1 | 81.036 | 0.716 | 0.400 |
| Group:partner_neurotype | 207.763 | 207.763 | 1 | 110.653 | 1.115 | 0.293 |
| S_Q_PC1:partner_neurotype | 99.413 | 99.413 | 1 | 115.254 | 0.534 | 0.467 |
| S_Q_PC2:partner_neurotype | 1653.394 | 1653.394 | 1 | 131.172 | 8.874 | 0.003 |
| S_Q_PC3:partner_neurotype | 291.337 | 291.337 | 1 | 99.296 | 1.564 | 0.214 |
| Group:S_Q_PC1 | 310.221 | 310.221 | 1 | 62.538 | 1.665 | 0.202 |
| Group:S_Q_PC2 | 2.522 | 2.522 | 1 | 69.207 | 0.014 | 0.908 |
| Group:S_Q_PC3 | 135.451 | 135.451 | 1 | 96.899 | 0.727 | 0.396 |
By DYAD TYPE
Rapport ~ predictors x dyad (LMM)
Model
Prepare the lmer formula:
main_predictors <- selected_predictor_summary$predictor
main_predictors <- main_predictors[!main_predictors == "Group"]
interaction_predictors <- main_predictors
fixed_formula_text <- paste0(
"R_Int_willingness ~ dyad + ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":dyad"), collapse = " + ")
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_willingness ~ dyad + S_Q_PC1 + S_Q_PC2 + S_Q_PC3 + S_Q_PC1:dyad +
## S_Q_PC2:dyad + S_Q_PC3:dyad + (1 | ID) + (1 | partner_ID) +
## (1 | dyad_indx)Fit the model:
m.will.select.dyad <- lmer(
formula = mixed_formula,
data = DFwill,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Performance
Singular fit?
print(isSingular(m.will.select.dyad, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.will.select.dyad), residuals(m.will.select.dyad), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.will.select.dyad), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.will.select.dyad), main='')
qqline(residuals(m.will.select.dyad))
layout(1)Variance components:
print(VarCorr(m.will.select.dyad))## Groups Name Std.Dev.
## dyad_indx (Intercept) 10.3035
## partner_ID (Intercept) 9.9845
## ID (Intercept) 14.8422
## Residual 14.0501Performance summary:
performance::model_performance(m.will.select.dyad)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ---------------------------------------------------------------------------
## 1780.5 | 1783.5 | 1833.2 | 0.730 | 0.146 | 0.683 | 9.905 | 14.050Interpretation
Anova:
a.m.will.select.dyad = anova(m.will.select.dyad)
anova_tab <- as.data.frame(a.m.will.select.dyad)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the social motivation model with dyad type",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| dyad | 731.692 | 365.846 | 2 | 134.936 | 1.853 | 0.161 |
| S_Q_PC1 | 457.999 | 457.999 | 1 | 102.160 | 2.320 | 0.131 |
| S_Q_PC2 | 1337.912 | 1337.912 | 1 | 84.564 | 6.777 | 0.011 |
| S_Q_PC3 | 822.746 | 822.746 | 1 | 99.360 | 4.168 | 0.044 |
| dyad:S_Q_PC1 | 129.826 | 64.913 | 2 | 147.261 | 0.329 | 0.720 |
| dyad:S_Q_PC2 | 1489.124 | 744.562 | 2 | 150.460 | 3.772 | 0.025 |
| dyad:S_Q_PC3 | 167.794 | 83.897 | 2 | 129.145 | 0.425 | 0.655 |
Plots:
Plot all significant results:
p_combined_ALL_motivation_dyad <- plot_lmm_significant(
model = m.will.select.dyad,
anova_obj = a.m.will.select.dyad,
df = DFwill,
grouping_var = "dyad",
DV = "R_Int_willingness",
n_cols_plotting = 3,
show_data = T,
colour_by = "dyad"
)
p_combined_ALL_motivation_dyad
Check the dyad comparisons for autism attitudes
Against 0:
# Estimated S_Q_PC2 slope within each dyad type
tr <- emtrends(
m.will.select.dyad,
~ dyad,
var = "S_Q_PC2"
)Against 0:
# Test whether each dyad-specific slope differs from zero
summary(tr, infer = c(TRUE, TRUE), adjust = "holm")## dyad S_Q_PC2.trend SE df lower.CL upper.CL t.ratio p.value
## AUT 12.755 5.65 197.2 -0.892 26.40 2.257 0.0502
## MIX 7.779 2.79 98.6 0.989 14.57 2.790 0.0190
## NAUT 0.312 3.41 126.4 -7.955 8.58 0.092 0.9272
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: holm method for 3 testsPairwise:
# Pairwise tests: do the S_Q_PC2 slopes differ across dyad types?
pairs(tr, adjust = "holm")## contrast estimate SE df t.ratio p.value
## AUT - MIX 4.98 5.74 207 0.867 0.3869
## AUT - NAUT 12.44 6.33 207 1.967 0.1010
## MIX - NAUT 7.47 3.15 141 2.371 0.0573
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: holm method for 3 testsNo pair-wise comparisons between the dyads survive corrections.
See effect sizes:
get_anova_effect_sizes(m.will.select.dyad)## Parameter Eta2_partial CI CI_low CI_high Omega2_partial Omega2_CI_low Omega2_CI_high
## 1 dyad 0.03 0.95 0.00 1 0.01 0.00 1
## 2 S_Q_PC1 0.02 0.95 0.00 1 0.01 0.00 1
## 3 S_Q_PC2 0.07 0.95 0.01 1 0.06 0.01 1
## 4 S_Q_PC3 0.04 0.95 0.00 1 0.03 0.00 1
## 5 dyad:S_Q_PC1 0.00 0.95 0.00 1 0.00 0.00 1
## 6 dyad:S_Q_PC2 0.05 0.95 0.00 1 0.03 0.00 1
## 7 dyad:S_Q_PC3 0.01 0.95 0.00 1 0.00 0.00 1#get_effect_sizes(m.will.select.dyad)See BF:
INTER: expressivity
model_null <- update(m.will.select.dyad, . ~ . - fam_synchrony_expressivity - fam_synchrony_expressivity:dyad)
compare_models_BF(model_null = model_null,
model_full = m.will.select.dyad,
model_null_name = "Without INTER: expressivity",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without INTER: expressivity model BIC: 1885.31
## - Full model model BIC: 1885.31
## - BIC difference (Full - Null): 0
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without INTER: expressivity vs Full model ): 1
## - BF₁₀ (Evidence for Full model vs Without INTER: expressivity ): 1
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Full model
## - Evidence strength: No evidence
## - Bayes Factor: 1 :1 in favor of Full model
##
## CONCLUSION: There is no evidence evidence ( 1 :1) in favor
## of the Full model model over the Without INTER: expressivity model.
## =======================================================================Q1:
model_null <- update(m.will.select.dyad, . ~ . - S_Q_PC1 - S_Q_PC1:dyad)
compare_models_BF(model_null = model_null,
model_full = m.will.select.dyad,
model_null_name = "Without autism-ike traits (PC1)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without autism-ike traits (PC1) model BIC: 1871.78
## - Full model model BIC: 1885.31
## - BIC difference (Full - Null): 13.53
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without autism-ike traits (PC1) vs Full model ): 868
## - BF₁₀ (Evidence for Full model vs Without autism-ike traits (PC1) ): 0.0012
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without autism-ike traits (PC1)
## - Evidence strength: Decisive/Extreme
## - Bayes Factor: 867.72 :1 in favor of Without autism-ike traits (PC1)
##
## CONCLUSION: There is decisive/extreme evidence ( 867.72 :1) in favor
## of the Without autism-ike traits (PC1) model over the Full model model.
## =======================================================================Q3:
model_null <- update(m.will.select.dyad, . ~ . - S_Q_PC3 - S_Q_PC3:dyad)
compare_models_BF(model_null = model_null,
model_full = m.will.select.dyad,
model_null_name = "Without cognitive engagement (PC3)",
model_full_name = "Full model")## =======================================================================
## BAYES FACTOR MODEL COMPARISON
## =======================================================================
##
## Model Information:
## - Without cognitive engagement (PC3) model BIC: 1873.46
## - Full model model BIC: 1885.31
## - BIC difference (Full - Null): 11.84
##
## ------------------------------------------------------------------------
## Bayes Factors:
## - BF₀₁ (Evidence for Without cognitive engagement (PC3) vs Full model ): 373
## - BF₁₀ (Evidence for Full model vs Without cognitive engagement (PC3) ): 0.0027
##
## ------------------------------------------------------------------------
## Interpretation:
## - Preferred model: Without cognitive engagement (PC3)
## - Evidence strength: Decisive/Extreme
## - Bayes Factor: 373.33 :1 in favor of Without cognitive engagement (PC3)
##
## CONCLUSION: There is decisive/extreme evidence ( 373.33 :1) in favor
## of the Without cognitive engagement (PC3) model over the Full model model.
## =======================================================================With gender + age + language
main_predictors <- selected_predictor_summary$predictor
main_predictors <- main_predictors[!main_predictors=="Group"]
interaction_predictors <- main_predictors
fixed_formula_text <- paste0(
"R_Int_willingness ~ dyad + ",
paste(main_predictors, collapse = " + "), " + ",
paste(paste0(interaction_predictors, ":dyad"), collapse = " + ")#,
# "+",
# paste(paste0(interaction_predictors, ":age_diff"), collapse = " + "),"+",
# paste(paste0(interaction_predictors, ":gender_match"), collapse = " + "),"+",
# paste(paste0(interaction_predictors, ":language"), collapse = " + ")
)
mixed_formula <- as.formula(
paste0(
fixed_formula_text,
" + age_diff + gender_match + language + (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)"
)
)
print(mixed_formula)## R_Int_willingness ~ dyad + S_Q_PC1 + S_Q_PC2 + S_Q_PC3 + S_Q_PC1:dyad +
## S_Q_PC2:dyad + S_Q_PC3:dyad + age_diff + gender_match + language +
## (1 | ID) + (1 | partner_ID) + (1 | dyad_indx)Fit the model:
m.will.select.dyad.agegender <- lmer(
formula = mixed_formula,
data = DFwill,
REML = FALSE,
control = lmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
)Singular fit?
print(isSingular(m.will.select.dyad.agegender, tol = 1e-4))## [1] FALSEAssumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(m.will.select.dyad.agegender), residuals(m.will.select.dyad.agegender), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(m.will.select.dyad.agegender), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(m.will.select.dyad.agegender), main='')
qqline(residuals(m.will.select.dyad.agegender))
layout(1)Variance components:
print(VarCorr(m.will.select.dyad.agegender))## Groups Name Std.Dev.
## dyad_indx (Intercept) 10.2636
## partner_ID (Intercept) 9.5478
## ID (Intercept) 14.4996
## Residual 14.1361Performance summary:
performance::model_performance(m.will.select.dyad.agegender)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ----------------------------------------------------------------------------
## 1773.3 | 1777.5 | 1835.9 | 0.725 | 0.166 | 0.671 | 10.038 | 14.136Did adding these factors improve the model?
comparison <- anova(m.will.select.dyad,m.will.select.dyad.agegender)
anova_tab <- as.data.frame(comparison)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "Comparison of social motivation models with dyad types: with and without age+gender predictors.",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | npar | AIC | BIC | logLik | -2*log(L) | Chisq | Df | Pr(>Chisq) |
|---|---|---|---|---|---|---|---|---|
| m.will.select.dyad | 16 | 1832.615 | 1885.308 | -900.308 | 1800.615 | NA | NA | NA |
| m.will.select.dyad.agegender | 19 | 1835.544 | 1898.116 | -898.772 | 1797.544 | 3.072 | 3 | 0.381 |
No.
Anova:
a.m.will.select.dyad.agegender = anova(m.will.select.dyad.agegender)
anova_tab <- as.data.frame(a.m.will.select.dyad.agegender)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for the age+gender+language social motivation model with dyad type",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| dyad | 610.303 | 305.151 | 2 | 138.964 | 1.527 | 0.221 |
| S_Q_PC1 | 374.902 | 374.902 | 1 | 104.473 | 1.876 | 0.174 |
| S_Q_PC2 | 1575.653 | 1575.653 | 1 | 85.468 | 7.885 | 0.006 |
| S_Q_PC3 | 781.066 | 781.066 | 1 | 100.556 | 3.909 | 0.051 |
| age_diff | 55.426 | 55.426 | 1 | 100.904 | 0.277 | 0.600 |
| gender_match | 329.034 | 329.034 | 1 | 91.164 | 1.647 | 0.203 |
| language | 303.981 | 303.981 | 1 | 77.251 | 1.521 | 0.221 |
| dyad:S_Q_PC1 | 83.824 | 41.912 | 2 | 146.301 | 0.210 | 0.811 |
| dyad:S_Q_PC2 | 1475.942 | 737.971 | 2 | 150.001 | 3.693 | 0.027 |
| dyad:S_Q_PC3 | 201.234 | 100.617 | 2 | 128.346 | 0.504 | 0.606 |
Additional
Neurotype guessed
Model
DFoutcomes$partner_neurotype <- ifelse(
DFoutcomes$dyad == "MIX" & DFoutcomes$Group == "NAUT", "AUT",
ifelse(
DFoutcomes$dyad == "MIX" & DFoutcomes$Group == "AUT", "NAUT",
as.character(DFoutcomes$Group)
)
)
DFoutcomes$partner_neurotype <- factor(
DFoutcomes$partner_neurotype,
levels = c("AUT", "NAUT")
)
neurotype.guess <- lmer(neurotype_guessed ~ Group * partner_neurotype + (1|ID) + (1|partner_ID), DFoutcomes, REML = F)
neurotype.guess@call## lmer(formula = neurotype_guessed ~ Group * partner_neurotype +
## (1 | ID) + (1 | partner_ID), data = DFoutcomes, REML = F)Performance
Singular fit?
print(isSingular(neurotype.guess, tol = 1e-4))## [1] TRUEThe fit is singular. We check the variance in random effects:
print(VarCorr(neurotype.guess))## Groups Name Std.Dev.
## ID (Intercept) 0.00000
## partner_ID (Intercept) 0.16253
## Residual 0.39480ID did not explain much variance. We remove it.
neurotype.guess <- lmer(neurotype_guessed ~ Group * partner_neurotype + (1|partner_ID), DFoutcomes, REML = F)
neurotype.guess@call## lmer(formula = neurotype_guessed ~ Group * partner_neurotype +
## (1 | partner_ID), data = DFoutcomes, REML = F)print(isSingular(neurotype.guess, tol = 1e-4))## [1] FALSEThis works.
Variance components again:
print(VarCorr(neurotype.guess))## Groups Name Std.Dev.
## partner_ID (Intercept) 0.16253
## Residual 0.39480Assumptions:
layout(matrix(c(1,1,2,3), 2, 2))
plot(fitted(neurotype.guess), residuals(neurotype.guess), main='', xlab='Fitted Values', ylab='Residuals', abline(h=0, lty=2))
hist(residuals(neurotype.guess), main='', xlab='Residuals', breaks = 30)
qqnorm(residuals(neurotype.guess), main='')
qqline(residuals(neurotype.guess))
layout(1)Performance summary:
performance::model_performance(neurotype.guess)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------
## 257.8 | 258.2 | 277.8 | 0.258 | 0.133 | 0.145 | 0.374 | 0.395Interpretation
Anova:
a.neurotype.guess <- anova(neurotype.guess)
anova_tab <- as.data.frame(a.neurotype.guess)
anova_tab$Effect <- rownames(anova_tab)
rownames(anova_tab) <- NULL
anova_tab <- anova_tab[, c("Effect", setdiff(names(anova_tab), "Effect"))]
kable(anova_tab,
digits = 3,
caption = "ANOVA table for guessing neurotype of interaction partner",
align = "lrrrrrr") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"))| Effect | Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) |
|---|---|---|---|---|---|---|
| Group | 0.540 | 0.540 | 1 | 158.900 | 3.463 | 0.065 |
| partner_neurotype | 0.070 | 0.070 | 1 | 57.831 | 0.448 | 0.506 |
| Group:partner_neurotype | 5.068 | 5.068 | 1 | 158.900 | 32.512 | 0.000 |
Interaction effect:
emm <- emmeans(neurotype.guess, ~ Group * partner_neurotype)
simple_effects <- rbind(
pairs(emm, by = "partner_neurotype", adjust = "none"),
pairs(emm, by = "Group", adjust = "none")
)
summary(simple_effects, adjust = "holm")## partner_neurotype Group contrast estimate SE df t.ratio p.value
## AUT . AUT - NAUT 0.416 0.0783 159 5.308 <0.0001
## NAUT . AUT - NAUT -0.211 0.0783 159 -2.696 0.0078
## . AUT AUT - NAUT 0.267 0.0903 130 2.953 0.0075
## . NAUT AUT - NAUT -0.360 0.0903 130 -3.991 0.0003
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: holm method for 4 testsPlot:
dodge <- position_dodge(width = 0.18)
# Model-estimated marginal means and 95% confidence intervals
emm_guess_df <- emmeans(
neurotype.guess,
~ Group * partner_neurotype
) %>%
summary(infer = c(TRUE, TRUE)) %>%
as.data.frame() %>%
dplyr::mutate(
Group = factor(Group, levels = c("AUT", "NAUT")),
partner_neurotype = factor(partner_neurotype, levels = c("AUT", "NAUT"))
)
p.neurotype_guess <- ggplot(
emm_guess_df,
aes(
x = partner_neurotype,
y = emmean,
group = Group,
colour = Group
)
) +
geom_line(
position = dodge,
linewidth = 1.2
) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = dodge,
width = 0.06,
linewidth = 0.9
) +
geom_point(
position = dodge,
size = 4.5
) +
scale_colour_manual(values = colours_neurotypes) +
scale_y_continuous(
breaks = c(0, 0.25, 0.50, 0.75, 1.00),
labels = scales::label_percent(accuracy = 1)
) +
coord_cartesian(ylim = c(0, 1.05)) +
labs(
x = "Partner neurotype\n",
y = "Correct guess",
colour = "Rater neurotype"
) +
my_theme
p.neurotype_guess
Plots
text_size <- 15Outcomes across dyad types
Rapport and social motivation across dyad types:
p.rapport + p.motivation + p.neurotype_guess
p.neurotype_guess_fixed <- p.neurotype_guess +
labs(color = "Participant", fill = "Participant")
combined_outcomes <- (
(p.rapport | p.motivation | p.neurotype_guess_fixed) / guide_area()
) +
plot_layout(
heights = c(1, 0.15),
widths = c(1.1, 1.1, 0.95),
guides = "collect"
) &
theme(
text = element_text(size = text_size),
axis.title = element_text(size = text_size),
axis.title.x = element_text(size = text_size, margin = margin(t = 6, b = 2)),
axis.title.y = element_text(size = text_size, margin = margin(r = 5, l = 2)), # Keeps Y-title close to axis
axis.text = element_text(size = text_size - 2, color = "black"),
axis.text.x = element_text(size = text_size - 2, color = "black"),
axis.text.y = element_text(size = text_size - 2, color = "black"),
legend.title = element_text(size = text_size),
legend.text = element_text(size = text_size - 2),
plot.margin = margin(t = 8, r = 8, b = 8, l = 8)
)
combined_outcomes
Predictors of rapport and social motivation:
rapp.dyad <- plot_lmm_significant(
model = m.rapp.select.dyad,
anova_obj = a.m.rapp.select.dyad,
df = DFrapp,
grouping_var = "dyad",
n_cols_plotting = 3,
show_data = T,
colour_by="dyad" #dyad, Group, none
)
rapp.neurotype <- plot_lmm_significant(
model = m.rapp.select.ALL,
anova_obj = a.m.rapp.select.ALL,
df = DFrapp,
grouping_var = "partner_neurotype",
n_cols_plotting = 3,
show_data = T,
colour_by = "Group"
)
will.dyad <- plot_lmm_significant(
model = m.will.select.dyad,
anova_obj = a.m.will.select.dyad,
df = DFwill,
grouping_var = "dyad",
DV = "R_Int_willingness",
n_cols_plotting = 3,
show_data = T,
colour_by = "dyad"
)
will.neurotype <- plot_lmm_significant(
model = m.will.select.ALL,
anova_obj = a.m.will.select.ALL,
df = DFwill,
grouping_var = "partner_neurotype",
DV = "R_Int_willingness",
n_cols_plotting = 3,
show_data = T,
colour_by = "Group"
)Rapport combined
# ── helper: strip y-axis title from all panels except the first ──────────
no_y <- theme(axis.title.y = element_blank())# Add model label to y-axis of first panel in each row
rapp.dyad[[1]] <- rapp.dyad[[1]] + labs(y = "Rapport\n(model with dyad type)")
rapp.neurotype[[1]] <- rapp.neurotype[[1]] + labs(y = "Rapport\n(model with neurotypes)")
rapp.dyad[[1]] <- rapp.dyad[[1]] + labs(x = "Non-verbal synchrony")
rapp.neurotype[[1]] <- rapp.neurotype[[1]] + labs(x = "Non-verbal synchrony")
# Remove y-axis title from all other panels
for (i in 2:3) rapp.dyad[[i]] <- rapp.dyad[[i]] + no_y
for (i in 2:3) rapp.neurotype[[i]] <- rapp.neurotype[[i]] + no_y
# Remove per-row legends, collect into one shared guide area
combined_rapport <- (rapp.dyad & theme(legend.position = "none")) /
(rapp.neurotype & theme(legend.position = "none")) /
guide_area() +
plot_layout(heights = c(1, 1, 0.15), guides = "collect") &
theme(
text = element_text(size = text_size), # all text
axis.title = element_text(size = text_size),
axis.text = element_text(size = text_size-2),
legend.text = element_text(size = text_size-2),
legend.title = element_text(size = text_size)
)
combined_rapport
Social motivation combined
will.dyad[[1]] <- will.dyad[[1]] + labs(y = "Social motivation\n(model with dyad type)")
will.neurotype[[1]] <- will.neurotype[[1]] + labs(y = "Social motivation\n(model with neurotypes)")
for (i in 2:3) will.dyad[[i]] <- will.dyad[[i]] + no_y
for (i in 2:3) will.neurotype[[i]] <- will.neurotype[[i]] + no_y
combined_motivation <- (will.dyad & theme(legend.position = "none")) /
(will.neurotype & theme(legend.position = "none")) /
guide_area() +
plot_layout(heights = c(1, 1, 0.15), guides = "collect") &
theme(
text = element_text(size = text_size),
axis.title = element_text(size = text_size),
axis.text = element_text(size = text_size-2),
legend.text = element_text(size = text_size-2),
legend.title = element_text(size = text_size)
)
combined_motivation
Predictor selection: bootstrap results
# 1. Processing and cleaning function
process_boot_data <- function(df, outcome_name) {
df %>%
mutate(
Outcome = outcome_name,
Sign_Stability = pmax(prop_positive, prop_negative) * 100,
Dominant_Sign = ifelse(prop_positive >= prop_negative, "+", "–"),
# Clean variable names with balanced multi-line wrapping
Var_Label = variable %>%
gsub(":GroupNAUT", "\n× neurotype", .) %>%
gsub("GroupNAUT", "Rater's neurotype", .) %>%
gsub("S_Q_PC1", "Autism-like traits", .) %>%
gsub("S_Q_PC2", "Attitudes towards\nautism", .) %>%
gsub("S_Q_PC3", "Cognitive engagement", .) %>%
gsub("fam_synchrony_expressivity", "Non-verbal\nsynchrony", .) %>%
gsub("fam_backchannel_tpt", "Back-\nchannelling", .) %>%
gsub("fam_intonation", "Intonation", .) %>%
gsub("fam_expressivity", "Expressivity", .) %>%
gsub("fam_turn_taking", "Turn taking", .),
# Categorize predictor type
Predictor_Type = case_when(
grepl("Autism-like|Attitudes|Cognitive", Var_Label) ~ "Disposition",
grepl("synchrony|channelling|Intonation|Expressivity|Turn taking", Var_Label) ~ "Behaviour",
TRUE ~ "Rater's neurotype"
)
)
}
# 2. Process datasets
df_rapport <- process_boot_data(boot_summary_rapport, "Rapport")
df_socmot <- process_boot_data(boot_summary_motivation, "Social Motivation")
df_all <- bind_rows(df_rapport, df_socmot) %>%
mutate(
Outcome = factor(Outcome, levels = c("Rapport", "Social Motivation")),
Predictor_Type = factor(Predictor_Type, levels = c("Disposition", "Behaviour", "Rater's neurotype"))
)
# 3. High-contrast palette
category_colors <- c(
"Disposition" = "#FFA941",
"Behaviour" = "#4EB3D3",
"Rater's neurotype" = "grey35"
)
text_colors <- c(
"Disposition" = "#B35A00", # Rich amber
"Behaviour" = "#136994", # Deep teal-blue
"Rater's neurotype" = "grey25"
)
# 4. Generate Plot
p_quadrant_modern <- ggplot(df_all, aes(x = mean_abs_coef, y = Sign_Stability)) +
# Pure white selected quadrant
annotate("rect", xmin = 1.5, xmax = Inf, ymin = 80, ymax = 101,
fill = "#FFFFFF", color = NA) +
# Selection threshold reference lines
geom_vline(xintercept = 1.5, linetype = "dashed", color = "grey30", linewidth = 0.7) +
geom_hline(yintercept = 80, linetype = "dashed", color = "grey30", linewidth = 0.7) +
# Point badges: colored discs
geom_point(
aes(fill = Predictor_Type),
shape = 21,
size = 7.0,
color = "white",
stroke = 1.2
) +
# Stamped bold '+' and '–' inside discs
geom_text(
aes(label = Dominant_Sign),
size = 4.6,
fontface = "bold",
color = "black",
vjust = 0.42
) +
# Repelled labels without connecting lines
geom_text_repel(
aes(label = Var_Label, color = Predictor_Type),
size = 3.6,
fontface = "bold",
lineheight = 0.85,
box.padding = 0.45,
point.padding = 0.40,
force = 6,
force_pull = 0.5,
min.segment.length = Inf, # Suppresses connecting lines entirely
max.overlaps = Inf,
show.legend = FALSE
) +
# Scales & axes
scale_fill_manual(values = category_colors, name = "Predictor Category") +
scale_color_manual(values = text_colors) +
scale_x_continuous(
breaks = seq(0.5, 3.5, by = 0.5),
expand = expansion(mult = c(0.04, 0.08))
) +
scale_y_continuous(
limits = c(50, 101),
breaks = seq(50, 100, by = 10),
expand = c(0, 0)
) +
facet_wrap(~ Outcome, nrow = 1) +
labs(
x = "Mean Absolute Coefficient (|β| across 500 resamples)",
y = "Sign Stability (% resamples in dominant direction)"
) +
theme_minimal(base_size = 13) +
theme(
# Clean floating titles
strip.text = element_text(face = "bold", size = 16, color = "grey15", margin = margin(b = 10)),
strip.background = element_blank(),
# Clean panels
panel.background = element_rect(fill = "#F2F4F7", color = NA),
panel.border = element_rect(color = "grey75", fill = NA, linewidth = 0.8),
panel.grid.major = element_line(color = "white", linewidth = 0.6),
panel.grid.minor = element_blank(),
panel.spacing = unit(1.4, "lines"),
# Typography & axes
axis.title = element_text(face = "bold", size = 12.5, color = "grey15"),
axis.text = element_text(size = 11, color = "grey25"),
axis.title.x = element_text(margin = margin(t = 10)),
axis.title.y = element_text(margin = margin(r = 10)),
# Legend
legend.position = "bottom",
legend.title = element_text(face = "bold", size = 12),
legend.text = element_text(size = 11),
legend.margin = margin(t = 12)
) +
guides(
fill = guide_legend(override.aes = list(size = 6, stroke = 1))
)
# 5. Display & Save
print(p_quadrant_modern)
if (!dir.exists("./figures")) dir.create("./figures", recursive = TRUE)
ggsave(
filename = "./figures/bootstrap_feature_screening.png",
plot = p_quadrant_modern,
width = 13.5,
height = 7.2,
units = "in",
dpi = 300
)All results in one plot
# ── 1. Global typography parameter ─────────────────────────────────────────
text_size <- 18
# ── 2. Helper functions ───────────────────────────────────────────────────
# Strip y-axis title from columns 2 & 3
# Update clean_col23 to enforce y-limits:
clean_col23 <- function(p) {
p +
coord_cartesian(ylim = c(0, 100)) +
labs(y = NULL, title = NULL, subtitle = NULL) +
theme(
axis.title.y = element_blank(),
axis.text.y = element_text(size = text_size - 3, color = "black")
)
}
# Add coord_cartesian(ylim = c(0, 100)) to column 1 plots (rapp.dyad[[1]], rapp.neurotype[[1]], will.dyad[[1]], will.neurotype[[1]]):
rapp.dyad[[1]] <- rapp.dyad[[1]] + coord_cartesian(ylim = c(0, 100))
rapp.neurotype[[1]] <- rapp.neurotype[[1]] + coord_cartesian(ylim = c(0, 100))
will.dyad[[1]] <- will.dyad[[1]] + coord_cartesian(ylim = c(0, 100))
will.neurotype[[1]] <- will.neurotype[[1]] + coord_cartesian(ylim = c(0, 100))
# Make significance labels (ns, **) smaller and prevent top clipping
fix_row1_signif <- function(p, target_text_size = 3.2) {
p <- p +
coord_cartesian(ylim = c(0, 122), clip = "off") +
theme(plot.margin = margin(t = 12, r = 6, b = 6, l = 6))
for (i in seq_along(p$layers)) {
if (any(c("GeomSignif", "GeomText") %in% class(p$layers[[i]]$geom))) {
p$layers[[i]]$aes_params$size <- target_text_size
}
}
return(p)
}
# ── 3. Block 1: Row 1 (Small "ns" + fixed margins) ────────────────────────
p.rapport_fixed <- fix_row1_signif(p.rapport, target_text_size = 4)
p.motivation_fixed <- fix_row1_signif(p.motivation, target_text_size = 4)
p.neurotype_guess_fixed <- p.neurotype_guess +
labs(colour = "Participant", fill = "Participant") +
theme(plot.margin = margin(t = 12, r = 6, b = 6, l = 6))
row1_outcomes <- (p.rapport_fixed | p.motivation_fixed | p.neurotype_guess_fixed)
# ── 4. Block 2: Rapport Models (Rows 2 & 3: Clean "Rapport" title) ─────────
rapp.dyad[[1]] <- rapp.dyad[[1]] +
labs(y = "Rapport", x = "Non-verbal synchrony", title = NULL) +
theme(axis.title.y = element_text(angle = 90, vjust = 0.5))
rapp.neurotype[[1]] <- rapp.neurotype[[1]] +
labs(y = "Rapport", x = "Non-verbal synchrony", title = NULL) +
theme(axis.title.y = element_text(angle = 90, vjust = 0.5))
for (i in 2:3) {
rapp.dyad[[i]] <- clean_col23(rapp.dyad[[i]])
rapp.neurotype[[i]] <- clean_col23(rapp.neurotype[[i]])
}
row2_rapp_dyad <- wrap_plots(rapp.dyad, nrow = 1)
row3_rapp_neuro <- wrap_plots(rapp.neurotype, nrow = 1)
# ── 5. Block 3: Social Motivation Models (Rows 4 & 5: Clean title) ────────
will.dyad[[1]] <- will.dyad[[1]] +
labs(y = "Social motivation", title = NULL) +
theme(axis.title.y = element_text(angle = 90, vjust = 0.5))
will.neurotype[[1]] <- will.neurotype[[1]] +
labs(y = "Social motivation", title = NULL) +
theme(axis.title.y = element_text(angle = 90, vjust = 0.5))
for (i in 2:3) {
will.dyad[[i]] <- clean_col23(will.dyad[[i]])
will.neurotype[[i]] <- clean_col23(will.neurotype[[i]])
}
row4_will_dyad <- wrap_plots(will.dyad, nrow = 1)
row5_will_neuro <- wrap_plots(will.neurotype, nrow = 1)
# ── 6. Assemble Main Mega-Figure WITHOUT Legends ──────────────────────────
mega_figure <- (
row1_outcomes /
plot_spacer() / # Clean gap between outcomes and rapport
row2_rapp_dyad /
row3_rapp_neuro /
plot_spacer() / # Clean gap between rapport and motivation
row4_will_dyad /
row5_will_neuro
) +
plot_layout(
heights = c(1.05, 0.08, 1, 1, 0.08, 1, 1)
) &
theme(
text = element_text(size = text_size),
axis.title = element_text(size = text_size, face = "bold"),
axis.title.x = element_text(size = text_size, margin = margin(t = 6, b = 2)),
axis.title.y = element_text(size = text_size, angle = 90, margin = margin(r = 6, l = 2)),
axis.text = element_text(size = text_size - 3, color = "black"),
axis.text.x = element_text(size = text_size - 3, color = "black"),
axis.text.y = element_text(size = text_size - 3, color = "black"),
# Hide legends on the main panel
legend.position = "none",
plot.margin = margin(t = 5, r = 5, b = 5, l = 5)
)
# ── 7. Build Standalone Legend Plot ───────────────────────────────────────
dummy_legend_plot <- (
p.rapport_fixed +
rapp.dyad[[1]] +
rapp.neurotype[[1]] +
rapp.dyad[[3]] # ensures all color/fill scales are represented
) +
plot_layout(guides = "collect") &
theme(
legend.position = "bottom",
legend.box = "vertical", # Two clean centered rows
legend.box.just = "center",
legend.title = element_text(size = text_size - 2, face = "bold"),
legend.text = element_text(size = text_size - 4),
legend.margin = margin(t = 5, b = 5)
)
standalone_legend <- cowplot::get_legend(dummy_legend_plot)
p_standalone_legend <- cowplot::ggdraw(standalone_legend)
# ── 8. Render & Save Both Objects Separately ──────────────────────────────
if (!dir.exists("./figures")) dir.create("./figures", recursive = TRUE)
# 1. Main 5x3 Mega-Figure
print(mega_figure)
ggsave(
filename = "./figures/mega_figure_interaction_models.png",
plot = mega_figure,
width = 12.0,
height = 15.5,
units = "in",
dpi = 300
)
# 2. Standalone Legend Figure
print(p_standalone_legend)
ggsave(
filename = "./figures/standalone_interaction_legend.png",
plot = p_standalone_legend,
width = 12.0,
height = 2.2,
units = "in",
dpi = 300
)#end
Social motivation across dyads
Is there a willingness difference for dyad types?
Pair-wise comparisons:
Plot:
Performance summary:
See BF:
See effect sizes:
So there is an effect of dyad type on willingness to meet again, but no pairwise comparisons survive the corrections for multiple comparisons.
SRM Model
We first fit a model with no predictors but with 4 random effects:
(1|ID) = random intercept for subject, i.e., the actor effect (1|partner_ID) = random intercept for the interaction partner of the subject, i.e., the partner effect (1|dyad_indx) = random intercept for the dyad, i.e., the relation effect (1|session) = random intercept for the experimental session (in case the groups in different sessions differed systematically)
Variance:
Let’s look at the explained variance in the model:
The random effects explained 68% of variance in the model. We use this SRM to estimate how much variation in willingness is a function of the individuals (autistic or neurotypical) in the role of the actor (the one rating willingness) and of the partner (the one with whom the actor is rating the willingness), and how much of the unique interactions across dyads. Let’s look at that in detail.
Out of the 68% variance explained by the random effects:
actor effects account for
ractor_v% of the variance, *partner effects account forrpartner_v% of the variance, dyad effects account forrdyad_v`% of the variance.Overall, the data suggest that the variance related to the unique relationship between two participants explains more variance than either of the individual effects of actor or partner. This supports the view that willingness is a resultant of a relational, unique combination of the two social agents i addition to individual tendencies to rate or be rated high on willingness.
Social motivation by neurotype in MIX
Performance summary:
See BF:
See effect sizes: