M2

Report

1 Define functions

## Define function that recodes to numeric, but watches out to coercion to not introduce NAs
colstonumeric <- function(df){
  tryCatch({
    df_num <- as.data.frame(
      lapply(df,
             function(x) { as.numeric(as.character(x))})) 
  },warning = function(stop_on_warning) {
    message("Stoped the execution of numeric conversion: ", conditionMessage(stop_on_warning))
  }) 
}
##
## Define function that reverse codes items
ReverseCode <- function(df, tonumeric = FALSE, min = NULL, max = NULL) {
  if(tonumeric) df <- colstonumeric(df)
  df <- (max + min) - df
}
##
## Define function that scores only rows with less than 10% NAs (returns NA if all or above threshold percentage of rows are NA); can reverse code if vector of column indexes and min, max are provided.
ScoreLikert <- function(df, napercent = .1, tonumeric = FALSE, reversecols = NULL, min = NULL, max = NULL) {
  reverse_list <- list(reversecols = reversecols, min = min, max = max)
  reverse_check <- !sapply(reverse_list, is.null)
  
  # Recode to numeric, but watch out to coercion to not introduce NAs
  colstonumeric <- function(df){
    tryCatch({
      df_num <- as.data.frame(
        lapply(df,
               function(x) { as.numeric(as.character(x))})) 
    },warning = function(stop_on_warning) {
      message("Stoped the execution of numeric conversion: ", conditionMessage(stop_on_warning))
    }) 
  }
  
  if(tonumeric) df <- colstonumeric(df)
  
  if(all(reverse_check)){
    df[ ,reversecols] <- (max + min) - df[ ,reversecols]
  }else if(any(reverse_check)){
    stop("Insuficient info for reversing. Please provide: ", paste(names(reverse_list)[!reverse_check], collapse = ", "))
  }
  
  ifelse(rowSums(is.na(df)) > ncol(df) * napercent,
         NA,
         rowSums(df, na.rm = TRUE) * NA ^ (rowSums(!is.na(df)) == 0)
  )
}
##

## Func t test si boxplot simplu
func_t_box <- function(df, ind, pre_var, post_var){
  vars <- c(ind, pre_var, post_var)                # to avoid new tidyverse error of ambiguity due to external vectos
  df_modif <-
    df %>%
    dplyr::select(tidyselect::all_of(vars)) %>% 
    tidyr::drop_na() %>%
    gather(tidyselect::all_of(pre_var), tidyselect::all_of(post_var), key = "Cond", value = "value") %>% 
    mutate_at(vars(c(1, 2)), as.factor) %>% 
    mutate(Cond = factor(Cond, levels = c(pre_var, post_var))) 
  
  stat_comp <- ggpubr::compare_means(value ~ Cond, data = df_modif, method = "t.test", paired = TRUE)
  
  stat_comp2 <-
    df_modif %>% 
    do(tidy(t.test(.$value ~ .$Cond,
                   paired = TRUE,
                   data=.)))
  
  plot <- 
    ggpubr::ggpaired(df_modif, x = "Cond", y = "value", id = ind, 
                     color = "Cond", line.color = "gray", line.size = 0.4,
                     palette = c("#00AFBB", "#FC4E07"), legend = "none") +
      stat_summary(fun.data = mean_se,  colour = "darkred") +
      ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label.x = as.numeric(df_modif$Cond)-0.4, label.y = max(df_modif$value)+0.5) + 
      ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label = "p.signif", comparisons = list(c(pre_var, post_var)))
  
  cat(paste0("#### ", pre_var, " ", post_var, "\n", "\n"))
  print(stat_comp)
  print(stat_comp2)
  print(plot)
}

# library(tidyverse)
# library(ggpubr)
# library(rstatix)
# library(broom)
# library(emmeans)
# library(rlang)


# Function ANCOVAPost
# Takes Long Data
ANCOVAPost_func <- 
  function(data, id_var, 
           time_var, pre_label, post_label,
           value_var = scores, cond_var, 
           assum_check = TRUE, posthoc = TRUE, 
           p_adjust_method = "bonferroni"){
  
  id_var_enq <- rlang::enquo(id_var)
  id_var_name <- rlang::as_name(id_var_enq)    
  time_var_enq <- rlang::enquo(time_var)
  time_var_name <- rlang::as_name(time_var_enq)
  pre_var_enq <- rlang::enquo(pre_label)
  pre_var_name <- rlang::as_name(pre_var_enq)  
  post_var_enq <- rlang::enquo(post_label)
  post_var_name <- rlang::as_name(post_var_enq)
  cond_var_enq <- rlang::enquo(cond_var)
  cond_var_name <- rlang::as_name(cond_var_enq)
  value_var_enq <- rlang::enquo(value_var)
  value_var_name <- rlang::as_name(value_var_enq) 
  
  data_wider <-
    data %>%
    dplyr::select(!!id_var_enq, !!time_var_enq, !!cond_var_enq, !!value_var_enq) %>%
    spread(key = time_var_name, value = value_var_name) # %>%     # if need to compute change score statistics go from here
    # mutate(difference = !!post_var_enq - !!pre_var_enq)  
    
  # Assumptions
  if(assum_check){
  cat("\n Linearity assumptionLinearity assumption (linear relationship between pre-test and post-test for each group) \n")
  # Create a scatter plot between the covariate (i.e., pretest) and the outcome variable (i.e., posttest)
  scatter_lin <-  
    ggscatter(data_wider, x = pre_var_name, y = post_var_name, color = cond_var_name, 
              add = "reg.line", title = "Linearity assumption") +
      stat_regline_equation(aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = !!cond_var_enq))
  
  cat("\n Homogeneity of regression slopes (interaction term is n.s.) \n")
  data_wider %>% 
    anova_test(as.formula(paste0(post_var_name, " ~ ", cond_var_name, " * ", pre_var_name))) %>% 
    print()   
  
  cat("\n Normality of residuals (Model diagnostics & Shapiro Wilk) \n")
  # Fit the model, the covariate goes first
  model <- 
    lm(as.formula(paste0(post_var_name, " ~ ", cond_var_name, " + ", pre_var_name)),
       data = data_wider)
  cat("\n Inspect the model diagnostic metrics \n")
  model.metrics <- 
    augment(model) %>%
    dplyr::select(-.hat, -.sigma, -.fitted, -.se.fit)  %>%    # Remove details
    print()
  cat("\n Normality of residuals (Shapiro Wilk p>.05) \n")
  shapiro_test(model.metrics$.resid) %>% 
    print()
  
  cat("\n Homogeneity of variances (Levene’s test p>.05) \n")
  model.metrics %>% 
    levene_test(as.formula(paste0(".resid", " ~ ", cond_var_name)) ) %>%     
    print()
  
  cat("\n Outliers (needs to be 0) \n")
  model.metrics %>% 
    filter(abs(.std.resid) > 3) %>%
    as.data.frame() %>% 
    print()
  
  }
  
  cat("\n ANCOVAPost \n")
  res_ancova <- 
    data_wider %>% 
    anova_test(as.formula(paste0(post_var_name, " ~ ",  pre_var_name, " + ", cond_var_name)))       # the covariate needs to be first term 
  get_anova_table(res_ancova) %>% print()
  
  cat("\n Pairwise comparisons \n")
  pwc <- 
    data_wider %>% 
    emmeans_test(as.formula(paste0(post_var_name, " ~ ", cond_var_name)),       
                 covariate = !!pre_var_enq,
                 p.adjust.method = p_adjust_method)
  pwc %>% print()
  cat("\n Display the adjusted means of each group, also called as the estimated marginal means (emmeans) \n")
  get_emmeans(pwc) %>% print()
  
  # Visualization: line plots with p-values
  pwc <- 
    pwc %>% 
    add_xy_position(x = cond_var_name, fun = "mean_se")
  
  line_plot <- 
    ggline(get_emmeans(pwc), x = cond_var_name, y = "emmean") +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) + 
      stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) +
      labs(subtitle = get_test_label(res_ancova, detailed = TRUE),
           caption = get_pwc_label(pwc))
  
  if(assum_check){
    list(scatter_lin, line_plot)
  }else{
    line_plot
  }
  
}

# ex.
# ANCOVAPost_func(new_anxiety, time_var = time, pre_label = pretest, post_label = posttest,
#           value_var = scores, cond_var = group, assum_check = TRUE, p_adjust_method = "bonferroni")

# library(tidyverse)
# library(ggpubr)
# library(rstatix)
# library(rlang)

# Define Function for Mixed Anova
tw_mixedANOVA_func <- 
  function(data, id_var, cond_var, time_var, value_var, 
           assum_check = TRUE, posthoc_sig_interac = FALSE, posthoc_ns_interac = FALSE,
           p_adjust_method = "bonferroni"){
    
    # input dataframe needs to have columns names diffrent from "variable" and "value" because it collides with rstatix::shapiro_test
    
    id_var_enq <- rlang::enquo(id_var)  
    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)
    
    data <-                  # need to subset becuase of strange contrasts error in anova_test()
      data %>%
      dplyr::select(!!id_var_enq, !!time_var_enq, !!cond_var_enq, !!value_var_enq) 
    
    # Assumptions
    if(assum_check){
      cat("\n Outliers \n")
      data %>%
        dplyr::group_by(!!time_var_enq, !!cond_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(!!time_var_enq, !!cond_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
      
      cat("\n Homogneity of variance assumption - Levene’s test (p>.05) \n")
      data %>%
        group_by(!!time_var_enq ) %>%
        levene_test(as.formula(paste0(value_var_name, " ~ ", cond_var_name))) %>%
        print()
      
      cat("\n Homogeneity of covariances assumption - Box’s test of equality of covariance matrices (p>.001) \n")
      box_m(data = data[, value_var_name, drop = FALSE], group = data[, cond_var_name, drop = TRUE]) %>%
        print
    }
    
    # Two-way rmANOVA - check for interaction (ex. F(2, 22) = 30.4, p < 0.0001)
    cat("\n Two-way rmANOVA \n")
    res_aov <- 
      anova_test(data = data, dv = !!value_var_enq, wid = !!id_var_enq,               # automatically does sphericity Mauchly’s test
                 within = !!time_var_enq, between = !!cond_var_enq)
    get_anova_table(res_aov) %>%  # ges: Greenhouse-Geisser sphericity correction is automatically applied to factors violating the sphericity assumption  
      print()
    
    
    # ------------------------------------------------------------------------
    
    #- Procedure for a significant two-way interaction -
    if(posthoc_sig_interac){
      cat("\n Effect of group at each time point - One-way ANOVA\n")
      one_way <- 
        data %>%
        group_by(!!time_var_enq) %>%
        anova_test(dv = !!value_var_enq, wid = !!id_var_enq, between = !!cond_var_enq) %>%
        get_anova_table() %>%
        adjust_pvalue(method = p_adjust_method)
      one_way %>% print()
      
      cat("\n Pairwise comparisons between group levels \n")
      pwc <- 
        data %>%
        group_by(!!time_var_enq) %>%
        pairwise_t_test(as.formula(paste0(value_var_name, " ~ ", cond_var_name)), 
                        p.adjust.method = p_adjust_method)
      pwc %>% print()
      cat("\n Effect of time at each level of exercises group  - One-way ANOVA \n")
      one_way2 <- 
        data %>%
        group_by(!!cond_var_enq) %>%
        anova_test(dv = !!value_var_enq, wid = !!id_var_enq, within = !!time_var_enq) %>%
        get_anova_table() %>%
        adjust_pvalue(method = p_adjust_method)
      one_way2 %>% print()
      
      cat("\n Pairwise comparisons between time points at each group levels (we have repeated measures by time) \n")
      pwc2 <-
        data %>%
        group_by(!!cond_var_enq) %>%
        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 %>%
        pairwise_t_test(
          as.formula(paste0(value_var_name, " ~ ", cond_var_name)),     # paste formula, not quosure             
          paired = FALSE,
          p.adjust.method = p_adjust_method
        )
      pwc_cond %>% print()
      cat("\n Comparisons for time variable \n")
      pwc_time <-
        data %>% 
        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()
    }
    
    # Visualization
    bx_plot <- 
      ggboxplot(data, x = time_var_name, y = value_var_name,
                color = cond_var_name, palette = "jco")
    pwc <- 
      pwc %>% 
      add_xy_position(x = time_var_name)
    bx_plot <- 
      bx_plot + 
      stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
      labs(
        subtitle = get_test_label(res_aov, detailed = TRUE),
        caption = get_pwc_label(pwc)
      )
    
    if(assum_check){ 
      list(qq_plot, bx_plot)
    }else{
      bx_plot
    } 
    
  }

# ex. - run on long format
# tw_mixedANOVA_func(data = anxiety, id_var = id, cond_var = group, time_var = time, value_var = score,
#                 posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

2 Read, Clean, Recode

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Read, Clean, Recode, Unite
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Read files
folder <- "C:/Users/Mihai/Desktop/R Notebooks/notebooks/M2-report"
file <- "DateM2 final.xlsx"

setwd(folder)
Data_pre <- rio::import(file.path(folder, file),
                           skip = 3, which = "M2 PRE")

Data_post <- rio::import(file.path(folder, file),
                           skip = 3, which = "M2 POST")



# PRE = colnames(Data_pre); POST = colnames(Data_post) 
# cbind(PRE, POST)                                             # "s1.7" var is missing form POST -- missmatch row 12
#                                                              # also missing "s2.7" and others


# index from "DateM2.xlsx"
# index_ryff_pre <- 94:135
# index_pss_pre <- 136:149
# index_ryff_post <- 94:135 - 5
# index_pss_post <- 136:149 - 5

index_ryff_pre <- 83:124
index_pss_pre <- 125:138
index_ryff_post <- 81:122
index_pss_post <- 123:136


colnames(Data_pre)[1] <- "ID"
colnames(Data_pre)[colnames(Data_pre) == "Stres pre"] <- "VAS_stress_pre"
colnames(Data_pre)[colnames(Data_pre) == "Stres post"] <- "VAS_stress_post"
colnames(Data_pre)[index_ryff_pre] <- sprintf("ryff_%01d", seq(1, 42))
colnames(Data_pre)[index_pss_pre] <- sprintf("pss_%01d", seq(1, 14))
Data_pre <-
  Data_pre %>%
  drop_na(ID) %>%
  dplyr::mutate_if(is.character, list(~dplyr::na_if(., "na"))) 
Data_pre[index_ryff_pre] <- colstonumeric(Data_pre[index_ryff_pre])
Data_pre$pss_8[Data_pre$pss_8 == "++++++"] <- NA                           # typo
Data_pre[index_pss_pre] <- colstonumeric(Data_pre[index_pss_pre])

colnames(Data_post)[1] <- "ID"
colnames(Data_post)[colnames(Data_post) == "Stres pre"] <- "VAS_stress_pre"
colnames(Data_post)[colnames(Data_post) == "Stres post"] <- "VAS_stress_post"
colnames(Data_post)[index_ryff_post] <- sprintf("ryff_%01d", seq(1, 42))
colnames(Data_post)[index_pss_post] <- sprintf("pss_%01d", seq(1, 14))
Data_post <-
  Data_post %>%
  dplyr::mutate_if(is.character, list(~dplyr::na_if(., "na"))) 
Data_post[index_ryff_post] <- colstonumeric(Data_post[index_ryff_post])
Data_post[index_pss_post] <- colstonumeric(Data_post[index_pss_post])

# typos
# check_numeric <- as.data.frame(sapply(Data_pre[index_pss_pre], varhandle::check.numeric)) 
# sapply(check_numeric, function(x) length(which(!x)))

colnames(Data_pre)[colnames(Data_pre) == "s1.1."] <- "s1.1"
colnames(Data_post)[colnames(Data_post) == "s1.1."] <- "s1.1"

Data_post$s2.10[which(Data_post$s2.10 == ".")] <- NA
Data_post$s2.10 <- as.numeric(Data_post$s2.10)
Data_post$s3.16[which(Data_post$s3.16 == "kq")] <- NA
Data_post$s3.16 <- as.numeric(Data_post$s3.16)

## PSS-SF 14 (likert 0-4)
# Items 4, 5, 6, 7, 9, 10, and 13 are scored in reverse direction.

indexitem_revPSS <- c(4, 5, 6, 7, 9, 10, 13)

Data_pre[, index_pss_pre][indexitem_revPSS] <- ReverseCode(Data_pre[, index_pss_pre][indexitem_revPSS], tonumeric = FALSE, min = 0, max = 4)
Data_post[, index_pss_post][indexitem_revPSS] <- ReverseCode(Data_post[, index_pss_post][indexitem_revPSS], tonumeric = FALSE, min = 0, max = 4)

Data_pre$PSS <- ScoreLikert(Data_pre[, index_pss_pre], napercent = .4)
Data_post$PSS <- ScoreLikert(Data_post[, index_pss_post], napercent = .4)
  

## Ryff (likert 1-6)
# Recode negative phrased items:  3,5,10,13,14,15,16,17,18,19,23,26,27,30,31,32,34,36,39,41.
# Autonomy: 1,7,13,19,25,31,37
# Environmental mastery: 2,8,14,20,26,32,38
# Personal Growth: 3,9,15,21,27,33,39
# Positive Relations: 4,10,16,22,28,34,40
# Purpose in life: 5,11,17,23,29,35,41
# Self-acceptance: 6,12,18,24,30,36,42

indexitem_revRYFF <- c(3,5,10,13,14,15,16,17,18,19,23,26,27,30,31,32,34,36,39,41)

Data_pre[, index_ryff_pre][indexitem_revRYFF] <- ReverseCode(Data_pre[, index_ryff_pre][indexitem_revRYFF], tonumeric = FALSE, min = 1, max = 6)
Data_post[, index_ryff_post][indexitem_revRYFF] <- ReverseCode(Data_post[, index_ryff_post][indexitem_revRYFF], tonumeric = FALSE, min = 1, max = 6)

indexitem_Auto <- c(1,7,13,19,25,31,37)
indexitem_EnvM <- c(2,8,14,20,26,32,38)
indexitem_PersG <- c(3,9,15,21,27,33,39)
indexitem_PosRel <- c(4,10,16,22,28,34,40)
indexitem_PurLif <- c(5,11,17,23,29,35,41)
indexitem_SelfAc <- c(6,12,18,24,30,36,42)


Data_pre$Auto <- ScoreLikert(Data_pre[, index_ryff_pre][indexitem_Auto], napercent = .4)
Data_pre$EnvM <- ScoreLikert(Data_pre[, index_ryff_pre][indexitem_EnvM], napercent = .4)
Data_pre$PersG <- ScoreLikert(Data_pre[, index_ryff_pre][indexitem_PersG], napercent = .4)
Data_pre$PosRel <- ScoreLikert(Data_pre[, index_ryff_pre][indexitem_PosRel], napercent = .4)
Data_pre$PurLif <- ScoreLikert(Data_pre[, index_ryff_pre][indexitem_PurLif], napercent = .4)
Data_pre$SelfAc <- ScoreLikert(Data_pre[, index_ryff_pre][indexitem_SelfAc], napercent = .4)

Data_post$Auto <- ScoreLikert(Data_post[, index_ryff_post][indexitem_Auto], napercent = .4)
Data_post$EnvM <- ScoreLikert(Data_post[, index_ryff_post][indexitem_EnvM], napercent = .4)
Data_post$PersG <- ScoreLikert(Data_post[, index_ryff_post][indexitem_PersG], napercent = .4)
Data_post$PosRel <- ScoreLikert(Data_post[, index_ryff_post][indexitem_PosRel], napercent = .4)
Data_post$PurLif <- ScoreLikert(Data_post[, index_ryff_post][indexitem_PurLif], napercent = .4)
Data_post$SelfAc <- ScoreLikert(Data_post[, index_ryff_post][indexitem_SelfAc], napercent = .4)

## Save Scores
# rio::export(Data_pre[, c(1, 150:156)])
# rio::export(Data_post[, c(1, 149:155)])
# nlastcol <- 7
# rio::export(list(PRE = Data_pre[, c(1, (ncol(Data_pre)-nlastcol+1):ncol(Data_pre))], POST = Data_post[, c(1, (ncol(Data_post)-nlastcol+1):ncol(Data_post))]), 
#             "M2 PSS Ryff final.xlsx")

Data_pre$S1_Mean <- rowMeans(Data_pre[, sprintf("s1.%d", c(1:6, 8:16))], na.rm = TRUE)  # Data_pre[, grep("s1.", colnames(Data_pre))]
Data_pre$S2_Mean <- rowMeans(Data_pre[, sprintf("s2.%d", c(1:6, 8:16))], na.rm = TRUE)  # Data_pre[, grep("s2.", colnames(Data_pre))]
Data_pre$S3_Mean <- rowMeans(Data_pre[, sprintf("s3.%d", c(1:6, 8:16))], na.rm = TRUE)  # Data_pre[, grep("s3.", colnames(Data_pre))]

Data_post$S1_Mean <- rowMeans(Data_post[, sprintf("s1.%d", c(1:6, 8:16))], na.rm = TRUE)  # Data_post[, grep("s1.", colnames(Data_post))]
Data_post$S2_Mean <- rowMeans(Data_post[, sprintf("s2.%d", c(1:6, 8:16))], na.rm = TRUE)  # Data_post[, grep("s2.", colnames(Data_post))]
Data_post$S3_Mean <- rowMeans(Data_post[, sprintf("s3.%d", c(1:6, 8:16))], na.rm = TRUE)  # Data_post[, grep("s3.", colnames(Data_post))]

3 Add TR & CTRL groups

tr_ids <- paste0(c(1, 2, 4, 5, 6, 10, 13, 14, 15, 16, 24, 25, 26, 32, 34, 35, 39, 40, 42, 46, 47, 48, 51, 53, 59, 60, 62, 63, 64), " M2") 
ctrl_ids <- paste0(c(3, 7, 8, 9, 11, 12, 17, 18, 19, 20, 21, 22, 23, 27, 28, 30, 33, 36, 37, 38, 41, 43, 45, 49, 52, 54, 55, 56, 57, 58), " M2") 
# Subj from TR: Data_pre[which(Data_pre$ID %in% tr_ids),]  
# Subj from CTRL: Data_pre[which(Data_pre$ID %in% ctrl_ids),] 

Data_pre$Cond <- dplyr::case_when(Data_pre$ID %in% tr_ids ~ "TR",
                                  Data_pre$ID %in% ctrl_ids ~ "CTRL",
                                  TRUE ~ NA_character_)
Data_post$Cond <- dplyr::case_when(Data_post$ID %in% tr_ids ~ "TR",
                                  Data_post$ID %in% ctrl_ids ~ "CTRL",
                                  TRUE ~ NA_character_)

4 Add Count-vars & Post-Vars

names(Data_post)[137:139] <- c("spontaneous", "voluntary", "dreams")

Data_pre$s1.count <- rowSums(!is.na(Data_pre[, grepl("s1.*", names(Data_pre))]))
Data_pre$s2.count <- rowSums(!is.na(Data_pre[, grepl("s1.*", names(Data_pre))]))   
Data_pre$s3.count <- rowSums(!is.na(Data_pre[, grepl("s3.*", names(Data_pre))]))

Data_post$s1.count <- rowSums(!is.na(Data_post[, grepl("s1.*", names(Data_post))]))
Data_post$s2.count <- rowSums(!is.na(Data_post[, grepl("s1.*", names(Data_post))]))    # for some reason it has +1 count from all rest
Data_post$s3.count <- rowSums(!is.na(Data_post[, grepl("s3.*", names(Data_post))]))

5 Unite data frames

cat("## Number of subjects in pre")

## Number of subjects in pre

Data_pre %>% 
 dplyr::summarise(count = dplyr::n_distinct(ID))

cat("## Number of subjects in post")

## Number of subjects in post

Data_post %>%
 dplyr::summarise(count = dplyr::n_distinct(ID))

Data_pre$PrePost <- rep("Pre", nrow(Data_pre))
Data_post$PrePost <- rep("Post", nrow(Data_post))

Data_pre_scales <- Data_pre[, c("ID", "Varsta",
                                "VAS_stress_pre", "VAS_stress_post", "PSS", 
                                "Auto", "EnvM", "PersG", "PosRel", "PurLif",  "SelfAc", "S1_Mean", "S2_Mean", "S3_Mean", 
                                "PrePost", "Cond",
                                "s1.count", "s2.count", "s3.count")]
Data_post_scales <- Data_post[, c("ID", "Varsta",
                                  "VAS_stress_pre", "VAS_stress_post", "PSS", 
                                  "Auto", "EnvM", "PersG", "PosRel", "PurLif",  "SelfAc", "S1_Mean", "S2_Mean", "S3_Mean", 
                                  "PrePost", "Cond",
                                  "s1.count", "s2.count", "s3.count",
                                  "spontaneous", "voluntary", "dreams")]

Data_unif_long <- dplyr::bind_rows(Data_pre_scales, Data_post_scales)
Data_unif_wide <-
  Data_unif_long %>%
  tidyr::pivot_wider(id_cols = ID, names_from = PrePost, values_from = c(VAS_stress_pre, VAS_stress_post, PSS, Auto, EnvM, PersG, PosRel, PurLif, SelfAc, Cond,
                                                           S1_Mean, S2_Mean, S3_Mean, s1.count, s2.count, s3.count, spontaneous, voluntary, dreams))

Data_unif_long <-
  Data_unif_long %>%
    dplyr::mutate(ID = as.factor(ID),
                  Cond = as.factor(Cond),
                  PrePost = as.factor(PrePost)) %>%
   dplyr::mutate(Vas_Diff = VAS_stress_post - VAS_stress_pre,
                 Vas_Mean = rowMeans(dplyr::select(.data = ., VAS_stress_post, VAS_stress_pre), na.rm = TRUE))




# Export data
# rio::export(Data_unif_wide, "M2 wide_format_2.xlsx")
# rio::export(Data_unif_long, "M2 long_format_2.xlsx")


# Export data for RMN in specified order
# rmn_order_ids <- paste0(c(1, 14, 13, 6, 4, 5, 10, 15, 2, 16, 42, 34, 24, 25, 26, 35, 39, 40, 32, 59, 53, 64, 46, 47, 62, 63, 48, 51, 60,
#                    8, 20, 12, 18, 9, 17, 19, 7, 3, 43, 41, 27, 23, 38, 22, 33, 28, 30, 21, 37, 49, 56, 58, 57, 52, 45, 55, 54), " ", "M2")
# 
# RMN_df <- Data_unif_wide[match(rmn_order_ids, Data_unif_wide$ID),]    # only 57 subjs for RMN
# 
# rio::export(RMN_df, "RMN ordered M2 wide_format_2.xlsx")

5.1 Read NFkB

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Read, Clean, Recode, Unite
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Read files
folder <- "C:/Users/Mihai/Desktop/R Notebooks/notebooks/M2-report"
file <- "lista_M2 cu rezultate NFkB.xlsx"

setwd(folder)

Data_nfkb_TR <-  readxl::read_xlsx(file.path(folder, file), 
                                   range = "A1:F30", sheet = "NFkB")
Data_nfkb_CTRL <-  readxl::read_xlsx(file.path(folder, file), 
                                     range = "H1:M31", sheet = "NFkB")

Data_nfkb_TR <- Data_nfkb_TR[-1, ]
colnames(Data_nfkb_TR)[c(1, 4, 6)] <- c("ID", "NFkB_Pre", "NFkB_Post")
Data_nfkb_TR <- 
  Data_nfkb_TR %>%
  dplyr::rename_all(~stringr::str_replace_all(., " ", "_")) %>%
  dplyr::select(-starts_with("ID_Parhon")) %>%
  dplyr::mutate(Cond = rep("TR", nrow(Data_nfkb_TR))) %>%
  dplyr::mutate(ID = stringr::str_remove(ID, "ID")) %>%
  dplyr::mutate(ID = stringr::str_remove(ID, " ")) %>%
  dplyr::mutate(ID = stringr::str_c(ID, " M2"))
  

Data_nfkb_CTRL <- Data_nfkb_CTRL[-1, ]
colnames(Data_nfkb_CTRL)[c(1, 4, 6)] <- c("ID", "NFkB_Pre", "NFkB_Post")
Data_nfkb_CTRL <- 
  Data_nfkb_CTRL %>%
  dplyr::rename_all(~stringr::str_replace_all(., " ", "_")) %>%
  dplyr::select(-starts_with("ID_Parhon")) %>%
  dplyr::mutate(Cond = rep("CTRL", nrow(Data_nfkb_CTRL))) %>%
  dplyr::mutate(ID = stringr::str_remove(ID, "ID")) %>%
  dplyr::mutate(ID = stringr::str_remove(ID, " ")) %>%
  dplyr::mutate(ID = stringr::str_c(ID, " M2"))

6 Unite with nfkb

# check
all.equal(colnames(Data_nfkb_TR), colnames(Data_nfkb_CTRL))

[1] TRUE

Data_unif_wide2 <- dplyr::left_join(Data_unif_wide, rbind(Data_nfkb_TR, Data_nfkb_CTRL), by = "ID")

# check
all.equal(Data_unif_wide2$Cond_Pre, Data_unif_wide2$Cond)

[1] "'is.NA' value mismatch: 2 in current 0 in target"

Data_unif_wide2 <- 
  Data_unif_wide2 %>%
  dplyr::select(-c(Nume_subiect, Cond))

7 NfKB

func_t_box(Data_nfkb_TR, ind = "ID", pre_var = "NFkB_Pre", post_var = "NFkB_Post")

#### NFkB_Pre NFkB_Post

Data_nfkb_TR %>%
  tidyr::pivot_longer(cols = c(NFkB_Pre, NFkB_Post), names_to = "PrePost", values_to = "NFkB") %>%
  rstatix::wilcox_test(NFkB ~ PrePost, paired = TRUE) %>%
  rstatix::add_significance()

func_t_box(Data_nfkb_CTRL, ind = "ID", pre_var = "NFkB_Pre", post_var = "NFkB_Post")

#### NFkB_Pre NFkB_Post

Data_nfkb_CTRL %>%
  tidyr::pivot_longer(cols = c(NFkB_Pre, NFkB_Post), names_to = "PrePost", values_to = "NFkB") %>%
  rstatix::wilcox_test(NFkB ~ PrePost, paired = TRUE) %>%
  rstatix::add_significance()

7.1 New Ideas: Correlations Number of memories & remembering

cat("PRE-POST Comparisons by Cond: TR")

PRE-POST Comparisons by Cond: TR

Data_unif_long %>%
  dplyr::filter(Cond == "TR") %>% 
  rstatix::wilcox_test(s1.count ~ PrePost, paired = TRUE) %>%
  rstatix::add_significance()

cat("## PRE-POST Comparisons by Cond: CTRL")

## PRE-POST Comparisons by Cond: CTRL

Data_unif_long %>%
  dplyr::filter(Cond == "CTRL") %>% 
  rstatix::wilcox_test(s1.count ~ PrePost, paired = TRUE) %>%
  rstatix::add_significance()

cat("## POST-vars Comparisons by Cond")

## POST-vars Comparisons by Cond

Data_unif_long %>%
  dplyr::filter(PrePost == "Post") %>%
  rstatix::wilcox_test(s1.count ~ Cond, paired = FALSE) %>%
  rstatix::add_significance()

Data_unif_long %>%
  dplyr::filter(PrePost == "Post") %>% 
    rstatix::t_test(s1.count ~ Cond, paired = FALSE) %>%
    rstatix::add_significance()

Data_unif_long %>%
  dplyr::filter(PrePost == "Post") %>%
  rstatix::wilcox_test(voluntary ~ Cond, paired = FALSE) %>%
  rstatix::add_significance()                   

Data_unif_long %>%
  dplyr::filter(PrePost == "Post") %>%
    rstatix::t_test(voluntary ~ Cond, paired = FALSE) %>%
    rstatix::add_significance()

cat("## Correlation - whole sample")

## Correlation - whole sample

Data_unif_long %>%
  dplyr::filter(PrePost == "Post") %>%
  dplyr::select(s1.count, spontaneous, voluntary, dreams) %>% 
  PerformanceAnalytics::chart.Correlation()

cat("## Correlation - TR")

## Correlation - TR

Data_unif_long %>%
  dplyr::filter(PrePost == "Post") %>%
  dplyr::filter(Cond == "TR") %>%
  dplyr::select(s1.count, spontaneous, voluntary, dreams) %>% 
  PerformanceAnalytics::chart.Correlation()

cat("## Correlation - CTRL")

## Correlation - CTRL

Data_unif_long %>%
  dplyr::filter(PrePost == "Post") %>%
  dplyr::filter(Cond == "CTRL") %>%
  dplyr::select(s1.count, spontaneous, voluntary, dreams) %>% 
  PerformanceAnalytics::chart.Correlation()

8 Data in wide format

Data_unif_wide2 %>%
    DT::datatable(
      extensions = 'Buttons',
      options = list(pageLength = 10,
                     scrollX='500px',
                     dom = 'Bfrtip',
                     buttons = c('excel', "csv")))

NA

9 Session Info

 

---

A work by Claudiu Papasteri

claudiu.papasteri@gmail.com