### R script 2 - Main analysis # Attractive Things Do Work Better: A Meta-Analysis on Visual Aesthetics and User Performance ### install and load packages # install.packages("rstudioapi") library(rstudioapi) # set working directory to active folder # install.packages("readxl") library(readxl) # read .xlsx datasheet # install.packages("dplyr") library(dplyr) # data manipulation and transformation (i.e., recoding variables, aggregating variables) # install.packages("purrr") library(purrr) # use map_dfr to return data frame from applied function # install.packages("metafor") library(metafor) # conduct meta-analysis, create forest plot # install.packages("clubSandwich") library(clubSandwich) # apply robust variance estimation (RVE) to meta-analytic model # install.packages("DescTools") library(DescTools) # calculate Cramer's V # install.packages("ggplot2") library(ggplot2) # visualize correlation matrix library(reshape2) # reorganize correlation matrix for visualization # install.packages("misty") library(misty) # conduct Grand Mean Centering for continuous moderators # install.packages("RColorBrewer") library(RColorBrewer) # color studies in funnel plot ### Set working directory and load data setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) Data_AestheticsPerformance_noNA <- read_excel("Data_AestheticsPerformance_noNA.xlsx", col_types = c("numeric", # study_id "numeric", # es_id "text", # authors "text", # title "text", # pub_status "text", # pub_type "numeric", # year "text", "text", "text", # country_author, country_study, continent_study_recoded "numeric", "numeric", # total_n, women_n "numeric", "numeric", "text", # age_mean, age_sd, age_range "text", # control_vision "text", "text", # sample_students, sample_compensation "text", "text", # mod_device, mod_interface_type "text", # mod_time_recoded "text", # mod_aesth_measure "numeric", "numeric", # rel_aesth "text", # mod_aesth_measure_time "text", # mod_aesth_measure_interface "text", # mod_aesth_reference "text", "text", "text", "text", "text", "text", "text", "text", # mod_aesth_manip "numeric", "numeric", # n_aesth "numeric", "numeric", # mean_aesth "numeric", "numeric", # sd_aesth "numeric", "numeric", "text", # g_diff_aesth, g_var_diff_aesth, mod_diff_aesth "text", # mod_task "text", "text", "text", # mod_actual_context, mod_intended_context, mod_mismatch_context "text", "numeric", # mod_perf_measure, direction_within_each_measure "text", # mod_confounders "text", # mod_usability "text", # setting_type "text", # design_type "text", # aesth_conditions "numeric", "numeric", "numeric", "numeric", # mean_perf_aesth, mean_perf_aesth_recoded "numeric", "numeric", # sd_perf_aesth "numeric", "numeric", # cor_perf_aesth "numeric", # eta_perf_aesth "numeric", "numeric", "numeric", #f_statistic_perf_aesth, df1_f_statistic_perf_aesth, df2_f_statistic_perf_aesth "numeric", "numeric", "numeric", "numeric", # d_perf_aesth, d_var_perf_aesth, g_perf_aesth, g_var_perf_aesth "text", # interpret_g_perf_aesth "text", # how_calculated "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", # diad_question1 "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", # diad_question2 "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", # diad_question3 "text", "text", "text", "text", "text", "text", "text", "text", "text", "text", "text"), # diad_question4 na = "NA") # View(Data_AestheticsPerformance_noNA) str(Data_AestheticsPerformance_noNA) ############################################# ### Descriptive statistics of final data set ############################################# # Number of studies, total number of effect sizes, and distribution of effect sizes summary_effects <- Data_AestheticsPerformance_noNA %>% group_by(study_id) %>% summarise(es_id = n()) %>% summarise( studies = n(), effects = sum(es_id), es_min = min(es_id), es_max = max(es_id), es_median = median(es_id), es_q1 = quantile(es_id, .25), es_q3 = quantile(es_id, .75) ) summary_effects # Total number of participants, mean age, and sd of age summary_subjects <- Data_AestheticsPerformance_noNA %>% summarise( total_n = sum(total_n, na.rm = T), age_mean = mean(age_mean, na.rm = T), age_sd = mean(age_sd, na.rm = T), total_n_women = sum(n_women, na.rm = T) ) summary_subjects # Percentage of women total_n_women <- summary_subjects$total_n_women total_n <- summary_subjects$total_n perc_women <- (total_n_women / total_n) * 100 perc_women # 51.12802 # Distribution of total_n summary_total_n <- Data_AestheticsPerformance_noNA %>% summarise( studies = n(), n_min = min(total_n), n_max = max(total_n), n_median = median(total_n), ) summary_total_n # Number of studies with n < 30 Data_AestheticsPerformance_noNA %>% group_by(study_id) %>% summarise( min_n_high = min(n_aesth_high, na.rm = TRUE), min_n_low = min(n_aesth_low, na.rm = TRUE) ) %>% filter(min_n_high < 30 | min_n_low < 30) %>% summarise(studies = n()) # 17 # Range of publication years year_summary <- Data_AestheticsPerformance_noNA %>% summarise( year_min = min(year, na.rm = T), year_max = max(year, na.rm = T)) year_summary # Number of author groups, continents/countries, publication status/type authors_summary <- Data_AestheticsPerformance_noNA %>% mutate(first_author = StrExtract(authors, "^[^, &]+")) %>% # Filtered by first author group_by(first_author) %>% summarise( study_count = n_distinct(study_id), effect_size_count = n() ) print(authors_summary, n = Inf) country_summary <- Data_AestheticsPerformance_noNA %>% group_by(continent_study_recoded, country_study) %>% summarise( study_count = n_distinct(study_id), effect_size_count = n() ) country_summary publication_summary <- Data_AestheticsPerformance_noNA %>% group_by(pub_status, pub_type) %>% summarise( study_count = n_distinct(study_id), effect_size_count = n() ) publication_summary ############################################# ### Main analysis ############################################# ## Multilevel meta-analysis model # Effect sizes (es_id) nested in primary studies (study_id) # Assumption Multilevel model: Correlation of true outcomes ≠ 0, Correlation of sampling errors = 0 within studies # cf. Viechtbauer (2010b): https://wviechtb.github.io/metafor/reference/misc-recs.html multi_model <- rma.mv(yi = g_perf_aesth, V = g_var_perf_aesth, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1| study_id/es_id, # multilevel model test = "t", # usage of t- and F-distributions (Note: standard errors are not adjusted); equivalent to the Knapp-Hartung method (Knapp & Hartung, 2003) (cf. Viechtbauer, 2010b) dfs = "contain", # improved method for approximating dfs of t- and F-distributions (Viechtbauer, 2010b) method = "REML") summary(multi_model) # Estimate Hedge's g = 0.2964; 95% CI [0.0826; 0.5102]; p-value = 0.0082 < 0.01 ## Correlated and hierarchical effects (CHE) model (Pustejovsky & Tipton, 2022) # Assumption CHE model: Correlation of true outcomes ≠ 0, Correlation of sampling errors ≠ 0 (but unknown) within studies # cf. Viechtbauer (2010b): https://wviechtb.github.io/metafor/reference/misc-recs.html # Sensitivity analysis for varying values of the constant sampling correlation (cf. Pustejovsky & Tipton, 2022) rhos <- seq(0, 0.95, 0.05) # values of rho ranging between .00 to .95, in increments of .05 names(rhos) <- rhos # Function to run meta-analysis for a specific value of rho che_model_analysis <- function(rho_value) { # Construct approximate variance-covariance matrix of dependent effect sizes / their sampling errors using vcalc function # cf. Viechtbauer (2010b): https://wviechtb.github.io/metafor/reference/misc-recs.html V <- vcalc(vi = g_var_perf_aesth, cluster = study_id, # cluster = clustering variable obs = es_id, # obs = distinguish different observed ES corresponding to the same construct data = Data_AestheticsPerformance_noNA, rho = rho_value) che_model <- rma.mv(yi = g_perf_aesth, V = V, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id/es_id, test = "t", dfs = "contain", method = "REML") # Calculate confidence intervals with RVE ci <- conf_int(che_model, vcov = "CR2") # Small-sample adjustment "CR2 method" (< 40 primary studies; cf. Tipton & Pustejovsky, 2015) # Extract results (e.g., average effect size, confidence interval) tibble( rho = rho_value, estimate = che_model$beta, se = ci$SE, ci.lb = ci$CI_L, ci.ub = ci$CI_U, tau2 = che_model$tau2 ) } sensitivity_results <- map_dfr(rhos, che_model_analysis) print(sensitivity_results) # Plot results for sensitivity analysis ggplot(sensitivity_results, aes(x = rho, y = estimate)) + geom_line(color = "blue", size = 1) + geom_point(color = "blue", size = 2) + geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = 0.2) + labs(title = "Average effect sizes for varying values of the assumed sampling correlation", x = "Assumed sampling correlation (ρ)", y = "Average effect size (g)") + theme_minimal() + theme( plot.title = element_text(hjust = 0.5), text = element_text(size = 10) ) # Estimate Hedge's g is nearly identical across the range of rho (range of g = 0.296 for rho = 0; g = 0.291 for rho = 0.95); # Calculated ES holds across different estimates for the constant sampling correlation ## CHE model for specific value of rho (rho = 0.5) V_rho_0.5 <- vcalc(vi = g_var_perf_aesth, cluster = study_id, # cluster = clustering variable obs = es_id, # obs = distinguish different observed ES corresponding to the same construct data = Data_AestheticsPerformance_noNA, rho = 0.5) che_model_rho_0.5 <- rma.mv(yi = g_perf_aesth, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id/es_id, test = "t", dfs = "contain", method = "REML") summary(che_model_rho_0.5) # Estimate Hedge's g = 0.2933; 95% CI [0.0772; 0.5095]; p-value = 0.0095 < 0.01 ## Usage of cluster-robust inference method RVE (Robust Variance Estimation) (cf. Assink & Wibbelink, 2016; Tipton & Pustejovsky, 2015; Viechtbauer, 2010b) # Assumption RVE: Exact dependence structure between effect size estimates are unknown # cf. Viechtbauer (2010b): https://wviechtb.github.io/metafor/reference/misc-recs.html# # Apply cluster-robust inference methods (RVE) che_model_RVE <- robust(che_model_rho_0.5, cluster = study_id, adjust = TRUE, # Small-sample adjustment "CR2 method" (< 40 primary studies; cf. Tipton & Pustejovsky, 2015) clubSandwich = TRUE) summary(che_model_RVE) # Estimate Hedge's g = 0.2933; 95% CI [0.0766; 0.5101]; p-value = 0.0097 < 0.01 predict(che_model_RVE) # cf. Viechtbauer (2024b, September 2): https://www.metafor-project.org/doku.php/faq#for_random-effects_models_fitt # 95% Prediction interval of true effect [-1.0739; 1.6606] # Variance components (between- / within-study variance) che_model_RVE$sigma2 # Sigma^2.1 = 0.2546275 # Sigma^2.2 = 0.1811502 # Within-cluster-correlation of true outcomes (= rho) che_model_RVE$sigma2[1] / sum(che_model_RVE$sigma2) # Rho = 0.5843059 # Total heterogeneity (= tau^2) sum(che_model_RVE$sigma2) # Tau^2 = 0.4357778 sqrt(sum(che_model_RVE$sigma2)) # Tau = 0.6601347 # Calculate I2 using formulas by Cheung (2014): https://doi.org/10.1037/a0032968 n <- length(Data_AestheticsPerformance_noNA$g_var_perf_aesth) list.inverse.variances <- 1 / (Data_AestheticsPerformance_noNA$g_var_perf_aesth) sum.inverse.variances <- sum(list.inverse.variances) squared.sum.inverse.variances <- (sum.inverse.variances) ^ 2 list.inverse.variances.square <- 1 / (Data_AestheticsPerformance_noNA$g_var_perf_aesth^2) sum.inverse.variances.square <- sum(list.inverse.variances.square) numerator <- (n - 1) * sum.inverse.variances denominator <- squared.sum.inverse.variances - sum.inverse.variances.square estimated.sampling.variance <- numerator / denominator I2_1 <- (estimated.sampling.variance) / (che_model_RVE$sigma2[1] + che_model_RVE$sigma2[2] + estimated.sampling.variance) # Level 1 variance I2_2 <- (che_model_RVE$sigma2[2]) / (che_model_RVE$sigma2[1] + che_model_RVE$sigma2[2] + estimated.sampling.variance) # Level 2 variance I2_3 <- (che_model_RVE$sigma2[1]) / (che_model_RVE$sigma2[1] + che_model_RVE$sigma2[2] + estimated.sampling.variance) # Level 3 variance amountvariancelevel1 <- I2_1 * 100 amountvariancelevel2 <- I2_2 * 100 amountvariancelevel3 <- I2_3 * 100 amountvariancelevel1 # 10.31644 (10.32%) # Level 1 (Sampling error) = 10.32% (< 75% substantial heterogeneity, cf. Hunter & Schmidt, 1990) amountvariancelevel2 # 37.28093 (37.28%) # Level 2 (Within-study heterogeneity) = 37.28% = medium (cf. Higgins & Thompson, 2002) amountvariancelevel3 # 52.40263 (52.40%) # Level 3 (Between-study heterogeneity) = 52.40% = medium (cf. Higgins & Thompson, 2002) # Total I2 = 37.28 + 52.40 = 89.68% ## Determining the significance of heterogeneity (cf. Harrer et al., 2021) # Perform two separate log-likelihood-ratio tests for within-study variance (level 2) and between-study variance (level 3) # 1) Between-study heterogeneity = 0 novar3_model <- rma.mv(yi = g_perf_aesth, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1| study_id/es_id, test = "t", dfs = "contain", method = "REML", sigma2 = c(NA, 0)) # sigma^2.1 = 0 summary(novar3_model) anova(che_model_rho_0.5,novar3_model) # Three-level model (full) shows better fit (i.e., AIC and BIC are lower; LRT is significant, pval <.0001); # Between-study heterogeneity is significant # 2) Within-study heterogeneity = 0 novar2_model <- rma.mv(yi = g_perf_aesth, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1| study_id/es_id, test = "t", dfs = "contain", method = "REML", sigma2 = c(0, NA)) # sigma^2.2 = 0 summary(novar2_model) anova(che_model_rho_0.5,novar2_model) # Three-level model (full) shows better fit (i.e., AIC and BIC are lower; LRT is significant, pval <.0001); # Within-study heterogeneity is significant ## Run model diagnostics: Detect Outliers and influential cases (cf. Viechtbauer & Cheung, 2010) par(mfrow = c(1, 2)) ## Check for influential cases # Cook's Distance cooksdistance_values <- cooks.distance(che_model_rho_0.5) print(cooksdistance_values) plot(cooksdistance_values, type = "o", pch = 19, xlab = "Observed Outcome", ylab = "Cook's Distance", main = "Cook's Distance") # Study 10, ES 3 = Lavie & Oron-Gilad (2013) # Study 31 = Zhai & Chen (2022) chi_value <- qchisq(0.5, 2) # high values following Cook & Weisberg (1987) chi_value # 1.386294 # Both studies are far below this cut-off ## Check for outliers # Studentized deleted (or externally studentized) residuals rstudent_values <- rstudent(che_model_rho_0.5) # rstudent = studentized deleted residuals plot(rstudent_values$slab, rstudent_values$resid, type = "o", pch = 19, xlab = "Observed Outcome", ylab = "Studentized Deleted Residuals", main = "Studentized Deleted Residuals") high_rstudent_values <- rstudent_values[abs(rstudent_values$resid) > 1.96, ] # residuals > |1.96| (cf. Viechtbauer & Cheung, 2010) high_rstudent_values # Study 23, ES 1-8 = Schnürer et al. (2015) # Study 12, ES 23 + 26 = Ling & van Schaik (2012) # Study 10, ES 1 = Lavie & Oron-Gilad (2013) # No overlap of influential cases and outliers # Therefore, no data exclusion is necessary (cf. Viechtbauer & Cheung, 2010) #### Forest plot #### # pdf(file=paste("ForestPlot.pdf",sep=""),width=10,height=60) cex_text <- 0.8 cex_header <- 0.8 cex_stats <- 0.65 forest(che_model_RVE, slab = paste(Data_AestheticsPerformance_noNA$authors, " (", Data_AestheticsPerformance_noNA$year, ")", sep = ""), at = c(-4, -3, -2, -1, 0, 1, 2, 3, 4), efac = .1, mlab = "CHE model with application of RVE for all studies", xlab = expression("Hedge's " * italic("g")), header= c("Author(s) and year", "Hedge's g [95% CI]"), shade = TRUE, cex = cex_text, cex.lab = cex_header, cex.axis = cex_text) mtext(bquote(paste("(Q = ", .(formatC(che_model_RVE$QE, digits=2, format="f")), ", df = ", .(che_model_RVE$k - che_model_RVE$p), ", p < .001", ")")), side = 1, line = -11, cex = cex_stats, adj = 0.01) # dev.off() ## Forest plot extension: Study DIAD rating # Function to convert the encoding of the DIAD questions convert_diad <- function(value) { if (value == 0) { return("++") # Yes } else if (value == 1) { return("+") # Maybe yes } else if (value == 2) { return("-") # Maybe no } else if (value == 3) { return("--") # No } } pdf(file = paste("StudyDIAD_ForestPlotextension.pdf", sep = ""), width = 5, height = 60) rows <- seq(from = length(che_model_RVE$yi), to = 1, by = -1) # Row positions for each study plot(1, type = "n", xlim = c(0, 5), ylim = c(0, length(rows) + 2), xlab = "", ylab = "", axes = FALSE) # Insert subheadings for diad_questions x_pos_base <- 1.5 # Starting X position for the colored boxes x_spacing <- 1 # Spacing between the boxes cex_heading <- 0.8 cex_text <- 0.8 text(x_pos_base, length(rows) + 2.1, "Q1", cex = cex_heading, font = 2) text(x_pos_base + x_spacing, length(rows) + 2.1, "Q2", cex = cex_heading, font = 2) text(x_pos_base + 2 * x_spacing, length(rows) + 2.1, "Q3", cex = cex_heading, font = 2) text(x_pos_base + 3 * x_spacing, length(rows) + 2.1, "Q4", cex = cex_heading, font = 2) # Horizontal line under the Q1 to Q4 headings line_y <- length(rows) + 1.0 segments(x_pos_base - 0.5, line_y, x_pos_base + 3 * x_spacing + 0.5, line_y, lwd = 1, col = "black") # Insert heading text(x_pos_base + 1.5 * x_spacing, length(rows) + 3.2, "Study DIAD rating", cex = cex_heading, font = 2) # Iterate over all studies and convert the DIAD questions into symbols for (i in seq_along(rows)) { # Convert the DIAD responses into symbols (DIAD questions 1-4) diad_information <- c( convert_diad(Data_AestheticsPerformance_noNA$diad_question1[i]), convert_diad(Data_AestheticsPerformance_noNA$diad_question2[i]), convert_diad(Data_AestheticsPerformance_noNA$diad_question3[i]), convert_diad(Data_AestheticsPerformance_noNA$diad_question4[i]) ) x_pos <- seq(x_pos_base, by = x_spacing, length.out = 4) # Draw color-coded rectangles for (j in 1:length(diad_information)) { diad_col <- match(diad_information[j], c("++", "+", "-", "--")) if (!is.na(diad_col)) { # Draw rectangle rect(x_pos[j] - 0.25, rows[i] - 0.46, x_pos[j] + 0.25, rows[i] + 0.46, density = NA, col = c("#5EB5188C", "#B2E8878C", "#E07B7D8C", "#B5181B8C")[diad_col]) # Insert text into the rectangle text(x_pos[j], rows[i], labels = diad_information[j], cex = cex_text, col = "black") } } } # dev.off() ## Legend for Study DIAD rating # pdf("StudyDIAD_Legend.pdf", width = 6, height = 6) plot(1, type = "n", xlim = c(0, 3), ylim = c(0, 5), xlab = "", ylab = "", axes = FALSE) diad_symbols <- c("++", "+", "-", "--") diad_labels <- c("Yes", "Maybe yes", "Maybe no", "No") diad_colors <- c("#5EB5188C", "#B2E8878C", "#E07B7D8C", "#B5181B8C") symbol_positions <- seq(4, 1, by = -1) # Positionen für die Symbole und Texte for (i in 1:length(diad_symbols)) { rect(0.5, symbol_positions[i] - 0.3, 1, symbol_positions[i] + 0.3, col = diad_colors[i], border = NA) text(0.75, symbol_positions[i], labels = diad_symbols[i], cex = cex_text, font = 2, col = "black") text(1.5, symbol_positions[i], labels = diad_labels[i], pos = 4, cex = cex_text) } legend_box_x <- c(0.2, 2.5, 2.5, 0.2) legend_box_y <- c(0.2, 5, 0.5, 0.5) rect(xleft = min(legend_box_x), ybottom = min(legend_box_y), xright = max(legend_box_x), ytop = max(legend_box_y), border = "black", lwd = 2) title_y_position <- 4.7 text(1.38, title_y_position, "Study DIAD rating", cex = cex_heading, font = 2) # dev.off() ############################################# ### Moderator analysis ############################################# #### Check for multi-collinearity #### categorical_mods <- Data_AestheticsPerformance_noNA[, c("mod_device", "mod_interface_type", # Device, interface type "mod_time_recoded", # Interaction time "mod_aesth_measure", "mod_aesth_measure_time", "mod_aesth_measure_interface", "mod_aesth_reference", # Aesthetic measurement "mod_aesth_manip_color", "mod_aesth_manip_texture", "mod_aesth_manip_layout", "mod_aesth_manip_typo", "mod_aesth_manip_graphics", "mod_aesth_manip_complexity", "mod_aesth_manip_shape", "mod_aesth_manip", # Aesthetic manipulations "mod_diff_aesth", # Difference between aesthetic conditions "mod_task", # Type of task "mod_actual_context", "mod_intended_context", "mod_mismatch_context", # Usage context "mod_perf_measure", # Performance measurement "mod_confounders", # Confounders "mod_usability")] # Usability cramers_v_matrix <- function(cramers_v_values) { mod_names <- colnames(cramers_v_values) n <- length(mod_names) matrix <- matrix(0, n, n, dimnames = list(mod_names, mod_names)) for (i in 1:n) { for (j in 1:n) { if (i == j) { matrix[i, j] <- 1 } else if (i < j) { cramer_v_value <- CramerV(cramers_v_values[[i]], cramers_v_values[[j]]) if (is.na(cramer_v_value)) { matrix[i, j] <- NA matrix[j, i] <- NA } else { matrix[i, j] <- cramer_v_value matrix[j, i] <- cramer_v_value } } } } return(as.data.frame(matrix)) } # Create correlation matrix cramer_matrix <- cramers_v_matrix(categorical_mods) print(cramer_matrix) cramer_matrix <- as.matrix(cramer_matrix) # Visualize correlation matrix (Heatmap) cramer_matrix_melted <- melt(cramer_matrix) ggplot(data = cramer_matrix_melted, aes(x = Var1, y = Var2, fill = value)) + geom_tile() + geom_text(aes(Var2, Var1, label = round(value, 2))) + scale_fill_gradient2(low = "blue", high = "red") # Check for multicollinearity (>= 0.8) high_cramers_v <- which(cramer_matrix >= 0.8, arr.ind = TRUE) high_cramers_v # Cramer's V = 1 for mod_interface_type and mod_aesth_manip_shape # Cramer's V = 0.85 for mod_aesth_manip and mod_aesth_manip_graphics # Cramer's V = 0.87 for mod_aesth_manip and mod_aesth_manip_typo # Therefore, mod_aesth_manip_shape and mod_aesth_manip are excluded from further analysis Data_AestheticsPerformance_noNA$mod_aesth_manip_shape <- NULL Data_AestheticsPerformance_noNA$mod_aesth_manip <- NULL #### Separate meta-regressions for moderator variables #### # Suggestion by Fu et al. (2011): For categorical variables, K = 4 primary studies are needed to include a moderator's category; For continous variables, K = 6 primary studies are needed # Therefore, moderators are excluded from further analysis if K < 4 or K < 6, respectively # For each moderator, separate meta-regression models are estimated using an intercept specification # cf. Viechtbauer (2024a, June 18): https://www.metafor-project.org/doku.php/tips:models_with_or_without_intercept ## mod_device # Levels: 0 = no information about device used, 1 = computer, 2 = mobile phone, 3 = tablet Data_AestheticsPerformance_noNA$mod_device <- factor(Data_AestheticsPerformance_noNA$mod_device) Data_AestheticsPerformance_noNA %>% group_by(mod_device, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Device2 (Mobile phone) with K = 2 primary studies, 6 ES is not further analyzed # Device3 (Tablet) with K = 1 primary studies, 3 ES is not further analyzed # Therefore all devices (computer, mobile phone, and tablet) are merged into one category Data_AestheticsPerformance_noNA <- Data_AestheticsPerformance_noNA %>% mutate(mod_device = recode(mod_device, `0` = 0, `1` = 1, # All devices are merged into one category `2` = 1, `3` = 1 )) Data_AestheticsPerformance_noNA$mod_device <- factor(Data_AestheticsPerformance_noNA$mod_device) levels(Data_AestheticsPerformance_noNA$mod_device) # New levels: 0 = no information about device used, 1 = information about device used model.mod_device <- rma.mv(g_perf_aesth ~ 1 + mod_device, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_device) # No significant results # Apply RVE model.mod_device_RVE <- robust(model.mod_device, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_device_RVE) # Estimate for Device1 (information about device) differs significantly from Device0 (no information about device) (Estimate = 0.3331, p-value = 0.0471 < 0.05) # Check associations with other control variables ## setting_type # Levels: 0 = in person, 1 = online, 2 = mix of settings CramerV(Data_AestheticsPerformance_noNA$mod_device, Data_AestheticsPerformance_noNA$setting_type) # 0.97 ## diad_questions CramerV(Data_AestheticsPerformance_noNA$mod_device, Data_AestheticsPerformance_noNA$diad_question1) # 0.36 CramerV(Data_AestheticsPerformance_noNA$mod_device, Data_AestheticsPerformance_noNA$diad_question2) # 0.37 CramerV(Data_AestheticsPerformance_noNA$mod_device, Data_AestheticsPerformance_noNA$diad_question3) # 0.14 CramerV(Data_AestheticsPerformance_noNA$mod_device, Data_AestheticsPerformance_noNA$diad_question4) # 0.51 ## mod_interface_type # Levels: 0 = software, 1 = website, 2 = prototype/mock-up, 3 = element of interface Data_AestheticsPerformance_noNA$mod_interface_type <- factor(Data_AestheticsPerformance_noNA$mod_interface_type) Data_AestheticsPerformance_noNA %>% group_by(mod_interface_type, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_interface_type <- rma.mv(g_perf_aesth ~ 1 + mod_interface_type, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_interface_type) # Estimate for Interface Type 2 (Prototype/Mock-up) differs significantly from Type0 (Software) (Estimate = 0.9915, p-value = 0.0024 < 0.01) anova(model.mod_interface_type, X=rbind(c(0,-1,1,0),c(0,0,-1,1))) # Estimate for Interface Type2 (Prototype/Mock-up) differs significantly from Type1 (Website) (Estimate = 0.9829, p-value = 0.0008 < 0.001) and # from Type3 (Element of Interface) (Estimate = -0.8963, p-value = 0.0042 < 0.01) # Apply RVE model.mod_interface_type_RVE <- robust(model.mod_interface_type, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_interface_type_RVE) # No significant results anova(model.mod_interface_type_RVE, X=rbind(c(0,-1,1,0),c(0,0,-1,1))) # No significant results ## mod_time_recoded # Levels: 0 = few seconds, 1 = few minutes to an hour, 2 = few hours/days/weeks Data_AestheticsPerformance_noNA$mod_time_recoded <- factor(Data_AestheticsPerformance_noNA$mod_time_recoded) Data_AestheticsPerformance_noNA %>% group_by(mod_time_recoded, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Time0 (few seconds) with K = 1 primary studies, 2 ES is not further analyzed; # Time2 (few hours/days/weeks) with K = 1 primary studies, 16 ES is not further analyzed; # Therefore, the moderator mod_time_recoded is dropped from the model Data_AestheticsPerformance_noNA$mod_time_recoded <- NULL ## mod_aesth_measure # Levels: 0 = single-item scales, 1 = multi-item scales with low-level validation, 2 = multi-item scales with high-level validation, # 3 = scale by Lavie & Tractinsky (2004), 4 = VisAWI, 5 = appeal items of AttrakDiff 1, 6 = other measurements Data_AestheticsPerformance_noNA$mod_aesth_measure <- factor(Data_AestheticsPerformance_noNA$mod_aesth_measure) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_measure, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Measure5 (AttrakDiff) with K = 2 primary studies, 14 ES is not further analyzed # Measure6 (Other) with K = 1 primary studies, 44 ES is not further analyzed Data_AestheticsPerformance_noNA$mod_aesth_measure[Data_AestheticsPerformance_noNA$mod_aesth_measure == 5] <- NA Data_AestheticsPerformance_noNA$mod_aesth_measure[Data_AestheticsPerformance_noNA$mod_aesth_measure == 6] <- NA Data_AestheticsPerformance_noNA$mod_aesth_measure <- factor(Data_AestheticsPerformance_noNA$mod_aesth_measure) # Measure1 (Multi-item scale with low-level validation) with K = 2 primary studies, 20 ES is not further analyzed # Measure2 (Multi-item scale with high-level validation) with K = 2 primary studies, 10 ES is not further analyzed # Therefore both moderator with multi-item scales are merged into one category Data_AestheticsPerformance_noNA <- Data_AestheticsPerformance_noNA %>% mutate(mod_aesth_measure = recode(mod_aesth_measure, `0` = 0, `1` = 1, # Multi-item scales (1, 2) are merged into one category `2` = 1, `3` = 2, `4` = 3 )) Data_AestheticsPerformance_noNA$mod_aesth_measure <- factor(Data_AestheticsPerformance_noNA$mod_aesth_measure) levels(Data_AestheticsPerformance_noNA$mod_aesth_measure) # New levels: 0 = single-item scales, 1 = multi-item scales (with low and high level validation), 2 = scale by Lavie & Tractinsky (2004), 3 = VisAWI model.mod_aesth_measure <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_measure, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_measure) # Intercept / Estimate for Measure0 (Single-item scales) differs significantly from zero (Estimate = 0.4020; p-value = 0.0186 < 0.05) anova(model.mod_aesth_measure, X=rbind(c(0,-1,1,0),c(0,-1,0,1),c(0,0,-1,1))) # No significant results # Multi-item scales are merged into one category to increase number of studies (K) Data_AestheticsPerformance_noNA <- Data_AestheticsPerformance_noNA %>% mutate(mod_aesth_measure = recode(mod_aesth_measure, `0` = 0, `1` = 1, # Multi-item scales, scale by Lavie & Tractinsky and VisAWI are merged into one subgroup `2` = 1, `3` = 1 )) Data_AestheticsPerformance_noNA$mod_aesth_measure <- factor(Data_AestheticsPerformance_noNA$mod_aesth_measure) # New levels: 0 = single-item scales, 1 = multi-item scales (low and high validation, scale by LAvie & Tractinsky, VisAWI) model.mod_aesth_measure <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_measure, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_measure) # Intercept / Estimate for Measure0 (Single-item scales) differs significantly from zero (Estimate = 0.4020; p-value = 0.0149 < 0.05) # Apply RVE model.mod_aesth_measure_RVE <- robust(model.mod_aesth_measure, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_measure_RVE) # No significant results ## mod_aesth_measure_time # Levels: 0 = in previous study, 1 = before interaction, 2 = after interaction Data_AestheticsPerformance_noNA$mod_aesth_measure_time <- factor(Data_AestheticsPerformance_noNA$mod_aesth_measure_time) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_measure_time, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_measure_time <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_measure_time, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_measure_time) # No significant results anova(model.mod_aesth_measure_time, X=rbind(c(0,-1,1))) # No significant results # Apply RVE model.mod_aesth_measure_time_RVE <- robust(model.mod_aesth_measure_time, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_measure_time_RVE) # Intercept / Estimate for Measure Time0 (in previous study) is significantly different from zero (Estimate = 0.3223; p-value = 0.0394 < 0.05) anova(model.mod_aesth_measure_time_RVE, X=rbind(c(0,-1,1))) # No significant results ## mod_aesth_measure_interface # Levels: 0 = no interface presented, 1 = interface presented Data_AestheticsPerformance_noNA$mod_aesth_measure_interface <- factor(Data_AestheticsPerformance_noNA$mod_aesth_measure_interface) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_measure_interface, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_measure_interface <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_measure_interface, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_measure_interface) # Intercept / Estimate for Interface0 (not presented) differs significantly from zero (Estimate = 0.3177; p-value = 0.0204 < 0.05) # Apply RVE model.mod_aesth_measure_interface_RVE <- robust(model.mod_aesth_measure_interface, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_measure_interface_RVE) # Intercept / Estimate for Interface0 (not presented) differs significantly from zero (Estimate = 0.3117; p-value = 0.0084 < 0.01) ## mod_aesth_reference # Levels: 0 = collective assessment, 1 = individual assessment Data_AestheticsPerformance_noNA$mod_aesth_reference <- factor(Data_AestheticsPerformance_noNA$mod_aesth_reference) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_reference, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Eval1 (Individual) with K = 3 primary studies, 37 ES is not further analyzed; # Therefore, the moderator mod_aesth_eval is dropped from the model Data_AestheticsPerformance_noNA$mod_aesth_reference <- NULL ## mod_aesth_manip # color # Levels: 0 = no color manipulation, 1 = color manipulation Data_AestheticsPerformance_noNA$mod_aesth_manip_color <- factor(Data_AestheticsPerformance_noNA$mod_aesth_manip_color) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_manip_color, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_manip_color <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_manip_color, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_manip_color) # Intercept / Estimate for Color0 (no manipulation) differs significantly from zero (Estimate = 0.5313; p-value = 0.0111 < 0.05) # Apply RVE model.mod_aesth_manip_color_RVE <- robust(model.mod_aesth_manip_color, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_manip_color_RVE) # No significant results # texture # Levels: 0 = no manipulation in texture, 1 = manipulation in texture Data_AestheticsPerformance_noNA$mod_aesth_manip_texture <- factor(Data_AestheticsPerformance_noNA$mod_aesth_manip_texture) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_manip_texture, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_manip_texture <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_manip_texture, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_manip_texture) # Intercept / Estimate for Texture0 (no manipulation) differs significantly from zero (Estimate = 0.3152; p-value = 0.0095 < 0.01) # Apply RVE model.mod_aesth_manip_texture_RVE <- robust(model.mod_aesth_manip_texture, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_manip_texture_RVE) # Intercept / Estimate for Texture0 (no manipulation) differs significantly from zero (Estimate = 0.3152; p-value = 0.0117) # layout # Levels: 0 = no manipulation in layout, 1 = manipulation in layout Data_AestheticsPerformance_noNA$mod_aesth_manip_layout <- factor(Data_AestheticsPerformance_noNA$mod_aesth_manip_layout) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_manip_layout, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_manip_layout <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_manip_layout, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_manip_layout) # No significant results # Apply RVE model.mod_aesth_manip_layout_RVE <- robust(model.mod_aesth_manip_layout, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_manip_layout_RVE) # Intercept / Estimate for Layout0 (no manipulation) differs significantly from zero (Estimate = 0.1902; p-value = 0.0057) # typo # Levels: 0 = no manipulation in typography, 1 = manipulation in typography Data_AestheticsPerformance_noNA$mod_aesth_manip_typo <- factor(Data_AestheticsPerformance_noNA$mod_aesth_manip_typo) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_manip_typo, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_manip_typo <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_manip_typo, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_manip_typo) # Intercept / Estimate for Typo0 (no manipulation) differs significantly from zero (Estimate = 0.3767; p-value = 0.0031 < 0.01) # Apply RVE model.mod_aesth_manip_typo_RVE <- robust(model.mod_aesth_manip_typo, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_manip_typo_RVE) # Intercept / Estimate for Typo0 (no manipulation) differs significantly from zero (Estimate = 0.3767; p-value = 0.0079) # Estimate for Typo1 (manipulation) differs significantly from Typo0 (no manipulation) (Estimate = -0.3977; p-value = 0.0191) # graphics # Levels: 0 = no manipulation in graphics, 1 = manipulation in graphics Data_AestheticsPerformance_noNA$mod_aesth_manip_graphics <- factor(Data_AestheticsPerformance_noNA$mod_aesth_manip_graphics) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_manip_graphics, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_manip_graphics <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_manip_graphics, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_manip_graphics) # Intercept / Estimate for Graphics0 (no manipulation) is significantly different from zero (Estimate = 0.3319; p-value = 0.0065 < 0.01) # Apply RVE model.mod_aesth_manip_graphics_RVE <- robust(model.mod_aesth_manip_graphics, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_manip_graphics_RVE) # Intercept / Estimate for Graphics0 (no manipulation) is significantly different from zero (Estimate = 0.3319; p-value = 0.0097 < 0.01) # complexity # Levels: 0 = no manipulation in complexity, 1 = manipulation in complexity Data_AestheticsPerformance_noNA$mod_aesth_manip_complexity <- factor(Data_AestheticsPerformance_noNA$mod_aesth_manip_complexity) Data_AestheticsPerformance_noNA %>% group_by(mod_aesth_manip_complexity, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_aesth_manip_complexity <- rma.mv(g_perf_aesth ~ 1 + mod_aesth_manip_complexity, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_aesth_manip_complexity) # Intercept / Estimate for Complexity0 (no manipulation) is significantly different from zero (Estimate = 0.3249; p-value = 0.0178 < 0.05) # Apply RVE model.mod_aesth_manip_complexity_RVE <- robust(model.mod_aesth_manip_complexity, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_aesth_manip_complexity_RVE) # Intercept / Estimate for Complexity0 (no manipulation) is significantly different from zero (Estimate = 0.3249; p-value = 0.0399 < 0.05) ## mod_diff_aesth # Levels: 0 = small, 1 = medium, 2 = large Data_AestheticsPerformance_noNA$mod_diff_aesth <- factor(Data_AestheticsPerformance_noNA$mod_diff_aesth) Data_AestheticsPerformance_noNA %>% group_by(mod_diff_aesth, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_diff_aesth <- rma.mv(g_perf_aesth ~ 1 + mod_diff_aesth, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_diff_aesth) # No significant results anova(model.mod_diff_aesth, X=rbind(c(0,-1,1))) # Estimate Diff_aesth2 (large) differs significantly from Diff_aesth1 (medium) (Estimate = 0.3314, p-value = 0.0250) # Apply RVE model.mod_diff_aesth_RVE <- robust(model.mod_diff_aesth, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_diff_aesth_RVE) # No significant results anova(model.mod_diff_aesth_RVE, X=rbind(c(0,-1,1))) # No significant results (Difference between Diff_aesth1 and Diff_aesth2 is significant on a trend level; Estimate = 0.3314, p-value = 0.0761) ## mod_task # Levels: 0 = free use, 1 = search task, 2 = use function Data_AestheticsPerformance_noNA$mod_task <- factor(Data_AestheticsPerformance_noNA$mod_task) Data_AestheticsPerformance_noNA %>% group_by(mod_task, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_task <- rma.mv(g_perf_aesth ~ 1 + mod_task, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_task) # No significant results anova(model.mod_task, X=rbind(c(0,-1,1))) # No significant results # Apply RVE model.mod_task_RVE <- robust(model.mod_task, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_task_RVE) # No significant results anova(model.mod_task_RVE, X=rbind(c(0,-1,1))) # No significant results ## mod_actual_context # Levels: 0 = leisure, 1 = learning Data_AestheticsPerformance_noNA$mod_actual_context <- factor(Data_AestheticsPerformance_noNA$mod_actual_context) Data_AestheticsPerformance_noNA %>% group_by(mod_actual_context, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_actual_context <- rma.mv(g_perf_aesth ~ 1 + mod_actual_context, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_actual_context) # No significant results # Apply RVE model.mod_actual_context_RVE <- robust(model.mod_actual_context, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_actual_context_RVE) # No significant results ## mod_intended_context # Levels: 0 = leisure, 1 = learning, 2 = work, 3 = other Data_AestheticsPerformance_noNA$mod_intended_context <- factor(Data_AestheticsPerformance_noNA$mod_intended_context) Data_AestheticsPerformance_noNA %>% group_by(mod_intended_context, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Context3 (Other) with K = 2 primary studies, 28 ES is not further analyzed Data_AestheticsPerformance_noNA$mod_intended_context[Data_AestheticsPerformance_noNA$mod_intended_context == 3] <- NA Data_AestheticsPerformance_noNA$mod_intended_context <- factor(Data_AestheticsPerformance_noNA$mod_intended_context) model.mod_intended_context <- rma.mv(g_perf_aesth ~ 1 + mod_intended_context, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_intended_context) # No significant results anova(model.mod_intended_context, X=rbind(c(0,-1,1))) # No significant results # Apply RVE model.mod_intended_context_RVE <- robust(model.mod_intended_context, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_intended_context_RVE) # Intercept / Estimate for Context0 (Leisure) differs significantly from zero (Estimate = 0.3135, p-value = 0.0264) anova(model.mod_intended_context_RVE, X=rbind(c(0,-1,1))) # No significant results ## mod_mismatch_context # Levels: 0 = no (i.e., the actual context of use matches the intended one), 1 = yes (i.e., the actual context of use does not match the intended one) Data_AestheticsPerformance_noNA$mod_mismatch_context <- factor(Data_AestheticsPerformance_noNA$mod_mismatch_context) Data_AestheticsPerformance_noNA %>% group_by(mod_mismatch_context, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_mismatch_context <- rma.mv(g_perf_aesth ~ 1 + mod_mismatch_context, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_mismatch_context) # Intercept / Estimate for Context0 (no mismatch) differs significantly from zero (Estimate = 0.3935, p-value = 0.0192) # Apply RVE model.mod_mismatch_context_RVE <- robust(model.mod_mismatch_context, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_mismatch_context_RVE) # No significant results ## mod_perf_measure # Levels: 0 = speed, 1 = accuracy, 2 = efficiency, 3 = learning output, 4 = depth of interaction, 5 = other measurements Data_AestheticsPerformance_noNA$mod_perf_measure <- factor(Data_AestheticsPerformance_noNA$mod_perf_measure) Data_AestheticsPerformance_noNA %>% group_by(mod_perf_measure, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Perf3 (Learning) with K = 3 primary studies, 10 ES is not further analyzed # Perf4 (Depth of interaction) with K = 2 primary studies, 10 ES is not further analyzed # Perf5 (Other) with K = 1 primary study, 1 ES is not further analyzed Data_AestheticsPerformance_noNA$mod_perf_measure[Data_AestheticsPerformance_noNA$mod_perf_measure == 3] <- NA Data_AestheticsPerformance_noNA$mod_perf_measure[Data_AestheticsPerformance_noNA$mod_perf_measure == 4] <- NA Data_AestheticsPerformance_noNA$mod_perf_measure[Data_AestheticsPerformance_noNA$mod_perf_measure == 5] <- NA Data_AestheticsPerformance_noNA$mod_perf_measure <- factor(Data_AestheticsPerformance_noNA$mod_perf_measure) model.mod_perf_measure <- rma.mv(g_perf_aesth ~ 1 + mod_perf_measure, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_perf_measure) # Intercept / Estimate for Measure0 (Speed) differs significantly from zero (Estimate = 0.3891; p-value = 0.0029 < 0.01) # Estimate for Measure1 (Accuracy) differs significantly from Measure0 (Speed) (Estimate = -0.2913; p-value = 0.0028 < 0.01) anova(model.mod_perf_measure, X=rbind(c(0,-1,1))) # Estimate for Measure1 (Accuracy) differs significantly from Measure2 (Efficiency) (Estimate = 0.3350; p-value = 0.0168) # Apply RVE model.mod_perf_measure_RVE <- robust(model.mod_perf_measure, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_perf_measure_RVE) # Intercept / Estimate for Measure0 (Speed) differs significantly from zero (Estimate = 0.3891; p-value = 0.0148 < 0.05) anova(model.mod_perf_measure_RVE, X=rbind(c(0,-1,1))) # No significant results ## mod_confounders # Levels: 0 = not controlled, 1 = controlled for the subjectively experienced usability when using the interface, # 2 = controlled for multiple variables, including experienced usability Data_AestheticsPerformance_noNA$mod_confounders <- factor(Data_AestheticsPerformance_noNA$mod_confounders) Data_AestheticsPerformance_noNA %>% group_by(mod_confounders, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.mod_confounders <- rma.mv(g_perf_aesth ~ 1 + mod_confounders, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.mod_confounders) # Intercept / Estimate for Confounders0 (not controlled) differs significantly from zero (Estimate = 0.3926; p-value = 0.0171 < 0.05) anova(model.mod_confounders, X=rbind(c(0,-1,1))) # No significant results # Apply RVE model.mod_confounders_RVE <- robust(model.mod_confounders, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.mod_confounders_RVE) # No significant results anova(model.mod_confounders_RVE, X=rbind(c(0,-1,1))) # No significant results ## mod_usability # Levels: 0 = good, 1 = poor Data_AestheticsPerformance_noNA$mod_usability <- factor(Data_AestheticsPerformance_noNA$mod_usability) Data_AestheticsPerformance_noNA %>% group_by(mod_usability, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Usability1 (poor) with K = 3 primary studies, 11 ES is not further analyzed # Therefore, the moderator mod_usability is dropped from the model Data_AestheticsPerformance_noNA$mod_usability <- NULL #### Model with all significant moderators #### CramerV(Data_AestheticsPerformance_noNA$mod_device, Data_AestheticsPerformance_noNA$mod_aesth_manip_typo) # 0.49 model.allmods <- rma.mv(g_perf_aesth ~ 1 + mod_device + mod_aesth_manip_typo, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.allmods) # Apply RVE model.allmods_RVE <- robust(model.allmods, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.allmods_RVE) # No significant results ## Cramer's V with final coding ## # mod_time_recoded, mod_aesth_reference, and mod_usability cannot be analyzed due to a lack of primary studies (k < 4) # mod_aesth_manip_shape, mod_aesth_manip are excluded due to multicollinearity (V > 0.8) categorical_mods_new <- Data_AestheticsPerformance_noNA[, c("mod_device", "mod_interface_type", # Device, interface type #"mod_time_recoded", # Interaction time "mod_aesth_measure", "mod_aesth_measure_time", "mod_aesth_measure_interface", #"mod_aesth_reference", # Aesthetic measurement "mod_aesth_manip_color", "mod_aesth_manip_texture", "mod_aesth_manip_layout", "mod_aesth_manip_typo", "mod_aesth_manip_graphics", "mod_aesth_manip_complexity", # "mod_aesth_manip_shape", "mod_aesth_manip", # Aesthetic manipulations "mod_diff_aesth", # Difference between aesthetic conditions "mod_task", # Type of task "mod_actual_context", "mod_intended_context", "mod_mismatch_context", # Usage context "mod_perf_measure", # Performance measurement "mod_confounders")] # Confounders #"mod_usability")] # Usability cramers_v_matrix_new <- function(cramers_v_values) { mod_names <- colnames(cramers_v_values) n <- length(mod_names) matrix <- matrix(0, n, n, dimnames = list(mod_names, mod_names)) for (i in 1:n) { for (j in 1:n) { if (i == j) { matrix[i, j] <- 1 } else if (i < j) { cramer_v_value <- CramerV(cramers_v_values[[i]], cramers_v_values[[j]]) if (is.na(cramer_v_value)) { matrix[i, j] <- NA matrix[j, i] <- NA } else { matrix[i, j] <- cramer_v_value matrix[j, i] <- cramer_v_value } } } } return(as.data.frame(matrix)) } # Create correlation matrix cramer_matrix_new <- cramers_v_matrix(categorical_mods_new) print(cramer_matrix_new) cramer_matrix_new <- as.matrix(cramer_matrix_new) # Visualize correlation matrix (Heatmap) cramer_matrix_melted_new <- melt(cramer_matrix_new) ggplot(data = cramer_matrix_melted_new, aes(x = Var1, y = Var2, fill = value)) + geom_tile() + geom_text(aes(Var2, Var1, label = round(value, 2))) + scale_fill_gradient2(low = "blue", high = "red") # Check for multicollinearity (>= 0.8) high_cramers_v_new <- which(cramer_matrix_new >= 0.8, arr.ind = TRUE) high_cramers_v_new #### Control variables #### #### Check for multi-collinearity #### categorical_controls <- Data_AestheticsPerformance_noNA[, c("pub_status", "pub_type", "continent_study_recoded", # Reference "control_vision", "sample_students", "sample_compensation", # Sample "design_type", "setting_type", "aesth_conditions", "how_calculated_recoded")] # Method cramers_v_matrix <- function(cramers_v_values) { mod_names <- colnames(cramers_v_values) n <- length(mod_names) matrix <- matrix(0, n, n, dimnames = list(mod_names, mod_names)) for (i in 1:n) { for (j in 1:n) { if (i == j) { matrix[i, j] <- 1 } else if (i < j) { cramer_v_value <- CramerV(cramers_v_values[[i]], cramers_v_values[[j]]) if (is.na(cramer_v_value)) { matrix[i, j] <- NA matrix[j, i] <- NA } else { matrix[i, j] <- cramer_v_value matrix[j, i] <- cramer_v_value } } } } return(as.data.frame(matrix)) } # Create correlation matrix cramer_matrix_controls <- cramers_v_matrix(categorical_controls) print(cramer_matrix_controls) cramer_matrix_controls <- as.matrix(cramer_matrix_controls) # Visualize correlation matrix (Heatmap) cramer_matrix_melted <- melt(cramer_matrix_controls) ggplot(data = cramer_matrix_melted, aes(x = Var1, y = Var2, fill = value)) + geom_tile() + geom_text(aes(Var2, Var1, label = round(value, 2))) + scale_fill_gradient2(low = "blue", high = "red") # Check for multicollinearity (>= 0.8) high_cramers_v <- which(cramer_matrix_controls >= 0.8, arr.ind = TRUE) high_cramers_v # No correlation >= 0.8 #### Separate meta-regressions for control variables #### ## Reference ## pub_status # Levels: 0 = published, 1 = unpublished Data_AestheticsPerformance_noNA$pub_status <- factor(Data_AestheticsPerformance_noNA$pub_status) Data_AestheticsPerformance_noNA %>% group_by(pub_status, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Pub_status1 (unpublished) with K = 1 primary study, 6 ES is not further analyzed # Therefore, the control variable is dropped from the model Data_AestheticsPerformance_noNA$pub_status <- NULL ## pub_type # Levels: 0 = peer-reviewed journal article, 1 = conference paper, 2 = book chapter, 3 = thesis/dissertation Data_AestheticsPerformance_noNA$pub_type <- factor(Data_AestheticsPerformance_noNA$pub_type) Data_AestheticsPerformance_noNA %>% group_by(pub_type, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Pub_type1 (conference paper) with K = 2 primary studies, 3 ES is not further analyzed # Pub_type2 (book chapter) with K = 2 primary studies, 10 ES is not further analyzed # Pub_type3 (thesis/dissertation) with K = 3 primary studies, 11 ES is not further analyzed # Therefore, the control variable is dropped from the model Data_AestheticsPerformance_noNA$pub_type <- NULL ## year # Grand Mean Centering ("CGM") # n_women = Level-1 predictor; i.e. varies between ES within a study sum(!is.na(Data_AestheticsPerformance_noNA$year)) studies_year <- Data_AestheticsPerformance_noNA %>% group_by(study_id) %>% summarise(noNA = any(!is.na(year)))%>% filter(noNA == TRUE) %>% select(study_id) nrow(studies_year) Data_AestheticsPerformance_noNA$year <- as.numeric(Data_AestheticsPerformance_noNA$year) round(mean(Data_AestheticsPerformance_noNA$year, na.rm = TRUE), 0) # Grand Mean = 2013 Data_AestheticsPerformance_noNA$year_centered <- center(Data_AestheticsPerformance_noNA$year, type = "CGM") model.year_centered <- rma.mv(g_perf_aesth ~ 1 + year_centered, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.year_centered) # Intercept / Estimate for average publication year (= 2013) differs significantly from zero (Estimate = 0.2966, p-value = 0.0124 < 0.05) # Apply RVE model.year_centered_RVE <- robust(model.year_centered, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.year_centered_RVE) # Intercept / Estimate for average publication year (= 2013) differs significantly from zero (Estimate = 0.2966, p-value = 0.0102 < 0.05) ## continent_study(_recoded) # Levels: 0 = Europe, 1 = North America, 2 = Asia Data_AestheticsPerformance_noNA$continent_study_recoded <- factor(Data_AestheticsPerformance_noNA$continent_study_recoded) Data_AestheticsPerformance_noNA %>% group_by(continent_study_recoded, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.continent_study_recoded <- rma.mv(g_perf_aesth ~ 1 + continent_study_recoded, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.continent_study_recoded) # Intercept / Estimate for Continent_study0 (Europe) differs significantly from zero (Estimate = 0.2716; p-value = 0.0392 < 0.05) anova(model.continent_study_recoded, X=rbind(c(0,-1,1))) # No significant results # Apply RVE model.continent_study_recoded_RVE <- robust(model.continent_study_recoded, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.continent_study_recoded_RVE) # No significant results anova(model.continent_study_recoded_RVE, X=rbind(c(0,-1,1))) # No significant results ## Sample ## n_women # Grand Mean Centering ("CGM") # n_women = Level-1 predictor; i.e. varies between ES within a study sum(!is.na(Data_AestheticsPerformance_noNA$n_women)) studies_n_women <- Data_AestheticsPerformance_noNA %>% group_by(study_id) %>% summarise(noNA = any(!is.na(n_women)))%>% filter(noNA == TRUE) %>% select(study_id) nrow(studies_n_women) Data_AestheticsPerformance_noNA$n_women <- as.numeric(Data_AestheticsPerformance_noNA$n_women) round(mean(Data_AestheticsPerformance_noNA$n_women, na.rm = TRUE), 0) # Grand Mean = 47 Data_AestheticsPerformance_noNA$n_women_centered <- center(Data_AestheticsPerformance_noNA$n_women, type = "CGM") model.n_women_centered <- rma.mv(g_perf_aesth ~ 1 + n_women_centered, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.n_women_centered) # Intercept / Estimate for average number of women (= 47) differs significantly from zero (Estimate = 0.2789, p-value = 0.0328 < 0.05) # Apply RVE model.n_women_centered_RVE <- robust(model.n_women_centered, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.n_women_centered_RVE) # Intercept / Estimate for average number of women (= 47) differs significantly from zero (Estimate = 0.2789, p-value = 0.0352 < 0.05) ## age_mean # Grand Mean Centering ("CGM") # age_mean = Level-1 predictor; i.e. varies between ES within a study sum(!is.na(Data_AestheticsPerformance_noNA$age_mean)) studies_age_mean <- Data_AestheticsPerformance_noNA %>% group_by(study_id) %>% summarise(noNA = any(!is.na(age_mean)))%>% filter(noNA == TRUE) %>% select(study_id) nrow(studies_age_mean) Data_AestheticsPerformance_noNA$age_mean <- as.numeric(Data_AestheticsPerformance_noNA$age_mean) round(mean(Data_AestheticsPerformance_noNA$age_mean, na.rm = TRUE), 0) # Grand Mean = 26 Data_AestheticsPerformance_noNA$age_mean_centered <- center(Data_AestheticsPerformance_noNA$age_mean, type = "CGM") model.age_mean_centered <- rma.mv(g_perf_aesth ~ 1 + age_mean_centered, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.age_mean_centered) # Intercept / Estimate for average age (= 26 years) differs significantly from zero (Estimate = 0.2890, p-value = 0.0150 < 0.05) # Apply RVE model.age_mean_centered_RVE <- robust(model.age_mean_centered, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.age_mean_centered_RVE) # Intercept / Estimate for average age (= 26 years) differs significantly from zero (Estimate = 0.2890, p-value = 0.0179 < 0.05) ## control_vision # Levels: 0 = no information about vision, 1 = sample with normal/fully corrected vision Data_AestheticsPerformance_noNA$control_vision <- factor(Data_AestheticsPerformance_noNA$control_vision) Data_AestheticsPerformance_noNA %>% group_by(control_vision, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.control_vision <- rma.mv(g_perf_aesth ~ 1 + control_vision, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.control_vision) # Intercept / Estimate for control_vision0 (no information) differs significantly from zero (Estimate = 0.3016, p-value = 0.0162 < 0.05) # Apply RVE model.control_vision_RVE <- robust(model.control_vision, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.control_vision_RVE) # Intercept / Estimate for control_vision0 (no information) differs significantly from zero (Estimate = 0.3016, p-value = 0.0486 < 0.05) ## sample_students # Levels: 0 = assured student sample, 1 = potential student sample, 2 = no student sample Data_AestheticsPerformance_noNA$sample_students <- factor(Data_AestheticsPerformance_noNA$sample_students) Data_AestheticsPerformance_noNA %>% group_by(sample_students, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.sample_students <- rma.mv(g_perf_aesth ~ 1 + sample_students, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.sample_students) # Intercept / Estimate for Students0 (assured student sample) differs significantly from zero (Estimate = 0.2796, p-value = 0.0340 < 0.05) # Apply RVE model.sample_students_RVE <- robust(model.sample_students, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.sample_students_RVE) # Intercept / Estimate for Students0 (assured student sample) differs significantly from zero (Estimate = 0.2796, p-value = 0.0049 < 0.01) ## sample_compensation # Levels: 0 = no information about compensation received, 1 = compensation received (e.g., monetary, course credits) Data_AestheticsPerformance_noNA$sample_compensation <- factor(Data_AestheticsPerformance_noNA$sample_compensation) Data_AestheticsPerformance_noNA %>% group_by(sample_compensation, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.sample_compensation <- rma.mv(g_perf_aesth ~ 1 + sample_compensation, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.sample_compensation) # No significant results # Apply RVE model.sample_compensation_RVE <- robust(model.sample_compensation, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.sample_compensation_RVE) # Intercept / Estimate for Compensation0 (assured student sample) differs significantly from zero (Estimate = 0.3218, p-value = 0.0054 < 0.01) ## Method ## design_type # Levels: 0 = between-subject, 1 = within-subject Data_AestheticsPerformance_noNA$design_type <- factor(Data_AestheticsPerformance_noNA$design_type) Data_AestheticsPerformance_noNA %>% group_by(design_type, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.design_type <- rma.mv(g_perf_aesth ~ 1 + design_type, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.design_type) # No significant results # Apply RVE model.design_type_RVE <- robust(model.design_type, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.design_type_RVE) # No significant results ## setting_type # Levels: 0 = in person setting, 1 = online setting, 2 = mix of both settings Data_AestheticsPerformance_noNA$setting_type <- factor(Data_AestheticsPerformance_noNA$setting_type) Data_AestheticsPerformance_noNA %>% group_by(setting_type, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Implementation2 (mix) with K = 1 primary study, 4 ES is not further analyzed Data_AestheticsPerformance_noNA$setting_type[Data_AestheticsPerformance_noNA$setting_type == 2] <- NA Data_AestheticsPerformance_noNA$setting_type <- factor(Data_AestheticsPerformance_noNA$setting_type) model.setting_type <- rma.mv(g_perf_aesth ~ 1 + setting_type, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.setting_type) # Intercept / Estimate for Implementation0 (in person) differs significantly from zero (Estimate = 0.3684, p-value = 0.0067) # Apply RVE model.setting_type_RVE <- robust(model.setting_type, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.setting_type_RVE) # Intercept / Estimate for Implementation0 (in person) differs significantly from zero (Estimate = 0.3684, p-value = 0.0203) ## aesth_conditions # Levels: 0 = a priori definition of conditions, 1 = a posteriori definition of conditions Data_AestheticsPerformance_noNA$aesth_conditions <- factor(Data_AestheticsPerformance_noNA$aesth_conditions) Data_AestheticsPerformance_noNA %>% group_by(aesth_conditions, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) model.aesth_conditions <- rma.mv(g_perf_aesth ~ 1 + aesth_conditions, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.aesth_conditions) # Estimate for aesth_conditions1 (a posteriori) differs significantly from aesth_conditions0 (a priori) (Estimate = 0.8486, p-value <.0001) # Apply RVE model.aesth_conditions_RVE <- robust(model.aesth_conditions, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.aesth_conditions_RVE) # Estimate for aesth_conditions1 (a posteriori) differs significantly from aesth_conditions0 (a priori) (Estimate = 0.8486, p-value = 0.0390 < 0.05) # Check associations with other moderator variables ## mod_diff_aesth # Levels: 0 = small, 1 = medium, 2 = large CramerV(Data_AestheticsPerformance_noNA$aesth_conditions, Data_AestheticsPerformance_noNA$mod_diff_aesth) # 0.08761828 ## mod_confounders # Levels: 0 = not controlled, 1 = controlled for the subjectively experienced usability when using the interface, # 2 = controlled for multiple variables, including experienced usability CramerV(Data_AestheticsPerformance_noNA$aesth_conditions, Data_AestheticsPerformance_noNA$mod_confounders) # 0.4564959 # rel_aesth_high, rel_aesth_low # Grand Mean Centering ("CGM") # rel_aesth = Level-1 predictor; i.e. varies between ES within a study sum(!is.na(Data_AestheticsPerformance_noNA$rel_aesth_high)) studies_rel_aesth_high <- Data_AestheticsPerformance_noNA %>% group_by(study_id) %>% summarise(noNA = any(!is.na(rel_aesth_high)))%>% filter(noNA == TRUE) %>% select(study_id) nrow(studies_rel_aesth_high) sum(!is.na(Data_AestheticsPerformance_noNA$rel_aesth_low)) studies_rel_aesth_low <- Data_AestheticsPerformance_noNA %>% group_by(study_id) %>% summarise(noNA = any(!is.na(rel_aesth_low)))%>% filter(noNA == TRUE) %>% select(study_id) nrow(studies_rel_aesth_low) Data_AestheticsPerformance_noNA$rel_aesth_high <- as.numeric(Data_AestheticsPerformance_noNA$rel_aesth_high) Data_AestheticsPerformance_noNA$rel_aesth_low <- as.numeric(Data_AestheticsPerformance_noNA$rel_aesth_low) mean(Data_AestheticsPerformance_noNA$rel_aesth_high, na.rm = TRUE) # Grand Mean = 0.8773469 mean(Data_AestheticsPerformance_noNA$rel_aesth_low, na.rm = TRUE) # Grand Mean = 0.8895918 Data_AestheticsPerformance_noNA$rel_aesth_high_centered <- center(Data_AestheticsPerformance_noNA$rel_aesth_high, type = "CGM") Data_AestheticsPerformance_noNA$rel_aesth_low_centered <- center(Data_AestheticsPerformance_noNA$rel_aesth_low, type = "CGM") model.rel_aesth_high_centered <- rma.mv(g_perf_aesth ~ 1 + rel_aesth_high_centered, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.rel_aesth_high_centered) # Intercept / Estimate for average reliability (= 0.877) differs significantly from zero (Estimate = 0.2132, p-value = 0.0462 < 0.05) # Apply RVE model.rel_aesth_high_centered_RVE <- robust(model.rel_aesth_high_centered, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.rel_aesth_high_centered_RVE) # No significant results model.rel_aesth_low_centered <- rma.mv(g_perf_aesth ~ 1 + rel_aesth_low_centered, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.rel_aesth_low_centered) # No significant results # Apply RVE model.rel_aesth_low_centered_RVE <- robust(model.rel_aesth_low_centered, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.rel_aesth_low_centered_RVE) # No significant results ## how_calculated(_recoded) # Levels: 0 = Calculation from raw data, 1 = Calculation from M, SD and N per condition, 2 = Calculation from correlation r, # 3 = Calculation from F-test statistics or eta squared, 4 = Calculation from Cohen’s d Data_AestheticsPerformance_noNA$how_calculated_recoded <- factor(Data_AestheticsPerformance_noNA$how_calculated_recoded) Data_AestheticsPerformance_noNA %>% group_by(how_calculated_recoded, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # How_calculated0 (raw data) with K = 1 primary study, 16 ES is not further analyzed Data_AestheticsPerformance_noNA$how_calculated_recoded[Data_AestheticsPerformance_noNA$how_calculated_recoded == 0] <- NA Data_AestheticsPerformance_noNA$how_calculated_recoded <- factor(Data_AestheticsPerformance_noNA$how_calculated_recoded) # How_calculated2 (correlation) with K = 3 primary studies, 8 ES is not further analyzed # How_calculated3 (F-statistics) with K = 1 primary study, 2 ES is not further analyzed # How_calculated4 (Cohen's d) with K = 3 primary studies, 18 ES is not further analyzed # Therefore, all control variables with higher aggregated data than M and SD are merged into one category Data_AestheticsPerformance_noNA <- Data_AestheticsPerformance_noNA %>% mutate(how_calculated_recoded = recode(how_calculated_recoded, `1` = 0, `2` = 1, # Calculation from r, F-Test or eta-squared, and d are merged into one category `3` = 1, `4` = 1 )) Data_AestheticsPerformance_noNA$how_calculated_recoded <- factor(Data_AestheticsPerformance_noNA$how_calculated_recoded) # New levels: 0 = Calculation from M, SD and N per condition, 1 = Calculation from correlation r, from F-test statistics or eta squared, or from Cohen’s d model.how_calculated_recoded <- rma.mv(g_perf_aesth ~ 1 + how_calculated_recoded, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.how_calculated_recoded) # Intercept / Estimate for how_calculated0 (in person) differs significantly from zero (Estimate = 0.2924, p-value = 0.0155 < 0.05) # Apply RVE model.how_calculated_recoded_RVE <- robust(model.how_calculated_recoded, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.how_calculated_recoded_RVE) # Intercept / Estimate for how_calculated0 (in person) differs significantly from zero (Estimate = 0.2924, p-value = 0.0240 < 0.05) #### Model with all significant control variables #### # Aesth_conditions is the only significant control variable (see model results above) #### Study Quality (Study DIAD) #### #### Check for multi-collinearity #### categorical_studydiad <- Data_AestheticsPerformance_noNA[, c("diad_question1", "diad_question2", "diad_question3", "diad_question4")] cramers_v_matrix <- function(cramers_v_values) { mod_names <- colnames(cramers_v_values) n <- length(mod_names) matrix <- matrix(0, n, n, dimnames = list(mod_names, mod_names)) for (i in 1:n) { for (j in 1:n) { if (i == j) { matrix[i, j] <- 1 } else if (i < j) { cramer_v_value <- CramerV(cramers_v_values[[i]], cramers_v_values[[j]]) if (is.na(cramer_v_value)) { matrix[i, j] <- NA matrix[j, i] <- NA } else { matrix[i, j] <- cramer_v_value matrix[j, i] <- cramer_v_value } } } } return(as.data.frame(matrix)) } # Create correlation matrix cramer_matrix_studydiad <- cramers_v_matrix(categorical_studydiad) print(cramer_matrix_studydiad) cramer_matrix_studydiad <- as.matrix(cramer_matrix_studydiad) # Visualize correlation matrix (Heatmap) cramer_matrix_melted <- melt(cramer_matrix_studydiad) ggplot(data = cramer_matrix_melted, aes(x = Var1, y = Var2, fill = value)) + geom_tile() + geom_text(aes(Var2, Var1, label = round(value, 2))) + scale_fill_gradient2(low = "blue", high = "red") # Check for multicollinearity (>= 0.8) high_cramers_v <- which(cramer_matrix_studydiad >= 0.8, arr.ind = TRUE) high_cramers_v # No correlation >= 0.8 #### Separate meta-regressions for study quality variables #### ## Construct validity: diad_question1 # Levels: 0 = Yes, 2 = Maybe no Data_AestheticsPerformance_noNA$diad_question1 <- factor(Data_AestheticsPerformance_noNA$diad_question1) Data_AestheticsPerformance_noNA %>% group_by(diad_question1, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Diad_question10 (Yes): K = 17, 145 ES # Diad_question12 (Maybe no): K = 14, 89 ES model.diad_question1 <- rma.mv(g_perf_aesth ~ 1 + diad_question1, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.diad_question1) # No significant results # Apply RVE model.diad_question1_RVE <- robust(model.diad_question1, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.diad_question1_RVE) # Intercept / Estimate for diad_question10 (yes) differs significantly from zero (Estimate = 0.1331, p-value = 0.0075 < 0.01) ## Internal validity: diad_question2 # Levels: 0 = Yes, 1 = Maybe yes, 2 = Maybe no Data_AestheticsPerformance_noNA$diad_question2 <- factor(Data_AestheticsPerformance_noNA$diad_question2) Data_AestheticsPerformance_noNA %>% group_by(diad_question2, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Diad_question20 (Yes): K = 29, 218 ES # Diad_question21 (Maybe yes) with K = 1 primary study, 10 ES is not further analyzed # Diad_question22 (Maybe no) with K = 1 primary studies, 6 ES is not further analyzed # Therefore, the control variable is dropped from the model Data_AestheticsPerformance_noNA$diad_question2 <- NULL ## External validity: diad_question3 # Levels: 2 = Maybe no, 3 = No Data_AestheticsPerformance_noNA$diad_question3 <- factor(Data_AestheticsPerformance_noNA$diad_question3) Data_AestheticsPerformance_noNA %>% group_by(diad_question3, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Diad_question32 (Maybe no): K = 27, 213 ES # Diad_question33 (No): K = 4, 21 ES model.diad_question3 <- rma.mv(g_perf_aesth ~ 1 + diad_question3, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.diad_question3) # Intercept / Estimate for diad_question32 (maybe no) differs significantly from zero (Estimate = 0.2736, p-value = 0.0238 < 0.05) # Apply RVE model.diad_question3_RVE <- robust(model.diad_question3, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.diad_question3_RVE) # Intercept / Estimate for diad_question32 (maybe no) differs significantly from zero (Estimate = 0.2736, p-value = 0.0286 < 0.05) ## Statistical validity: diad_question4 # Levels: 0 = Yes, 1 = Maybe yes, 2 = Maybe no, 3 = No Data_AestheticsPerformance_noNA$diad_question4 <- factor(Data_AestheticsPerformance_noNA$diad_question4) Data_AestheticsPerformance_noNA %>% group_by(diad_question4, study_id) %>% summarise(effects = n()) %>% summarise( studies = n(), effects = sum(effects)) # Diad_question40 (Yes): K = 9, 59 ES # Diad_question41 (Maybe yes): K = 4, 51 ES # Diad_question42 (Maybe no): K = 7, 40 ES # Diad_question43 (No): K = 11, 84 ES model.diad_question4 <- rma.mv(g_perf_aesth ~ 1 + diad_question4, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.diad_question4) # No significant results # "Yes" and "Maybe yes" as well as "Maybe no" and "No" are merged into one category each to increase number of studies (K) Data_AestheticsPerformance_noNA <- Data_AestheticsPerformance_noNA %>% mutate(diad_question4 = recode(diad_question4, `0` = 0, # # "Yes" and "Maybe yes" are merged into one category `1` = 0, `2` = 1, # "Maybe no" and "No" are merged into one category `3` = 1 )) Data_AestheticsPerformance_noNA$diad_question4 <- factor(Data_AestheticsPerformance_noNA$diad_question4) # New levels: 0 = Yes or Maybe yes, 1 = Maybe no or No model.diad_question4 <- rma.mv(g_perf_aesth ~ 1 + diad_question4, V = V_rho_0.5, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1 | study_id / es_id, test = "t", dfs = "contain", method = "REML") summary(model.diad_question4) # No significant results # Apply RVE model.diad_question4_RVE <- robust(model.diad_question4, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(model.diad_question4_RVE) # Intercept / Estimate for diad_question40 (yes + maybe yes) differs significantly from zero (Estimate = 0.1805, p-value = 0.0152 < 0.05) #### Check associations with other control / moderator variables #### ## setting_type # Levels: 0 = in person setting, 1 = online setting CramerV(Data_AestheticsPerformance_noNA$setting_type, Data_AestheticsPerformance_noNA$diad_question1) # 0.3755184 CramerV(Data_AestheticsPerformance_noNA$setting_type, Data_AestheticsPerformance_noNA$diad_question1.1.2) # 0.3755184 ## mod_mismatch_context # Levels: 0 = no, 1 = yes CramerV(Data_AestheticsPerformance_noNA$mod_mismatch_context, Data_AestheticsPerformance_noNA$diad_question3.1.1) # 0.2942598 ## year lm_diad_question1 <- lm(year ~ diad_question1, data = Data_AestheticsPerformance_noNA) sqrt(summary(lm_diad_question1)$r.squared) # 0.0316588 # diad_question2 cannot be analyzed to due lack of studies in categories lm_diad_question3 <- lm(year ~ diad_question3, data = Data_AestheticsPerformance_noNA) sqrt(summary(lm_diad_question3)$r.squared) # 0.03185143 lm_diad_question4 <- lm(year ~ diad_question4, data = Data_AestheticsPerformance_noNA) sqrt(summary(lm_diad_question4)$r.squared) # 0.243198 ############################################# ### Publication Bias ############################################# ## Funnel Plot # Color coding of effect sizes nested in studies author_year_data <- Data_AestheticsPerformance_noNA %>% select(authors, year) %>% distinct() legend_labels <- paste(author_year_data$authors, author_year_data$year, sep = " (") legend_labels <- paste(legend_labels, ")", sep = "") num_studies <- length(unique(Data_AestheticsPerformance_noNA$study_id)) # 31 studies custom_colors <- colorRampPalette(brewer.pal(12, "Set1"))(num_studies) # define color palette using RColorBrewer # Alternative: Paired, Accent study_colors <- custom_colors[as.factor(Data_AestheticsPerformance_noNA$study_id)] # assign unique color to each study # png(filename = "FunnelPlot_colored.png", units = "in", width = 10, height = 8.0, res = 300, pointsize = 16) par(mfrow = c(1, 2)) funnel(che_model_RVE, yaxis = "sei", xlab = expression("Hedge's"~italic(g)), main = "(a) Standard Error", col = study_colors, bg = study_colors) funnel(che_model_RVE, yaxis = "seinv", xlab = expression("Hedge's"~italic(g)), main = "(b) Inverse Standard Error", col = study_colors, bg = study_colors) # dev.off() # png(filename = "FunnelPlot_colored_legend.png", units = "in", width = 12, height = 8.0, res = 300, pointsize = 16) par(mfrow = c(1, 1)) plot.new() legend("center", legend = legend_labels, fill = custom_colors, xpd = TRUE, cex = 0.8) # dev.off() ## Egger's test for asymmetry # Egger's regression test for multilevel models (Egger MLMA): Including the standard error of the effect size (or a related measure of precision, e.g. sample size) as a moderator (cf. Rodgers & Pustejovsky, 2021) # When SMD as effect size is calculated: sample size or "effective sample size" should be used as moderator (cf. Nakagawa et al., 2021) # cf. Viechtbauer (2017, June 25): https://stats.stackexchange.com/questions/155693/metafor-package-bias-and-sensitivity-diagnostics ## Egger's test for asymmetry using standard errors of effect sizes eggertest <- rma.mv(yi = g_perf_aesth, V = V_rho_0.5, mods = ~ sqrt(g_var_perf_aesth), # standard error = square root of sampling variances slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1| study_id/es_id, test = "t", dfs = "contain", method = "REML") eggertest_RVE <- robust(eggertest, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(eggertest_RVE) # Moderator sqrt(g_var_perf_aesth) does not have a significant influence on the model (p = 0.4752); # Test for asymmetry is negative ## Egger's test for asymmetry using the "effective sample size" (cf. Nakagawa et al., 2021) # calculating "effective sample size" (esz) to account for unbalanced sampling (see Equation 25) Data_AestheticsPerformance_noNA$SMD.esz <- (4*Data_AestheticsPerformance_noNA$n_aesth_high*Data_AestheticsPerformance_noNA$n_aesth_low) / (Data_AestheticsPerformance_noNA$n_aesth_high + Data_AestheticsPerformance_noNA$n_aesth_low) # creating "effective sample size" based "sampling variance" (see Equation 26) Data_AestheticsPerformance_noNA$SMD.esz.sv <- 4/Data_AestheticsPerformance_noNA$SMD.esz # creating corresponding "standard error" (i.e. the square root of the sampling variance) Data_AestheticsPerformance_noNA$SMD.esz.sei <- sqrt(Data_AestheticsPerformance_noNA$SMD.esz.sv) eggertest.SMD.se <- rma.mv(yi = g_perf_aesth, V = V_rho_0.5, mods = ~ SMD.esz.sei, slab = study_id, data = Data_AestheticsPerformance_noNA, random = ~ 1| study_id/es_id, test = "t", dfs = "contain", method = "REML") eggertest.SMD.se_RVE <- robust(eggertest.SMD.se, cluster = study_id, adjust = TRUE, clubSandwich = TRUE) summary(eggertest.SMD.se_RVE) # Moderator SMD.esz.sei does not have a significant influence on the model (p = 0.6969); # Test for asymmetry is negative