Overview

This document provides a comprehensive analysis of emotional responses in two instructional conditions:

PF (Productive Failure): Students work on problems before receiving instruction
PS (Productive Success): Students receive instruction before working on problems

Emotions are measured using the VAD (Valence-Arousal-Dominance) framework:

Valence: Negative (1) to Positive (9)
Arousal: Calm (1) to Excited (9)
Dominance: Controlled (1) to In Control (9)

1. Setup & Data Processing

1.1 Load Libraries

# Data manipulation
library(tidyverse)
library(readxl)
library(writexl)
library(stringr)
library(purrr)
library(tools)

# Statistical analysis
library(lme4)        # Mixed models
library(lmerTest)    # P-values for mixed models
library(emmeans)     # Post-hoc comparisons
library(effectsize)  # Effect sizes

# Visualization
library(ggpubr)      # Publication-ready plots
library(patchwork)   # Combining plots
library(scales)      # Formatting

1.2 Read Raw Data

# File list for emotion survey data
file_list <- c(
  "surveys/ES-Baseline-Open.xlsx",
  "surveys/ES-PF1-Open.xlsx",
  "surveys/ES-PF2-Open.xlsx",
  "surveys/ES-PF3-Open.xlsx",
  "surveys/ES-PS1-Open.xlsx",
  "surveys/ES-PS2-Open.xlsx",
  "surveys/ES-PS3-Open.xlsx"
)

# Function to process each file
process_file <- function(file_path) {
  file_name <- file_path_sans_ext(basename(file_path))
  
  # Extract condition and timepoint
  if (grepl("Baseline", file_name, ignore.case = TRUE)) {
    condition <- NA
    timepoint <- "Baseline"
  } else {
    condition <- str_extract(file_name, "PF|PS")
    timepoint <- str_extract(file_name, "\\d+") %>% paste0("T", .)
  }
  
  # Read the file
  df <- read_excel(file_path)
  names(df)[1] <- "Email_Address"
  
  # Pivot emotion words
  word_df <- df %>%
    pivot_longer(
      cols = ends_with(".word"),
      names_to = "Context",
      values_to = "Emotion_Word"
    ) %>%
    mutate(Context = str_remove(Context, "\\.word$"))
  
  # Pivot explanations
  why_df <- df %>%
    pivot_longer(
      cols = ends_with(".why"),
      names_to = "Context",
      values_to = "Explanation"
    ) %>%
    mutate(Context = str_remove(Context, "\\.why$"))
  
  # Join and add metadata
  df_long <- left_join(word_df, why_df, by = c("Email_Address", "Context")) %>%
    mutate(
      Condition = condition,
      Timepoint = timepoint
    )
  
  return(df_long)
}

# Combine all files
all_data <- map_dfr(file_list, process_file)

# Factor the timepoint
all_data <- all_data %>%
  mutate(Timepoint = factor(Timepoint, levels = c("Baseline", "T1", "T2", "T3")))

# Create participant ID lookup
id_lookup <- all_data %>%
  distinct(Email_Address) %>%
  arrange(Email_Address) %>%
  mutate(ID = paste0("P", row_number()))

# Add ID and remove Email_Address
all_data <- all_data %>%
  left_join(id_lookup, by = "Email_Address") %>%
  select(ID, Context, Emotion_Word, Explanation, Condition, Timepoint)

cat(sprintf("Data loaded: %d observations from %d participants\n", 
            nrow(all_data), n_distinct(all_data$ID)))

## Data loaded: 483 observations from 45 participants

2. Emotion Word Consolidation

2.1 Define Emotion Lexicon

# Comprehensive emotion word dictionary
emotion_lexicon <- list(
  # Positive emotions
  positive = c("happy", "excited", "pleased", "satisfied", "joyful", "glad",
               "delighted", "grateful", "gratified", "proud", "confident",
               "comfortable", "content", "relieved", "encouraged", "hopeful",
               "energized", "enlightened", "good", "better", "fine", "okay",
               "nice", "reassured", "alleviated", "comforted"),
  
  # Negative emotions
  negative = c("sad", "unhappy", "disappointed", "discouraged", "frustrated",
               "angry", "annoyed", "upset", "stressed", "anxious", "worried",
               "scared", "nervous", "overwhelmed", "confused", "lost",
               "embarrassed", "ashamed", "uncomfortable", "tense", "silly"),
  
  # Neutral/Mixed
  neutral = c("challenged", "engaged", "focused", "busy", "indifferent",
              "uncertain", "ambivalent", "okay", "fine", "meh"),
  
  # Cognitive states
  cognitive = c("confused", "clarity", "understanding", "enlightened",
                "aware", "lost", "uncertain", "struggling")
)

# Flatten to single vector for matching
all_emotion_words <- unlist(emotion_lexicon)

2.2 Consolidation Function

# Function to consolidate multi-word entries to single words
consolidate_to_single_word <- function(text) {
  if (is.na(text) || text == "") return(NA)
  
  # Rule 0: If already single word, return it
  if (!str_detect(text, " |,|;|/|\\(|\\)")) {
    return(str_to_title(str_trim(text)))
  }
  
  # Convert to lowercase for processing
  text_lower <- tolower(str_trim(text))
  
  # Rule 1: Handle "and"/"but" patterns
  if (str_detect(text_lower, "and|but")) {
    parts <- str_split(text_lower, " and | but ")[[1]]
    words <- str_extract_all(parts[1], "\\w+")[[1]]
    for (word in words) {
      if (word %in% all_emotion_words) return(str_to_title(word))
    }
  }
  
  # Rule 2: Handle comma-separated lists
  if (str_detect(text_lower, ",")) {
    parts <- str_split(text_lower, ",")[[1]]
    first_part <- str_trim(parts[1])
    words <- str_extract_all(first_part, "\\w+")[[1]]
    for (word in words) {
      if (word %in% all_emotion_words) return(str_to_title(word))
    }
  }
  
  # Rule 3: Extract all words and find emotion words
  all_words <- str_extract_all(text_lower, "\\w+")[[1]]
  
  # Remove stop words
  stop_words <- c("i", "me", "my", "felt", "feel", "feeling", "was", "were",
                  "am", "is", "are", "the", "a", "an", "and", "or", "but",
                  "very", "really", "quite", "somewhat", "kind", "sort", "of",
                  "little", "bit", "more", "less", "it", "that", "this", "to",
                  "would", "could", "should", "like", "as", "when", "about",
                  "with", "working", "worked", "work", "classes", "class",
                  "classmates", "classmate", "teacher", "teachers", "help",
                  "helped", "helping", "problems", "problem", "difficult")
  
  meaningful_words <- all_words[!all_words %in% stop_words]
  
  # Rule 4: Find matching emotion words
  for (word in meaningful_words) {
    if (word %in% all_emotion_words) {
      return(str_to_title(word))
    }
  }
  
  # Rule 5: Look for emotion-indicating suffixes
  for (word in meaningful_words) {
    if (str_detect(word, "ed$|ing$|ful$")) {
      return(str_to_title(word))
    }
  }
  
  # Rule 6: Special case handling
  special_cases <- list(
    "a-ha moment" = "Enlightened",
    "ah-ha" = "Enlightened",
    "struggle bus" = "Struggling",
    "meh" = "Indifferent",
    "not applicable" = NA,
    "pretty good" = "Good",
    "pretty easy" = "Comfortable",
    "less confused" = "Clarity",
    "still confused" = "Confused",
    "more reassured" = "Reassured"
  )
  
  for (phrase in names(special_cases)) {
    if (str_detect(text_lower, fixed(phrase))) {
      return(special_cases[[phrase]])
    }
  }
  
  # Rule 7: Sentiment-based fallback
  positive_indicators <- c("good", "better", "best", "easy", "easier", "helpful",
                          "happy", "glad", "nice", "great", "relief", "relieved")
  negative_indicators <- c("hard", "harder", "difficult", "struggle", "struggling",
                          "bad", "worse", "worst", "scary", "scared")
  
  has_positive <- any(meaningful_words %in% positive_indicators)
  has_negative <- any(meaningful_words %in% negative_indicators)
  
  if (has_positive && !has_negative) return("Positive")
  if (has_negative && !has_positive) return("Challenged")
  
  # Rule 8: Use first meaningful word
  if (length(meaningful_words) > 0) {
    return(str_to_title(meaningful_words[1]))
  }
  
  # Ultimate fallback
  return(str_to_title(str_trim(text)))
}

2.3 Apply Consolidation

# Apply consolidation to the data
all_data <- all_data %>%
  mutate(
    Emotion_Word_Original = Emotion_Word,
    Is_Multi_Word = str_detect(Emotion_Word, " |,|;|/"),
    Emotion_Word_Single = sapply(Emotion_Word, consolidate_to_single_word),
    Was_Consolidated = Emotion_Word_Original != Emotion_Word_Single
  )

# Summary
cat("\n", rep("=", 60), "\n", sep="")

## 
## ============================================================

cat("EMOTION WORD CONSOLIDATION SUMMARY\n")

## EMOTION WORD CONSOLIDATION SUMMARY

cat(rep("=", 60), "\n", sep="")

## ============================================================

cat(sprintf("Total rows: %d\n", nrow(all_data)))

## Total rows: 483

cat(sprintf("Rows with emotion words: %d\n", sum(!is.na(all_data$Emotion_Word_Original))))

## Rows with emotion words: 483

cat(sprintf("Rows with multiple words: %d\n", sum(all_data$Is_Multi_Word, na.rm = TRUE)))

## Rows with multiple words: 118

cat(sprintf("Rows consolidated: %d\n", sum(all_data$Was_Consolidated, na.rm = TRUE)))

## Rows consolidated: 260

# Show top emotions
cat("\n", rep("=", 60), "\n", sep="")

## 
## ============================================================

cat("TOP 15 MOST FREQUENT CONSOLIDATED EMOTIONS\n")

## TOP 15 MOST FREQUENT CONSOLIDATED EMOTIONS

cat(rep("=", 60), "\n\n", sep="")

## ============================================================

top_emotions <- all_data %>%
  filter(!is.na(Emotion_Word_Single)) %>%
  count(Emotion_Word_Single, sort = TRUE) %>%
  head(15)
print(as.data.frame(top_emotions), row.names = FALSE)

##  Emotion_Word_Single  n
##                 Good 49
##             Confused 36
##             Relieved 24
##             Stressed 18
##              Helpful 17
##           Frustrated 15
##                 Fine 13
##                Great 13
##            Confident 12
##               Better 11
##                Happy 11
##                   Na 11
##                 Okay 11
##                 Calm  8
##          Indifferent  7

3. VAD Score Assignment

3.1 Load VAD Lexicon

# Read the VAD lexicon
vad_lexicon <- read_excel("surveys/VAD_Cleaned.xlsx") %>%
  select(
    Word,
    Valence = V.Mean.Sum,
    Arousal = A.Mean.Sum,
    Dominance = D.Mean.Sum
  ) %>%
  mutate(Word_Lower = tolower(Word))

cat(sprintf("VAD lexicon loaded: %d words\n", nrow(vad_lexicon)))

## VAD lexicon loaded: 13915 words

3.2 Manual VAD Mappings

# Manual mappings for emotions not in the VAD lexicon
manual_vad_mappings <- tibble(
  Word_Lower = c(
    "relieved", "stressed", "better", "okay", "challenged",
    "frstrated", "frusterated", "frustrated.", "releived", "relived",
    "encouraged", "energized", "enlightened", "gratified", "helped",
    "reassured", "comforted", "alleviated", "struggling", "struggled",
    "supported", "accomplished", "intellectually-stimulated", "intimidated",
    "panicked", "trepidatious", "cathartic", "collaborative"
  ),
  Valence = c(
    7.24, 2.25, 6.50, 5.50, 5.10,
    3.00, 3.00, 3.00, 7.24, 7.24,
    7.00, 7.20, 7.26, 7.00, 6.80, 7.00, 7.24, 7.24, 3.50, 3.50,
    6.80, 7.20, 7.00, 3.00, 2.00, 3.50, 6.50, 6.50
  ),
  Arousal = c(
    3.76, 6.68, 3.50, 3.00, 5.10,
    5.50, 5.50, 5.50, 3.76, 3.76,
    4.50, 6.00, 3.64, 4.00, 3.50, 3.00, 2.76, 3.76, 5.50, 5.50,
    3.50, 5.00, 6.00, 5.00, 7.00, 6.00, 4.00, 4.00
  ),
  Dominance = c(
    6.41, 3.19, 6.00, 5.50, 5.60,
    3.78, 3.78, 3.78, 6.41, 6.41,
    6.50, 6.50, 6.59, 6.50, 6.00, 6.50, 7.18, 6.41, 3.50, 3.50,
    6.50, 6.80, 6.50, 3.50, 2.50, 3.50, 6.00, 6.00
  ),
  Source = "Manual"
)

# Combine VAD lexicon with manual mappings
vad_complete <- vad_lexicon %>%
  select(Word_Lower, Valence, Arousal, Dominance) %>%
  mutate(Source = "Lexicon") %>%
  bind_rows(manual_vad_mappings %>% select(Word_Lower, Valence, Arousal, Dominance, Source))

cat(sprintf("Total VAD mappings available: %d\n", nrow(vad_complete)))

## Total VAD mappings available: 13943

3.3 Map VAD to Emotions

# Map VAD scores to emotion words
all_data <- all_data %>%
  mutate(Emotion_Word_Single_Lower = tolower(Emotion_Word_Single)) %>%
  left_join(
    vad_complete,
    by = c("Emotion_Word_Single_Lower" = "Word_Lower")
  ) %>%
  select(-Emotion_Word_Single_Lower)

# Check mapping success
mapping_summary <- all_data %>%
  filter(!is.na(Emotion_Word_Single)) %>%
  summarise(
    Total_Emotions = n(),
    With_VAD = sum(!is.na(Valence)),
    Without_VAD = sum(is.na(Valence)),
    Percent_Mapped = round(100 * sum(!is.na(Valence)) / n(), 1)
  )

cat("\n", rep("=", 60), "\n", sep="")

## 
## ============================================================

cat("VAD MAPPING SUMMARY\n")

## VAD MAPPING SUMMARY

cat(rep("=", 60), "\n", sep="")

## ============================================================

cat(sprintf("Total emotion entries: %d\n", mapping_summary$Total_Emotions))

## Total emotion entries: 482

cat(sprintf("Successfully mapped to VAD: %d (%.1f%%)\n",
            mapping_summary$With_VAD, mapping_summary$Percent_Mapped))

## Successfully mapped to VAD: 422 (87.6%)

cat(sprintf("Missing VAD scores: %d (%.1f%%)\n",
            mapping_summary$Without_VAD, 100 - mapping_summary$Percent_Mapped))

## Missing VAD scores: 60 (12.4%)

# Save processed data
write_xlsx(all_data, "tables/all_data_with_VAD.xlsx")

4. Descriptive Statistics

4.1 Prepare Analysis Data

# Filter and prepare data for analysis
analysis_data <- all_data %>%
  filter(!is.na(Valence) & !is.na(Condition)) %>%
  mutate(
    Timepoint = factor(Timepoint, levels = c("Baseline", "T1", "T2", "T3")),
    Condition = factor(Condition, levels = c("PF", "PS")),
    Context = factor(Context, levels = c("activity", "group", "social")),
    ID = factor(ID)
  )

cat("\n", rep("=", 80), "\n", sep="")

## 
## ================================================================================

cat("ANALYSIS SAMPLE\n")

## ANALYSIS SAMPLE

cat(rep("=", 80), "\n\n", sep="")

## ================================================================================

cat(sprintf("Total observations: %d\n", nrow(analysis_data)))

## Total observations: 328

cat(sprintf("Unique participants: %d\n", n_distinct(analysis_data$ID)))

## Unique participants: 41

cat(sprintf("PF observations: %d\n", sum(analysis_data$Condition == "PF")))

## PF observations: 172

cat(sprintf("PS observations: %d\n", sum(analysis_data$Condition == "PS")))

## PS observations: 156

4.2 Overall Descriptives

# Overall by Condition
overall_stats <- analysis_data %>%
  group_by(Condition) %>%
  summarise(
    n = n(),
    Valence_M = mean(Valence, na.rm = TRUE),
    Valence_SD = sd(Valence, na.rm = TRUE),
    Arousal_M = mean(Arousal, na.rm = TRUE),
    Arousal_SD = sd(Arousal, na.rm = TRUE),
    Dominance_M = mean(Dominance, na.rm = TRUE),
    Dominance_SD = sd(Dominance, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(overall_stats, digits = 2, caption = "Overall VAD by Condition")

Overall VAD by Condition
Condition	n	Valence_M	Valence_SD	Arousal_M	Arousal_SD	Dominance_M	Dominance_SD
PF	172	5.63	2.15	4.33	1.08	5.52	1.34
PS	156	6.43	1.79	4.06	1.04	6.04	1.11

4.3 By Condition × Timepoint

condition_time_stats <- analysis_data %>%
  group_by(Condition, Timepoint) %>%
  summarise(
    n = n(),
    Valence_M = mean(Valence, na.rm = TRUE),
    Valence_SD = sd(Valence, na.rm = TRUE),
    Arousal_M = mean(Arousal, na.rm = TRUE),
    Arousal_SD = sd(Arousal, na.rm = TRUE),
    Dominance_M = mean(Dominance, na.rm = TRUE),
    Dominance_SD = sd(Dominance, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(condition_time_stats, digits = 2, caption = "VAD by Condition × Timepoint")

VAD by Condition × Timepoint
Condition	Timepoint	n	Valence_M	Valence_SD	Arousal_M	Arousal_SD	Dominance_M	Dominance_SD
PF	T1	68	6.00	2.10	4.52	1.12	5.76	1.30
PF	T2	51	5.30	2.13	4.29	0.97	5.32	1.33
PF	T3	53	5.47	2.22	4.11	1.11	5.40	1.37
PS	T1	66	6.14	1.97	4.19	1.02	5.91	1.28
PS	T2	55	6.85	1.50	3.99	1.05	6.23	0.88
PS	T3	35	6.33	1.76	3.94	1.06	5.99	1.11

4.4 By Condition × Context

condition_context_stats <- analysis_data %>%
  group_by(Condition, Context) %>%
  summarise(
    n = n(),
    Valence_M = mean(Valence, na.rm = TRUE),
    Valence_SD = sd(Valence, na.rm = TRUE),
    Arousal_M = mean(Arousal, na.rm = TRUE),
    Arousal_SD = sd(Arousal, na.rm = TRUE),
    Dominance_M = mean(Dominance, na.rm = TRUE),
    Dominance_SD = sd(Dominance, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(condition_context_stats, digits = 2, caption = "VAD by Condition × Context")

VAD by Condition × Context
Condition	Context	n	Valence_M	Valence_SD	Arousal_M	Arousal_SD	Dominance_M	Dominance_SD
PF	activity	62	4.32	2.10	4.93	1.22	4.68	1.34
PF	group	55	6.57	1.70	4.04	0.86	6.07	1.05
PF	social	55	6.17	1.91	3.93	0.80	5.91	1.13
PS	activity	54	5.69	2.06	4.09	1.07	5.56	1.33
PS	group	56	6.89	1.49	4.12	1.20	6.38	0.87
PS	social	46	6.74	1.51	3.96	0.75	6.20	0.90

4.5 Three-Way Interaction

three_way_stats <- analysis_data %>%
  group_by(Condition, Timepoint, Context) %>%
  summarise(
    n = n(),
    Valence_M = mean(Valence, na.rm = TRUE),
    Valence_SD = sd(Valence, na.rm = TRUE),
    Arousal_M = mean(Arousal, na.rm = TRUE),
    Arousal_SD = sd(Arousal, na.rm = TRUE),
    Dominance_M = mean(Dominance, na.rm = TRUE),
    Dominance_SD = sd(Dominance, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(head(three_way_stats, 20), digits = 2, 
             caption = "VAD by Condition × Timepoint × Context (first 20 rows)")

VAD by Condition × Timepoint × Context (first 20 rows)
Condition	Timepoint	Context	n	Valence_M	Valence_SD	Arousal_M	Arousal_SD	Dominance_M	Dominance_SD
PF	T1	activity	23	5.03	2.41	5.22	1.10	5.09	1.52
PF	T1	group	23	6.54	1.88	4.37	0.97	6.12	1.01
PF	T1	social	22	6.45	1.65	3.95	0.92	6.09	1.07
PF	T2	activity	19	4.10	1.91	4.84	1.18	4.58	1.34
PF	T2	group	16	6.18	1.78	3.87	0.70	5.82	1.11
PF	T2	social	16	5.85	2.16	4.05	0.57	5.71	1.18
PF	T3	activity	20	3.70	1.73	4.69	1.37	4.32	1.03
PF	T3	group	16	7.00	1.33	3.72	0.68	6.26	1.04
PF	T3	social	17	6.10	2.03	3.80	0.84	5.86	1.19
PS	T1	activity	23	5.37	2.10	4.22	1.05	5.32	1.49
PS	T1	group	26	6.61	1.84	4.30	1.14	6.32	1.03
PS	T1	social	17	6.47	1.75	3.96	0.78	6.09	1.06
PS	T2	activity	19	6.21	1.94	3.90	1.11	5.92	1.15
PS	T2	group	18	7.26	1.08	4.12	1.23	6.42	0.66
PS	T2	social	18	7.11	1.13	3.97	0.81	6.36	0.66
PS	T3	activity	12	5.49	2.17	4.15	1.11	5.43	1.26
PS	T3	group	12	6.94	1.08	3.72	1.31	6.44	0.82
PS	T3	social	11	6.58	1.68	3.95	0.68	6.12	1.02

5. Inferential Statistics

5.1 Valence Analysis

# Mixed effects model for Valence
valence_model <- lmer(Valence ~ Condition * Timepoint * Context + (1|ID),
                      data = analysis_data)

cat("Mixed Effects Model: Valence ~ Condition × Timepoint × Context\n")

## Mixed Effects Model: Valence ~ Condition × Timepoint × Context

cat(rep("=", 80), "\n")

## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

print(anova(valence_model))

## Type III Analysis of Variance Table with Satterthwaite's method
##                              Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Condition                    42.705  42.705     1 309.99 13.5309 0.0002766 ***
## Timepoint                     1.313   0.656     2 307.98  0.2080 0.8123343    
## Context                     190.070  95.035     2 279.08 30.1113 1.442e-12 ***
## Condition:Timepoint          23.929  11.965     2 309.88  3.7909 0.0236282 *  
## Condition:Context            16.825   8.413     2 278.79  2.6655 0.0713392 .  
## Timepoint:Context             9.592   2.398     4 280.83  0.7598 0.5522455    
## Condition:Timepoint:Context   5.394   1.349     4 280.14  0.4273 0.7889047    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Main effect of Condition
valence_condition_emm <- emmeans(valence_model, ~ Condition)
valence_condition_contrast <- pairs(valence_condition_emm)

cat("\n\nMain Effect of Condition on Valence:\n")

## 
## 
## Main Effect of Condition on Valence:

cat(rep("-", 80), "\n")

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

print(valence_condition_emm)

##  Condition emmean    SE   df lower.CL upper.CL
##  PF          5.66 0.155 79.5     5.35     5.97
##  PS          6.42 0.167 84.3     6.09     6.75
## 
## Results are averaged over the levels of: Timepoint, Context 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

print(valence_condition_contrast)

##  contrast estimate    SE  df t.ratio p.value
##  PF - PS    -0.759 0.207 310  -3.658  0.0003
## 
## Results are averaged over the levels of: Timepoint, Context 
## Degrees-of-freedom method: kenward-roger

# Effect size
valence_condition_d <- cohens_d(Valence ~ Condition, data = analysis_data)
cat("\nCohen's d:\n")

## 
## Cohen's d:

print(valence_condition_d)

## Cohen's d |         95% CI
## --------------------------
## -0.41     | [-0.62, -0.19]
## 
## - Estimated using pooled SD.

5.2 Arousal Analysis

# Mixed effects model for Arousal
arousal_model <- lmer(Arousal ~ Condition * Timepoint * Context + (1|ID),
                      data = analysis_data)

cat("Mixed Effects Model: Arousal ~ Condition × Timepoint × Context\n")

## Mixed Effects Model: Arousal ~ Condition × Timepoint × Context

cat(rep("=", 80), "\n")

## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

print(anova(arousal_model))

## Type III Analysis of Variance Table with Satterthwaite's method
##                              Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Condition                    4.4953  4.4953     1 310.00  4.6181  0.032410 *  
## Timepoint                    5.8013  2.9006     2 308.37  2.9799  0.052263 .  
## Context                     19.2993  9.6496     2 283.10  9.9133 6.899e-05 ***
## Condition:Timepoint          0.8118  0.4059     2 309.88  0.4170  0.659402    
## Condition:Context           13.6779  6.8389     2 282.85  7.0258  0.001052 ** 
## Timepoint:Context            3.6508  0.9127     4 284.61  0.9376  0.442461    
## Condition:Timepoint:Context  0.9786  0.2447     4 284.01  0.2513  0.908713    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Main effect of Condition
arousal_condition_emm <- emmeans(arousal_model, ~ Condition)
arousal_condition_contrast <- pairs(arousal_condition_emm)

cat("\n\nMain Effect of Condition on Arousal:\n")

## 
## 
## Main Effect of Condition on Arousal:

cat(rep("-", 80), "\n")

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

print(arousal_condition_emm)

##  Condition emmean     SE   df lower.CL upper.CL
##  PF          4.28 0.0864 78.9     4.11     4.45
##  PS          4.03 0.0935 84.0     3.85     4.22
## 
## Results are averaged over the levels of: Timepoint, Context 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

print(arousal_condition_contrast)

##  contrast estimate    SE  df t.ratio p.value
##  PF - PS     0.246 0.115 310   2.137  0.0334
## 
## Results are averaged over the levels of: Timepoint, Context 
## Degrees-of-freedom method: kenward-roger

# Effect size
arousal_condition_d <- cohens_d(Arousal ~ Condition, data = analysis_data)
cat("\nCohen's d:\n")

## 
## Cohen's d:

print(arousal_condition_d)

## Cohen's d |       95% CI
## ------------------------
## 0.25      | [0.03, 0.47]
## 
## - Estimated using pooled SD.

5.3 Dominance Analysis

# Mixed effects model for Dominance
dominance_model <- lmer(Dominance ~ Condition * Timepoint * Context + (1|ID),
                        data = analysis_data)

cat("Mixed Effects Model: Dominance ~ Condition × Timepoint × Context\n")

## Mixed Effects Model: Dominance ~ Condition × Timepoint × Context

cat(rep("=", 80), "\n")

## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

print(anova(dominance_model))

## Type III Analysis of Variance Table with Satterthwaite's method
##                             Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Condition                   18.050  18.050     1 309.62 14.5787 0.0001624 ***
## Timepoint                    0.782   0.391     2 307.42  0.3159 0.7293624    
## Context                     77.113  38.556     2 282.69 31.1418 5.986e-13 ***
## Condition:Timepoint          7.216   3.608     2 309.82  2.9140 0.0557449 .  
## Condition:Context            6.139   3.069     2 282.43  2.4791 0.0856387 .  
## Timepoint:Context            3.085   0.771     4 284.45  0.6229 0.6465503    
## Condition:Timepoint:Context  2.221   0.555     4 283.75  0.4485 0.7734979    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Main effect of Condition
dominance_condition_emm <- emmeans(dominance_model, ~ Condition)
dominance_condition_contrast <- pairs(dominance_condition_emm)

cat("\n\nMain Effect of Condition on Dominance:\n")

## 
## 
## Main Effect of Condition on Dominance:

cat(rep("-", 80), "\n")

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

print(dominance_condition_emm)

##  Condition emmean    SE   df lower.CL upper.CL
##  PF          5.53 0.093 83.3     5.35     5.72
##  PS          6.02 0.101 86.0     5.82     6.22
## 
## Results are averaged over the levels of: Timepoint, Context 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

print(dominance_condition_contrast)

##  contrast estimate    SE  df t.ratio p.value
##  PF - PS     -0.49 0.129 310  -3.796  0.0002
## 
## Results are averaged over the levels of: Timepoint, Context 
## Degrees-of-freedom method: kenward-roger

# Effect size
dominance_condition_d <- cohens_d(Dominance ~ Condition, data = analysis_data)
cat("\nCohen's d:\n")

## 
## Cohen's d:

print(dominance_condition_d)

## Cohen's d |         95% CI
## --------------------------
## -0.42     | [-0.64, -0.20]
## 
## - Estimated using pooled SD.

6. Visualizations

library(ggplot2)
library(patchwork)
library(ggthemes)  # For Tableau color schemes

# Use Tableau color palette
condition_colors <- tableau_color_pal("Tableau 10")(2)
names(condition_colors) <- c("PF", "PS")

# OPTIONAL: Set default plot size for your R session (comment out if not needed)
# This provides a default size but plots will still resize with window
# options(repr.plot.width=12, repr.plot.height=4)  # For Jupyter notebooks
# If using RStudio, you can resize the Plots pane to your preference

6.1 Main Effects by Condition

# Helper function for cleaner bar plots
create_bar_plot <- function(data, y_var, y_label, title, y_limits) {
  sd_var <- paste0(sub("_M$", "_SD", y_var))
  
  ggplot(data, aes(x = Condition, y = .data[[y_var]], fill = Condition)) +
    geom_col(width = 0.6) +  # Slightly wider bars
    geom_errorbar(aes(ymin = .data[[y_var]] - .data[[sd_var]],
                      ymax = .data[[y_var]] + .data[[sd_var]]),
                  width = 0.2, linewidth = 0.6) +  # Slightly wider error bars
    scale_fill_manual(values = condition_colors) +
    scale_y_continuous(limits = y_limits, expand = expansion(mult = c(0, 0.1))) +
    labs(title = title, y = y_label, x = NULL) +
    theme_minimal(base_size = 12) +
    theme(legend.position = "none",
          plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
          panel.grid.minor = element_blank(),
          panel.grid.major.x = element_blank(),
          axis.text = element_text(size = 11),
          axis.title.y = element_text(size = 11, face = "bold"))
}

fig1a <- create_bar_plot(overall_stats, "Valence_M", "Valence", 
                         "Valence by Condition", c(0, 8))
fig1b <- create_bar_plot(overall_stats, "Arousal_M", "Arousal", 
                         "Arousal by Condition", c(0, 8))
fig1c <- create_bar_plot(overall_stats, "Dominance_M", "Dominance", 
                         "Dominance by Condition", c(0, 8))

figure1 <- fig1a + fig1b + fig1c +
  plot_annotation(
    title = "Figure 1: Main Effects of Condition on VAD Dimensions",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

print(figure1)

ggsave("figures/VAD_Analysis_Figure1_MainEffects.png", figure1,
       width = 12, height = 4, dpi = 300)

6.2 Condition × Timepoint Interaction

# Helper function for cleaner line plots
create_line_plot <- function(data, y_var, y_label, title, y_limits, legend_pos = NULL) {
  sd_var <- paste0(sub("_M$", "_SD", y_var))
  
  p <- ggplot(data, aes(x = Timepoint, y = .data[[y_var]],
                        color = Condition, group = Condition)) +
    geom_line(linewidth = 1.2) +  # Thicker lines
    geom_point(size = 3) +  # Larger points
    geom_errorbar(aes(ymin = .data[[y_var]] - .data[[sd_var]] / sqrt(n),
                      ymax = .data[[y_var]] + .data[[sd_var]] / sqrt(n)),
                  width = 0.15, linewidth = 0.6) +
    scale_color_manual(values = condition_colors) +
    scale_y_continuous(limits = y_limits) +
    labs(title = title, y = y_label, x = "Timepoint") +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
          panel.grid.minor = element_blank(),
          axis.text = element_text(size = 11),
          axis.title = element_text(size = 11, face = "bold"))
  
  if (!is.null(legend_pos)) {
    p <- p + theme(legend.position = legend_pos,
                   legend.background = element_rect(fill = "white", color = "gray80", linewidth = 0.3),
                   legend.title = element_text(size = 10, face = "bold"),
                   legend.text = element_text(size = 10),
                   legend.key.size = unit(0.8, "cm"))
  } else {
    p <- p + theme(legend.position = "none")
  }
  
  return(p)
}

# IMPROVED: Better legend placement to avoid covering data
fig2a <- create_line_plot(condition_time_stats, "Valence_M", "Valence",
                          "Valence Over Time", c(3, 7.5), legend_pos = c(0.15, 0.2))
fig2b <- create_line_plot(condition_time_stats, "Arousal_M", "Arousal",
                          "Arousal Over Time", c(3, 6))
fig2c <- create_line_plot(condition_time_stats, "Dominance_M", "Dominance",
                          "Dominance Over Time", c(3, 7))

figure2 <- fig2a + fig2b + fig2c +
  plot_annotation(
    title = "Figure 2: Condition × Timepoint Interaction",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

print(figure2)

ggsave("figures/VAD_Analysis_Figure2_TimeInteraction.png", figure2,
       width = 12, height = 4, dpi = 300)

6.3 Condition × Context Interaction

# Helper function for grouped bar plots
create_grouped_bar <- function(data, y_var, y_label, title, y_limits, legend_pos = NULL) {
  sd_var <- paste0(sub("_M$", "_SD", y_var))
  
  p <- ggplot(data, aes(x = Context, y = .data[[y_var]], fill = Condition)) +
    geom_col(position = position_dodge(0.8), width = 0.7) +  # Better spacing
    geom_errorbar(aes(ymin = .data[[y_var]] - .data[[sd_var]] / sqrt(n),
                      ymax = .data[[y_var]] + .data[[sd_var]] / sqrt(n)),
                  position = position_dodge(0.8), width = 0.2, linewidth = 0.6) +
    scale_fill_manual(values = condition_colors) +
    scale_y_continuous(limits = y_limits, expand = expansion(mult = c(0, 0.1))) +
    labs(title = title, y = y_label, x = "Context") +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
          panel.grid.minor = element_blank(),
          panel.grid.major.x = element_blank(),
          axis.text = element_text(size = 11),
          axis.title = element_text(size = 11, face = "bold"))
  
  if (!is.null(legend_pos)) {
    p <- p + theme(legend.position = legend_pos,
                   legend.background = element_rect(fill = "white", color = "gray80", linewidth = 0.3),
                   legend.title = element_text(size = 10, face = "bold"),
                   legend.text = element_text(size = 10),
                   legend.key.size = unit(0.6, "cm"))
  } else {
    p <- p + theme(legend.position = "none")
  }
  
  return(p)
}

# IMPROVED: Better legend placement
fig3a <- create_grouped_bar(condition_context_stats, "Valence_M", "Valence",
                            "Valence by Context", c(0, 7.5), legend_pos = c(0.15, 0.85))
fig3b <- create_grouped_bar(condition_context_stats, "Arousal_M", "Arousal",
                            "Arousal by Context", c(0, 6.5))
fig3c <- create_grouped_bar(condition_context_stats, "Dominance_M", "Dominance",
                            "Dominance by Context", c(0, 7))

figure3 <- fig3a + fig3b + fig3c +
  plot_annotation(
    title = "Figure 3: Condition × Context Interaction",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

print(figure3)

ggsave("figures/VAD_Analysis_Figure3_ContextInteraction.png", figure3,
       width = 12, height = 4, dpi = 300)

6.4 Emotional Space (Valence-Arousal)

# IMPROVED: Better legend placement and clearer quadrant lines
figure4 <- ggplot(analysis_data, aes(Valence, Arousal, color = Condition)) +
  geom_point(alpha = 0.20, size = 2) +
  stat_ellipse(level = 0.95, linewidth = 1.2) +
  scale_color_manual(values = condition_colors) +

  geom_vline(xintercept = 5, linetype = "dashed", alpha = 0.35) +
  geom_hline(yintercept = 5, linetype = "dashed", alpha = 0.35) +

  # Annotations in the four quadrants with proper alignment
  annotate("text", x = 0.5, y = 9.5, label = "High Arousal\nLow Valence", hjust = 0, vjust = 1, size = 4, alpha = 0.7) +
  annotate("text", x = 0.5, y = 0.5, label = "Low Arousal\nLow Valence", hjust = 0, vjust = 0, size = 4, alpha = 0.7) +
  annotate("text", x = 9.5, y = 9.5, label = "High Arousal\nHigh Valence", hjust = 1, vjust = 1, size = 4, alpha = 0.7) +
  annotate("text", x = 9.5, y = 0.5, label = "Low Arousal\nHigh Valence", hjust = 1, vjust = 0, size = 4, alpha = 0.7) +

  scale_x_continuous(expand = expansion(mult = 0.02)) +
  scale_y_continuous(expand = expansion(mult = 0.02)) +

  theme_minimal(base_size = 14) +
  theme(
    axis.title = element_text(face = "bold"),
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    legend.position = "right"
  )


print(figure4)

# Save with fixed size for publication, but display is responsive
ggsave("figures/VAD_Analysis_Figure4_EmotionalSpace.png", figure4,
       width = 8, height = 8, dpi = 300)

6.5 Three-Way Interaction Heatmaps

# Helper function for heatmaps with improved text visibility
create_heatmap <- function(data, fill_var, title, color_low, color_high, midpoint, limits) {
  ggplot(data, aes(x = Timepoint, y = Context, fill = .data[[fill_var]])) +
    geom_tile(color = "white", linewidth = 1) +  # Thicker tile borders
    # IMPROVED: Dynamic text color for better contrast
    geom_text(aes(label = sprintf("%.1f", .data[[fill_var]]),
                  color = ifelse(.data[[fill_var]] > midpoint + 0.8 | 
                                .data[[fill_var]] < midpoint - 0.8, 
                                "white", "black")),
              fontface = "bold", size = 3.5, show.legend = FALSE) +
    scale_color_identity() +
    scale_fill_gradient2(low = color_low, mid = "grey80", high = color_high,
                        midpoint = midpoint, limits = limits) +
    facet_wrap(~ Condition, ncol = 2) +
    labs(title = title, x = "Timepoint", y = "Context", fill = title) +
    theme_minimal(base_size = 11) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5, size = 10),  # Horizontal text
          axis.text.y = element_text(size = 10),
          axis.title = element_text(size = 11, face = "bold"),
          panel.grid = element_blank(),
          strip.text = element_text(size = 11, face = "bold"),
          strip.background = element_rect(fill = "grey95", color = "grey70"),
          legend.title = element_text(size = 10, face = "bold"),
          legend.text = element_text(size = 9))
}

fig5a <- create_heatmap(three_way_stats, "Valence_M", "Valence",
                        "#d7191c", "#2b83ba", 5, c(3, 7))
fig5b <- create_heatmap(three_way_stats, "Arousal_M", "Arousal",
                        "#2b83ba", "#d7191c", 4.5, c(3, 6))
fig5c <- create_heatmap(three_way_stats, "Dominance_M", "Dominance",
                        "#d7191c", "#2b83ba", 5, c(3, 7))

figure5 <- fig5a / fig5b / fig5c +
  plot_annotation(
    title = "Figure 5: Three-Way Interaction (Condition × Timepoint × Context)",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

print(figure5)

ggsave("figures/VAD_Analysis_Figure5_ThreeWayInteraction.png", figure5,
       width = 10, height = 11, dpi = 300)

6.6 Valence Trajectories by Context

# Create plots for each context with improved spacing and legend
contexts <- c("activity", "group", "social")
context_titles <- c("Activity", "Group", "Social")

plots <- lapply(1:3, function(i) {
  y_label <- if (i == 1) "Valence" else ""
  
  p <- ggplot(three_way_stats %>% filter(Context == contexts[i]),
              aes(x = Timepoint, y = Valence_M, color = Condition, group = Condition)) +
    geom_line(linewidth = 1.2) +  # Thicker lines
    geom_point(size = 3) +  # Larger points
    geom_errorbar(aes(ymin = Valence_M - Valence_SD / sqrt(n),
                      ymax = Valence_M + Valence_SD / sqrt(n)),
                  width = 0.15, linewidth = 0.6) +
    scale_color_manual(values = condition_colors) +
    scale_y_continuous(limits = c(3, 8), breaks = seq(3, 7, 1)) +  # More breaks
    labs(title = context_titles[i], y = y_label, x = "Timepoint") +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
          panel.grid.minor = element_blank(),
          axis.text = element_text(size = 10),
          axis.title = element_text(size = 11, face = "bold"))
  
  # IMPROVED: Better legend placement for middle panel
  if (i == 2) {
    p <- p + theme(legend.position = c(0.85, 0.25),  # Lower right to avoid data
                   legend.background = element_rect(fill = "white", color = "gray80", linewidth = 0.3),
                   legend.title = element_text(size = 10, face = "bold"),
                   legend.text = element_text(size = 10),
                   legend.key.size = unit(0.7, "cm"))
  } else {
    p <- p + theme(legend.position = "none")
  }
  
  return(p)
})

figure6 <- plots[[1]] + plots[[2]] + plots[[3]] +
  plot_annotation(
    title = "Figure 6: Valence Trajectories by Context",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

print(figure6)

ggsave("figures/VAD_Analysis_Figure6_ValenceTrajectories.png", figure6,
       width = 12, height = 4.5, dpi = 300)

6.7 Arousal Trajectories by Context

# Create plots for each context - Arousal
contexts <- c("activity", "group", "social")
context_titles <- c("Activity", "Group", "Social")

plots_arousal <- lapply(1:3, function(i) {
  y_label <- if (i == 1) "Arousal" else ""
  
  p <- ggplot(three_way_stats %>% filter(Context == contexts[i]),
              aes(x = Timepoint, y = Arousal_M, color = Condition, group = Condition)) +
    geom_line(linewidth = 1.2) +  # Thicker lines
    geom_point(size = 3) +  # Larger points
    geom_errorbar(aes(ymin = Arousal_M - Arousal_SD / sqrt(n),
                      ymax = Arousal_M + Arousal_SD / sqrt(n)),
                  width = 0.15, linewidth = 0.6) +
    scale_color_manual(values = condition_colors) +
    scale_y_continuous(limits = c(3, 6), breaks = seq(3, 6, 1)) +  # Arousal range
    labs(title = context_titles[i], y = y_label, x = "Timepoint") +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
          panel.grid.minor = element_blank(),
          axis.text = element_text(size = 10),
          axis.title = element_text(size = 11, face = "bold"))
  
  # Better legend placement for middle panel
  if (i == 2) {
    p <- p + theme(legend.position = c(0.85, 0.25),
                   legend.background = element_rect(fill = "white", color = "gray80", linewidth = 0.3),
                   legend.title = element_text(size = 10, face = "bold"),
                   legend.text = element_text(size = 10),
                   legend.key.size = unit(0.7, "cm"))
  } else {
    p <- p + theme(legend.position = "none")
  }
  
  return(p)
})

figure7 <- plots_arousal[[1]] + plots_arousal[[2]] + plots_arousal[[3]] +
  plot_annotation(
    title = "Figure 7: Arousal Trajectories by Context",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

print(figure7)

ggsave("figures/VAD_Analysis_Figure7_ArousalTrajectories.png", figure7,
       width = 12, height = 4.5, dpi = 300)

6.8 Dominance Trajectories by Context

# Create plots for each context - Dominance
contexts <- c("activity", "group", "social")
context_titles <- c("Activity", "Group", "Social")

plots_dominance <- lapply(1:3, function(i) {
  y_label <- if (i == 1) "Dominance" else ""
  
  p <- ggplot(three_way_stats %>% filter(Context == contexts[i]),
              aes(x = Timepoint, y = Dominance_M, color = Condition, group = Condition)) +
    geom_line(linewidth = 1.2) +  # Thicker lines
    geom_point(size = 3) +  # Larger points
    geom_errorbar(aes(ymin = Dominance_M - Dominance_SD / sqrt(n),
                      ymax = Dominance_M + Dominance_SD / sqrt(n)),
                  width = 0.15, linewidth = 0.6) +
    scale_color_manual(values = condition_colors) +
    scale_y_continuous(limits = c(3, 7.5), breaks = seq(3, 7, 1)) +  # Dominance range
    labs(title = context_titles[i], y = y_label, x = "Timepoint") +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
          panel.grid.minor = element_blank(),
          axis.text = element_text(size = 10),
          axis.title = element_text(size = 11, face = "bold"))
  
  # Better legend placement for middle panel
  if (i == 2) {
    p <- p + theme(legend.position = c(0.85, 0.25),
                   legend.background = element_rect(fill = "white", color = "gray80", linewidth = 0.3),
                   legend.title = element_text(size = 10, face = "bold"),
                   legend.text = element_text(size = 10),
                   legend.key.size = unit(0.7, "cm"))
  } else {
    p <- p + theme(legend.position = "none")
  }
  
  return(p)
})

figure8 <- plots_dominance[[1]] + plots_dominance[[2]] + plots_dominance[[3]] +
  plot_annotation(
    title = "Figure 8: Dominance Trajectories by Context",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

print(figure8)

ggsave("figures/VAD_Analysis_Figure8_DominanceTrajectories.png", figure8,
       width = 12, height = 4.5, dpi = 300)

# Save descriptive statistics
write_xlsx(list(
  Overall = overall_stats,
  ConditionByTime = condition_time_stats,
  ConditionByContext = condition_context_stats,
  ThreeWay = three_way_stats
), "tables/VAD_Analysis_DescriptiveStats.xlsx")

# Extract and save model results
valence_anova <- anova(valence_model)
arousal_anova <- anova(arousal_model)
dominance_anova <- anova(dominance_model)

# Make sure ANOVA results are data frames
valence_df <- as.data.frame(valence_anova)
arousal_df <- as.data.frame(arousal_anova)
dominance_df <- as.data.frame(dominance_anova)

# Example: assuming "Condition" is the effect of interest
# You can adjust the row name ("Condition") to match your actual ANOVA term
get_f <- function(df, term) {
  if (term %in% rownames(df)) df[term, "F value"] else NA
}
get_p <- function(df, term) {
  if (term %in% rownames(df)) df[term, "Pr(>F)"] else NA
}

# Build summary table
model_summary <- tibble(
  Model = c("Valence", "Arousal", "Dominance"),
  Condition_F = c(
    get_f(valence_df, "Condition"),
    get_f(arousal_df, "Condition"),
    get_f(dominance_df, "Condition")
  ),
  Condition_p = c(
    get_p(valence_df, "Condition"),
    get_p(arousal_df, "Condition"),
    get_p(dominance_df, "Condition")
  )
)

# Save ANOVA summaries to Excel
write_xlsx(list(
  DescriptiveStats = overall_stats,
  ConditionByTime = condition_time_stats,
  ConditionByContext = condition_context_stats,
  ThreeWay = three_way_stats,
  ANOVA_Summary = model_summary
), "tables/VAD_Analysis_Results.xlsx")

VAD Emotion Analysis: Productive Failure and Productive Success

Zack Carpenter

2025-10-28

Overview

1. Setup & Data Processing

1.1 Load Libraries

1.2 Read Raw Data

2. Emotion Word Consolidation

2.1 Define Emotion Lexicon

2.2 Consolidation Function

2.3 Apply Consolidation

3. VAD Score Assignment

3.1 Load VAD Lexicon

3.2 Manual VAD Mappings

3.3 Map VAD to Emotions

4. Descriptive Statistics

4.1 Prepare Analysis Data

4.2 Overall Descriptives

4.3 By Condition × Timepoint

4.4 By Condition × Context

4.5 Three-Way Interaction

5. Inferential Statistics

5.1 Valence Analysis

5.2 Arousal Analysis

5.3 Dominance Analysis

6. Visualizations

6.1 Main Effects by Condition

6.2 Condition × Timepoint Interaction

6.3 Condition × Context Interaction

6.4 Emotional Space (Valence-Arousal)

6.5 Three-Way Interaction Heatmaps

6.6 Valence Trajectories by Context

6.7 Arousal Trajectories by Context

6.8 Dominance Trajectories by Context