O1B Report - behavioral

Initial Analysis

1 Define functions

1.1 Quick ggwithinstats function

my_ggwithinstats <- function(data, title, x, y, outlier.label, xlab, ylab) {
  x <- rlang::enquo(x)
  y <- rlang::enquo(y)
  outlier.label <- rlang::enquo(outlier.label)
  
  data %>%
    ggstatsplot::ggwithinstats(
      x = !!x,
      y = !!y,
      title = title,
      xlab = xlab,
      ylab = ylab,
      outlier.tagging = TRUE,            # whether outliers need to be tagged
      outlier.label = !!outlier.label,   # variable to be used for tagging outliers
      outlier.coef = 2,
      pairwise.comparisons = TRUE,
      pairwise.display = "significant",
      results.subtitle = TRUE,
      type = "parametric",
      bf.message = FALSE, 
      p.adjust.method = "none",
      point.path = TRUE,
      ggtheme = ggprism::theme_prism(),
      # package = "RColorBrewer",  # "ggsci",
      # palette = "Dark",         # "default_jco",
      violin.args = list(width = 0.9, alpha = 0.2, size = 1, color = "black"),
      centrality.point.args = list(size = 5, color = "darkred"),
      centrality.label.args = list(size = 3, nudge_x = 0.2, segment.linetype = 5, fill = "#FFF8E7"), 
      ggplot.component = list(
        theme(
          plot.title = element_text(hjust = 0, size = 16),
          plot.subtitle = element_text(hjust = 0, size = 12), 
          plot.caption = element_text(hjust = 0, size = 12), 
          text = element_text(size = 14)
      ))
    ) + scale_colour_grey(start = 0.2, end = 0.2)   # hacky way to change point color
}

# For publication
my_ggwithinstats2 <- function(data, title, x, y, outlier.label, xlab, ylab, 
                              outlier.tagging = FALSE, results.subtitle = TRUE, 
                              centrality.label.args = TRUE, point.path = TRUE,
                              power = TRUE,
                              ...) {  # ... for limits and breaks
  x <- rlang::enquo(x)
  y <- rlang::enquo(y)
  outlier.label <- rlang::enquo(outlier.label)
  
  if(centrality.label.args){
    centrality.label.args <- list(size = 3, nudge_x = 0.2, segment.linetype = 5, fill = "#FFF8E7")
  }else{
    centrality.label.args <- list(size = 0, nudge_x = 10, segment.linetype = 0, alpha = 0) # very hacky way of not showing label
  }
  
  plot <- 
    data %>%
      ggstatsplot::ggwithinstats(
        x = !!x,
        y = !!y,
        title = title,
        xlab = xlab,
        ylab = ylab,
        outlier.tagging = outlier.tagging,                    # whether outlines need to be tagged
        outlier.label = !!outlier.label,                      # variable to be used for tagging outliers
        outlier.coef = 2,
        pairwise.comparisons = TRUE,
        pairwise.display = "all",
        results.subtitle = results.subtitle,
        type = "parametric",
        bf.message = FALSE, 
        p.adjust.method = "none",
        point.path = point.path,
        ggtheme = ggprism::theme_prism(palette = "black_and_white"),
        # package = "RColorBrewer",  # "ggsci",
        # palette = "Dark",         # "default_jco",
        violin.args = list(width = 0.9, alpha = 0.2, size = 1, color = "black"),
        centrality.plotting = TRUE,
        centrality.type = "parameteric",
        centrality.point.args = list(size = 5, color = "black"),
        centrality.path.args = list(color = "black", size = 1, alpha = 1),
        centrality.label.args = centrality.label.args,
        ggplot.component = list(
          theme(
            plot.title = element_text(hjust = 0, size = 16),
            plot.subtitle = element_text(hjust = 0, size = 12), 
            plot.caption = element_text(hjust = 0, size = 12), 
            text = element_text(family = "Sans", size = 14)
        ))
      ) + scale_colour_grey(start = 0.2, end = 0.2) +  # hacky way to change point color
      scale_y_continuous(...)
  
  if(power) {
  stats_df <- ggstatsplot::extract_stats(plot)$subtitle_data
  ttest_power <- WebPower::wp.t(
    n1 = stats_df$df.error + 1,         # it's df + 1 
    d = stats_df$estimate,              # this is Hedges'g, very close to d, but less biased
    type = "paired"
  )
  ttest_power %>% print()
  }  
  
  plot
  
}

# Fast ggsave - saves plot with filename of R plot object
fast_ggsave <- function(plot, device = "tiff", path = NULL,
                        units = "in", dpi = 300, width = 5, height = 5, ...){ 
  plot_name <- deparse(substitute(plot))
  ggplot2::ggsave(filename = paste0(plot_name, ".", device), plot = plot,
                  device = device, path = path,
                  units = units, dpi = dpi,
                  width = width, height = height,
                  ...
  )
  
} # use: fast_ggsave(jrad_ox_p, path = savefolder)

# Fast tiff save
fast_tiffsave <- function(plot, path = NULL,
                          units = "in", res = 300, width = 5, height = 5, ...){ 
  plot_name <- deparse(substitute(plot))
  tiff(filename = file.path(path, paste0(plot_name, ".", "tiff")),
       units = units, res = res,
       width = width, height = height,
       ...
  )
  plot(plot)
  dev.off()
}  # use: fast_tiffsave(jrad_ox_p, path = savefolder)

# Errors with ggplot2  --- can use this to save:
# Cairo::Cairo(
#   width = 5, 
#   height = 5, 
#   file = file.path(savefolder, paste0("jrad_ox_p", ".", "tiff")),
#   type = "png", 
#   bg =  "white",    # "transparent" 
#   dpi = 300,
#   units = "in"
# )
# plot(jrad_ox_p) 
# dev.off()

1.2 rmANOVA

1.2.1 helper for mixed type to long format

select_to_longer <- function(data, var_prefix = NULL) {
  var_regex <- paste0("(", var_prefix, ")_(pre|post)")
  
  data %>%
    dplyr::select(ID, Conditie, dplyr::matches(var_regex)) %>%
    tidyr::pivot_longer(cols = dplyr::matches(var_regex), 
                 names_to = c("Variable", "PrePost"), 
                 names_pattern = "(.*)_(pre|post)",   
                 values_to = var_prefix
    ) %>% 
    dplyr::mutate(
      ID = as.factor(ID),
      Cond = factor(Conditie, levels = c("ES", "JRAD")),
      PrePost = factor(PrePost, levels = c("pre", "post"))
    ) %>%
    dplyr::select(-c(Conditie, Variable))
}  # select_to_longer(long_df, "IRI_PD")

1.2.2 rmANOVA

# Define Function for Two-way rmANOVA
# Check with https://www.datanovia.com/en/lessons/repeated-measures-anova-in-r/

tw_rmANOVA_func <- 
  function(data, id_var, cond_var, time_var, value_var, 
           assum_check = TRUE, p_adjust_method = "bonferroni") {
  
  # input dataframe needs to have columns names different from "variable" and "value" because it collides with rstatix::shapiro_test
  # github issue was raised: https://github.com/kassambara/rstatix/issues/52
  
  # Dependencies    
  library(dplyr)
  library(rlang)
  library(ggpubr)
  library(ggplot2)
  library(rstatix)
  library(WebPower)
  library(afex)  
  library(effectsize)
  
  # Defaults    
  posthoc_sig_interac <- FALSE
  posthoc_ns_interac <- FALSE  
      
  id_var_enq <- rlang::enquo(id_var)
  id_var_name <- rlang::as_name(id_var_enq)
  cond_var_enq <- rlang::enquo(cond_var)
  cond_var_name <- rlang::as_name(cond_var_enq)
  time_var_enq <- rlang::enquo(time_var)
  time_var_name <- rlang::as_name(time_var_enq)
  value_var_enq <- rlang::enquo(value_var)
  value_var_name <- rlang::as_name(value_var_enq)
  
  anova_formula_str <- paste0(value_var_name, " ~ ",
                        cond_var_name, " * ", time_var_name, " + ",
                        "Error(", id_var_name, " / ",
                        "(", cond_var_name, " * ", time_var_name, "))") 
  anova_formula <- as.formula(anova_formula_str)
    
  # Assumptions
  if(assum_check){
    cat("\n Outliers \n")
    data %>%
      dplyr::group_by(!!cond_var_enq, !!time_var_enq) %>%
      rstatix::identify_outliers(!!value_var_enq) %>%                                  # outliers (needs to be 0)
      print()
    
    cat("\n Normality assumption (p >.05) \n")
    data %>%
      dplyr::group_by(!!cond_var_enq, !!time_var_enq) %>%
      rstatix::shapiro_test(!!value_var_enq) %>%                                        # normality assumption (p>.05)
      print()
    
    qq_plot <- 
      ggpubr::ggqqplot(data = data, value_var_name, ggtheme = theme_bw(), title = "QQ Plot") +
        ggplot2::facet_grid(vars(!!time_var_enq), vars(!!cond_var_enq), labeller = "label_both")    # QQ plot
    
    qq_plot %>% print()
  }
  
  
  # Two-way rmANOVA - check for interaction (ex. F(2, 22) = 30.4, p < 0.0001)
  cat("\n Two-way rmANOVA \n")
  res_aov <- rstatix::anova_test(                             # automatically does sphericity Mauchly’s test
    data = data, dv = !!value_var_enq, wid = !!id_var_enq,    # also could have used anova_formula                         
    within = c(!!cond_var_enq, !!time_var_enq),
    type = 2, effect.size = "ges", detailed = TRUE
  )
  
  res_aov %>% print()
  # DFn = Degrees of Freedom in the numerator (a.k.a. DFeffect)
  # DFd = Degrees of Freedom in the denominator (a.k.a. DFerror)
  # SSn = Sum of Squares in the numerator (a.k.a. SSeffect)
  # SSd = Sum of Squares in the denominator (a.k.a. SSerror)
  
    # "auto": apply automatically GG correction to only within-subjects factors violating the sphericity assumption (i.e., Mauchly's test p-value is significant, p <= 0.05).
  res_aov_table <- get_anova_table(res_aov, correction = c("auto"))
  
  # Get Cohen's f partial for power analysis
  cat("\n Cohen's f (partial) \n")
  cat("Used formula: ", anova_formula_str)
  
  f_cohen <- 
    data %>%
    afex::aov_car(anova_formula, data = ., type = 2) %>%    # aov(anova_formula, data = .) is the same but with type = 1 
    effectsize::cohens_f(partial = TRUE) 
  f_cohen %>% print()
  
  f_cohen_cond <- f_cohen$Cohens_f_partial[[1]]
  f_cohen_time <- f_cohen$Cohens_f_partial[[2]]
  f_cohen_interac <- f_cohen$Cohens_f_partial[[3]]
  
  # Sphericity tests and Sphericity Correction (e.g. GG) make sens only for factors with more than 2 levels
  if(length(unique(data[, cond_var_name])) > 2 | length(unique(data[, time_var_name])) > 2) {
    # res_aov$`Mauchly's Test for Sphericity`
    # res_aov$`Sphericity Corrections`
    gge_time <- res_aov$`Sphericity Corrections`$GGe[[1]]   # time
    gge_interac <- res_aov$`Sphericity Corrections`$GGe[[2]]   # interaction
  }else {
    gge_time <- 1
    gge_interac <- 1
  } 
  

  # default p-value sig level for interaction term is 0.05
  if(res_aov_table[, "p"][[4]] <= 0.05) {
    cat("\n Following reporting procedure for significant two-way interaction \n")
    posthoc_sig_interac <- TRUE
  } else {
    cat("\n Following reporting procedure for non-significant two-way interaction \n")
    posthoc_ns_interac <- TRUE
  }
  
  # Pairwise comparisons between treatment groups
  # Compute it for plot, regardless of interaction significance, print it only for significant interaction
  pwc <-
    data %>%
    dplyr::group_by(!!time_var_enq) %>%
    rstatix::pairwise_t_test(
      as.formula(paste0(value_var_name, " ~ ", cond_var_name)),
      paired = TRUE,
      p.adjust.method = p_adjust_method
    )
  
  
  #- Procedure for a significant two-way interaction -
  if(posthoc_sig_interac) {
    
    cat("\n Effect of treatment at each time point \n")
    one_way <- 
      data %>%
      dplyr::group_by(!!time_var_enq) %>%
      rstatix::anova_test(dv = !!value_var_enq, wid = !!id_var_enq, within = !!cond_var_enq) %>%
      rstatix::get_anova_table() %>%
      rstatix::adjust_pvalue(method = "bonferroni")
    one_way %>% print()
    
    cat("\n Pairwise comparisons between treatment groups \n")
    pwc %>% print()

    cat("\n Effect of time at each level of treatment - One-way ANOVA \n")
    one_way2 <-
      data %>%
      dplyr::group_by(!!cond_var_enq) %>%
      rstatix::anova_test(dv = !!value_var_enq, wid = !!id_var_enq, within = !!time_var_enq) %>%
      rstatix::get_anova_table() %>%
      adjust_pvalue(method = p_adjust_method)
    one_way2  %>% print()
    
    cat("\n Pairwise comparisons between time points \n")
    pwc2 <-
      data %>%
      dplyr::group_by(!!cond_var_enq) %>%
      rstatix::pairwise_t_test(
        as.formula(paste0(value_var_name, " ~ ", time_var_name)),     # paste formula, not quosure
        paired = TRUE,
        p.adjust.method = p_adjust_method
      )
    pwc2  %>% print()
    
  }
  
  #- Procedure for non-significant two-way interaction- 
  # If the interaction is not significant, you need to interpret the main effects for each of the two variables: treatment and time.
  if(posthoc_ns_interac){
    
    cat("\n Comparisons for treatment variable \n")
    pwc_cond <-
      data %>%
      rstatix::pairwise_t_test(
        as.formula(paste0(value_var_name, " ~ ", cond_var_name)),     # paste formula, not quosure             
        paired = TRUE,
        p.adjust.method = p_adjust_method
      )
    pwc_cond %>% print()
    
    cat("\n Comparisons for time variable \n")
    pwc_time <-
      data %>% 
      rstatix::pairwise_t_test(
        as.formula(paste0(value_var_name, " ~ ", time_var_name)),     # paste formula, not quosure
        paired = TRUE,
        p.adjust.method = p_adjust_method
      )
    pwc_time %>% print()
    
  }
  
  # Post-hoc Power Analysis
  cat("\n Post-hoc Power Analysis for Repeated-measures ANOVA \n")
  pwr_rmanova_cond <-
    WebPower::wp.rmanova(n = length(unique(data[[id_var_name]])),
                         ng = length(unique(data[[cond_var_name]])),
                         nm = length(unique(data[[time_var_name]])),
                         f = f_cohen_cond,
                         nscor = gge_time,          # don't have gge for cond, this should be close
                         type = 1                   # for within-effect; 1 is for between-effect
    )
  pwr_rmanova_cond$note %>% print()
  pwr_rmanova_cond[1:7] %>%
    as.data.frame() %>% print()

  pwr_rmanova_time <-
    WebPower::wp.rmanova(n = length(unique(data[[id_var_name]])),
                         ng = length(unique(data[[cond_var_name]])),
                         nm = length(unique(data[[time_var_name]])),
                         f = f_cohen_time,
                         nscor = gge_time,
                         type = 1                  # for within-effect
    )
  pwr_rmanova_time$note %>% print()
  pwr_rmanova_time[1:7] %>%
    as.data.frame() %>% print()

  pwr_rmanova_interac <-
    WebPower::wp.rmanova(n = length(unique(data[[id_var_name]])),
                         ng = length(unique(data[[cond_var_name]])),
                         nm = length(unique(data[[time_var_name]])),
                         f = f_cohen_interac,
                         nscor = gge_interac,
                         type = 2                  # for interaction
    )
  pwr_rmanova_interac$note %>% print()
  pwr_rmanova_interac[1:7] %>%
    as.data.frame() %>% print()

  
  # Visualization
  violin_plot <- 
    ggpubr::ggviolin(data, x = time_var_name, y = value_var_name, color = cond_var_name, 
                     add = c("boxplot", "mean_se", "dotplot"), 
                     add.params = list(group = cond_var_name, fill = "white")  
                     ) +
    ggplot2::scale_color_grey(start = 0.7, end = 0.3) + 
    ggplot2::theme_classic()
  
  pwc <- 
    pwc %>% 
    rstatix::add_xy_position(x = time_var_name)
  
  violin_plot <- 
    violin_plot + 
    ggpubr::stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
    labs(
      subtitle = get_test_label(res_aov, detailed = TRUE),
      caption = get_pwc_label(pwc)
    )
  
    violin_plot %>% print()
  
}

# ex. - run on long format
# selfesteem2 <- datarium::selfesteem2
# selfesteem2 <- selfesteem2 %>%                         # interac, within, between
#   tidyr::gather(key = "time", value = "score", t1, t2, t3) %>%
#   rstatix::convert_as_factor(id, time)
# selfesteem3 <- selfesteem2 %>%         # sig interact, sig within
#   dplyr::filter(time %in% c("t1", "t2")) %>% 
#   dplyr::mutate(time = factor(time, levels = c("t1", "t2"))) 
# selfesteem4 <- selfesteem2 %>%        # sig interact, sig within, sig between
#   dplyr::filter(time %in% c("t2", "t3")) %>% 
#   dplyr::mutate(time = factor(time, levels = c("t2", "t3"))) 
# tw_rmANOVA_func(data = selfesteem2, id_var = id, cond_var = treatment, time_var = time, value_var = score)
# 
# weightloss <- datarium::weightloss
# weightloss2 <- weightloss %>%
#   dplyr::mutate(treatment = paste0(diet, "_diet", "/", exercises, "_exerc")) %>%
#   dplyr::filter(treatment %in% c("yes_diet/no_exerc", "no_diet/yes_exerc")) %>%
#   tidyr::gather(key = "time", value = "score", t1, t2, t3) %>%
#   rstatix::convert_as_factor(id, time) 
# weightloss3 <- weightloss2 %>%
#   dplyr::filter(time %in% c("t1", "t2")) %>% 
#   dplyr::mutate(time = factor(time, levels = c("t1", "t2")))
# weightloss4 <- weightloss %>%
#   dplyr::mutate(treatment = paste0(diet, "_diet", "/", exercises, "_exerc")) %>%
#   dplyr::filter(treatment %in% c("no_diet/no_exerc", "yes_diet/no_exerc")) %>%
#   tidyr::gather(key = "time", value = "score", t1, t2, t3) %>%
#   rstatix::convert_as_factor(id, time) 
# weightloss5 <- weightloss4 %>%
#   dplyr::filter(time %in% c("t1", "t2")) %>% 
#   dplyr::mutate(time = factor(time, levels = c("t1", "t2")))
# weightloss6 <- weightloss4 %>%
#   dplyr::filter(time %in% c("t2", "t3")) %>% 
#   dplyr::mutate(time = factor(time, levels = c("t2", "t3")))
# tw_rmANOVA_func(data = weightloss3, id_var = id, cond_var = treatment, time_var = time, value_var = score) # within & interac
# tw_rmANOVA_func(data = weightloss4, id_var = id, cond_var = treatment, time_var = time, value_var = score) # between & within
# tw_rmANOVA_func(data = weightloss5, id_var = id, cond_var = treatment, time_var = time, value_var = score) # between 

2 Read data

folder <- "C:/Users/Mihai/Desktop/R Notebooks/notebooks/o1b-report-behavior"
file <- "O1b data.RDS"

savefolder <- "C:/Users/Mihai/Desktop/R Notebooks/notebooks/o1b-report-behavior/Art"

# Data with all Participants
long_df_all <- readRDS(file.path(folder, file))


# Correct a typo on participant 11
long_df_all[(long_df_all$ID == 11 & long_df_all$Conditie == "JRAD"), "APS_post"] <- 28
long_df_all <- 
  long_df_all %>% 
  dplyr::mutate(Diff_APS = APS_post - APS_pre)   # need to redo this because of correction

# Remove a participant with odd data and no OXT data
long_df_all <-
  long_df_all %>%
  dplyr::filter(ID != 23)

# Data only with participants with both Cond 
ids_single_Cond <- 
  long_df_all %>%
  dplyr::count(ID) %>%
  dplyr::filter(n < 2) %>%
  dplyr::pull(ID)

long_df <- 
  long_df_all %>%
  dplyr::filter(!ID %in% ids_single_Cond)

## Probably leave NAs in
# long_df_nona <- 
#   long_df %>%
#   group_by(ID) %>%
#   tidyr::drop_na(PA_pre:APS_post, IOS:OX_post) %>%
#   ungroup()

jrad_df <-
  long_df_all %>%
  dplyr::filter(Conditie == "JRAD")

es_df <-
  long_df_all %>%
  dplyr::filter(Conditie == "ES")

3 Demografics

length(unique(long_df_all$ID))  # 41 + id 21 twice in JRAD

[1] 41

gender_stat <- 
  long_df_all %>%
    pivot_wider(id_cols = ID, values_from = Gen, names_from = Conditie) %>%    # here is a problem because participant 21 is in JRAD twice
    unnest(ES, JRAD) %>%
    dplyr::mutate(Gen = coalesce(JRAD, ES))

gender_stat %>%
  dplyr::group_by(ID) %>%
  dplyr::summarise(n = n()) %>%          # 41 particpants but id 21 appears twice in JRAD
  dplyr::summarise(n_subj = sum(n))

gender_stat %>%
  dplyr::group_by(Gen) %>%
  dplyr::summarise(n = n()) %>%
  dplyr::mutate(freq = n / sum(n))

 

age_stat <- 
  long_df_all %>%
    pivot_wider(id_cols = ID, values_from = Varsta, names_from = Conditie) %>%    # here is a problem because participant 21 is in JRAD twice
    unnest(ES, JRAD) %>%
    dplyr::mutate(Age = coalesce(JRAD, ES))
age_stat %>%
  rstatix::get_summary_stats(Age, type = "common")

jrad_df %>%
  dplyr::group_by(Gen) %>%
  dplyr::summarise(n = n()) %>%
  dplyr::mutate(freq = n / sum(n))

jrad_df %>%
  rstatix::get_summary_stats(Varsta, type = "common")

ggplot(data = jrad_df, aes(x = Varsta)) + geom_histogram(aes(fill = Gen))

t.test(jrad_df$Varsta ~ jrad_df$Gen) %>% broom::tidy()

es_df %>%
  dplyr::group_by(Gen) %>%
  dplyr::summarise(n = n()) %>%
  dplyr::mutate(freq = n / sum(n))

es_df %>%
  rstatix::get_summary_stats(Varsta, type = "common")

ggplot(data = es_df, aes(x = Varsta)) + geom_histogram(aes(fill = Gen))

t.test(es_df$Varsta ~ es_df$Gen) %>% broom::tidy()

4 Some new resutls

lm(Diff_PA ~ PA_pre + StaiS_pre, data = jrad_df) %>% summary()

Call:
lm(formula = Diff_PA ~ PA_pre + StaiS_pre, data = jrad_df)

Residuals:
   Min     1Q Median     3Q    Max 
-7.756 -3.325 -0.620  2.034 12.687 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 11.64786    5.30764   2.195   0.0351 *
PA_pre      -0.16476    0.12600  -1.308   0.1998  
StaiS_pre   -0.13068    0.07355  -1.777   0.0845 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.526 on 34 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.1025,    Adjusted R-squared:  0.04969 
F-statistic: 1.941 on 2 and 34 DF,  p-value: 0.1591

lm(Diff_PA ~ PA_pre + StaiS_pre, data = es_df) %>% summary()

Call:
lm(formula = Diff_PA ~ PA_pre + StaiS_pre, data = es_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4592 -1.8704 -0.8811  2.1809  8.9339 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 11.64405    5.75538   2.023   0.0512 .
PA_pre      -0.15939    0.13910  -1.146   0.2601  
StaiS_pre   -0.16008    0.06569  -2.437   0.0204 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.707 on 33 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.1533,    Adjusted R-squared:  0.1019 
F-statistic: 2.987 on 2 and 33 DF,  p-value: 0.06425

lm(Diff_PA ~ PA_pre + OX_post, data = jrad_df) %>% summary()

Call:
lm(formula = Diff_PA ~ PA_pre + OX_post, data = jrad_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.0453 -2.9057 -0.4404  1.6999 11.4906 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -13.87760    8.06331  -1.721   0.0946 .
PA_pre       -0.02761    0.12067  -0.229   0.8204  
OX_post      16.21583    6.43889   2.518   0.0168 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.318 on 33 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.1674,    Adjusted R-squared:  0.117 
F-statistic: 3.318 on 2 and 33 DF,  p-value: 0.04864

lm(Diff_PA ~ PA_pre + OX_post, data = es_df) %>% summary()

Call:
lm(formula = Diff_PA ~ PA_pre + OX_post, data = es_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.3405 -2.0382 -0.2352  2.2273  9.0547 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -7.9484     6.8726  -1.157    0.256
PA_pre        0.0730     0.1289   0.566    0.575
OX_post       7.0514     5.0267   1.403    0.170

Residual standard error: 3.625 on 32 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.06283,   Adjusted R-squared:  0.004253 
F-statistic: 1.073 on 2 and 32 DF,  p-value: 0.3541

# lm(Diff_APS ~ APS_pre + OX_post, data = jrad_df) %>% summary()
# lm(Diff_PA ~ PA_pre + OX_post, data = es_df) %>% summary()

# lm(Diff_NA ~ NA_pre + StaiS_pre, data = jrad_df) %>% summary()
# lm(Diff_NA ~ NA_pre + StaiS_pre, data = es_df) %>% summary()


jrad_df %>%
  dplyr::select(-c(ID, Conditie, contains("IRI"))) %>%
  correlation::correlation(p_adjust = "none") %>%
  dplyr::filter(p < .05) 

Registered S3 methods overwritten by 'parameters':
  method                           from      
  as.double.parameters_kurtosis    datawizard
  as.double.parameters_skewness    datawizard
  as.double.parameters_smoothness  datawizard
  as.numeric.parameters_kurtosis   datawizard
  as.numeric.parameters_skewness   datawizard
  as.numeric.parameters_smoothness datawizard
  print.parameters_distribution    datawizard
  print.parameters_kurtosis        datawizard
  print.parameters_skewness        datawizard
  summary.parameters_kurtosis      datawizard
  summary.parameters_skewness      datawizard

# Correlation Matrix (pearson-method)

Parameter1    |     Parameter2 |     r |         95% CI |     t | df |         p
--------------------------------------------------------------------------------
StaiS_pre     |     StaiS_post |  0.92 | [ 0.83,  0.96] | 11.98 | 27 | < .001***
StaiS_pre     |        PA_post | -0.39 | [-0.63, -0.08] | -2.51 | 35 | 0.017*   
StaiS_pre     |         NA_pre |  0.74 | [ 0.54,  0.85] |  6.52 | 36 | < .001***
StaiS_pre     |        NA_post |  0.68 | [ 0.46,  0.82] |  5.53 | 35 | < .001***
StaiS_pre     |  Rosenberg_pre | -0.53 | [-0.73, -0.25] | -3.75 | 36 | < .001***
StaiS_pre     | Rosenberg_post | -0.56 | [-0.75, -0.30] | -4.11 | 36 | < .001***
StaiS_post    |         PA_pre | -0.52 | [-0.74, -0.19] | -3.15 | 27 | 0.004**  
StaiS_post    |        PA_post | -0.65 | [-0.83, -0.37] | -4.40 | 26 | < .001***
StaiS_post    |         NA_pre |  0.67 | [ 0.41,  0.83] |  4.75 | 27 | < .001***
StaiS_post    |        NA_post |  0.68 | [ 0.41,  0.84] |  4.70 | 26 | < .001***
StaiS_post    |  Rosenberg_pre | -0.55 | [-0.76, -0.23] | -3.40 | 27 | 0.002**  
StaiS_post    | Rosenberg_post | -0.68 | [-0.84, -0.41] | -4.76 | 27 | < .001***
PA_pre        |        PA_post |  0.77 | [ 0.60,  0.88] |  7.19 | 35 | < .001***
PA_post       |         NA_pre | -0.33 | [-0.59, -0.01] | -2.07 | 35 | 0.046*   
PA_post       |            IOS |  0.36 | [ 0.04,  0.61] |  2.26 | 35 | 0.030*   
PA_post       |        Diff_PA |  0.52 | [ 0.24,  0.72] |  3.63 | 35 | < .001***
NA_pre        |        NA_post |  0.88 | [ 0.78,  0.94] | 11.00 | 35 | < .001***
NA_pre        |       APS_post | -0.33 | [-0.58, -0.01] | -2.06 | 36 | 0.046*   
NA_pre        |  Rosenberg_pre | -0.38 | [-0.62, -0.07] | -2.45 | 36 | 0.019*   
NA_pre        | Rosenberg_post | -0.42 | [-0.65, -0.12] | -2.77 | 36 | 0.009**  
NA_pre        |        Diff_PA | -0.42 | [-0.66, -0.12] | -2.77 | 35 | 0.009**  
NA_pre        |        Diff_NA | -0.49 | [-0.70, -0.20] | -3.36 | 35 | 0.002**  
NA_post       |  Rosenberg_pre | -0.40 | [-0.64, -0.08] | -2.56 | 35 | 0.015*   
NA_post       | Rosenberg_post | -0.47 | [-0.69, -0.17] | -3.17 | 35 | 0.003**  
NA_post       |        Diff_PA | -0.40 | [-0.64, -0.09] | -2.58 | 35 | 0.014*   
APS_pre       |       APS_post |  0.93 | [ 0.88,  0.97] | 15.65 | 36 | < .001***
APS_pre       |        OX_post |  0.43 | [ 0.13,  0.66] |  2.83 | 35 | 0.008**  
APS_pre       |        Diff_OX |  0.36 | [ 0.04,  0.61] |  2.27 | 35 | 0.029*   
APS_post      |        OX_post |  0.45 | [ 0.15,  0.68] |  2.98 | 35 | 0.005**  
APS_post      |        Diff_OX |  0.37 | [ 0.06,  0.62] |  2.38 | 35 | 0.023*   
APS_post      |        Diff_PA |  0.37 | [ 0.05,  0.62] |  2.35 | 35 | 0.025*   
Rosenberg_pre | Rosenberg_post |  0.88 | [ 0.77,  0.93] | 10.87 | 36 | < .001***
OX_pre        |        Diff_OX | -0.64 | [-0.80, -0.39] | -4.87 | 35 | < .001***
OX_post       |        Diff_OX |  0.66 | [ 0.43,  0.81] |  5.20 | 35 | < .001***
OX_post       |        Diff_PA |  0.41 | [ 0.09,  0.65] |  2.60 | 34 | 0.014*   
Diff_PA       |       Diff_APS |  0.33 | [ 0.01,  0.59] |  2.07 | 35 | 0.046*   

p-value adjustment method: none
Observations: 28-38

es_df %>%
  dplyr::select(-c(ID, Conditie, contains("IRI"))) %>%
  correlation::correlation(p_adjust = "none") %>%
  dplyr::filter(p < .05)

# Correlation Matrix (pearson-method)

Parameter1    |     Parameter2 |     r |         95% CI |     t | df |         p
--------------------------------------------------------------------------------
StaiS_pre     |     StaiS_post |  0.64 | [ 0.34,  0.82] |  4.09 | 24 | < .001***
StaiS_pre     |         PA_pre | -0.40 | [-0.64, -0.09] | -2.60 | 35 | 0.014*   
StaiS_pre     |        PA_post | -0.54 | [-0.74, -0.25] | -3.70 | 34 | < .001***
StaiS_pre     |         NA_pre |  0.58 | [ 0.32,  0.76] |  4.23 | 35 | < .001***
StaiS_pre     |        NA_post |  0.53 | [ 0.24,  0.73] |  3.60 | 34 | 0.001**  
StaiS_pre     |  Rosenberg_pre | -0.41 | [-0.65, -0.10] | -2.65 | 35 | 0.012*   
StaiS_pre     |        Diff_PA | -0.35 | [-0.61, -0.02] | -2.15 | 34 | 0.039*   
StaiS_post    |        PA_post | -0.46 | [-0.72, -0.09] | -2.56 | 24 | 0.017*   
StaiS_post    |         NA_pre |  0.45 | [ 0.07,  0.71] |  2.46 | 24 | 0.022*   
StaiS_post    |        NA_post |  0.77 | [ 0.54,  0.89] |  5.91 | 24 | < .001***
StaiS_post    |        APS_pre | -0.39 | [-0.68,  0.00] | -2.07 | 24 | 0.049*   
StaiS_post    |  Rosenberg_pre | -0.46 | [-0.72, -0.09] | -2.52 | 24 | 0.019*   
StaiS_post    | Rosenberg_post | -0.40 | [-0.68, -0.02] | -2.15 | 24 | 0.042*   
StaiS_post    |        Diff_PA | -0.50 | [-0.74, -0.14] | -2.84 | 24 | 0.009**  
PA_pre        |        PA_post |  0.78 | [ 0.60,  0.88] |  7.16 | 34 | < .001***
PA_pre        |  Rosenberg_pre |  0.47 | [ 0.18,  0.69] |  3.18 | 35 | 0.003**  
PA_pre        | Rosenberg_post |  0.39 | [ 0.07,  0.64] |  2.46 | 34 | 0.019*   
PA_post       |  Rosenberg_pre |  0.45 | [ 0.14,  0.68] |  2.93 | 34 | 0.006**  
PA_post       | Rosenberg_post |  0.44 | [ 0.13,  0.67] |  2.88 | 34 | 0.007**  
PA_post       |        Diff_PA |  0.61 | [ 0.35,  0.78] |  4.47 | 34 | < .001***
NA_pre        |        NA_post |  0.77 | [ 0.60,  0.88] |  7.13 | 34 | < .001***
NA_pre        |        Diff_NA | -0.63 | [-0.80, -0.38] | -4.75 | 34 | < .001***
APS_pre       |       APS_post |  0.87 | [ 0.76,  0.93] | 10.24 | 33 | < .001***
APS_post      |       Diff_APS |  0.62 | [ 0.36,  0.79] |  4.56 | 33 | < .001***
Rosenberg_pre | Rosenberg_post |  0.89 | [ 0.79,  0.94] | 11.26 | 34 | < .001***
OX_pre        |        Diff_OX | -0.70 | [-0.84, -0.49] | -5.74 | 34 | < .001***
OX_post       |        Diff_OX |  0.61 | [ 0.35,  0.78] |  4.43 | 34 | < .001***

p-value adjustment method: none
Observations: 26-37

PerformanceAnalytics::chart.Correlation(jrad_df[, c("APS_pre", "APS_post", "Diff_APS", "OX_pre", "OX_post", "Diff_OX")])

PerformanceAnalytics::chart.Correlation(jrad_df[, c("PA_pre", "PA_post", "Diff_PA", "OX_pre", "OX_post", "Diff_OX")])

PerformanceAnalytics::chart.Correlation(jrad_df[, c("NA_pre", "NA_post", "Diff_NA", "OX_pre", "OX_post", "Diff_OX")])

PerformanceAnalytics::chart.Correlation(jrad_df[, c("NA_pre", "NA_post", "Diff_NA", "APS_pre", "APS_post", "Diff_APS")])

PerformanceAnalytics::chart.Correlation(jrad_df[, c("StaiS_pre", "StaiS_post", "OX_pre", "OX_post", "Diff_OX")])

PerformanceAnalytics::chart.Correlation(jrad_df[, c("StaiS_pre", "StaiS_post", "APS_pre", "APS_post", "Diff_APS")])

PerformanceAnalytics::chart.Correlation(jrad_df[, c("StaiS_pre", "StaiS_post", "PA_pre", "PA_post", "Diff_PA")])

PerformanceAnalytics::chart.Correlation(jrad_df[, c("StaiS_pre", "StaiS_post", "NA_pre", "NA_post", "Diff_NA")])

library(MASS)

Attaching package: ‘MASS’

The following object is masked from ‘package:rstatix’:

    select

The following object is masked from ‘package:dplyr’:

    select

bla <- jrad_df %>%
  dplyr::select(-c(ID, Conditie, contains("IRI"), NA_pre, NA_post))

full_model <- lm(Diff_NA ~., data = bla)

step_model <- stepAIC(full_model, direction = "both", trace = FALSE)
summary(step_model)

Call:
lm(formula = Diff_NA ~ StaiS_pre + StaiS_post + PA_post + OX_post, 
    data = bla)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.9041 -0.9345  0.2900  1.2005  2.9167 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -10.98693    5.41805  -2.028   0.0549 .
StaiS_pre    -0.25666    0.09528  -2.694   0.0133 *
StaiS_post    0.27020    0.11288   2.394   0.0256 *
PA_post       0.11133    0.07187   1.549   0.1356  
OX_post       5.56450    3.65230   1.524   0.1419  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.028 on 22 degrees of freedom
  (11 observations deleted due to missingness)
Multiple R-squared:  0.4127,    Adjusted R-squared:  0.3059 
F-statistic: 3.865 on 4 and 22 DF,  p-value: 0.01588

5 OXT

jrad_df %>%
  rstatix::get_summary_stats(OX_pre, OX_post, type = "common")

# Difference in results for OXT is caused by subject excluded in left_join (they had OXT data but not in psychological data)
jrad_ox_p <- 
  jrad_df %>%
    dplyr::select(ID, OX_pre, OX_post) %>%
    tidyr::pivot_longer(cols = c(OX_pre, OX_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "RP",
      ylab = "OXT (pg/ml total protein)",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = .8, to = 1.3, by = .1),
      limits = c(.75, 1.32)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha      power
    37 0.0797434  0.05 0.07591611

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

jrad_ox_p

# fast_tiffsave(jrad_ox_p, path = savefolder)



es_df %>%
  rstatix::get_summary_stats(OX_pre, OX_post, type = "common")

es_ox_p <- 
  es_df %>%
    dplyr::select(ID, OX_pre, OX_post) %>%
    tidyr::pivot_longer(cols = c(OX_pre, OX_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "MS",
      ylab = "OXT (pg/ml total protein)",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = .8, to = 1.3, by = .1),
      limits = c(.75, 1.32)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n        d alpha     power
    36 0.363632  0.05 0.5642651

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

es_ox_p

# fast_tiffsave(es_ox_p, path = savefolder)

6 Behavioral

# STAI
jrad_df %>%
  rstatix::get_summary_stats(StaiS_pre, StaiS_post, type = "common")

jrad_stai_p <- 
  jrad_df %>%
    dplyr::select(ID, StaiS_pre, StaiS_post) %>%
    tidyr::pivot_longer(cols = c(StaiS_pre, StaiS_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "RP",
      ylab = "Anxiety",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 20, to = 65, by = 5),
      limits = c(19, 65)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha   power
    29 0.5823418  0.05 0.85698

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

jrad_stai_p

# fast_tiffsave(jrad_stai_p, path = savefolder)

es_df %>%
  rstatix::get_summary_stats(StaiS_pre, StaiS_post, type = "common")

es_stai_p <-
  es_df %>%
    dplyr::select(ID, StaiS_pre, StaiS_post) %>%
    tidyr::pivot_longer(cols = c(StaiS_pre, StaiS_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "MS",
      ylab = "Anxiety",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 20, to = 65, by = 5),
      limits = c(19, 65)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha     power
    26 0.3049887  0.05 0.3214605

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

es_stai_p

# fast_tiffsave(es_stai_p, path = savefolder)

# PA
jrad_df %>%
  rstatix::get_summary_stats(PA_pre, PA_post, type = "common")

jrad_pa_p <- 
  jrad_df %>%
    dplyr::select(ID, PA_pre, PA_post) %>%
    tidyr::pivot_longer(cols = c(PA_pre, PA_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "RP",
      ylab = "Positive Affect",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 20, to = 50, by = 5),
      limits = c(17, 51)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha     power
    37 0.4217271  0.05 0.7041574

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

jrad_pa_p

# fast_tiffsave(jrad_pa_p, path = savefolder)

es_df %>%
  rstatix::get_summary_stats(PA_pre, PA_post, type = "common")

es_pa_p <- 
  es_df %>%
    dplyr::select(ID, PA_pre, PA_post) %>%
    tidyr::pivot_longer(cols = c(PA_pre, PA_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "MS",
      ylab = "Positive Affect",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 20, to = 50, by = 5),
      limits = c(17, 51)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha    power
    36 0.2987284  0.05 0.414293

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

es_pa_p

# fast_tiffsave(es_pa_p, path = savefolder)

# NA  
jrad_df %>%
  rstatix::get_summary_stats(NA_pre, NA_post, type = "common")

jrad_na_p <- 
  jrad_df %>%
    dplyr::select(ID, NA_pre, NA_post) %>%
    tidyr::pivot_longer(cols = c(NA_pre, NA_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "RP",
      ylab = "Negative Affect",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 10, to = 30, by = 5),
      limits = c(8, 30)
    ) 

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha    power
    37 0.6723652  0.05 0.978232

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

jrad_na_p

# fast_tiffsave(jrad_na_p, path = savefolder)

es_df %>%
  rstatix::get_summary_stats(NA_pre, NA_post, type = "common")

es_na_p <- 
  es_df %>%
    dplyr::select(ID, NA_pre, NA_post) %>%
    tidyr::pivot_longer(cols = c(NA_pre, NA_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "MS",
      ylab = "Negative Affect",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 10, to = 30, by = 5),
      limits = c(8, 30)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha     power
    36 0.4970078  0.05 0.8262263

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

es_na_p

# fast_tiffsave(es_na_p, path = savefolder)

# APS
jrad_df %>%
  rstatix::get_summary_stats(APS_pre, APS_post, type = "common")

jrad_aps_p <- 
  jrad_df %>%
    dplyr::select(ID, APS_pre, APS_post) %>%
    tidyr::pivot_longer(cols = c(APS_pre, APS_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "RP",
      ylab = "Prosocial Attitudes",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 20, to = 80, by = 10),
      limits = c(22, 80)
    ) 

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha     power
    38 0.5433895  0.05 0.9034706

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

jrad_aps_p

# fast_tiffsave(jrad_aps_p, path = savefolder)

es_df %>%
  rstatix::get_summary_stats(APS_pre, APS_post, type = "common")

es_aps_p <- 
  es_df %>%                                        # id 28 is a clear outlier
    dplyr::select(ID, APS_pre, APS_post) %>%
    tidyr::pivot_longer(cols = c(APS_pre, APS_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "",
      xlab = "MS",
      ylab = "Prosocial Attitudes",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 20, to = 80, by = 10),
      limits = c(22, 80)
    ) 

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha     power
    35 0.1239919  0.05 0.1099659

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

es_aps_p

# fast_tiffsave(es_aps_p, path = savefolder)


es_df %>%                                        # id 28 is a clear outlier
    dplyr::select(ID, APS_pre, APS_post) %>%
    dplyr::filter(ID != 28) %>%
    tidyr::pivot_longer(cols = c(APS_pre, APS_post), names_to = "Time", values_to = "value") %>%
    dplyr::mutate(
      Time = stringr::str_replace(Time, "pre", "Pre"),
      Time = stringr::str_replace(Time, "post", "Post"), 
      Time = stringr::str_remove(Time, ".*_"),
      Time = stringr::str_remove(Time, "_"), 
      Time = factor(Time, levels = c("Pre", "Post"))) %>%
    group_by(ID) %>%
    my_ggwithinstats2(
      x = Time,
      y = value,
      outlier.label = ID, 
      title = "ES APS - Excluded outlier 28",
      xlab = "MS",
      ylab = "Prosocial Attitudes",
      results.subtitle = TRUE,
      centrality.label.args = TRUE,
      breaks = seq(from = 20, to = 80, by = 10),
      limits = c(22, 80)
    )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Paired t-test

     n         d alpha     power
    34 0.4505963  0.05 0.7225248

NOTE: n is number of *pairs*
URL: http://psychstat.org/ttest

7 Within-participant mediation APS

7.1 JRAD

library(JSmediation)

set.seed(101)
withinmed_jrad <- 
  mdt_within_wide(data = jrad_df,
                  DV_B = APS_pre, DV_A = APS_post,     # mdt_within_wide uses M_A - M_B and DV_A - DV_B in these models
                  M_B = PA_pre, M_A = PA_post)
add_index(withinmed_jrad)

Test of mediation (within-participant mediation)
==============================================

Variables:

- DV difference: APS_post - APS_pre 
- M difference: PA_post - PA_pre 

Paths:

==== ============== ===== ====================== Path Point estimate SE APA
==== ============== ===== ====================== a 2.000 0.763 t(36) = 2.62, p = .013 b 0.363 0.160 t(34) = 2.26, p = .030 c 2.500 0.731 t(37) = 3.42, p = .002 c’ 1.842 0.783 t(34) = 2.35, p = .025 ==== ============== ===== ======================

Indirect effect index:

- type: Within-participant indirect effect 
- point estimate: 0.726 
- confidence interval:
  - method: Monte Carlo (5000 iterations)
  - level: 0.05 
  - CI: [0.0415; 1.69]

Fitted models:

- 1 -> DV_diff 
- 1 -> M_diff 
- 1 + M_diff + M_mean -> DV_diff 

7.2 ES

set.seed(101)
# remember that id 28 is outlier on aps ... mediation does not hold even when 28 is excluded
withinmed_es <- 
  mdt_within_wide(data = es_df,   # %>% dplyr::filter(ID != 28)
                  DV_B = APS_pre, DV_A = APS_post, 
                  M_B = PA_pre, M_A = PA_post)
add_index(withinmed_es) 

Test of mediation (within-participant mediation)
==============================================

Variables:

- DV difference: APS_post - APS_pre 
- M difference: PA_post - PA_pre 

Paths:

==== ============== ===== ====================== Path Point estimate SE APA
==== ============== ===== ====================== a 1.194 0.652 t(35) = 1.83, p = .075 b 0.065 0.292 t(32) = 0.22, p = .827 c 0.771 1.028 t(34) = 0.75, p = .458 c’ 0.700 1.100 t(32) = 0.64, p = .529 ==== ============== ===== ======================

Indirect effect index:

- type: Within-participant indirect effect 
- point estimate: 0.077 
- confidence interval:
  - method: Monte Carlo (5000 iterations)
  - level: 0.05 
  - CI: [-0.683; 0.933]

Fitted models:

- 1 -> DV_diff 
- 1 -> M_diff 
- 1 + M_diff + M_mean -> DV_diff 

8 rmAnova

# id 21 twice in JRAD
# set second round to "ES"
long_df_new <- long_df
long_df_new[long_df_new$ID == 21,]$Conditie[2] <- "ES"

# Analyses
long_df_new %>%
  select_to_longer("OX") %>%
  tw_rmANOVA_func(id_var = ID, cond_var = Cond, time_var = PrePost, value_var = OX)

 Outliers 

 Normality assumption (p >.05) 

 Two-way rmANOVA 
ANOVA Table (type II tests)

        Effect DFn DFd     SSn   SSd        F                                              p p<.05   ges
1  (Intercept)   1  33 142.539 0.513 9160.579 0.00000000000000000000000000000000000000000628     * 0.987
2         Cond   1  33   0.009 0.566    0.535 0.46999999999999997335464740899624302983283997       0.005
3      PrePost   1  33   0.017 0.422    1.330 0.25700000000000000621724893790087662637233734       0.009
4 Cond:PrePost   1  33   0.045 0.428    3.472 0.07099999999999999367172875963660771958529949       0.023

 Cohen's f (partial) 
Used formula:  OX ~ Cond * PrePost + Error(ID / (Cond * PrePost))# Effect Size for ANOVA (Type II)

Parameter    | Cohen's f (partial) |       90% CI
-------------------------------------------------
Cond         |                0.13 | [0.00, 0.41]
PrePost      |                0.20 | [0.00, 0.49]
Cond:PrePost |                0.32 | [0.00, 0.62]

 Following reporting procedure for non-significant two-way interaction 

 Comparisons for treatment variable 

 Comparisons for time variable 

 Post-hoc Power Analysis for Repeated-measures ANOVA 
[1] "Power analysis for within-effect test"

[1] "Power analysis for within-effect test"

[1] "Power analysis for interaction-effect test"

long_df_new %>%
  select_to_longer("StaiS") %>%
  tw_rmANOVA_func(id_var = ID, cond_var = Cond, time_var = PrePost, value_var = StaiS)

 Outliers 

 Normality assumption (p >.05) 

 Two-way rmANOVA 
ANOVA Table (type II tests)

        Effect DFn DFd       SSn      SSd       F                p p<.05      ges
1  (Intercept)   1  19 80327.812 5788.438 263.668 0.00000000000135     * 0.914000
2         Cond   1  19    59.512  603.737   1.873 0.18700000000000       0.008000
3      PrePost   1  19   201.612  677.638   5.653 0.02800000000000     * 0.026000
4 Cond:PrePost   1  19     2.813  481.438   0.111 0.74300000000000       0.000372

 Cohen's f (partial) 
Used formula:  StaiS ~ Cond * PrePost + Error(ID / (Cond * PrePost))# Effect Size for ANOVA (Type II)

Parameter    | Cohen's f (partial) |       90% CI
-------------------------------------------------
Cond         |                0.31 | [0.00, 0.70]
PrePost      |                0.55 | [0.13, 0.94]
Cond:PrePost |                0.08 | [0.00, 0.42]

 Following reporting procedure for non-significant two-way interaction 

 Comparisons for treatment variable 

 Comparisons for time variable 

 Post-hoc Power Analysis for Repeated-measures ANOVA 
[1] "Power analysis for within-effect test"

[1] "Power analysis for within-effect test"

[1] "Power analysis for interaction-effect test"

long_df_new %>%
  select_to_longer("PA") %>%
  tw_rmANOVA_func(id_var = ID, cond_var = Cond, time_var = PrePost, value_var = PA)

 Outliers 

 Normality assumption (p >.05) 

 Two-way rmANOVA 
ANOVA Table (type II tests)

        Effect DFn DFd        SSn      SSd        F                               p p<.05       ges
1  (Intercept)   1  31 138206.531 3177.469 1348.370 0.00000000000000000000000000406     * 0.9700000
2         Cond   1  31      7.031  604.969    0.360 0.55300000000000004707345624411       0.0020000
3      PrePost   1  31     87.781  287.219    9.474 0.00400000000000000008326672685     * 0.0200000
4 Cond:PrePost   1  31      0.281  236.719    0.037 0.84899999999999997690736108780       0.0000653

 Cohen's f (partial) 
Used formula:  PA ~ Cond * PrePost + Error(ID / (Cond * PrePost))# Effect Size for ANOVA (Type II)

Parameter    | Cohen's f (partial) |       90% CI
-------------------------------------------------
Cond         |                0.11 | [0.00, 0.40]
PrePost      |                0.55 | [0.23, 0.87]
Cond:PrePost |                0.03 | [0.00, 0.27]

 Following reporting procedure for non-significant two-way interaction 

 Comparisons for treatment variable 

 Comparisons for time variable 

 Post-hoc Power Analysis for Repeated-measures ANOVA 
[1] "Power analysis for within-effect test"

[1] "Power analysis for within-effect test"

[1] "Power analysis for interaction-effect test"

long_df_new %>%
  dplyr::rename(NegAf_pre = NA_pre, NegAf_post = NA_post) %>%
  select_to_longer("NegAf") %>%
  tw_rmANOVA_func(id_var = ID, cond_var = Cond, time_var = PrePost, value_var = NegAf)

 Outliers 

 Normality assumption (p >.05) 

 Two-way rmANOVA 
ANOVA Table (type II tests)

        Effect DFn DFd       SSn      SSd       F                       p p<.05       ges
1  (Intercept)   1  31 22578.125 1662.875 420.911 0.000000000000000000135     * 0.8990000
2         Cond   1  31     0.125  611.875   0.006 0.937000000000000055067       0.0000493
3      PrePost   1  31    94.531  168.469  17.395 0.000226999999999999988     * 0.0360000
4 Cond:PrePost   1  31     0.781   93.219   0.260 0.613999999999999990230       0.0003080

 Cohen's f (partial) 
Used formula:  NegAf ~ Cond * PrePost + Error(ID / (Cond * PrePost))# Effect Size for ANOVA (Type II)

Parameter    | Cohen's f (partial) |       90% CI
-------------------------------------------------
Cond         |                0.01 | [0.00, 0.12]
PrePost      |                0.75 | [0.41, 1.08]
Cond:PrePost |                0.09 | [0.00, 0.38]

 Following reporting procedure for non-significant two-way interaction 

 Comparisons for treatment variable 

 Comparisons for time variable 

 Post-hoc Power Analysis for Repeated-measures ANOVA 
[1] "Power analysis for within-effect test"

[1] "Power analysis for within-effect test"

[1] "Power analysis for interaction-effect test"

long_df_new %>%
  select_to_longer("APS") %>%
  tw_rmANOVA_func(id_var = ID, cond_var = Cond, time_var = PrePost, value_var = APS)

 Outliers 

 Normality assumption (p >.05) 

 Two-way rmANOVA 
ANOVA Table (type II tests)

        Effect DFn DFd        SSn       SSd        F                              p p<.05      ges
1  (Intercept)   1  31 488566.125 11756.375 1288.284 0.0000000000000000000000000081     * 0.969000
2         Cond   1  31      5.281  2681.219    0.061 0.8060000000000000497379915032       0.000343
3      PrePost   1  31     75.031   529.469    4.393 0.0439999999999999974464870434     * 0.005000
4 Cond:PrePost   1  31     32.000   416.500    2.382 0.1330000000000000071054273576       0.002000

 Cohen's f (partial) 
Used formula:  APS ~ Cond * PrePost + Error(ID / (Cond * PrePost))# Effect Size for ANOVA (Type II)

Parameter    | Cohen's f (partial) |       90% CI
-------------------------------------------------
Cond         |                0.04 | [0.00, 0.30]
PrePost      |                0.38 | [0.04, 0.68]
Cond:PrePost |                0.28 | [0.00, 0.58]

 Following reporting procedure for non-significant two-way interaction 

 Comparisons for treatment variable 

 Comparisons for time variable 

 Post-hoc Power Analysis for Repeated-measures ANOVA 
[1] "Power analysis for within-effect test"

[1] "Power analysis for within-effect test"

[1] "Power analysis for interaction-effect test"

library(lavaan)
library(semPlot)

# model=
# "
#   # Regressions
#   m1 ~ a*iv
#   m2 ~ b*m1 + iv
#   dv ~ c*m2 + m1 + d*iv + cv1 + cv2
#   
#   # Defined Parameters:
#   ie := a*b*c
#   de := d
# "


model <- "
  # Regressions
  PA_pre ~ a*OX_post
  PA_post ~ b*PA_pre + OX_post
  APS_post ~ c*PA_post + PA_pre + d*OX_post + StaiS_post
  
  # Defined Parameters:
  ie := a*b*c
  de := d
"

fit <- lavaan::sem(model, data = jrad_df)  # long_df_all %>% dplyr::filter(ID != 28)
summary(fit)

lavaanPlot::lavaanPlot(model = fit, 
                       graph_options = list(layout = "dot", rankdir = "LR"), 
                       node_options = list(shape = "box", fontname = "Helvetica"), 
                       edge_options = list(color = "grey"), 
                       coefs = TRUE, stand = TRUE, covs = TRUE, 
                       sig = 0.05, stars = c("regress", "cov"))  # unstandardized: sig = 1.00 

–>

source("C:/Users/Mihai/Desktop/LCS mediation/Function LCS ANCOVA mediation.R")
mod <- lcs_ancova_med(df = long_df_all, x = Conditie, y1 = APS_pre, y2 = APS_post, m1 = PA_pre, m2 = PA_post)
summary(mod)
semPlot::semPaths(mod, layout = "spring",  nCharNodes = 0, nCharEdges = 0, what = "path", whatLabels = "path", edge.label.cex = 0.8)

–>

9 Session Info

R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8.1 x64 (build 9600)

Matrix products: default

locale:
[1] LC_COLLATE=Romanian_Romania.1250  LC_CTYPE=Romanian_Romania.1250    LC_MONETARY=Romanian_Romania.1250
[4] LC_NUMERIC=C                      LC_TIME=Romanian_Romania.1250    
system code page: 1252

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] effectsize_0.4.5           afex_1.1-1                 WebPower_0.8.6             PearsonDS_1.2             
 [5] lavaan_0.6-9               lme4_1.1-29                Matrix_1.3-4               rlang_1.0.2               
 [9] JSmediation_0.2.0          MASS_7.3-54                rio_0.5.29                 scales_1.2.0              
[13] ggpubr_0.4.0               summarytools_1.0.0         rstatix_0.7.0              broom_0.7.11              
[17] PerformanceAnalytics_2.0.4 xts_0.12.1                 zoo_1.8-9                  psych_2.1.9               
[21] plyr_1.8.6                 forcats_0.5.1              stringr_1.4.0              dplyr_1.0.9               
[25] purrr_0.3.4                readr_2.0.1                tidyr_1.1.3                tibble_3.1.7              
[29] ggplot2_3.3.5              tidyverse_1.3.1            papaja_0.1.0.9997          pacman_0.5.1              

loaded via a namespace (and not attached):
  [1] readxl_1.3.1              pairwiseComparisons_3.1.6 backports_1.2.1           splines_4.1.0            
  [5] gmp_0.6-2                 kSamples_1.2-9            ipmisc_6.0.2              TH.data_1.0-10           
  [9] pryr_0.1.5                digest_0.6.28             SuppDists_1.1-9.7         htmltools_0.5.2          
 [13] magick_2.7.3              lmerTest_3.1-3            fansi_0.5.0               magrittr_2.0.1           
 [17] checkmate_2.0.0           memoise_2.0.0             paletteer_1.4.0           tzdb_0.1.2               
 [21] openxlsx_4.2.4            modelr_0.1.8              matrixStats_0.60.1        sandwich_3.0-1           
 [25] colorspace_2.0-3          rvest_1.0.2               ggrepel_0.9.1             haven_2.4.3              
 [29] xfun_0.30                 tcltk_4.1.0               crayon_1.5.1              jsonlite_1.8.0           
 [33] zeallot_0.1.0             survival_3.2-13           glue_1.6.2                gtable_0.3.0             
 [37] emmeans_1.6.3             MatrixModels_0.5-0        statsExpressions_1.1.0    car_3.0-11               
 [41] Rmpfr_0.8-4               abind_1.4-5               rapportools_1.0           mvtnorm_1.1-2            
 [45] DBI_1.1.1                 PMCMRplus_1.9.0           Rcpp_1.0.7                performance_0.7.3        
 [49] xtable_1.8-4              tmvnsim_1.0-2             foreign_0.8-81            stats4_4.1.0             
 [53] datawizard_0.2.0.1        httr_1.4.3                ellipsis_0.3.2            pkgconfig_2.0.3          
 [57] reshape_0.8.8             farver_2.1.0              multcompView_0.1-8        dbplyr_2.1.1             
 [61] utf8_1.2.2                reshape2_1.4.4            tidyselect_1.1.2          labeling_0.4.2           
 [65] munsell_0.5.0             cellranger_1.1.0          tools_4.1.0               cachem_1.0.6             
 [69] ggprism_1.0.3             cli_3.0.1                 generics_0.1.2            evaluate_0.15            
 [73] fastmap_1.1.0             BWStest_0.2.2             rematch2_2.1.2            knitr_1.39               
 [77] fs_1.5.2                  zip_2.2.0                 pander_0.6.5              WRS2_1.1-3               
 [81] pbapply_1.4-3             nlme_3.1-152              xml2_1.3.3                correlation_0.7.0        
 [85] compiler_4.1.0            rstudioapi_0.13           curl_4.3.2                ggsignif_0.6.2           
 [89] reprex_2.0.1              pbivnorm_0.6.0            stringi_1.7.4             highr_0.9                
 [93] parameters_0.14.0         lattice_0.20-44           nloptr_2.0.0              vctrs_0.4.1              
 [97] pillar_1.7.0              lifecycle_1.0.1           mc2d_0.1-21               estimability_1.3         
[101] data.table_1.14.0         insight_0.14.4            patchwork_1.1.1           R6_2.5.1                 
[105] BayesFactor_0.9.12-4.2    codetools_0.2-18          boot_1.3-28               gtools_3.9.2             
[109] assertthat_0.2.1          withr_2.5.0               mnormt_2.0.2              multcomp_1.4-17          
[113] bayestestR_0.11.0         hms_1.1.1                 quadprog_1.5-8            grid_4.1.0               
[117] minqa_1.2.4               coda_0.19-4               rmarkdown_2.14            carData_3.0-4            
[121] numDeriv_2016.8-1.1       lubridate_1.7.10          base64enc_0.1-3           ggstatsplot_0.8.0        

 

---

A work by Claudiu Papasteri

claudiu.papasteri@gmail.com