pacman::p_load(tidyr,dplyr,ggplot2,gridExtra,psych,stats,osfr,rio)Dyadic Analysis workshop, Zurich Nov. 24: Dyadic Sync
Dyadic Sync
Libraries
Loading the data from OSF
file <- osf_retrieve_file("6731ddb79a41a311408eaa69")
setwd("C:\\Users\\97254\\Downloads")
osf_download(file, progress = T,conflicts = "overwrigt")
temp2 <- import("SyncData.csv")
head(temp2)Example data
foo <- temp2 %>% filter(Couple_Int=="104_1")
#plot of one interaction
cor <- round(cor(foo$M_Close,foo$W_Close,use="pairwise.complete.obs"),3)
# Plot using ggplot2
plot1 <- ggplot(foo) +
geom_line(aes(x = Time, y = M_Close), color = "blue") +
geom_line(aes(x = Time, y = W_Close), color = "red") +
geom_hline(yintercept = 50, color = "black", linetype = "dashed") +
labs(
x = "Time (Sec)",
y = "Closeness",
title = paste0("Correlation= ", cor)
) +
ylim(0, 100) +
theme_minimal()
plot1
Cross Correlation Function
Global CCF
data <- foo
window_size <- 50 #(about +/-3 sec)
ccf_res_glob <- ccf(data$M_Close, data$W_Close, plot = T, lag.max = window_size)
ccfs_glob <- ccf_res_glob$acf
(ccfs_glob), , 1
[,1]
[1,] 0.5305032
[2,] 0.5303950
[3,] 0.5302753
[4,] 0.5301365
[5,] 0.5299895
[6,] 0.5298391
[7,] 0.5296871
[8,] 0.5295293
[9,] 0.5293592
[10,] 0.5291766
[11,] 0.5289915
[12,] 0.5287899
[13,] 0.5285792
[14,] 0.5283702
[15,] 0.5281595
[16,] 0.5279488
[17,] 0.5277373
[18,] 0.5275216
[19,] 0.5273035
[20,] 0.5270853
[21,] 0.5268498
[22,] 0.5266003
[23,] 0.5263375
[24,] 0.5260723
[25,] 0.5258078
[26,] 0.5255327
[27,] 0.5252426
[28,] 0.5249509
[29,] 0.5246542
[30,] 0.5243559
[31,] 0.5240410
[32,] 0.5237237
[33,] 0.5234006
[34,] 0.5230692
[35,] 0.5227628
[36,] 0.5224506
[37,] 0.5221335
[38,] 0.5218172
[39,] 0.5214901
[40,] 0.5211573
[41,] 0.5208252
[42,] 0.5204916
[43,] 0.5201263
[44,] 0.5197585
[45,] 0.5193850
[46,] 0.5190007
[47,] 0.5186064
[48,] 0.5182114
[49,] 0.5178056
[50,] 0.5173957
[51,] 0.5169882
[52,] 0.5159829
[53,] 0.5149692
[54,] 0.5139449
[55,] 0.5129114
[56,] 0.5118812
[57,] 0.5108445
[58,] 0.5098060
[59,] 0.5087593
[60,] 0.5077102
[61,] 0.5066569
[62,] 0.5056011
[63,] 0.5045271
[64,] 0.5034465
[65,] 0.5023659
[66,] 0.5012944
[67,] 0.5002180
[68,] 0.4991415
[69,] 0.4980659
[70,] 0.4969828
[71,] 0.4958989
[72,] 0.4948167
[73,] 0.4937278
[74,] 0.4926414
[75,] 0.4915460
[76,] 0.4904381
[77,] 0.4893236
[78,] 0.4882042
[79,] 0.4870641
[80,] 0.4859207
[81,] 0.4847706
[82,] 0.4836123
[83,] 0.4824532
[84,] 0.4812924
[85,] 0.4801259
[86,] 0.4789535
[87,] 0.4777812
[88,] 0.4765989
[89,] 0.4754092
[90,] 0.4742120
[91,] 0.4730157
[92,] 0.4718227
[93,] 0.4706313
[94,] 0.4694350
[95,] 0.4682387
[96,] 0.4670481
[97,] 0.4658609
[98,] 0.4646836
[99,] 0.4635104
[100,] 0.4623339
[101,] 0.4611591Time de-trending
# Compute global CCF for detrended time series across all data
data_detrended_glob <- data %>%
mutate(
M_Close_detrended = residuals(lm(M_Close ~ Time, data = data)),
W_Close_detrended = residuals(lm(W_Close ~ Time, data = data))
)
ccf_res_glob_detrended <- ccf(data_detrended_glob$M_Close_detrended, data_detrended_glob$W_Close_detrended, plot = T, lag.max = window_size)
ccfs_glob_detrended <- ccf_res_glob_detrended$acf
(ccfs_glob_detrended), , 1
[,1]
[1,] 0.6166339
[2,] 0.6161443
[3,] 0.6156404
[4,] 0.6151130
[5,] 0.6145755
[6,] 0.6140340
[7,] 0.6134906
[8,] 0.6129402
[9,] 0.6123744
[10,] 0.6117934
[11,] 0.6112094
[12,] 0.6106050
[13,] 0.6099894
[14,] 0.6093760
[15,] 0.6087607
[16,] 0.6081456
[17,] 0.6075295
[18,] 0.6069084
[19,] 0.6062843
[20,] 0.6056604
[21,] 0.6050151
[22,] 0.6043524
[23,] 0.6036734
[24,] 0.6029915
[25,] 0.6023107
[26,] 0.6016167
[27,] 0.6009043
[28,] 0.6001900
[29,] 0.5994697
[30,] 0.5987474
[31,] 0.5980047
[32,] 0.5972591
[33,] 0.5965064
[34,] 0.5957436
[35,] 0.5949967
[36,] 0.5942428
[37,] 0.5934829
[38,] 0.5927241
[39,] 0.5919521
[40,] 0.5911731
[41,] 0.5903952
[42,] 0.5896154
[43,] 0.5888117
[44,] 0.5880050
[45,] 0.5871912
[46,] 0.5863642
[47,] 0.5855251
[48,] 0.5846850
[49,] 0.5838317
[50,] 0.5829734
[51,] 0.5821183
[52,] 0.5805356
[53,] 0.5789425
[54,] 0.5773360
[55,] 0.5757180
[56,] 0.5741041
[57,] 0.5724818
[58,] 0.5708574
[59,] 0.5692226
[60,] 0.5675847
[61,] 0.5659415
[62,] 0.5642951
[63,] 0.5626260
[64,] 0.5609486
[65,] 0.5592710
[66,] 0.5576047
[67,] 0.5559321
[68,] 0.5542594
[69,] 0.5525877
[70,] 0.5509066
[71,] 0.5492243
[72,] 0.5475440
[73,] 0.5458554
[74,] 0.5441698
[75,] 0.5424727
[76,] 0.5407602
[77,] 0.5390393
[78,] 0.5373122
[79,] 0.5355592
[80,] 0.5338021
[81,] 0.5320366
[82,] 0.5302608
[83,] 0.5284839
[84,] 0.5267047
[85,] 0.5249183
[86,] 0.5231246
[87,] 0.5213308
[88,] 0.5195246
[89,] 0.5177090
[90,] 0.5158841
[91,] 0.5140600
[92,] 0.5122400
[93,] 0.5104220
[94,] 0.5085977
[95,] 0.5067732
[96,] 0.5049559
[97,] 0.5031426
[98,] 0.5013415
[99,] 0.4995454
[100,] 0.4977451
[101,] 0.4959468Windowed CCF
data <- foo
window_size <- 50 #(about +/-3 sec)
slide_step <- 10
max_lag <- 5
window_size_samples <- window_size
slide_step_samples <- slide_step
Min_Time <- min(data$Time,na.rm = T)
Max_Time <- max(data$Time,na.rm = T)
# Initialize a list to store results
ccf_results <- list()
# Loop over each window
for (start_idx in seq(Min_Time, Max_Time- window_size_samples, by = slide_step_samples)) {
# start_idx <- seq(Min_Time, Max_Time- window_size_samples, by = slide_step_samples)[1]
# Define the end index of the window
end_idx <- start_idx + window_size_samples
# Extract the window for each time series
datai <- data %>% filter(Time>=start_idx,Time<end_idx) %>% select(M_Close,W_Close)
# Calculate the CCF for this window (lag.max can be adjusted as needed)
ccf_res <- ccf(datai$M_Close, datai$W_Close, plot = F, lag.max = max_lag)
ccfs <- ccf_res$acf
#handling constant values (no correlation without variability)
if (sd(datai$M_Close)==0|sd(datai$W_Close)==0){
datai <- data.frame(Start_Time=start_idx,End_Time=end_idx, CCF_0=NA,CCF_Max=NA,CCF_Max_Lag=NA,
CCF_avg=NA,CCF_Max_Abs=NA)
}else{datai <- data.frame(Start_Time=start_idx,End_Time=end_idx, CCF_0=ccfs[max_lag+1],CCF_Max=max(ccfs),CCF_Max_Lag=(which(ccfs==max(ccfs))-(max_lag+1)),
CCF_avg=mean(ccfs),CCF_Max_Abs=max(abs(ccfs)))
}
# Store the CCF values for this window
ccf_results[[paste("Time_",start_idx)]] <- datai # Only storing CCF values
}
CCF <- dplyr::bind_rows(ccf_results)
CCF Start_Time End_Time CCF_0 CCF_Max CCF_Max_Lag CCF_avg
1 0.12 50.12 0.83543584 0.83543584 0 0.82800771
2 10.12 60.12 0.82566474 0.82566474 0 0.81976897
3 20.12 70.12 0.86673782 0.86673782 0 0.86011078
4 30.12 80.12 0.93034784 0.93034784 0 0.92128800
5 40.12 90.12 0.72713115 0.72861092 -5 0.72035403
6 50.12 100.12 0.41447335 0.41953555 -5 0.40260388
7 60.12 110.12 0.75070890 0.75276066 -5 0.74624940
8 70.12 120.12 0.61908522 0.63082429 -5 0.61779350
9 80.12 130.12 0.01410805 0.02690825 -5 0.01005804
10 90.12 140.12 -0.57606685 -0.56550164 -5 -0.57334008
11 100.12 150.12 -0.64651572 -0.64109016 -5 -0.64932979
12 110.12 160.12 0.59582358 0.59582358 0 0.58952276
13 120.12 170.12 0.93237971 0.93298863 5 0.92897883
14 130.12 180.12 0.91841700 0.91911485 5 0.91618236
15 140.12 190.12 0.32375964 0.33094438 -5 0.31583975
16 150.12 200.12 0.74676192 0.76447800 -5 0.74441452
17 160.12 210.12 0.62364593 0.64412717 -5 0.62249043
18 170.12 220.12 0.29574476 0.32152533 -5 0.29421422
19 180.12 230.12 -0.24929176 -0.21900901 -5 -0.25388506
20 190.12 240.12 0.52490306 0.53273790 -5 0.52337615
21 200.12 250.12 0.57146016 0.57146016 0 0.56817736
22 210.12 260.12 0.40995485 0.40995485 0 0.40768161
23 220.12 270.12 0.14470088 0.14957294 5 0.14562040
24 230.12 280.12 0.07709515 0.07974169 5 0.07575467
25 240.12 290.12 0.47552090 0.47576777 1 0.46604372
26 250.12 300.12 0.73141356 0.73141356 0 0.72586752
27 260.12 310.12 0.66815032 0.67484318 5 0.66783906
28 270.12 320.12 0.35929131 0.35929131 0 0.35704009
29 280.12 330.12 0.70499453 0.70499453 0 0.70253969
30 290.12 340.12 0.66852851 0.67511304 -5 0.66305347
31 300.12 350.12 0.59786067 0.60205108 -5 0.59181611
CCF_Max_Abs
1 0.83543584
2 0.82566474
3 0.86673782
4 0.93034784
5 0.72861092
6 0.41953555
7 0.75276066
8 0.63082429
9 0.02690825
10 0.57629092
11 0.66099354
12 0.59582358
13 0.93298863
14 0.91911485
15 0.33094438
16 0.76447800
17 0.64412717
18 0.32152533
19 0.29646227
20 0.53273790
21 0.57146016
22 0.40995485
23 0.14957294
24 0.07974169
25 0.47576777
26 0.73141356
27 0.67484318
28 0.35929131
29 0.70499453
30 0.67511304
31 0.60205108Plotting Results
# Add global CCF indices for both original and detrended time series to CCF data frame
CCF <- CCF %>%
mutate(
# Original time series global CCF
CCF_glob_0 = ccfs_glob[max_lag + 1],
CCF_glob_Max = max(ccfs_glob),
CCF_glob_Max_Lag = which(ccfs_glob == max(ccfs_glob)) - (max_lag + 1),
CCF_glob_avg = mean(ccfs_glob),
CCF_glob_Max_Abs = max(abs(ccfs_glob)),
# Detrended time series global CCF
CCF_glob_0_detrended = ccfs_glob_detrended[max_lag + 1],
CCF_glob_Max_detrended = max(ccfs_glob_detrended),
CCF_glob_Max_Lag_detrended = which(ccfs_glob_detrended == max(ccfs_glob_detrended)) - (max_lag + 1),
CCF_glob_avg_detrended = mean(ccfs_glob_detrended),
CCF_glob_Max_Abs_detrended = max(abs(ccfs_glob_detrended))
)
graph.1 <- ggplot(data) +
geom_line(aes(x= Time , y= M_Close),color='blue',linewidth=1) +
geom_line(aes(x= Time , y= W_Close),color='orange',linewidth=1) +
# scale_x_continuous(name = "Time (seconds)", breaks =seq(0,300,30)) +
# guides(colour = guide_legend(override.aes = list(size=3)))+
ylab('Closeness')+xlab('Time') +
ggtitle(paste0("Couple_Int: ",data$Couple_Int,"\n",
"CCF_glob_0: ", round(ccfs_glob[max_lag+1],3)))+
theme_bw()+geom_vline(xintercept = seq(Min_Time,Max_Time,window_size_samples),
linetype="dotted",
color = "red", linewidth=1)
graph.2<-ggplot(data=CCF,aes(x=(End_Time-window_size/2) ,
y=CCF_0)) +
geom_line(size=1)+
geom_point(size=2,color="red")+
scale_x_continuous(name = "Time") +
scale_y_continuous(name = "Max CCF",limits = c(-1,1),breaks =seq(-1,1,1))+
geom_vline(xintercept = seq(Min_Time,Max_Time,window_size_samples),
linetype="dotted",
color = "red", linewidth=1)+
theme_bw()
gridExtra::grid.arrange(graph.1, graph.2, ncol = 1) 
Wrapping everything within a function
compute_ccf <- function(data, window_size, slide_step, max_lag) {
# Convert to number of samples (assuming 1 sample per second)
window_size_samples <- window_size
slide_step_samples <- slide_step
# Determine minimum and maximum time in the dataset
Min_Time <- min(data$Time, na.rm = TRUE)
Max_Time <- max(data$Time, na.rm = TRUE)
# Initialize a list to store results
ccf_results <- list()
# Loop over each window
for (start_idx in seq(Min_Time, Max_Time - window_size_samples, by = slide_step_samples)) {
# Define the end index of the window
end_idx <- start_idx + window_size_samples
# Extract the window for each time series
datai <- data %>% filter(Time >= start_idx, Time < end_idx) %>% select(Time, M_Close, W_Close)
# Detrend each time series by fitting a linear model and using residuals
datai <- datai %>%
mutate(
M_Close_detrended = residuals(lm(M_Close ~ Time, data = datai)),
W_Close_detrended = residuals(lm(W_Close ~ Time, data = datai))
)
# Calculate the CCF for the original data in this window
ccf_res <- ccf(datai$M_Close, datai$W_Close, plot = FALSE, lag.max = max_lag)
ccfs <- ccf_res$acf
# Calculate the CCF for the detrended data in this window
ccf_res_detrended <- ccf(datai$M_Close_detrended, datai$W_Close_detrended, plot = FALSE, lag.max = max_lag)
ccfs_detrended <- ccf_res_detrended$acf
# Handle constant values (no correlation without variability)
if (sd(datai$M_Close) == 0 | sd(datai$W_Close) == 0) {
datai_result <- data.frame(
Start_Time = start_idx,
End_Time = end_idx,
avg_M_Close = mean(datai$M_Close, na.rm = TRUE),
avg_W_Close = mean(datai$W_Close, na.rm = TRUE),
CCF_0 = NA,
CCF_Max = NA,
CCF_Max_Lag = NA,
CCF_avg = NA,
CCF_Max_Abs = NA,
CCF_0_detrended = NA,
CCF_Max_detrended = NA,
CCF_Max_Lag_detrended = NA,
CCF_avg_detrended = NA,
CCF_Max_Abs_detrended = NA
)
} else {
datai_result <- data.frame(
Start_Time = start_idx,
End_Time = end_idx,
avg_M_Close = mean(datai$M_Close, na.rm = TRUE),
avg_W_Close = mean(datai$W_Close, na.rm = TRUE),
CCF_0 = ccfs[max_lag + 1],
CCF_Max = max(ccfs),
CCF_Max_Lag = which(ccfs == max(ccfs)) - (max_lag + 1),
CCF_avg = mean(ccfs),
CCF_Max_Abs = max(abs(ccfs)),
CCF_0_detrended = ccfs_detrended[max_lag + 1],
CCF_Max_detrended = max(ccfs_detrended),
CCF_Max_Lag_detrended = which(ccfs_detrended == max(ccfs_detrended)) - (max_lag + 1),
CCF_avg_detrended = mean(ccfs_detrended),
CCF_Max_Abs_detrended = max(abs(ccfs_detrended))
)
}
# Store the CCF values for this window
ccf_results[[paste("Time_", start_idx)]] <- datai_result
}
# Combine results into a single data frame
CCF <- dplyr::bind_rows(ccf_results)
# Compute global CCF for original time series
ccf_res_glob <- ccf(data$M_Close, data$W_Close, plot = FALSE, lag.max = max_lag)
ccfs_glob <- ccf_res_glob$acf
# Compute global CCF for detrended time series
data_detrended_glob <- data %>%
mutate(
M_Close_detrended = residuals(lm(M_Close ~ Time, data = data)),
W_Close_detrended = residuals(lm(W_Close ~ Time, data = data))
)
ccf_res_glob_detrended <- ccf(data_detrended_glob$M_Close_detrended, data_detrended_glob$W_Close_detrended, plot = FALSE, lag.max = max_lag)
ccfs_glob_detrended <- ccf_res_glob_detrended$acf
# Add global CCF indices for both original and detrended time series to CCF data frame
CCF <- CCF %>%
mutate(
avg_glob_M_Close=mean(data$M_Close,na.rm=T),
avg_glob_W_Close=mean(data$W_Close,na.rm=T),
CCF_glob_0 = ccfs_glob[max_lag + 1],
CCF_glob_Max = max(ccfs_glob),
CCF_glob_Max_Lag = which(ccfs_glob == max(ccfs_glob)) - (max_lag + 1),
CCF_glob_avg = mean(ccfs_glob),
CCF_glob_Max_Abs = max(abs(ccfs_glob)),
CCF_glob_0_detrended = ccfs_glob_detrended[max_lag + 1],
CCF_glob_Max_detrended = max(ccfs_glob_detrended),
CCF_glob_Max_Lag_detrended = which(ccfs_glob_detrended == max(ccfs_glob_detrended)) - (max_lag + 1),
CCF_glob_avg_detrended = mean(ccfs_glob_detrended),
CCF_glob_Max_Abs_detrended = max(abs(ccfs_glob_detrended))
)
return(CCF)
}Using the function
CCF <- compute_ccf(data=foo, window_size=20, slide_step=20, max_lag=50)
CCF Start_Time End_Time avg_M_Close avg_W_Close CCF_0 CCF_Max
1 0.12 20.12 50.00000 49.57485 NA NA
2 20.12 40.12 53.54655 62.28228 0.8987949 0.89879487
3 40.12 60.12 78.62462 78.31832 0.9454401 0.94544005
4 60.12 80.12 96.01198 89.83234 0.2795919 0.27959189
5 80.12 100.12 81.17417 86.79580 0.4056604 0.56629427
6 100.12 120.12 75.76276 72.01802 -0.7003380 -0.08900336
7 120.12 140.12 86.11677 64.50000 -0.3184014 -0.13081794
8 140.12 160.12 89.36336 80.67267 0.9179189 0.91791895
9 160.12 180.12 95.66366 96.00000 NA NA
10 180.12 200.12 66.46707 79.07784 0.4158366 0.67013318
11 200.12 220.12 82.94294 66.11411 0.4044937 0.42625411
12 220.12 240.12 93.84685 71.93093 0.3033036 0.37617907
13 240.12 260.12 92.08683 66.73054 -0.4220888 -0.04367101
14 260.12 280.12 94.11712 65.82883 -0.4635128 0.07651078
15 280.12 300.12 97.94595 89.74174 0.3765581 0.46698199
16 300.12 320.12 93.17365 87.89222 0.5949246 0.59492456
17 320.12 340.12 86.01802 73.24625 0.2945438 0.49707332
CCF_Max_Lag CCF_avg CCF_Max_Abs CCF_0_detrended CCF_Max_detrended
1 NA NA NA NA NA
2 0 0.74082429 0.8987949 0.24350097 0.272588802
3 0 0.76026355 0.9454401 0.64198534 0.645868346
4 0 0.06492932 0.2795919 0.09090530 0.120988390
5 -50 0.30525483 0.5662943 -0.75846215 -0.003184127
6 -50 -0.55420524 0.7159108 -0.45053248 0.127557504
7 -50 -0.29860048 0.3693053 0.08225300 0.104999197
8 0 0.78860265 0.9179189 0.32813572 0.632996229
9 NA NA NA NA NA
10 -50 0.30326510 0.6701332 0.15447547 0.627933291
11 -50 0.33275245 0.4262541 -0.55774619 0.368185855
12 7 0.26638247 0.3761791 0.04533115 0.217779925
13 50 -0.30188284 0.5184244 -0.40036801 0.372552251
14 50 -0.32844588 0.5331094 -0.04415145 -0.027293094
15 50 0.21755191 0.4669820 0.01773973 0.422855037
16 0 0.40247190 0.5949246 0.09739412 0.298984894
17 -26 0.27670706 0.4970733 0.32604827 0.442633965
CCF_Max_Lag_detrended CCF_avg_detrended CCF_Max_Abs_detrended
1 NA NA NA
2 6 0.12592858 0.2725888
3 2 0.36563625 0.6458683
4 -3 -0.06733808 0.2167454
5 50 -0.43222512 0.7584621
6 -50 -0.28816176 0.4543157
7 -13 0.05228603 0.1049992
8 37 0.24854320 0.6329962
9 NA NA NA
10 -50 0.16376178 0.6279333
11 50 -0.35261232 0.6362429
12 -50 0.03181253 0.2177799
13 50 -0.28199643 0.6931544
14 22 -0.12450449 0.3884496
15 50 0.05113643 0.4228550
16 -48 0.09691391 0.2989849
17 -19 0.16260746 0.4426340
avg_glob_M_Close avg_glob_W_Close CCF_glob_0 CCF_glob_Max CCF_glob_Max_Lag
1 83.24004 74.4319 0.5169882 0.5305032 -50
2 83.24004 74.4319 0.5169882 0.5305032 -50
3 83.24004 74.4319 0.5169882 0.5305032 -50
4 83.24004 74.4319 0.5169882 0.5305032 -50
5 83.24004 74.4319 0.5169882 0.5305032 -50
6 83.24004 74.4319 0.5169882 0.5305032 -50
7 83.24004 74.4319 0.5169882 0.5305032 -50
8 83.24004 74.4319 0.5169882 0.5305032 -50
9 83.24004 74.4319 0.5169882 0.5305032 -50
10 83.24004 74.4319 0.5169882 0.5305032 -50
11 83.24004 74.4319 0.5169882 0.5305032 -50
12 83.24004 74.4319 0.5169882 0.5305032 -50
13 83.24004 74.4319 0.5169882 0.5305032 -50
14 83.24004 74.4319 0.5169882 0.5305032 -50
15 83.24004 74.4319 0.5169882 0.5305032 -50
16 83.24004 74.4319 0.5169882 0.5305032 -50
17 83.24004 74.4319 0.5169882 0.5305032 -50
CCF_glob_avg CCF_glob_Max_Abs CCF_glob_0_detrended CCF_glob_Max_detrended
1 0.5073075 0.5305032 0.5821183 0.6166339
2 0.5073075 0.5305032 0.5821183 0.6166339
3 0.5073075 0.5305032 0.5821183 0.6166339
4 0.5073075 0.5305032 0.5821183 0.6166339
5 0.5073075 0.5305032 0.5821183 0.6166339
6 0.5073075 0.5305032 0.5821183 0.6166339
7 0.5073075 0.5305032 0.5821183 0.6166339
8 0.5073075 0.5305032 0.5821183 0.6166339
9 0.5073075 0.5305032 0.5821183 0.6166339
10 0.5073075 0.5305032 0.5821183 0.6166339
11 0.5073075 0.5305032 0.5821183 0.6166339
12 0.5073075 0.5305032 0.5821183 0.6166339
13 0.5073075 0.5305032 0.5821183 0.6166339
14 0.5073075 0.5305032 0.5821183 0.6166339
15 0.5073075 0.5305032 0.5821183 0.6166339
16 0.5073075 0.5305032 0.5821183 0.6166339
17 0.5073075 0.5305032 0.5821183 0.6166339
CCF_glob_Max_Lag_detrended CCF_glob_avg_detrended CCF_glob_Max_Abs_detrended
1 -50 0.5703339 0.6166339
2 -50 0.5703339 0.6166339
3 -50 0.5703339 0.6166339
4 -50 0.5703339 0.6166339
5 -50 0.5703339 0.6166339
6 -50 0.5703339 0.6166339
7 -50 0.5703339 0.6166339
8 -50 0.5703339 0.6166339
9 -50 0.5703339 0.6166339
10 -50 0.5703339 0.6166339
11 -50 0.5703339 0.6166339
12 -50 0.5703339 0.6166339
13 -50 0.5703339 0.6166339
14 -50 0.5703339 0.6166339
15 -50 0.5703339 0.6166339
16 -50 0.5703339 0.6166339
17 -50 0.5703339 0.6166339CCF2 <- compute_ccf(data=foo, window_size=30, slide_step=10, max_lag=50)
CCF2 Start_Time End_Time avg_M_Close avg_W_Close CCF_0 CCF_Max
1 0.12 30.12 50.15600 50.14200 0.64504668 0.64504668
2 10.12 40.12 52.36200 58.11600 0.92443544 0.92443544
3 20.12 50.12 58.80400 66.69600 0.76621451 0.76621451
4 30.12 60.12 71.30938 76.62675 0.92018724 0.92018724
5 40.12 70.12 83.93600 82.10800 0.89612944 0.89612944
6 50.12 80.12 93.35800 86.95200 0.81769783 0.81769783
7 60.12 90.12 91.97600 90.10000 -0.25161111 -0.04230794
8 70.12 100.12 86.62600 87.86600 0.49357663 0.58700700
9 80.12 110.12 78.43800 82.53200 0.71251097 0.71251097
10 90.12 120.12 76.68000 75.67800 0.07216259 0.30759085
11 100.12 130.12 79.15200 70.00800 -0.76192983 -0.58858010
12 110.12 140.12 83.60800 66.32400 -0.65149018 -0.30221544
13 120.12 150.12 85.69400 65.47400 0.24307661 0.24307661
14 130.12 160.12 88.34800 74.77000 0.88853928 0.88853928
15 140.12 170.12 91.78400 85.79200 0.93873734 0.93873734
16 150.12 180.12 95.06000 95.27600 0.59377358 0.59501692
17 160.12 190.12 86.57000 94.45400 0.67378310 0.67382561
18 170.12 200.12 75.84400 84.69600 0.67761922 0.83318599
19 180.12 210.12 70.26400 74.56400 -0.02140251 0.31149603
20 190.12 220.12 76.78200 66.33800 -0.10071951 0.11648992
21 200.12 230.12 86.53400 67.07400 0.58283606 0.60420822
22 210.12 240.12 91.87800 70.19800 0.60715489 0.60715489
23 220.12 250.12 93.11800 70.35200 0.48708094 0.50461956
24 230.12 260.12 92.72200 69.44000 0.35606102 0.45126037
25 240.12 270.12 92.80000 66.30000 -0.57017305 -0.46273851
26 250.12 280.12 93.58000 65.97200 -0.67288489 -0.21338831
27 260.12 290.12 94.71400 73.43000 0.56305740 0.66393663
28 270.12 300.12 96.62800 81.88600 0.70596199 0.78528150
29 280.12 310.12 97.25000 89.51200 0.43554881 0.43554881
30 290.12 320.12 95.44000 88.89200 0.66177260 0.66177260
31 300.12 330.12 90.84800 85.58000 0.81124762 0.81124762
32 310.12 340.12 87.51000 77.75000 0.56292909 0.64253992
33 320.12 350.12 86.29800 65.48200 -0.16949470 -0.01587322
CCF_Max_Lag CCF_avg CCF_Max_Abs CCF_0_detrended CCF_Max_detrended
1 0 0.17828102 0.6450467 0.6059374521 0.605937452
2 0 0.79518863 0.9244354 0.7223978738 0.738872887
3 0 0.66933473 0.7662145 -0.5692526406 -0.401270161
4 0 0.78312696 0.9201872 0.3603235790 0.546554429
5 0 0.77650168 0.8961294 -0.0064767531 0.001443623
6 0 0.65828454 0.8176978 0.2789038367 0.347836395
7 -50 -0.29641740 0.4857263 -0.0713887940 -0.026387210
8 -50 0.42236939 0.5870070 -0.4542785425 -0.082742443
9 0 0.60600486 0.7125110 -0.6739303159 -0.037059151
10 -50 0.11407217 0.3075908 0.4086647303 0.473689540
11 -50 -0.68353449 0.7619298 -0.0007614398 0.143622339
12 50 -0.53907191 0.6514902 -0.0148782804 0.092204763
13 0 0.04670751 0.3482530 0.2838223734 0.283822373
14 0 0.79518211 0.8885393 0.7338534024 0.839202879
15 0 0.84370167 0.9387373 0.7302853011 0.814478626
16 1 0.31302077 0.5950169 0.6059238769 0.615196094
17 -2 0.46507664 0.6738256 0.5049700903 0.609498251
18 -50 0.62998534 0.8331860 -0.2367955851 0.335786939
19 -50 -0.08144585 0.5580327 0.5751055269 0.775403262
20 -50 -0.09210565 0.2578701 -0.6068891094 -0.006724633
21 -50 0.49609047 0.6042082 -0.4968499530 0.202811158
22 0 0.53869577 0.6071549 -0.1668403174 -0.083419286
23 50 0.44461711 0.5046196 0.4722671984 0.583099526
24 50 0.37532075 0.4512604 0.2776691166 0.710321825
25 -50 -0.54112619 0.6370538 -0.3673819260 0.072279691
26 50 -0.51658447 0.7237890 -0.7579062261 -0.131597110
27 50 0.44249983 0.6639366 0.2947119785 0.546553441
28 50 0.65125965 0.7852815 -0.1863179627 0.362714723
29 0 0.29823987 0.4355488 0.4364282718 0.436428272
30 0 0.53291113 0.6617726 0.3022391756 0.347799355
31 0 0.70758924 0.8112476 0.2029165898 0.253224437
32 -21 0.47489991 0.6425399 -0.2173358193 -0.011521201
33 -50 -0.16364167 0.3877495 0.3097968334 0.427452199
CCF_Max_Lag_detrended CCF_avg_detrended CCF_Max_Abs_detrended
1 0 0.16648482 0.60593745
2 7 0.60742161 0.73887289
3 50 -0.52182534 0.57253156
4 -16 0.34049086 0.54655443
5 -3 -0.13819514 0.34730309
6 -50 0.15214617 0.34783640
7 -50 -0.15056664 0.22429321
8 -50 -0.37415717 0.46401856
9 50 -0.43954036 0.67538627
10 50 0.35163377 0.47368954
11 -50 0.04417425 0.14362234
12 50 0.03038122 0.09220476
13 0 0.07876288 0.28382237
14 34 0.66697697 0.83920288
15 33 0.64770101 0.81447863
16 7 0.32315199 0.61519609
17 -50 0.36996102 0.60949825
18 -50 -0.23484143 0.76418734
19 -49 0.49967618 0.77540326
20 50 -0.46980192 0.72688924
21 50 -0.24759119 0.49789354
22 50 -0.16586031 0.22463998
23 50 0.43077765 0.58309953
24 50 0.27024188 0.71032183
25 50 -0.31062171 0.58151187
26 50 -0.58192598 0.84954245
27 50 0.26525012 0.54655344
28 50 -0.17611302 0.58643875
29 0 0.29928123 0.43642827
30 -48 0.21954383 0.34779935
31 -47 0.16849211 0.25322444
32 50 -0.13553777 0.22319408
33 -23 0.16060725 0.42745220
avg_glob_M_Close avg_glob_W_Close CCF_glob_0 CCF_glob_Max CCF_glob_Max_Lag
1 83.24004 74.4319 0.5169882 0.5305032 -50
2 83.24004 74.4319 0.5169882 0.5305032 -50
3 83.24004 74.4319 0.5169882 0.5305032 -50
4 83.24004 74.4319 0.5169882 0.5305032 -50
5 83.24004 74.4319 0.5169882 0.5305032 -50
6 83.24004 74.4319 0.5169882 0.5305032 -50
7 83.24004 74.4319 0.5169882 0.5305032 -50
8 83.24004 74.4319 0.5169882 0.5305032 -50
9 83.24004 74.4319 0.5169882 0.5305032 -50
10 83.24004 74.4319 0.5169882 0.5305032 -50
11 83.24004 74.4319 0.5169882 0.5305032 -50
12 83.24004 74.4319 0.5169882 0.5305032 -50
13 83.24004 74.4319 0.5169882 0.5305032 -50
14 83.24004 74.4319 0.5169882 0.5305032 -50
15 83.24004 74.4319 0.5169882 0.5305032 -50
16 83.24004 74.4319 0.5169882 0.5305032 -50
17 83.24004 74.4319 0.5169882 0.5305032 -50
18 83.24004 74.4319 0.5169882 0.5305032 -50
19 83.24004 74.4319 0.5169882 0.5305032 -50
20 83.24004 74.4319 0.5169882 0.5305032 -50
21 83.24004 74.4319 0.5169882 0.5305032 -50
22 83.24004 74.4319 0.5169882 0.5305032 -50
23 83.24004 74.4319 0.5169882 0.5305032 -50
24 83.24004 74.4319 0.5169882 0.5305032 -50
25 83.24004 74.4319 0.5169882 0.5305032 -50
26 83.24004 74.4319 0.5169882 0.5305032 -50
27 83.24004 74.4319 0.5169882 0.5305032 -50
28 83.24004 74.4319 0.5169882 0.5305032 -50
29 83.24004 74.4319 0.5169882 0.5305032 -50
30 83.24004 74.4319 0.5169882 0.5305032 -50
31 83.24004 74.4319 0.5169882 0.5305032 -50
32 83.24004 74.4319 0.5169882 0.5305032 -50
33 83.24004 74.4319 0.5169882 0.5305032 -50
CCF_glob_avg CCF_glob_Max_Abs CCF_glob_0_detrended CCF_glob_Max_detrended
1 0.5073075 0.5305032 0.5821183 0.6166339
2 0.5073075 0.5305032 0.5821183 0.6166339
3 0.5073075 0.5305032 0.5821183 0.6166339
4 0.5073075 0.5305032 0.5821183 0.6166339
5 0.5073075 0.5305032 0.5821183 0.6166339
6 0.5073075 0.5305032 0.5821183 0.6166339
7 0.5073075 0.5305032 0.5821183 0.6166339
8 0.5073075 0.5305032 0.5821183 0.6166339
9 0.5073075 0.5305032 0.5821183 0.6166339
10 0.5073075 0.5305032 0.5821183 0.6166339
11 0.5073075 0.5305032 0.5821183 0.6166339
12 0.5073075 0.5305032 0.5821183 0.6166339
13 0.5073075 0.5305032 0.5821183 0.6166339
14 0.5073075 0.5305032 0.5821183 0.6166339
15 0.5073075 0.5305032 0.5821183 0.6166339
16 0.5073075 0.5305032 0.5821183 0.6166339
17 0.5073075 0.5305032 0.5821183 0.6166339
18 0.5073075 0.5305032 0.5821183 0.6166339
19 0.5073075 0.5305032 0.5821183 0.6166339
20 0.5073075 0.5305032 0.5821183 0.6166339
21 0.5073075 0.5305032 0.5821183 0.6166339
22 0.5073075 0.5305032 0.5821183 0.6166339
23 0.5073075 0.5305032 0.5821183 0.6166339
24 0.5073075 0.5305032 0.5821183 0.6166339
25 0.5073075 0.5305032 0.5821183 0.6166339
26 0.5073075 0.5305032 0.5821183 0.6166339
27 0.5073075 0.5305032 0.5821183 0.6166339
28 0.5073075 0.5305032 0.5821183 0.6166339
29 0.5073075 0.5305032 0.5821183 0.6166339
30 0.5073075 0.5305032 0.5821183 0.6166339
31 0.5073075 0.5305032 0.5821183 0.6166339
32 0.5073075 0.5305032 0.5821183 0.6166339
33 0.5073075 0.5305032 0.5821183 0.6166339
CCF_glob_Max_Lag_detrended CCF_glob_avg_detrended CCF_glob_Max_Abs_detrended
1 -50 0.5703339 0.6166339
2 -50 0.5703339 0.6166339
3 -50 0.5703339 0.6166339
4 -50 0.5703339 0.6166339
5 -50 0.5703339 0.6166339
6 -50 0.5703339 0.6166339
7 -50 0.5703339 0.6166339
8 -50 0.5703339 0.6166339
9 -50 0.5703339 0.6166339
10 -50 0.5703339 0.6166339
11 -50 0.5703339 0.6166339
12 -50 0.5703339 0.6166339
13 -50 0.5703339 0.6166339
14 -50 0.5703339 0.6166339
15 -50 0.5703339 0.6166339
16 -50 0.5703339 0.6166339
17 -50 0.5703339 0.6166339
18 -50 0.5703339 0.6166339
19 -50 0.5703339 0.6166339
20 -50 0.5703339 0.6166339
21 -50 0.5703339 0.6166339
22 -50 0.5703339 0.6166339
23 -50 0.5703339 0.6166339
24 -50 0.5703339 0.6166339
25 -50 0.5703339 0.6166339
26 -50 0.5703339 0.6166339
27 -50 0.5703339 0.6166339
28 -50 0.5703339 0.6166339
29 -50 0.5703339 0.6166339
30 -50 0.5703339 0.6166339
31 -50 0.5703339 0.6166339
32 -50 0.5703339 0.6166339
33 -50 0.5703339 0.6166339Looping across interactions:
# Initialize a list to store results
all_ccf_results <- list()
for (i in unique(temp2$Couple_Int)){
print(i)
datai <- temp2 %>% filter(Couple_Int==i)
foo <- compute_ccf(data=datai, window_size=30, slide_step=10, max_lag=5)
foo <- foo %>% mutate(Couple_Int=datai$Couple_Int[1], Positive=datai$Positive[1])
all_ccf_results[[i]] <- foo
}[1] "110_4"
[1] "110_3"
[1] "110_2"
[1] "110_1"
[1] "109_4"
[1] "109_3"
[1] "109_2"
[1] "109_1"
[1] "107_4"
[1] "107_3"
[1] "107_2"
[1] "107_1"
[1] "106_4"
[1] "106_3"
[1] "106_2"
[1] "106_1"
[1] "105_4"
[1] "105_3"
[1] "105_2"
[1] "105_1"
[1] "104_4"
[1] "104_3"
[1] "104_2"
[1] "104_1"
[1] "103_4"
[1] "103_3"
[1] "103_2"
[1] "103_1"
[1] "102_4"
[1] "102_3"
[1] "102_2"
[1] "102_1"
[1] "101_4"
[1] "101_3"
[1] "101_2"
[1] "101_1"
[1] "100_4"
[1] "100_3"
[1] "100_2"
[1] "100_1"all_CCF <- dplyr::bind_rows(all_ccf_results) %>% relocate(Couple_Int)
head(all_CCF) Couple_Int Start_Time End_Time avg_M_Close avg_W_Close CCF_0
1 110_4 0.12 30.12 55.09000 84.42600 0.29658123
2 110_4 10.12 40.12 60.27000 88.03000 -0.47217721
3 110_4 20.12 50.12 64.79800 64.69600 -0.63436021
4 110_4 30.12 60.12 65.02994 38.72854 0.05732798
5 110_4 40.12 70.12 63.31400 16.30000 -0.16826473
6 110_4 50.12 80.12 62.32600 17.79400 -0.24743372
CCF_Max CCF_Max_Lag CCF_avg CCF_Max_Abs CCF_0_detrended
1 0.30202921 5 0.28779708 0.30202921 -0.5708657
2 -0.39569591 5 -0.44259816 0.47217721 0.2016695
3 -0.60492360 5 -0.62609686 0.63436021 0.7791491
4 0.07389114 5 0.05452046 0.07389114 -0.7009876
5 -0.15799856 5 -0.17417548 0.19965522 -0.2889539
6 -0.22656025 5 -0.24168106 0.24763517 -0.9133842
CCF_Max_detrended CCF_Max_Lag_detrended CCF_avg_detrended
1 -0.5491948 -5 -0.5592911
2 0.2026526 -2 0.1699840
3 0.7791491 0 0.7739165
4 -0.6743289 -5 -0.6913130
5 -0.2805551 5 -0.2858395
6 -0.8764301 -5 -0.9030037
CCF_Max_Abs_detrended avg_glob_M_Close avg_glob_W_Close CCF_glob_0
1 0.5708657 62.639 46.35615 -0.1511389
2 0.2026526 62.639 46.35615 -0.1511389
3 0.7791491 62.639 46.35615 -0.1511389
4 0.7009876 62.639 46.35615 -0.1511389
5 0.2889539 62.639 46.35615 -0.1511389
6 0.9133842 62.639 46.35615 -0.1511389
CCF_glob_Max CCF_glob_Max_Lag CCF_glob_avg CCF_glob_Max_Abs
1 -0.1452664 5 -0.1501755 0.1532556
2 -0.1452664 5 -0.1501755 0.1532556
3 -0.1452664 5 -0.1501755 0.1532556
4 -0.1452664 5 -0.1501755 0.1532556
5 -0.1452664 5 -0.1501755 0.1532556
6 -0.1452664 5 -0.1501755 0.1532556
CCF_glob_0_detrended CCF_glob_Max_detrended CCF_glob_Max_Lag_detrended
1 -0.150714 -0.1449197 5
2 -0.150714 -0.1449197 5
3 -0.150714 -0.1449197 5
4 -0.150714 -0.1449197 5
5 -0.150714 -0.1449197 5
6 -0.150714 -0.1449197 5
CCF_glob_avg_detrended CCF_glob_Max_Abs_detrended Positive
1 -0.1497976 0.152927 0
2 -0.1497976 0.152927 0
3 -0.1497976 0.152927 0
4 -0.1497976 0.152927 0
5 -0.1497976 0.152927 0
6 -0.1497976 0.152927 0GGally::ggpairs(all_CCF %>%select(CCF_0,CCF_Max,CCF_0_detrended,CCF_Max_detrended) , title="Segment-level indices")
GGally::ggpairs(all_CCF %>% distinct(Couple_Int,.keep_all = T) %>%select(CCF_glob_0,CCF_Max,CCF_glob_0_detrended,CCF_glob_Max_detrended) , title="Global indices")
Plots (positive vs. conflict)
# Create the boxplot
# Reshape data to long format for combined plot
all_CCF_long <- all_CCF %>%
mutate(Type = ifelse(Positive == 1, "Pos", "Neg")) %>%
pivot_longer(cols = c(avg_glob_M_Close, avg_glob_W_Close),
names_to = "Variable",
values_to = "Value")
ggplot(all_CCF_long, aes(x = Type, y = Value, fill = Variable)) +
geom_boxplot() +
labs( title = "Level of Warmth", x = "Positive", y = "Value", fill = "Variable") +
theme_minimal()
ggplot(all_CCF %>% mutate(Type=ifelse(Positive==1,"Pos","Neg")), aes(x=Type,y=CCF_glob_0)) +
geom_boxplot() +
labs(title= "Synchrony",x = "Positive", y = "CCF_glob_0") +
theme_minimal()
describeBy(all_CCF %>% select(CCF_glob_0,avg_glob_W_Close,avg_glob_M_Close),all_CCF$Positive)
Descriptive statistics by group
group: 0
vars n mean sd median trimmed mad min max range
CCF_glob_0 1 643 0.11 0.29 0.18 0.13 0.32 -0.56 0.56 1.12
avg_glob_W_Close 2 643 57.96 10.94 55.68 58.44 12.11 32.01 75.65 43.64
avg_glob_M_Close 3 643 65.01 12.53 62.64 63.57 11.68 49.28 92.98 43.71
skew kurtosis se
CCF_glob_0 -0.60 -0.17 0.01
avg_glob_W_Close -0.26 -0.29 0.43
avg_glob_M_Close 0.78 -0.32 0.49
------------------------------------------------------------
group: 1
vars n mean sd median trimmed mad min max range
CCF_glob_0 1 644 0.20 0.26 0.22 0.24 0.24 -0.60 0.52 1.13
avg_glob_W_Close 2 644 74.28 12.43 74.07 75.10 15.29 50.47 92.84 42.37
avg_glob_M_Close 3 644 76.34 12.65 78.36 76.75 16.70 54.86 96.28 41.42
skew kurtosis se
CCF_glob_0 -1.18 1.87 0.01
avg_glob_W_Close -0.49 -0.89 0.49
avg_glob_M_Close -0.15 -1.23 0.50Pseudo couples
This is the main function that creates the combinations of all possible matching
N_Pseudo_Couples <- 1000
UniqueInt <- unique(temp2$Couple_Int)
combinations<-combn(UniqueInt,2)
#printing the first 10
combinations[,1:10] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] "110_4" "110_4" "110_4" "110_4" "110_4" "110_4" "110_4" "110_4" "110_4"
[2,] "110_3" "110_2" "110_1" "109_4" "109_3" "109_2" "109_1" "107_4" "107_3"
[,10]
[1,] "110_4"
[2,] "107_2"Creating the random sync
simulated_data <- list()
for (i in 1:N_Pseudo_Couples){
print(paste0(i, " out of ",N_Pseudo_Couples))
sample.combinations<-sample(length(combinations[1,]),1)#here we randomly sample one combination
a<-temp2 %>% filter(Couple_Int==combinations[1,sample.combinations]) %>% select(Time,M_Close)
b<-temp2 %>% filter(Couple_Int==combinations[2,sample.combinations]) %>% select(Time,W_Close)
a.b <- full_join(a,b) %>% filter(!is.na(M_Close) & !is.na(W_Close))
CCFi <- compute_ccf(data=a.b, window_size=60, slide_step=60, max_lag=5)
Num_simulation <- i
simulated_data[[paste("sim",i,sep="_")]] <- data.frame(Num_simulation,
CCFi)
}
Pseudo_data <- dplyr::bind_rows(simulated_data)
head(Pseudo_data)Creating Sampling distribution of the means (H0)
#Creating sampling distribution of the means for N=the current sample
Sampling_Dist <- list()
for (i in 1:N_Pseudo_Couples){
print(i)
sampled <- sample(1:N_Pseudo_Couples,
length(unique(temp2$Couple_Int)))
sampled_data <- Pseudo_data %>% filter(Num_simulation%in%sampled)
sampled_data_int <- sampled_data %>% distinct(Num_simulation,.keep_all = T)
foo <- data.frame(Iter=i,
Avg_MaxCCF=mean(sampled_data$CCF_Max,na.rm = T),
Glob_Avg_MaxCCF=mean(sampled_data_int$CCF_glob_Max,na.rm = T),
Avg_MaxCCF_detrended=mean(sampled_data$CCF_Max_detrended,na.rm = T),
Glob_Avg_MaxCCF_detrended=mean(sampled_data_int$CCF_glob_Max_detrended,na.rm = T)
)
Sampling_Dist[[paste("sim",i,sep="_")]] <- foo
}
Sampling_Dist <- dplyr::bind_rows(Sampling_Dist)
head(Sampling_Dist)Comparing observed sync mean to pseudo sync means
Segment-level
#hist sampling dist data (randomly combination paired of partners):
(qts<-quantile(Sampling_Dist$Avg_MaxCCF_detrended,probs = c(.025,.975))) 2.5% 97.5%
-0.02488224 0.07553681 hist(Sampling_Dist$Avg_MaxCCF_detrended,breaks = 30,
#xlim= c(0.01,0.5),
, col="grey",xlab="CCF",
main = "CCF H0 Distribution- Segment")
abline(v=qts[1],col="black",lwd=2)
abline(v=qts[2],col="black",lwd=2)
abline(v=mean(all_CCF$CCF_Max_detrended, na.rm=T),col="red",lty=2,lwd=2)
f<-ecdf(Sampling_Dist$Avg_MaxCCF_detrended)
1-(f(mean(all_CCF$CCF_Max_detrended, na.rm=T)))[1] 0.058Interaction level
(qts<-quantile(Sampling_Dist$Avg_MaxCCF_detrended,probs = c(.025,.975))) 2.5% 97.5%
-0.02488224 0.07553681 hist(Sampling_Dist$Avg_MaxCCF_detrended,breaks = 30,
#xlim= c(0.01,0.5),
, col="grey",xlab="CCF",
main = "CCF H0 Distribution- Segment")
abline(v=qts[1],col="black",lwd=2)
abline(v=qts[2],col="black",lwd=2)
abline(v=mean(all_CCF$CCF_Max_detrended, na.rm=T),col="red",lty=2,lwd=2)
f<-ecdf(Sampling_Dist$Avg_MaxCCF_detrended)
1-(f(mean(all_CCF$CCF_Max_detrended, na.rm=T)))[1] 0.058(qts<-quantile(Sampling_Dist$Glob_Avg_MaxCCF_detrended,probs = c(.025,.975))) 2.5% 97.5%
-0.01872703 0.11264097 hist(Sampling_Dist$Avg_MaxCCF_detrended,breaks = 30,
xlim= c(-0.20,0.25),
, col="grey",xlab="CCF",
main = "CCF H0 Distribution- Interaction")
abline(v=qts[1],col="black",lwd=2)
abline(v=qts[2],col="black",lwd=2)
interaction_level <- all_CCF %>% distinct(Couple_Int,.keep_all =T) %>% select(CCF_glob_Max_detrended)
abline(v=mean(interaction_level$CCF_glob_Max_detrended, na.rm=T),col="red",lty=2,lwd=2)
f<-ecdf(Sampling_Dist$Glob_Avg_MaxCCF_detrended)
1-(f(mean(interaction_level$CCF_glob_Max_detrended, na.rm=T)))[1] 0Cross-spectral analysis: Coherence
Provides information about relation between two signals in their frequency domain
1. Simulated signals
Simulate two related time series with most variance in low frequencies
# 1. Generate simulated data
set.seed(123) # Set random seed for reproducibility
minutes <- 60 # Total duration of the simulated time series
samples_per_minute <- 1000 # Sampling rate
n <- minutes * samples_per_minute # Total number of samples
t <- seq(0, minutes, length.out = n) # Time vector in minutesEach sine wave represents a different frequency component
The general form: A * sin(2π * f * t)
A: Amplitude- The amplitudes (2, 1, 0.5, 0.2 in our case) determine the strength or prominence of each frequency component in the overall signal.
f: Frequency (cycles per unit time)
t: Time
2π: Used to convert from cycles to radians (a full circle is 2π radians)
x <- 2*sin(2*pi*0.1*t) + # 0.1 cycles per minute (10-minute cycle); This completes one full cycle every 10 minutes (1/0.1 = 10)
1*sin(2*pi*0.3*t) + # 0.3 cycles per minute (3.33-minute cycle)
0.5*sin(2*pi*0.8*t) + # 0.8 cycles per minute (1.25-minute cycle)
0.2*sin(2*pi*3*t) + # 3 cycles per minute (0.33-minute cycle, small contribution)
rnorm(n, 0, 0.1) # Add small amount of random noise
plot(x)
y is similar to x but with phase shifts
The ‘y’ variable is constructed similarly, but with phase shifts (the + pi/6, + pi/4, etc.).
This phase angle represents how much one wave (or signal) is shifted horizontally compared to another.
These phase shifts create a slight offset between x and y,
Radians: One full cycle of a sine wave corresponds to 2π radians.
These values are fractions of𝜋because 𝜋 radians is equivalent to 180 degrees, which is half of a full cycle
phase_shift = pi/12 # Small phase shift
# phase_shift = pi/2 # Large phase shift
y <- 2*sin(2*pi*0.1*t + phase_shift) +
1*sin(2*pi*0.3*t + phase_shift) +
0.5*sin(2*pi*0.8*t + phase_shift) +
0.2*sin(2*pi*3*t + phase_shift) +
rnorm(n, 0, 0.1)
plot(y)
2. Detrend the time series
Remove any linear trend from the data
x_detrended <- residuals(lm(x ~ t))
y_detrended <- residuals(lm(y ~ t))3. Perform cross-spectral analysis
spans=c(3,3) applies some smoothing to the spectral estimates
spec <- spectrum(cbind(x_detrended, y_detrended), spans = c(3,3), plot = T)
4. Compute average weighted coherence
Setting the cutoff
freq <- spec$freq * samples_per_minute # spec$freq is in cycles per sample
coherence <- spec$coh # Coherence between x and y at each frequency
spec_x <- spec$spec[,1] # Power spectrum of x
spec_y <- spec$spec[,2] # Power spectrum of y
# Find frequencies corresponding to cycles longer than minumum
# cutoff of frequency, 2 cycles per minute
# The value 0.5 represents the period in minutes for one cycle, meaning one complete oscillation takes 0.5 minutes.
# Frequency is the number of cycles per minute, so to convert a period (in minutes per cycle)
# to frequency (cycles per minute), we take the reciprocal of the period: Freq=1/period
freq_cutoff <- 1 / 0.5
low_freq_idx <- which(freq <= freq_cutoff) # Indices of frequencies below cutoff
# Compute weighted coherence
# Weight the coherence by the average power at each frequency
weighted_coh <- coherence[low_freq_idx] *
((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
# Compute average weighted coherence
# This is the final measure of coordination between x and y- standardization
avg_weighted_coh <- sum(weighted_coh) / sum((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
# Print result
cat("Average Weighted Coherence:", avg_weighted_coh, "\n")Average Weighted Coherence: 0.9996846 Making sure that the cutoff was set correctly
# Calculate variance explained by low frequencies
total_var_x <- sum(spec_x) # Total variance in x
total_var_y <- sum(spec_y) # Total variance in y
low_freq_var_x <- sum(spec_x[low_freq_idx]) # Variance in x from low frequencies
low_freq_var_y <- sum(spec_y[low_freq_idx]) # Variance in y from low frequencies
# Print percentage of variance explained by low frequencies
cat("Variance explained by low frequencies (X):", low_freq_var_x / total_var_x * 100, "%\n")Variance explained by low frequencies (X): 98.90183 %cat("Variance explained by low frequencies (Y):", low_freq_var_y / total_var_y * 100, "%\n")Variance explained by low frequencies (Y): 98.90043 %5. Calculate average phase relationship
Extract phase difference from cross-spectrum phase
phase_diff <- spec$phase # Phase difference in radians
phase_fraction <- phase_diff / (2 * pi) # Convert to fraction of a cycle- 2 pi is one full cycle
# Focus on low-frequency components
# low_freq_phase_fraction <- phase_fraction[low_freq_idx]
# Weight the phase_fraction by the average power at each frequency
weighted_low_freq_phase_fraction <- phase_fraction[low_freq_idx] *
((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
# Compute average phase relationship over the low-frequency range
# avg_phase_fraction <- mean(low_freq_phase_fraction)
#standardized
avg_phase_fraction <- sum(weighted_low_freq_phase_fraction) / sum((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
# Print average phase fraction
cat("Average Phase Difference (fraction of cycle):", avg_phase_fraction, "\n")Average Phase Difference (fraction of cycle): -0.04100541 Now with absolute phase diff
abs_weighted_low_freq_phase_fraction <- abs(phase_fraction[low_freq_idx]) *
((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
abs_avg_phase_fraction <- sum(abs_weighted_low_freq_phase_fraction) / sum((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
cat("Average ABS Phase Difference (fraction of cycle):", abs_avg_phase_fraction, "\n")Average ABS Phase Difference (fraction of cycle): 0.04100741 Compute time-domain indices for comparision
ccfi <- ccf(x_detrended,y_detrended,lag.max =50,plot=F)
avg_abs_diff=mean(abs(x_detrended-y_detrended))
trend_x <- cor(t,x);trend_y <- cor(t,y)
ar_x<- ar(x,order.max = 2)$resid;ar_y<- ar(y,order.max = 2)$resid
cor_white <- corr.test(as.numeric(ar_x),as.numeric(ar_y))$rPlots
# Create a 2x2 plot layout
par(mfrow = c(2, 2))
# par(mfrow = c(1, 1))
# Plot the original time series
plot(t, x, type = "l", col = "blue", xlab = "Time (minutes)", ylab = "Value",#ylim = c(0, 100),
main = paste0(paste0(" simulated Time Series \n CCF= ",round(max(ccfi$acf),3))))
lines(t, y, col = "red")
abline(h = 50, col = "black", lty = 2)
# Plot the power spectrum of x (black) and y (red)
plot(freq, spec_x, type = "l", col = "blue",xlab = "Frequency (cycles/min)", ylab = "Power",
main = "Power Spectrum of X(blue) and Y(red)", xlim = c(0, 5))
abline(v = freq_cutoff, col = "red", lty = 2)
lines(freq,spec_y,col="red")
# Plot the coherence
plot(freq, coherence, type = "l", xlab = "Frequency (cycles/min)\n size of red dots represents power strength ", ylab = "Coherence",
main = paste0("Weighted Coherence:", round(avg_weighted_coh,3)), xlim = c(0, 5))
abline(v = freq_cutoff, col = "red", lty = 2) # Add vertical line at cutoff frequency
#adding points representing power strength
#points are drawn larger, highlighting those frequencies where the power is greater.
# Calculate the power
power <- (spec_x + spec_y) / 2
# Scale the point size according to power
# You may need to adjust the scaling factor (e.g., 5) for appropriate visibility
point_size <- power / sum(power) * 7
# Add points with sizes proportional to the power
points(freq, coherence,
pch = 16, col = "red", cex = point_size)
#plot phase_fraction
# Create the plot
plot(freq, phase_fraction, type = "h", lwd = 2,xlim = c(0, 5),
xlab = "Frequency (cycles/min)\n size of red dots represents power strength", ylab = "Phase Difference (fraction of cycle)",
main = paste0("Weighted Phase Difference: ", round(avg_phase_fraction, 3)," \n Abs = ",round(abs_avg_phase_fraction,3)),
col = "blue")
abline(v = freq_cutoff, col = "red", lty = 2) # Add vertical line at cutoff frequency
# Add points with sizes proportional to the power
points(freq, phase_fraction,
pch = 16, col = "red", cex = point_size)
Wrap into a function:
compute_cross_spectral_analysis <- function(data, freq_cutoff, samples_per_minute) {
# Extract x, y, and t from the data
x <- data$x
y <- data$y
t <- data$t
ID <-data$ID[1]
# 2. Detrend the time series
x_detrended <- residuals(lm(x ~ t))
y_detrended <- residuals(lm(y ~ t))
# 3. Perform cross-spectral analysis
spec <- spectrum(cbind(x_detrended, y_detrended), spans = c(3, 3), plot = FALSE)
# Convert frequencies to cycles per minute
freq <- spec$freq * samples_per_minute
coherence <- spec$coh
spec_x <- spec$spec[, 1]
spec_y <- spec$spec[, 2]
# Identify low-frequency indices
low_freq_idx <- which(freq <= freq_cutoff)
# Compute weighted coherence
weighted_coh <- coherence[low_freq_idx] * ((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
avg_weighted_coh <- sum(weighted_coh) / sum((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
# Calculate variance explained by low frequencies
total_var_x <- sum(spec_x)
total_var_y <- sum(spec_y)
low_freq_var_x <- sum(spec_x[low_freq_idx])
low_freq_var_y <- sum(spec_y[low_freq_idx])
# Compute average phase relationship
phase_diff <- spec$phase
phase_fraction <- phase_diff / (2 * pi)
weighted_low_freq_phase_fraction <- phase_fraction[low_freq_idx] * ((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
avg_phase_fraction <- sum(weighted_low_freq_phase_fraction) / sum((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
# Compute absolute average phase difference
abs_weighted_low_freq_phase_fraction <- abs(phase_fraction[low_freq_idx]) * ((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
abs_avg_phase_fraction <- sum(abs_weighted_low_freq_phase_fraction) / sum((spec_x[low_freq_idx] + spec_y[low_freq_idx]) / 2)
# Time-domain indices
ccfi <- ccf(x_detrended, y_detrended, lag.max = 50, plot = FALSE)
avg_abs_diff <- mean(abs(x_detrended - y_detrended))
trend_x <- cor(t, x)
trend_y <- cor(t, y)
ar_x <- ar(x, order.max = 2)$resid
ar_y <- ar(y, order.max = 2)$resid
cor_white <- cor(ar_x, ar_y,use="pairwise")
# Data frame to store results
datai <- data.frame(
avg_weighted_coh = avg_weighted_coh,
avg_phase_fraction = avg_phase_fraction,
cor = round(cor(x_detrended, y_detrended), 3),
ccf = round(max(ccfi$acf), 3),
abs_ccf = round(max(abs(ccfi$acf)), 3),
avg_abs_diff = round(avg_abs_diff, 3),
abs_avg_phase_fraction = round(abs_avg_phase_fraction, 3),
cor_white = cor_white,
trend_x = trend_x,
trend_y = trend_y
)
# Plot results
par(mfrow = c(2, 2))
# Plot the original time series
plot(t, x, type = "l", col = "blue", xlab = "Time (minutes)", ylab = "Value",
main = paste0("Time Series of: ",ID,"\nCCF = ", round(max(ccfi$acf), 3)))
lines(t, y, col = "red")
abline(h = 50, col = "black", lty = 2)
# Plot the power spectrum of x and y
plot(freq, spec_x, type = "l", col = "blue", xlab = "Frequency (cycles/min)", ylab = "Power",
main = "Power Spectrum of X (blue) and Y (red)", xlim = c(0, 5))
abline(v = freq_cutoff, col = "red", lty = 2)
lines(freq, spec_y, col = "red")
# Plot coherence
plot(freq, coherence, type = "l", xlab = "Frequency (cycles/min)\nsize of red dots represents power strength",
ylab = "Coherence", main = paste0("Weighted Coherence:", round(avg_weighted_coh, 3)), xlim = c(0, 5))
abline(v = freq_cutoff, col = "red", lty = 2)
points(freq, coherence, pch = 16, col = "red", cex = (spec_x + spec_y) / sum(spec_x + spec_y) * 7)
# Plot phase fraction
plot(freq, phase_fraction, type = "h", lwd = 2, xlim = c(0, 5),
xlab = "Frequency (cycles/min)\nsize of red dots represents power strength",
ylab = "Phase Difference (fraction of cycle)", col = "blue",
main = paste0("Weighted Phase Difference: ", round(avg_phase_fraction, 3),
"\nAbs = ", round(abs_avg_phase_fraction, 3)))
abline(v = freq_cutoff, col = "red", lty = 2)
points(freq, phase_fraction, pch = 16, col = "red", cex = (spec_x + spec_y) / sum(spec_x + spec_y) * 7)
return(datai)
}With large pahse shift
x <- 2*sin(2*pi*0.1*t) + # 0.1 cycles per minute (10-minute cycle); This completes one full cycle every 10 minutes (1/0.1 = 10)
1*sin(2*pi*0.3*t) + # 0.3 cycles per minute (3.33-minute cycle)
0.5*sin(2*pi*0.8*t) + # 0.8 cycles per minute (1.25-minute cycle)
0.2*sin(2*pi*3*t) + # 3 cycles per minute (0.33-minute cycle, small contribution)
rnorm(n, 0, 0.1) # Add small amount of random noise
# y is similar to x but with phase shifts
# The 'y' variable is constructed similarly, but with phase shifts (the + pi/6, + pi/4, etc.).
#This phase angle represents how much one wave (or signal) is shifted horizontally compared to another.
# These phase shifts create a slight offset between x and y,
# Radians: One full cycle of a sine wave corresponds to 2π radians.
# These values are fractions of𝜋because 𝜋 radians is equivalent to 180 degrees, which is half of a full cycle
phase_shift = pi/2 # Large phase shift
y <- 2*sin(2*pi*0.1*t + phase_shift) +
1*sin(2*pi*0.3*t + phase_shift) +
0.5*sin(2*pi*0.8*t + phase_shift) +
0.2*sin(2*pi*3*t + phase_shift) +
rnorm(n, 0, 0.1)
datai <- data.frame(x=x,
y=y,
t=t,
ID="Simulated")
compute_cross_spectral_analysis(data=datai, freq_cutoff=(1 / 0.5 ), samples_per_minute=1000) 
avg_weighted_coh avg_phase_fraction cor ccf abs_ccf avg_abs_diff
1 0.9916892 -0.249985 0 0.059 0.059 1.947
abs_avg_phase_fraction cor_white trend_x trend_y
1 0.25 0.002291427 -0.1353387 -0.0001717559Different frequencies low phase shift
y2 <- 2*sin(2*pi*0.2*t + phase_shift) +
1*sin(2*pi*0.4*t + phase_shift) +
0.5*sin(2*pi*0.8*t + phase_shift) +
0.2*sin(2*pi*3*t + phase_shift) +
rnorm(n, 0, 0.1)
datai <- data.frame(x=x,
y=y2,
t=t)
compute_cross_spectral_analysis(data=datai, freq_cutoff=(1 / 0.5 ), samples_per_minute=1000) 
avg_weighted_coh avg_phase_fraction cor ccf abs_ccf avg_abs_diff
1 0.1996135 0.1692063 0 0.018 0.018 1.971
abs_avg_phase_fraction cor_white trend_x trend_y
1 0.251 -0.0009013889 -0.1353387 0.0002472748Real data
Couple_Int <- unique(temp2$Couple_Int)
agg_data <- list()
# pdf("G:\\My Drive\\Research\\Conferences & Workshops\\2024-Zurich\\plot1.pdf")
for (i in Couple_Int){
# i <- Couple_Int[1]
print(which(i==Couple_Int))
fooi <- temp2 %>% filter(Couple_Int==i) %>% select(Entry,M_Close,W_Close,Positive)
fooi <- na.omit(fooi)
if (length(fooi[,1])!=0&var(fooi$W_Close)!=0&var(fooi$M_Close)!=0){
t = fooi$Entry; x = fooi$M_Close; y = fooi$W_Close
datai <- data.frame(x=x,
y=y,
t=t,
ID=i)
resulti<- compute_cross_spectral_analysis(data=datai, freq_cutoff=(1 / 0.5 ), samples_per_minute=1000)
agg_data[[paste("data",i,sep="_")]] <- data.frame(i=i,Positive=fooi[1,"Positive"],
resulti)
}else{
print(paste0("problem with data of ",i))
}
}[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
Associations between time-domain and frequency-domain indices
my_df2 <- as.data.frame(dplyr::bind_rows(agg_data)) %>% mutate(abs_cor=abs(cor),
abs2_avg_phase_fraction=abs(avg_phase_fraction))
GGally::ggpairs(my_df2 %>% select(cor,abs_cor,ccf,abs_ccf,avg_weighted_coh,avg_phase_fraction,avg_abs_diff))