General Plots for M.1. (Autobiographical Memories)

Initial Dataset

1 Read

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Read
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Read files
file <- "Date Complete M1 v.13 siPPGGSRamilaza.sav"

# setwd(folder)
Data <- rio::import(file.path(folder, file))

2 Make data frame

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

2.1 Exclude Protocol 8 (mother)

Data <- 
  Data %>%
  filter(P != 8)

3 Define Functions

## Func t test si boxplot simplu
func_t_box <- function(df, ind, pre_var, post_var, facet = FALSE, xlab = ""){  
  if(facet){
    facet <- "Protocol"
  }else{
    facet <- NULL
  }
  
  df_modif <-
    df %>%
    select(ind, P, 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))) 
  
  if(!is.null(facet)){
    df_modif <-
      df_modif %>%
      group_by(P) %>%
      mutate(Protocol = paste0("Protocol = ", P, ", n = ", n()))
  }
  
  stat_comp <-
    df_modif %>% 
    do(tidy(t.test(.$value ~ .$PrePost,
                   paired = TRUE,
                   data=.)))
  
  plot <- 
    ggpubr::ggpaired(df_modif, x = "PrePost", y = "value", id = ind, 
                     color = "PrePost", line.color = "gray", line.size = 0.4,
                     palette = c("#00AFBB", "#FC4E07"), legend = "none",
                     facet.by = facet, ncol = 3, 
                     xlab = xlab) +
    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.4, label.y = max(df_modif$value)+1) + 
    ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label = "p.signif", comparisons = list(c(pre_var, post_var)))
  
  print(stat_comp)
  cat("\n")                      
  print(plot)
  cat("\n")
  plot.new()                     # Need this workaround for interleaving tables and plots in R Markdown, within loop
  dev.off()
}

heat_cor_plotly <- function(df, x_vars = NULL, y_vars = NULL, low_color = "cyan",  high_color = "red",  ...){   
  # inherit type = c("pearson","spearman") from Hmisc::rcorr() 
  library(ggplot2)
  library(plotly)
  library(reshape2)
  library(Hmisc)
  
  # use all numeric columns only, print message if non-numeric are found
  numeric_cols <- unlist(lapply(df, is.numeric))
  if(!all(numeric_cols)) message("Warning: Non-numeric columns were excluded!")
  df <- df[, numeric_cols]
  
  df_mat <- as.matrix(df)
  rt <- Hmisc::rcorr(df_mat, ...)
  
  # extract correlations, p-values and merge into another dataframe
  mtlr <- reshape2::melt(rt$r, value.name = "Correlation")
  mtlp <- reshape2::melt(rt$P, value.name = "P-Value")
  
  mtl <- merge(mtlr, mtlp)
  
  # give possibility to prune the correlation matrix
  if(!is.null(x_vars)){
    mtl <- mtl[(mtl$Var1 %in% x_vars), ]
  }
  if(!is.null(x_vars)){
    mtl <- mtl[(mtl$Var2 %in% y_vars), ]
  }
  
  # want to avoid scientific notetion, but this doesnt work as numeric
  # mtl$Correlation <- as.numeric(format(mtl$Correlation, digits = 4, scientific = FALSE))  # doesnt work
  # mtl$`P-Value` <- as.numeric(format(mtl$`P-Value`, digits = 4, scientific = FALSE)) 
  options(scipen = 999)
  mtl$Correlation <- round(mtl$Correlation, 3)
  mtl$`P-Value` <- round(mtl$`P-Value`, 3)

  gx <-
    ggplot2::ggplot(mtl, 
           aes(Var1, Var2, 
               fill = Correlation,  
               text = paste("P-val = ", `P-Value`))) +
    ggplot2::geom_tile() + 
    ggplot2::scale_fill_gradient(low = low_color,  high = high_color, limits = c(-1, 1), breaks = c(-1, -.5, 0, .5, 1)) +
    ggplot2::theme_minimal() +
    {if(any(nchar(names(df)) > 6)) ggplot2::theme(axis.text.x = element_text(angle = 90, hjust = 1))}  # vertical x axis labels if lenghty
  plotly::ggplotly(gx)  
}

4 Plot Age

## Dodged Bar plot of Age and Gender
Data  %>%
  mutate(Varta_categ = cut(Varsta, 
                           breaks=c(-Inf, 25, 30, 35, 40, 45, 50, 55, 60, Inf), 
                           labels=c("<25","25-29","30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60>"), 
                           right = FALSE)) %>%  
  mutate(Varsta = as.factor(Varsta),
         Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "1", "masculin" = "2")) %>%
  dplyr::count(Varta_categ, Gen, .drop = FALSE) %>%         # Group by, then count number in each group (dont drop 0 counts)
  mutate(pct = prop.table(n)) %>%                           # Calculate percent within each var
    ggplot(aes(x = Varta_categ, y = pct, fill = Gen, label = scales::percent(pct))) + 
      geom_col(position = position_dodge(preserve = "single"), stat = "identity",) +    # Don't drop zero count
      geom_text(position = position_dodge(width = .9),      # move to center of bars
                vjust = -0.5,                               # nudge above top of bar
                size = 3) + 
      scale_y_continuous(labels = scales::percent) +
      ggtitle("") +
      xlab("Varsta") + ylab("Percentage %") + 
      guides(fill = guide_legend(title = "Gen", ncol = 1)) + 
      scale_fill_grey(start = 0.8, end = 0.2, na.value = "red", aesthetics = "fill") +
      theme(legend.position = "right", legend.direction = "vertical", 
            legend.justification = c(0, 1), panel.border = element_rect(fill = NA, colour = "black"))

4.1 By Protocol

## Dodged Bar plot of Age and Gender by Protocol
Data  %>%
  mutate(Varta_categ = cut(Varsta, 
                           breaks=c(-Inf, 25, 30, 35, 40, 45, 50, 55, 60, Inf), 
                           labels=c("<25","25-29","30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60>"), 
                           right = FALSE)) %>%  
  mutate(Varsta = as.factor(Varsta),
         Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "1", "masculin" = "2")) %>%
  group_by(P) %>%
  dplyr::count(Varta_categ, Gen, .drop = FALSE) %>%         # Group by, then count number in each group (dont drop 0 counts)
  mutate(pct = prop.table(n)) %>%                           # Calculate percent within each var
    ggplot(aes(x = Varta_categ, y = pct, fill = Gen, label = scales::percent(pct))) +
      facet_wrap(~P, scales = "free", ncol = 1) +
      geom_col(position = position_dodge(preserve = "single"), stat = "identity",) +    # Don't drop zero count
      geom_text(position = position_dodge(width = .9),      # move to center of bars
                vjust = -0.5,                               # nudge above top of bar
                size = 3) + 
      scale_y_continuous(labels = scales::percent) +
      ggtitle("") +
      xlab("Varsta") + ylab("Percentage %") + 
      guides(fill = guide_legend(title = "Gen", ncol = 1)) + 
      scale_fill_grey(start = 0.8, end = 0.2, na.value = "red", aesthetics = "fill") +
      theme(legend.position = "right", legend.direction = "vertical", 
            legend.justification = c(0, 1), panel.border = element_rect(fill = NA, colour = "black"))

## Pie chart
Data  %>%
  mutate(Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "1", "masculin" = "2")) %>%
  group_by(Gen) %>%
  dplyr::summarise(counts = n()) %>%
  mutate(prop = round(counts*100/sum(counts), 1),
         lab.ypos = cumsum(prop) - .5*prop,
         Percent = paste0(prop, " %")) %>% 
  ggpubr::ggpie(x = "prop", label = "Percent",
                fill = "Gen", color = "white", 
                lab.pos = "in", lab.font = list(color = "white"),
                palette = "grey")

5 Analyses

5.1 Simple before-after analyses with t test

## Simple before-after analyses with t test
cat("#### VAS Stress")

5.1.0.1 VAS Stress

func_t_box(Data, ind = "ID", "Stres_pre", "Stres_post", facet = FALSE) 

NANA

null device 1

## Simple before-after analyses with t test
cat("#### VAS Stress")

5.1.0.2 VAS Stress

func_t_box(Data, ind = "ID", "Stres_pre", "Stres_post", facet = TRUE) 

NANA

null device 1

5.2 Correlations: Anotimpuri - Calitate Amintiri (without P6, P7)

dateplot1 <- Data[, c("P", "Primavara", "Vara", "Toamna", "Iarna", "Media_s1", "Media_s2", "Media_s3",  "SocDih_Part",  "SocDih_FamN",  "SocDih_FamInd",  "SocDih_Priet",  "SocDih_Amici",  "SocDih_Necun",  "SocDih_Antag",  "SocDih_TotAprop",  "SocDih_TotNeaprop", "STAI_T")] 
names(dateplot1) <- c("P", "Primavara", "Vara", "Toamna", "Iarna", "S1- Valenta", "S2 - Vividness", "S3 - Relevanta",  "Partener",  "Familie nucleu",  "Familie extinsa",  "Prieteni",  "Amici",  "Necunoscuti",  "Antagonisti",  "Toti Apropiatii",  "Toti Neapropiatii", "STAI_T")
dateplot1 <- subset(dateplot1, P!=6 & P!=7)

COR <- Hmisc::rcorr(as.matrix(dateplot1[,-1]))   
M <- COR$r
P_MAT <- COR$P
corrplot::corrplot(M, method = "number", type = "upper", p.mat = P_MAT, sig.level = 0.05, insig = "blank", tl.col = "black", tl.cex = .9, tl.srt = 45)  

heat_cor_plotly(dateplot1[,-1])

5.3 Correlations: Personality - Qualities of Memories (without P6, P7)

dateplot2 <- Data[, c(24, 40, 56, 87:121, 126)] 
names(dateplot2)[1:3] <- c("S1- Valenta", "S2 - Vividness", "S3 - Relevanta")

COR <- Hmisc::rcorr(as.matrix(dateplot2))   
M <- COR$r
P_MAT <- COR$P
corrplot::corrplot(M, type = "upper", p.mat = P_MAT, sig.level = 0.05, insig = "blank", tl.col = "black", tl.cex = .7, cl.pos = "b", tl.srt = 45)

heat_cor_plotly(dateplot2, x_vars = names(dateplot2)[1:3], y_vars = names(dateplot2)[-(1:3)])

5.4 Correlations: Social - Personality

dateplot3 <- Data[, c(131:139, 87:121)]
names(dateplot3)[1:9] <- c("Partener",  "Familie nucleu",  "Familie extinsa",  "Prieteni",  "Amici",  "Necunoscuti",  "Antagonisti",  "Toti Apropiatii",  "Toti Neapropiatii")

COR <- Hmisc::rcorr(as.matrix(dateplot3))   
M <- COR$r
P_MAT <- COR$P
corrplot::corrplot(M, type = "upper", p.mat = P_MAT, sig.level = 0.05, insig = "blank", tl.col = "black", tl.cex = .7, cl.pos = "b", tl.srt = 45)

heat_cor_plotly(dateplot3, x_vars = names(dateplot3)[1:9], y_vars = names(dateplot3)[-(1:9)])

6 Social

# names(Data[str_detect(colnames(Data), fixed("SocDih", ignore_case=TRUE))])
Data_soc <-
  Data %>%
  dplyr::select("ID", "P", starts_with("SocDih")) %>%
  haven::zap_formats(.) %>%                # to not get warning  
  haven::zap_labels(.) %>%                 # "attributes are not identical across measure variables"
  haven::zap_widths(.) %>%                 # on gather() because of SPSS
  dplyr::rename_all(list(~stringr::str_replace(., "SocDih_", ""))) %>% 
  gather(key = "Variable", value = "Value", -ID, -P) 
  
# Create a custom color scale for all ASCQ graphs
library(RColorBrewer)
myColors <- brewer.pal(10,"Set1")
names(myColors) <- levels(Data_soc$Variable)
colScale <- scale_colour_manual(name = "Variable", values = myColors)

# Plot
ggpubr::ggviolin(data = Data_soc, x = "Variable", y = "Value", fill = "Variable",
      add = "boxplot", add.params = list(fill = "white"),
      xlab = "", legend = "none") +
  colScale +                                                   # color scale here keep consistency of color with factor level
  stat_summary(fun.data = mean_se,  colour = "darkred") +         
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Plot faceted by Protocol
p <- 
  ggpubr::ggviolin(data = Data_soc, x = "Variable", y = "Value", fill = "Variable",
        add = "boxplot", add.params = list(fill = "white"),
        xlab = "", legend = "none") +
         
    colScale +                                                   # color scale here keep consistency of color with factor level
    stat_summary(fun.data = mean_se,  colour = "darkred") +         
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

ggpubr::facet(p, facet.by = "P", ncol = 1, scales = "free_x")

6.1 Social - this doesnt make sense

# PerformanceAnalytics::chart.Correlation(Data[, c(8, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140)])

# Example of Simpsons Paradox
coplot(DifStres ~ SocDih_Amici | Media_s1,
       data = Data,
       rows = 1,
     panel = function(x, y, ...) {
          panel.smooth(x, y, span = .8, iter = 5,...)
          abline(lm(y ~ x), col = "blue")})

7 Varsta Amint - P1,P2,P3

Data_P1P2P3 <- 
  Data %>%
  filter(P %in% c("1", "2", "3")) %>%
  mutate(Med_amintvarsta = as.numeric(as.character(Med_amintvarsta)),
         Dif_Med_amintvarsta = Varsta - Med_amintvarsta)

PerformanceAnalytics::chart.Correlation(Data_P1P2P3[, c(8, 85, 122, 148)])

coplot(DifStres ~ Dif_Med_amintvarsta | Media_s1,
       data = Data_P1P2P3,
       rows = 1,
     panel = function(x, y, ...) {
          panel.smooth(x, y, span = .8, iter = 5,...)
          abline(lm(y ~ x), col = "blue")})

8 Protocol 3 - Social

Data_P3 <- 
  Data %>%
  filter(P == "3")

PerformanceAnalytics::chart.Correlation(Data_P3[, c(8, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140)])

# coplot(SocDih_Amici ~ DifStres | Media_s1, 
#        data = Data,
#        columns = 3,
#      panel = function(x, y, ...) {
#           panel.smooth(x, y, span = .8, iter = 5,...)
#           abline(lm(y ~ x), col = "blue") } )

9 Protocol 3 - Varsta Amint

Data_P3 <- 
  Data %>%
  filter(P == "3") %>%
  mutate(Med_amintvarsta = as.numeric(as.character(Med_amintvarsta)),
         Dif_Med_amintvarsta = Varsta - Med_amintvarsta)

PerformanceAnalytics::chart.Correlation(Data_P3[, c(8, 85, 122, 148)])

10 Session Info

R version 3.6.1 (2019-07-05)
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 LC_NUMERIC=C                     
[5] LC_TIME=Romanian_Romania.1250    

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

other attached packages:
 [1] RColorBrewer_1.1-2         Hmisc_4.1-1                Formula_1.2-3              survival_2.44-1.1          lattice_0.20-38           
 [6] reshape2_1.4.3             plotly_4.9.0               rio_0.5.16                 scales_1.0.0               ggpubr_0.2                
[11] magrittr_1.5               tadaatoolbox_0.16.1        summarytools_0.8.8         rstatix_0.2.0              broom_0.5.2               
[16] PerformanceAnalytics_1.5.2 xts_0.11-2                 zoo_1.8-4                  psych_1.8.12               plyr_1.8.4                
[21] forcats_0.4.0              stringr_1.4.0              dplyr_0.8.3                purrr_0.3.2                readr_1.3.1               
[26] tidyr_1.0.0                tibble_2.1.3               ggplot2_3.2.1              tidyverse_1.2.1            papaja_0.1.0.9842         
[31] pacman_0.5.1              

loaded via a namespace (and not attached):
 [1] colorspace_1.4-1    ggsignif_0.4.0      pryr_0.1.4          ellipsis_0.3.0      htmlTable_1.12      base64enc_0.1-3     rstudioapi_0.8     
 [8] DT_0.5              mvtnorm_1.0-11      lubridate_1.7.4     xml2_1.2.0          codetools_0.2-16    splines_3.6.1       mnormt_1.5-5       
[15] knitr_1.25          zeallot_0.1.0       pixiedust_0.8.6     jsonlite_1.6        cluster_2.1.0       shiny_1.2.0         compiler_3.6.1     
[22] httr_1.4.0          backports_1.1.4     assertthat_0.2.1    Matrix_1.2-17       lazyeval_0.2.2      cli_1.1.0           later_0.7.5        
[29] acepack_1.4.1       htmltools_0.3.6     tools_3.6.1         gtable_0.3.0        glue_1.3.1          Rcpp_1.0.2          carData_3.0-2      
[36] cellranger_1.1.0    vctrs_0.2.0         nlme_3.1-140        crosstalk_1.0.0     xfun_0.9            openxlsx_4.1.0      rvest_0.3.2        
[43] mime_0.7            lifecycle_0.1.0     MASS_7.3-51.4       hms_0.5.1           promises_1.0.1      parallel_3.6.1      expm_0.999-3       
[50] pwr_1.2-2           yaml_2.2.0          curl_3.2            gridExtra_2.3       pander_0.6.3        rpart_4.1-15        latticeExtra_0.6-28
[57] stringi_1.4.3       corrplot_0.84       nortest_1.0-4       checkmate_1.8.5     boot_1.3-22         zip_1.0.0           rlang_0.4.0        
[64] pkgconfig_2.0.3     matrixStats_0.54.0  bitops_1.0-6        rapportools_1.0     htmlwidgets_1.3     labeling_0.3        tidyselect_0.2.5   
[71] R6_2.4.0            DescTools_0.99.29   generics_0.0.2      pillar_1.4.2        haven_2.1.1         foreign_0.8-71      withr_2.1.2        
[78] nnet_7.3-12         abind_1.4-5         RCurl_1.95-4.11     modelr_0.1.5        crayon_1.3.4        car_3.0-2           viridis_0.5.1      
[85] grid_3.6.1          readxl_1.1.0        data.table_1.11.8   digest_0.6.21       xtable_1.8-4        httpuv_1.4.5        munsell_0.5.0      
[92] viridisLite_0.3.0   quadprog_1.5-5     

 

---

A work by Claudiu Papasteri

claudiu.papasteri@gmail.com