PRE_anxiety_state_MEAN   = scale(PRE_anxiety_state_MEAN),
PRE_anxiety_trait_MEAN   = scale(PRE_anxiety_trait_MEAN),
PRE_we_PROD              = scale(PRE_we_PROD),
PRE_dfi_PROD             = scale(PRE_dfi_PROD),
PRE_anxiety_state_DIFF   = scale(PRE_anxiety_state_DIFF),
PRE_anxiety_trait_DIFF   = scale(PRE_anxiety_trait_DIFF),
sam_arousal_MEAN         = scale(sam_arousal_MEAN),
sam_arousal_DIFF         = scale(sam_arousal_DIFF),
sam_valence_MEAN         = scale(sam_valence_MEAN),
sam_valence_DIFF         = scale(sam_valence_DIFF),
ecg_IBI_MEAN             = scale(ecg_IBI_MEAN),
rsp_Rate_MEAN            = scale(rsp_Rate_MEAN),
eda_Amplitude_MEAN       = scale(eda_Amplitude_MEAN),
ecg_WTC_LF               = scale(ecg_WTC_LF),
ecg_WTC_HF               = scale(ecg_WTC_HF),
rsp_WTC                  = scale(rsp_WTC),
eda_WTC                  = scale(eda_WTC),
ecg_rC                   = scale(ecg_rC),
ecg_rCC                  = scale(ecg_rCC),
rsp_rC                   = scale(rsp_rC),
rsp_rCC                  = scale(rsp_rCC),
eda_rC                   = scale(eda_rC),
eda_rCC                  = scale(eda_rCC))
# LOAD PACKAGES ####
library("tidyverse")
library("rlang")
library("MASS")
library("rstatix")
library("psych")
library("ggthemes")
library("ggpubr")
library("ggpp")
library("reshape2")
library("Hmisc")
library("data.table")
library("lme4")
library("nlme")
library("ordinal")
library("lmerTest")
library('buildmer')
library("interactions")
library('emmeans')
library('optimx')
library("sjPlot")
library("jtools")
library("performance")
library("reghelper")
library("ltm")
library("corrplot")
library("car")
library("coin")
# ADDITIONAL FUNCTIONS ####
# TO ADJUST ERROR BARS FOR REPEATED MEASURES
rmMeanAdjustLong <- function(dataframe, var_id, var_fact, var_value)
{
id         <- dataframe %>%  dplyr::pull(var_id)
df_values     <- dataframe %>%
pivot_wider(id_cols= all_of(var_id),
names_from = var_fact,
values_from = var_value)
df_id         <- df_values %>%  dplyr::select(Dyad)
df_values     <-  df_values %>%  dplyr::select(!Dyad)
pMean         <- rowSums(df_values)/length(df_values)
grandMean     <- mean(colMeans(df_values, na.rm  = TRUE), na.rm  = TRUE)
adj           <- grandMean - pMean
for(i in colnames(df_values)){
df_values[[i]]  <- df_values[[i]] + adj
}
output        <- df_values %>%
pivot_longer(everything(),
names_to = var_fact,
values_to = var_value) %>%
add_column(id, .before = var_fact) %>%
dplyr::rename(Dyad = id)
return(output)
}
# TO MERGE ALL FILES INTO A DATA FRAME
unite <- function(path, pattern, ignore_cols)
{
filesList <- list.files(path, pattern, full.names = TRUE)
for (file in filesList){
# if the merged dataset doesn't exist, create it
if (!exists("output")){
output <- read.table(file,
header=TRUE,
sep=";",
quote = "",
encoding = "latin1")
next
}else if (exists("output")){
temp_dataframe <- read.table(file, header=TRUE, sep=";", quote = "", encoding = "latin1")
output         <- cbind(output, temp_dataframe[-c(ignore_cols)])
rm(temp_dataframe)
}
}
return(output)
}
# LOAD GENERAL DATA FRAME ####
df_wide <- read.csv(paste0(path, 'data_table.csv'),
sep = ";", header = TRUE, encoding = "latin1")
# COMPUTE DYAD-LEVEL INDICES ####
# DYAD-LEVEL DIFFERENCES
df_wide_dyad              <- df_wide[seq(1, nrow(df_wide), 2), ]
Age_DIFF_tmp               = abs(diff(df_wide$Age))
PRE_anxiety_trait_DIFF_tmp = abs(diff(df_wide$PRE_anxiety_trait))
PRE_anxiety_state_DIFF_tmp = abs(diff(df_wide$PRE_anxiety_state))
POS_sam_arousal_DIFF_tmp   = abs(diff(df_wide$POS_sam_arousal))
NEG_sam_arousal_DIFF_tmp   = abs(diff(df_wide$NEG_sam_arousal))
NEU_sam_arousal_DIFF_tmp   = abs(diff(df_wide$NEU_sam_arousal))
POS_sam_valence_DIFF_tmp   = abs(diff(df_wide$POS_sam_valence))
NEG_sam_valence_DIFF_tmp   = abs(diff(df_wide$NEG_sam_valence))
NEU_sam_valence_DIFF_tmp   = abs(diff(df_wide$NEU_sam_valence))
df_wide_dyad               <- df_wide_dyad %>% mutate(
Age_DIFF                 = Age_DIFF_tmp[seq(1, length(Age_DIFF_tmp), 2)],
PRE_anxiety_trait_DIFF   = PRE_anxiety_trait_DIFF_tmp[seq(1, length(PRE_anxiety_trait_DIFF_tmp), 2)],
PRE_anxiety_state_DIFF   = PRE_anxiety_state_DIFF_tmp[seq(1, length(PRE_anxiety_state_DIFF_tmp), 2)],
POS_sam_arousal_DIFF     = POS_sam_arousal_DIFF_tmp[seq(1, length(POS_sam_arousal_DIFF_tmp), 2)],
NEG_sam_arousal_DIFF     = NEG_sam_arousal_DIFF_tmp[seq(1, length(NEG_sam_arousal_DIFF_tmp), 2)],
NEU_sam_arousal_DIFF     = NEU_sam_arousal_DIFF_tmp[seq(1, length(NEU_sam_arousal_DIFF_tmp), 2)],
POS_sam_valence_DIFF     = POS_sam_valence_DIFF_tmp[seq(1, length(POS_sam_valence_DIFF_tmp), 2)],
NEG_sam_valence_DIFF     = NEG_sam_valence_DIFF_tmp[seq(1, length(NEG_sam_valence_DIFF_tmp), 2)],
NEU_sam_valence_DIFF     = NEU_sam_valence_DIFF_tmp[seq(1, length(NEU_sam_valence_DIFF_tmp), 2)])
# DYAD-LEVEL MEANS
agg_tbl <- df_wide %>% group_by(Dyad) %>%
summarise(across(everything(), mean),
.groups = 'drop') %>%
dplyr::select(-c(Dyad, Participant, AttentionType, Gender, Stimuli_Sequence,ecg_clean, rsp_clean, eda_clean)) %>%
rename_with( ~ paste0(.x, "_MEAN"))
df_wide_dyad <- cbind(df_wide_dyad, agg_tbl)
rm(agg_tbl)
# DYAD-LEVEL PRODUCTS
j <- 0
for(i in seq(1,nrow(df_wide), 2)){
j <- j+1
# PRE_DFI
ppt1 <- as.numeric(df_wide$PRE_dfi[i])
ppt2 <- as.numeric(df_wide$PRE_dfi[i+1])
df_wide_dyad$PRE_dfi_PROD[j] <- ppt1*ppt2
# PRE_WE
ppt1 <- as.numeric(df_wide$PRE_we[i])
ppt2 <- as.numeric(df_wide$PRE_we[i+1])
df_wide_dyad$PRE_we_PROD[j] <- ppt1*ppt2
# POST_DFI
ppt1 <- as.numeric(df_wide$POST_dfi[i])
ppt2 <- as.numeric(df_wide$POST_dfi[i+1])
df_wide_dyad$POST_dfi_PROD[j] <- ppt1*ppt2
# POST_WE
ppt1 <- as.numeric(df_wide$POST_we[i])
ppt2 <- as.numeric(df_wide$POST_we[i+1])
df_wide_dyad$POST_we_PROD[j] <- ppt1*ppt2
# DELTA_DFI
df_wide_dyad$DELTA_dfi_PROD[j] <- df_wide_dyad$POST_dfi_PROD[j] -
df_wide_dyad$PRE_dfi_PROD[j]
# DELTA_WE
df_wide_dyad$DELTA_we_PROD[j] <- df_wide_dyad$POST_we_PROD[j] -
df_wide_dyad$PRE_we_PROD[j]
# CONNECTEDNESS
# Fixation
ppt1 <- as.numeric(df_wide$fix_connectedness[i])
ppt2 <- as.numeric(df_wide$fix_connectedness[i+1])
df_wide_dyad$fix_connectedness_PROD[j] <- ppt1*ppt2
# VIDEO POSITIVE
ppt1 <- as.numeric(df_wide$POS_connectedness[i])
ppt2 <- as.numeric(df_wide$POS_connectedness[i+1])
df_wide_dyad$POS_connectedness_PROD[j] <- ppt1*ppt2
# VIDEO NEGATIVE
ppt1 <- as.numeric(df_wide$NEG_connectedness[i])
ppt2 <- as.numeric(df_wide$NEG_connectedness[i+1])
df_wide_dyad$NEG_connectedness_PROD[j] <- ppt1*ppt2
# VIDEO NEUTRAL
ppt1 <- as.numeric(df_wide$NEU_connectedness[i])
ppt2 <- as.numeric(df_wide$NEU_connectedness[i+1])
df_wide_dyad$NEU_connectedness_PROD[j] <- ppt1*ppt2
}
# INFORMATION ABOUT PHYSIOLOGICAL ARTEFACTS
j = 0
for(i in seq(1,nrow(df_wide), 2)){
j = j+1
df_wide_dyad$ecg_clean[j] <- if_else(df_wide$ecg_clean[i] == 'y' && df_wide$ecg_clean[i+1] == 'y', 'y', 'n')
df_wide_dyad$rsp_clean[j] <- if_else(df_wide$rsp_clean[i] == 'y' && df_wide$rsp_clean[i+1] == 'y', 'y', 'n')
df_wide_dyad$eda_clean[j] <- if_else(df_wide$eda_clean[i] == 'y' && df_wide$eda_clean[i+1] == 'y', 'y', 'n')
}
# CREATE A LONGITUDINAL DATA FRAME FOR REGRESSION ANALYSES #####
# SELECT THE VARIABLES OF INTEREST
# Non-repeated measures and factors
myVarOne <- c("Dyad",
"AttentionType",
"Gender",
"Stimuli_Sequence",
"ecg_clean",
"rsp_clean",
"eda_clean",
"Age_MEAN",
"Age_DIFF",
"PRE_anxiety_trait_MEAN",
"PRE_anxiety_state_MEAN",
"PRE_anxiety_trait_DIFF",
"PRE_anxiety_state_DIFF",
"PRE_we_MEAN",
"PRE_dfi_MEAN",
"PRE_we_PROD",
"PRE_dfi_PROD",
"POST_we_MEAN",
"POST_dfi_MEAN",
"DELTA_we_MEAN",
"DELTA_we_PROD",
"DELTA_dfi_MEAN",
"DELTA_dfi_PROD")
# Repeated measures
myVarRep  <- c("sam_valence_MEAN",
"sam_valence_DIFF",
"sam_arousal_MEAN",
"sam_arousal_DIFF",
"ecg_IBI_MEAN",
"rsp_Rate_MEAN",
"eda_Amplitude_MEAN",
"eda_PeakNumber_MEAN",
"connectedness_MEAN",
"connectedness_PROD",
"exposure_MEAN",
"WTC",
"rC",
"rCC",
"PLV",
"WTC",
"CRQA")
# FIND THE CORRESPONDING COLUMNS IN THE INPUT DATA FRAME
myVarList <- c(myVarOne, myVarRep)
myColList <- c()
for(myVar in 1:length(myVarList)){
myCol <- colnames(df_wide_dyad[grepl(myVarList[myVar], names(df_wide_dyad))])
myColList <- unique(c(myColList, myCol))
}
# Exclude the baseline and average index of collective effervescence
myColList <- myColList[!grepl("FIX", myColList)]
myColList <- myColList[!grepl(".1", myColList)]
# SUBSET THE INPUT DATA FRAME ACCORDINGLY
df_regression  <- subset(df_wide_dyad, select = myColList)
# TRANSFORM FROM WIDE TO LONGITUDINAL FORMAT
myColListRep  <- setdiff(names(df_regression), myVarOne)
df_regression <- df_regression %>% pivot_longer(cols = all_of(myColListRep)) %>%
mutate(VideoValence = str_extract(name, "[POSNEGNEU]_"),
suffix = str_sub(str_extract(name, "_\\w+"),2)) %>%
dplyr::select(-name) %>%
pivot_wider(names_from =  suffix, values_from = value)
# RENAME THE CONDITION FACTOR
df_regression$VideoValence[grepl("S_",df_regression$VideoValence)] <- 'positive'
df_regression$VideoValence[grepl("G_",df_regression$VideoValence)] <- 'negative'
df_regression$VideoValence[grepl("U_",df_regression$VideoValence)] <- 'neutral'
# ADD A COLUMN TO SPECIFY THE ORDER OF VIDEO PRESENTATION
for(i in 1:nrow(df_regression)){
if(df_regression$VideoValence[i] == 'positive'){
if(df_regression$Stimuli_Sequence[i] == 'PosNegNeu'){
df_regression$VideoIndex[i] = 'first'
}else if(df_regression$Stimuli_Sequence[i] == 'PosNeuNeg'){
df_regression$VideoIndex[i] = 'first'
}else if(df_regression$Stimuli_Sequence[i] == 'NegPosNeu'){
df_regression$VideoIndex[i] = 'second'
}else if(df_regression$Stimuli_Sequence[i] == 'NegNeuPos'){
df_regression$VideoIndex[i] = 'third'
}else if(df_regression$Stimuli_Sequence[i] == 'NeuPosNeg'){
df_regression$VideoIndex[i] = 'second'
}else if(df_regression$Stimuli_Sequence[i] == 'NeuNegPos'){
df_regression$VideoIndex[i] = 'third'
}
}else if(df_regression$VideoValence[i] == 'negative'){
if(df_regression$Stimuli_Sequence[i] == 'PosNegNeu'){
df_regression$VideoIndex[i] = 'second'
}else if(df_regression$Stimuli_Sequence[i] == 'PosNeuNeg'){
df_regression$VideoIndex[i] = 'third'
}else if(df_regression$Stimuli_Sequence[i] == 'NegPosNeu'){
df_regression$VideoIndex[i] = 'first'
}else if(df_regression$Stimuli_Sequence[i] == 'NegNeuPos'){
df_regression$VideoIndex[i] = 'first'
}else if(df_regression$Stimuli_Sequence[i] == 'NeuPosNeg'){
df_regression$VideoIndex[i] = 'third'
}else if(df_regression$Stimuli_Sequence[i] == 'NeuNegPos'){
df_regression$VideoIndex[i] = 'second'
}
}else if(df_regression$VideoValence[i] == 'neutral'){
if(df_regression$Stimuli_Sequence[i] == 'PosNegNeu'){
df_regression$VideoIndex[i] = 'third'
}else if(df_regression$Stimuli_Sequence[i] == 'PosNeuNeg'){
df_regression$VideoIndex[i] = 'second'
}else if(df_regression$Stimuli_Sequence[i] == 'NegPosNeu'){
df_regression$VideoIndex[i] = 'third'
}else if(df_regression$Stimuli_Sequence[i] == 'NegNeuPos'){
df_regression$VideoIndex[i] = 'second'
}else if(df_regression$Stimuli_Sequence[i] == 'NeuPosNeg'){
df_regression$VideoIndex[i] = 'first'
}else if(df_regression$Stimuli_Sequence[i] == 'NeuNegPos'){
df_regression$VideoIndex[i] = 'first'
}
}
}
# ORDER FACTORS
df_regression$Participant        <- factor(df_regression$Participant)
df_regression$Dyad               <- factor(df_regression$Dyad)
df_regression$Stimuli_Sequence   <- factor(df_regression$Stimuli_Sequence)
df_regression$AttentionType      <- factor(df_regression$AttentionType, levels = c('disjoint', 'joint'))
df_regression$Gender             <- factor(df_regression$Gender, levels = c('male', 'female'))
df_regression$VideoValence       <- factor(df_regression$VideoValence, levels = c('neutral', 'negative', 'positive'))
df_regression$VideoIndex         <- factor(df_regression$VideoIndex, levels = c('first', 'second', 'third'))
df_regression$JointExposure      <- factor(df_regression$exposure_MEAN, levels = c(0.5, 0, 1))
df_regression                    <- subset(df_regression, select = -c(exposure_MEAN))
# REPLACE ARTEFACTUAL PHYSIOLOGY VALUES BY NA
df_regression  <- df_regression %>%
group_by(Dyad, VideoValence)%>% mutate(
ecg_IBI_MEAN        = if_else(ecg_clean == 'n', NA_real_, ecg_IBI_MEAN),
ecg_WTC_LF          = if_else(ecg_clean == 'n', NA_real_, ecg_WTC_LF),
ecg_WTC_HF          = if_else(ecg_clean == 'n', NA_real_, ecg_WTC_HF),
rsp_Rate_MEAN       = if_else(rsp_clean == 'n', NA_real_, rsp_Rate_MEAN),
rsp_WTC             = if_else(rsp_clean == 'n', NA_real_, rsp_WTC),
eda_Amplitude_MEAN  = if_else(eda_clean == 'n', NA_real_, eda_Amplitude_MEAN),
eda_PeakNumber_MEAN = if_else(eda_clean == 'n', NA_real_, eda_PeakNumber_MEAN),
eda_WTC             = if_else(eda_clean == 'n', NA_real_, eda_WTC)) %>%
ungroup()
# CREATE A DATA FRAME WITH DIFFERENCES POS-NEU AND NEG-NEU ####
df_regression_baseline <- df_regression %>% group_by(Dyad) %>%
mutate(ecg_IBI_MEAN         = ecg_IBI_MEAN - ecg_IBI_MEAN[VideoValence == 'neutral'],
ecg_WTC_LF           = ecg_WTC_LF - ecg_WTC_LF[VideoValence == 'neutral'],
ecg_WTC_HF           = ecg_WTC_HF - ecg_WTC_HF[VideoValence == 'neutral'],
rsp_Rate_MEAN        = rsp_Rate_MEAN - rsp_Rate_MEAN[VideoValence == 'neutral'],
rsp_WTC              = rsp_WTC - rsp_WTC[VideoValence == 'neutral'],
eda_Amplitude_MEAN   = eda_Amplitude_MEAN - eda_Amplitude_MEAN[VideoValence == 'neutral'],
eda_WTC              = eda_WTC - eda_WTC[VideoValence == 'neutral'],
sam_valence_MEAN     = sam_valence_MEAN - sam_valence_MEAN[VideoValence == 'neutral'],
sam_valence_DIFF     = sam_valence_DIFF - sam_valence_DIFF[VideoValence == 'neutral'],
sam_arousal_MEAN     = sam_arousal_MEAN - sam_arousal_MEAN[VideoValence == 'neutral'],
sam_arousal_DIFF     = sam_arousal_DIFF - sam_arousal_DIFF[VideoValence == 'neutral'],
connectedness_MEAN   = connectedness_MEAN - connectedness_MEAN[VideoValence == 'neutral'],
connectedness_PROD   = connectedness_PROD - connectedness_PROD[VideoValence == 'neutral']
)  %>% subset(VideoValence != 'neutral')  %>% ungroup()
# PREPARE DATA ####
myDF <- df_regression_baseline %>%
mutate(
Age_MEAN                 = scale(Age_MEAN),
Age_DIFF                 = scale(Age_DIFF),
PRE_anxiety_state_MEAN   = scale(PRE_anxiety_state_MEAN),
PRE_anxiety_trait_MEAN   = scale(PRE_anxiety_trait_MEAN),
PRE_we_PROD              = scale(PRE_we_PROD),
PRE_dfi_PROD             = scale(PRE_dfi_PROD),
PRE_anxiety_state_DIFF   = scale(PRE_anxiety_state_DIFF),
PRE_anxiety_trait_DIFF   = scale(PRE_anxiety_trait_DIFF),
sam_arousal_MEAN         = scale(sam_arousal_MEAN),
sam_arousal_DIFF         = scale(sam_arousal_DIFF),
sam_valence_MEAN         = scale(sam_valence_MEAN),
sam_valence_DIFF         = scale(sam_valence_DIFF),
ecg_IBI_MEAN             = scale(ecg_IBI_MEAN),
rsp_Rate_MEAN            = scale(rsp_Rate_MEAN),
eda_Amplitude_MEAN       = scale(eda_Amplitude_MEAN),
ecg_WTC_LF               = scale(ecg_WTC_LF),
ecg_WTC_HF               = scale(ecg_WTC_HF),
rsp_WTC                  = scale(rsp_WTC),
eda_WTC                  = scale(eda_WTC),
ecg_rC                   = scale(ecg_rC),
ecg_rCC                  = scale(ecg_rCC),
rsp_rC                   = scale(rsp_rC),
rsp_rCC                  = scale(rsp_rCC),
eda_rC                   = scale(eda_rC),
eda_rCC                  = scale(eda_rCC))
# FULL MODEL (NO MODEL SELECTION) ####
mod_connectedness_full <- lmer(connectedness_PROD ~
AttentionType +
Gender +
Age_MEAN +
Age_DIFF +
PRE_we_PROD*AttentionType +
PRE_dfi_PROD*AttentionType +
sam_arousal_MEAN*AttentionType  +
sam_valence_MEAN*AttentionType   +
ecg_IBI_MEAN*AttentionType +
rsp_Rate_MEAN*AttentionType +
eda_Amplitude_MEAN*AttentionType +
(1|Dyad),
myDF,
REML = F)
summary(mod_connectedness_full)
check_model(mod_connectedness_full)
vif(mod_connectedness_full)
model_performance(mod_connectedness_full)
tab_model(mod_connectedness_full,
df.method = "satterthwaite",
show.df = TRUE)
# SELECTED MODEL (BACKWARD) ####
allFit(mod_connectedness_full)
build_connectedness_select <- buildmer(connectedness_PROD ~
AttentionType +
Gender +
Age_MEAN +
Age_DIFF +
PRE_we_PROD*AttentionType +
PRE_dfi_PROD*AttentionType +
sam_arousal_MEAN*AttentionType  +
sam_valence_MEAN*AttentionType   +
ecg_IBI_MEAN *AttentionType  +
rsp_Rate_MEAN*AttentionType +
eda_Amplitude_MEAN*AttentionType +
(1+sam_arousal_MEAN+sam_valence_MEAN+ecg_IBI_MEAN+rsp_Rate_MEAN+eda_Amplitude_MEAN|Dyad) +
(1+sam_arousal_MEAN+sam_valence_MEAN+ecg_IBI_MEAN+rsp_Rate_MEAN+eda_Amplitude_MEAN|VideoIndex) +
(1+sam_arousal_MEAN+sam_valence_MEAN+ecg_IBI_MEAN+rsp_Rate_MEAN+eda_Amplitude_MEAN|VideoValence),
myDF,
buildmerControl= buildmerControl(
direction = c('order',"backward"),
crit='LRT',
args =list(control=lmerControl(optimizer="optimx",
optCtrl = list(method = 'Nelder_Mead ')))))
mod_connectedness_select <- lmer(build_connectedness_select@model@call[["formula"]],
REML = F,
myDF)
summary(mod_connectedness_select)
check_model(mod_connectedness_select)
vif(mod_connectedness_select)
model_performance(mod_connectedness_select)
tab_model(mod_connectedness_select,
df.method = "satterthwaite",
show.df = TRUE)
# CARDIAC SYNCHRONY ####
# PREPARE DATA
myDF <- df_regression %>% filter(ecg_clean == "y") %>%
dplyr::select(Dyad,
VideoValence,
VideoIndex,
ecg_rC,
ecg_rCC,
ecg_WTC_LF,
ecg_WTC_HF,
ecg_PLV,
ecg_CRQA_DET,
ecg_CRQA_ADL,
ecg_CRQA_MDL,
ecg_CRQA_ENT,
ecg_CRQA_LAM)  %>%
group_by(Dyad) %>%
dplyr::summarize(ecg_rC = mean(ecg_rC, na.rm = T),
ecg_rCC = mean(ecg_rCC, na.rm = T),
ecg_WTC_LF = mean(ecg_WTC_LF, na.rm = T),
ecg_WTC_HF = mean(ecg_WTC_HF, na.rm = T),
ecg_PLV = mean(ecg_PLV, na.rm = T),
ecg_CRQA_DET = mean(ecg_CRQA_DET, na.rm = T),
ecg_CRQA_ADL = mean(ecg_CRQA_ADL, na.rm = T),
ecg_CRQA_MDL = mean(ecg_CRQA_MDL, na.rm = T),
ecg_CRQA_ENT = mean(ecg_CRQA_ENT, na.rm = T),
ecg_CRQA_LAM = mean(ecg_CRQA_LAM, na.rm = T)) %>%
ungroup()
# REMOVE OUTLIERS
extreme_rC    <- myDF %>% identify_outliers(ecg_rC)
outliers_rC   <- unique(as.character(extreme_rC$Dyad))
no_outlier_rC <- unique(setdiff(myDF$Dyad, outliers_rC))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_rC)
extreme_rCC    <- myDF %>% identify_outliers(ecg_rCC)
outliers_rCC   <- unique(as.character(extreme_rCC$Dyad))
no_outlier_rCC <- unique(setdiff(myDF$Dyad, outliers_rCC))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_rCC)
extreme_WTC_LF    <- myDF %>% identify_outliers(ecg_WTC_LF)
outliers_WTC_LF   <- unique(as.character(extreme_WTC_LF$Dyad))
no_outlier_WTC_LF <- unique(setdiff(myDF$Dyad, outliers_WTC_LF))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_WTC_LF)
extreme_WTC_HF    <- myDF %>% identify_outliers(ecg_WTC_HF)
outliers_WTC_HF   <- unique(as.character(extreme_WTC_HF$Dyad))
no_outlier_WTC_HF <- unique(setdiff(myDF$Dyad, outliers_WTC_HF))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_WTC_HF)
extreme_PLV    <- myDF %>% identify_outliers(ecg_PLV)
outliers_PLV   <- unique(as.character(extreme_PLV$Dyad))
no_outlier_PLV <- unique(setdiff(myDF$Dyad, outliers_PLV))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_PLV)
extreme_CRQA_DET    <- myDF %>% identify_outliers(ecg_CRQA_DET)
outliers_CRQA_DET   <- unique(as.character(extreme_CRQA_DET$Dyad))
no_outlier_CRQA_DET <- unique(setdiff(myDF$Dyad, outliers_CRQA_DET))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_CRQA_DET)
myDF       <- myDF %>% filter(Dyad %in% no_outlier_PLV)
extreme_CRQA_ADL    <- myDF %>% identify_outliers(ecg_CRQA_ADL)
outliers_CRQA_ADL   <- unique(as.character(extreme_CRQA_ADL$Dyad))
no_outlier_CRQA_ADL <- unique(setdiff(myDF$Dyad, outliers_CRQA_ADL))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_CRQA_ADL)
extreme_CRQA_MDL    <- myDF %>% identify_outliers(ecg_CRQA_MDL)
outliers_CRQA_MDL   <- unique(as.character(extreme_CRQA_MDL$Dyad))
no_outlier_CRQA_MDL <- unique(setdiff(myDF$Dyad, outliers_CRQA_MDL))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_CRQA_MDL)
extreme_CRQA_ENT    <- myDF %>% identify_outliers(ecg_CRQA_ENT)
outliers_CRQA_ENT   <- unique(as.character(extreme_CRQA_ENT$Dyad))
no_outlier_CRQA_ENT <- unique(setdiff(myDF$Dyad, outliers_CRQA_ENT))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_CRQA_ENT)
extreme_CRQA_LAM    <- myDF %>% identify_outliers(ecg_CRQA_LAM)
outliers_CRQA_LAM   <- unique(as.character(extreme_CRQA_LAM$Dyad))
no_outlier_CRQA_LAM <- unique(setdiff(myDF$Dyad, outliers_CRQA_LAM))
myDF       <- myDF %>% filter(Dyad %in% no_outlier_CRQA_LAM)
# DESCRIPTIVE STATISTICS
stats_desc <- myDF %>%
get_summary_stats(everything(), type = 'mean_sd')
stats_desc
# PLOT
pairs.panels(myDF,
smooth = TRUE,      # If TRUE, draws loess smooths
scale = FALSE,      # If TRUE, scales the correlation text font
density = TRUE,     # If TRUE, adds density plots and histograms
ellipses = FALSE,    # If TRUE, draws ellipses
method = "spearman", # Correlation method (also "spearman" or "kendall")
pch = 21,           # pch symbol
lm = TRUE,         # If TRUE, plots linear fit rather than the LOESS (smoothed) fit
cor = TRUE,         # If TRUE, reports correlations
jiggle = FALSE,     # If TRUE, data points are jittered
factor = 2,         # Jittering factor
hist.col = 4,       # Histograms color
stars = TRUE,       # If TRUE, adds significance level with stars
ci = TRUE)          # If TRUE, adds confidence intervals
# SPEARMAN RANK CORRELATIONS ACROSS CONDITIONS
myDFcorr <- myDF %>% dplyr::select(-c(1))
rcor.test(data.matrix(myDFcorr),
method = "spearman",
p.adjust = FALSE,
p.adjust.method = "holm")
