kid_expt1_main

library(pacman)
pacman::p_load(tidyverse,
               conflicted,
               here,
               lme4,
               lmerTest,
               effects,
               performance,
               effectsize,
               lsmeans,
               patchwork,
               rmarkdown,
               psych,
               ggpubr,
               ggplot2,
               ggplotly,
               tidyr,
               purrr,
               tidytable,
               stringr,
               emmeans,
               plotly
               )

Warning: package 'ggplotly' is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

Warning: 'BiocManager' not available.  Could not check Bioconductor.

Please use `install.packages('BiocManager')` and then retry.

Warning in p_install(package, character.only = TRUE, ...):

Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
logical.return = TRUE, : there is no package called 'ggplotly'

The downloaded binary packages are in
    /var/folders/92/xm82pvms01l67cph1f7lch_80000gn/T//RtmpvAO91F/downloaded_packages

tidytable installed

The downloaded binary packages are in
    /var/folders/92/xm82pvms01l67cph1f7lch_80000gn/T//RtmpvAO91F/downloaded_packages

plotly installed

Warning: package 'plotly' was built under R version 4.4.3

Warning in pacman::p_load(tidyverse, conflicted, here, lme4, lmerTest, effects, : Failed to install/load:
ggplotly, tidytable, plotly

conflict_prefer("select", "dplyr")

[conflicted] Will prefer dplyr::select over any other package.

conflict_prefer("rename", "dplyr")

[conflicted] Will prefer dplyr::rename over any other package.

conflict_prefer("lmer", "lmerTest")

[conflicted] Will prefer lmerTest::lmer over any other package.

conflict_prefer("here", "here")

[conflicted] Will prefer here::here over any other package.

conflicted::conflicts_prefer(psych::alpha)

[conflicted] Will prefer psych::alpha over any other package.

conflicts_prefer(dplyr::arrange)

[conflicted] Will prefer dplyr::arrange over any other package.

conflicts_prefer(tidyr::pivot_longer)

[conflicted] Will prefer tidyr::pivot_longer over any other package.

conflicts_prefer(tidyr::pivot_wider)

[conflicted] Will prefer tidyr::pivot_wider over any other package.

conflicts_prefer(tidytable::`%in%`)

[conflicted] Will prefer tidytable::`%in%` over any other package.

conflicts_prefer(dplyr::group_by)

[conflicted] Will prefer dplyr::group_by over any other package.

conflicts_prefer(dplyr::summarize)

[conflicted] Will prefer dplyr::summarize over any other package.

conflicts_prefer(dplyr::mutate)

[conflicted] Will prefer dplyr::mutate over any other package.

conflicts_prefer(dplyr::filter)

[conflicted] Will prefer dplyr::filter over any other package.

conflicts_prefer(dplyr::row_number)

[conflicted] Will prefer dplyr::row_number over any other package.

conflicts_prefer(tidytable::ungroup)

[conflicted] Will prefer tidytable::ungroup over any other package.

conflicts_prefer(dplyr::n)

[conflicted] Will prefer dplyr::n over any other package.

conflicts_prefer(dplyr::case_when)

[conflicted] Will prefer dplyr::case_when over any other package.

Wrangling

data_path <- "/Users/apple/Desktop/JHU/Energy_Kid_Analysis/Main/expt1_main_raw"
file_paths <- list.files(path = data_path, full.names = TRUE) 
kid_expt1_main <- tidytable::map_df(file_paths, ~ read.csv(.x) %>% 
              mutate(filename = basename(.x))) 
write.csv(kid_expt1_main, "kid_expt1_raw.csv")

vital_verbal <- read.csv("/Users/apple/Desktop/JHU/Energy_Kid_Analysis/Main/vital_verbal.csv")

Compute Vitalism score

#Judgment

vital_judge_data <- kid_expt1_main %>%
  filter(task == "vital_judge") %>%
  filter(filename != "fkk50flduj.csv") %>%  # Remove participant who didn't complete vital task
  select(filename, item, response) %>%
  mutate(response = as.numeric(response),
         response_reversed = ifelse(response == 0, 1, 0)) # Reverse: 0 (yes) → 1, 1 (no) → 0

vital_judge_data <- vital_judge_data %>%
  mutate(
    category = case_when(
      item %in% c("bear1", "bluejay1", "bee1") ~ "animal",
      item %in% c("sunflower3", "cranberry_bush3", "maple_tree3") ~ "plant",
      item %in% c("sun2", "cloud2", "rock2") ~ "natural",
      item %in% c("bicyle4", "lamp4", "pencil4") ~ "artifact",
      TRUE ~ NA_character_
    )
  )

compute_vital_judge_score <- function(df) {
  
  # Summarize by participant (filename) and category
  summary_df <- df %>%
    group_by(filename, category) %>%
    summarize(yes_count = sum(response_reversed), .groups = "drop") %>%
    pivot_wider(names_from = category, values_from = yes_count, names_prefix = "yes_")
  
  # Compute no counts (3 items per category)
  summary_df <- summary_df %>%
    mutate(
      no_natural = 3 - yes_natural,
      no_artifact = 3 - yes_artifact,
      no_plant = 3 - yes_plant
    )
  
  # Compute score based on criteria
  summary_df <- summary_df %>%
    mutate(
      judge_score = case_when(
        # 0 points: Some animals no
        yes_animal < 3 ~ 0,
        
        # 4 points: All animals yes + All plants yes + All natural no + All artifacts no
        yes_animal == 3 & 
        yes_plant == 3 & 
        no_natural == 3 & 
        no_artifact == 3 ~ 4,
        
        # 3 points: All animals yes + Some/all plants no + All natural no + All artifacts no
        yes_animal == 3 & 
        yes_plant < 3 & 
        no_natural == 3 & 
        no_artifact == 3 ~ 3,
        
        # 2 points: All animals yes + Some natural yes + All artifacts no
        yes_animal == 3 & 
        yes_natural >= 1 & 
        no_artifact == 3 ~ 2,
        
        # 1 point: All animals yes + Some artifacts yes
        yes_animal == 3 & 
        yes_artifact >= 1 ~ 1,
        
        # Default
        TRUE ~ NA_real_
      )
    )
  return(summary_df)
}
vital_judge_score <- compute_vital_judge_score(vital_judge_data)

#Verbal
vital_verbal[2:8] <- lapply(vital_verbal[2:8], as.numeric)
vital_verbal$verbal_total <- rowSums(vital_verbal[, 2:8], na.rm = FALSE)

#Combine
vital_judge_score$id <- gsub("\\.csv$", "", vital_judge_score$filename)

# Full join and select columns
vital_score_combined <- vital_verbal %>%
  select(id = participant.ID..osf., verbal_total) %>%
  full_join(
    vital_judge_score %>% select(id, judge_score),
    by = "id"
  ) %>%
  select(id, judge_score, verbal_total)

#compute final score
vital_score_combined$vital_total <- rowSums(vital_score_combined[, 2:3], na.rm = FALSE)

kid_expt1_main_filtered <- kid_expt1_main |>
  dplyr::filter(task %in% c("energy","fun","difficulty")) |> #subset to relevant trials 
  mutate(osfID = gsub("\\.csv$", "", filename))

length(unique(kid_expt1_main_filtered$filename))

[1] 64

#Import subject info (excluding 3 compcheck failers, 1 scammer)
kid_expt1_main_subjects <- read.csv("/Users/apple/Desktop/JHU/Energy_Kid_Analysis/Main/expt1_main_subjects.csv")

length(unique(kid_expt1_main_subjects$participant.ID..osf.))

[1] 60

#merge 
kid_expt1_main_full <- merge(kid_expt1_main_filtered, kid_expt1_main_subjects,
                     by.x = "osfID",
                     by.y = "participant.ID..osf.")

length(unique(kid_expt1_main_full$osfID)) #missing data file from e9ord0enzf

[1] 59

#reformat reponse
kid_expt1_main_clean <- kid_expt1_main_full |> 
  select(osfID, task,response, activity,X,Gender,Month,Age,Age..in.years.,task.Set) |>
   mutate(response = as.numeric(response)) |>
  rename(
    Age_in_years = Age..in.years.,
    subjectID = X,
    Age_in_months = Month
  )

#label activity type
kid_expt1_main_labeled <- kid_expt1_main_clean %>%
  mutate(act_type = case_when(
    str_detect(activity, "run") | str_detect(activity, "clean") | str_detect(activity, "hike") | str_detect(activity, "backpack") | str_detect(activity, "climb") | str_detect(activity, "swim") ~ "physical",
    str_detect(activity, "handwriting") | str_detect(activity, "math") | str_detect(activity, "music") | str_detect(activity, "homework") | str_detect(activity, "flashcards") | str_detect(activity, "alphabet") ~ "mental",
    TRUE ~ "restorative"
  ))

kid_expt1_main_labeled_wide <- kid_expt1_main_labeled |>
  pivot_wider(
    id_cols = c(osfID, subjectID, activity, act_type, 
                Gender, Age_in_months, Age, Age_in_years, task.Set),
    names_from = task,
    values_from = response
  )

#ids_in_data <- unique(kid_expt1_main_filtered$osfID)
#ids_in_subjects <- unique(kid_expt1_main_subjects$participant.ID..osf.)
#setdiff(ids_in_data, ids_in_subjects) #missing ids

A total of 64 recorded responses. Three failed comprehension checks, 1 was confirmed scammer, and 1 was test (cn5ryt825h).

Demographic

kid_expt1_main_labeled_wide %>%
  distinct(subjectID, Age_in_years) %>%
  count(Age_in_years)

# A tibble: 4 × 2
  Age_in_years     n
         <int> <int>
1            4    14
2            5    16
3            6    15
4            7    14

Adult vs. Kid Energy Rating - not correlated

#loading adult data
#energy_adult <- read.csv("/Users/apple/Desktop/JHU/Energy_Adult_Analysis/mean_energy_adult.csv")

Visualization

#Energy by activity
ggplot(kid_expt1_main_labeled_wide, 
       aes(y = fct_reorder(activity, energy, .fun = median),
           x=energy,
           fill=act_type)) + 
  stat_summary(fun = median)+
#  stat_boxplot(geom='errorbar', linetype=1, width=0.5) +
  geom_boxplot() +
#  facet_wrap(~Age_in_years) +
  theme_minimal() +
  scale_fill_manual(values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  labs(title="Energy Change by Activity (Kids)",y="Activity",x="Energy Change")+
    theme(plot.title = element_text(size=40))+
    theme(
      axis.text.y = element_text(size = 17),
      axis.text.x = element_text(size = 15), 
    axis.title.x = element_text(size = 18),
axis.title.y = element_text(size = 18), 
legend.text = element_text(size = 16),
legend.title = element_text(size = 16),
    legend.key.size = unit(2, "lines"))

Warning: Removed 18 rows containing missing values or values outside the scale range
(`geom_segment()`).

#Energy by activity type and age
kid_expt1_main_labeled_wide$act_type <- factor(
  kid_expt1_main_labeled_wide$act_type,
  levels = c("restorative", "physical", "mental")
)

ggplot(kid_expt1_main_labeled_wide, 
           aes(y=energy, x=as.factor(Age_in_years), fill=act_type)) + 
  geom_boxplot() +
  stat_summary(
  fun = mean, 
  geom = "point", 
  shape = 19, 
  size = 2, 
  color = "black",
  aes(group = Age_in_years)
) +
stat_summary(
  fun.data = mean_cl_normal, 
  geom = "errorbar", 
  width = 0.2, 
  color = "black", 
  linewidth = 0.8,
  aes(group = Age_in_years)
)+
  facet_wrap(~act_type) +
  theme_minimal(base_size = 16) +
  scale_fill_manual(values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  labs(title="Energy Change by Age",y="Energy Rating", x="Age")+
  theme(legend.position = "none")

# Difficulty by activity
ggplot(kid_expt1_main_labeled_wide,
       aes(y = fct_reorder(activity, difficulty, .fun = median),
           x=difficulty,
           fill=act_type)) + 
  stat_summary(fun = median)+
  geom_boxplot() +
  #facet_wrap(~age_in_years) +
  theme_minimal() +
  scale_fill_manual(values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  labs(title="Difficulty by Activity",y="Activity",x="Difficulty")

Warning: Removed 18 rows containing missing values or values outside the scale range
(`geom_segment()`).

# Difficulty by activity type
ggplot(kid_expt1_main_labeled_wide, 
       aes(x=difficulty, y=act_type, fill=act_type)) + 
  stat_summary(fun = median)+
  geom_boxplot() +
#  facet_wrap(~age_in_years) +
  theme_minimal() +
  scale_fill_manual(values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  labs(title="Difficulty by Activity Type",y="Activity",x="Difficulty")

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

#Fun by activity
ggplot(kid_expt1_main_labeled_wide, 
       aes(y = fct_reorder(activity, fun, .fun = median),
           x=fun,
           fill=act_type)) + 
  stat_summary(fun = median)+
  geom_boxplot() +
 # facet_wrap(~age_in_years) +
  theme_minimal() +
  scale_fill_manual(values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  labs(title="Fun by Activity",y="Activity",x="Fun")

Warning: Removed 18 rows containing missing values or values outside the scale range
(`geom_segment()`).

#Fun by activity type
ggplot(kid_expt1_main_labeled_wide, 
       aes(x=fun, y=act_type, fill=act_type)) + 
  stat_summary(fun = median)+
  geom_boxplot() +
  facet_wrap(~Age_in_years) +
  theme_minimal() +
  scale_fill_manual(values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  labs(title="Fun by Activity",y="Activity",x="Fun")

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).
Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).
Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Kids dont’ generally find restorative activities as energy-restoring, in fact it is similar to mental activities - because they involve similar amount of movement? physical activities are energy-depleting because they involve more movement?

#pretty plot
p <- ggplot(kid_expt1_main_labeled_wide, aes(y = reorder(activity, energy, FUN = median),x=energy,fill=act_type)) +
  geom_boxplot(outlier.shape = NA, width = 0.6) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, 
               fill = "black", color = "black",
               aes(group = activity)) +  # diamond for mean
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", 
               width = 0.2, color = "black", linewidth = 0.8,
               aes(group = activity)) +  # 95% CI
  scale_fill_manual(
    name = "Action Type",  
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
#  geom_jitter(width = 0.1, alpha = 0.3, size = 1) +
  theme_classic() +
  scale_x_continuous(
    breaks = c(-50, -25, 0, 25, 50),
   labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic")
  )+
  theme(
    axis.title = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(margin = margin(r = 0)),  # Add this line
    axis.text = element_text(size = 17),
    axis.text.x = element_text(angle = 15, hjust = 0.6),
    panel.grid.major.x = element_line(color = "gray80"),
    legend.title = element_text(size = 16, face = "bold"),
    legend.text = element_text(size = 14),
    theme(axis.text.x = element_text(lineheight = 0.1)) #secondary label spacing

  ) +
  labs(
    x = "Energy Rating",
    y = "Action",
  )

p

ggsave("experiment1_kids_ci.png", width = 14, height = 7, dpi = 300)

#pretty plot by activity type
ggplot(kid_expt1_main_labeled_wide, aes(x = energy, y = act_type, fill = act_type)) +
  geom_boxplot(outlier.shape = NA, width = 0.6) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, 
               fill = "black", color = "black",
               aes(group = act_type)) +  # diamond for mean
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", 
               width = 0.2, color = "black", linewidth = 0.8,
               aes(group = act_type)) +  # 95% CI
  scale_fill_manual(
    name = "Action Type",  
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
#  geom_jitter(width = 0.1, alpha = 0.3, size = 1) +
  facet_wrap(~Age_in_years) +
  theme_classic() +
  scale_x_continuous(
    breaks = c(-50, -25, 0, 25, 50),
   labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic")
  )+
  theme(
    axis.title = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(margin = margin(r = 0)),  # Add this line
    axis.text = element_text(size = 17),
    axis.text.x = element_text(angle = 15, hjust = 0.6),
    panel.grid.major.x = element_line(color = "gray80"),
    legend.title = element_text(size = 16, face = "bold"),
    legend.text = element_text(size = 14),
    theme(axis.text.x = element_text(lineheight = 0.1)) #secondary label spacing

  ) +
  labs(
    x = "Energy Rating",
    y = "Action",
  )

Kid & Adult comparison

# Prepare adult data
adult_data <- read.csv("/Users/apple/Desktop/JHU/Energy_Adult_Analysis/adult_expt1_final.csv")

adult_data <- adult_data |>
  mutate(
    group = "Adult",
    unified_activity = case_when(
      activity == "practice violin" ~ "practice violin / music*",
      activity == "jump up and down" ~ "jump up and down / hold a backpack full of books*",
      TRUE ~ activity
    )
  )

# Prepare child data
child_data <- kid_expt1_main_labeled_wide |>
  mutate(
    group = "Child",
    unified_activity = case_when(
      activity == "practice music" ~ "practice violin / music*",
      activity == "hold a backpack full of books" ~ "jump up and down / hold a backpack full of books*",
      activity == "do their homework" ~ "do homework",
      activity == "watch their favorite show on the couch" ~ "watch favorite show on the couch",
      activity == "listen to their favorite stories" ~ "listen to favorite stories",
      TRUE ~ activity
    )
  )

combined_data <- bind_rows(adult_data, child_data)

# Set activity order (by median of adult data)
activity_order <- combined_data |>
  filter(group == "Adult") |>
  group_by(unified_activity) |>
  summarise(median_val = median(energy, na.rm = TRUE)) |>
  arrange(median_val) |>
  pull(unified_activity)

combined_data$unified_activity <- factor(combined_data$unified_activity, levels = activity_order)

combined_data <- combined_data |>
  mutate(group = factor(group, levels = c("Child", "Adult"))) #change order of boxplot display

# Create the plot
p <- ggplot(combined_data, aes(x = energy, y = unified_activity, fill = act_type, alpha=group)) +
  geom_boxplot(
    position = position_dodge(width = 0.8),
    width = 0.7,
    outlier.shape = NA
  ) +
 stat_summary(
  fun = mean, 
  geom = "point", 
  shape = 19, 
  size = 2, 
  color = "black",
  aes(group = group),
  position = position_dodge(width = 0.8)
) +
stat_summary(
  fun.data = mean_cl_normal, 
  geom = "errorbar", 
  width = 0.2, 
  color = "black", 
  linewidth = 0.8,
  aes(group = group),
  position = position_dodge(width = 0.8)
)+
   scale_fill_manual(
    name = "Action Type",
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c")
  ) +
  scale_alpha_manual(
    name = "Group",
    values = c("Adult" = 1, "Child" = 0.4)
  ) +
  scale_x_continuous(
    limits = c(-50, 50),
    breaks = c(-50, -25, 0, 25, 50),
    labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral/Just okay", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic")
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.text = element_text(size = 14),
    axis.text.x = element_text(angle = 0, hjust = 0.6),
    panel.grid.major.x = element_line(color = "gray80"),
    legend.position = "top",
    legend.title = element_text(size = 16, face = "bold"),
    legend.text = element_text(size = 14),
    plot.caption = element_text(hjust = 0, size = 12, face = "italic"),
    theme(axis.text.x = element_text(lineheight = 0.1))) + #secondary label spacing
  labs(
    x = "Energy Rating",
    y = "Action",
    caption = "* Action differed slightly between adults and children"
  )+
  theme(
  plot.margin = margin(t = 10, r = 50, b = 10, l = 10)  # increase right margin
)

p

ggsave("experiment1_comparison_indiv_action.png", width = 14, height = 10)

ggplot(combined_data, aes(x = energy, y = act_type, fill = act_type, alpha=group)) +
  geom_boxplot(
    position = position_dodge(width = 0.8),
    width = 0.7,
    outlier.shape = NA
  ) +
  scale_y_discrete(limits = c("physical", "mental", "restorative")) +
 stat_summary(
  fun = mean, 
  geom = "point", 
  shape = 19, 
  size = 2, 
  color = "black",
  aes(group = group),
  position = position_dodge(width = 0.8)
) +
stat_summary(
  fun.data = mean_cl_normal, 
  geom = "errorbar", 
  width = 0.2, 
  color = "black", 
  linewidth = 0.8,
  aes(group = group),
  position = position_dodge(width = 0.8)
)+
   scale_fill_manual(
    name = "Action Type",
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c")
  ) +
  scale_alpha_manual(
    name = "Group",
    values = c("Adult" = 1, "Child" = 0.4)
  ) +
  scale_x_continuous(
    limits = c(-50, 50),
    breaks = c(-50, -25, 0, 25, 50),
    labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral/Just okay", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic")
  ) +
   theme_classic() +
  theme(
    axis.title = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.text = element_text(size = 14),
    axis.text.x = element_text(angle = 0, hjust = 0.6),
    panel.grid.major.x = element_line(color = "gray80"),
    legend.position = "top",
    legend.title = element_text(size = 16, face = "bold"),
    legend.text = element_text(size = 14),
    plot.caption = element_text(hjust = 0, size = 12, face = "italic"),
    theme(axis.text.x = element_text(lineheight = 0.1))) + #secondary label spacing
  labs(
    x = "Energy Rating",
    y = "Action Type"
  )+
    theme(
  plot.margin = margin(t = 10, r = 50, b = 10, l = 10)  # increase right margin
)

ggsave("experiment1_comparison_act_type.png", width = 10, height = 6)

ggsave("experiment1_comparison_act_type_adult.png", width = 10, height = 6)

p_adult <- ggplot(adult_data, aes(x = energy, y = reorder(activity, energy, FUN = median), fill = act_type)) +
  geom_boxplot(
    width = 0.7,
    outlier.shape = NA
  ) +
 stat_summary(
  fun = mean, 
  geom = "point", 
  shape = 19, 
  size = 2, 
  color = "black"
) +
stat_summary(
  fun.data = mean_cl_normal, 
  geom = "errorbar", 
  width = 0.2, 
  color = "black", 
  linewidth = 0.8
)+
   scale_fill_manual(
    name = "Action Type",
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c")
  ) +
  scale_x_continuous(
    limits = c(-50, 50),
    breaks = c(-50, -25, 0, 25, 50),
    labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic")
  ) +
  theme_classic() +
  theme(
    axis.title = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.text = element_text(size = 14),
    axis.text.x = element_text(angle = 0, hjust = 0.6),
    panel.grid.major.x = element_line(color = "gray80"),
    legend.position = "top",
    legend.title = element_text(size = 16, face = "bold"),
    legend.text = element_text(size = 14),
    plot.caption = element_text(hjust = 0, size = 12, face = "italic"),
    theme(axis.text.x = element_text(lineheight = 0.1))) + #secondary label spacing
  labs(
    x = "Energy Rating",
    y = "Action"
  )+
  theme(
  plot.margin = margin(t = 10, r = 50, b = 10, l = 10)  # increase right margin
)

p_adult

ggsave("experiment1_adult_indiv_action.png", width = 14, height = 10)

ggplot(adult_data, aes(x = energy, y = act_type, fill = act_type)) +
  geom_boxplot(
    width = 0.7,
    outlier.shape = NA
  ) +
    scale_y_discrete(limits = c("physical", "mental", "restorative")) +
 stat_summary(
  fun = mean, 
  geom = "point", 
  shape = 19, 
  size = 2, 
  color = "black"
) +
stat_summary(
  fun.data = mean_cl_normal, 
  geom = "errorbar", 
  width = 0.2, 
  color = "black", 
  linewidth = 0.8
)+
   scale_fill_manual(
    name = "Action Type",
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c")
  ) +
  scale_x_continuous(
    limits = c(-50, 50),
    breaks = c(-50, -25, 0, 25, 50),
    labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic")
  ) +
   theme_classic() +
  theme(
    axis.title = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.text = element_text(size = 14),
    axis.text.x = element_text(angle = 0, hjust = 0.6),
    panel.grid.major.x = element_line(color = "gray80"),
    legend.position = "top",
    legend.title = element_text(size = 16, face = "bold"),
    legend.text = element_text(size = 14),
    plot.caption = element_text(hjust = 0, size = 12, face = "italic"),
    theme(axis.text.x = element_text(lineheight = 0.1))) + #secondary label spacing
  labs(
    x = "Energy Rating",
    y = "Action Type"
  )+
    theme(
  plot.margin = margin(t = 10, r = 50, b = 10, l = 10)  # increase right margin
)

ggsave("experiment1_adult_act_type.png", width = 10, height = 6)

Correlation

library(ggcorrplot)

kid_expt1_main_wide <- kid_expt1_main_labeled |>
  pivot_wider(
    id_cols = c(subjectID, activity, act_type, Age, Age_in_years),
    names_from = task,     
    values_from = response
  )

cor_wide <- kid_expt1_main_wide |>
  select(subjectID, activity, energy) |>
  pivot_wider(
    names_from = activity,
    values_from = energy
  )

activity_matrix <- cor_wide |>
  select(-subjectID) |>
  t()

activity_order <- kid_expt1_main_wide |>
  select(activity, act_type) |>
  distinct() |>
  arrange(act_type, activity) |>
  pull(activity)

# Reorder RSM
#rsm <- rsm[activity_order, activity_order]
#ggcorrplot(rsm)

# Function to compute RSM
get_rsm <- function(data, activity_order) {
  cor_wide <- data |>
    select(subjectID, activity, energy) |>
    pivot_wider(names_from = activity, values_from = energy)
  
  rsm <- cor(cor_wide |> select(-subjectID), use = "pairwise.complete.obs")
  rsm[activity_order, activity_order]
}

# Create individual plots
p1 <- ggcorrplot(get_rsm(filter(kid_expt1_main_wide, Age_in_years == 4), activity_order)) +
  labs(title = "4 years")

Warning in cor(select(cor_wide, -subjectID), use = "pairwise.complete.obs"):
the standard deviation is zero

p2 <- ggcorrplot(get_rsm(filter(kid_expt1_main_wide, Age_in_years == 5), activity_order)) +
  labs(title = "5 years")

p3 <- ggcorrplot(get_rsm(filter(kid_expt1_main_wide, Age_in_years == 6), activity_order)) +
  labs(title = "6 years")

p4 <- ggcorrplot(get_rsm(filter(kid_expt1_main_wide, Age_in_years == 7), activity_order)) +
  labs(title = "7 years")

# Combine
(p1 + p2) / (p3 + p4)

corr_phys_older <- kid_expt1_main_labeled_wide %>%
    filter(act_type == "physical" & Age_in_years >=6) %>%
    select(energy, difficulty, fun) %>%
    cor(use = "pairwise.complete.obs")

corr_ment_older <- kid_expt1_main_labeled_wide %>%
    filter(act_type == "mental" & Age_in_years >=6) %>%
    select(energy, difficulty, fun) %>%
    cor(use = "pairwise.complete.obs")

corr_rest_older <- kid_expt1_main_labeled_wide %>%
    filter(act_type == "restorative" & Age_in_years >=6) %>%
    select(energy, difficulty, fun) %>%
    cor(use = "pairwise.complete.obs")

fig_phys_older <- ggcorrplot(corr_phys_older, lab = TRUE, lab_col = "black", title = "Physical Activities")
fig_ment_older <- ggcorrplot(corr_ment_older, lab = TRUE, lab_col = "black", title = "Mental Activities")
fig_rest_older <- ggcorrplot(corr_rest_older, lab = TRUE, lab_col = "black", title = "Restorative Activities")

fig_phys_older

fig_ment_older

fig_rest_older

younger kids (3-5): the more difficult, the less energy-depleting/more restorative for all activities (esp mental ones). For mental activities only, the more difficult more fun.

older kids (6-7): Like younger kids, for mental activities, they say the more difficult the less energy-depleting. Like adults (unlike younger kids), they say the less difficult the more fun for mental and physical activities, and the more energy depleting the less fun for restorative activities.

kids below 5 have trouble with difficulty sliders (numbers??)

Modeling

library(lme4)
library(lmerTest)
#comparing across types
kid_expt1_main_labeled$act_type <- factor(kid_expt1_main_labeled$act_type, levels = c("physical", "mental", "restorative"))

#Primary H
phys_energy <- kid_expt1_main_labeled |> filter(act_type=="physical" & task=="energy")
summary(lmer(formula = response ~ 1 + (1 |subjectID), data=filter(phys_energy)))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: response ~ 1 + (1 | subjectID)
   Data: filter(phys_energy)

REML criterion at convergence: 2334.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9466 -0.6818 -0.1220  0.7054  2.2060 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept) 316.0    17.78   
 Residual              956.9    30.93   
Number of obs: 236, groups:  subjectID, 59

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)   
(Intercept)   -9.309      3.068 58.000  -3.035   0.0036 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ment_energy <- kid_expt1_main_labeled |> filter(act_type=="mental" & task=="energy")
summary(lmer(formula = response ~ 1 + (1 |subjectID), data=filter(ment_energy))) # not different from 0

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: response ~ 1 + (1 | subjectID)
   Data: filter(ment_energy)

REML criterion at convergence: 2287.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.42225 -0.65831 -0.04397  0.72149  2.18744 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept) 314.2    17.72   
 Residual              760.3    27.57   
Number of obs: 236, groups:  subjectID, 59

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)
(Intercept)    3.216      2.923 58.000     1.1    0.276

rest_energy <- kid_expt1_main_labeled |> filter(act_type=="restorative" & task=="energy")
summary(lmer(formula = response ~ 1 + (1 |subjectID), data=rest_energy))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: response ~ 1 + (1 | subjectID)
   Data: rest_energy

REML criterion at convergence: 2338.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9538 -0.7463  0.2411  0.7576  1.6438 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept)  173.1   13.16   
 Residual              1061.0   32.57   
Number of obs: 236, groups:  subjectID, 59

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   10.466      2.726 58.000    3.84 0.000307 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model1_primary <- lmer(energy ~ act_type + (1|subjectID), data=kid_expt1_main_labeled_wide)
emmeans(model1_primary, ~ act_type) 

 act_type    emmean   SE  df lower.CL upper.CL
 restorative  10.47 2.71 151     5.12    15.82
 physical     -9.31 2.71 151   -14.66    -3.96
 mental        3.22 2.71 151    -2.13     8.57

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

#Exploratory act_type*age
kid_expt1_main_labeled_wide$Age_in_years <- as.factor(kid_expt1_main_labeled_wide$Age_in_years)

model1 <- lmer(energy ~ act_type*Age_in_years + (1|subjectID), data=kid_expt1_main_labeled_wide)
summary(model1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ act_type * Age_in_years + (1 | subjectID)
   Data: kid_expt1_main_labeled_wide

REML criterion at convergence: 6900.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.51868 -0.69214 -0.03846  0.77517  2.23224 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept)  178.8   13.37   
 Residual              1008.2   31.75   
Number of obs: 708, groups:  subjectID, 59

Fixed effects:
                               Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)                     14.5357     5.5475 142.7961   2.620  0.00974 **
act_typephysical               -16.4643     6.0005 641.0000  -2.744  0.00624 **
act_typemental                  -5.8393     6.0005 641.0000  -0.973  0.33085   
Age_in_years5                   -3.1138     7.5963 142.7961  -0.410  0.68248   
Age_in_years6                   -4.5690     7.7135 142.7961  -0.592  0.55456   
Age_in_years7                   -8.6964     7.8454 142.7961  -1.108  0.26952   
act_typephysical:Age_in_years5  -6.0357     8.2165 641.0000  -0.735  0.46286   
act_typemental:Age_in_years5    -0.7232     8.2165 641.0000  -0.088  0.92989   
act_typephysical:Age_in_years6  -8.9357     8.3433 641.0000  -1.071  0.28457   
act_typemental:Age_in_years6     3.7726     8.3433 641.0000   0.452  0.65130   
act_typephysical:Age_in_years7   2.5179     8.4860 641.0000   0.297  0.76678   
act_typemental:Age_in_years7    -9.1607     8.4860 641.0000  -1.080  0.28076   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
              (Intr) act_typp act_typm Ag_n_5 Ag_n_6 Ag_n_7 act_typp:A__5
act_typphys   -0.541                                                     
act_typmntl   -0.541  0.500                                              
Age_in_yrs5   -0.730  0.395    0.395                                     
Age_in_yrs6   -0.719  0.389    0.389    0.525                            
Age_in_yrs7   -0.707  0.382    0.382    0.516  0.509                     
act_typp:A__5  0.395 -0.730   -0.365   -0.541 -0.284 -0.279              
act_typm:A__5  0.395 -0.365   -0.730   -0.541 -0.284 -0.279  0.500       
act_typp:A__6  0.389 -0.719   -0.360   -0.284 -0.541 -0.275  0.525       
act_typm:A__6  0.389 -0.360   -0.719   -0.284 -0.541 -0.275  0.263       
act_typp:A__7  0.382 -0.707   -0.354   -0.279 -0.275 -0.541  0.516       
act_typm:A__7  0.382 -0.354   -0.707   -0.279 -0.275 -0.541  0.258       
              act_typm:A__5 act_typp:A__6 act_typm:A__6 act_typp:A__7
act_typphys                                                          
act_typmntl                                                          
Age_in_yrs5                                                          
Age_in_yrs6                                                          
Age_in_yrs7                                                          
act_typp:A__5                                                        
act_typm:A__5                                                        
act_typp:A__6  0.263                                                 
act_typm:A__6  0.525         0.500                                   
act_typp:A__7  0.258         0.509         0.254                     
act_typm:A__7  0.516         0.254         0.509         0.500       

plot(allEffects(model1))

kid_expt1_main_labeled_wide$act_type <- factor(kid_expt1_main_labeled_wide$act_type, levels = c("physical", "mental", "restorative"))

model1_kids <- lmer(energy ~ difficulty+fun + (1|subjectID), data=kid_expt1_main_labeled_wide)
model2_kids <- lmer(energy ~ act_type + difficulty+fun + (1|subjectID), data=kid_expt1_main_labeled_wide)
model3_kids <- lmer(energy ~ Age*act_type + difficulty+fun + (1|subjectID), data=kid_expt1_main_labeled_wide)


summary(model1_kids)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ difficulty + fun + (1 | subjectID)
   Data: kid_expt1_main_labeled_wide

REML criterion at convergence: 6977.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.48551 -0.75489 -0.04874  0.77680  2.36598 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept)  142.3   11.93   
 Residual              1022.4   31.97   
Number of obs: 708, groups:  subjectID, 59

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)   1.52047    3.25444 296.30707   0.467    0.641    
difficulty   -0.17396    0.03475 704.42266  -5.006 7.02e-07 ***
fun           0.13928    0.03462 702.10827   4.023 6.38e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) dffclt
difficulty -0.564       
fun        -0.654  0.174

summary(model2_kids)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ act_type + difficulty + fun + (1 | subjectID)
   Data: kid_expt1_main_labeled_wide

REML criterion at convergence: 6944

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.66060 -0.72024 -0.05817  0.75896  2.24779 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept) 148.3    12.18   
 Residual              985.0    31.39   
Number of obs: 708, groups:  subjectID, 59

Fixed effects:
                     Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)          -8.47434    3.79461 423.95214  -2.233 0.026053 *  
act_typemental       11.67769    2.90758 645.48896   4.016 6.61e-05 ***
act_typerestorative  14.42913    3.07694 657.64659   4.689 3.33e-06 ***
difficulty           -0.12664    0.03588 702.94396  -3.530 0.000443 ***
fun                   0.12508    0.03456 700.08542   3.620 0.000316 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) act_typm act_typr dffclt
act_typmntl -0.452                         
act_typrstr -0.458  0.496                  
difficulty  -0.579  0.110    0.301         
fun         -0.514  0.036   -0.129    0.119

summary(model3_kids)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ Age * act_type + difficulty + fun + (1 | subjectID)
   Data: kid_expt1_main_labeled_wide

REML criterion at convergence: 6931.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.65823 -0.72116 -0.05724  0.73387  2.15983 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept) 147.9    12.16   
 Residual              986.3    31.41   
Number of obs: 708, groups:  subjectID, 59

Fixed effects:
                         Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)              -3.50681   14.44392 165.44838  -0.243 0.808471    
Age                      -0.75616    2.32400 158.52852  -0.325 0.745331    
act_typemental           27.70072   15.97764 642.91547   1.734 0.083448 .  
act_typerestorative      19.77513   15.93925 642.53801   1.241 0.215186    
difficulty               -0.13060    0.03607 699.98763  -3.621 0.000315 ***
fun                       0.12131    0.03468 697.01804   3.499 0.000498 ***
Age:act_typemental       -2.66339    2.60996 643.24255  -1.020 0.307888    
Age:act_typerestorative  -0.89511    2.60296 642.89829  -0.344 0.731048    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Age    act_typm act_typr dffclt fun    Ag:ct_typm
Age         -0.965                                                  
act_typmntl -0.536  0.550                                           
act_typrstr -0.544  0.549  0.501                                    
difficulty  -0.113 -0.042 -0.067   -0.012                           
fun         -0.136  0.000 -0.056   -0.044    0.124                  
Ag:ct_typmn  0.523 -0.559 -0.983   -0.492    0.089  0.063           
Ag:ct_typrs  0.531 -0.560 -0.494   -0.981    0.072  0.020  0.503    

r2(model1_kids)

# R2 for Mixed Models

  Conditional R2: 0.180
     Marginal R2: 0.066

r2(model2_kids)

# R2 for Mixed Models

  Conditional R2: 0.210
     Marginal R2: 0.091

anova(model2_kids,model3_kids)

refitting model(s) with ML (instead of REML)

Data: kid_expt1_main_labeled_wide
Models:
model2_kids: energy ~ act_type + difficulty + fun + (1 | subjectID)
model3_kids: energy ~ Age * act_type + difficulty + fun + (1 | subjectID)
            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
model2_kids    7 6959.1 6991.0 -3472.6   6945.1                     
model3_kids   10 6962.8 7008.4 -3471.4   6942.8 2.3082  3     0.5109

Explorary Analysis

correlation matrix (energy, fun, difficulty)

# Correlation heatmaps by activity
library(ggcorrplot)

corr_phys <- kid_expt1_main_labeled_wide %>%
    filter(act_type == "physical") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

corr_ment <- kid_expt1_main_labeled_wide %>%
    filter(act_type == "mental") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

corr_rest <- kid_expt1_main_labeled_wide %>%
    filter(act_type == "restorative") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

fig_phys <- ggcorrplot(corr_phys, lab = TRUE, lab_col = "black", title = "Physical Activities")
fig_ment <- ggcorrplot(corr_ment, lab = TRUE, lab_col = "black", title = "Mental Activities")
fig_rest <- ggcorrplot(corr_rest, lab = TRUE, lab_col = "black", title = "Restorative Activities")

fig_phys

fig_ment

fig_rest

Vitalism

kid_expt1_main_fullvital <- full_join(
  vital_score_combined,
  kid_expt1_main_labeled_wide,
  by = c("id" = "osfID")
)
kid_expt1_main_fullvital <- kid_expt1_main_fullvital |>
  filter(!is.na(vital_total)) |>
  filter(!is.na(energy)) #exclude NAs, "cn5ryt825h", "elyvcsz1xu", "k65b39fgz6", "ahb8sd3jmx"

kid_expt1_main_fullvital$vital_total <- as.numeric(kid_expt1_main_fullvital$vital_total)

# is vitalism score related to age?
ggplot(kid_expt1_main_fullvital, aes(x = Age, y = vital_total, color=as.factor(Age_in_years))) +
  geom_point()

cor(kid_expt1_main_fullvital$vital_total, kid_expt1_main_fullvital$Age, use = "complete.obs")

[1] 0.4360212

# group by vital score
vital_score_combined |>
  group_by(vital_total) |>
  summarize(n =n())

# A tibble: 17 × 2
   vital_total     n
         <dbl> <int>
 1           0     5
 2           1     1
 3           2     4
 4           3     4
 5           4     7
 6           5     7
 7           6     7
 8           7     5
 9           8     4
10           9     6
11          10     1
12          11     2
13          12     3
14          13     2
15          14     1
16          18     1
17          NA     5

median(vital_score_combined$vital_total, na.rm = TRUE) #median = 6

[1] 6

kid_expt1_main_fullvital$vital_median_split <- cut(
  kid_expt1_main_fullvital$vital_total,
  breaks = c(-Inf, 5, Inf),
  labels = c("<6", ">=6"), # 28 in <6, 32 in >=6
  include.lowest = TRUE
)

kid_expt1_main_fullvital |>
  group_by(vital_median_split) |>
  summarize(sd = sd(energy)) # vitalistic group showed less variability in rating

# A tibble: 2 × 2
  vital_median_split    sd
  <fct>              <dbl>
1 <6                  38.3
2 >=6                 32.5

#plot energy rating by vitalism median split
ggplot(kid_expt1_main_fullvital, aes(x = energy, y = act_type, fill = act_type)) +
  geom_boxplot(outlier.shape = NA, width = 0.6) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, 
               fill = "black", color = "black",
               aes(group = act_type)) +  # diamond for mean
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", 
               width = 0.2, color = "black", linewidth = 0.8,
               aes(group = act_type)) +  # 95% CI
  scale_fill_manual(
    name = "Action Type",  
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  facet_wrap(~vital_median_split) +
  theme_classic() +
  scale_x_continuous(
    breaks = c(-50, -25, 0, 25, 50),
   labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic"))+
    theme_minimal()+
   theme(axis.text.x = element_text(angle = 15, hjust = 1)) +
  labs(x = "Energy Rating", y = "Action")

possible range -6 ~ 19

actual range 0 ~ 18

#difficulty rating by vitalism median split
ggplot(kid_expt1_main_fullvital, aes(x = difficulty, y = act_type, fill = act_type)) +
  geom_boxplot(outlier.shape = NA, width = 0.6) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, 
               fill = "black", color = "black",
               aes(group = act_type)) +  # diamond for mean
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", 
               width = 0.2, color = "black", linewidth = 0.8,
               aes(group = act_type)) +  # 95% CI
  scale_fill_manual(
    name = "Action Type",  
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  facet_wrap(~vital_median_split) +
  theme_classic() +
  scale_x_continuous(
    breaks = c(0, 25, 50,75,100),
   labels = c("0\nExtremely easy", 
             "25\nSomewhat easy", 
             "50\nNeutral", 
             "75\nSomewhat hard", 
             "100\nExtremely hard"))+
    theme_minimal()+
   theme(axis.text.x = element_text(angle = 15, hjust = 1)) +
  labs(x = "Difficulty Rating", y = "Action")

#fun rating by vitalism median split
ggplot(kid_expt1_main_fullvital, aes(x = fun, y = act_type, fill = act_type)) +
  geom_boxplot(outlier.shape = NA, width = 0.6) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, 
               fill = "black", color = "black",
               aes(group = act_type)) +  # diamond for mean
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", 
               width = 0.2, color = "black", linewidth = 0.8,
               aes(group = act_type)) +  # 95% CI
  scale_fill_manual(
    name = "Action Type",  
    values = c("mental" = "#1f77b4", "physical" = "#ff7f0e", "restorative" = "#2ca02c"))+
  facet_wrap(~vital_median_split) +
  theme_classic() +
  scale_x_continuous(
    breaks = c(0, 25, 50,75,100),
   labels = c("0\nExtremely boring", 
             "25\nSomewhat boring", 
             "50\nNeutral", 
             "75\nSomewhat fun", 
             "100\nExtremely fun"))+
    theme_minimal()+
   theme(axis.text.x = element_text(angle = 15, hjust = 1)) +
  labs(x = "Fun Rating", y = "Action")

#plot continuous vitalism score and energy rating, by type
ggplot(kid_expt1_main_fullvital, aes(x = vital_total, y = energy, color=vital_median_split)) +
  facet_wrap(~act_type)+
  geom_point(alpha=0.5)+
  geom_smooth(method = "lm", se = TRUE)

`geom_smooth()` using formula = 'y ~ x'

#individual action comparison

ggplot(subset(kid_expt1_main_fullvital, act_type=="mental"), 
       aes(x = energy, y = activity)) +
  geom_boxplot(outlier.shape = NA, width = 0.6) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, 
               fill = "black", color = "black",
               aes(group = act_type)) +  # diamond for mean
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", 
               width = 0.2, color = "black", linewidth = 0.8,
               aes(group = act_type)) +  # 95% CI
  facet_wrap(~vital_median_split) +
  theme_classic() +
  scale_x_continuous(
    breaks = c(-50, -25, 0, 25, 50),
   labels = c("-50\nExtremely tired", 
             "-25\nSomewhat tired", 
             "0\nNeutral", 
             "25\nSomewhat energetic", 
             "50\nExtremely energetic"))+
    theme_minimal()+
   theme(axis.text.x = element_text(angle = 15, hjust = 1)) +
  labs(x = "Energy Rating", y = "Action")

correlation matrix

# Correlation heatmaps by activity
corr_phys_lowvital <- kid_expt1_main_fullvital %>%
    filter(act_type == "physical" & vital_median_split == "<6") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

corr_phys_highvital <- kid_expt1_main_fullvital %>%
    filter(act_type == "physical" & vital_median_split == ">=6") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

corr_ment_lowvital <- kid_expt1_main_fullvital %>%
    filter(act_type == "mental"& vital_median_split == "<6") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

corr_ment_highvital <- kid_expt1_main_fullvital %>%
    filter(act_type == "mental"& vital_median_split == ">=6") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

corr_rest_lowvital <- kid_expt1_main_fullvital %>%
    filter(act_type == "restorative" & vital_median_split == "<6") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

corr_rest_highvital <- kid_expt1_main_fullvital %>%
    filter(act_type == "restorative" & vital_median_split == ">=6") %>%
    select(energy, difficulty, fun) %>%
    cor(method = "spearman",use = "pairwise.complete.obs")

fig_phys_low <- ggcorrplot(corr_phys_lowvital, lab = TRUE, lab_col = "black", title = "Physical Activities - low vital")
fig_phys_high <- ggcorrplot(corr_phys_highvital, lab = TRUE, lab_col = "black", title = "Physical Activities - high vital")
fig_ment_low <- ggcorrplot(corr_ment_lowvital, lab = TRUE, lab_col = "black", title = "Mental Activities - low vital")
fig_ment_high <- ggcorrplot(corr_ment_highvital, lab = TRUE, lab_col = "black", title = "Mental Activities - high vital")
fig_rest_low <- ggcorrplot(corr_rest_lowvital, lab = TRUE, lab_col = "black", title = "Restorative Activities - low vital")
fig_rest_high <- ggcorrplot(corr_rest_highvital, lab = TRUE, lab_col = "black", title = "Restorative Activities - high vital")

fig_phys_low

fig_phys_high

fig_ment_low

fig_ment_high

fig_rest_low

fig_rest_high

We don’t see that much of a difference in two vital groups in how energy, difficulty and fun correlate. For physical actions, high vital group showed less positive correlation for fun and energy (more fun less depleting), and more negative cor for difficulty and energy (more difficult more depleting). For mental actions only high vital group showed neg cor for fun and difficulty. For restorative actions only high vital group showed neg cor for difficulty and energy (the more difficult the less restorative)

Modeling

#vitalism as continuous measure
m1.1 <- lmer(energy ~ act_type + difficulty + fun + Age + (1|subjectID), data=kid_expt1_main_fullvital)

m1.2 <- lmer(energy ~ act_type + difficulty + fun + vital_total + (1|subjectID), data=kid_expt1_main_fullvital)

m1.5 <- lmer(energy ~ act_type + difficulty + fun + Age + vital_total + (1|subjectID), data=kid_expt1_main_fullvital)

m1.3 <- lmer(energy ~ act_type*Age + difficulty + fun + vital_total + (1|subjectID), data=kid_expt1_main_fullvital)

m1.4 <- lmer(energy ~ act_type*vital_median_split + difficulty + fun + Age + (1|subjectID), data=kid_expt1_main_fullvital)

BIC(m1.1,m1.2,m1.3,m1.4, m1.5) #m1.2 has lowest BIC

     df      BIC
m1.1  8 6750.531
m1.2  8 6746.765
m1.3 11 6754.827
m1.4 11 6751.124
m1.5  9 6750.143

summary(m1.4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ act_type * vital_median_split + difficulty + fun + Age +  
    (1 | subjectID)
   Data: kid_expt1_main_fullvital

REML criterion at convergence: 6679.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.65001 -0.70064 -0.05266  0.72709  2.17141 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept) 151.0    12.29   
 Residual              974.1    31.21   
Number of obs: 684, groups:  subjectID, 57

Fixed effects:
                                           Estimate Std. Error        df
(Intercept)                                 0.05861   12.19921  65.34472
act_typemental                             13.50825    4.41699 620.67603
act_typerestorative                        12.03482    4.48802 623.81978
vital_median_split>=6                      -6.49068    5.72018 125.81593
difficulty                                 -0.13122    0.03664 673.40468
fun                                         0.11901    0.03478 671.41156
Age                                        -0.68243    2.03785  53.34951
act_typemental:vital_median_split>=6       -3.87834    5.90300 620.93946
act_typerestorative:vital_median_split>=6   3.39208    5.91619 621.35629
                                          t value Pr(>|t|)    
(Intercept)                                 0.005 0.996181    
act_typemental                              3.058 0.002322 ** 
act_typerestorative                         2.682 0.007522 ** 
vital_median_split>=6                      -1.135 0.258659    
difficulty                                 -3.581 0.000367 ***
fun                                         3.421 0.000661 ***
Age                                        -0.335 0.739030    
act_typemental:vital_median_split>=6       -0.657 0.511417    
act_typerestorative:vital_median_split>=6   0.573 0.566611    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
               (Intr) act_typm act_typr v__>=6 dffclt fun    Age   
act_typmntl    -0.190                                              
act_typrstr    -0.187  0.494                                       
vtl_mdn_>=6     0.112  0.385    0.369                              
difficulty     -0.166  0.032    0.154   -0.042                     
fun            -0.187  0.024   -0.077    0.041  0.114              
Age            -0.917  0.000   -0.002   -0.369 -0.004  0.014       
act_typm:__>=6  0.125 -0.745   -0.357   -0.517  0.064  0.002  0.000
act_typr:__>=6  0.123 -0.370   -0.718   -0.517  0.091 -0.007 -0.001
               act_typm:__>=6
act_typmntl                  
act_typrstr                  
vtl_mdn_>=6                  
difficulty                   
fun                          
Age                          
act_typm:__>=6               
act_typr:__>=6  0.503        

plot(allEffects(m1.4))

car::vif(m1.2) #no multicollinearity issue

                GVIF Df GVIF^(1/(2*Df))
act_type    1.160678  2        1.037954
difficulty  1.145675  1        1.070362
fun         1.063073  1        1.031054
vital_total 1.005592  1        1.002792

The higher the vitalism score, the lower the energy rating, and no difference across action type. Age does not explain energy rating at all. That is, kids who are more vitalistic tend to rate physical, mental, restorative actions as more depleting.

#Lets break down by action type (depleting vs. restorative)

mental_phys_fullvital <- kid_expt1_main_fullvital |> filter(act_type!="restorative")
restorative_fullvital <- kid_expt1_main_fullvital |> filter(act_type=="restorative")

#depleting
m3.1 <- lmer(energy ~ act_type + difficulty + fun + vital_total + Age + (1|subjectID), data=mental_phys_fullvital)

summary(m3.1) 

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ act_type + difficulty + fun + vital_total + Age + (1 |  
    subjectID)
   Data: mental_phys_fullvital

REML criterion at convergence: 4434.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.51240 -0.69299 -0.08143  0.69472  2.47754 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept) 177.5    13.32   
 Residual              884.5    29.74   
Number of obs: 456, groups:  subjectID, 57

Fixed effects:
                Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)     -0.07763   13.15314  63.24373  -0.006 0.995310    
act_typemental  11.28826    2.81712 396.49487   4.007 7.35e-05 ***
difficulty      -0.14164    0.04230 441.71530  -3.348 0.000883 ***
fun              0.14452    0.04062 440.02973   3.558 0.000415 ***
vital_total     -1.43913    0.63754  53.59516  -2.257 0.028086 *  
Age              0.12212    2.23762  53.06740   0.055 0.956682    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) act_ty dffclt fun    vtl_tt
act_typmntl -0.140                            
difficulty  -0.201  0.146                     
fun         -0.223  0.054  0.162              
vital_total  0.135 -0.001 -0.025  0.067       
Age         -0.903  0.002  0.007  0.019 -0.434

plot(allEffects(m3.1))

#restorative
m3.2 <- lmer(energy ~ difficulty + fun + vital_total + (1|subjectID), data=restorative_fullvital)

summary(m3.2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ difficulty + fun + vital_total + (1 | subjectID)
   Data: restorative_fullvital

REML criterion at convergence: 2257.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1353 -0.7151  0.2140  0.7617  1.6208 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept)  143.7   11.99   
 Residual              1055.8   32.49   
Number of obs: 228, groups:  subjectID, 57

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)  
(Intercept)  16.44473    7.11302 119.07490   2.312   0.0225 *
difficulty   -0.07483    0.06886 223.54685  -1.087   0.2783  
fun           0.05666    0.06526 220.54278   0.868   0.3862  
vital_total  -1.22157    0.68655  55.67267  -1.779   0.0807 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) dffclt fun   
difficulty  -0.324              
fun         -0.622 -0.004       
vital_total -0.674  0.119  0.058

plot(allEffects(m3.2))

For depleting actions (mental and physical), kids with higher vital understanding tend to rate them as more depleting. For restorative actions this effect is marginally present (higher vital understanding leads to lower restoration rating)

#vitalism as categorical level (median split)
 
m2.1 <- lmer(energy ~ act_type + vital_median_split + difficulty + fun + (1|subjectID), data=kid_expt1_main_fullvital)

summary(m2.1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: energy ~ act_type + vital_median_split + difficulty + fun + (1 |  
    subjectID)
   Data: kid_expt1_main_fullvital

REML criterion at convergence: 6694.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.71661 -0.72103 -0.04387  0.73006  2.21580 

Random effects:
 Groups    Name        Variance Std.Dev.
 subjectID (Intercept) 147.2    12.13   
 Residual              973.4    31.20   
Number of obs: 684, groups:  subjectID, 57

Fixed effects:
                       Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)            -3.57459    4.52262 223.57933  -0.790 0.430143    
act_typemental         11.32210    2.94411 623.69177   3.846 0.000133 ***
act_typerestorative    13.89407    3.12220 634.34769   4.450 1.01e-05 ***
vital_median_split>=6  -7.35608    4.04150  54.72447  -1.820 0.074209 .  
difficulty             -0.13236    0.03645 676.23624  -3.631 0.000303 ***
fun                     0.11989    0.03475 674.42860   3.450 0.000595 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) act_typm act_typr v__>=6 dffclt
act_typmntl -0.390                                
act_typrstr -0.394  0.498                         
vtl_mdn_>=6 -0.533  0.003   -0.005                
difficulty  -0.494  0.119    0.316    0.014       
fun         -0.467  0.039   -0.118    0.064  0.116

plot(allEffects(m2.1))

m2.2 <- lmer(energy ~ act_type*vital_median_split + difficulty + fun + (1|subjectID), data=kid_expt1_main_fullvital)

m2.3 <- lmer(energy ~ act_type + vital_total + difficulty + fun + (1|subjectID), data=kid_expt1_main_fullvital)

BIC(m2.1,m2.2,m2.3)

     df      BIC
m2.1  8 6746.921
m2.2 10 6747.962
m2.3  8 6746.765

plot(allEffects(m2.2))

emmeans(m2.2, ~ act_type*vital_median_split)

 act_type    vital_median_split emmean   SE  df lower.CL upper.CL
 physical    <6                 -2.796 3.97 155   -10.63     5.04
 mental      <6                 10.713 3.96 154     2.89    18.53
 restorative <6                  9.231 4.00 159     1.33    17.13
 physical    >=6                -9.991 3.55 161   -17.00    -2.98
 mental      >=6                -0.361 3.51 156    -7.30     6.57
 restorative >=6                 5.425 3.57 164    -1.62    12.47

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

By median split, <6 group rated mental and restorative actions as energy-increasing, and physical actions as not changing energy; >=6 group rated physical actions as energy-decreasing, and mental and restorative actions as not changing energy.

#Variability

indiv_var <- kid_expt1_main_fullvital %>%
  group_by(subjectID, act_type, vital_median_split, vital_total, Age) %>%
  summarise(
    sd_energy = sd(energy, na.rm = TRUE),
    sd_difficulty = sd(difficulty, na.rm = TRUE),
    sd_fun = sd(fun, na.rm = TRUE),
    .groups = "drop"
  )

model1 <- lm(sd_energy ~ vital_median_split + act_type, data = indiv_var)
model2 <- lm(sd_energy ~ vital_total + act_type, data = indiv_var)

BIC(model1, model2)

       df      BIC
model1  5 1431.350
model2  5 1424.361

summary(model2)

Call:
lm(formula = sd_energy ~ vital_total + act_type, data = indiv_var)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.089  -9.711   0.502  10.754  26.001 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          32.0578     2.6420  12.134  < 2e-16 ***
vital_total          -0.8489     0.2845  -2.984  0.00328 ** 
act_typemental       -3.5301     2.7389  -1.289  0.19923    
act_typerestorative   2.0312     2.7389   0.742  0.45936    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.62 on 167 degrees of freedom
Multiple R-squared:  0.07287,   Adjusted R-squared:  0.05621 
F-statistic: 4.375 on 3 and 167 DF,  p-value: 0.005402

next steps:

• variability

• adult vs. kids

for adult: focus on energy ratings (use fun and difficulty as covariates, controlled for in model)

for kids - look at relationships between three measures and also compare adult