ML Temporal nets CZ code.R

library(tidyverse)
library(lubridate)
library(here)

#load data
load(here("data", "adult_inters_mx.RData")) #load basic adult interaction data, symmetrical list of edgelists, and data formatted into muxViz style
#all generated by "Social inters to muxviz Functions" code

# Reducibility Analysis ---------------------------------------------------
source("muxLib DF.R") #set source to the location of the "muxLib DF.R" file
#note this is the original with the muxViz master file except we made some edits with "Re" for the calculation of eigenvalues to only obtain real numbers. 
#See L 413, 431, & 461 of "muxLib DF.R". This solution was provided by Nitika Sharma

# Build supra-adjacency matrices
build_supra_AM <- function(edgelist){
  BuildSupraAdjacencyMatrixFromExtendedEdgelist(
    mEdges = as.data.frame(edgelist[,1:5]),
    Layers = length(unique(edgelist$layer1)),
    Nodes = length(unique(edgelist$node1)), isDirected=F)
}

adult_inters_SAM = lapply(adult_inters_mx, build_supra_AM) # apply the function to build the supra AM's


#Colony temporal reducibility to attack speed####

adult_inters_reduce = list() # create an empty list to hold the reduced data
#takes a minute or so. Have added a "tryCatch" function to keep working if an error comes up
for (i in 1:length(unique(adult_inters$Colony))) {
  tryCatch({
    adult_inters_reduce[[i]] = GetMultilayerReducibility(
      as(adult_inters_SAM[[i]], "matrix"), # David Fisher's hack to make this run without errors
      Layers = length(unique(adult_inters_mx[[i]]$layer1)),
      Nodes = length(unique(adult_inters_mx[[i]]$node1)),
      Method="single",
      Type = "Categorical")
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
} # we get debug messages but the code does run. 

lapply(adult_inters_reduce, '[[', "gQualityFunction") # pull out the gQualityFunction element of each
#10 and 24 show peaks not at max layers

# Create some functions for plotting
plot_ML_reduce <- function(reduce_list) { # for line plot
  points(1:19, reduce_list$gQualityFunction, type="l")}

plot_ML_reduce2 = function(reduce_list) { # for dashed line plot
  points(1:19, reduce_list$gQualityFunction, type="l", lty=2)}

# Make the initial plot (just the axes and axis labels)
plot(1, type = "n", xlab = "Amount of aggregation", 
     ylab = "Difference in relative entropy",
     xlim = c(0, 20), ylim = c(0, 14), frame = F, las = 1, yaxs = "i", xaxs = "i")

# Add lines for the interactions besides 10 and 24, i.e. for the ones that don't peak somewhere besides at the max # of layers.
lapply(adult_inters_reduce[c(-10,-24)], plot_ML_reduce) #2 weirdo (10 & 24 colony, BR3-5 & BR10-9), rest OK

# Add dashed lines for 10 and 24 because they peak somewhere besides at the max # of layers
lapply(adult_inters_reduce[c(10,24)], plot_ML_reduce2)

# Make a data frame to contain the data about reduction in information
uniqueColonies <- unique(adult_inters$Colony)
nLayers <- length(unique(adult_inters$layer1))
col_info_data <- data.frame(colony = rep(uniqueColonies, each = nLayers),
                            layer = 1:nLayers,
                            information = unlist(lapply(adult_inters_reduce, 
                                                        "[[", "gQualityFunction"))) %>%
  mutate(colony = factor(colony))

# Summarize
col_info_summ = col_info_data %>% 
  group_by(colony) %>%
  summarise(max_info = as.numeric(max(information))) %>%
  mutate(colony = unique(adult_inters$Colony))

# Load adult prey attack speed data
# Download from figshare
download.file("https://ndownloader.figshare.com/files/21743820", 
              destfile = here("data", "adult_prey_raw.csv"))
adult_prey <- read.csv(here("data", "adult_prey_raw.csv"), header = T)

# fix the column names (replace full stops with underscores)
names(adult_prey) <- gsub("\\.", "_", names(adult_prey)) 
names(adult_prey)[names(adult_prey) == "latancy_to_attack"] <- "latency_to_attack"

#shape into useful format & calculate mean & sd of latency to attack
# colony collective behavior
col_collect_beh <- adult_prey %>% # 
  mutate(latency_to_attack = as.numeric(hms(latency_to_attack))) %>%
  filter(!is.na(latency_to_attack)) %>%
  group_by(Colony) %>%
  summarise(m_latency_att = mean(latency_to_attack),
            sd_latency_att = sd(latency_to_attack)) %>%
  mutate(cv_latency_att = sd_latency_att/m_latency_att) %>%
  select(colony=Colony,  m_latency_att, sd_latency_att, cv_latency_att)

# add the max info information to the summary stats
attack_to_reduce_data <- col_info_summ %>%
  left_join(col_collect_beh)

# Prepare to plot: latency to attack vs. variability of network
par(mfrow = c(1,2), mar = c(5,6,2,2), xpd = NA)

with(attack_to_reduce_data, 
     plot(max_info, m_latency_att, pch = 16, frame = F,
          ylim = c(0,500), xlim = c(5,15),
          xlab = "Variability of network", ylab = "Mean latency to attack",
          las = 1))
text(2,600,"a",font=2) #whether this is in the correct place will depend on the size of your plotting window. 
#If you cannot see the "a" try reducing the height.
box(which = "plot", bty = "l")
with(attack_to_reduce_data, cor.test(max_info, m_latency_att, method="spearman")) #no relationship

#par(mar=c(5,6,1,1), xpd=NA)

with(attack_to_reduce_data, 
     plot(max_info, cv_latency_att, 
          pch = 16, frame = F, ylim = c(0,1.5), xlim = c(5,15), 
          xlab = "Variability of network", 
          ylab ="Coefficient of variation\nin latency to attack",
          las = 1))
text(2,1.8,"b",font=2)
box(which = "plot", bty = "l")
with(attack_to_reduce_data, cor.test(max_info, cv_latency_att, method="spearman")) #variance in max info scores not related to latency to attack

# Attack speed/variation vs.  keystone ------------------------------------
#Does attack speed/variation relate to number of times the same individual was keystone?####
#keystone is defined as having the highest boldness in the colony that week
#We don't allow anyone to be keystone if all indivs did not react (got 0s)

download.file("https://ndownloader.figshare.com/files/21743826", 
              destfile = here("data", "adult_bold_raw.csv"))
adult_bold_raw = read.csv(here("data", "adult_bold_raw.csv"), header=T)

adult_bold = adult_bold_raw %>%
  mutate(spider = factor(paste0(adult_bold_raw$Individual, "_",
                                adult_bold_raw$colony))) %>%
  gather(week, shyness, X1:X7) %>%
  mutate(boldness = 600-shyness) %>%
  select(spider, colony, treatment, week, boldness)

adult_keystone = adult_bold %>%
  group_by(colony, week) %>%
  filter(boldness == max(boldness, na.rm=T) , boldness > 0) %>%
  #note we retain 2 equally bold spiders so they both count as "keystone" for that week
  group_by(colony) %>%
  count(spider) %>% #count how often each ID appears per colony
  filter(n == max(n)) %>% #only keep largest
  top_n(-1,spider) #if there are two that were keystone for the same number of weeks, then pick one arbitrarily based on ID

attack_to_keys_data = col_collect_beh %>%
  left_join(adult_keystone) 

par(mar=c(5,5,2,2),mfrow=c(1,2), xpd=NA)
#on mean attack speed
with(attack_to_keys_data, boxplot(m_latency_att ~ factor(n), frame=F, ylim=c(0,500),las=1,
                                  xlab="Number of weeks the same individual was keystone",
                                  ylab="Mean latency to attack"))
text(-0.5,570,"a",font=2)
box(which = "plot", bty = "l")
with(attack_to_keys_data, kruskal.test(m_latency_att ~ factor(n))) #no

#on CV of attack speed
with(attack_to_keys_data, boxplot(cv_latency_att ~ factor(n), frame=F,las=1,
                                  xlab="Number of weeks the same individual was keystone",
                                  ylab="Coefficient of variation\nin latency to attack")) 
text(-0.5,1.66,"b",font=2)
box(which = "plot", bty = "l")
with(attack_to_keys_data, kruskal.test(cv_latency_att ~ factor(n))) #not diff


# Multilayer connectivity and indiv. boldness -----------------------------
#Connectivity in ML net to individual boldness####

#calculate some ML eigenvector centrality using muxViz function "GetMultiEigenvectorCentrality" within a loop for each group
adult_inters_eigenvec = list()
#Have added a "tryCatch" function to keep working if an error comes up
for (i in 1:length(unique(adult_inters$Colony))) {
  tryCatch({
    adult_inters_eigenvec[[i]] = GetMultiEigenvectorCentrality(
      adult_inters_SAM[[i]],
      Layers = length(unique(adult_inters_mx[[i]]$layer1)),
      Nodes = length(unique(adult_inters_mx[[i]]$node1)))
  }, error = function(e){cat("ERROR :",conditionMessage(e), "\n")})
}

#Also want some measure of connectivity based on non-ML approach, e.g. mean/median and sd/cv of degree across the 19 time points

get_degree <- function(edgelist) {
  edgelist %>% group_by(node1, layer1) %>%
    summarise(degree = length(unique(node2))-1) #-1 is to not count the interlayer edges to itself for degree
}

adult_inters_degree <- lapply(adult_inters_all, get_degree) #note using non-muxViz version of edgelist as that had kept the original spiders IDs

summarise_degree <- function(degree_list) {
  degree_list %>% group_by(node1) %>%
    summarise(mean_degree = mean(degree),
              median_degree = median(degree),
              sd_degree = sd(degree),
              cv_degree = sd_degree/mean_degree)
}

adult_inters_sumdegree <- lapply(adult_inters_degree, summarise_degree) %>%
  bind_rows( .id = "column_label") %>%
  mutate(evc = unlist(adult_inters_eigenvec)) %>%
  separate(node1, into=c(NA, "colony"), sep="_", remove=F) %>%
  group_by(colony) %>%
  mutate(col_m_mdeg = mean(mean_degree, na.rm=T),
         col_m_cvdeg = mean(cv_degree, na.rm=T),
         col_m_evc = mean(evc, na.rm=T),
         #create delta mean degree etc to allow more sensible comaprisons
         d_mdeg = mean_degree-col_m_mdeg,
         d_cvdeg = cv_degree-col_m_cvdeg,
         d_evc = evc - col_m_evc)

#Plots and randomisation test of correlations among centrality measures

#Mean and CV of degree
par(xpd = NA)
with(adult_inters_sumdegree, plot(d_cvdeg, d_mdeg, las = 1,
                                  pch=16, frame = F, 
                                  xlab = "Relative coefficient of variation in degree",
                                  ylab = "Relative mean degree"))
text(-1.4,1.45,"a",font=2)
box(which = "plot", bty = "l")
par(xpd = F)
abline(lm(d_mdeg~d_cvdeg, data = adult_inters_sumdegree))
with(adult_inters_sumdegree, cor.test(d_cvdeg, d_mdeg, method="spearman"))
#strong negative corr, -0.806

dcvdeg_dmdeg_obv = cor.test(adult_inters_sumdegree$d_cvdeg, adult_inters_sumdegree$d_mdeg, method="spearman")$estimate #extract observed value
#then shuffle degree scores and calculate correlation 1000 times
dcvdeg_dmdeg_r = numeric(1000)
for (i in 1:1000) {
  test = adult_inters_sumdegree %>% group_by(colony) %>%
    mutate(rand_d_cvdeg = sample(d_cvdeg, replace=F))
  dcvdeg_dmdeg_r[i] = cor.test(test$rand_d_cvdeg, test$d_mdeg, method="spearman")$estimate
}
#calculate p value
min(sum(dcvdeg_dmdeg_obv>dcvdeg_dmdeg_r),sum(dcvdeg_dmdeg_obv<dcvdeg_dmdeg_r)) /500 #0

#Mean degree and EVC
par(xpd = NA)
with(adult_inters_sumdegree, plot(d_mdeg, d_evc,  pch=16, las = 1,
                                  frame = F, ylim=c(-3,2),
                                  xlab = "Relative mean degree", 
                                  ylab = "Relative eigenvector centrality"))
text(-3,2.8,"b",font=2)
box(which = "plot", bty = "l")
par(xpd = F)
abline(lm(d_evc~d_mdeg, data = adult_inters_sumdegree))
with(adult_inters_sumdegree, cor.test(d_mdeg, d_evc, method="spearman"))


dmdeg_devc_obv = cor.test(adult_inters_sumdegree$d_mdeg, adult_inters_sumdegree$d_evc, method="spearman")$estimate 
dmdeg_devc_r = numeric(1000)
for (i in 1:1000) {
  test = adult_inters_sumdegree %>% group_by(colony) %>%
    mutate(rand_d_mdeg = sample(d_mdeg, replace=F))
  dmdeg_devc_r[i] = cor.test(test$rand_d_mdeg, test$d_evc, method="spearman")$estimate
}

#p value
min(sum(dmdeg_devc_obv>dmdeg_devc_r),sum(dmdeg_devc_obv<dmdeg_devc_r)) /500 #0

#CV of degree and EVC
par(xpd = NA)
with(adult_inters_sumdegree, plot(d_cvdeg, d_evc,  pch=16, las=1,
                                  frame = F,  ylim=c(-2,2), 
                                  xlab = "Relative coefficient of variation in degree",
                                  ylab = "Relative eigenvector centrality"))
text(-1.4,2.7,"c",font=2)
box(which = "plot", bty = "l")
with(adult_inters_sumdegree, cor.test(d_cvdeg, d_evc, method="spearman"))

dcvdeg_devc_obv = cor.test(adult_inters_sumdegree$d_cvdeg, adult_inters_sumdegree$d_evc, method="spearman")$estimate
dcvdeg_devc_r = numeric(1000)
for (i in 1:1000) {
  test = adult_inters_sumdegree %>% group_by(colony) %>%
    mutate(rand_d_cvdeg = sample(d_cvdeg, replace=F))
  dcvdeg_devc_r[i] = cor.test(test$rand_d_cvdeg, test$d_evc, method="spearman")$estimate
}

#p value
min(sum(dcvdeg_devc_obv>dcvdeg_devc_r),sum(dcvdeg_devc_obv<dcvdeg_devc_r)) /500 #0.080

#Plot the distributions of the randomised correlation coefs for the supp mat
par(mfrow=c(1,3), xpd=NA)
hist(dmdeg_devc_r, xlim=c(-0.2, 0.5), las=1, 
     xlab = "Spearman correlation coefficient",
     main=NULL, breaks=15)
text(-0.45,135,"a",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dmdeg_devc_obv, col="red")

par(xpd=NA)
hist(dcvdeg_devc_r, xlim=c(-0.2,0.2),las=1,
     xlab = "Spearman correlation coefficient", 
     main=NULL, breaks = 15)
text(-0.35,150,"b",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dcvdeg_devc_obv, col="red")

par(xpd=NA)
hist(dcvdeg_dmdeg_r, xlim=c(-0.8,0.2),las=1,
     xlab = "Spearman correlation coefficient",
     main=NULL, breaks=15)
text(-1.2,277,"c",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dcvdeg_dmdeg_obv, col="red")


# Boldness and centrality -------------------------------------------------
#Boldness and centrality measures####
#Have 3 measures of an individuals network position, eigenvector centrality, mean degree, and cv of degree.
#see how each of these 3 relate to mean boldness and cv of boldness

adult_sumbold = adult_bold %>%
  group_by(spider) %>%
  summarise(m_bold = mean(boldness, na.rm=T),
            sd_bold = sd(boldness, na.rm=T),
            cv_bold = sd_bold/m_bold)

adult_bold_inters = adult_sumbold %>%
  left_join(adult_inters_sumdegree, by=c("spider" = "node1" )) %>%
  filter(!is.na(colony)) %>%
  group_by(colony) %>%
  mutate(col_m_mbold = mean(m_bold, na.rm=T),
         col_m_cvbold = mean(cv_bold, na.rm=T),
         d_mbold = m_bold-col_m_mbold,
         d_cvbold = cv_bold - col_m_cvbold)

#Plotting and conducting randomisation tests
#Note the exact p values for all randomisation tests will be different to the published version (and each tie it is run) as the randomisations are run anew each time
#the overall conclusions should not change however.

par(mar=c(5,6,2,1),mfrow=c(3,2), xpd = NA)

#EVC
with(adult_bold_inters, plot(d_mbold, d_evc,  pch=16, frame=F, las = 1,
                             xlim=c(-200,300),ylim=c(-2,2),
                             xlab="Relative mean boldness", 
                             ylab="Relative eigenvector centrality")) 
text(-300,2.5,"a",font=2)
box(which = "plot", bty = "l")
with(adult_bold_inters, cor.test(d_mbold, d_evc, method="spearman")) 

dmbold_devc_obv = cor.test(adult_bold_inters$d_mbold, adult_bold_inters$d_evc, method="spearman")$estimate
dmbold_devc_r = numeric(1000)
for (i in 1:1000) {
  test = adult_bold_inters %>% group_by(colony) %>%
    mutate(rand_d_mbold = sample(d_mbold, replace=F))
  dmbold_devc_r[i] = cor.test(test$rand_d_mbold, test$d_evc, method="spearman")$estimate
}

min(sum(dmbold_devc_obv>dmbold_devc_r),sum(dmbold_devc_obv<dmbold_devc_r)) /500 


with(adult_bold_inters, plot(d_cvbold, d_evc, pch=16, frame=F,las = 1,
                             ylim=c(-3,2), xlim=c(-1.5,1.5),
                             xlab="Relative coefficient of variation in boldness", 
                             ylab="Relative eigenvector centrality")) 
text(-2.1,2.8,"b",font=2)
box(which = "plot", bty = "l")

with(adult_bold_inters, cor.test(d_cvbold, d_evc, method="spearman")) 

dcvbold_devc_obv = cor.test(adult_bold_inters$d_cvbold, adult_bold_inters$d_evc, method="spearman")$estimate
dcvbold_devc_r = numeric(1000)
for (i in 1:1000) {
  test = adult_bold_inters %>% group_by(colony) %>%
    mutate(rand_d_cvbold = sample(d_cvbold, replace=F))
  dcvbold_devc_r[i] = cor.test(test$rand_d_cvbold, test$d_evc, method="spearman")$estimate
}

min(sum(dcvbold_devc_obv>dcvbold_devc_r),sum(dcvbold_devc_obv<dcvbold_devc_r)) /500 #0.866

#Mean degree
with(adult_bold_inters, plot(d_mbold, d_mdeg, pch=16, frame=F,las = 1,
                             ylim=c(-3,1), xlim=c(-200,300), 
                             xlab="Relative mean boldness", 
                             ylab="Relative mean degree"))
text(-300,1.5,"c",font=2)
box(which = "plot", bty = "l")
par(xpd = F)
abline(lm(d_mdeg~d_mbold, data = adult_bold_inters))
with(adult_bold_inters, cor.test(d_mbold, d_mdeg, method="spearman"))

dmbold_dmdeg_obv = cor.test(adult_bold_inters$d_mbold, adult_bold_inters$d_mdeg, method="spearman")$estimate
dmbold_dmdeg_r = numeric(1000)
for (i in 1:1000) {
  test = adult_bold_inters %>% group_by(colony) %>%
    mutate(rand_d_mbold = sample(d_mbold, replace=F))
  dmbold_dmdeg_r[i] = cor.test(test$rand_d_mbold, test$d_mdeg, method="spearman")$estimate
}

min(sum(dmbold_dmdeg_obv>dmbold_dmdeg_r),sum(dmbold_dmdeg_obv<dmbold_dmdeg_r)) /500 

with(adult_bold_inters, plot(d_cvbold, d_mdeg, pch=16, frame=F,las = 1,
                             ylim=c(-2,1), xlim=c(-1.5,1.5), 
                             xlab="Relative coefficient of variation in boldness", 
                             ylab="Relative mean degree"))
par(xpd = NA)
text(-2.1,1.5,"d",font=2)
box(which = "plot", bty = "l")

with(adult_bold_inters, cor.test(d_cvbold, d_mdeg, method="spearman")) 

dcvbold_dmdeg_obv = cor.test(adult_bold_inters$d_cvbold, adult_bold_inters$d_mdeg, method="spearman")$estimate
dcvbold_dmdeg_r = numeric(1000)
for (i in 1:1000) {
  test = adult_bold_inters %>% group_by(colony) %>%
    mutate(rand_d_cvbold = sample(d_cvbold, replace=F))
  dcvbold_dmdeg_r[i] = cor.test(test$rand_d_cvbold, test$d_mdeg, method="spearman")$estimate
}

min(sum(dcvbold_dmdeg_obv>dcvbold_dmdeg_r),sum(dcvbold_dmdeg_obv<dcvbold_dmdeg_r)) /500 #0.140

#CV of degree
par(xpd = NA)
with(adult_bold_inters, plot(d_mbold, d_cvdeg, pch=16, frame=F,las = 1,
                             ylim=c(-1,2.5), xlim=c(-200,300), 
                             xlab="Relative mean boldness", 
                             ylab="Relative coefficient\nof variation in degree"))
text(-300,3,"e",font=2)
box(which = "plot", bty = "l")
par(xpd = F)
abline(lm(d_cvdeg~d_mbold, data = adult_bold_inters))
with(adult_bold_inters, cor.test(d_mbold, d_cvdeg, method="spearman")) 

dmbold_dcvdeg_obv = cor.test(adult_bold_inters$d_mbold, adult_bold_inters$d_cvdeg, method="spearman")$estimate
dmbold_dcvdeg_r = numeric(1000)
for (i in 1:1000) {
  test = adult_bold_inters %>% group_by(colony) %>%
    mutate(rand_d_mbold = sample(d_mbold, replace=F))
  dmbold_dcvdeg_r[i] = cor.test(test$rand_d_mbold, test$d_cvdeg, method="spearman")$estimate
}

min(sum(dmbold_dcvdeg_obv>dmbold_dcvdeg_r),sum(dmbold_dcvdeg_obv<dmbold_dcvdeg_r)) /500 

par(xpd = NA)
with(adult_bold_inters, plot(d_cvbold, d_cvdeg, pch=16, frame=F,las = 1,
                             ylim=c(-1,1.5), xlim=c(-1.5,1.5), 
                             xlab="Relative coefficient of variation in boldness", 
                             ylab="Relative coefficient\nof variation in degree"))
text(-2.1,2,"f",font=2)
box(which = "plot", bty = "l")
with(adult_bold_inters, cor.test(d_cvbold, d_cvdeg, method="spearman"))

dcvbold_dcvdeg_obv = cor.test(adult_bold_inters$d_cvbold, adult_bold_inters$d_cvdeg, method="spearman")$estimate
dcvbold_dcvdeg_r = numeric(1000)
for (i in 1:1000) {
  test = adult_bold_inters %>% group_by(colony) %>%
    mutate(rand_d_cvbold = sample(d_cvbold, replace=F))
  dcvbold_dcvdeg_r[i] = cor.test(test$rand_d_cvbold, test$d_cvdeg, method="spearman")$estimate
}

min(sum(dcvbold_dcvdeg_obv>dcvbold_dcvdeg_r),sum(dcvbold_dcvdeg_obv<dcvbold_dcvdeg_r)) /500 


#Plot distributions of correlation coefficients for the supp mat
par(mfrow=c(3,2), xpd=NA)

hist(dmbold_devc_r, xlim=c(-0.3,0.3), las = 1,
     xlab = "Spearman correlation coefficient",
     main=NULL,breaks=20)
text(-0.45,140,"a",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dmbold_devc_obv, col="red")

par(xpd=NA)
hist(dcvbold_devc_r,xlim=c(-0.3,0.3),las = 1,
     xlab = "Spearman correlation coefficient",
     main=NULL, breaks=20)
text(-0.45,133,"b",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dcvbold_devc_obv, col="red")

par(xpd=NA)
hist(dmbold_dmdeg_r, xlim=c(-0.2,0.2),las = 1,
     xlab = "Spearman correlation coefficient",
     main=NULL, breaks=20)
text(-0.3,140,"c",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dmbold_dmdeg_obv, col="red")

par(xpd=NA)
hist(dcvbold_dmdeg_r, xlim=c(-0.3,0.3),las = 1,
     xlab = "Spearman correlation coefficient",
     main=NULL, breaks=25)
text(-0.47,110,"d",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dcvbold_dmdeg_obv, col="red")

par(xpd=NA)
hist(dmbold_dcvdeg_r,xlim=c(-0.2,0.2),las = 1,
     xlab = "Spearman correlation coefficient",
     main=NULL, breaks=20)
text(-0.3,145,"e",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dmbold_dcvdeg_obv, col="red")

par(xpd=NA)
hist(dcvbold_dcvdeg_r,xlim=c(-0.3,0.3),las = 1,
     xlab = "Spearman correlation coefficient",
     main=NULL, breaks=20)
text(-0.45,120,"f",font=2)
box(which = "plot", bty = "l")

par(xpd=F)
abline(v = dcvbold_dcvdeg_obv, col="red")


####END####