O1A Report

Oxytocin and Cortisol

1 Load data

## Read
Data <- rio::import("rezultate_oxitocina_pre_post_O1a_id28corectat.xlsx", sheet = "oxitocina+cortizol")
Data_PsySOP <- readRDS("Date_merged_Psy&SOP.RDS")

2 Processing data

## Clean
Data <-
  Data %>%
    dplyr::filter(rowSums(is.na(.)) < 8) %>%               # filter out rows (no more than 7 NA on row)
    dplyr::select(which(colMeans(is.na(.)) < 0.5))         # filter out columns with Catalina's statistics (no more than 1/2 NA on column values)
## Transform in oder to be consistent with past behavioral data 
oldnames = colnames(Data)
newnames = c("ID_pre","Cortizol_pre", "Oxitocina_pre", "ID_post","Cortizol_post", "Oxitocina_post", "Conditie")
Data <-
  Data %>% 
  dplyr::rename_at(vars(oldnames), ~ newnames) %>%
  dplyr::select(-ID_post) %>%
  dplyr::mutate(ID_pre = stringr::str_remove_all(ID_pre, c(" proba A|/proba A"))) %>%     # small inconsistency "/proba"
  tidyr::separate(ID_pre,  c("ID", "Ziua", "Nr_zi"), "\\s+") %>%                          # split on white space
  dplyr::mutate(Ziua = rep("zi", length.out = n())) %>%
  tidyr::unite("Zi", c("Ziua", "Nr_zi"), sep = "", remove = TRUE) %>%
  mutate(ID = as.numeric(str_extract(ID, "[^/]+"))) %>%                                   # [^/]+ matches 1 or more chars other than /
  mutate(ID = as.character(ID))                                                             # ID Psy&SOP are char, not numeric...for merge
## Melt -- not needed here
Data_melt <-
  Data %>% 
  gather(variable, value, -ID, -Zi, -Conditie) %>%
  tidyr::separate(variable,  c("variable", "PrePost"), "_") %>%
  spread(variable, value) %>%
  mutate(PrePost = factor(PrePost, levels = c("pre","post"))) %>%
  dplyr::arrange(ID)

3 Merge with SOP

## Merge with Psy&SOP data
Data_merged <- left_join(Data_PsySOP, Data, by = c("ID", "Conditie"))             # all good: ID 1,4,36,42 dont have Oxy&Cort 
varnottable <- c("Nume_Prenume", "NA_per_row",
                 "IOS_mama", "IOS_tata", "IOS_iubit", "IOS_prieten", "IOS_personalitate",
                 sprintf("Stais_pre_%01d", seq(1,20)), 
                 sprintf("Stais_post_%01d", seq(1,20)))
Data_merged <- 
  Data_merged %>%
  rename_at(vars(ends_with("Pre")), funs(gsub("Pre", "pre", .))) %>%              # _Pre and _Post -> tolower ... helps for automation
  rename_at(vars(ends_with("Post")), funs(gsub("Post", "post", .))) %>%
  select(-varnottable)

3.1 Table of Merged O.1.A (Psy&SOP&Oxy)

Data_merged %>%                              
    DT::datatable(                                  # excel downloadable  DT table
      extensions = 'Buttons',
      options = list(pageLength = 20,
                     scrollX='500px', 
                     dom = 'Bfrtip', 
                     buttons = c('excel', "csv")))

4 Define Functions

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Define Function for mining correlations
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Function for p-value significance -- both for func_ancova_multibox(), Get_Top_Relationships() and Correlations_With_One()
stars_signif <- function(pval) {
  stars = "ns"
  if(pval <= 0.001)
    stars = "***"
  if(pval > 0.001 & pval <= 0.01)
    stars = "**"
  if(pval > 0.01 & pval <= 0.05)
    stars = "*"
  if(pval > 0.05 & pval <= 0.1)
    stars = "."
  stars
}
## Function that returns correlations of all variables in descending order.
# Arg for threshold with default at .3 will keep only correlantions above .3 and below -.3. Also has threshhold for p-value. 
Get_Top_Relationships <- function(data_set, 
                                  correlation_abs_threshold=0.3,
                                  pvalue_threshold=0.05) {
  require(psych)
  require(dplyr)
  feature_names <- names(data_set)
  # strip var names to index for pair-wise identification
  names(data_set) <- seq(1:ncol(data_set))
  # calculate correlation and significance numbers
  cor_data_df <- psych::corr.test(data_set)
  # apply var names to correlation matrix over index
  rownames(cor_data_df$r) <- feature_names
  colnames(cor_data_df$r) <- feature_names
  # top cor and sig
  relationships_set <- cor_data_df$ci[,c('r','p')]
  # apply var names to data over index pairs
  relationships_set$feature_1 <- feature_names[as.numeric(sapply(strsplit(rownames(relationships_set), "-"), `[`, 1))]
  relationships_set$feature_2 <- feature_names[as.numeric(
    sapply(strsplit(rownames(relationships_set), "-"), `[`, 2))]
  relationships_set <- dplyr::select(relationships_set, feature_1, feature_2, r, p) %>% dplyr::rename(correlation = r, p.value = p)
  # return only the most insteresting relationships
  return(relationships_set %>%
         filter(abs(correlation) > correlation_abs_threshold &
                  p.value < pvalue_threshold) %>%
         arrange(p.value) %>%
         mutate(p.signif = sapply(p.value, function(x) stars_signif(x)))) %>%
         mutate(p.value = round(p.value, 3)) 
}
## Function for ploting correlation data frames resulting from Get_Top_Relationships and Correlations_With_One()
func_dotplot_cor <- function(df){                                        # https://www.r-pkg.org/pkg/ggpubr
  dotplotcor_scale_fill <- function(...){                                # Fix colors to signif factor levels even if missing
    ggplot2:::manual_scale(                                   
      'color', 
      values = setNames(
        c("darkgreen", "green3", "lawngreen", "yellow", "red"), 
        c("***", "**", "*", ".", "ns")), 
      ...
    )
  }                                           
  
  dtoplot_theme <- 
    ggpubr::theme_pubr() +
    theme(axis.text.y = element_text(size = 10))
  
  if(!"Variable" %in% colnames(df)){                                             # in oder to work for both Get_Top_Relationships and Correlations_With_One()
  df <- 
    df %>%                                            
      unite(cor_between, c("feature_1", "feature_2"), sep = " X ")               # unite 2 columns to x name from plot
  }else df <- df %>% dplyr::rename(cor_between = Variable)                       # change Variable to x name from plot
  
  df %>%
    ggpubr::ggdotchart(x = "cor_between", y = "correlation",
                       color = "p.signif",                                       # Color by sig
                       #   palette = c("#00AFBB", "#E7B800", "#FC4E07"),         # Custom color palette
                       sorting = "descending",                                   # Sort value in descending order
                       add = "segments",                                         # Add segments from y = 0 to dots
                       add.params = list(color = "lightgray", size = 2),         # Change segment color and size
                       group = "p.signif",                                       # Order by groups
                       dot.size = 8,                                             # Large dot size
                       xlab = "",
                       rotate = TRUE,                                            # Rotate vertically
                       label = round(.$correlation, 1),                          # Add mpg values as dot labels
                       font.label = list(color = "white", size = 9, 
                                         vjust = 0.5),                           # Adjust label parameters
                       ggtheme = dtoplot_theme) +                                # ggplot2 theme
    dotplotcor_scale_fill() +                                            # Fix colors to signif factor levels even if missing
    geom_hline(yintercept = 0, linetype = 2, color = "lightgray")
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Function for t test and boxplot
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
func_t_box <- function(df, ID, Conditie, pre_var, post_var){
  df_modif <-
    df %>%
    select(ID, Conditie, pre_var, post_var) %>% 
    tidyr::drop_na() %>%
    gather(pre_var, post_var, key = "PrePost", value = "value") %>% 
    mutate_at(vars(c(1, 2)), funs(as.factor)) %>% 
    mutate(PrePost = factor(PrePost, levels = c(pre_var, post_var))) 
  
  stat_comp <-
    df_modif %>% 
    group_by(Conditie) %>% 
    do(tidy(t.test(.$value ~ .$PrePost,
                   paired = TRUE,
                   data=.)))
  
  plot <- 
    ggpubr::ggpaired(df_modif, x = "PrePost", y = "value", id = ID, 
                     color = "PrePost", line.color = "gray", line.size = 0.4,
                     palette = c("#00AFBB", "#FC4E07"), legend = "none") +
    facet_wrap(~Conditie) +
    stat_summary(fun.data = mean_se,  colour = "darkred") +
    ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label.x = as.numeric(df_modif$PrePost) * 0.90, label.y = max(df_modif$value) * 1.15) + 
    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(plot)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

5 Analyses

5.1 Pre-Post differences on Oxytocin and Cortisol

func_t_box(Data, "ID", "Conditie", "Cortizol_pre", "Cortizol_post")  

5.1.0.1 Cortizol_pre Cortizol_post

func_t_box(Data, "ID", "Conditie", "Oxitocina_pre", "Oxitocina_post") 

5.1.0.2 Oxitocina_pre Oxitocina_post

5.2 Correlations between Difference Scores

## Automatic Post-Pre Difference Scores
var_prepost <- colnames(Data_merged[, grepl(".*pre|.*post", colnames(Data_merged))])    # find variables that have "pre" and "post"
var_prepost <-  gsub("_pre", "", var_prepost)                                           # delete pre and post from colname
var_prepost <-  gsub("_post", "", var_prepost)
var_prepost <- unique(var_prepost)                                                      # keep unique -- all good         
new_diffvar_name <- vector(mode="character", length = length(var_prepost))              # initialize vector for for()
for(i in seq_along(var_prepost)){
    var_name_pre <- paste0(var_prepost[i], "_pre")
    var_name_post <- paste0(var_prepost[i], "_post")
    new_diffvar_name[i] <- paste0(var_prepost[i], "_Diff")
    Data_merged[, new_diffvar_name[i]] <- Data_merged[, var_name_post] - Data_merged[, var_name_pre]
}
# Select variables
vars_cor <- new_diffvar_name[!new_diffvar_name %in% c("AllButtons_Diff", "AllButtons45_Diff", 
                                                      "Ones_Diff", "Ones45_Diff",
                                                      "Threes_Diff", "Threes45_Diff", 
                                                      "Twos_Diff", "Twos45_Diff")]       # exclude these
vars_cor <- c(vars_cor,
              c("Vas_rel_global", "Vas_rel_arousal", "CRQ_1", "CRQ_2", "CRQ_3", "CRQ_4", "CRQ_5", "CRQ_6"))
## Correlations
Get_Top_Relationships(Data_merged[, vars_cor])

Get_Top_Relationships(Data_merged[, vars_cor]) %>% func_dotplot_cor()

6 Conditioning on VAS global

6.1 Pre-Post differences on Oxytocin and Cortisol

Data_merged %>%
  filter(Vas_rel_global >= 7) %>%
  func_t_box(., "ID", "Conditie", "Cortizol_pre", "Cortizol_post") 

6.1.0.1 Cortizol_pre Cortizol_post

Data_merged %>%
  filter(Vas_rel_global >= 7) %>%
  func_t_box(., "ID", "Conditie", "Oxitocina_pre", "Oxitocina_post") 

6.1.0.2 Oxitocina_pre Oxitocina_post

7 Moderation Analysis

## Define Function
func_moderation <- function(Data, dep, mod, pred){
  moderation <- 
    Data %>%
      medmod::mod(., dep = dep, mod = mod, pred = pred,
                  ci = TRUE, estMethod = 'standard', test = TRUE, simpleSlopeEst = FALSE, simpleSlopePlot = TRUE)
  
  cat(paste("<b> Moderation: ", "Dep = ", dep, "Pred = ", pred, "Mod = ", mod, "</b>"))
  moderation$mod %>% 
    knitr::kable(caption = "Moderation", digits = 3) %>%
    print()
  moderation$simpleSlope$plot %>%
    print()
}
## Apply Function
func_moderation(Data = Data_merged, dep = "Oxitocina_post", mod = "Vas_rel_arousal", pred = "Oxitocina_pre")

Moderation: Dep = Oxitocina_post Pred = Oxitocina_pre Mod = Vas_rel_arousal

term	est	se	lower	upper	z	p
Oxitocina_pre	0.727	0.057	0.615	0.840	12.669	0.000
Vas_rel_arousal	0.000	0.001	-0.003	0.003	-0.171	0.864
Oxitocina_pre:Vas_rel_arousal	0.004	0.002	0.000	0.008	2.038	0.042

func_moderation(Data = Data_merged, dep = "Cortizol_post", mod = "Vas_rel_arousal", pred = "Cortizol_pre")

Moderation: Dep = Cortizol_post Pred = Cortizol_pre Mod = Vas_rel_arousal

term	est	se	lower	upper	z	p
Cortizol_pre	0.483	0.025	0.434	0.533	19.142	0.000
Vas_rel_arousal	0.011	0.005	0.002	0.021	2.382	0.017
Cortizol_pre:Vas_rel_arousal	0.003	0.001	0.000	0.005	2.237	0.025

# bla <-
#   Data_merged %>%
#   mutate(Conditie = as.numeric(as.factor(Conditie)))
# 
# psych::mediate(data = bla, Oxitocina_post ~ Oxitocina_pre + Conditie + Vas_rel_global)   # moderation with covariate and diagram

8 Mediation Analysis

9 Session Info

R version 3.5.2 (2018-12-20)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=Romanian_Romania.1250  LC_CTYPE=Romanian_Romania.1250    LC_MONETARY=Romanian_Romania.1250 LC_NUMERIC=C                     
[5] LC_TIME=Romanian_Romania.1250    

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

other attached packages:
 [1] bindrcpp_0.2.2     rio_0.5.16         plyr_1.8.4         summarytools_0.9.3 DT_0.5             ggpubr_0.2         magrittr_1.5      
 [8] broom_0.5.1        papaja_0.1.0.9842  psych_1.8.10       forcats_0.3.0      stringr_1.3.1      dplyr_0.7.8        purrr_0.2.5       
[15] readr_1.3.0        tidyr_0.8.2        tibble_1.4.2       ggplot2_3.1.0      tidyverse_1.2.1    pacman_0.5.0      

loaded via a namespace (and not attached):
 [1] nlme_3.1-137       bitops_1.0-6       matrixStats_0.54.0 lubridate_1.7.4    httr_1.4.0         tools_3.5.2        backports_1.1.3   
 [8] R6_2.4.0           lazyeval_0.2.1     colorspace_1.3-2   withr_2.1.2        tidyselect_0.2.5   mnormt_1.5-5       curl_3.3          
[15] compiler_3.5.2     cli_1.0.1          rvest_0.3.2        xml2_1.2.0         labeling_0.3       scales_1.0.0       checkmate_1.8.5   
[22] pbivnorm_0.6.0     digest_0.6.18      foreign_0.8-71     rmarkdown_1.11     base64enc_0.1-3    pkgconfig_2.0.2    htmltools_0.3.6   
[29] highr_0.7          htmlwidgets_1.3    rlang_0.3.0.1      readxl_1.1.0       rstudioapi_0.8     pryr_0.1.4         jmvcore_0.9.5.2   
[36] shiny_1.2.0        bindr_0.1.1        generics_0.0.2     jsonlite_1.6       crosstalk_1.0.0    zip_1.0.0          RCurl_1.95-4.11   
[43] rapportools_1.0    Rcpp_1.0.1         munsell_0.5.0      stringi_1.2.4      yaml_2.2.0         MASS_7.3-51.1      lavaan_0.6-3      
[50] grid_3.5.2         parallel_3.5.2     promises_1.0.1     crayon_1.3.4       lattice_0.20-38    haven_2.1.0        pander_0.6.3      
[57] hms_0.4.2          magick_2.0         knitr_1.21         pillar_1.3.1       tcltk_3.5.2        rjson_0.2.20       ggsignif_0.4.0    
[64] stats4_3.5.2       codetools_0.2-15   glue_1.3.1         evaluate_0.12      data.table_1.12.2  modelr_0.1.2       httpuv_1.4.5      
[71] cellranger_1.1.0   gtable_0.2.0       assertthat_0.2.1   xfun_0.4           openxlsx_4.1.0     mime_0.6           xtable_1.8-3      
[78] later_0.7.5        medmod_1.0.0      

 

---

A work by Claudiu Papasteri

claudiu.papasteri@gmail.com