R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

rm(list = ls())
gc()

          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  585283 31.3    1238594 66.2  1238594 66.2
Vcells 1457144 11.2    8388608 64.0  7583308 57.9

# si no andan los paths, setear esa siguiente linea
#setwd("acá va el path hasta el directorio 'neuro'")

library(brainGraph)

Loading required package: igraph

Attaching package: 㤼㸱igraph㤼㸲

The following objects are masked from 㤼㸱package:stats㤼㸲:

    decompose, spectrum

The following object is masked from 㤼㸱package:base㤼㸲:

    union

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

library(igraph)
library(ggplot2)
library(tidyr)

Attaching package: 㤼㸱tidyr㤼㸲

The following object is masked from 㤼㸱package:igraph㤼㸲:

    crossing

source(paste0(getwd(),"/red_por_densidad.R"))

# se generan todos los grafos individuales para cada sujeto en cada estadio de sue?o para cada delta

seq <- seq(from = 0.025, to = 0.15, by = 0.005)
suj.count <- 18

N1 <- array(dim = length(seq))
N2 <- array(dim = length(seq))
N3 <- array(dim = length(seq))
W <- array(dim = length(seq))

for(i in 1:length(seq)) {
  tmpN1 <- array(dim = suj.count)
  tmpN2 <- array(dim = suj.count)
  tmpN3 <- array(dim = suj.count)
  tmpW <- array(dim = suj.count)
  for (j in 1:18) {
    tmp <- read.table( paste0(getwd(),paste0("/datasets/DataSujetos/N1_suj",j,".csv")), header=F, sep = ",")
    result <- red.por.densidad(tmp, seq[i])
    tmpN1[j] <- list(grafo = result$red)
    
    tmp <- read.table( paste0(getwd(),paste0("/datasets/DataSujetos/N2_suj",j,".csv")), header=F, sep = ",")
    result <- red.por.densidad(tmp, seq[i])
    tmpN2[j] <- list(grafo = result$red)
    
    tmp <- read.table( paste0(getwd(),paste0("/datasets/DataSujetos/N3_suj",j,".csv")), header=F, sep = ",")
    result <- red.por.densidad(tmp, seq[i])
    tmpN3[j] <- list(grafo = result$red)

    tmp <- read.table( paste0(getwd(),paste0("/datasets/DataSujetos/W_suj",j,".csv")), header=F, sep = ",")
    result <- red.por.densidad(tmp, seq[i])
    tmpW[j] <- list(grafo = result$red)

  }
  N1[i] <- list(tmpN1)
  N2[i] <- list(tmpN2)
  N3[i] <- list(tmpN3)
  W[i] <- list(tmpW)
}

[1] "Warning!, density 0.025 aproximada con error -0.000187406296851572"
[1] "Warning!, density 0.025 aproximada con error -0.000187406296851572"
[1] "Warning!, density 0.03 aproximada con error 0.000164917541229388"
[1] "Warning!, density 0.03 aproximada con error 0.000164917541229388"
[1] "Warning!, density 0.04 aproximada con error 0.000119940029985006"
[1] "Warning!, density 0.04 aproximada con error -0.000179910044977506"
[1] "Warning!, density 0.055 aproximada con error -0.000172413793103447"
[1] "Warning!, density 0.06 aproximada con error 0.000179910044977513"
[1] "Warning!, density 0.07 aproximada con error -0.000164917541229381"
[1] "Warning!, density 0.07 aproximada con error 0.000134932533733131"
[1] "Warning!, density 0.07 aproximada con error 0.000134932533733131"
[1] "Warning!, density 0.075 aproximada con error 0.000187406296851586"
[1] "Warning!, density 0.075 aproximada con error -0.000112443778110927"
[1] "Warning!, density 0.085 aproximada con error 0.00014242878560719"
[1] "Warning!, density 0.085 aproximada con error -0.000157421289355322"
[1] "Warning!, density 0.085 aproximada con error -0.000157421289355322"
[1] "Warning!, density 0.085 aproximada con error -0.000157421289355322"
[1] "Warning!, density 0.085 aproximada con error 0.00014242878560719"
[1] "Warning!, density 0.085 aproximada con error 0.000292353823088454"
[1] "Warning!, density 0.085 aproximada con error -0.000157421289355322"
[1] "Warning!, density 0.09 aproximada con error 0.000194902548725631"
[1] "Warning!, density 0.1 aproximada con error 0.000149925037481263"
[1] "Warning!, density 0.1 aproximada con error 0.000149925037481263"
[1] "Warning!, density 0.1 aproximada con error 0.000149925037481263"
[1] "Warning!, density 0.1 aproximada con error -0.000149925037481263"
[1] "Warning!, density 0.11 aproximada con error 0.000104947526236895"
[1] "Warning!, density 0.11 aproximada con error 0.000104947526236895"
[1] "Warning!, density 0.125 aproximada con error 0.000112443778110941"
[1] "Warning!, density 0.125 aproximada con error 0.000112443778110941"
[1] "Warning!, density 0.125 aproximada con error -0.000337331334332835"
[1] "Warning!, density 0.125 aproximada con error 0.000112443778110941"
[1] "Warning!, density 0.13 aproximada con error -0.000284857571214381"
[1] "Warning!, density 0.13 aproximada con error -0.000134932533733118"
[1] "Warning!, density 0.13 aproximada con error -0.000134932533733118"
[1] "Warning!, density 0.14 aproximada con error -0.000179910044977499"
[1] "Warning!, density 0.14 aproximada con error 0.000119940029985027"
[1] "Warning!, density 0.14 aproximada con error 0.000119940029985027"
[1] "Warning!, density 0.14 aproximada con error -0.000179910044977499"
[1] "Warning!, density 0.14 aproximada con error -0.000179910044977499"
[1] "Warning!, density 0.14 aproximada con error 0.000119940029985027"
[1] "Warning!, density 0.145 aproximada con error -0.000127436281859072"
[1] "Warning!, density 0.145 aproximada con error 0.000172413793103454"

louvain_n1<-louvain_n2<-louvain_n3<-louvain_w<-list()
degree_n1<-degree_n2<-degree_n3<-degree_w<-list()
z_score_n1<-z_score_n2<-z_score_n3<-z_score_w<-list()
p_factor_n1<-p_factor_n2<-p_factor_n3<-p_factor_w<-list()



for(i in 1:length(seq)){ # itero sobre los deltas
  
  louvain_n1[[i]] <- louvain_n2[[i]] <- louvain_n3[[i]] <- louvain_w[[i]] <- list()
  degree_n1[[i]] <- degree_n2[[i]] <- degree_n3[[i]] <- degree_w[[i]] <-  list()
  z_score_n1[[i]] <- z_score_n2[[i]] <- z_score_n3[[i]] <- z_score_w[[i]] <- list()
  p_factor_n1[[i]] <- p_factor_n2[[i]] <- p_factor_n3[[i]] <- p_factor_w[[i]]<-list()

  for (j in 1:suj.count) {#itero sobre los sujetos
    louvain_n1[[i]][[j]] <- cluster_louvain(graph = N1[[i]][[j]], weights = NULL)
    louvain_n2[[i]][[j]] <- cluster_louvain(graph = N2[[i]][[j]], weights = NULL)
    louvain_n3[[i]][[j]] <- cluster_louvain(graph = N3[[i]][[j]], weights = NULL)
    louvain_w[[i]][[j]] <- cluster_louvain(graph = W[[i]][[j]], weights = NULL)
  
    degree_n1[[i]][[j]] <- degree(N1[[i]][[j]])
    degree_n2[[i]][[j]]<- degree(N3[[i]][[j]])
    degree_n3[[i]][[j]]<- degree(N3[[i]][[j]])
    degree_w[[i]][[j]] <- degree(W[[i]][[j]])
  
    z_score_n1[[i]][[j]] <- within_module_deg_z_score(N1[[i]][[j]], louvain_n1[[i]][[j]]$membership)
    z_score_n2[[i]][[j]] <- within_module_deg_z_score(N2[[i]][[j]], louvain_n2[[i]][[j]]$membership)
    z_score_n3[[i]][[j]] <- within_module_deg_z_score(N3[[i]][[j]], louvain_n3[[i]][[j]]$membership)
    z_score_w[[i]][[j]] <- within_module_deg_z_score(W[[i]][[j]], louvain_w[[i]][[j]]$membership)
  
    p_factor_n1[[i]][[j]] <- part_coeff(N1[[i]][[j]], louvain_n1[[i]][[j]]$membership)
    p_factor_n2[[i]][[j]] <- part_coeff(N2[[i]][[j]], louvain_n2[[i]][[j]]$membership)
    p_factor_n3[[i]][[j]] <- part_coeff(N3[[i]][[j]], louvain_n3[[i]][[j]]$membership)
    p_factor_w[[i]][[j]] <- part_coeff(W[[i]][[j]], louvain_w[[i]][[j]]$membership)
  
  }
}

# conversion de z_score y part_coef en rol

z0<- 1
p0 <- 0.05


node_rol <- function(z,p){
  
  result<- array(dim = length(z)) 
  
  for(i in 1:length(z)){
    #print(result)
    #print(i)
    if (is.na(z[i]) | is.na(p[i])){}
    else if (z[i]>z0 & p[i]>p0 ){result[i] <- 'Hub'}
    else if (z[i]>z0 & p[i]<p0 ) {result[i] <- 'Provincial Hub'}
    else if (z[i]<z0 & p[i]<p0 ) {result[i] <- 'Provincial Node'}
    else if (z[i]<z0 & p[i]>p0 ) {result[i] <- 'Connector Node'}
    
  }
  return(result)
}

role_count_n1 <- list()
role_count_n2 <- list()
role_count_n3 <- list()
role_count_w <- list()


for(i in 1:length(seq)){ # itero sobre los deltas
  role_count_n1_temp_H <- role_count_n2_temp_H <- role_count_n3_temp_H <- role_count_w_temp_H <- NA
  role_count_n1_temp_PH <- role_count_n2_temp_PH <- role_count_n3_temp_PH <- role_count_w_temp_PH <- NA
  role_count_n1_temp_PN <- role_count_n2_temp_PN <- role_count_n3_temp_PN <- role_count_w_temp_PN <- NA
  role_count_n1_temp_CN <- role_count_n2_temp_CN <- role_count_n3_temp_CN <- role_count_w_temp_CN <- NA

  for (j in 1:suj.count) {#itero sobre los sujetos
     
    temp <- node_rol(z_score_n1[[i]][[j]] , p_factor_n1[[i]][[j]])
    role_count_n1_temp_H[[j]] <- length(which(temp == "Hub"))
    role_count_n1_temp_PH[[j]] <- length(which(temp == "Provincial Hub"))
    role_count_n1_temp_PN[[j]] <- length(which(temp == "Provincial Node"))
    role_count_n1_temp_CN[[j]] <- length(which(temp == "Connector Node"))
    
    temp <- node_rol(z_score_n2[[i]][[j]] , p_factor_n2[[i]][[j]])
    role_count_n2_temp_H[[j]] <- length(which(temp == "Hub"))
    role_count_n2_temp_PH[[j]] <- length(which(temp == "Provincial Hub"))
    role_count_n2_temp_PN[[j]] <- length(which(temp == "Provincial Node"))
    role_count_n2_temp_CN[[j]] <- length(which(temp == "Connector Node"))
    
    temp <- node_rol(z_score_n3[[i]][[j]] , p_factor_n3[[i]][[j]])
    role_count_n3_temp_H[[j]] <- length(which(temp == "Hub"))
    role_count_n3_temp_PH[[j]] <- length(which(temp == "Provincial Hub"))
    role_count_n3_temp_PN[[j]] <- length(which(temp == "Provincial Node"))
    role_count_n3_temp_CN[[j]] <- length(which(temp == "Connector Node"))
    
    temp <- node_rol(z_score_w[[i]][[j]] , p_factor_w[[i]][[j]])
    role_count_w_temp_H[[j]] <- length(which(temp == "Hub"))
    role_count_w_temp_PH[[j]] <- length(which(temp == "Provincial Hub"))
    role_count_w_temp_PN[[j]] <- length(which(temp == "Provincial Node"))
    role_count_w_temp_CN[[j]] <- length(which(temp == "Connector Node"))
    
     
   }

  role_count_n1[[i]] <- data.frame(H =role_count_n1_temp_H, 
                                   PH =role_count_n1_temp_PH,
                                   PN =role_count_n1_temp_PN, 
                                   CN =role_count_n1_temp_CN)
  role_count_n2[[i]] <- data.frame(H =role_count_n2_temp_H, 
                                   PH =role_count_n2_temp_PH,
                                   PN =role_count_n2_temp_PN, 
                                   CN =role_count_n2_temp_CN)
  role_count_n3[[i]] <- data.frame(H =role_count_n3_temp_H, 
                                   PH =role_count_n3_temp_PH,
                                   PN =role_count_n3_temp_PN, 
                                   CN =role_count_n3_temp_CN)
  role_count_w[[i]] <- data.frame(H =role_count_w_temp_H, 
                                  PH =role_count_w_temp_PH,
                                  PN =role_count_w_temp_PN, 
                                  CN =role_count_w_temp_CN)
  
  
}

#### Creacion de Dataframe con comparaciones


rol_vector<- c("H","PH", "PN", "CN")
comparison_results_df <- NULL







for(i in 1:length(seq)){ # itero sobre los deltas
 
  
  
  
  for(r in 1:length(rol_vector)){# itero sobre los roles H, PH, PN, CN
  rol <- rol_vector[r]
  # t.test
  #p.valor.n1 <- t.test(role_count_n1[[i]][,rol],role_count_w[[i]][,rol])$p.value
  #p.valor.n2 <- t.test(role_count_n2[[i]][,rol],role_count_w[[i]][,rol])$p.value
  #p.valor.n3 <- t.test(role_count_n3[[i]][,rol],role_count_w[[i]][,rol])$p.value
  #p.valor.w <- NA
  
  
  #ANOVA
  p.valor.n1 <- unlist(summary(aov(role_count_n1[[i]][,rol]~role_count_w[[i]][,rol]))[[1]])[[9]]
  p.valor.n2 <- unlist(summary(aov(role_count_n2[[i]][,rol]~role_count_w[[i]][,rol]))[[1]])[[9]]
  p.valor.n3 <- unlist(summary(aov(role_count_n3[[i]][,rol]~role_count_w[[i]][,rol]))[[1]])[[9]]
  p.valor.w <- NA
  
  #Wilcoxon
  #p.valor.n1 <- wilcox.test(role_count_n1[[i]][,rol],role_count_w[[i]][,rol], conf.level = 0.95)$p.value
  #p.valor.n2 <- wilcox.test(role_count_n2[[i]][,rol],role_count_w[[i]][,rol], conf.level = 0.95)$p.value
  #p.valor.n3 <- wilcox.test(role_count_n3[[i]][,rol],role_count_w[[i]][,rol], conf.level = 0.95)$p.value
  #p.valor.w <- NA
  
  
  N1_mean_temp <- mean(role_count_n1[[i]][,rol])
  N2_mean_temp <- mean(role_count_n2[[i]][,rol])
  N3_mean_temp <- mean(role_count_n3[[i]][,rol])
  W_mean_temp <- mean(role_count_w[[i]][,rol])
  N1_se_temp <- sd(role_count_n1[[i]][,rol])/sqrt(length(role_count_n1[[i]][,rol]))
  N2_se_temp <- sd(role_count_n2[[i]][,rol])/sqrt(length(role_count_n2[[i]][,rol]))
  N3_se_temp <- sd(role_count_n3[[i]][,rol])/sqrt(length(role_count_n3[[i]][,rol]))
  W_se_temp <- sd(role_count_w[[i]][,rol])/sqrt(length(role_count_w[[i]][,rol]))
  
  
  
  
  
  comparison_results_df <- rbind(comparison_results_df,
                             data.frame(estadio = "N1",
                                        rol = rol, 
                                        delta=seq[i],
                                        mean=N1_mean_temp, 
                                        se= N1_se_temp,
                                        p_valor = p.valor.n1) )
  
  comparison_results_df <- rbind(comparison_results_df,
                             data.frame(estadio = "N2",
                                        rol = rol, 
                                        delta=seq[i],
                                        mean=N2_mean_temp, 
                                        se= N2_se_temp,
                                        p_valor = p.valor.n2) )
  comparison_results_df <- rbind(comparison_results_df,
                             data.frame(estadio = "N3",
                                        rol = rol, 
                                        delta=seq[i],
                                        mean=N3_mean_temp, 
                                        se= N3_se_temp,
                                        p_valor = p.valor.n3) )
  comparison_results_df <- rbind(comparison_results_df,
                             data.frame(estadio = "W",
                                        rol = rol, 
                                        delta=seq[i],
                                        mean=W_mean_temp, 
                                        se= W_se_temp,
                                        p_valor = p.valor.w) )
  
  

}
}

#### graficos
pv_ref <- 0.05 #significancia



i<- 1
j<- 2

estadio_ref_vector<- c("N1","N2","N3")
role_ref_vector <- c("H","PH","PN", "CN")
role_name_ref_vector <- c("Hubs","Provincial Hubs","Povincial Nodes", "Connection Nodes")


for(i in 1:length(role_ref_vector)){
  for(j in 1:length(estadio_ref_vector)){



estadio_ref<- estadio_ref_vector[j]
role_ref <- role_ref_vector[i]
role_name_ref <-role_name_ref_vector[i]
y_label <- paste(role_name_ref)
title_label <- paste("WAKE vs.",estadio_ref)

#title_label <- paste(role_ref,"mean count")
#y_label <- role_ref



test_df <- comparison_results_df[comparison_results_df$estadio %in% c(estadio_ref,"W") &
comparison_results_df$rol %in% c(role_ref),]
test_df<- transform(test_df, p_valor = ifelse(p_valor < pv_ref, 0, NA))

pd <- position_dodge(0.001) # move them .05 to the left and right


if (all(is.na(test_df$p_valor))) {
   p<- ggplot(test_df, aes(x=delta, y=mean, colour=estadio)) + 
        geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.005,position=pd) +
        geom_line(position=pd) +
        geom_point(position=pd) +
        #geom_point(aes(y=p_valor),position=pd) +
        labs(title=title_label, x = "delta", y=y_label)
} else {
   p<-ggplot(test_df, aes(x=delta, y=mean, colour=estadio)) + 
        geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.005,position=pd) +
        geom_line(position=pd) +
        geom_point(position=pd) +
        geom_point(aes(y=p_valor),position=pd) +
        labs(title=title_label, x = "delta", y=y_label)
}
print(p)
}
}

#grid plot

library("cowplot")

plot_grid(plots[[1]],
          plots[[2]],
          plots[[3]],
          plots[[4]],
          plots[[5]],
          plots[[6]],
          plots[[7]],
          plots[[8]],
          plots[[9]],
          plots[[10]],
          plots[[11]],
          plots[[12]],ncol = 3)

Cannot convert object of class logical into a grob.Removed 49 rows containing missing values (geom_point).Removed 51 rows containing missing values (geom_point).Removed 47 rows containing missing values (geom_point).Removed 40 rows containing missing values (geom_point).Removed 50 rows containing missing values (geom_point).Removed 51 rows containing missing values (geom_point).

### Graficos interesantes

#forma copara de dibujar los clusters
#clp <- cluster_louvain(graph = N2_no_pesados[[1]], weights = NULL)
i<-1 # delta
j<-1 # sujeto

    g<- N1[[i]][[j]] 
    roles <- as.factor(node_rol(z_score_n1[[i]][[j]] , p_factor_n1[[i]][[j]]))
    colrs <- c("olivedrab3", "tomato","darkcyan", "gold")
    V(g)$color <- colrs[roles]
    clp <- cluster_louvain(graph = g, weights = NULL)
    
    l <- layout_with_fr(g)
    
    plot(g, layout=l,vertex.label = "" , vertex.size = 7 )
legend(x=1.2, y=0.5, levels(as.factor(roles)), pch=21,col="#777777", pt.cex=2, cex=.8, bty="n", ncol=1, pt.bg=colrs)

title("Roles en N1 - sujeto 1 -  delta = 0.025", cex.main=.8)

NA
NA
NA

#Cambio de roles de los nodos

i<-1 # delta
j<-1 # sujeto
delta <- 0.065
i <- which(seq == delta)

    g<- N1[[i]][[j]] 
    roles_n1 <- node_rol(z_score_n1[[i]][[j]] , p_factor_n1[[i]][[j]])
    roles_n2 <- node_rol(z_score_n2[[i]][[j]] , p_factor_n2[[i]][[j]])
    roles_n3 <-node_rol(z_score_n3[[i]][[j]] , p_factor_n3[[i]][[j]])
    roles_w <- node_rol(z_score_w[[i]][[j]] , p_factor_w[[i]][[j]])
    
    roles <- rbind(roles_n1, roles_n2)
    roles <- rbind(roles, roles_n3)
    roles <- rbind(roles, roles_w)
    roles <- as.data.frame(roles)
    rownames(roles) <- c("N1","N2","N3","W")
    regions <- read.table( paste0(getwd(),"/datasets/aal_no_rep.csv"), header=F, sep = ",")
    colnames(roles) <- regions[,1]
    #roles[is.na(roles)] <- ""
    roles <-as.data.frame(t(roles))
    str(roles)

'data.frame':   116 obs. of  4 variables:
 $ N1: Factor w/ 4 levels "Connector Node",..: 1 2 2 2 4 4 2 2 4 4 ...
  ..- attr(*, "names")= chr  "Precentral_L" "Precentral_R" "Frontal_Sup_L" "Frontal_Sup_R" ...
 $ N2: Factor w/ 4 levels "Connector Node",..: 1 1 2 2 4 4 2 2 4 3 ...
  ..- attr(*, "names")= chr  "Precentral_L" "Precentral_R" "Frontal_Sup_L" "Frontal_Sup_R" ...
 $ N3: Factor w/ 4 levels "Connector Node",..: 1 2 2 2 4 1 2 1 4 1 ...
  ..- attr(*, "names")= chr  "Precentral_L" "Precentral_R" "Frontal_Sup_L" "Frontal_Sup_R" ...
 $ W : Factor w/ 4 levels "Connector Node",..: 1 4 3 2 4 1 2 1 4 NA ...
  ..- attr(*, "names")= chr  "Precentral_L" "Precentral_R" "Frontal_Sup_L" "Frontal_Sup_R" ...

    plot(roles)

    View(roles)
    roles_sin_nas <- roles[(!is.na(roles$N1) & !is.na(roles$N2) & !is.na(roles$N3) & !is.na(roles$W)),]
    roles_sin_repetidos <- roles[(roles$N1 != roles$N2) 
                                 | (roles$N1 != roles$N3) 
                                 | (roles$N1 != roles$W) 
                                 | (roles$N2 != roles$N3) 
                                 | (roles$N2 != roles$W) 
                                 | (roles$N3 != roles$W),]
    roles_sin_repetidos <- na.omit(roles_sin_repetidos)
    roles_sin_repetidos$region<-rownames(roles_sin_repetidos)
    View(roles_sin_repetidos)

    roles_sin_repetidos <- as.data.frame(roles_sin_repetidos)

    # separamos en 3 porque no se ven bien todos en 1 gráfico
    roles_sin_repetidos_1 <- roles_sin_repetidos[1:26,]
    roles_sin_repetidos_2 <- roles_sin_repetidos[27:52,]
    roles_sin_repetidos_3 <- roles_sin_repetidos[53:78,]
    
    roles_sin_repetidos_1 %>% gather(delta,value, N1,N2,N3,W) %>%
      ggplot() +
      geom_point(aes(x=delta, y=region, colour=value), size=3) +
      ylab(label="Región del cerebro") +
      xlab(label = "Estadío del sueño")

    roles_sin_repetidos_2 %>% gather(delta,value, N1,N2,N3,W) %>%
      ggplot() +
      geom_point(aes(x=delta, y=region, colour=value), size=3) +
      ylab(label="Región del cerebro") +
      xlab(label = "Estadío del sueño")

    roles_sin_repetidos_3 %>% gather(delta,value, N1,N2,N3,W) %>%
      ggplot() +
      geom_point(aes(x=delta, y=region, colour=value), size=3) +
      ylab(label="Región del cerebro") +
      xlab(label = "Estadío del sueño")

NA
NA

IGNORAR CODIGO A PARTIR DE ACA

#### graficos
pv_ref <- 0.05 #significancia

plots <- NA

i<- 1
j<- 2

estadio_ref_vector<- c("N1","N2","N3")
role_ref_vector <- c("H","PH","PN", "CN")
role_name_ref_vector <- c("Hubs","Provincial Hubs","Povincial Nodes", "Connection Nodes")


for(i in 1:length(role_ref_vector)){
  for(j in 1:length(estadio_ref_vector)){



estadio_ref<- estadio_ref_vector[j]
role_ref <- role_ref_vector[i]
role_name_ref <-role_name_ref_vector[i]
y_label <- paste(role_name_ref)
title_label <- paste("WAKE vs.",estadio_ref)

#title_label <- paste(role_ref,"mean count")
#y_label <- role_ref



test_df <- comparison_results_df[comparison_results_df$estadio %in% c(estadio_ref,"W") &
comparison_results_df$rol %in% c(role_ref),]
test_df<- transform(test_df, p_valor = ifelse(p_valor < pv_ref, 0, NA))

pd <- position_dodge(0.001) # move them .05 to the left and right


if (all(is.na(test_df$p_valor))) {
   p<- ggplot(test_df, aes(x=delta, y=mean, colour=estadio)) + 
        geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.005,position=pd) +
        geom_line(position=pd) +
        geom_point(position=pd) +
        #geom_point(aes(y=p_valor),position=pd) +
        labs(title=title_label, x = "delta", y=y_label)
} else {
   p<-ggplot(test_df, aes(x=delta, y=mean, colour=estadio)) + 
        geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.005,position=pd) +
        geom_line(position=pd) +
        geom_point(position=pd) +
        geom_point(aes(y=p_valor),position=pd) +
        labs(title=title_label, x = "delta", y=y_label)
}
plots <-cbind(plots,list(p))
}
}

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuXG5cblxuXG5wbG90X3JvbGVzKHNlcSwgTjFfbWVhbl9zZXJpZSAsIE4xX3N0ZF9zZXJpZSwgV19tZWFuX3NlcmllICwgV19zdGRfc2VyaWUsIGVzdGFkaW89IFwiTjFcIiwgcm9sID1yb2xfbmFtZSlcbnBsb3Rfcm9sZXMoc2VxLCBOMl9tZWFuX3NlcmllICwgTjJfc3RkX3NlcmllLCBXX21lYW5fc2VyaWUgLCBXX3N0ZF9zZXJpZSwgZXN0YWRpbz0gXCJOMlwiLCByb2wgPXJvbF9uYW1lKVxucGxvdF9yb2xlcyhzZXEsIE4zX21lYW5fc2VyaWUgLCBOM19zdGRfc2VyaWUsIFdfbWVhbl9zZXJpZSAsIFdfc3RkX3NlcmllLCBlc3RhZGlvPSBcIk4zXCIsIHJvbCA9cm9sX25hbWUpXG5gYGAifQ== -->

```r




plot_roles(seq, N1_mean_serie , N1_std_serie, W_mean_serie , W_std_serie, estadio= "N1", rol =rol_name)
plot_roles(seq, N2_mean_serie , N2_std_serie, W_mean_serie , W_std_serie, estadio= "N2", rol =rol_name)
plot_roles(seq, N3_mean_serie , N3_std_serie, W_mean_serie , W_std_serie, estadio= "N3", rol =rol_name)

for(i in 1:length(seq)){ # itero sobre los deltas
  df <- role_count_n1[[i]]
  plot(df)
}

for(i in 1:length(seq)){ # itero sobre los deltas
  df <- role_count_n2[[i]]
  plot(df)
}

for(i in 1:length(seq)){ # itero sobre los deltas
  df <- role_count_n3[[i]]
  print(paste0("Delta ",seq[i]," - Mean H ",mean(df$H)))
  print(paste0("Delta ",seq[i]," - Mean PH ",mean(df$PH)))
  print(paste0("Delta ",seq[i]," - Mean PN ",mean(df$PN)))
  print(paste0("Delta ",seq[i]," - Mean CN ",mean(df$CN)))
  plot(df)
}

for(i in 1:length(seq)){ # itero sobre los deltas
  df <- role_count_w[[i]]
  plot(df)
}

role_count_n1_temp_H <- NA
  role_count_n2_temp_H <- NA
  role_count_n3_temp_H <- NA
  role_count_w_temp_H <- NA
  
  role_count_n1_temp_PH <- NA
  role_count_n2_temp_PH <- NA
  role_count_n3_temp_PH <- NA
  role_count_w_temp_PH <- NA
  
  role_count_n1_temp_PN <- NA
  role_count_n2_temp_PN <- NA
  role_count_n3_temp_PN <- NA
  role_count_w_temp_PN <- NA
  
  role_count_n1_temp_CN <- NA
  role_count_n2_temp_CN <- NA
  role_count_n3_temp_CN <- NA
  role_count_w_temp_CN <- NA
  
  
  
role_count_n1 <- list()
role_count_n2 <- list()
role_count_n3 <- list()
role_count_w <- list()

louvain_n1<-list()
louvain_n2<-list()
louvain_n3<-list()
louvain_w<-list()

degree_n1<-list()
degree_n2<-list()
degree_n3<-list()
degree_w<-list()

z_score_n1<-list()
z_score_n2<-list()
z_score_n3<-list()
z_score_w<-list()

p_factor_n1<-list()
p_factor_n2<-list()
p_factor_n3<-list()
p_factor_w<-list()