library(nlme) library(plyr) library(dplyr) library(reshape2) library(stringr) library(effsize) library(multgee) library(lmtest) library(modelr) library(rr2) library(parameters) library(effectsize) library(repolr) library(ARTool) library(ggplot2) library(psych) library(qgcomp) library(MASS, include.only = c("polr")) library(ordinal) library(car) # For Levene's test library(nortest) # For Anderson-Darling test as an alternative to Shapiro-Wilk library(PMCMRplus) # For post-hoc tests like Nemenyi test library(FSA) # For Dunn's test library(purrr) library(fmsb) library(multcomp) library(RColorBrewer) library(data.table) sink("RetroSketch-Results.txt") dl = read.csv("RetroSketch-Data.csv", fileEncoding="UTF-8-BOM") dl$Condition = as.factor(dl$Condition) dl$Game = as.factor(dl$Game) dl$Video_game_genre = as.factor(dl$Video_game_genre) # Calculate Valence and Arousal scores from subjective ratings based on Russel's Circumplex Model dl$Retro_Valence4 = ( .592 * dl$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating -.333 * dl$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating -.350 * dl$Exposure_Windowed_60s_Retro_Mean_Fear_Rating +.729 * dl$Exposure_Windowed_60s_Retro_Mean_Joy_Rating) / (dl$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating + dl$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating + dl$Exposure_Windowed_60s_Retro_Mean_Fear_Rating + dl$Exposure_Windowed_60s_Retro_Mean_Joy_Rating) dl$Retro_Arousal4 = ( -.577 * dl$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating -.603 * dl$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating +.679 * dl$Exposure_Windowed_60s_Retro_Mean_Fear_Rating +.108 * dl$Exposure_Windowed_60s_Retro_Mean_Joy_Rating) / (dl$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating + dl$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating + dl$Exposure_Windowed_60s_Retro_Mean_Fear_Rating + dl$Exposure_Windowed_60s_Retro_Mean_Joy_Rating) dl$ESM_Valence4 = ( .592 * dl$ESM_Relaxation_Value -.333 * dl$ESM_Bored_Value -.350 * dl$ESM_Fear_Value +.729 * dl$ESM_Joy_Value) / (dl$ESM_Relaxation_Value + dl$ESM_Bored_Value + dl$ESM_Fear_Value + dl$ESM_Joy_Value) dl$ESM_Arousal4 = ( -.577 * dl$ESM_Relaxation_Value -.603 * dl$ESM_Bored_Value +.679 * dl$ESM_Fear_Value +.108 * dl$ESM_Joy_Value) / (dl$ESM_Relaxation_Value + dl$ESM_Bored_Value + dl$ESM_Fear_Value + dl$ESM_Joy_Value) # Average between left and right side pupil and smile measures dl$Exposure_Windowed_60s_Mean_Pupil_Dilation = (dl$Exposure_Windowed_60s_Mean_Pupil_Dilation_Right + dl$Exposure_Windowed_60s_Mean_Pupil_Dilation_Left) /2 dl$Exposure_Windowed_60s_SD_Pupil_Dilation = (dl$Exposure_Windowed_60s_SD_Pupil_Dilation_Right + dl$Exposure_Windowed_60s_SD_Pupil_Dilation_Left) /2 dl$Exposure_Windowed_60s_Mean_Smile = (dl$Exposure_Windowed_60s_Mean_9_Mouth_Smile_Right + dl$Exposure_Windowed_60s_Mean_10_Mouth_Smile_Left) / 2 # Assign unique id to all rows of a participant based on time dl$Global_Sample_ID = dl$Retro_Sample_Frame + dl$Session_Number * 1801 # Recode VR Experience as ordinal variable dl <- dl %>% mutate(VR_Experience_Rating = case_when( VR_Experience == "Never" ~ 0, VR_Experience == "Occasionally" ~ 1, VR_Experience == "Weekly" ~ 2, VR_Experience == "Daily" ~ 3, TRUE ~ NA_real_ # Default case for any unmatched values )) # Quantize dependent variables to enable ordinal regression qdata <- quantize(data=dl, expnms=c("ESM_Joy_Value"), q=11) dl$ESM_Joy_Value_quant <- qdata$data$ESM_Joy_Value qdata <- quantize(data=dl, expnms=c("ESM_Fear_Value"), q=11) dl$ESM_Fear_Value_quant <- qdata$data$ESM_Fear_Value qdata <- quantize(data=dl, expnms=c("ESM_Relaxation_Value"), q=11) dl$ESM_Relaxation_Value_quant <- qdata$data$ESM_Relaxation_Value qdata <- quantize(data=dl, expnms=c("ESM_Bored_Value"), q=11) dl$ESM_Bored_Value_quant <- qdata$data$ESM_Bored_Value qdata <- quantize(data=dl, expnms=c("ESM_Presence_Value"), q=11) dl$ESM_Presence_Value_quant <- qdata$data$ESM_Presence_Value qdata <- quantize(data=dl, expnms=c("ESM_Arousal4"), q=11) dl$ESM_Arousal4_quant <- qdata$data$ESM_Arousal4 qdata <- quantize(data=dl, expnms=c("ESM_Valence4"), q=11) dl$ESM_Valence4_quant <- qdata$data$ESM_Valence4 qdata <- quantize(data=dl, expnms=c("Exposure_Windowed_60s_Retro_Mean_Joy_Rating"), q=11) dl$Exposure_Windowed_60s_Retro_Mean_Joy_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Joy_Rating qdata <- quantize(data=dl, expnms=c("Exposure_Windowed_60s_Retro_Mean_Fear_Rating"), q=11) dl$Exposure_Windowed_60s_Retro_Mean_Fear_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Fear_Rating qdata <- quantize(data=dl, expnms=c("Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating"), q=11) dl$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating qdata <- quantize(data=dl, expnms=c("Exposure_Windowed_60s_Retro_Mean_Boredom_Rating"), q=11) dl$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating qdata <- quantize(data=dl, expnms=c("Exposure_Windowed_60s_Retro_Mean_Presence_Rating"), q=11) dl$Exposure_Windowed_60s_Retro_Mean_Presence_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Presence_Rating qdata <- quantize(data=dl, expnms=c("Retro_Valence4"), q=11) dl$Retro_Valence4_quant <- qdata$data$Retro_Valence4 qdata <- quantize(data=dl, expnms=c("Retro_Arousal4"), q=11) dl$Retro_Arousal4_quant <- qdata$data$Retro_Arousal4 # Filter the data into subsets d_ESM_Retro = filter(dl, dl$Condition=="ESM_Retro") d_Retro = filter(dl, dl$Condition=="Retro") d_ESM = filter(d_ESM_Retro, !is.na(d_ESM_Retro$ESM_Fear_Value)) # Set correlation structure for linear regression DefaultCorr = corCompSymm(, form= ~ 1 | Participant_ID) # Create data frames that aggregate variables per participant using mean columns_to_aggregate <- c("Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "Retro_Valence4", "Retro_Arousal4", "ESM_Joy_Value", "ESM_Fear_Value", "ESM_Relaxation_Value", "ESM_Bored_Value", "ESM_Presence_Value", "ESM_Valence4", "ESM_Arousal4") # Step 1: Aggregate the columns over Participant_ID aggregated_data <- d_ESM %>% group_by(Participant_ID) %>% summarize(across(all_of(columns_to_aggregate), mean, na.rm = TRUE)) # Step 2: Filter the data to retain one row per Participant_ID (with Retro_Sample_Frame == 300) filtered_data <- d_ESM %>% filter(Retro_Sample_Frame == 300) %>% select(-c(ESM_Fear_Value, all_of(columns_to_aggregate))) %>% distinct(Participant_ID, .keep_all = TRUE) # Ensure only one row per Participant_ID # Step 3: Merge the aggregated data with the filtered data d_ESM_Aggregated <- left_join(filtered_data, aggregated_data, by = "Participant_ID") # Quantize the data qdata <- quantize(data=d_ESM_Aggregated, expnms=c("ESM_Joy_Value"), q=11) d_ESM_Aggregated$ESM_Joy_Value_quant <- qdata$data$ESM_Joy_Value qdata <- quantize(data=d_ESM_Aggregated, expnms=c("ESM_Fear_Value"), q=11) d_ESM_Aggregated$ESM_Fear_Value_quant <- qdata$data$ESM_Fear_Value qdata <- quantize(data=d_ESM_Aggregated, expnms=c("ESM_Relaxation_Value"), q=11) d_ESM_Aggregated$ESM_Relaxation_Value_quant <- qdata$data$ESM_Relaxation_Value qdata <- quantize(data=d_ESM_Aggregated, expnms=c("ESM_Bored_Value"), q=11) d_ESM_Aggregated$ESM_Bored_Value_quant <- qdata$data$ESM_Bored_Value qdata <- quantize(data=d_ESM_Aggregated, expnms=c("ESM_Presence_Value"), q=11) d_ESM_Aggregated$ESM_Presence_Value_quant <- qdata$data$ESM_Presence_Value qdata <- quantize(data=d_ESM_Aggregated, expnms=c("ESM_Arousal4"), q=11) d_ESM_Aggregated$ESM_Arousal4_quant <- qdata$data$ESM_Arousal4 qdata <- quantize(data=d_ESM_Aggregated, expnms=c("ESM_Valence4"), q=11) d_ESM_Aggregated$ESM_Valence4_quant <- qdata$data$ESM_Valence4 qdata <- quantize(data=d_ESM_Aggregated, expnms=c("Exposure_Windowed_60s_Retro_Mean_Joy_Rating"), q=11) d_ESM_Aggregated$Exposure_Windowed_60s_Retro_Mean_Joy_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Joy_Rating qdata <- quantize(data=d_ESM_Aggregated, expnms=c("Exposure_Windowed_60s_Retro_Mean_Fear_Rating"), q=11) d_ESM_Aggregated$Exposure_Windowed_60s_Retro_Mean_Fear_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Fear_Rating qdata <- quantize(data=d_ESM_Aggregated, expnms=c("Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating"), q=11) d_ESM_Aggregated$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating qdata <- quantize(data=d_ESM_Aggregated, expnms=c("Exposure_Windowed_60s_Retro_Mean_Boredom_Rating"), q=11) d_ESM_Aggregated$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating qdata <- quantize(data=d_ESM_Aggregated, expnms=c("Exposure_Windowed_60s_Retro_Mean_Presence_Rating"), q=11) d_ESM_Aggregated$Exposure_Windowed_60s_Retro_Mean_Presence_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Presence_Rating qdata <- quantize(data=d_ESM_Aggregated, expnms=c("Retro_Valence4"), q=11) d_ESM_Aggregated$Retro_Valence4_quant <- qdata$data$Retro_Valence4 qdata <- quantize(data=d_ESM_Aggregated, expnms=c("Retro_Arousal4"), q=11) d_ESM_Aggregated$Retro_Arousal4_quant <- qdata$data$Retro_Arousal4 # Filter binary genders for testing covariate interactions (as not enough data for others) d_ESM_Main_Gender = filter(d_ESM_Aggregated, d_ESM_Aggregated$Gender == "Male" | d_ESM_Aggregated$Gender == "Female") d_ESM_Main_Gender$VR_Experience = as.factor(d_ESM_Main_Gender$VR_Experience) # Step 1: Aggregate the columns over Participant_ID aggregated_data <- d_ESM_Main_Gender %>% group_by(Participant_ID) %>% summarize(across(all_of(columns_to_aggregate), mean, na.rm = TRUE)) # Step 2: Filter the data to retain one row per Participant_ID (with Retro_Sample_Frame == 300) filtered_data <- d_ESM_Main_Gender %>% filter(Retro_Sample_Frame == 300) %>% select(-c(ESM_Fear_Value, all_of(columns_to_aggregate))) %>% distinct(Participant_ID, .keep_all = TRUE) # Ensure only one row per Participant_ID # Step 3: Merge the aggregated data with the filtered data d_ESM_Main_Gender_Aggr <- left_join(filtered_data, aggregated_data, by = "Participant_ID") # # Quantize the data qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("ESM_Joy_Value"), q=11) d_ESM_Main_Gender_Aggr$ESM_Joy_Value_quant <- qdata$data$ESM_Joy_Value qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("ESM_Fear_Value"), q=11) d_ESM_Main_Gender_Aggr$ESM_Fear_Value_quant <- qdata$data$ESM_Fear_Value qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("ESM_Relaxation_Value"), q=11) d_ESM_Main_Gender_Aggr$ESM_Relaxation_Value_quant <- qdata$data$ESM_Relaxation_Value qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("ESM_Bored_Value"), q=11) d_ESM_Main_Gender_Aggr$ESM_Bored_Value_quant <- qdata$data$ESM_Bored_Value qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("ESM_Presence_Value"), q=11) d_ESM_Main_Gender_Aggr$ESM_Presence_Value_quant <- qdata$data$ESM_Presence_Value qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("ESM_Arousal4"), q=11) d_ESM_Main_Gender_Aggr$ESM_Arousal4_quant <- qdata$data$ESM_Arousal4 qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("ESM_Valence4"), q=11) d_ESM_Main_Gender_Aggr$ESM_Valence4_quant <- qdata$data$ESM_Valence4 qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("Exposure_Windowed_60s_Retro_Mean_Joy_Rating"), q=11) d_ESM_Main_Gender_Aggr$Exposure_Windowed_60s_Retro_Mean_Joy_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Joy_Rating qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("Exposure_Windowed_60s_Retro_Mean_Fear_Rating"), q=11) d_ESM_Main_Gender_Aggr$Exposure_Windowed_60s_Retro_Mean_Fear_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Fear_Rating qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating"), q=11) d_ESM_Main_Gender_Aggr$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("Exposure_Windowed_60s_Retro_Mean_Boredom_Rating"), q=11) d_ESM_Main_Gender_Aggr$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Boredom_Rating qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("Exposure_Windowed_60s_Retro_Mean_Presence_Rating"), q=11) d_ESM_Main_Gender_Aggr$Exposure_Windowed_60s_Retro_Mean_Presence_Rating_quant <- qdata$data$Exposure_Windowed_60s_Retro_Mean_Presence_Rating qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("Retro_Valence4"), q=11) d_ESM_Main_Gender_Aggr$Retro_Valence4_quant <- qdata$data$Retro_Valence4 qdata <- quantize(data=d_ESM_Main_Gender_Aggr, expnms=c("Retro_Arousal4"), q=11) d_ESM_Main_Gender_Aggr$Retro_Arousal4_quant <- qdata$data$Retro_Arousal4 # Helper method for pretty-printing p-values with asterisks print_p <- function(p) { starsv <- c("^{***}", "^{**}", "^*", "") vec <- c(0, 0.001, 0.01, 0.05, 1) star_i <- findInterval(p, vec) if (p < 0.001) pstr = "p<.001^{***}" else pstr = sprintf("p=%s%s", str_remove(sprintf("%.3f", p), "^0+"), starsv[star_i]) pstr } ######## Emotion Manipulation Check ######## check_assumptions <- function(dep_var, data) { # Perform ANOVA aov_result <- aov(as.formula(paste(dep_var, "~ Game")), data = data) # Get residuals residuals_aov <- residuals(aov_result) # Check normality with Shapiro-Wilk if sample size is appropriate if (length(residuals_aov) <= 5000) { shapiro_test <- shapiro.test(residuals_aov) cat("\nShapiro-Wilk Normality Test for", dep_var, ":\n") print(shapiro_test) } else { # Use Anderson-Darling test for larger sample sizes ad_test <- ad.test(residuals_aov) cat("\nAnderson-Darling Normality Test for", dep_var, ":\n") print(ad_test) } # Homogeneity of variances check with Levene's test levene_test <- leveneTest(as.formula(paste(dep_var, "~ Game")), data = data) cat("\nLevene's Test for Homogeneity of Variance for", dep_var, ":\n") print(levene_test) } # Function to run non-parametric tests run_nonparametric_tests <- function(dep_var, data) { print("--------------------------------------------") cat("\nDescriptive statistics for", dep_var, ":\n") desc_stats <- aggregate(as.formula(paste(dep_var, "~ Game")), data = data, FUN = mean) print(desc_stats) # Kruskal-Wallis test kruskal_result <- kruskal.test(as.formula(paste(dep_var, "~ Game")), data = data) cat("\nKruskal-Wallis Test for", dep_var, ":\n") print(kruskal_result) # If Kruskal-Wallis test is significant, run Dunn's test for pairwise comparisons if (kruskal_result$p.value < 0.05) { dunn_result <- dunnTest(as.formula(paste(dep_var, "~ Game")), data = data, method = "holm") cat("\nDunn's Test for pairwise comparisons for", dep_var, ":\n") # Round p.adj to three decimal places and add significance stars dunn_result$res$P.adj <- formatC(dunn_result$res$P.adj, format = "f", digits = 3) dunn_result$res$Significance <- ifelse(dunn_result$res$P.adj < 0.001, "***", ifelse(dunn_result$res$P.adj < 0.01, "**", ifelse(dunn_result$res$P.adj < 0.05, "*", ""))) print(dunn_result$res) } } # List of dependent variables (emotion measures) emotion_measures <- c( "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating" ) for (measure in emotion_measures) { check_assumptions(measure, dl) run_nonparametric_tests(measure, dl) } # Calculate mean values for each emotion measure across all participants for each game mean_vals <- aggregate(cbind(Exposure_Windowed_60s_Retro_Mean_Joy_Rating, Exposure_Windowed_60s_Retro_Mean_Fear_Rating, Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating, Exposure_Windowed_60s_Retro_Mean_Boredom_Rating, Exposure_Windowed_60s_Retro_Mean_Presence_Rating) ~ Game, data = dl, FUN = mean) # Define the maximum and minimum values for the radar chart scale max_min <- data.frame( Exposure_Windowed_60s_Retro_Mean_Joy_Rating = c(8, 0), Exposure_Windowed_60s_Retro_Mean_Fear_Rating = c(8, 0), Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating = c(8, 0), Exposure_Windowed_60s_Retro_Mean_Boredom_Rating = c(8, 0), Exposure_Windowed_60s_Retro_Mean_Presence_Rating = c(8, 0) ) rownames(max_min) <- c("Max", "Min") # Remove the 'Game' column from mean_vals mean_vals_clean <- mean_vals[, -1] # Combine the max/min data with the cleaned mean values radar_data <- rbind(max_min, mean_vals_clean) # Set the rownames to be the games rownames(radar_data)[3:nrow(radar_data)] <- mean_vals$Game # Define a mapping for readable column names readable_names <- c( Exposure_Windowed_60s_Retro_Mean_Joy_Rating = "Joy", Exposure_Windowed_60s_Retro_Mean_Fear_Rating = "Fear", Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating = "Relax", Exposure_Windowed_60s_Retro_Mean_Boredom_Rating = "Bored", Exposure_Windowed_60s_Retro_Mean_Presence_Rating = "Presence" ) # Rename columns in radar_data colnames(radar_data) <- readable_names[colnames(radar_data)] # Print the radar_data to verify the new column names print(radar_data) # Set colors for the radar chart colors_border <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728","#FF00FF") # Adjust colors as needed colors_fill <- adjustcolor(colors_border, alpha.f = 0.2) # Plot the radar chart radarchart(radar_data, axistype = 1, pcol = colors_border, pfcol = NA, plwd = 2, plty = 1, cglcol = "grey", cglty = 1, cglwd = 0.8, caxislabels = seq(0, 10, 2), axislabcol = "grey", title = "Emotion Spread Across Games") # Create the game_names vector from the mean_vals$Game column game_names <- mean_vals$Game # Add a legend with proper game names legend(x = "topright", legend = game_names, bty = "n", pch = 20, col = colors_border, text.col = "black", cex = 1.2, pt.cex = 4) ######## RQ1: How do Retro measures relate to ESM measures? ######## # RQ1 correlations between Retro and ESM # Interpret tau correlation effect size interpret_kendalls_tau <- function(kendall_tau){ if (kendall_tau < 0.20) { return("Weak correlation") } else if (kendall_tau < 0.40) { return("Moderate correlation") } else if (kendall_tau < 0.60) { return("Strong correlation") } else if (kendall_tau < 0.80) { return("Very strong correlation") } else { return("Almost perfect correlation") } } run_correlations_with_interval <- function(d1, interval, dv_name, iv_name){ print("############################") cat("Runnning Correlations for Interval:", interval, "\n") d4 = filter(d1, d1$Retro_Sample_Frame == interval) p_ct = cor.test(d4[, dv_name], d4[, iv_name], method="pearson") k_ct = cor.test(d4[, dv_name], d4[, iv_name], method="kendall") print(p_ct) print(interpret_r(p_ct$estimate)) print(k_ct) print(interpret_kendalls_tau(cor(d4[, dv_name], d4[, iv_name], method="kendall"))) ## test assumptions for Pearson correlation tests qqnorm(d4[, dv_name]) qqnorm(d4[, iv_name]) sw_a = shapiro.test(d4[, dv_name]) sw_b = shapiro.test(d4[, iv_name]) plot(d4[, dv_name], d4[, iv_name]) # Decide which correlation to use based on Shapiro-Wilk test if (sw_a$p.value <= 0.05 || sw_b$p.value <= 0.05) { # If data is not normally distributed, use Kendall's correlation print("Returning Kendall's") print("############################") return(list(value = cor(d4[, dv_name], d4[, iv_name], method = "kendall"), method = "kendall")) } else { # If data is normally distributed, use Pearson's correlation print("Returning Pearson's") print("############################") return(list(value = p_ct$estimate, method = "pearson")) } } run_ordinal_regression <- function(formula, d1){ fit <- tryCatch({ ordLORgee( formula = formula, link = "logit", id = d1$Participant_ID, repeated = d1$Retro_Sample_Frame, data = d1, LORstr = "uniform") }, error = function(e) { message("Error in ordLORgee: ", e$message) return(NULL) }) if (!is.null(fit)) { s1 <- summary(fit) return(s1) } else{ return(NULL) } } run_ordinal_regression_using_global_sample <- function(formula, d1){ fit <- tryCatch({ ordLORgee( formula = formula, link = "logit", id = d1$Participant_ID, repeated = d1$Global_Sample_ID, data = d1, LORstr = "uniform") }, error = function(e) { message("Error in ordLORgee: ", e$message) return(NULL) }) if (!is.null(fit)) { s1 <- summary(fit) return(s1) } else{ return(NULL) } } corr_regression <- function(dv_name, iv_name, d1, regtest = TRUE) { # workaround to fix scoping problem f1 <<- as.formula(paste(dv_name, "~", iv_name)) d1 <<- d1 iv_name <<- iv_name pearson_list <- list() kendall_list <- list() if(regtest){ if(dv_name == "ESM_Arousal4" || dv_name == "ESM_Valence4"){ if (dv_name == "ESM_Arousal4") { # Quantize the data qdata <- quantize(data=d1, expnms=c("ESM_Arousal4"), q=11) # Assign the quantized data to the new column in d1 d1$ESM_Arousal4_quant <- qdata$data$ESM_Arousal4 formula <- as.formula(paste("ESM_Arousal4_quant", "~", iv_name, "* Game")) } else{ # Quantize the data qdata <- quantize(data=d1, expnms=c("ESM_Valence4"), q=11) # Assign the quantized data to the new column in d1 d1$ESM_Valence4_quant <- qdata$data$ESM_Valence4 formula <- as.formula(paste("ESM_Valence4_quant", "~", iv_name, "* Game")) } } else{ formula <- as.formula(paste(dv_name, "~", iv_name, "* Game")) } # Initialize an empty data frame to store interactions and reference levels interaction_results <- data.frame(BaseGame = character(), Interaction = character(), PValue = numeric(), stringsAsFactors = FALSE) # First analysis s1 <- run_ordinal_regression(formula, d1) interactions <- s1$coefficients[c(16:19), 4] interaction_names <- rownames(s1$coefficients)[c(16:19)] base_game <- "AC" # Store the results interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Second analysis d1$Game <- relevel(d1$Game, ref = "GotS") s1 <- run_ordinal_regression(formula, d1) #print(s1) interactions <- s1$coefficients[c(17:19), 4] interaction_names <- rownames(s1$coefficients)[c(17:19)] base_game <- "GotS" # Store the results interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Third analysis d1$Game <- relevel(d1$Game, ref = "HLA") s1 <- run_ordinal_regression(formula, d1) interactions <- s1$coefficients[c(18:19), 4] interaction_names <- rownames(s1$coefficients)[c(18:19)] base_game <- "HLA" # Store the results interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Fourth analysis d1$Game <- relevel(d1$Game, ref = "IEYTD") s1 <- run_ordinal_regression(formula, d1) interactions <- s1$coefficients[c(19), 4] interaction_names <- rownames(s1$coefficients)[c(19)] base_game <- "IEYTD" # Store the results interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Adjust p-values interaction_results$AdjustedPValue <- p.adjust(interaction_results$PValue, method = "holm") # Print the results print("ORDINAL REGRESSION INTERACTIONS: ") rownames(interaction_results) <- NULL print(interaction_results, row.names = FALSE) } # correlation overall d2 = aggregate(as.formula(paste(dv_name, "~ Participant_ID")), data = d1, FUN = mean) d3 = aggregate(as.formula(paste(iv_name, "~ Participant_ID")), data = d1, FUN = mean) p_ct = cor.test(d2[, dv_name], d3[, iv_name], method="pearson") k_ct = cor.test(d2[, dv_name], d3[, iv_name], method="kendall") print("Running Correlation overall") print(p_ct) print(interpret_r(p_ct$estimate)) print(k_ct) print(interpret_kendalls_tau(cor(d2[, dv_name], d3[, iv_name], method="kendall"))) ## test assumptions for Pearson correlation tests qqnorm(d2[, dv_name]) qqnorm(d3[, iv_name]) print(shapiro.test(d2[, dv_name])) print(shapiro.test(d3[, iv_name])) plot(d2[, dv_name], d3[, iv_name]) # correlation for ESM Intervals for (interval in seq(300, 1800, by = 300)) { result <- run_correlations_with_interval(d1, interval, dv_name, iv_name) # Store the result in the appropriate list based on the method if (result$method == "pearson") { pearson_list[[length(pearson_list) + 1]] <- result$value } else if (result$method == "kendall") { kendall_list[[length(kendall_list) + 1]] <- result$value } } # Convert lists to numeric vectors pearson_list <- unlist(pearson_list) kendall_list <- unlist(kendall_list) # Print the minimum value from each list if (length(pearson_list) > 0) { cat("Minimum Pearson's r:", min(pearson_list), "\n") print(interpret_r(min(pearson_list))) } else { cat("No Pearson's correlations calculated.\n") } if (length(kendall_list) > 0) { cat("Minimum Kendall's tau:", min(kendall_list), "\n") print(interpret_kendalls_tau(min(kendall_list))) } else { cat("No Kendall's correlations calculated.\n") } print(paste0(dv_name, " ")) } corr_regression("ESM_Joy_Value", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", d_ESM) corr_regression("ESM_Fear_Value", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", d_ESM, regtest = FALSE) corr_regression("ESM_Relaxation_Value", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", d_ESM, regtest = FALSE) corr_regression("ESM_Bored_Value", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", d_ESM, regtest = FALSE) corr_regression("ESM_Presence_Value", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", d_ESM, regtest = FALSE) corr_regression("ESM_Arousal4", "Retro_Arousal4", d_ESM, regtest = FALSE) corr_regression("ESM_Valence4", "Retro_Valence4", d_ESM, regtest = FALSE) # RQ1 differences in distributions between Retro and ESM distr_tests <- function(dv_name, iv_name, d) { calculate_cohen_d <- function(data, iv_name, dv_name) { m1 <- mean(filter(data, Method == iv_name)$Rating) m2 <- mean(filter(data, Method == dv_name)$Rating) sd_pooled <- sqrt(((var(filter(data, Method == iv_name)$Rating) + var(filter(data, Method == dv_name)$Rating)) / 2)) return((m1 - m2) / sd_pooled) } # transform data into long format d_iv = melt(setDT(d), id.vars=c("Participant_ID"), measure.vars=c(iv_name), variable.name="Method", value.name="Rating") d_iv$Game = d$Game d_iv$Retro_Sample_Frame = d$Retro_Sample_Frame d_dv = melt(setDT(d), id.vars=c("Participant_ID"), measure.vars=c(dv_name), variable.name="Method", value.name="Rating") d_dv$Game = d$Game d_dv$Retro_Sample_Frame = d$Retro_Sample_Frame d_long = rbind(d_iv, d_dv) # Mean d_mean = aggregate(Rating ~ Participant_ID + Method + Game, data = d_long, FUN = median) model = art(Rating ~ Method * Game + (1|Participant_ID), data = d_long) # use all data for better power res = anova(model) res$part.eta.sq = res$F * res$Df / (res$F * res$Df + res$Df.res) res$cohen.d = sqrt(4 * res$part.eta.sq / (1 - res$part.eta.sq )) res$effect = interpret_omega_squared(res$part.eta.sq, rules = "cohen1992") res$effect.d = interpret_cohens_d(res$cohen.d, rules = "cohen1988") print(res) print(art.con(model, "Method:Game", adjust="holm")) res_mean = res model_mean = model # SD d_sd = aggregate(Rating ~ Participant_ID + Method + Game, data = d_long, FUN = mad) model = art(Rating ~ Method * Game + (1|Participant_ID), data = d_sd) res = anova(model) res$part.eta.sq = res$F * res$Df / (res$F * res$Df + res$Df.res) res$cohen.d = sqrt(4 * res$part.eta.sq / (1 - res$part.eta.sq )) res$effect = interpret_omega_squared(res$part.eta.sq, rules = "cohen1992") res$effect.d = interpret_cohens_d(res$cohen.d, rules = "cohen1988") print(res) print(art.con(model, "Method:Game", adjust="holm")) res_sd = res model_sd = model # min d_min = aggregate(Rating ~ Participant_ID + Method + Game, data = d_long, FUN = min) model = art(Rating ~ Method * Game + (1|Participant_ID), data = d_min) res = anova(model) res$part.eta.sq = res$F * res$Df / (res$F * res$Df + res$Df.res) res$cohen.d = sqrt(4 * res$part.eta.sq / (1 - res$part.eta.sq )) res$effect = interpret_omega_squared(res$part.eta.sq, rules = "cohen1992") res$effect.d = interpret_cohens_d(res$cohen.d, rules = "cohen1988") print(res) print(art.con(model, "Method:Game", adjust="holm")) res_min = res model_min = model # max d_max = aggregate(Rating ~ Participant_ID + Method + Game, data = d_long, FUN = max) model = art(Rating ~ Method * Game + (1|Participant_ID), data = d_max) res = anova(model) res$part.eta.sq = res$F * res$Df / (res$F * res$Df + res$Df.res) res$cohen.d = sqrt(4 * res$part.eta.sq / (1 - res$part.eta.sq )) res$effect = interpret_omega_squared(res$part.eta.sq, rules = "cohen1992") res$effect.d = interpret_cohens_d(res$cohen.d, rules = "cohen1988") print(res) print(art.con(model, "Method:Game", adjust="holm")) res_max = res model_max = model ## post-hoc correct p-values p_main = p.adjust(c(res_mean$`Pr(>F)`[1], res_sd$`Pr(>F)`[1], res_min$`Pr(>F)`[1], res_max$`Pr(>F)`[1])) p_int = p.adjust(c(res_mean$`Pr(>F)`[3], res_sd$`Pr(>F)`[3], res_min$`Pr(>F)`[3], res_max$`Pr(>F)`[3])) print(paste( "Median ", "p_main=", round(p_main[1], 3), "d=", round(res_mean$cohen.d[1], 3), res_mean$effect.d[1], iv_name, "=", round(mean(filter(d_mean, d_mean$Method == iv_name)$Rating), 3), dv_name, "=", round(mean(filter(d_mean, d_mean$Method == dv_name)$Rating), 3), "p_int=", round(p_int[1], 3))) if (p_int[1] <= .05) { ints <- summary(art.con(model_mean, "Method:Game", adjust="holm")) # Adjust p-values and print contrasts p_ints <- p.adjust(c(ints$p.value[5], ints$p.value[14], ints$p.value[22], ints$p.value[29], ints$p.value[35])) for (i in 1:5) { game <- c("AC", "GotS", "HLA", "IEYTD", "RM")[i] cohen_d <- calculate_cohen_d(filter(d_long, Game == game), iv_name, dv_name) cat(paste( paste0(game, "_Retro="), mean(filter(d_long, Method == iv_name & Game == game)$Rating), paste0(game, "_ESM="), mean(filter(d_long, Method == dv_name & Game == game)$Rating), "p=", p_ints[i], "Cohen's d=", round(cohen_d, 3), "\n" )) } } print(paste( "MAD ", "p_main=", round(p_main[2], 3), "d=", round(res_sd$cohen.d[1], 3), res_sd$effect.d[1], iv_name, "=", round(mean(filter(d_sd, d_sd$Method == iv_name)$Rating), 3), dv_name, "=", round(mean(filter(d_sd, d_sd$Method == dv_name)$Rating), 3), "p_int=", round(p_int[2], 3))) if (p_int[2] <= .05) { ints <- summary(art.con(model_sd, "Method:Game")) # Adjust p-values and print contrasts p_ints <- p.adjust(c(ints$p.value[5], ints$p.value[14], ints$p.value[22], ints$p.value[29], ints$p.value[35])) for (i in 1:5) { game <- c("AC", "GotS", "HLA", "IEYTD", "RM")[i] cohen_d <- calculate_cohen_d(filter(d_long, Game == game), iv_name, dv_name) cat(paste( paste0(game, "_Retro="), mean(filter(d_long, Method == iv_name & Game == game)$Rating), paste0(game, "_ESM="), mean(filter(d_long, Method == dv_name & Game == game)$Rating), "p=", p_ints[i], "Cohen's d=", round(cohen_d, 3), "\n" )) } } print(paste( "Min", "p_main=", round(p_main[3], 3), "d=", round(res_min$cohen.d[1], 3), res_min$effect.d[1], iv_name, "=", round(mean(filter(d_min, d_min$Method == iv_name)$Rating), 3), dv_name, "=", round(mean(filter(d_min, d_min$Method == dv_name)$Rating), 3), "p_int=", round(p_int[3], 3))) if (p_int[3] <= .05) { ints <- summary(art.con(model_min, "Method:Game")) # Adjust p-values and print contrasts p_ints <- p.adjust(c(ints$p.value[5], ints$p.value[14], ints$p.value[22], ints$p.value[29], ints$p.value[35])) for (i in 1:5) { game <- c("AC", "GotS", "HLA", "IEYTD", "RM")[i] cohen_d <- calculate_cohen_d(filter(d_long, Game == game), iv_name, dv_name) cat(paste( paste0(game, "_Retro="), mean(filter(d_long, Method == iv_name & Game == game)$Rating), paste0(game, "_ESM="), mean(filter(d_long, Method == dv_name & Game == game)$Rating), "p=", p_ints[i], "Cohen's d=", round(cohen_d, 3), "\n" )) } } print(paste( "Max", "p_main=", round(p_main[4], 3), "d=", round(res_max$cohen.d[1], 3), res_max$effect.d[1], iv_name, "=", round(mean(filter(d_max, d_max$Method == iv_name)$Rating), 3), dv_name, "=", round(mean(filter(d_max, d_max$Method == dv_name)$Rating), 3), "p_int=", round(p_int[4], 3))) if (p_int[4] <= .05) { ints <- summary(art.con(model_max, "Method:Game")) # Adjust p-values and print contrasts p_ints <- p.adjust(c(ints$p.value[5], ints$p.value[14], ints$p.value[22], ints$p.value[29], ints$p.value[35])) for (i in 1:5) { game <- c("AC", "GotS", "HLA", "IEYTD", "RM")[i] cohen_d <- calculate_cohen_d(filter(d_long, Game == game), iv_name, dv_name) cat(paste( paste0(game, "_Retro="), mean(filter(d_long, Method == iv_name & Game == game)$Rating), paste0(game, "_ESM="), mean(filter(d_long, Method == dv_name & Game == game)$Rating), "p=", p_ints[i], "Cohen's d=", round(cohen_d, 3), "\n" )) } } } distr_tests("Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) distr_tests("Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) distr_tests("Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) distr_tests("Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) distr_tests("Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) distr_tests("Retro_Valence4", "ESM_Valence4", d_ESM) distr_tests("Retro_Arousal4", "ESM_Arousal4", d_ESM) ######## RQ2: How reliable is RetroSketch and ESM across different VR experiences and users? ######## run_covar_ordinal_regression <- function(formula, data, numOf_interaction_Terms){ localOR_to_cohens_D <-function(local_or){ # Convert log odds ratio to Cohen's d cohens_d <- oddsratio_to_d(local_or, log=TRUE) # Return Cohen's d return(cohens_d) } m <- polr(formula, data = data, Hess = TRUE) s1 <- summary(m) ctable <- coef(s1) p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2 (ctable <- cbind(ctable, "p value" = p)) print(ctable) if(numOf_interaction_Terms > 0){ adjusted_ps <- p.adjust(ctable[c(numOf_interaction_Terms+2):(numOf_interaction_Terms*2+1), "p value"], method = "holm") print(adjusted_ps) } } # Ensure data types are correct d_ESM_Main_Gender_Aggr$ESM_Joy_Value_quant = factor(d_ESM_Main_Gender_Aggr$ESM_Joy_Value_quant) d_ESM_Main_Gender_Aggr$ESM_Fear_Value_quant = factor(d_ESM_Main_Gender_Aggr$ESM_Fear_Value_quant) d_ESM_Main_Gender_Aggr$ESM_Relaxation_Value_quant = factor(d_ESM_Main_Gender_Aggr$ESM_Relaxation_Value_quant) d_ESM_Main_Gender_Aggr$ESM_Bored_Value_quant = factor(d_ESM_Main_Gender_Aggr$ESM_Bored_Value_quant) d_ESM_Main_Gender_Aggr$ESM_Presence_Value_quant = factor(d_ESM_Main_Gender_Aggr$ESM_Presence_Value_quant) d_ESM_Main_Gender_Aggr$ESM_Arousal4_quant = factor(d_ESM_Main_Gender_Aggr$ESM_Arousal4_quant) d_ESM_Main_Gender_Aggr$ESM_Valence4_quant = factor(d_ESM_Main_Gender_Aggr$ESM_Valence4_quant) run_covar_ordinal_regression(as.formula("ESM_Joy_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Age + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Gender + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * VR_Experience_Rating + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Qual7_Preference"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Fear_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Age + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Gender + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * VR_Experience_Rating + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Qual7_Preference"), d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Relaxation_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Age + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Gender + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * VR_Experience_Rating+ Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Qual7_Preference"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Bored_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Age + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Gender + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * VR_Experience_Rating+ Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Qual7_Preference"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Presence_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Age + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Gender + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * VR_Experience_Rating+ Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Qual7_Preference"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Arousal4_quant ~ Retro_Arousal4 * Age + Retro_Arousal4 * Gender + Retro_Arousal4 * VR_Experience_Rating+ Retro_Arousal4 * Qual7_Preference"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Valence4_quant ~ Retro_Valence4 * Age + Retro_Valence4 * Gender + Retro_Valence4 * VR_Experience_Rating+ Retro_Valence4 * Qual7_Preference"),d_ESM_Main_Gender_Aggr, 5) ############# PLAYER TYPE run_covar_ordinal_regression(as.formula("ESM_Joy_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Tondello_Action + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Tondello_Aesthetic + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Tondello_Immersion + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Tondello_Goal + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * Tondello_Social"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Fear_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Tondello_Action + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Tondello_Aesthetic + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Tondello_Immersion + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Tondello_Goal + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * Tondello_Social"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Relaxation_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Tondello_Action + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Tondello_Aesthetic + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Tondello_Immersion + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Tondello_Goal + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * Tondello_Social"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Bored_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Tondello_Action + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Tondello_Aesthetic + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Tondello_Immersion + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Tondello_Goal + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * Tondello_Social"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Presence_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Tondello_Action + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Tondello_Aesthetic + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Tondello_Immersion + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Tondello_Goal + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * Tondello_Social"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Arousal4_quant ~ Retro_Arousal4 * Tondello_Action + Retro_Arousal4 * Tondello_Aesthetic + Retro_Arousal4 * Tondello_Immersion + Retro_Arousal4 * Tondello_Goal + Retro_Arousal4 * Tondello_Social"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Valence4_quant ~ Retro_Valence4 * Tondello_Action + Retro_Valence4 * Tondello_Aesthetic + Retro_Valence4 * Tondello_Immersion + Retro_Valence4 * Tondello_Goal + Retro_Valence4 * Tondello_Social"),d_ESM_Main_Gender_Aggr, 5) ############# Personality Type run_covar_ordinal_regression(as.formula("ESM_Joy_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Joy_Rating * B5_Extroversion_Score + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * B5_Agreeableness_Score + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * B5_Conscientiousness_Score + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * B5_Neuroticism_Score + Exposure_Windowed_60s_Retro_Mean_Joy_Rating * B5_Openness_Score"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Fear_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Fear_Rating * B5_Extroversion_Score + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * B5_Agreeableness_Score + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * B5_Conscientiousness_Score + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * B5_Neuroticism_Score + Exposure_Windowed_60s_Retro_Mean_Fear_Rating * B5_Openness_Score"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Relaxation_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * B5_Extroversion_Score + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * B5_Agreeableness_Score + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * B5_Conscientiousness_Score + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * B5_Neuroticism_Score + Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * B5_Openness_Score"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Bored_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * B5_Extroversion_Score + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * B5_Agreeableness_Score + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * B5_Conscientiousness_Score + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * B5_Neuroticism_Score + Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * B5_Openness_Score"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Presence_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Presence_Rating * B5_Extroversion_Score + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * B5_Agreeableness_Score + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * B5_Conscientiousness_Score + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * B5_Neuroticism_Score + Exposure_Windowed_60s_Retro_Mean_Presence_Rating * B5_Openness_Score"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Arousal4_quant ~ Retro_Arousal4 * B5_Extroversion_Score + Retro_Arousal4 * B5_Agreeableness_Score + Retro_Arousal4 * B5_Conscientiousness_Score + Retro_Arousal4 * B5_Neuroticism_Score + Retro_Arousal4 * B5_Openness_Score"),d_ESM_Main_Gender_Aggr, 5) run_covar_ordinal_regression(as.formula("ESM_Valence4_quant ~ Retro_Valence4 * B5_Extroversion_Score + Retro_Valence4 * B5_Agreeableness_Score + Retro_Valence4 * B5_Conscientiousness_Score + Retro_Valence4 * B5_Neuroticism_Score + Retro_Valence4 * B5_Openness_Score"),d_ESM_Main_Gender_Aggr, 5) ############# Immersive Tendencies run_covar_ordinal_regression(as.formula("ESM_Joy_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Joy_Rating * ITQ_Total_Score"),d_ESM_Main_Gender_Aggr, 0) run_covar_ordinal_regression(as.formula("ESM_Fear_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Fear_Rating * ITQ_Total_Score"),d_ESM_Main_Gender_Aggr, 0) run_covar_ordinal_regression(as.formula("ESM_Relaxation_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating * ITQ_Total_Score"),d_ESM_Main_Gender_Aggr, 0) run_covar_ordinal_regression(as.formula("ESM_Bored_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Boredom_Rating * ITQ_Total_Score"),d_ESM_Main_Gender_Aggr, 0) run_covar_ordinal_regression(as.formula("ESM_Presence_Value_quant ~ Exposure_Windowed_60s_Retro_Mean_Presence_Rating * ITQ_Total_Score"),d_ESM_Main_Gender_Aggr, 0) run_covar_ordinal_regression(as.formula("ESM_Arousal4_quant ~ Retro_Arousal4 * ITQ_Total_Score"),d_ESM_Main_Gender_Aggr, 0) run_covar_ordinal_regression(as.formula("ESM_Valence4_quant ~ Retro_Valence4 * ITQ_Total_Score"),d_ESM_Main_Gender_Aggr, 0) # RQ2 influence of covariates on RetroSketch and ESM measures, considering interactions with Game run_covar_interactions <- function(formula, d1, global_sample=FALSE){ # Initialize an empty data frame to store interactions and reference levels interaction_results <- data.frame(BaseGame = character(), Interaction = character(), PValue = numeric(), stringsAsFactors = FALSE) # First analysis if(global_sample){ s1 <- run_ordinal_regression_using_global_sample(formula, d1) } else{ s1 <- run_ordinal_regression(formula, d1) } interactions <- s1$coefficients[c(16:19), 4] interaction_names <- rownames(s1$coefficients)[c(16:19)] base_game <- "AC" # Store the results interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Second analysis d1$Game <- relevel(d1$Game, ref = "GotS") if(global_sample){ s1 <- run_ordinal_regression_using_global_sample(formula, d1) } else{ s1 <- run_ordinal_regression(formula, d1) } interactions <- s1$coefficients[c(17:19), 4] interaction_names <- rownames(s1$coefficients)[c(17:19)] base_game <- "GotS" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Third analysis d1$Game <- relevel(d1$Game, ref = "HLA") if(global_sample){ s1 <- run_ordinal_regression_using_global_sample(formula, d1) } else{ s1 <- run_ordinal_regression(formula, d1) } interactions <- s1$coefficients[c(18:19), 4] interaction_names <- rownames(s1$coefficients)[c(18:19)] base_game <- "HLA" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Fourth analysis d1$Game <- relevel(d1$Game, ref = "IEYTD") if(global_sample){ s1 <- run_ordinal_regression_using_global_sample(formula, d1) } else{ s1 <- run_ordinal_regression(formula, d1) } interactions <- s1$coefficients[c(19), 4] interaction_names <- rownames(s1$coefficients)[c(19)] base_game <- "IEYTD" # Store the results interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions)) # Adjust p-values interaction_results$AdjustedPValue <- p.adjust(interaction_results$PValue, method = "holm") # Print the results print("ORDINAL REGRESSION INTERACTIONS: ") rownames(interaction_results) <- NULL print(interaction_results, row.names = FALSE) } covar_corr <- function(dv_name, iv_name, d1, regtest = TRUE, global_sample = FALSE) { #Interaction tests if(regtest){ if(dv_name == "ESM_Arousal4"){ formula <- as.formula(paste("ESM_Arousal4_quant", "~", iv_name, "* Game")) } else if(dv_name == "ESM_Valence4") { formula <- as.formula(paste("ESM_Valence4_quant", "~", iv_name, "* Game")) } else{ formula <- as.formula(paste(dv_name, "~", iv_name, "* Game")) } run_covar_interactions(formula, d1, global_sample = global_sample) } # Aggregate data for each variable d2 = aggregate(as.formula(paste(dv_name, "~ Participant_ID")), data = d1, FUN = mean) d3 = aggregate(as.formula(paste(iv_name, "~ Participant_ID")), data = d1, FUN = mean) # Merge datasets on Participant_ID to ensure matching entries merged_data = merge(d2, d3, by = "Participant_ID") # Check lengths of merged vectors if (nrow(merged_data) == 0) { cat("Warning: No matching data found for", dv_name, "and", iv_name, "\n") return(NULL) # Exit the function if no matching data } # Extract the variables of interest dv_values = merged_data[, dv_name] iv_values = merged_data[, iv_name] # Check for missing values missing_dv = sum(is.na(dv_values)) missing_iv = sum(is.na(iv_values)) if (missing_dv > 0 | missing_iv > 0) { cat("Warning: There are missing values in the data.\n") cat("Missing values in", dv_name, ":", missing_dv, "\n") cat("Missing values in", iv_name, ":", missing_iv, "\n") } # Perform the correlation test p_ct = cor.test(dv_values, iv_values, method="pearson") k_ct = cor.test(dv_values, iv_values, method="kendall") cat("Running Correlation overall: ", dv_name, "~", iv_name, "\n") print(p_ct) print(k_ct) print(interpret_kendalls_tau(cor(dv_values, iv_values, method="kendall"))) # Plot the results plot(dv_values, iv_values, xlab=dv_name, ylab=iv_name, main=paste("Scatter plot of", dv_name, "vs", iv_name)) } run_Retro_ESM_covar_corr <-function(covar) { d_retro_filtered <- dl %>% filter(Retro_Sample_Frame %in% seq(300, 1800, by = 300)) covar_corr("ESM_Joy_Value", covar, d_ESM, regtest=FALSE) #covar_corr("Exposure_Windowed_60s_Retro_Mean_Joy_Rating_quant", covar, d_ESM) covar_corr("Exposure_Windowed_60s_Retro_Mean_Joy_Rating_quant", covar, d_retro_filtered, regtest=FALSE, global_sample = TRUE) print("----------------------") covar_corr("ESM_Fear_Value", covar, d_ESM, regtest=FALSE) covar_corr("Exposure_Windowed_60s_Retro_Mean_Fear_Rating_quant", covar, d_retro_filtered, regtest=FALSE, global_sample = TRUE) print("----------------------") covar_corr("ESM_Relaxation_Value", covar, d_ESM,regtest=FALSE) covar_corr("Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating_quant", covar, d_retro_filtered, regtest=FALSE, global_sample = TRUE) print("----------------------") covar_corr("ESM_Bored_Value", covar, d_ESM, regtest=FALSE) covar_corr("Exposure_Windowed_60s_Retro_Mean_Boredom_Rating_quant", covar, d_retro_filtered, regtest=FALSE, global_sample = TRUE) print("----------------------") covar_corr("ESM_Presence_Value", covar, d_ESM, regtest=FALSE) covar_corr("Exposure_Windowed_60s_Retro_Mean_Presence_Rating_quant", covar, d_retro_filtered, regtest=FALSE, global_sample = TRUE) print("----------------------") covar_corr("ESM_Arousal4", covar, d_ESM, regtest=FALSE) covar_corr("Retro_Arousal4_quant", covar, d_retro_filtered, regtest=FALSE, global_sample = TRUE) print("----------------------") covar_corr("ESM_Valence4", covar, d_ESM, regtest=FALSE) covar_corr("Retro_Valence4_quant", covar, d_retro_filtered, regtest=FALSE, global_sample = TRUE) print("----------------------") } run_Retro_ESM_covar_corr("Age") run_Retro_ESM_covar_corr("VR_Experience_Rating") run_Retro_ESM_covar_corr("Tondello_Action") run_Retro_ESM_covar_corr("Tondello_Aesthetic") run_Retro_ESM_covar_corr("Tondello_Immersion") run_Retro_ESM_covar_corr("Tondello_Goal") run_Retro_ESM_covar_corr("Tondello_Social") run_Retro_ESM_covar_corr("B5_Extroversion_Score") run_Retro_ESM_covar_corr("B5_Agreeableness_Score") run_Retro_ESM_covar_corr("B5_Conscientiousness_Score") run_Retro_ESM_covar_corr("B5_Neuroticism_Score") run_Retro_ESM_covar_corr("B5_Openness_Score") run_Retro_ESM_covar_corr("ITQ_Total_Score") run_art_gender <- function(dv1, dv2, data){ d_iv = melt(setDT(data), id.vars=c("Participant_ID"), measure.vars=c(dv1, dv2), variable.name="Method", value.name="Rating") d_iv <- left_join(d_iv, data, by = "Participant_ID") d_iv <- filter(d_iv, d_iv$Gender == "Male" | d_iv$Gender == "Female") d_iv$Gender <- factor(d_iv$Gender) d_iv$Method <- factor(d_iv$Method) d_iv$Game <- factor(d_iv$Game) model = art(Rating ~ Gender * Game * Method, data = d_iv) # use all data for better power res = anova(model) res$part.eta.sq = res$F * res$Df / (res$F * res$Df + res$Df.res) res$effect = interpret_omega_squared(res$part.eta.sq, rules = "cohen1992") res # Perform contrasts for significant interactions contrast_results <- list() if (res$`Pr(>F)`[which(rownames(res) == "Gender:Game")] < 0.05) { contrast_results[["Gender:Game"]] <- summary(art.con(model, "Gender:Game")) } if (res$`Pr(>F)`[which(rownames(res) == "Gender:Method")] < 0.05) { contrast_results[["Gender:Method"]] <- summary(art.con(model, "Gender:Method")) } if (res$`Pr(>F)`[which(rownames(res) == "Game:Method")] < 0.05) { contrast_results[["Game:Method"]] <- summary(art.con(model, "Game:Method")) } if (res$`Pr(>F)`[which(rownames(res) == "Gender:Game:Method")] < 0.05) { contrast_results[["Gender:Game:Method"]] <- summary(art.con(model, "Gender:Game:Method")) } return(list(ANOVA = res, Contrasts = contrast_results)) } run_art_gender("ESM_Joy_Value", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", d_ESM_Aggregated) run_art_gender("ESM_Fear_Value", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", d_ESM_Aggregated) run_art_gender("ESM_Relaxation_Value", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", d_ESM_Aggregated) run_art_gender("ESM_Bored_Value", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", d_ESM_Aggregated) run_art_gender("ESM_Presence_Value", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", d_ESM_Aggregated) run_art_gender("ESM_Arousal4", "Retro_Arousal4", d_ESM_Aggregated) run_art_gender("ESM_Valence4", "Retro_Valence4", d_ESM_Aggregated) run_art_preference <- function(dv1, dv2, data){ d_iv = melt(setDT(data), id.vars=c("Participant_ID"), measure.vars=c(dv1, dv2), variable.name="Method", value.name="Rating") d_iv <- left_join(d_iv, data, by = "Participant_ID") d_iv$Qual7_Preference <- factor(d_iv$Qual7_Preference) d_iv$Method <- factor(d_iv$Method) d_iv$Game <- factor(d_iv$Game) model = art(Rating ~ Qual7_Preference * Game * Method, data = d_iv) # use all data for better power res = anova(model) res$part.eta.sq = res$F * res$Df / (res$F * res$Df + res$Df.res) res$effect = interpret_omega_squared(res$part.eta.sq, rules = "cohen1992") print(res) print(summary(art.con(model, "Qual7_Preference"))) # Perform contrasts for significant interactions contrast_results <- list() if (res$`Pr(>F)`[which(rownames(res) == "Qual7_Preference:Game")] < 0.05) { contrast_results[["Qual7_Preference:Game"]] <- summary(art.con(model, "Qual7_Preference:Game")) } if (res$`Pr(>F)`[which(rownames(res) == "Qual7_Preference:Method")] < 0.05) { contrast_results[["Qual7_Preference:Method"]] <- summary(art.con(model, "Qual7_Preference:Method")) } if (res$`Pr(>F)`[which(rownames(res) == "Qual7_Preference:Method")] < 0.05) { contrast_results[["Qual7_Preference:Method"]] <- summary(art.con(model, "Qual7_Preference:Method")) } if (res$`Pr(>F)`[which(rownames(res) == "Qual7_Preference:Game:Method")] < 0.05) { contrast_results[["Qual7_Preference:Game:Method"]] <- summary(art.con(model, "Qual7_Preference:Game:Method")) } return(list(ANOVA = res, Contrasts = contrast_results)) } run_art_preference("ESM_Joy_Value", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", d_ESM_Aggregated) run_art_preference("ESM_Fear_Value", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", d_ESM_Aggregated) run_art_preference("ESM_Relaxation_Value", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", d_ESM_Aggregated) run_art_preference("ESM_Bored_Value", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", d_ESM_Aggregated) run_art_preference("ESM_Presence_Value", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", d_ESM_Aggregated) run_art_preference("ESM_Arousal4", "Retro_Arousal4", d_ESM_Aggregated) run_art_preference("ESM_Valence4", "Retro_Valence4", d_ESM_Aggregated) ######## RQ3: How does ESM influence the VR user experience? ######## # RQ3: Compare post-session questionnaire scores diff_within <- function(dv_name, d, d_A, d_B) { # Filter the data d1_NoESM <- filter(d, d$Condition == "Retro") d1_ESM <- filter(d, d$Condition == "ESM_Retro") # Shapiro-Wilk test for normality sw_A <- shapiro.test(d1_NoESM[, dv_name]) sw_B <- shapiro.test(d1_ESM[, dv_name]) # Print the Shapiro-Wilk test result cat("Shapiro-Wilk Test for Normality:\n") print(sw_A) print(sw_B) # ART ANOVA if data is not normally distributed model <- art(as.formula(paste(dv_name, "~ Condition * Game + (1|Participant_ID)")), data = d) res <- anova(model) res$part.eta.sq <- res$F * res$Df / (res$F * res$Df + res$Df.res) res$cohen.d = sqrt(4 * res$part.eta.sq / (1 - res$part.eta.sq )) res$effect = interpret_omega_squared(res$part.eta.sq, rules = "cohen1992") res$effect.d = interpret_cohens_d(res$cohen.d, rules = "cohen1988") print(res) if (res$`Pr(>F)`[3] <= .05) { ints <- summary(art.con(model, "Condition:Game")) # Adjust p-values p_ints <- p.adjust(c(ints$p.value[5], ints$p.value[14], ints$p.value[22], ints$p.value[29], ints$p.value[35])) # Calculate means for contrasts means <- data.frame( Game = c("AC", "GotS", "HLA", "IEYTD", "RM"), NoESM = sapply(c("AC", "GotS", "HLA", "IEYTD", "RM"), function(game) mean(filter(d1_NoESM, Game == game)[[dv_name]])), ESM = sapply(c("AC", "GotS", "HLA", "IEYTD", "RM"), function(game) mean(filter(d1_ESM, Game == game)[[dv_name]])), p_value = p_ints ) # Print contrasts with p-values and effect sizes cat("Contrasts:\n") for (i in 1:nrow(means)) { # Calculate effect sizes sd_NoESM <- sd(filter(d1_NoESM, Game == means$Game[i])[[dv_name]]) sd_ESM <- sd(filter(d1_ESM, Game == means$Game[i])[[dv_name]]) pooled_sd <- sqrt((sd_NoESM^2 + sd_ESM^2) / 2) effect_size <- (means$ESM[i] - means$NoESM[i]) / pooled_sd # Print with effect sizes cat(sprintf("%s_delta= %.3f \ p=%.3f \t d=%.3f\n", means$Game[i], (means$NoESM[i] - means$ESM[i]), means$p_value[i], effect_size)) } } cat(paste("\n", dv_name, " MEAN Delta =", round((mean(filter(d, d$Condition == "Retro")[, dv_name]) - mean(filter(d, d$Condition == "ESM_Retro"))), 3))) } # Scale psychometric variables so scores are expressed in the range of the item scale d1 = dl d1$IMI_Competence_Score_AVG = d1$IMI_Competence_Score / 6 d1$IMI_Interest_Score_AVG = d1$IMI_Interest_Score / 7 d1$IMI_Pressure_Score_AVG = d1$IMI_Pressure_Score / 5 d1$IMI_Effort_Score_AVG = d1$IMI_Effort_Score / 5 d1$PPL_FSQ_Absorption_Score_AVG= d1$PPL_FSQ_Absorption_Score / 9 d1$PPL_FSQ_Challenge_Skill_Score_AVG = d1$PPL_FSQ_Challenge_Skill_Score / 11 d1$MPS_Phys_Presence_Score_AVG= d1$MPS_Phys_Presence_Score / 5 d1$MPS_Self_Presence_Score_AVG = d1$MPS_Self_Presence_Score / 5 d1$MPS_Social_Presence_Score_AVG = d1$MPS_Social_Presence_Score / 5 diff_within("IMI_Competence_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("IMI_Interest_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("IMI_Pressure_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("IMI_Effort_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("MPS_Phys_Presence_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("MPS_Self_Presence_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("MPS_Social_Presence_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("PPL_FSQ_Absorption_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("PPL_FSQ_Challenge_Skill_Score_AVG", d1, d_Retro, d_ESM_Retro) diff_within("SSQ_Nausea", d1, d_Retro, d_ESM_Retro) diff_within("SSQ_Occulomotor_Disturbance", d1, d_Retro, d_ESM_Retro) diff_within("SSQ_Disorientation", d1, d_Retro, d_ESM_Retro) diff_within("SSQ_Total", d1, d_Retro, d_ESM_Retro) diff_within("RPE_BORG_CR100", d1, d_Retro, d_ESM_Retro) diff_within("FSS_Absorption_Score", filter(d1, !is.na(d1$FSS_Absorption_Score)), d_Retro, d_ESM_Retro) diff_within("FSS_Demand_Score", filter(d1, !is.na(d1$FSS_Demand_Score)), d_Retro, d_ESM_Retro) diff_within("FSS_Performance_Fluency_Score", filter(d1, !is.na(d1$FSS_Performance_Fluency_Score)), d_Retro, d_ESM_Retro) diff_within("FSS_Total_Score", filter(d1, !is.na(d1$FSS_Total_Score)), d_Retro, d_ESM_Retro) # RQ3: Test effects of ESM by comparing RetroSketch scores between within-participant conditions # Consider data points around ESM measures, 3 minutes from ESM start d1 = filter(dl, dl$Retro_Sample_Frame %% 300 >= 0 & dl$Retro_Sample_Frame %% 300 < 180 & dl$Retro_Sample_Frame != 0) diff_within("Retro_Sample_Joy_Rating", d1, d_Retro, d_ESM_Retro) diff_within("Retro_Sample_Fear_Rating", d1, d_Retro, d_ESM_Retro) diff_within("Retro_Sample_Relaxation_Rating", d1, d_Retro, d_ESM_Retro) diff_within("Retro_Sample_Boredom_Rating", d1, d_Retro, d_ESM_Retro) diff_within("Retro_Sample_Presence_Rating", d1, d_Retro, d_ESM_Retro) diff_within("Retro_Valence4", d1, d_Retro, d_ESM_Retro) diff_within("Retro_Arousal4", d1, d_Retro, d_ESM_Retro) ######## RQ4: How do RetroSketch and ESM relate to physiological measures? ######## # Use ordinal regressions to test interactions with Game run_interactions <- function(formula, d1, buckets){ localOR_to_cohens_D <-function(local_or){ # Convert log odds ratio to Cohen's d cohens_d <- oddsratio_to_d(local_or, log=TRUE) # Return Cohen's d return(cohens_d) } # Initialize an empty data frame to store interactions and reference levels interaction_results <- data.frame(BaseGame = character(), Interaction = character(), estimate=numeric(), PValue = numeric(), cohens.d = numeric(), stringsAsFactors = FALSE) # First analysis s1 <- run_ordinal_regression(formula, d1) print(s1) interactions <- s1$coefficients[c((buckets+5):(buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+5):(buckets+8))] estimates <- s1$coefficients[c((buckets+5):(buckets+8)), 1] base_game <- "AC" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Second analysis d1$Game <- relevel(d1$Game, ref = "GotS") s1 <- run_ordinal_regression(formula, d1) interactions <- s1$coefficients[c((buckets+6):(buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+6):(buckets+8))] estimates <- s1$coefficients[c((buckets+6):(buckets+8)), 1] base_game <- "GotS" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Third analysis d1$Game <- relevel(d1$Game, ref = "HLA") s1 <- run_ordinal_regression(formula, d1) interactions <- s1$coefficients[c((buckets+7):(buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+7):(buckets+8))] estimates <- s1$coefficients[c((buckets+7):(buckets+8)), 1] base_game <- "HLA" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Fourth analysis d1$Game <- relevel(d1$Game, ref = "IEYTD") s1 <- run_ordinal_regression(formula, d1) interactions <- s1$coefficients[c((buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+8))] estimates <- s1$coefficients[c((buckets+8)), 1] base_game <- "IEYTD" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Adjust p-values interaction_results$AdjustedPValue <- p.adjust(interaction_results$PValue, method = "holm") interaction_results$cohens.d <- sapply(interaction_results$estimate, localOR_to_cohens_D) return(interaction_results) } # Create linear models and compare their fit using encompassing test corr_regression_comparison <- function (dv_name, iv1_name, iv2_name, d1, regtest = TRUE, buckets=30) { # workaround to fix scoping problem dv_name <<- dv_name iv1_name <<- iv1_name iv2_name <<- iv2_name d1 <<- d1 model_retro = gls(as.formula(paste("scale(",dv_name, ")~scale(", iv1_name, ")")), data = d1, na.action=na.exclude, corr = DefaultCorr) stdpar_retro = standardize_parameters(model_retro) r2_retro = R2_lik(model_retro) summary_retro = summary(model_retro) if (regtest) { # Use ordinal GEE regression to test main effect of IV and interactions with Game # Quantize physiological measure to enable use of ordinal regression qdata <- quantize(data=d1, expnms=c(dv_name), q=buckets) d1$dv_quant <- qdata$data[[dv_name]] formula = as.formula(paste("dv_quant", "~", iv1_name, "* Game")) if(grepl("Smile", dv_name, fixed = TRUE)){ interaction_results <- run_interactions(formula, d1, buckets = buckets-10) } else if(grepl("O_Shape", dv_name, fixed = TRUE)){ interaction_results <- run_interactions(formula, d1, buckets = buckets-5) } else{ interaction_results <- run_interactions(formula, d1, buckets) } # Print the results print(paste("ORDINAL REGRESSION INTERACTIONS FOR ", iv1_name, ":")) rownames(interaction_results) <- NULL print(interaction_results, row.names = FALSE) } model_esm = gls(as.formula(paste("scale(",dv_name, ")~", "scale(",iv2_name,")")), data = d1, na.action=na.exclude, corr = DefaultCorr) r2_esm = R2_lik(model_esm) stdpar_esm = standardize_parameters(model_esm) summary_esm = summary(model_esm) if (regtest) { # Use ordinal GEE regression to test main effect of IV and interactions with Game qdata <- quantize(data=d1, expnms=c(dv_name), q=buckets) d1$dv_quant <- qdata$data[[dv_name]] formula = as.formula(paste("dv_quant", "~", iv2_name, "* Game")) if(grepl("Smile", dv_name, fixed = TRUE)){ interaction_results <- run_interactions(formula, d1, buckets = buckets-10) } else if(grepl("O_Shape", dv_name, fixed = TRUE)){ interaction_results <- run_interactions(formula, d1, buckets = buckets-5) } else{ interaction_results <- run_interactions(formula, d1, buckets) } # Print the results print(paste("ORDINAL REGRESSION INTERACTIONS FOR ", iv2_name, ":")) rownames(interaction_results) <- NULL print(interaction_results, row.names = FALSE) } model_enc = gls(as.formula(paste("scale(", dv_name, ") ~"," scale(" , iv1_name, ") + scale(", iv2_name, ")")), data = d1, na.action=na.exclude, corr = DefaultCorr) r2_enc = R2_lik(model_enc) f2_retro = (r2_retro - r2_enc) / (1 - r2_enc) f2_esm = (r2_esm - r2_enc) / (1 - r2_enc) cohen_d_retro = 2*sqrt(abs(as.numeric(f2_retro))) cohen_d_esm = 2*sqrt(abs(as.numeric(f2_esm))) print(cohen_d_retro) print(cohen_d_esm) cohen_d_delta = cohen_d_retro - cohen_d_esm enc = encomptest(model_retro, model_esm) print(enc) cat("DV Method Coeff R2 r AIC p_encompass f2 d \n") cat(paste0( dv_name, " Retro ", round(summary_retro$coefficients[2], 3), " ", round(r2_retro, 3), " ", round(sqrt(r2_retro), 3), " ", round(summary_retro$AIC, 3), " ", print_p(enc$`Pr(>Chisq)`[1]), " ", round(f2_retro, 3), " ",round(cohen_d_retro, 3), "\n")) cat(paste0( dv_name, " ESM ", round(summary_esm$coefficients[2], 3), " ", round(r2_esm, 3), " ", round(sqrt(r2_esm), 3), " ", round(summary_esm$AIC, 3), " ", print_p(enc$`Pr(>Chisq)`[2]), " ", round(f2_esm, 3), " ", round(cohen_d_esm, 3), "\n")) cat("Delta F^2= ", round((f2_retro - f2_esm), 3), "\n", "Delta Cohen's d= ", round(cohen_d_delta, 3), "\n") } ########## Fear corr_regression_comparison("Exposure_Windowed_60s_Mean_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) # Randomly downsample data to make some ordinal regressions computationally tractable set.seed(123) d_ESM_downsampled <- d_ESM[sample(1:nrow(d_ESM), size = 0.9 * nrow(d_ESM)), ] corr_regression_comparison("Exposure_Windowed_60s_Mean_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM_downsampled) corr_regression_comparison("Exposure_Windowed_60s_SD_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_HearRate_BPM", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_RR_Interval_CLEANED_ABS_RMSSD", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_InterBlinkInterval", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_blinkDuration", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Smile", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_5_Mouth_O_Shape", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_ESM) ########## Joy corr_regression_comparison("Exposure_Windowed_60s_Mean_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_HearRate_BPM", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_RR_Interval_CLEANED_ABS_RMSSD", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_InterBlinkInterval", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_blinkDuration", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Smile", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_5_Mouth_O_Shape", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_ESM) ## Boredom corr_regression_comparison("Exposure_Windowed_60s_Mean_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_HearRate_BPM", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_RR_Interval_CLEANED_ABS_RMSSD", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_InterBlinkInterval", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_blinkDuration", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Smile", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_5_Mouth_O_Shape", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_ESM) ## Relaxation corr_regression_comparison("Exposure_Windowed_60s_Mean_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_HearRate_BPM", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_RR_Interval_CLEANED_ABS_RMSSD", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_InterBlinkInterval", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_blinkDuration", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Smile", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_5_Mouth_O_Shape", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_ESM) ## Presence corr_regression_comparison("Exposure_Windowed_60s_Mean_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Pupil_Dilation", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_HearRate_BPM", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_RR_Interval_CLEANED_ABS_RMSSD", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_InterBlinkInterval", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_blinkDuration", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Smile", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_5_Mouth_O_Shape", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_ESM) ## Valence corr_regression_comparison("Exposure_Windowed_60s_Mean_Pupil_Dilation", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Pupil_Dilation", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_HearRate_BPM", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_RR_Interval_CLEANED_ABS_RMSSD", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_InterBlinkInterval", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_blinkDuration", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Smile", "Retro_Valence4", "ESM_Valence4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_5_Mouth_O_Shape", "Retro_Valence4", "ESM_Valence4", d_ESM) ## Arousal corr_regression_comparison("Exposure_Windowed_60s_Mean_Pupil_Dilation", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Pupil_Dilation", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_SD_Skin_Conductance_uS_CLEANED_ABS_Threshold", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_HearRate_BPM", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_RR_Interval_CLEANED_ABS_RMSSD", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_InterBlinkInterval", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_blinkDuration", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_Smile", "Retro_Arousal4", "ESM_Arousal4", d_ESM) corr_regression_comparison("Exposure_Windowed_60s_Mean_5_Mouth_O_Shape", "Retro_Arousal4", "ESM_Arousal4", d_ESM) ######## RQ5: How do RetroSketch and ESM relate to qualitative measures? ######## d_keypoint = filter(dl, !is.na(dl$Exposure_Windowed_60s_Retro_Keypoint_Mean_Sentiment_Positive)) # Calculate overall sentiment by combining positive, negative and neutral sentiment sub-scores d_keypoint$sentiment = (d_keypoint$Exposure_Windowed_60s_Retro_Keypoint_Mean_Sentiment_Positive - d_keypoint$Exposure_Windowed_60s_Retro_Keypoint_Mean_Sentiment_Negative) / (d_keypoint$Exposure_Windowed_60s_Retro_Keypoint_Mean_Sentiment_Positive + d_keypoint$Exposure_Windowed_60s_Retro_Keypoint_Mean_Sentiment_Negative + d_keypoint$Exposure_Windowed_60s_Retro_Keypoint_Mean_Sentiment_Neutral) d_keypoint_ESM = filter(d_keypoint, !is.na(d_keypoint$ESM_Fear_Value)) d_keypoint_filtered <- d_keypoint %>% filter(Retro_Sample_Frame %in% seq(300, 1800, by = 300)) run_sentiment_ordinal <-function(iv_name, formula, d1){ if(grepl("ESM", iv_name, fixed = TRUE)){ return(run_ordinal_regression(formula, d1)) } else{ return(run_ordinal_regression_using_global_sample(formula, d1)) } } # Test correlations and regressions with sentiment scores, including their interactions with Game. corr_regression_sentiment <- function(dv_name, iv_name, d1, regtest = TRUE, buckets=30) { localOR_to_cohens_D <-function(local_or){ # Convert log odds ratio to Cohen's d cohens_d <- oddsratio_to_d(local_or, log=TRUE) # Return Cohen's d return(cohens_d) } # Workaround to fix scoping problem f1 <<- as.formula(paste(dv_name, "~", iv_name)) d1 <<- d1 iv_name <<- iv_name pearson_list <- list() kendall_list <- list() if(regtest){ # Quantize the data qdata <- quantize(data=d1, expnms=c("sentiment"), q=buckets) d1$sentiment_quant <- qdata$data$sentiment formula <- as.formula(paste("sentiment_quant", "~", iv_name, "* Game")) # Initialize an empty data frame to store interactions and reference levels interaction_results <- data.frame(BaseGame = character(), Interaction = character(), PValue = numeric(), stringsAsFactors = FALSE) # First analysis s1 <- run_sentiment_ordinal(iv_name, formula, d1) interactions <- s1$coefficients[c((buckets+5):(buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+5):(buckets+8))] estimates <- s1$coefficients[c((buckets+5):(buckets+8)), 1] base_game <- "AC" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Second analysis d1$Game <- relevel(d1$Game, ref = "GotS") s1 <- run_sentiment_ordinal(iv_name, formula, d1) interactions <- s1$coefficients[c((buckets+6):(buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+6):(buckets+8))] base_game <- "GotS" estimates <- s1$coefficients[c((buckets+6):(buckets+8)), 1] interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Third analysis d1$Game <- relevel(d1$Game, ref = "HLA") s1 <- run_sentiment_ordinal(iv_name, formula, d1) interactions <- s1$coefficients[c((buckets+7):(buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+7):(buckets+8))] estimates <- s1$coefficients[c((buckets+7):(buckets+8)), 1] base_game <- "HLA" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Fourth analysis d1$Game <- relevel(d1$Game, ref = "IEYTD") s1 <- run_sentiment_ordinal(iv_name, formula, d1) interactions <- s1$coefficients[c((buckets+8)), 4] interaction_names <- rownames(s1$coefficients)[c((buckets+8))] estimates <- s1$coefficients[c((buckets+8)), 1] base_game <- "IEYTD" interaction_results <- rbind(interaction_results, data.frame(BaseGame = base_game, Interaction = interaction_names, PValue = interactions, estimate = estimates)) # Adjust p-values interaction_results$AdjustedPValue <- p.adjust(interaction_results$PValue, method = "holm") # Calculate eff size interaction_results$cohens.d <- sapply(interaction_results$estimate, localOR_to_cohens_D) # Print the results print("ORDINAL REGRESSION INTERACTIONS: ") rownames(interaction_results) <- NULL print(interaction_results, row.names = FALSE) } # Correlation overall d2 = aggregate(as.formula(paste(dv_name, "~ Participant_ID")), data = d1, FUN = mean) d3 = aggregate(as.formula(paste(iv_name, "~ Participant_ID")), data = d1, FUN = mean) p_ct = cor.test(d2[, dv_name], d3[, iv_name], method="pearson") k_ct = cor.test(d2[, dv_name], d3[, iv_name], method="kendall") print("Running Correlation overall") print(p_ct) print(interpret_r(p_ct$estimate)) # https://easystats.github.io/effectsize/articles/interpret.html print(k_ct) print(interpret_kendalls_tau(cor(d2[, dv_name], d3[, iv_name], method="kendall"))) ## Test assumptions for Pearson correlation tests qqnorm(d2[, dv_name]) qqnorm(d3[, iv_name]) print(shapiro.test(d2[, dv_name])) print(shapiro.test(d3[, iv_name])) plot(d2[, dv_name], d3[, iv_name]) # Correlation for ESM Intervals for (interval in seq(300, 1800, by = 300)) { result <- run_correlations_with_interval(d1, interval, dv_name, iv_name) # Store the result in the appropriate list based on the method if (result$method == "pearson") { pearson_list[[length(pearson_list) + 1]] <- result$value } else if (result$method == "kendall") { kendall_list[[length(kendall_list) + 1]] <- result$value } } # Convert lists to numeric vectors pearson_list <- unlist(pearson_list) kendall_list <- unlist(kendall_list) # Print the minimum value from each list if (length(pearson_list) > 0) { cat("Minimum Pearson's r:", min(pearson_list), "\n") print(interpret_r(min(pearson_list))) } else { cat("No Pearson's correlations calculated.\n") } if (length(kendall_list) > 0) { cat("Minimum Kendall's tau:", min(kendall_list), "\n") print(interpret_kendalls_tau(min(kendall_list))) } else { cat("No Kendall's correlations calculated.\n") } print(paste0(dv_name, " ")) if (regtest) { print(paste0(print_p(s1$coefficients[buckets, 4]))) } } corr_regression_sentiment("sentiment", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", d_keypoint_filtered, buckets=11) corr_regression_sentiment("sentiment", "ESM_Joy_Value", d_keypoint_ESM, buckets=11) corr_regression_comparison("sentiment", "Exposure_Windowed_60s_Retro_Mean_Joy_Rating", "ESM_Joy_Value", d_keypoint_ESM, regtest = FALSE) corr_regression_sentiment("sentiment", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", d_keypoint_filtered, buckets=11) corr_regression_sentiment("sentiment", "ESM_Fear_Value_quant", d_keypoint_ESM, buckets=11) corr_regression_comparison("sentiment", "Exposure_Windowed_60s_Retro_Mean_Fear_Rating", "ESM_Fear_Value", d_keypoint_ESM, regtest = FALSE) corr_regression_sentiment("sentiment", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", d_keypoint_filtered, buckets=11) corr_regression_sentiment("sentiment", "ESM_Relaxation_Value", d_keypoint_ESM, buckets=11) corr_regression_comparison("sentiment", "Exposure_Windowed_60s_Retro_Mean_Relaxation_Rating", "ESM_Relaxation_Value", d_keypoint_ESM, regtest = FALSE) corr_regression_sentiment("sentiment", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", d_keypoint_filtered, buckets=11) corr_regression_sentiment("sentiment", "ESM_Bored_Value", d_keypoint_ESM, buckets=11) corr_regression_comparison("sentiment", "Exposure_Windowed_60s_Retro_Mean_Boredom_Rating", "ESM_Bored_Value", d_keypoint_ESM, regtest = FALSE) corr_regression_sentiment("sentiment", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", d_keypoint_filtered, buckets=11) corr_regression_sentiment("sentiment", "ESM_Presence_Value", d_keypoint_ESM, buckets=11) corr_regression_comparison("sentiment", "Exposure_Windowed_60s_Retro_Mean_Presence_Rating", "ESM_Presence_Value", d_keypoint_ESM, regtest = FALSE) corr_regression_sentiment("sentiment", "Retro_Valence4", d_keypoint_filtered, buckets=11) corr_regression_sentiment("sentiment", "ESM_Valence4", d_keypoint_ESM, buckets=11) corr_regression_comparison("sentiment", "Retro_Valence4", "ESM_Valence4", d_keypoint_ESM, regtest = FALSE) corr_regression_sentiment("sentiment", "Retro_Arousal4", d_keypoint_filtered, buckets=11) corr_regression_sentiment("sentiment", "ESM_Arousal4", d_keypoint_ESM, buckets=11) corr_regression_comparison("sentiment", "Retro_Arousal4", "ESM_Arousal4", d_keypoint_ESM, regtest = FALSE)