diff --git a/.gitignore b/.gitignore index 791e731..416c4ff 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ .history/ -eohi3_2.csv \ No newline at end of file +eohi3_2.csv +*~ +.~lock* \ No newline at end of file diff --git a/eohi3/DA00_fmaxVALS.r b/eohi3/DA00_fmaxVALS.r new file mode 100644 index 0000000..f3fd22e --- /dev/null +++ b/eohi3/DA00_fmaxVALS.r @@ -0,0 +1,149 @@ +library(SuppDists) +library(dplyr) +library(tidyr) + +setwd("/home/ladmin/Documents/DND/EOHI/eohi3") + +between_vars <- c("perspective", "temporalDO") +within_vars_MEAN <- c( + "past_pref_MEAN", "past_pers_MEAN", "past_val_MEAN", + "fut_pref_MEAN", "fut_pers_MEAN", "fut_val_MEAN" +) +within_vars_DGEN <- c( + "past_pref_DGEN", "past_pers_DGEN", "past_val_DGEN", + "fut_pref_DGEN", "fut_pers_DGEN", "fut_val_DGEN" +) + +df <- read.csv("eohi3.csv", stringsAsFactors = FALSE, check.names = FALSE, na.strings = "NA") + +anova_data_MEAN <- df %>% + select(pID, all_of(between_vars), all_of(within_vars_MEAN)) %>% + filter(!is.na(perspective), perspective != "", + !is.na(temporalDO), temporalDO != "") + +long_data_MEAN <- anova_data_MEAN %>% + pivot_longer( + cols = all_of(within_vars_MEAN), + names_to = "variable", + values_to = "MEAN_SCORE" + ) %>% + mutate( + time = ifelse(grepl("^past_", variable), "past", "fut"), + domain = case_when( + grepl("_pref_MEAN$", variable) ~ "pref", + grepl("_pers_MEAN$", variable) ~ "pers", + grepl("_val_MEAN$", variable) ~ "val", + TRUE ~ NA_character_ + ) + ) %>% + mutate( + TIME = factor(time, levels = c("past", "fut")), + DOMAIN = factor(domain, levels = c("pref", "pers", "val")), + perspective = factor(perspective), + temporalDO = factor(temporalDO) + ) %>% + select(pID, perspective, temporalDO, TIME, DOMAIN, MEAN_SCORE) %>% + filter(!is.na(MEAN_SCORE)) + +cell_vars_MEAN <- long_data_MEAN %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + variance = var(MEAN_SCORE, na.rm = TRUE), + .groups = "drop" + ) + +fmax_by_cell_MEAN <- cell_vars_MEAN %>% + group_by(TIME, DOMAIN) %>% + summarise( + Fmax_observed = max(variance, na.rm = TRUE) / min(variance, na.rm = TRUE), + df_min = min(n) - 1L, + .groups = "drop" + ) + +k <- 4 + +fmax_table_MEAN <- fmax_by_cell_MEAN %>% + rowwise() %>% + mutate( + alpha_0.05 = SuppDists::qmaxFratio(0.95, df = df_min, k = k), + alpha_0.01 = SuppDists::qmaxFratio(0.99, df = df_min, k = k) + ) %>% + ungroup() %>% + mutate( + Fmax_observed = round(Fmax_observed, 4), + alpha_0.05 = round(alpha_0.05, 4), + alpha_0.01 = round(alpha_0.01, 4) + ) %>% + select(TIME, DOMAIN, Fmax_observed, alpha_0.05, alpha_0.01) + +# ---- MEAN: Print observed Hartley ratios ---- +cat("\n--- Hartley ratios (MEAN) ---\n") +fmax_table_MEAN %>% + mutate(across(where(is.numeric), ~ format(round(., 4), nsmall = 4))) %>% + print() + +# ---- DGEN: Observed Hartley ratios ---- +anova_data_DGEN <- df %>% + select(pID, all_of(between_vars), all_of(within_vars_DGEN)) %>% + filter(!is.na(perspective), perspective != "", + !is.na(temporalDO), temporalDO != "") + +long_data_DGEN <- anova_data_DGEN %>% + pivot_longer( + cols = all_of(within_vars_DGEN), + names_to = "variable", + values_to = "DGEN_SCORE" + ) %>% + mutate( + time = ifelse(grepl("^past_", variable), "past", "fut"), + domain = case_when( + grepl("_pref_DGEN$", variable) ~ "pref", + grepl("_pers_DGEN$", variable) ~ "pers", + grepl("_val_DGEN$", variable) ~ "val", + TRUE ~ NA_character_ + ) + ) %>% + mutate( + TIME = factor(time, levels = c("past", "fut")), + DOMAIN = factor(domain, levels = c("pref", "pers", "val")), + perspective = factor(perspective), + temporalDO = factor(temporalDO) + ) %>% + select(pID, perspective, temporalDO, TIME, DOMAIN, DGEN_SCORE) %>% + filter(!is.na(DGEN_SCORE)) + +cell_vars_DGEN <- long_data_DGEN %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + variance = var(DGEN_SCORE, na.rm = TRUE), + .groups = "drop" + ) + +fmax_by_cell_DGEN <- cell_vars_DGEN %>% + group_by(TIME, DOMAIN) %>% + summarise( + Fmax_observed = max(variance, na.rm = TRUE) / min(variance, na.rm = TRUE), + df_min = min(n) - 1L, + .groups = "drop" + ) + +fmax_table_DGEN <- fmax_by_cell_DGEN %>% + rowwise() %>% + mutate( + alpha_0.05 = SuppDists::qmaxFratio(0.95, df = df_min, k = k), + alpha_0.01 = SuppDists::qmaxFratio(0.99, df = df_min, k = k) + ) %>% + ungroup() %>% + mutate( + Fmax_observed = round(Fmax_observed, 4), + alpha_0.05 = round(alpha_0.05, 4), + alpha_0.01 = round(alpha_0.01, 4) + ) %>% + select(TIME, DOMAIN, Fmax_observed, alpha_0.05, alpha_0.01) + +cat("\n--- Hartley ratios (DGEN) ---\n") +fmax_table_DGEN %>% + mutate(across(where(is.numeric), ~ format(round(., 4), nsmall = 4))) %>% + print() \ No newline at end of file diff --git a/eohi3/DA01_anova_DS.r b/eohi3/DA01_anova_DS.r new file mode 100644 index 0000000..b6f9fad --- /dev/null +++ b/eohi3/DA01_anova_DS.r @@ -0,0 +1,235 @@ +library(tidyverse) +library(rstatix) +library(emmeans) +library(effectsize) +library(afex) +library(car) + +options(scipen = 999) +afex::set_sum_contrasts() + +setwd("/home/ladmin/Documents/DND/EOHI/eohi3") + +df <- read.csv( + "eohi3.csv", + stringsAsFactors = FALSE, + check.names = FALSE, + na.strings = "NA" +) + +between_vars <- c("perspective", "temporalDO") +within_vars <- c( + "past_pref_MEAN", "past_pers_MEAN", "past_val_MEAN", + "fut_pref_MEAN", "fut_pers_MEAN", "fut_val_MEAN" +) + +missing_vars <- setdiff(c(between_vars, within_vars, "pID"), names(df)) +if (length(missing_vars) > 0) { + stop(paste("Missing required variables:", paste(missing_vars, collapse = ", "))) +} + +anova_data <- df %>% + select(pID, all_of(between_vars), all_of(within_vars)) %>% + filter( + !is.na(perspective), perspective != "", + !is.na(temporalDO), temporalDO != "" + ) + +long_data <- anova_data %>% + pivot_longer( + cols = all_of(within_vars), + names_to = "variable", + values_to = "MEAN_SCORE" + ) %>% + mutate( + time = case_when( + grepl("^past_", variable) ~ "past", + grepl("^fut_", variable) ~ "fut", + TRUE ~ NA_character_ + ), + domain = case_when( + grepl("_pref_MEAN$", variable) ~ "pref", + grepl("_pers_MEAN$", variable) ~ "pers", + grepl("_val_MEAN$", variable) ~ "val", + TRUE ~ NA_character_ + ) + ) %>% + mutate( + TIME = factor(time, levels = c("past", "fut")), + DOMAIN = factor(domain, levels = c("pref", "pers", "val")), + perspective = factor(perspective), + temporalDO = factor(temporalDO), + pID = factor(pID) + ) %>% + select(pID, perspective, temporalDO, TIME, DOMAIN, MEAN_SCORE) %>% + filter(!is.na(MEAN_SCORE)) + +desc_stats <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + mean = round(mean(MEAN_SCORE), 5), + variance = round(var(MEAN_SCORE), 5), + sd = round(sd(MEAN_SCORE), 5), + median = round(median(MEAN_SCORE), 5), + q1 = round(quantile(MEAN_SCORE, 0.25), 5), + q3 = round(quantile(MEAN_SCORE, 0.75), 5), + min = round(min(MEAN_SCORE), 5), + max = round(max(MEAN_SCORE), 5), + .groups = "drop" + ) + +print(desc_stats, n = Inf) + +missing_summary <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n_total = n(), + n_missing = sum(is.na(MEAN_SCORE)), + pct_missing = round(100 * n_missing / n_total, 2), + .groups = "drop" + ) + +print(missing_summary, n = Inf) + +outlier_summary <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + n_outliers = sum(abs(scale(MEAN_SCORE)) > 3), + .groups = "drop" + ) + +print(outlier_summary, n = Inf) + +homogeneity_between <- long_data %>% + group_by(TIME, DOMAIN) %>% + rstatix::levene_test(MEAN_SCORE ~ perspective * temporalDO) + +print(homogeneity_between, n = Inf) + +# Normality: within-subjects residuals (deviation from each participant's mean) +resid_within <- long_data %>% + group_by(pID) %>% + mutate(person_mean = mean(MEAN_SCORE, na.rm = TRUE)) %>% + ungroup() %>% + mutate(resid = MEAN_SCORE - person_mean) %>% + pull(resid) +resid_within <- resid_within[!is.na(resid_within)] +n_resid <- length(resid_within) +if (n_resid < 3L) { + message("Too few within-subjects residuals (n < 3); skipping Shapiro-Wilk.") +} else { + resid_for_shapiro <- if (n_resid > 5000L) { + set.seed(1L) + sample(resid_within, 5000L) + } else { + resid_within + } + print(shapiro.test(resid_for_shapiro)) +} +# qqnorm(resid_within) +# qqline(resid_within) + +aov_afex <- aov_ez( + id = "pID", + dv = "MEAN_SCORE", + data = long_data, + between = c("perspective", "temporalDO"), + within = c("TIME", "DOMAIN"), + type = 3 +) + +# ANOVA table: uncorrected and Greenhouse–Geisser +cat("\n--- ANOVA Table (Type 3, uncorrected) ---\n") +print(nice(aov_afex, correction = "none")) +cat("\n--- ANOVA Table (Type 3, Greenhouse–Geisser correction) ---\n") +print(nice(aov_afex, correction = "GG")) + +# Mauchly's test of sphericity and epsilon (via car::Anova on wide data) +anova_wide <- anova_data %>% + select(pID, perspective, temporalDO, all_of(within_vars)) %>% + filter(if_all(all_of(within_vars), ~ !is.na(.))) +response_matrix <- as.matrix(anova_wide[, within_vars]) +rm_model <- lm(response_matrix ~ perspective * temporalDO, data = anova_wide) +idata <- data.frame( + TIME = factor(rep(c("past", "fut"), each = 3), levels = c("past", "fut")), + DOMAIN = factor(rep(c("pref", "pers", "val"), 2), levels = c("pref", "pers", "val")) +) +rm_anova <- car::Anova(rm_model, idata = idata, idesign = ~ TIME * DOMAIN, type = 3) +rm_summary <- summary(rm_anova, multivariate = FALSE) +if (!is.null(rm_summary$sphericity.tests)) { + cat("\nMauchly's Test of Sphericity:\n") + print(rm_summary$sphericity.tests) +} +if (!is.null(rm_summary$epsilon)) { + cat("\nEpsilon (GG, HF):\n") + print(rm_summary$epsilon) +} + +# Within-subjects residuals: deviation from each participant's mean (one per observation) +resid_within <- long_data %>% + group_by(pID) %>% + mutate(person_mean = mean(MEAN_SCORE, na.rm = TRUE)) %>% + ungroup() %>% + mutate(resid = MEAN_SCORE - person_mean) %>% + pull(resid) +resid_within <- resid_within[!is.na(resid_within)] +# R's shapiro.test() allows 3 <= n <= 5000; use a random sample of 5000 if we have more +n_resid <- length(resid_within) +if (n_resid < 3L) { + message("Too few within-subjects residuals (n < 3); skipping Shapiro-Wilk.") +} else { + resid_for_shapiro <- if (n_resid > 5000L) { + set.seed(1L) + sample(resid_within, 5000L) + } else { + resid_within + } + print(shapiro.test(resid_for_shapiro)) +} + +# qqnorm(resid_within) +# qqline(resid_within) + +# POST-HOC COMPARISONS (significant effects only) + +# TIME (main effect) +emm_TIME <- emmeans(aov_afex, ~ TIME) +print(pairs(emm_TIME, adjust = "bonferroni")) + +# temporalDO:TIME — ~TIME and ~temporalDO +emm_temporalDO_TIME <- emmeans(aov_afex, ~ TIME | temporalDO) +print(pairs(emm_temporalDO_TIME, adjust = "bonferroni")) +emm_temporalDO_temporalDO <- emmeans(aov_afex, ~ temporalDO | TIME) +print(pairs(emm_temporalDO_temporalDO, adjust = "bonferroni")) + +# perspective:temporalDO:TIME — ~TIME, ~perspective, ~temporalDO +emm_pt_TIME <- emmeans(aov_afex, ~ TIME | perspective + temporalDO) +print(pairs(emm_pt_TIME, adjust = "bonferroni")) +emm_pt_perspective <- emmeans(aov_afex, ~ perspective | temporalDO + TIME) +print(pairs(emm_pt_perspective, adjust = "bonferroni")) +emm_pt_temporalDO <- emmeans(aov_afex, ~ temporalDO | perspective + TIME) +print(pairs(emm_pt_temporalDO, adjust = "bonferroni")) + +# perspective:DOMAIN — ~perspective and ~DOMAIN +emm_perspective_DOMAIN <- emmeans(aov_afex, ~ perspective | DOMAIN) +print(pairs(emm_perspective_DOMAIN, adjust = "bonferroni")) +emm_perspective_DOMAIN_domain <- emmeans(aov_afex, ~ DOMAIN | perspective) +print(pairs(emm_perspective_DOMAIN_domain, adjust = "bonferroni")) + +# perspective:TIME:DOMAIN — ~TIME, ~perspective, ~DOMAIN +emm_pt_TIME_domain <- emmeans(aov_afex, ~ TIME | perspective + DOMAIN) +print(pairs(emm_pt_TIME_domain, adjust = "bonferroni")) +emm_pt_domain_perspective <- emmeans(aov_afex, ~ perspective | TIME + DOMAIN) +print(pairs(emm_pt_domain_perspective, adjust = "bonferroni")) +emm_pt_domain_domain <- emmeans(aov_afex, ~ DOMAIN | perspective + TIME) +print(pairs(emm_pt_domain_domain, adjust = "bonferroni")) + +# perspective:temporalDO:TIME:DOMAIN — ~TIME, ~perspective, ~temporalDO +emm_ptt_TIME <- emmeans(aov_afex, ~ TIME | perspective + temporalDO + DOMAIN) +print(pairs(emm_ptt_TIME, adjust = "bonferroni")) +emm_ptt_perspective <- emmeans(aov_afex, ~ perspective | temporalDO + TIME + DOMAIN) +print(pairs(emm_ptt_perspective, adjust = "bonferroni")) +emm_ptt_temporalDO <- emmeans(aov_afex, ~ temporalDO | perspective + TIME + DOMAIN) +print(pairs(emm_ptt_temporalDO, adjust = "bonferroni")) diff --git a/eohi3/knit/DA01_anova_DS.Rmd b/eohi3/knit/DA01_anova_DS.Rmd new file mode 100644 index 0000000..f8870c4 --- /dev/null +++ b/eohi3/knit/DA01_anova_DS.Rmd @@ -0,0 +1,501 @@ +--- +title: "Mixed ANOVA - Domain Specific Means (DA01)" +author: "" +date: "`r Sys.Date()`" +output: + html_document: + toc: true + toc_float: true + code_folding: hide +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = TRUE) +``` + +# Setup + +```{r libraries} +library(tidyverse) +library(rstatix) +library(emmeans) +library(effectsize) +library(afex) +library(car) + +options(scipen = 999) +afex::set_sum_contrasts() +``` + +# Data + +```{r read-data} +# Data file is in parent of knit/ (eohi3/eohi3.csv) +df <- read.csv( + "/home/ladmin/Documents/DND/EOHI/eohi3/eohi3.csv", + stringsAsFactors = FALSE, + check.names = FALSE, + na.strings = "NA" +) + +between_vars <- c("perspective", "temporalDO") +within_vars <- c( + "past_pref_MEAN", "past_pers_MEAN", "past_val_MEAN", + "fut_pref_MEAN", "fut_pers_MEAN", "fut_val_MEAN" +) + +missing_vars <- setdiff(c(between_vars, within_vars, "pID"), names(df)) +if (length(missing_vars) > 0) { + stop(paste("Missing required variables:", paste(missing_vars, collapse = ", "))) +} + +anova_data <- df %>% + select(pID, all_of(between_vars), all_of(within_vars)) %>% + filter( + !is.na(perspective), perspective != "", + !is.na(temporalDO), temporalDO != "" + ) +``` + +# Long format + +```{r long-format} +long_data <- anova_data %>% + pivot_longer( + cols = all_of(within_vars), + names_to = "variable", + values_to = "MEAN_SCORE" + ) %>% + mutate( + time = case_when( + grepl("^past_", variable) ~ "past", + grepl("^fut_", variable) ~ "fut", + TRUE ~ NA_character_ + ), + domain = case_when( + grepl("_pref_MEAN$", variable) ~ "pref", + grepl("_pers_MEAN$", variable) ~ "pers", + grepl("_val_MEAN$", variable) ~ "val", + TRUE ~ NA_character_ + ) + ) %>% + mutate( + TIME = factor(time, levels = c("past", "fut")), + DOMAIN = factor(domain, levels = c("pref", "pers", "val")), + perspective = factor(perspective), + temporalDO = factor(temporalDO), + pID = factor(pID) + ) %>% + select(pID, perspective, temporalDO, TIME, DOMAIN, MEAN_SCORE) %>% + filter(!is.na(MEAN_SCORE)) +``` + +# Descriptive statistics + +```{r desc-stats} +desc_stats <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + mean = round(mean(MEAN_SCORE), 5), + variance = round(var(MEAN_SCORE), 5), + sd = round(sd(MEAN_SCORE), 5), + median = round(median(MEAN_SCORE), 5), + q1 = round(quantile(MEAN_SCORE, 0.25), 5), + q3 = round(quantile(MEAN_SCORE, 0.75), 5), + min = round(min(MEAN_SCORE), 5), + max = round(max(MEAN_SCORE), 5), + .groups = "drop" + ) + +# Show all rows and columns (no truncation) +options(tibble.width = Inf) +print(desc_stats, n = Inf) +``` + +Interpretations: +1. Mean and median values are similar, indicating distribution is relatively symmetric and any skew is minimal. Any outliers are not extreme. +2. Highest to lowest group n size ratio is 1.14 (139/122). Acceptable ratio for ANOVA (under 1.5). +3. Highest to lowest overall group variance ratio is 1.67 (7.93/4.74). Acceptable ratio for ANOVA (under 4). + For the sake of consistency w/ the other EHI studies, I also calculated Hartley's F-max ratio. + The conservative F-max critical value is 1.60, which is still higher than the highest observed F-max ratio of 1.53. + +# Assumption checks + +## Missing values + +```{r missing} +missing_summary <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n_total = n(), + n_missing = sum(is.na(MEAN_SCORE)), + pct_missing = round(100 * n_missing / n_total, 2), + .groups = "drop" + ) + +print(missing_summary, n = Inf) +``` + +No missing values. As expected. + +## Outliers + +```{r outliers} +outlier_summary <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + n_outliers = sum(abs(scale(MEAN_SCORE)) > 3), + .groups = "drop" + ) + +print(outlier_summary, n = Inf) +``` + +No outliers present. Good. + +## Homogeneity of variance + +```{r homogeneity} +homogeneity_between <- long_data %>% + group_by(TIME, DOMAIN) %>% + rstatix::levene_test(MEAN_SCORE ~ perspective * temporalDO) + +print(homogeneity_between, n = Inf) +``` + +Levene's test is sigifnicant for two cells: fut-pers and fut-val. +However, variance ratios and F-max tests show that any heteroscedasticity is mild. + +## Normality (within-subjects residuals) + +```{r normality} +resid_within <- long_data %>% + group_by(pID) %>% + mutate(person_mean = mean(MEAN_SCORE, na.rm = TRUE)) %>% + ungroup() %>% + mutate(resid = MEAN_SCORE - person_mean) %>% + pull(resid) +resid_within <- resid_within[!is.na(resid_within)] + +n_resid <- length(resid_within) +if (n_resid < 3L) { + message("Too few within-subjects residuals (n < 3); skipping Shapiro-Wilk.") +} else { + resid_for_shapiro <- if (n_resid > 5000L) { + set.seed(1L) + sample(resid_within, 5000L) + } else { + resid_within + } + print(shapiro.test(resid_for_shapiro)) +} +``` + +### Q-Q plot + +```{r qqplot, fig.height = 4} +qqnorm(resid_within) +qqline(resid_within) +``` + +Shapiro-Wilk is significant but is sensitive to large sample size. +QQ plot shows that centre residuals are normally distributed, with some tail heaviness. +ANOVA is robust to violations of normality w/ large sample size. + +Overall, ANOVA can proceed. + +# Mixed ANOVA + +```{r anova} +aov_afex <- aov_ez( + id = "pID", + dv = "MEAN_SCORE", + data = long_data, + between = c("perspective", "temporalDO"), + within = c("TIME", "DOMAIN"), + type = 3 +) + +cat("\n--- ANOVA Table (Type 3, uncorrected) ---\n") +print(nice(aov_afex, correction = "none")) +cat("\n--- ANOVA Table (Type 3, Greenhouse–Geisser correction) ---\n") +print(nice(aov_afex, correction = "GG")) +``` + +Mauchly's test of sphericity is sig for DOMAIN main effect and interaction (except w/ TIME). Use GG correction for these: +- 8 DOMAIN 1.94, 1004.66 1.21 0.63 <.001, p = .529 +## 9 perspective:DOMAIN 1.94, 1004.66 1.21 7.79 *** <.001, p <.001 +## 10 temporalDO:DOMAIN 1.94, 1004.66 1.21 0.76 <.001, p = .466 +## 11 perspective:temporalDO:DOMAIN 1.94, 1004.66 1.21 0.17 <.001, p = .837 + + +The following are significant main effects and interactions: +## 15 perspective:temporalDO:TIME:DOMAIN 2, 1036 0.75 3.11 * <.001 .045 +## 13 perspective:TIME:DOMAIN 2, 1036 0.75 3.58 * <.001 .028 +## 9 perspective:DOMAIN 1.94, 1004.66 1.21 7.79 *** <.001, p <.001 (GG corrected) +## 6 temporalDO:TIME 1, 518 1.86 9.81 ** <.001 .002 +## 7 perspective:temporalDO:TIME 1, 518 1.86 7.91 ** <.001 .005 +## 4 TIME 1, 518 1.86 10.05 ** <.001 .002 + + +# Mauchly and epsilon + +```{r mauchly} +anova_wide <- anova_data %>% + select(pID, perspective, temporalDO, all_of(within_vars)) %>% + filter(if_all(all_of(within_vars), ~ !is.na(.))) +response_matrix <- as.matrix(anova_wide[, within_vars]) +rm_model <- lm(response_matrix ~ perspective * temporalDO, data = anova_wide) +idata <- data.frame( + TIME = factor(rep(c("past", "fut"), each = 3), levels = c("past", "fut")), + DOMAIN = factor(rep(c("pref", "pers", "val"), 2), levels = c("pref", "pers", "val")) +) +rm_anova <- car::Anova(rm_model, idata = idata, idesign = ~ TIME * DOMAIN, type = 3) +rm_summary <- summary(rm_anova, multivariate = FALSE) +if (!is.null(rm_summary$sphericity.tests)) { + cat("\nMauchly's Test of Sphericity:\n") + print(rm_summary$sphericity.tests) +} +if (!is.null(rm_summary$epsilon)) { + cat("\nEpsilon (GG, HF):\n") + print(rm_summary$epsilon) +} +``` + +# Post-hoc comparisons + +## TIME (main effect) + +```{r posthoc-TIME} +emm_TIME <- emmeans(aov_afex, ~ TIME) +print(pairs(emm_TIME, adjust = "bonferroni")) +``` + +Pairwise comparison provide supprot for EHI effect. + +## temporalDO:TIME + +```{r posthoc-temporalDO-TIME} +emm_temporalDO_TIME <- emmeans(aov_afex, ~ TIME | temporalDO) +print(pairs(emm_temporalDO_TIME, adjust = "bonferroni")) +``` + +Contrast significant only for temporal display order of past first, then future. +When grouped by time instead of temporalDO, no contrasts are significant. + +## perspective:temporalDO:TIME + +```{r posthoc-pt-TIME} +emm_pt_TIME <- emmeans(aov_afex, ~ TIME | perspective + temporalDO) +print(pairs(emm_pt_TIME, adjust = "bonferroni")) +``` + +EHI is significant only for self perspective and past first temporal display order. + +When grouped by perspective or temporalDO instead of TIME, no contrasts are significant. + +## perspective:DOMAIN + +```{r posthoc-perspective-DOMAIN} +emm_perspective_DOMAIN <- emmeans(aov_afex, ~ perspective | DOMAIN) +print(pairs(emm_perspective_DOMAIN, adjust = "bonferroni")) +emm_perspective_DOMAIN_domain <- emmeans(aov_afex, ~ DOMAIN | perspective) +print(pairs(emm_perspective_DOMAIN_domain, adjust = "bonferroni")) +``` + +significance is driven by the change from preferences to values in the "other" perspective. + +## perspective:TIME:DOMAIN + +```{r posthoc-pt-DOMAIN} +emm_pt_TIME_domain <- emmeans(aov_afex, ~ TIME | perspective + DOMAIN) +print(pairs(emm_pt_TIME_domain, adjust = "bonferroni")) +``` + +EHI effects present for other-perspective in the preferences domain and for self-perspective in the values domain. +Estimate is higher in the self-perspective than in the other-perspective. + +```{r posthoc-pt-DOMAIN-domain} +emm_pt_domain_domain <- emmeans(aov_afex, ~ DOMAIN | perspective + TIME) +print(pairs(emm_pt_domain_domain, adjust = "bonferroni")) +``` + +Significant contrasts are driven by domain changes from preferences to values in the self vs other perspectives, in the past-oriented questions. +Trends reverse depending on perspective, where values have higher estimates than preferences in the self-perspective, but lower estimates than preferences in the other-perspective. + +## perspective:temporalDO:TIME:DOMAIN + +```{r posthoc-ptt-TIME} +emm_ptt_TIME <- emmeans(aov_afex, ~ TIME | perspective + temporalDO + DOMAIN) +print(pairs(emm_ptt_TIME, adjust = "bonferroni")) +``` +EHI effects are present for three contrasts: +- past - fut 0.2806 0.118 518 2.380 0.0177 (other-perspective, preferences domain, past-first temporal display order) +- past - fut 0.4358 0.138 518 3.156 0.0017 (self-perspective, personality domain, past-first temporal display order) +- past - fut 0.7276 0.141 518 5.169 <0.0001 (self-perspective, values domain, past-first temporal display order) + +A reverse EHI effect is present for 1 contrast: +- past - fut -0.2367 0.118 518 -2.001 0.0459 (self-personality, preferences domain, future-first temporal display order) + +```{r posthoc-ptt-perspective} +emm_ptt_perspective <- emmeans(aov_afex, ~ perspective | temporalDO + TIME + DOMAIN) +print(pairs(emm_ptt_perspective, adjust = "bonferroni")) +``` +1 significant contrast: +- other - self -0.6972 0.314 518 -2.220 0.0268 (values domain, past-oriented questions, past-first temporal display order) + + +not really of theoretical interest but speaks to the perspective:TIME:DOMAIN interaction. + +no significant contrasts when grouped by temporalDO instead of TIME or perspective. + +## Cohen's d (significant contrasts only) + +```{r cohens-d-significant} +d_data <- anova_data %>% + mutate( + past_mean = (past_pref_MEAN + past_pers_MEAN + past_val_MEAN) / 3, + fut_mean = (fut_pref_MEAN + fut_pers_MEAN + fut_val_MEAN) / 3, + pref_mean = (past_pref_MEAN + fut_pref_MEAN) / 2, + pers_mean = (past_pers_MEAN + fut_pers_MEAN) / 2, + val_mean = (past_val_MEAN + fut_val_MEAN) / 2 + ) + +cohens_d_results <- tibble( + contrast = character(), + condition = character(), + d = double() +) + +# TIME main: past vs fut +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "overall", + d = suppressMessages(effectsize::cohens_d(d_data$past_mean, d_data$fut_mean, paired = TRUE)$Cohens_d) + ) + +# temporalDO:TIME: past vs fut for temporalDO = past +d_past_tdo <- d_data %>% filter(temporalDO == "past") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "temporalDO = past", + d = suppressMessages(effectsize::cohens_d(d_past_tdo$past_mean, d_past_tdo$fut_mean, paired = TRUE)$Cohens_d) + ) + +# perspective:temporalDO:TIME: past vs fut for self, past temporalDO +d_self_past <- d_data %>% filter(perspective == "self", temporalDO == "past") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "self, temporalDO = past", + d = suppressMessages(effectsize::cohens_d(d_self_past$past_mean, d_self_past$fut_mean, paired = TRUE)$Cohens_d) + ) + +# perspective:DOMAIN: pref vs val for perspective = other +d_other <- d_data %>% filter(perspective == "other") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - val)", + condition = "perspective = other", + d = suppressMessages(effectsize::cohens_d(d_other$pref_mean, d_other$val_mean, paired = TRUE)$Cohens_d) + ) + +# perspective:TIME:DOMAIN - TIME: other, pref +d_other_pref <- d_data %>% filter(perspective == "other") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "other, pref domain", + d = suppressMessages(effectsize::cohens_d(d_other_pref$past_pref_MEAN, d_other_pref$fut_pref_MEAN, paired = TRUE)$Cohens_d) + ) + +# perspective:TIME:DOMAIN - TIME: self, val +d_self_val <- d_data %>% filter(perspective == "self") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "self, val domain", + d = suppressMessages(effectsize::cohens_d(d_self_val$past_val_MEAN, d_self_val$fut_val_MEAN, paired = TRUE)$Cohens_d) + ) + +# perspective:TIME:DOMAIN - DOMAIN: other, past TIME +d_other_past <- d_data %>% filter(perspective == "other") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - val)", + condition = "other, past TIME", + d = suppressMessages(effectsize::cohens_d(d_other_past$past_pref_MEAN, d_other_past$past_val_MEAN, paired = TRUE)$Cohens_d) + ) + +# perspective:TIME:DOMAIN - DOMAIN: self, past TIME: pref - pers +d_self_past_t <- d_data %>% filter(perspective == "self") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - pers)", + condition = "self, past TIME", + d = suppressMessages(effectsize::cohens_d(d_self_past_t$past_pref_MEAN, d_self_past_t$past_pers_MEAN, paired = TRUE)$Cohens_d) + ) + +# perspective:TIME:DOMAIN - DOMAIN: self, past TIME: pref - val +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - val)", + condition = "self, past TIME", + d = suppressMessages(effectsize::cohens_d(d_self_past_t$past_pref_MEAN, d_self_past_t$past_val_MEAN, paired = TRUE)$Cohens_d) + ) + +# 4-way TIME: self, future temporalDO, pref (reverse EHI) +d_self_fut_pref <- d_data %>% filter(perspective == "self", temporalDO == "future") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut) [reverse]", + condition = "self, future temporalDO, pref domain", + d = suppressMessages(effectsize::cohens_d(d_self_fut_pref$past_pref_MEAN, d_self_fut_pref$fut_pref_MEAN, paired = TRUE)$Cohens_d) + ) + +# 4-way TIME: other, past temporalDO, pref +d_other_past_pref <- d_data %>% filter(perspective == "other", temporalDO == "past") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "other, past temporalDO, pref domain", + d = suppressMessages(effectsize::cohens_d(d_other_past_pref$past_pref_MEAN, d_other_past_pref$fut_pref_MEAN, paired = TRUE)$Cohens_d) + ) + +# 4-way TIME: self, past temporalDO, pers +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "self, past temporalDO, pers domain", + d = suppressMessages(effectsize::cohens_d(d_self_past$past_pers_MEAN, d_self_past$fut_pers_MEAN, paired = TRUE)$Cohens_d) + ) + +# 4-way TIME: self, past temporalDO, val +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "self, past temporalDO, val domain", + d = suppressMessages(effectsize::cohens_d(d_self_past$past_val_MEAN, d_self_past$fut_val_MEAN, paired = TRUE)$Cohens_d) + ) + +# 4-way perspective: past temporalDO, past TIME, val domain (other - self, between-subjects) +d_ptt_val <- d_data %>% + filter(temporalDO == "past") %>% + select(perspective, past_val_MEAN) +d_other_ptt <- d_ptt_val %>% filter(perspective == "other") %>% pull(past_val_MEAN) +d_self_ptt <- d_ptt_val %>% filter(perspective == "self") %>% pull(past_val_MEAN) +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "perspective (other - self)", + condition = "past temporalDO, past TIME, val domain", + d = suppressMessages(effectsize::cohens_d(d_other_ptt, d_self_ptt, paired = FALSE)$Cohens_d) + ) + +cohens_d_results %>% + mutate(d = round(d, 3)) %>% + print(n = Inf) +``` diff --git a/eohi3/knit/DA01_anova_DS.html b/eohi3/knit/DA01_anova_DS.html new file mode 100644 index 0000000..3481599 --- /dev/null +++ b/eohi3/knit/DA01_anova_DS.html @@ -0,0 +1,2535 @@ + + + + + + + + + + + + + + + +Mixed ANOVA - Domain Specific Means (DA01) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + +
+
+
+
+
+ +
+ + + + + + + +
+

Setup

+
library(tidyverse)
+library(rstatix)
+library(emmeans)
+library(effectsize)
+library(afex)
+library(car)
+
+options(scipen = 999)
+afex::set_sum_contrasts()
+
+
+

Data

+
# Data file is in parent of knit/ (eohi3/eohi3.csv)
+df <- read.csv(
+  "/home/ladmin/Documents/DND/EOHI/eohi3/eohi3.csv",
+  stringsAsFactors = FALSE,
+  check.names = FALSE,
+  na.strings = "NA"
+)
+
+between_vars <- c("perspective", "temporalDO")
+within_vars <- c(
+  "past_pref_MEAN", "past_pers_MEAN", "past_val_MEAN",
+  "fut_pref_MEAN",  "fut_pers_MEAN",  "fut_val_MEAN"
+)
+
+missing_vars <- setdiff(c(between_vars, within_vars, "pID"), names(df))
+if (length(missing_vars) > 0) {
+  stop(paste("Missing required variables:", paste(missing_vars, collapse = ", ")))
+}
+
+anova_data <- df %>%
+  select(pID, all_of(between_vars), all_of(within_vars)) %>%
+  filter(
+    !is.na(perspective), perspective != "",
+    !is.na(temporalDO),  temporalDO  != ""
+  )
+
+
+

Long format

+
long_data <- anova_data %>%
+  pivot_longer(
+    cols = all_of(within_vars),
+    names_to = "variable",
+    values_to = "MEAN_SCORE"
+  ) %>%
+  mutate(
+    time = case_when(
+      grepl("^past_", variable) ~ "past",
+      grepl("^fut_",  variable) ~ "fut",
+      TRUE ~ NA_character_
+    ),
+    domain = case_when(
+      grepl("_pref_MEAN$", variable) ~ "pref",
+      grepl("_pers_MEAN$", variable) ~ "pers",
+      grepl("_val_MEAN$",  variable) ~ "val",
+      TRUE ~ NA_character_
+    )
+  ) %>%
+  mutate(
+    TIME       = factor(time,   levels = c("past", "fut")),
+    DOMAIN     = factor(domain, levels = c("pref", "pers", "val")),
+    perspective = factor(perspective),
+    temporalDO  = factor(temporalDO),
+    pID         = factor(pID)
+  ) %>%
+  select(pID, perspective, temporalDO, TIME, DOMAIN, MEAN_SCORE) %>%
+  filter(!is.na(MEAN_SCORE))
+
+
+

Descriptive statistics

+
desc_stats <- long_data %>%
+  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
+  summarise(
+    n        = n(),
+    mean     = round(mean(MEAN_SCORE), 5),
+    variance = round(var(MEAN_SCORE), 5),
+    sd       = round(sd(MEAN_SCORE), 5),
+    median   = round(median(MEAN_SCORE), 5),
+    q1       = round(quantile(MEAN_SCORE, 0.25), 5),
+    q3       = round(quantile(MEAN_SCORE, 0.75), 5),
+    min      = round(min(MEAN_SCORE), 5),
+    max      = round(max(MEAN_SCORE), 5),
+    .groups  = "drop"
+  )
+
+# Show all rows and columns (no truncation)
+options(tibble.width = Inf)
+print(desc_stats, n = Inf)
+
## # A tibble: 24 × 13
+##    perspective temporalDO TIME  DOMAIN     n  mean variance    sd median    q1
+##    <fct>       <fct>      <fct> <fct>  <int> <dbl>    <dbl> <dbl>  <dbl> <dbl>
+##  1 other       future     past  pref     122  3.73     5.75  2.40   3.67  1.67
+##  2 other       future     past  pers     122  3.60     5.65  2.38   3.6   1.4 
+##  3 other       future     past  val      122  3.60     6.27  2.50   3.5   1.33
+##  4 other       future     fut   pref     122  3.65     4.76  2.18   3.67  1.83
+##  5 other       future     fut   pers     122  3.50     4.74  2.18   3.4   1.65
+##  6 other       future     fut   val      122  3.38     5.17  2.27   3     1.5 
+##  7 other       past       past  pref     139  3.95     5.86  2.42   4     2   
+##  8 other       past       past  pers     139  3.82     6.58  2.56   3.6   1.6 
+##  9 other       past       past  val      139  3.60     5.75  2.40   3.33  1.83
+## 10 other       past       fut   pref     139  3.67     6.24  2.50   3.33  1.75
+## 11 other       past       fut   pers     139  3.65     6.86  2.62   3.4   1.5 
+## 12 other       past       fut   val      139  3.56     7.14  2.67   3.17  1.5 
+## 13 self        future     past  pref     138  3.60     5.77  2.40   3.42  1.5 
+## 14 self        future     past  pers     138  3.74     6.00  2.45   3.6   1.65
+## 15 self        future     past  val      138  3.79     6.52  2.55   3.67  1.67
+## 16 self        future     fut   pref     138  3.83     5.36  2.32   3.75  2   
+## 17 self        future     fut   pers     138  3.83     6.24  2.50   3.9   1.4 
+## 18 self        future     fut   val      138  3.86     6.15  2.48   4     1.83
+## 19 self        past       past  pref     123  3.84     5.82  2.41   3.83  1.83
+## 20 self        past       past  pers     123  4.19     6.60  2.57   4.2   1.9 
+## 21 self        past       past  val      123  4.30     7.29  2.70   4.33  1.83
+## 22 self        past       fut   pref     123  3.65     6.76  2.60   3.33  1.33
+## 23 self        past       fut   pers     123  3.75     7.09  2.66   3.4   1.6 
+## 24 self        past       fut   val      123  3.57     7.93  2.82   3.17  1   
+##       q3   min   max
+##    <dbl> <dbl> <dbl>
+##  1  5.5      0 10   
+##  2  5        0 10   
+##  3  5.29     0 10   
+##  4  5.25     0 10   
+##  5  5.2      0  8   
+##  6  5.17     0  8.67
+##  7  5.5      0 10   
+##  8  5.8      0 10   
+##  9  5.33     0 10   
+## 10  5.17     0 10   
+## 11  5.6      0 10   
+## 12  5.25     0 10   
+## 13  5.29     0 10   
+## 14  5.75     0  8.6 
+## 15  5.62     0  9   
+## 16  5.5      0 10   
+## 17  5.95     0  8.6 
+## 18  5.62     0  9.17
+## 19  5.58     0 10   
+## 20  6.2      0 10   
+## 21  6.5      0 10   
+## 22  5.67     0 10   
+## 23  5.5      0 10   
+## 24  5.67     0 10
+

Interpretations:
+1. Mean and median values are similar, indicating distribution is +relatively symmetric and any skew is minimal. Any outliers are not +extreme.
+2. Highest to lowest group n size ratio is 1.14 (139/122). Acceptable +ratio for ANOVA (under 1.5).
+3. Highest to lowest overall group variance ratio is 1.67 (7.93/4.74). +Acceptable ratio for ANOVA (under 4).
+For the sake of consistency w/ the other EHI studies, I also calculated +Hartley’s F-max ratio.
+The conservative F-max critical value is 1.60, which is still higher +than the highest observed F-max ratio of 1.53.

+
+
+

Assumption checks

+
+

Missing values

+
missing_summary <- long_data %>%
+  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
+  summarise(
+    n_total     = n(),
+    n_missing   = sum(is.na(MEAN_SCORE)),
+    pct_missing = round(100 * n_missing / n_total, 2),
+    .groups     = "drop"
+  )
+
+print(missing_summary, n = Inf)
+
## # A tibble: 24 × 7
+##    perspective temporalDO TIME  DOMAIN n_total n_missing pct_missing
+##    <fct>       <fct>      <fct> <fct>    <int>     <int>       <dbl>
+##  1 other       future     past  pref       122         0           0
+##  2 other       future     past  pers       122         0           0
+##  3 other       future     past  val        122         0           0
+##  4 other       future     fut   pref       122         0           0
+##  5 other       future     fut   pers       122         0           0
+##  6 other       future     fut   val        122         0           0
+##  7 other       past       past  pref       139         0           0
+##  8 other       past       past  pers       139         0           0
+##  9 other       past       past  val        139         0           0
+## 10 other       past       fut   pref       139         0           0
+## 11 other       past       fut   pers       139         0           0
+## 12 other       past       fut   val        139         0           0
+## 13 self        future     past  pref       138         0           0
+## 14 self        future     past  pers       138         0           0
+## 15 self        future     past  val        138         0           0
+## 16 self        future     fut   pref       138         0           0
+## 17 self        future     fut   pers       138         0           0
+## 18 self        future     fut   val        138         0           0
+## 19 self        past       past  pref       123         0           0
+## 20 self        past       past  pers       123         0           0
+## 21 self        past       past  val        123         0           0
+## 22 self        past       fut   pref       123         0           0
+## 23 self        past       fut   pers       123         0           0
+## 24 self        past       fut   val        123         0           0
+

No missing values. As expected.

+
+
+

Outliers

+
outlier_summary <- long_data %>%
+  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
+  summarise(
+    n          = n(),
+    n_outliers = sum(abs(scale(MEAN_SCORE)) > 3),
+    .groups    = "drop"
+  )
+
+print(outlier_summary, n = Inf)
+
## # A tibble: 24 × 6
+##    perspective temporalDO TIME  DOMAIN     n n_outliers
+##    <fct>       <fct>      <fct> <fct>  <int>      <int>
+##  1 other       future     past  pref     122          0
+##  2 other       future     past  pers     122          0
+##  3 other       future     past  val      122          0
+##  4 other       future     fut   pref     122          0
+##  5 other       future     fut   pers     122          0
+##  6 other       future     fut   val      122          0
+##  7 other       past       past  pref     139          0
+##  8 other       past       past  pers     139          0
+##  9 other       past       past  val      139          0
+## 10 other       past       fut   pref     139          0
+## 11 other       past       fut   pers     139          0
+## 12 other       past       fut   val      139          0
+## 13 self        future     past  pref     138          0
+## 14 self        future     past  pers     138          0
+## 15 self        future     past  val      138          0
+## 16 self        future     fut   pref     138          0
+## 17 self        future     fut   pers     138          0
+## 18 self        future     fut   val      138          0
+## 19 self        past       past  pref     123          0
+## 20 self        past       past  pers     123          0
+## 21 self        past       past  val      123          0
+## 22 self        past       fut   pref     123          0
+## 23 self        past       fut   pers     123          0
+## 24 self        past       fut   val      123          0
+

No outliers present. Good.

+
+
+

Homogeneity of variance

+
homogeneity_between <- long_data %>%
+  group_by(TIME, DOMAIN) %>%
+  rstatix::levene_test(MEAN_SCORE ~ perspective * temporalDO)
+
+print(homogeneity_between, n = Inf)
+
## # A tibble: 6 × 6
+##   TIME  DOMAIN   df1   df2 statistic      p
+##   <fct> <fct>  <int> <int>     <dbl>  <dbl>
+## 1 past  pref       3   518    0.0857 0.968 
+## 2 past  pers       3   518    0.974  0.405 
+## 3 past  val        3   518    1.81   0.144 
+## 4 fut   pref       3   518    1.68   0.170 
+## 5 fut   pers       3   518    2.82   0.0383
+## 6 fut   val        3   518    2.87   0.0360
+

Levene’s test is sigifnicant for two cells: fut-pers and fut-val. +However, variance ratios and F-max tests show that any +heteroscedasticity is mild.

+
+
+

Normality (within-subjects residuals)

+
resid_within <- long_data %>%
+  group_by(pID) %>%
+  mutate(person_mean = mean(MEAN_SCORE, na.rm = TRUE)) %>%
+  ungroup() %>%
+  mutate(resid = MEAN_SCORE - person_mean) %>%
+  pull(resid)
+resid_within <- resid_within[!is.na(resid_within)]
+
+n_resid <- length(resid_within)
+if (n_resid < 3L) {
+  message("Too few within-subjects residuals (n < 3); skipping Shapiro-Wilk.")
+} else {
+  resid_for_shapiro <- if (n_resid > 5000L) {
+    set.seed(1L)
+    sample(resid_within, 5000L)
+  } else {
+    resid_within
+  }
+  print(shapiro.test(resid_for_shapiro))
+}
+
## 
+##  Shapiro-Wilk normality test
+## 
+## data:  resid_for_shapiro
+## W = 0.9431, p-value < 0.00000000000000022
+
+

Q-Q plot

+
qqnorm(resid_within)
+qqline(resid_within)
+

+

Shapiro-Wilk is significant but is sensitive to large sample size. QQ +plot shows that centre residuals are normally distributed, with some +tail heaviness. ANOVA is robust to violations of normality w/ large +sample size.

+

Overall, ANOVA can proceed.

+
+
+
+
+

Mixed ANOVA

+
aov_afex <- aov_ez(
+  id      = "pID",
+  dv      = "MEAN_SCORE",
+  data    = long_data,
+  between = c("perspective", "temporalDO"),
+  within  = c("TIME", "DOMAIN"),
+  type    = 3
+)
+
+cat("\n--- ANOVA Table (Type 3, uncorrected) ---\n")
+
## 
+## --- ANOVA Table (Type 3, uncorrected) ---
+
print(nice(aov_afex, correction = "none"))
+
## Anova Table (Type 3 tests)
+## 
+## Response: MEAN_SCORE
+##                                Effect      df   MSE        F   ges p.value
+## 1                         perspective  1, 518 31.39     0.88  .001    .349
+## 2                          temporalDO  1, 518 31.39     0.36 <.001    .551
+## 3              perspective:temporalDO  1, 518 31.39     0.00 <.001    .955
+## 4                                TIME  1, 518  1.86 10.05 ** <.001    .002
+## 5                    perspective:TIME  1, 518  1.86     0.02 <.001    .895
+## 6                     temporalDO:TIME  1, 518  1.86  9.81 ** <.001    .002
+## 7         perspective:temporalDO:TIME  1, 518  1.86  7.91 ** <.001    .005
+## 8                              DOMAIN 2, 1036  1.17     0.63 <.001    .534
+## 9                  perspective:DOMAIN 2, 1036  1.17 7.79 *** <.001   <.001
+## 10                  temporalDO:DOMAIN 2, 1036  1.17     0.76 <.001    .470
+## 11      perspective:temporalDO:DOMAIN 2, 1036  1.17     0.17 <.001    .843
+## 12                        TIME:DOMAIN 2, 1036  0.75     2.00 <.001    .136
+## 13            perspective:TIME:DOMAIN 2, 1036  0.75   3.58 * <.001    .028
+## 14             temporalDO:TIME:DOMAIN 2, 1036  0.75     0.01 <.001    .987
+## 15 perspective:temporalDO:TIME:DOMAIN 2, 1036  0.75   3.11 * <.001    .045
+## ---
+## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
+
cat("\n--- ANOVA Table (Type 3, Greenhouse–Geisser correction) ---\n")
+
## 
+## --- ANOVA Table (Type 3, Greenhouse–Geisser correction) ---
+
print(nice(aov_afex, correction = "GG"))
+
## Anova Table (Type 3 tests)
+## 
+## Response: MEAN_SCORE
+##                                Effect            df   MSE        F   ges
+## 1                         perspective        1, 518 31.39     0.88  .001
+## 2                          temporalDO        1, 518 31.39     0.36 <.001
+## 3              perspective:temporalDO        1, 518 31.39     0.00 <.001
+## 4                                TIME        1, 518  1.86 10.05 ** <.001
+## 5                    perspective:TIME        1, 518  1.86     0.02 <.001
+## 6                     temporalDO:TIME        1, 518  1.86  9.81 ** <.001
+## 7         perspective:temporalDO:TIME        1, 518  1.86  7.91 ** <.001
+## 8                              DOMAIN 1.94, 1004.66  1.21     0.63 <.001
+## 9                  perspective:DOMAIN 1.94, 1004.66  1.21 7.79 *** <.001
+## 10                  temporalDO:DOMAIN 1.94, 1004.66  1.21     0.76 <.001
+## 11      perspective:temporalDO:DOMAIN 1.94, 1004.66  1.21     0.17 <.001
+## 12                        TIME:DOMAIN 2.00, 1034.51  0.75     2.00 <.001
+## 13            perspective:TIME:DOMAIN 2.00, 1034.51  0.75   3.58 * <.001
+## 14             temporalDO:TIME:DOMAIN 2.00, 1034.51  0.75     0.01 <.001
+## 15 perspective:temporalDO:TIME:DOMAIN 2.00, 1034.51  0.75   3.11 * <.001
+##    p.value
+## 1     .349
+## 2     .551
+## 3     .955
+## 4     .002
+## 5     .895
+## 6     .002
+## 7     .005
+## 8     .529
+## 9    <.001
+## 10    .466
+## 11    .837
+## 12    .136
+## 13    .028
+## 14    .987
+## 15    .045
+## ---
+## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
+## 
+## Sphericity correction method: GG
+

Mauchly’s test of sphericity is sig for DOMAIN main effect and +interaction (except w/ TIME). Use GG correction for these:
+- 8 DOMAIN 1.94, 1004.66 1.21 0.63 <.001, p = .529
+## 9 perspective:DOMAIN 1.94, 1004.66 1.21 7.79 *** <.001, p +<.001
+## 10 temporalDO:DOMAIN 1.94, 1004.66 1.21 0.76 <.001, p = .466
+## 11 perspective:temporalDO:DOMAIN 1.94, 1004.66 1.21 0.17 <.001, p += .837

+

The following are significant main effects and interactions:
+## 15 perspective:temporalDO:TIME:DOMAIN 2, 1036 0.75 3.11 * <.001 +.045
+## 13 perspective:TIME:DOMAIN 2, 1036 0.75 3.58 * <.001 .028
+## 9 perspective:DOMAIN 1.94, 1004.66 1.21 7.79 *** <.001, p <.001 +(GG corrected)
+## 6 temporalDO:TIME 1, 518 1.86 9.81 ** <.001 .002
+## 7 perspective:temporalDO:TIME 1, 518 1.86 7.91 ** <.001 .005
+## 4 TIME 1, 518 1.86 10.05 ** <.001 .002

+
+
+

Mauchly and epsilon

+
anova_wide <- anova_data %>%
+  select(pID, perspective, temporalDO, all_of(within_vars)) %>%
+  filter(if_all(all_of(within_vars), ~ !is.na(.)))
+response_matrix <- as.matrix(anova_wide[, within_vars])
+rm_model <- lm(response_matrix ~ perspective * temporalDO, data = anova_wide)
+idata <- data.frame(
+  TIME   = factor(rep(c("past", "fut"), each = 3), levels = c("past", "fut")),
+  DOMAIN = factor(rep(c("pref", "pers", "val"), 2), levels = c("pref", "pers", "val"))
+)
+rm_anova <- car::Anova(rm_model, idata = idata, idesign = ~ TIME * DOMAIN, type = 3)
+rm_summary <- summary(rm_anova, multivariate = FALSE)
+
## Warning in summary.Anova.mlm(rm_anova, multivariate = FALSE): HF eps > 1
+## treated as 1
+
if (!is.null(rm_summary$sphericity.tests)) {
+  cat("\nMauchly's Test of Sphericity:\n")
+  print(rm_summary$sphericity.tests)
+}
+
## 
+## Mauchly's Test of Sphericity:
+##                                    Test statistic p-value
+## DOMAIN                                    0.96880 0.00028
+## perspective:DOMAIN                        0.96880 0.00028
+## temporalDO:DOMAIN                         0.96880 0.00028
+## perspective:temporalDO:DOMAIN             0.96880 0.00028
+## TIME:DOMAIN                               0.99856 0.68858
+## perspective:TIME:DOMAIN                   0.99856 0.68858
+## temporalDO:TIME:DOMAIN                    0.99856 0.68858
+## perspective:temporalDO:TIME:DOMAIN        0.99856 0.68858
+
if (!is.null(rm_summary$epsilon)) {
+  cat("\nEpsilon (GG, HF):\n")
+  print(rm_summary$epsilon)
+}
+
+
+

Post-hoc comparisons

+
+

TIME (main effect)

+
emm_TIME <- emmeans(aov_afex, ~ TIME)
+print(pairs(emm_TIME, adjust = "bonferroni"))
+
##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut    0.155 0.0488 518   3.170  0.0016
+## 
+## Results are averaged over the levels of: perspective, temporalDO, DOMAIN
+

Pairwise comparison provide supprot for EHI effect.

+
+
+

temporalDO:TIME

+
emm_temporalDO_TIME <- emmeans(aov_afex, ~ TIME | temporalDO)
+print(pairs(emm_temporalDO_TIME, adjust = "bonferroni"))
+
## temporalDO = future:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut  0.00186 0.0692 518   0.027  0.9786
+## 
+## temporalDO = past:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut  0.30772 0.0689 518   4.465 <0.0001
+## 
+## Results are averaged over the levels of: perspective, DOMAIN
+

Contrast significant only for temporal display order of past first, +then future. When grouped by time instead of temporalDO, no contrasts +are significant.

+
+
+

perspective:temporalDO:TIME

+
emm_pt_TIME <- emmeans(aov_afex, ~ TIME | perspective + temporalDO)
+print(pairs(emm_pt_TIME, adjust = "bonferroni"))
+
## perspective = other, temporalDO = future:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut    0.133 0.1010 518   1.317  0.1883
+## 
+## perspective = self, temporalDO = future:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut   -0.129 0.0948 518  -1.362  0.1739
+## 
+## perspective = other, temporalDO = past:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut    0.164 0.0944 518   1.736  0.0831
+## 
+## perspective = self, temporalDO = past:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut    0.451 0.1000 518   4.497 <0.0001
+## 
+## Results are averaged over the levels of: DOMAIN
+

EHI is significant only for self perspective and past first temporal +display order.

+

When grouped by perspective or temporalDO instead of TIME, no +contrasts are significant.

+
+
+

perspective:DOMAIN

+
emm_perspective_DOMAIN <- emmeans(aov_afex, ~ perspective | DOMAIN)
+print(pairs(emm_perspective_DOMAIN, adjust = "bonferroni"))
+
## DOMAIN = pref:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self   0.0191 0.202 518   0.094  0.9248
+## 
+## DOMAIN = pers:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.2379 0.208 518  -1.141  0.2543
+## 
+## DOMAIN = val:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.3459 0.213 518  -1.622  0.1053
+## 
+## Results are averaged over the levels of: temporalDO, TIME
+
emm_perspective_DOMAIN_domain <- emmeans(aov_afex, ~ DOMAIN | perspective)
+print(pairs(emm_perspective_DOMAIN_domain, adjust = "bonferroni"))
+
## perspective = other:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers  0.11064 0.0696 518   1.591  0.3369
+##  pref - val   0.21702 0.0706 518   3.074  0.0067
+##  pers - val   0.10638 0.0610 518   1.744  0.2454
+## 
+## perspective = self:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers -0.14631 0.0695 518  -2.104  0.1075
+##  pref - val  -0.14800 0.0706 518  -2.097  0.1093
+##  pers - val  -0.00169 0.0610 518  -0.028  1.0000
+## 
+## Results are averaged over the levels of: temporalDO, TIME 
+## P value adjustment: bonferroni method for 3 tests
+

significance is driven by the change from preferences to values in +the “other” perspective.

+
+
+

perspective:TIME:DOMAIN

+
emm_pt_TIME_domain <- emmeans(aov_afex, ~ TIME | perspective + DOMAIN)
+print(pairs(emm_pt_TIME_domain, adjust = "bonferroni"))
+
## perspective = other, DOMAIN = pref:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut   0.1799 0.0862 518   2.087  0.0374
+## 
+## perspective = self, DOMAIN = pref:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut  -0.0228 0.0862 518  -0.265  0.7912
+## 
+## perspective = other, DOMAIN = pers:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut   0.1351 0.0950 518   1.422  0.1555
+## 
+## perspective = self, DOMAIN = pers:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut   0.1759 0.0949 518   1.852  0.0645
+## 
+## perspective = other, DOMAIN = val:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut   0.1301 0.0968 518   1.344  0.1797
+## 
+## perspective = self, DOMAIN = val:
+##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut   0.3306 0.0968 518   3.416  0.0007
+## 
+## Results are averaged over the levels of: temporalDO
+

EHI effects present for other-perspective in the preferences domain +and for self-perspective in the values domain.
+Estimate is higher in the self-perspective than in the +other-perspective.

+
emm_pt_domain_domain <- emmeans(aov_afex, ~ DOMAIN | perspective + TIME)
+print(pairs(emm_pt_domain_domain, adjust = "bonferroni"))
+
## perspective = other, TIME = past:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers   0.1330 0.0894 518   1.488  0.4122
+##  pref - val    0.2419 0.0936 518   2.584  0.0301
+##  pers - val    0.1089 0.0791 518   1.377  0.5072
+## 
+## perspective = self, TIME = past:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers  -0.2457 0.0894 518  -2.749  0.0186
+##  pref - val   -0.3247 0.0936 518  -3.470  0.0017
+##  pers - val   -0.0791 0.0790 518  -1.001  0.9526
+## 
+## perspective = other, TIME = fut:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers   0.0882 0.0850 518   1.038  0.8993
+##  pref - val    0.1921 0.0843 518   2.280  0.0690
+##  pers - val    0.1039 0.0839 518   1.239  0.6480
+## 
+## perspective = self, TIME = fut:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers  -0.0470 0.0850 518  -0.553  1.0000
+##  pref - val    0.0287 0.0842 518   0.341  1.0000
+##  pers - val    0.0757 0.0838 518   0.903  1.0000
+## 
+## Results are averaged over the levels of: temporalDO 
+## P value adjustment: bonferroni method for 3 tests
+

Significant contrasts are driven by domain changes from preferences +to values in the self vs other perspectives, in the past-oriented +questions.
+Trends reverse depending on perspective, where values have higher +estimates than preferences in the self-perspective, but lower estimates +than preferences in the other-perspective.

+
+
+

perspective:temporalDO:TIME:DOMAIN

+
emm_ptt_TIME <- emmeans(aov_afex, ~ TIME | perspective + temporalDO + DOMAIN)
+print(pairs(emm_ptt_TIME, adjust = "bonferroni"))
+
## perspective = other, temporalDO = future, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.0792 0.126 518   0.630  0.5292
+## 
+## perspective = self, temporalDO = future, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut  -0.2367 0.118 518  -2.001  0.0459
+## 
+## perspective = other, temporalDO = past, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.2806 0.118 518   2.380  0.0177
+## 
+## perspective = self, temporalDO = past, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.1911 0.125 518   1.525  0.1280
+## 
+## perspective = other, temporalDO = future, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.1033 0.139 518   0.745  0.4566
+## 
+## perspective = self, temporalDO = future, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut  -0.0841 0.130 518  -0.645  0.5193
+## 
+## perspective = other, temporalDO = past, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.1669 0.130 518   1.285  0.1993
+## 
+## perspective = self, temporalDO = past, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.4358 0.138 518   3.156  0.0017
+## 
+## perspective = other, temporalDO = future, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.2158 0.141 518   1.527  0.1273
+## 
+## perspective = self, temporalDO = future, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut  -0.0664 0.133 518  -0.500  0.6174
+## 
+## perspective = other, temporalDO = past, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.0444 0.132 518   0.335  0.7377
+## 
+## perspective = self, temporalDO = past, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.7276 0.141 518   5.169 <0.0001
+

EHI effects are present for three contrasts:
+- past - fut 0.2806 0.118 518 2.380 0.0177 (other-perspective, +preferences domain, past-first temporal display order)
+- past - fut 0.4358 0.138 518 3.156 0.0017 (self-perspective, +personality domain, past-first temporal display order)
+- past - fut 0.7276 0.141 518 5.169 <0.0001 (self-perspective, values +domain, past-first temporal display order)

+

A reverse EHI effect is present for 1 contrast:
+- past - fut -0.2367 0.118 518 -2.001 0.0459 (self-personality, +preferences domain, future-first temporal display order)

+
emm_ptt_perspective <- emmeans(aov_afex, ~ perspective | temporalDO + TIME + DOMAIN)
+print(pairs(emm_ptt_perspective, adjust = "bonferroni"))
+
## temporalDO = future, TIME = past, DOMAIN = pref:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self   0.1343 0.299 518   0.449  0.6540
+## 
+## temporalDO = past, TIME = past, DOMAIN = pref:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self   0.1067 0.298 518   0.358  0.7207
+## 
+## temporalDO = future, TIME = fut, DOMAIN = pref:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.1817 0.299 518  -0.608  0.5436
+## 
+## temporalDO = past, TIME = fut, DOMAIN = pref:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self   0.0171 0.298 518   0.058  0.9541
+## 
+## temporalDO = future, TIME = past, DOMAIN = pers:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.1451 0.310 518  -0.469  0.6396
+## 
+## temporalDO = past, TIME = past, DOMAIN = pers:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.3714 0.309 518  -1.204  0.2292
+## 
+## temporalDO = future, TIME = fut, DOMAIN = pers:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.3325 0.311 518  -1.070  0.2852
+## 
+## temporalDO = past, TIME = fut, DOMAIN = pers:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.1025 0.310 518  -0.331  0.7407
+## 
+## temporalDO = future, TIME = past, DOMAIN = val:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.1951 0.315 518  -0.619  0.5362
+## 
+## temporalDO = past, TIME = past, DOMAIN = val:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.6972 0.314 518  -2.220  0.0268
+## 
+## temporalDO = future, TIME = fut, DOMAIN = val:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.4774 0.319 518  -1.495  0.1355
+## 
+## temporalDO = past, TIME = fut, DOMAIN = val:
+##  contrast     estimate    SE  df t.ratio p.value
+##  other - self  -0.0139 0.318 518  -0.044  0.9650
+

1 significant contrast:
+- other - self -0.6972 0.314 518 -2.220 0.0268 (values domain, +past-oriented questions, past-first temporal display order)

+

not really of theoretical interest but speaks to the +perspective:TIME:DOMAIN interaction.

+

no significant contrasts when grouped by temporalDO instead of TIME +or perspective.

+
+
+

Cohen’s d (significant contrasts only)

+
d_data <- anova_data %>%
+  mutate(
+    past_mean = (past_pref_MEAN + past_pers_MEAN + past_val_MEAN) / 3,
+    fut_mean  = (fut_pref_MEAN + fut_pers_MEAN + fut_val_MEAN) / 3,
+    pref_mean = (past_pref_MEAN + fut_pref_MEAN) / 2,
+    pers_mean = (past_pers_MEAN + fut_pers_MEAN) / 2,
+    val_mean  = (past_val_MEAN + fut_val_MEAN) / 2
+  )
+
+cohens_d_results <- tibble(
+  contrast = character(),
+  condition = character(),
+  d = double()
+)
+
+# TIME main: past vs fut
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "overall",
+    d = suppressMessages(effectsize::cohens_d(d_data$past_mean, d_data$fut_mean, paired = TRUE)$Cohens_d)
+  )
+
+# temporalDO:TIME: past vs fut for temporalDO = past
+d_past_tdo <- d_data %>% filter(temporalDO == "past")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "temporalDO = past",
+    d = suppressMessages(effectsize::cohens_d(d_past_tdo$past_mean, d_past_tdo$fut_mean, paired = TRUE)$Cohens_d)
+  )
+
+# perspective:temporalDO:TIME: past vs fut for self, past temporalDO
+d_self_past <- d_data %>% filter(perspective == "self", temporalDO == "past")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "self, temporalDO = past",
+    d = suppressMessages(effectsize::cohens_d(d_self_past$past_mean, d_self_past$fut_mean, paired = TRUE)$Cohens_d)
+  )
+
+# perspective:DOMAIN: pref vs val for perspective = other
+d_other <- d_data %>% filter(perspective == "other")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - val)",
+    condition = "perspective = other",
+    d = suppressMessages(effectsize::cohens_d(d_other$pref_mean, d_other$val_mean, paired = TRUE)$Cohens_d)
+  )
+
+# perspective:TIME:DOMAIN - TIME: other, pref
+d_other_pref <- d_data %>% filter(perspective == "other")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "other, pref domain",
+    d = suppressMessages(effectsize::cohens_d(d_other_pref$past_pref_MEAN, d_other_pref$fut_pref_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# perspective:TIME:DOMAIN - TIME: self, val
+d_self_val <- d_data %>% filter(perspective == "self")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "self, val domain",
+    d = suppressMessages(effectsize::cohens_d(d_self_val$past_val_MEAN, d_self_val$fut_val_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# perspective:TIME:DOMAIN - DOMAIN: other, past TIME
+d_other_past <- d_data %>% filter(perspective == "other")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - val)",
+    condition = "other, past TIME",
+    d = suppressMessages(effectsize::cohens_d(d_other_past$past_pref_MEAN, d_other_past$past_val_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# perspective:TIME:DOMAIN - DOMAIN: self, past TIME: pref - pers
+d_self_past_t <- d_data %>% filter(perspective == "self")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - pers)",
+    condition = "self, past TIME",
+    d = suppressMessages(effectsize::cohens_d(d_self_past_t$past_pref_MEAN, d_self_past_t$past_pers_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# perspective:TIME:DOMAIN - DOMAIN: self, past TIME: pref - val
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - val)",
+    condition = "self, past TIME",
+    d = suppressMessages(effectsize::cohens_d(d_self_past_t$past_pref_MEAN, d_self_past_t$past_val_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way TIME: self, future temporalDO, pref (reverse EHI)
+d_self_fut_pref <- d_data %>% filter(perspective == "self", temporalDO == "future")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut) [reverse]",
+    condition = "self, future temporalDO, pref domain",
+    d = suppressMessages(effectsize::cohens_d(d_self_fut_pref$past_pref_MEAN, d_self_fut_pref$fut_pref_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way TIME: other, past temporalDO, pref
+d_other_past_pref <- d_data %>% filter(perspective == "other", temporalDO == "past")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "other, past temporalDO, pref domain",
+    d = suppressMessages(effectsize::cohens_d(d_other_past_pref$past_pref_MEAN, d_other_past_pref$fut_pref_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way TIME: self, past temporalDO, pers
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "self, past temporalDO, pers domain",
+    d = suppressMessages(effectsize::cohens_d(d_self_past$past_pers_MEAN, d_self_past$fut_pers_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way TIME: self, past temporalDO, val
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "self, past temporalDO, val domain",
+    d = suppressMessages(effectsize::cohens_d(d_self_past$past_val_MEAN, d_self_past$fut_val_MEAN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way perspective: past temporalDO, past TIME, val domain (other - self, between-subjects)
+d_ptt_val <- d_data %>%
+  filter(temporalDO == "past") %>%
+  select(perspective, past_val_MEAN)
+d_other_ptt <- d_ptt_val %>% filter(perspective == "other") %>% pull(past_val_MEAN)
+d_self_ptt <- d_ptt_val %>% filter(perspective == "self") %>% pull(past_val_MEAN)
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "perspective (other - self)",
+    condition = "past temporalDO, past TIME, val domain",
+    d = suppressMessages(effectsize::cohens_d(d_other_ptt, d_self_ptt, paired = FALSE)$Cohens_d)
+  )
+
+cohens_d_results %>%
+  mutate(d = round(d, 3)) %>%
+  print(n = Inf)
+
## # A tibble: 14 × 3
+##    contrast                    condition                                   d
+##    <chr>                       <chr>                                   <dbl>
+##  1 TIME (past - fut)           overall                                 0.13 
+##  2 TIME (past - fut)           temporalDO = past                       0.249
+##  3 TIME (past - fut)           self, temporalDO = past                 0.413
+##  4 DOMAIN (pref - val)         perspective = other                     0.185
+##  5 TIME (past - fut)           other, pref domain                      0.126
+##  6 TIME (past - fut)           self, val domain                        0.189
+##  7 DOMAIN (pref - val)         other, past TIME                        0.159
+##  8 DOMAIN (pref - pers)        self, past TIME                        -0.171
+##  9 DOMAIN (pref - val)         self, past TIME                        -0.218
+## 10 TIME (past - fut) [reverse] self, future temporalDO, pref domain   -0.174
+## 11 TIME (past - fut)           other, past temporalDO, pref domain     0.172
+## 12 TIME (past - fut)           self, past temporalDO, pers domain      0.309
+## 13 TIME (past - fut)           self, past temporalDO, val domain       0.435
+## 14 perspective (other - self)  past temporalDO, past TIME, val domain -0.274
+
+
+ + + +
+
+ +
+ + + + + + + + + + + + + + + + + diff --git a/eohi3/knit/DA02_anova_DGEN.html b/eohi3/knit/DA02_anova_DGEN.html new file mode 100644 index 0000000..8343800 --- /dev/null +++ b/eohi3/knit/DA02_anova_DGEN.html @@ -0,0 +1,2345 @@ + + + + + + + + + + + + + + + +Mixed ANOVA - Domain General Vars + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + +
+
+
+
+
+ +
+ + + + + + + +
+

Setup

+
library(tidyverse)
+library(rstatix)
+library(emmeans)
+library(effectsize)
+library(afex)
+library(car)
+
+options(scipen = 999)
+afex::set_sum_contrasts()
+
+
+

Data

+
# Data file is in parent of knit/ (eohi3/eohi3.csv)
+df <- read.csv(
+  "/home/ladmin/Documents/DND/EOHI/eohi3/eohi3.csv",
+  stringsAsFactors = FALSE,
+  check.names = FALSE,
+  na.strings = "NA"
+)
+
+between_vars <- c("perspective", "temporalDO")
+within_vars <- c(
+  "past_pref_DGEN", "past_pers_DGEN", "past_val_DGEN",
+  "fut_pref_DGEN",  "fut_pers_DGEN",  "fut_val_DGEN"
+)
+
+missing_vars <- setdiff(c(between_vars, within_vars, "pID"), names(df))
+if (length(missing_vars) > 0) {
+  stop(paste("Missing required variables:", paste(missing_vars, collapse = ", ")))
+}
+
+anova_data <- df %>%
+  select(pID, all_of(between_vars), all_of(within_vars)) %>%
+  filter(
+    !is.na(perspective), perspective != "",
+    !is.na(temporalDO),  temporalDO  != ""
+  )
+
+
+

Long format

+
long_data <- anova_data %>%
+  pivot_longer(
+    cols = all_of(within_vars),
+    names_to = "variable",
+    values_to = "DGEN_SCORE"
+  ) %>%
+  mutate(
+    time = case_when(
+      grepl("^past_", variable) ~ "past",
+      grepl("^fut_",  variable) ~ "fut",
+      TRUE ~ NA_character_
+    ),
+    domain = case_when(
+      grepl("_pref_DGEN$", variable) ~ "pref",
+      grepl("_pers_DGEN$", variable) ~ "pers",
+      grepl("_val_DGEN$",  variable) ~ "val",
+      TRUE ~ NA_character_
+    )
+  ) %>%
+  mutate(
+    TIME       = factor(time,   levels = c("past", "fut")),
+    DOMAIN     = factor(domain, levels = c("pref", "pers", "val")),
+    perspective = factor(perspective),
+    temporalDO  = factor(temporalDO),
+    pID         = factor(pID)
+  ) %>%
+  select(pID, perspective, temporalDO, TIME, DOMAIN, DGEN_SCORE) %>%
+  filter(!is.na(DGEN_SCORE))
+
+
+

Descriptive statistics

+
desc_stats <- long_data %>%
+  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
+  summarise(
+    n        = n(),
+    mean     = round(mean(DGEN_SCORE), 5),
+    variance = round(var(DGEN_SCORE), 5),
+    sd       = round(sd(DGEN_SCORE), 5),
+    median   = round(median(DGEN_SCORE), 5),
+    q1       = round(quantile(DGEN_SCORE, 0.25), 5),
+    q3       = round(quantile(DGEN_SCORE, 0.75), 5),
+    min      = round(min(DGEN_SCORE), 5),
+    max      = round(max(DGEN_SCORE), 5),
+    .groups  = "drop"
+  )
+
+# Show all rows and columns (no truncation)
+options(tibble.width = Inf)
+print(desc_stats, n = Inf)
+
## # A tibble: 24 × 13
+##    perspective temporalDO TIME  DOMAIN     n  mean variance    sd median    q1
+##    <fct>       <fct>      <fct> <fct>  <int> <dbl>    <dbl> <dbl>  <dbl> <dbl>
+##  1 other       future     past  pref     122  3.74     7.19  2.68    4       1
+##  2 other       future     past  pers     122  3.49     7.56  2.75    3       1
+##  3 other       future     past  val      122  3.72     7.64  2.76    3       1
+##  4 other       future     fut   pref     122  3.69     6.65  2.58    3       2
+##  5 other       future     fut   pers     122  3.29     7.03  2.65    3       1
+##  6 other       future     fut   val      122  3.43     7.59  2.75    3       1
+##  7 other       past       past  pref     139  3.93     7.56  2.75    4       2
+##  8 other       past       past  pers     139  3.79     8.12  2.85    3       1
+##  9 other       past       past  val      139  3.30     7.26  2.69    3       1
+## 10 other       past       fut   pref     139  3.60     7.53  2.74    3       1
+## 11 other       past       fut   pers     139  3.50     8.31  2.88    3       1
+## 12 other       past       fut   val      139  3.42     9.24  3.04    3       1
+## 13 self        future     past  pref     138  3.88     7.39  2.72    3.5     2
+## 14 self        future     past  pers     138  3.73     7.48  2.74    3       1
+## 15 self        future     past  val      138  3.65     7.75  2.78    3       1
+## 16 self        future     fut   pref     138  3.92     6.86  2.62    4       2
+## 17 self        future     fut   pers     138  3.72     7.66  2.77    3       1
+## 18 self        future     fut   val      138  3.73     8.04  2.83    4       1
+## 19 self        past       past  pref     123  3.96     7.51  2.74    4       2
+## 20 self        past       past  pers     123  4.24     8.30  2.88    4       2
+## 21 self        past       past  val      123  3.88     9.32  3.05    4       1
+## 22 self        past       fut   pref     123  3.77     7.87  2.80    4       1
+## 23 self        past       fut   pers     123  3.71     8.55  2.92    3       1
+## 24 self        past       fut   val      123  3.34     8.19  2.86    2       1
+##       q3   min   max
+##    <dbl> <dbl> <dbl>
+##  1  6        0    10
+##  2  5.75     0    10
+##  3  5.75     0    10
+##  4  5        0    10
+##  5  5        0    10
+##  6  5        0    10
+##  7  6        0    10
+##  8  6        0    10
+##  9  5        0    10
+## 10  5.5      0    10
+## 11  5        0    10
+## 12  5        0    10
+## 13  6        0    10
+## 14  6        0    10
+## 15  6        0    10
+## 16  6        0    10
+## 17  6        0    10
+## 18  6        0    10
+## 19  6        0    10
+## 20  6        0    10
+## 21  6        0    10
+## 22  6        0    10
+## 23  6        0    10
+## 24  6        0    10
+

Interpretation:
+1. Mean and median values are similar w/ slightly more variation than in +the domain specific anova.
+2. Highest to lowest group n size ratio is 1.14 (139/122). Acceptable +ratio for ANOVA (under 1.5).
+3. Highest to lowest overall group variance ratio is 1.40 (9.32/6.65). +Acceptable ratio for ANOVA (under 4).
+For the sake of consistency w/ the other EHI studies, I also calculated +Hartley’s F-max ratio.
+The conservative F-max critical value is 1.60 (same as DS anova since +number of groups and n sizes doesn’t change), which is still higher than +the highest observed F-max ratio of 1.28.

+
+
+

Assumption checks

+
+

Missing values

+
missing_summary <- long_data %>%
+  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
+  summarise(
+    n_total     = n(),
+    n_missing   = sum(is.na(DGEN_SCORE)),
+    pct_missing = round(100 * n_missing / n_total, 2),
+    .groups     = "drop"
+  )
+
+print(missing_summary, n = Inf)
+
## # A tibble: 24 × 7
+##    perspective temporalDO TIME  DOMAIN n_total n_missing pct_missing
+##    <fct>       <fct>      <fct> <fct>    <int>     <int>       <dbl>
+##  1 other       future     past  pref       122         0           0
+##  2 other       future     past  pers       122         0           0
+##  3 other       future     past  val        122         0           0
+##  4 other       future     fut   pref       122         0           0
+##  5 other       future     fut   pers       122         0           0
+##  6 other       future     fut   val        122         0           0
+##  7 other       past       past  pref       139         0           0
+##  8 other       past       past  pers       139         0           0
+##  9 other       past       past  val        139         0           0
+## 10 other       past       fut   pref       139         0           0
+## 11 other       past       fut   pers       139         0           0
+## 12 other       past       fut   val        139         0           0
+## 13 self        future     past  pref       138         0           0
+## 14 self        future     past  pers       138         0           0
+## 15 self        future     past  val        138         0           0
+## 16 self        future     fut   pref       138         0           0
+## 17 self        future     fut   pers       138         0           0
+## 18 self        future     fut   val        138         0           0
+## 19 self        past       past  pref       123         0           0
+## 20 self        past       past  pers       123         0           0
+## 21 self        past       past  val        123         0           0
+## 22 self        past       fut   pref       123         0           0
+## 23 self        past       fut   pers       123         0           0
+## 24 self        past       fut   val        123         0           0
+

No missing values. As expected.

+
+
+

Outliers

+
outlier_summary <- long_data %>%
+  group_by(perspective, temporalDO, TIME, DOMAIN) %>%
+  summarise(
+    n          = n(),
+    n_outliers = sum(abs(scale(DGEN_SCORE)) > 3),
+    .groups    = "drop"
+  )
+
+print(outlier_summary, n = Inf)
+
## # A tibble: 24 × 6
+##    perspective temporalDO TIME  DOMAIN     n n_outliers
+##    <fct>       <fct>      <fct> <fct>  <int>      <int>
+##  1 other       future     past  pref     122          0
+##  2 other       future     past  pers     122          0
+##  3 other       future     past  val      122          0
+##  4 other       future     fut   pref     122          0
+##  5 other       future     fut   pers     122          0
+##  6 other       future     fut   val      122          0
+##  7 other       past       past  pref     139          0
+##  8 other       past       past  pers     139          0
+##  9 other       past       past  val      139          0
+## 10 other       past       fut   pref     139          0
+## 11 other       past       fut   pers     139          0
+## 12 other       past       fut   val      139          0
+## 13 self        future     past  pref     138          0
+## 14 self        future     past  pers     138          0
+## 15 self        future     past  val      138          0
+## 16 self        future     fut   pref     138          0
+## 17 self        future     fut   pers     138          0
+## 18 self        future     fut   val      138          0
+## 19 self        past       past  pref     123          0
+## 20 self        past       past  pers     123          0
+## 21 self        past       past  val      123          0
+## 22 self        past       fut   pref     123          0
+## 23 self        past       fut   pers     123          0
+## 24 self        past       fut   val      123          0
+

No outliers present. Good. Same as DS anova.

+
+
+

Homogeneity of variance

+
homogeneity_between <- long_data %>%
+  group_by(TIME, DOMAIN) %>%
+  rstatix::levene_test(DGEN_SCORE ~ perspective * temporalDO)
+
+print(homogeneity_between, n = Inf)
+
## # A tibble: 6 × 6
+##   TIME  DOMAIN   df1   df2 statistic      p
+##   <fct> <fct>  <int> <int>     <dbl>  <dbl>
+## 1 past  pref       3   518    0.0902 0.965 
+## 2 past  pers       3   518    0.407  0.748 
+## 3 past  val        3   518    2.69   0.0460
+## 4 fut   pref       3   518    0.824  0.481 
+## 5 fut   pers       3   518    0.876  0.454 
+## 6 fut   val        3   518    0.238  0.870
+

Levene’s test is signiicant for 1 cell only: past-val. However, +variance ratios and F-max tests show that any heteroscedasticity is +mild.

+
+
+

Normality (within-subjects residuals)

+
resid_within <- long_data %>%
+  group_by(pID) %>%
+  mutate(person_mean = mean(DGEN_SCORE, na.rm = TRUE)) %>%
+  ungroup() %>%
+  mutate(resid = DGEN_SCORE - person_mean) %>%
+  pull(resid)
+resid_within <- resid_within[!is.na(resid_within)]
+
+n_resid <- length(resid_within)
+if (n_resid < 3L) {
+  message("Too few within-subjects residuals (n < 3); skipping Shapiro-Wilk.")
+} else {
+  resid_for_shapiro <- if (n_resid > 5000L) {
+    set.seed(1L)
+    sample(resid_within, 5000L)
+  } else {
+    resid_within
+  }
+  print(shapiro.test(resid_for_shapiro))
+}
+
## 
+##  Shapiro-Wilk normality test
+## 
+## data:  resid_for_shapiro
+## W = 0.94192, p-value < 0.00000000000000022
+
+

Q-Q plot

+
qqnorm(resid_within)
+qqline(resid_within)
+

+

Shapiro-Wilk test is significant but is sensitive to large sample +size.
+QQ plot shows that strict centre residuals are normally distributed, but +there is some deviation from normality. ANOVA is robust to violations of +normality w/ large sample size.

+

Overall, ANOVA can proceed.

+
+
+
+
+

Mixed ANOVA

+
aov_afex <- aov_ez(
+  id      = "pID",
+  dv      = "DGEN_SCORE",
+  data    = long_data,
+  between = c("perspective", "temporalDO"),
+  within  = c("TIME", "DOMAIN"),
+  type    = 3,
+  anova_table = list(correction = "none")
+)
+
+print(aov_afex)
+
## Anova Table (Type 3 tests)
+## 
+## Response: DGEN_SCORE
+##                                Effect      df   MSE        F   ges p.value
+## 1                         perspective  1, 518 36.26     1.04  .002    .309
+## 2                          temporalDO  1, 518 36.26     0.03 <.001    .862
+## 3              perspective:temporalDO  1, 518 36.26     0.00 <.001    .981
+## 4                                TIME  1, 518  3.11  8.39 **  .001    .004
+## 5                    perspective:TIME  1, 518  3.11     0.02 <.001    .890
+## 6                     temporalDO:TIME  1, 518  3.11   2.94 + <.001    .087
+## 7         perspective:temporalDO:TIME  1, 518  3.11   3.44 + <.001    .064
+## 8                              DOMAIN 2, 1036  2.12 7.85 ***  .001   <.001
+## 9                  perspective:DOMAIN 2, 1036  2.12     1.17 <.001    .312
+## 10                  temporalDO:DOMAIN 2, 1036  2.12  5.00 ** <.001    .007
+## 11      perspective:temporalDO:DOMAIN 2, 1036  2.12     0.39 <.001    .680
+## 12                        TIME:DOMAIN 2, 1036  1.52     0.77 <.001    .461
+## 13            perspective:TIME:DOMAIN 2, 1036  1.52     0.67 <.001    .513
+## 14             temporalDO:TIME:DOMAIN 2, 1036  1.52     0.44 <.001    .643
+## 15 perspective:temporalDO:TIME:DOMAIN 2, 1036  1.52   3.12 * <.001    .045
+## ---
+## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
+

Mauchly’s test of sphericity is not significant. Using uncorrected +values for interpretation and analysis.

+

Significant main effects and interactions:
+Effect df MSE F ges p 4 TIME 1, 518 3.11 8.39 ** .001 .004
+8 DOMAIN 2, 1036 2.13 7.85 *** .001 <.001
+10 temporalDO:DOMAIN 2, 1036 2.13 5.00 ** <.001 .007
+15 perspective:temporalDO:TIME:DOMAIN 2, 1036 1.52 3.12 * <.001 +.045

+
+
+

Mauchly and epsilon

+
anova_wide <- anova_data %>%
+  select(pID, perspective, temporalDO, all_of(within_vars)) %>%
+  filter(if_all(all_of(within_vars), ~ !is.na(.)))
+response_matrix <- as.matrix(anova_wide[, within_vars])
+rm_model <- lm(response_matrix ~ perspective * temporalDO, data = anova_wide)
+idata <- data.frame(
+  TIME   = factor(rep(c("past", "fut"), each = 3), levels = c("past", "fut")),
+  DOMAIN = factor(rep(c("pref", "pers", "val"), 2), levels = c("pref", "pers", "val"))
+)
+rm_anova <- car::Anova(rm_model, idata = idata, idesign = ~ TIME * DOMAIN, type = 3)
+rm_summary <- summary(rm_anova, multivariate = FALSE)
+
## Warning in summary.Anova.mlm(rm_anova, multivariate = FALSE): HF eps > 1
+## treated as 1
+
if (!is.null(rm_summary$sphericity.tests)) {
+  print(rm_summary$sphericity.tests)
+}
+
##                                    Test statistic p-value
+## DOMAIN                                    0.99806 0.60476
+## perspective:DOMAIN                        0.99806 0.60476
+## temporalDO:DOMAIN                         0.99806 0.60476
+## perspective:temporalDO:DOMAIN             0.99806 0.60476
+## TIME:DOMAIN                               0.99654 0.40772
+## perspective:TIME:DOMAIN                   0.99654 0.40772
+## temporalDO:TIME:DOMAIN                    0.99654 0.40772
+## perspective:temporalDO:TIME:DOMAIN        0.99654 0.40772
+
if (!is.null(rm_summary$epsilon)) {
+  print(rm_summary$epsilon)
+}
+
+
+

Post hoc comparisons

+
+

TIME (main effect)

+
emm_TIME <- emmeans(aov_afex, ~ TIME)
+print(pairs(emm_TIME, adjust = "none"))
+
##  contrast   estimate     SE  df t.ratio p.value
+##  past - fut    0.183 0.0632 518   2.897  0.0039
+## 
+## Results are averaged over the levels of: perspective, temporalDO, DOMAIN
+

Supports presence of EHI effect.

+
+
+

DOMAIN (main effect)

+
emm_DOMAIN <- emmeans(aov_afex, ~ DOMAIN)
+print(pairs(emm_DOMAIN, adjust = "tukey"))
+
##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers    0.129 0.0629 518   2.043  0.1032
+##  pref - val     0.253 0.0634 518   3.991  0.0002
+##  pers - val     0.124 0.0652 518   1.908  0.1375
+## 
+## Results are averaged over the levels of: perspective, temporalDO, TIME 
+## P value adjustment: tukey method for comparing a family of 3 estimates
+

Only preference to values contrast is significant.

+
+
+

temporalDO:DOMAIN

+
emmeans(aov_afex, pairwise ~ temporalDO | DOMAIN, adjust = "none")$contrasts
+
## DOMAIN = pref:
+##  contrast      estimate    SE  df t.ratio p.value
+##  future - past -0.00838 0.220 518  -0.038  0.9697
+## 
+## DOMAIN = pers:
+##  contrast      estimate    SE  df t.ratio p.value
+##  future - past -0.25252 0.230 518  -1.096  0.2735
+## 
+## DOMAIN = val:
+##  contrast      estimate    SE  df t.ratio p.value
+##  future - past  0.14817 0.233 518   0.636  0.5248
+## 
+## Results are averaged over the levels of: perspective, TIME
+
emmeans(aov_afex, pairwise ~ DOMAIN | temporalDO, adjust = "tukey")$contrasts
+
## temporalDO = future:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers  0.25065 0.0892 518   2.810  0.0142
+##  pref - val   0.17474 0.0898 518   1.945  0.1271
+##  pers - val  -0.07591 0.0924 518  -0.821  0.6900
+## 
+## temporalDO = past:
+##  contrast    estimate     SE  df t.ratio p.value
+##  pref - pers  0.00651 0.0888 518   0.073  0.9970
+##  pref - val   0.33129 0.0895 518   3.702  0.0007
+##  pers - val   0.32478 0.0921 518   3.527  0.0013
+## 
+## Results are averaged over the levels of: perspective, TIME 
+## P value adjustment: tukey method for comparing a family of 3 estimates
+

When grouped by domain, no contrasts are significant.

+

When grouped by temporalDO, some contrasts are significant:

+

Future-first temporal display order:
+contrast estimate SE df t.ratio p.value
+pref - pers 0.25065 0.0892 518 2.810 0.0142

+

Past-first temporal display order:
+contrast estimate SE df t.ratio p.value
+pref - val 0.33129 0.0895 518 3.702 0.0007
+pers - val 0.32478 0.0921 518 3.527 0.0013

+
+
+

perspective:temporalDO:TIME:DOMAIN

+
+

contrasts for TIME grouped by perspective, temporalDO, and +DOMAIN

+
emm_fourway <- emmeans(aov_afex, pairwise ~ TIME | perspective * temporalDO * DOMAIN, adjust = "tukey")
+print(emm_fourway$contrasts)
+
## perspective = other, temporalDO = future, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.0492 0.182 518   0.270  0.7873
+## 
+## perspective = self, temporalDO = future, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut  -0.0362 0.171 518  -0.212  0.8326
+## 
+## perspective = other, temporalDO = past, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.3237 0.171 518   1.897  0.0584
+## 
+## perspective = self, temporalDO = past, DOMAIN = pref:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.1870 0.181 518   1.031  0.3032
+## 
+## perspective = other, temporalDO = future, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.2049 0.179 518   1.142  0.2540
+## 
+## perspective = self, temporalDO = future, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.0145 0.169 518   0.086  0.9316
+## 
+## perspective = other, temporalDO = past, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.2878 0.168 518   1.712  0.0875
+## 
+## perspective = self, temporalDO = past, DOMAIN = pers:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.5285 0.179 518   2.957  0.0032
+## 
+## perspective = other, temporalDO = future, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.2951 0.188 518   1.568  0.1175
+## 
+## perspective = self, temporalDO = future, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut  -0.0797 0.177 518  -0.450  0.6526
+## 
+## perspective = other, temporalDO = past, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut  -0.1151 0.176 518  -0.653  0.5141
+## 
+## perspective = self, temporalDO = past, DOMAIN = val:
+##  contrast   estimate    SE  df t.ratio p.value
+##  past - fut   0.5366 0.187 518   2.863  0.0044
+

Significant contrasts:

+

contrast estimate SE df t.ratio p.value
+past - fut 0.5285 0.179 518 2.957 0.0032 (self-perspective, personality +domain, past-first temporal display order)
+past - fut 0.5366 0.187 518 2.863 0.0044 (self-perspective, values +domain, past-first temporal display order)

+

### contrasts for DOMAIN grouped by perspective, TIME, and +temporalDO

+
emm_fourway2 <- emmeans(aov_afex, pairwise ~ DOMAIN | perspective * TIME * temporalDO, adjust = "tukey")
+print(emm_fourway2$contrasts)
+
## perspective = other, TIME = past, temporalDO = future:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers   0.2459 0.168 518   1.468  0.3073
+##  pref - val    0.0164 0.177 518   0.093  0.9953
+##  pers - val   -0.2295 0.171 518  -1.343  0.3718
+## 
+## perspective = self, TIME = past, temporalDO = future:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers   0.1522 0.158 518   0.966  0.5986
+##  pref - val    0.2319 0.166 518   1.395  0.3445
+##  pers - val    0.0797 0.161 518   0.496  0.8732
+## 
+## perspective = other, TIME = fut, temporalDO = future:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers   0.4016 0.177 518   2.268  0.0613
+##  pref - val    0.2623 0.169 518   1.552  0.2678
+##  pers - val   -0.1393 0.175 518  -0.798  0.7046
+## 
+## perspective = self, TIME = fut, temporalDO = future:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers   0.2029 0.167 518   1.219  0.4427
+##  pref - val    0.1884 0.159 518   1.185  0.4625
+##  pers - val   -0.0145 0.164 518  -0.088  0.9957
+## 
+## perspective = other, TIME = past, temporalDO = past:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers   0.1367 0.157 518   0.871  0.6589
+##  pref - val    0.6259 0.166 518   3.778  0.0005
+##  pers - val    0.4892 0.160 518   3.056  0.0066
+## 
+## perspective = self, TIME = past, temporalDO = past:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers  -0.2764 0.167 518  -1.657  0.2230
+##  pref - val    0.0813 0.176 518   0.462  0.8892
+##  pers - val    0.3577 0.170 518   2.102  0.0903
+## 
+## perspective = other, TIME = fut, temporalDO = past:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers   0.1007 0.166 518   0.607  0.8163
+##  pref - val    0.1871 0.158 518   1.181  0.4650
+##  pers - val    0.0863 0.164 518   0.528  0.8579
+## 
+## perspective = self, TIME = fut, temporalDO = past:
+##  contrast    estimate    SE  df t.ratio p.value
+##  pref - pers   0.0650 0.176 518   0.369  0.9278
+##  pref - val    0.4309 0.168 518   2.559  0.0290
+##  pers - val    0.3659 0.174 518   2.103  0.0902
+## 
+## P value adjustment: tukey method for comparing a family of 3 estimates
+

Significant contrasts:

+

contrast estimate SE df t.ratio p.value
+pref - val 0.6259 0.166 518 3.778 0.0005 (other-perspective, +past-directed questions, past-first temporal display order)
+pers - val 0.4892 0.160 518 3.056 0.0066 (other-perspective, +past-directed questions, past-first temporal display order)
+pref - val 0.4309 0.168 518 2.559 0.0290 (self-perspective, +future-directed questions, past-first temporal display order)

+
+
+
+

Cohen’s d (significant contrasts only)

+
d_data <- anova_data %>%
+  mutate(
+    past_mean = (past_pref_DGEN + past_pers_DGEN + past_val_DGEN) / 3,
+    fut_mean  = (fut_pref_DGEN + fut_pers_DGEN + fut_val_DGEN) / 3,
+    pref_mean = (past_pref_DGEN + fut_pref_DGEN) / 2,
+    pers_mean = (past_pers_DGEN + fut_pers_DGEN) / 2,
+    val_mean  = (past_val_DGEN + fut_val_DGEN) / 2
+  )
+
+cohens_d_results <- tibble(
+  contrast = character(),
+  condition = character(),
+  d = double()
+)
+
+# TIME main: past vs fut
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "overall",
+    d = suppressMessages(effectsize::cohens_d(d_data$past_mean, d_data$fut_mean, paired = TRUE)$Cohens_d)
+  )
+
+# DOMAIN main: pref vs val
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - val)",
+    condition = "overall",
+    d = suppressMessages(effectsize::cohens_d(d_data$pref_mean, d_data$val_mean, paired = TRUE)$Cohens_d)
+  )
+
+# temporalDO:DOMAIN - future: pref vs pers
+d_fut <- d_data %>% filter(temporalDO == "future")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - pers)",
+    condition = "temporalDO = future",
+    d = suppressMessages(effectsize::cohens_d(d_fut$pref_mean, d_fut$pers_mean, paired = TRUE)$Cohens_d)
+  )
+
+# temporalDO:DOMAIN - past: pref vs val, pers vs val
+d_past <- d_data %>% filter(temporalDO == "past")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - val)",
+    condition = "temporalDO = past",
+    d = suppressMessages(effectsize::cohens_d(d_past$pref_mean, d_past$val_mean, paired = TRUE)$Cohens_d)
+  ) %>%
+  add_row(
+    contrast = "DOMAIN (pers - val)",
+    condition = "temporalDO = past",
+    d = suppressMessages(effectsize::cohens_d(d_past$pers_mean, d_past$val_mean, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way TIME: self, past temporalDO, pers
+d_self_past <- d_data %>% filter(perspective == "self", temporalDO == "past")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "self, past temporalDO, pers domain",
+    d = suppressMessages(effectsize::cohens_d(d_self_past$past_pers_DGEN, d_self_past$fut_pers_DGEN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way TIME: self, past temporalDO, val
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "TIME (past - fut)",
+    condition = "self, past temporalDO, val domain",
+    d = suppressMessages(effectsize::cohens_d(d_self_past$past_val_DGEN, d_self_past$fut_val_DGEN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way DOMAIN: other, past TIME, past temporalDO - pref vs val
+d_other_past_tpast <- d_data %>% filter(perspective == "other", temporalDO == "past")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - val)",
+    condition = "other, past TIME, past temporalDO",
+    d = suppressMessages(effectsize::cohens_d(d_other_past_tpast$past_pref_DGEN, d_other_past_tpast$past_val_DGEN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way DOMAIN: other, past TIME, past temporalDO - pers vs val
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pers - val)",
+    condition = "other, past TIME, past temporalDO",
+    d = suppressMessages(effectsize::cohens_d(d_other_past_tpast$past_pers_DGEN, d_other_past_tpast$past_val_DGEN, paired = TRUE)$Cohens_d)
+  )
+
+# 4-way DOMAIN: self, fut TIME, past temporalDO - pref vs val
+d_self_fut_tpast <- d_data %>% filter(perspective == "self", temporalDO == "past")
+cohens_d_results <- cohens_d_results %>%
+  add_row(
+    contrast = "DOMAIN (pref - val)",
+    condition = "self, fut TIME, past temporalDO",
+    d = suppressMessages(effectsize::cohens_d(d_self_fut_tpast$fut_pref_DGEN, d_self_fut_tpast$fut_val_DGEN, paired = TRUE)$Cohens_d)
+  )
+
+cohens_d_results %>%
+  mutate(d = round(d, 3)) %>%
+  print(n = Inf)
+
## # A tibble: 10 × 3
+##    contrast             condition                              d
+##    <chr>                <chr>                              <dbl>
+##  1 TIME (past - fut)    overall                            0.122
+##  2 DOMAIN (pref - val)  overall                            0.178
+##  3 DOMAIN (pref - pers) temporalDO = future                0.174
+##  4 DOMAIN (pref - val)  temporalDO = past                  0.242
+##  5 DOMAIN (pers - val)  temporalDO = past                  0.239
+##  6 TIME (past - fut)    self, past temporalDO, pers domain 0.26 
+##  7 TIME (past - fut)    self, past temporalDO, val domain  0.254
+##  8 DOMAIN (pref - val)  other, past TIME, past temporalDO  0.348
+##  9 DOMAIN (pers - val)  other, past TIME, past temporalDO  0.241
+## 10 DOMAIN (pref - val)  self, fut TIME, past temporalDO    0.257
+

Size d Interpretation
+Small 0.2 Weak effect
+Medium 0.5 Moderate effect
+Large 0.8 Strong effect

+
+
+ + + +
+
+ +
+ + + + + + + + + + + + + + + + + diff --git a/eohi3/knit/DA02_anova_DGEN.rmd b/eohi3/knit/DA02_anova_DGEN.rmd new file mode 100644 index 0000000..657e09c --- /dev/null +++ b/eohi3/knit/DA02_anova_DGEN.rmd @@ -0,0 +1,434 @@ +--- +title: "Mixed ANOVA - Domain General Vars" +author: "" +date: "`r Sys.Date()`" +output: + html_document: + toc: true + toc_float: true + code_folding: hide +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = TRUE) +``` + +# Setup + +```{r libraries} +library(tidyverse) +library(rstatix) +library(emmeans) +library(effectsize) +library(afex) +library(car) + +options(scipen = 999) +afex::set_sum_contrasts() +``` + +# Data + +```{r read-data} +# Data file is in parent of knit/ (eohi3/eohi3.csv) +df <- read.csv( + "/home/ladmin/Documents/DND/EOHI/eohi3/eohi3.csv", + stringsAsFactors = FALSE, + check.names = FALSE, + na.strings = "NA" +) + +between_vars <- c("perspective", "temporalDO") +within_vars <- c( + "past_pref_DGEN", "past_pers_DGEN", "past_val_DGEN", + "fut_pref_DGEN", "fut_pers_DGEN", "fut_val_DGEN" +) + +missing_vars <- setdiff(c(between_vars, within_vars, "pID"), names(df)) +if (length(missing_vars) > 0) { + stop(paste("Missing required variables:", paste(missing_vars, collapse = ", "))) +} + +anova_data <- df %>% + select(pID, all_of(between_vars), all_of(within_vars)) %>% + filter( + !is.na(perspective), perspective != "", + !is.na(temporalDO), temporalDO != "" + ) +``` + +# Long format + +```{r long-format} +long_data <- anova_data %>% + pivot_longer( + cols = all_of(within_vars), + names_to = "variable", + values_to = "DGEN_SCORE" + ) %>% + mutate( + time = case_when( + grepl("^past_", variable) ~ "past", + grepl("^fut_", variable) ~ "fut", + TRUE ~ NA_character_ + ), + domain = case_when( + grepl("_pref_DGEN$", variable) ~ "pref", + grepl("_pers_DGEN$", variable) ~ "pers", + grepl("_val_DGEN$", variable) ~ "val", + TRUE ~ NA_character_ + ) + ) %>% + mutate( + TIME = factor(time, levels = c("past", "fut")), + DOMAIN = factor(domain, levels = c("pref", "pers", "val")), + perspective = factor(perspective), + temporalDO = factor(temporalDO), + pID = factor(pID) + ) %>% + select(pID, perspective, temporalDO, TIME, DOMAIN, DGEN_SCORE) %>% + filter(!is.na(DGEN_SCORE)) +``` + +# Descriptive statistics + +```{r desc-stats} +desc_stats <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + mean = round(mean(DGEN_SCORE), 5), + variance = round(var(DGEN_SCORE), 5), + sd = round(sd(DGEN_SCORE), 5), + median = round(median(DGEN_SCORE), 5), + q1 = round(quantile(DGEN_SCORE, 0.25), 5), + q3 = round(quantile(DGEN_SCORE, 0.75), 5), + min = round(min(DGEN_SCORE), 5), + max = round(max(DGEN_SCORE), 5), + .groups = "drop" + ) + +# Show all rows and columns (no truncation) +options(tibble.width = Inf) +print(desc_stats, n = Inf) +``` + +Interpretation: +1. Mean and median values are similar w/ slightly more variation than in the domain specific anova. +2. Highest to lowest group n size ratio is 1.14 (139/122). Acceptable ratio for ANOVA (under 1.5). +3. Highest to lowest overall group variance ratio is 1.40 (9.32/6.65). Acceptable ratio for ANOVA (under 4). + For the sake of consistency w/ the other EHI studies, I also calculated Hartley's F-max ratio. + The conservative F-max critical value is 1.60 (same as DS anova since number of groups and n sizes doesn't change), which is still higher than the highest observed F-max ratio of 1.28. + +# Assumption checks + +## Missing values + +```{r missing} +missing_summary <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n_total = n(), + n_missing = sum(is.na(DGEN_SCORE)), + pct_missing = round(100 * n_missing / n_total, 2), + .groups = "drop" + ) + +print(missing_summary, n = Inf) +``` + +No missing values. As expected. + +## Outliers + +```{r outliers} +outlier_summary <- long_data %>% + group_by(perspective, temporalDO, TIME, DOMAIN) %>% + summarise( + n = n(), + n_outliers = sum(abs(scale(DGEN_SCORE)) > 3), + .groups = "drop" + ) + +print(outlier_summary, n = Inf) +``` + +No outliers present. Good. Same as DS anova. + +## Homogeneity of variance + +```{r homogeneity} +homogeneity_between <- long_data %>% + group_by(TIME, DOMAIN) %>% + rstatix::levene_test(DGEN_SCORE ~ perspective * temporalDO) + +print(homogeneity_between, n = Inf) +``` + +Levene's test is signiicant for 1 cell only: past-val. However, variance ratios and F-max tests show that any heteroscedasticity is mild. + +## Normality (within-subjects residuals) + +```{r normality} +resid_within <- long_data %>% + group_by(pID) %>% + mutate(person_mean = mean(DGEN_SCORE, na.rm = TRUE)) %>% + ungroup() %>% + mutate(resid = DGEN_SCORE - person_mean) %>% + pull(resid) +resid_within <- resid_within[!is.na(resid_within)] + +n_resid <- length(resid_within) +if (n_resid < 3L) { + message("Too few within-subjects residuals (n < 3); skipping Shapiro-Wilk.") +} else { + resid_for_shapiro <- if (n_resid > 5000L) { + set.seed(1L) + sample(resid_within, 5000L) + } else { + resid_within + } + print(shapiro.test(resid_for_shapiro)) +} +``` + +### Q-Q plot + +```{r qqplot, fig.height = 4} +qqnorm(resid_within) +qqline(resid_within) +``` + +Shapiro-Wilk test is significant but is sensitive to large sample size. +QQ plot shows that strict centre residuals are normally distributed, but there is some deviation from normality. +ANOVA is robust to violations of normality w/ large sample size. + +Overall, ANOVA can proceed. + +# Mixed ANOVA + +```{r anova} +aov_afex <- aov_ez( + id = "pID", + dv = "DGEN_SCORE", + data = long_data, + between = c("perspective", "temporalDO"), + within = c("TIME", "DOMAIN"), + type = 3, + anova_table = list(correction = "none") +) + +print(aov_afex) +``` + +Mauchly's test of sphericity is not significant. Using uncorrected values for interpretation and analysis. + + +Significant main effects and interactions: + Effect df MSE F ges p +4 TIME 1, 518 3.11 8.39 ** .001 .004 +8 DOMAIN 2, 1036 2.13 7.85 *** .001 <.001 +10 temporalDO:DOMAIN 2, 1036 2.13 5.00 ** <.001 .007 +15 perspective:temporalDO:TIME:DOMAIN 2, 1036 1.52 3.12 * <.001 .045 + + +# Mauchly and epsilon + +```{r mauchly} +anova_wide <- anova_data %>% + select(pID, perspective, temporalDO, all_of(within_vars)) %>% + filter(if_all(all_of(within_vars), ~ !is.na(.))) +response_matrix <- as.matrix(anova_wide[, within_vars]) +rm_model <- lm(response_matrix ~ perspective * temporalDO, data = anova_wide) +idata <- data.frame( + TIME = factor(rep(c("past", "fut"), each = 3), levels = c("past", "fut")), + DOMAIN = factor(rep(c("pref", "pers", "val"), 2), levels = c("pref", "pers", "val")) +) +rm_anova <- car::Anova(rm_model, idata = idata, idesign = ~ TIME * DOMAIN, type = 3) +rm_summary <- summary(rm_anova, multivariate = FALSE) +if (!is.null(rm_summary$sphericity.tests)) { + print(rm_summary$sphericity.tests) +} +if (!is.null(rm_summary$epsilon)) { + print(rm_summary$epsilon) +} +``` + +# Post hoc comparisons + +## TIME (main effect) + +```{r posthoc-TIME} +emm_TIME <- emmeans(aov_afex, ~ TIME) +print(pairs(emm_TIME, adjust = "none")) +``` + +Supports presence of EHI effect. + +## DOMAIN (main effect) + +```{r posthoc-domain} +emm_DOMAIN <- emmeans(aov_afex, ~ DOMAIN) +print(pairs(emm_DOMAIN, adjust = "tukey")) +``` + +Only preference to values contrast is significant. + +## temporalDO:DOMAIN + +```{r posthoc-temporaldo-domain} +emmeans(aov_afex, pairwise ~ temporalDO | DOMAIN, adjust = "none")$contrasts +emmeans(aov_afex, pairwise ~ DOMAIN | temporalDO, adjust = "tukey")$contrasts +``` + +When grouped by domain, no contrasts are significant. + + +When grouped by temporalDO, some contrasts are significant: + +Future-first temporal display order: +contrast estimate SE df t.ratio p.value +pref - pers 0.25065 0.0892 518 2.810 0.0142 + + +Past-first temporal display order: +contrast estimate SE df t.ratio p.value +pref - val 0.33129 0.0895 518 3.702 0.0007 +pers - val 0.32478 0.0921 518 3.527 0.0013 + +## perspective:temporalDO:TIME:DOMAIN + +### contrasts for TIME grouped by perspective, temporalDO, and DOMAIN +```{r posthoc-fourway} +emm_fourway <- emmeans(aov_afex, pairwise ~ TIME | perspective * temporalDO * DOMAIN, adjust = "tukey") +print(emm_fourway$contrasts) +``` + +Significant contrasts: + +contrast estimate SE df t.ratio p.value + past - fut 0.5285 0.179 518 2.957 0.0032 (self-perspective, personality domain, past-first temporal display order) + past - fut 0.5366 0.187 518 2.863 0.0044 (self-perspective, values domain, past-first temporal display order) + + ### contrasts for DOMAIN grouped by perspective, TIME, and temporalDO +```{r posthoc-fourway2} +emm_fourway2 <- emmeans(aov_afex, pairwise ~ DOMAIN | perspective * TIME * temporalDO, adjust = "tukey") +print(emm_fourway2$contrasts) +``` + +Significant contrasts: + +contrast estimate SE df t.ratio p.value + pref - val 0.6259 0.166 518 3.778 0.0005 (other-perspective, past-directed questions, past-first temporal display order) + pers - val 0.4892 0.160 518 3.056 0.0066 (other-perspective, past-directed questions, past-first temporal display order) + pref - val 0.4309 0.168 518 2.559 0.0290 (self-perspective, future-directed questions, past-first temporal display order) + +## Cohen's d (significant contrasts only) + +```{r cohens-d-significant} +d_data <- anova_data %>% + mutate( + past_mean = (past_pref_DGEN + past_pers_DGEN + past_val_DGEN) / 3, + fut_mean = (fut_pref_DGEN + fut_pers_DGEN + fut_val_DGEN) / 3, + pref_mean = (past_pref_DGEN + fut_pref_DGEN) / 2, + pers_mean = (past_pers_DGEN + fut_pers_DGEN) / 2, + val_mean = (past_val_DGEN + fut_val_DGEN) / 2 + ) + +cohens_d_results <- tibble( + contrast = character(), + condition = character(), + d = double() +) + +# TIME main: past vs fut +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "overall", + d = suppressMessages(effectsize::cohens_d(d_data$past_mean, d_data$fut_mean, paired = TRUE)$Cohens_d) + ) + +# DOMAIN main: pref vs val +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - val)", + condition = "overall", + d = suppressMessages(effectsize::cohens_d(d_data$pref_mean, d_data$val_mean, paired = TRUE)$Cohens_d) + ) + +# temporalDO:DOMAIN - future: pref vs pers +d_fut <- d_data %>% filter(temporalDO == "future") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - pers)", + condition = "temporalDO = future", + d = suppressMessages(effectsize::cohens_d(d_fut$pref_mean, d_fut$pers_mean, paired = TRUE)$Cohens_d) + ) + +# temporalDO:DOMAIN - past: pref vs val, pers vs val +d_past <- d_data %>% filter(temporalDO == "past") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - val)", + condition = "temporalDO = past", + d = suppressMessages(effectsize::cohens_d(d_past$pref_mean, d_past$val_mean, paired = TRUE)$Cohens_d) + ) %>% + add_row( + contrast = "DOMAIN (pers - val)", + condition = "temporalDO = past", + d = suppressMessages(effectsize::cohens_d(d_past$pers_mean, d_past$val_mean, paired = TRUE)$Cohens_d) + ) + +# 4-way TIME: self, past temporalDO, pers +d_self_past <- d_data %>% filter(perspective == "self", temporalDO == "past") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "self, past temporalDO, pers domain", + d = suppressMessages(effectsize::cohens_d(d_self_past$past_pers_DGEN, d_self_past$fut_pers_DGEN, paired = TRUE)$Cohens_d) + ) + +# 4-way TIME: self, past temporalDO, val +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "TIME (past - fut)", + condition = "self, past temporalDO, val domain", + d = suppressMessages(effectsize::cohens_d(d_self_past$past_val_DGEN, d_self_past$fut_val_DGEN, paired = TRUE)$Cohens_d) + ) + +# 4-way DOMAIN: other, past TIME, past temporalDO - pref vs val +d_other_past_tpast <- d_data %>% filter(perspective == "other", temporalDO == "past") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - val)", + condition = "other, past TIME, past temporalDO", + d = suppressMessages(effectsize::cohens_d(d_other_past_tpast$past_pref_DGEN, d_other_past_tpast$past_val_DGEN, paired = TRUE)$Cohens_d) + ) + +# 4-way DOMAIN: other, past TIME, past temporalDO - pers vs val +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pers - val)", + condition = "other, past TIME, past temporalDO", + d = suppressMessages(effectsize::cohens_d(d_other_past_tpast$past_pers_DGEN, d_other_past_tpast$past_val_DGEN, paired = TRUE)$Cohens_d) + ) + +# 4-way DOMAIN: self, fut TIME, past temporalDO - pref vs val +d_self_fut_tpast <- d_data %>% filter(perspective == "self", temporalDO == "past") +cohens_d_results <- cohens_d_results %>% + add_row( + contrast = "DOMAIN (pref - val)", + condition = "self, fut TIME, past temporalDO", + d = suppressMessages(effectsize::cohens_d(d_self_fut_tpast$fut_pref_DGEN, d_self_fut_tpast$fut_val_DGEN, paired = TRUE)$Cohens_d) + ) + +cohens_d_results %>% + mutate(d = round(d, 3)) %>% + print(n = Inf) +``` + +Size d Interpretation +Small 0.2 Weak effect +Medium 0.5 Moderate effect +Large 0.8 Strong effect \ No newline at end of file