Network analysis workshop
Adrián Alacreu Crespo
2022-07-14
1. ¿Qué es el network analysis?
El network analysis consiste en la visualización gráfica de una gran cantidad de asociaciones entre variables. En el network analysis tendremos por un lado los nodos, cada variable es un nodo, y las edges o líneas, que son las conexiones (asociaciones) entre dichas variables. Es por tanto, un método que permite, por un lado, visualizar estas asociaciones entre variables, pero también permite ver que variables son las más centrales en la red, si hay variables “puente”, que variables se pueden agrupar como un Análisis factorial exploratorio, e incluso podemos ver que variables influyen sobre otras cuando obtenemos muestras longitudinales con Ecological momentary assesment
1.1. Ventajas del network analysis
Control de variables: Todas las variables que se introducen en el modelo están controladas por el resto de variables del modelo
Complejidad real de los fenomenos psicológicos: Es un método que puede ayudar a comprender la complejidad real de los fenomenos psicológicos, en muchos estudios se observa una o pocas variables sin tener en cuenta el rico mundo de variables que pueden predecir un mismo fenómeno
Diferencias entre sub-grupos de pacientes: Puede ayudar a comprender diferencias en sub-grupos de pacientes, si sus networks son efectivamente diferentes
Tratamiento con datos individuales: Para un mismo paciente, si se recogen una gran cantidad de datos durante un largo periodo de tiempo podría ayudar a diferenciar que síntoma es el más importante o el que lleva a empeorar el resto de síntomas
1.2. Inconvenientes del network analysis
Cantidad de muestra: Para 10 síntomas, se necesitan al menos 55 observaciones (10 nodos + 10 × 9/2 posibles interacciones), para 20 se necesitan 210 y para 50 se necesitan 1250
Exploración: En muchos casos es díficil hacer hipótesis sobre como va a aparecer una red compleja de síntomas o variables unidos unos con otros, por ello, muchas veces se usa como método exploratorio y no confirmatorio
Variables sin sentido: Ojo con meter variables por meter, esto es síntoma también de que es un método que se usa más bien de forma exploratoria, hay que tener coherencia con que variables se incluyen en el modelo teniendo en cuenta que es un analisis de asociaciones entre variables, por ejemplo, si introducimos insomnio e hipersomnia, estás variables van a estar muy relacionadas (de forma negativa entre si), otro caso sería usar el score total de la depresión y los items de la BDI…
2. Preparación
2.1 Instalación de librerias
Para empezar vamos a instalar una serie de librerias que utilizaremos a posteriori en el network analysis y cargarlas en R
install.packages("mgm")
install.packages("qgraph")
install.packages("bootnet")
install.packages("dplyr")
install.packages("ggplot2")
install.packages("readxl")
install.packages("NetworkComparisonTest")library(mgm) # libreria de network analysis para incluir variables cualitativas
library(qgraph) # libreria de network analysis
library(bootnet) # libreria para hacer analysis de estabilidad de las network
library(dplyr)
library(ggplot2)
library(readxl)
library(NetworkComparisonTest)2.2 El estudio (los datos que vamos a usar)
Vamos a utilizar los datos del estudio Alacreu-Crespo et al., 2020
El estudio pretendía identificar si el dolor psicológico máximo en los ultimos 15 días (evaluado con escala analogo visual), era una variable predictora de un evento suicida a 1 año desde la hospitalización tanto en pacientes con intento de suicidio previo como en pacientes hospitalizados en emergencias psiquiatricas por otra causa relacionada con un trastorno del estado de ánimo,
Se hizo una evaluación al ser hospitalizados, incluyendo varias escalas de depresión, sociodemográficos, medicación y la conducta suicida previa,
A continuacion se realizó un seguimiento que duró 1 año donde se evaluó los eventos suicidas, incluyendo:
Hospitalización por Ideación suicida
Hospitalización por Intento de suicidio
Muerte por suicidio
Uno de los análisis que se utilizaron fue el network analysis, para observar, entre los items de depresión de la BDI y el dolor psicológico (en la inclusión), cuales se encontraban más relacionadas con el evento suicida futuro
2.3 Cargar los datos
Vamos a cargar los siguientes datos del fichero “Psych_Posp.xlsx”, con el comando importar, para ello teneis que tener los datos en la misma carpeta que el proyecto R, recordad que el cambio de una sola coma puede significar que el programa no funcione:
BBDD <- read_excel("Psych_Posp.xlsx")Ahora vamos a seleccionar los datos que necesitamos para hacer el primer network analysis, seleccionaremos la variables conducta suicida, llamada SB, la variable dolor psicológico y todos los ítems de la depresión:
BBDDNet1 <- BBDD[, c("Psych_Pain", "BDI_1", "BDI_2", "BDI_3", "BDI_4", "BDI_5", "BDI_6",
"BDI_7", "BDI_8", "BDI_9", "BDI_10", "BDI_11", "BDI_12","BDI_13", "SB_0_12")]Antes de empezar a calcular el network analysis, vamos a renombrar las variables y simplificarlas de cara a la figura, como mucho 3 carácteres, si ponemos más no aparecerá bien en la figura, hay que renombrarlas por orden:
names(BBDDNet1) <-c("Psy", "Sad", "Pes", "Fai", "Ple", "Cul","Hat", "Sui", "Rel", "Ind", "Bod", "Ene","Fat", "App", "SE")Además, les vamos a asignar un grupo a cada variable, esto lo que hará que en el gráfico aparezca el nodo asignado a una variable de cada grupo con un color diferente, por ejemplo, tenemos la BDI por un lado (grupo 1), el dolor psicológico que queremos destacarlo (grupo 2) y los eventos suicidas que sería nuestra variable de interés (grupo 3):
groups <- structure(list(BDI = c(2,3,4,5,6,7,8,9,10,11,12,13, 14),
EVA = c(1),
SE = c(15)),
Names = c("BDI",
"EVA", "SE"))En list indicaremos las variables (por órden) asignadas a cada grupo, y en names el nombre de cada grupo.
Para que no haya problemas a la hora de calcular la network, es recomendable eliminar los NAs
BBDDNet1 <- na.omit(BBDDNet1) Por último, para entender bien el gráfico que va a salir, será necesario especificar los nombres al completo de cada una de las variables, para ello vamos a cargar un archivo de texto .txt con los nombres de las variables ordenadas de la misma forma que tenemos nuestra base de datos BBDDnet1:
Recordad que el archivo debe estar en la misma carpeta que vuestro proyecto R
NamesNet <- scan("NamesNet.txt",what = "character", sep = "\n")3. Network analysis
Ahora estamos preparados para hacer el network analysis, vamos a empezar realizando un network analysis básico con correlaciones policóricas, ya que nuestras variables son ordinales (los items de BDI y la escala analogo visual), para eso vamos a empezar calculando las correlaciones:
cormatrix <- cor_auto(BBDDNet1)## Variables detected as ordinal: Sad; Pes; Fai; Ple; Cul; Hat; Sui; Rel; Ind; Bod; Ene; Fat; App; SEcormatrix## Psy Sad Pes Fai Ple Cul Hat
## Psy 1.0000000 0.48075333 0.35029296 0.3515318 0.3082138 0.3894204 0.3858792
## Sad 0.4807533 1.00000000 0.66530434 0.4918420 0.5353452 0.4794683 0.5283717
## Pes 0.3502930 0.66530434 1.00000000 0.5059793 0.5290393 0.5932658 0.5971790
## Fai 0.3515318 0.49184205 0.50597934 1.0000000 0.5301659 0.5399505 0.5770538
## Ple 0.3082138 0.53534515 0.52903932 0.5301659 1.0000000 0.4774027 0.5185411
## Cul 0.3894204 0.47946826 0.59326579 0.5399505 0.4774027 1.0000000 0.7011537
## Hat 0.3858792 0.52837173 0.59717900 0.5770538 0.5185411 0.7011537 1.0000000
## Sui 0.4834388 0.55973063 0.56148922 0.4436743 0.4317125 0.4667335 0.5212247
## Rel 0.1663819 0.38713090 0.47050824 0.3080353 0.5332418 0.1867429 0.2644579
## Ind 0.3058797 0.51320783 0.58945757 0.3863370 0.4926662 0.4900396 0.4669761
## Bod 0.2332223 0.34449742 0.44028901 0.4229285 0.2901075 0.4621750 0.5816694
## Ene 0.2975452 0.45874901 0.52159863 0.3137974 0.4307041 0.4069514 0.4368635
## Fat 0.3202780 0.47273948 0.44555845 0.3157568 0.4588779 0.3449310 0.4204393
## App 0.3715288 0.33425228 0.31507468 0.2194203 0.3154568 0.3011412 0.3449086
## SE 0.3593113 0.08952825 0.08001154 0.3139916 0.1131460 0.1664143 0.1247404
## Sui Rel Ind Bod Ene Fat App
## Psy 0.4834388 0.1663819 0.30587975 0.2332223 0.297545194 0.32027795 0.37152877
## Sad 0.5597306 0.3871309 0.51320783 0.3444974 0.458749012 0.47273948 0.33425228
## Pes 0.5614892 0.4705082 0.58945757 0.4402890 0.521598628 0.44555845 0.31507468
## Fai 0.4436743 0.3080353 0.38633699 0.4229285 0.313797369 0.31575678 0.21942028
## Ple 0.4317125 0.5332418 0.49266621 0.2901075 0.430704101 0.45887792 0.31545680
## Cul 0.4667335 0.1867429 0.49003962 0.4621750 0.406951408 0.34493097 0.30114123
## Hat 0.5212247 0.2644579 0.46697613 0.5816694 0.436863535 0.42043935 0.34490856
## Sui 1.0000000 0.3128572 0.36103527 0.3936327 0.337081702 0.36909032 0.30817407
## Rel 0.3128572 1.0000000 0.34081202 0.1710022 0.366901016 0.35168901 0.27808704
## Ind 0.3610353 0.3408120 1.00000000 0.3808378 0.611956589 0.57541265 0.23330623
## Bod 0.3936327 0.1710022 0.38083779 1.0000000 0.282259832 0.29614342 0.25859107
## Ene 0.3370817 0.3669010 0.61195659 0.2822598 1.000000000 0.71408126 0.31061807
## Fat 0.3690903 0.3516890 0.57541265 0.2961434 0.714081256 1.00000000 0.33056465
## App 0.3081741 0.2780870 0.23330623 0.2585911 0.310618066 0.33056465 1.00000000
## SE 0.2557605 0.1433200 0.02199074 0.1490440 -0.007013743 0.03331642 0.08145881
## SE
## Psy 0.359311325
## Sad 0.089528246
## Pes 0.080011542
## Fai 0.313991630
## Ple 0.113145969
## Cul 0.166414302
## Hat 0.124740385
## Sui 0.255760506
## Rel 0.143320031
## Ind 0.021990740
## Bod 0.149044029
## Ene -0.007013743
## Fat 0.033316423
## App 0.081458812
## SE 1.000000000Ahora ya tendriamos las correlaciones, y vamos a calcular la network, tenemos que tener en cuenta que estamos ante variables que son síntomas psiquiatricos en pacientes con depresión mayor, y esto probablemente de lugar a muchas relaciones entre variables (muchas conexiones) pero muchas de ellas serán relaciones falsas (es decir relaciones que violen el error de tipo 1, falsos positivos), y tenemos que intentar tratar con ello, para ello el paquete qgraph aplica un umbral (el cual podemos decidir activarlo o no) que dicho de forma burda es como un ajuste post-hoc de las relaciones, donde las relaciones muy pequeñas serán eliminadas de la network:
Net1 <- qgraph(cor_auto(BBDDNet1), graph = "glasso", sampleSize = nrow(BBDDNet1), layout = "spring", threshold = T, tuning = 0.5, cut = 0, palette = "colorblind", theme = "Hollywood", esize = 10, groups = groups, vsize = 7.5, nodeNames = NamesNet, legend = T, legend.cex = 0.375, layoutScale = c(1,1), layoutOffset= c(0,0), GLratio = 1.25, title = "Network spring")## Variables detected as ordinal: Sad; Pes; Fai; Ple; Cul; Hat; Sui; Rel; Ind; Bod; Ene; Fat; App; SE## Note: Network with lowest lambda selected as best network: assumption of sparsity might be violated.
Los comandos que hemos utilizado son los siguientes:
cor_auto: indica la matriz de correlaciones que se va a usar
qgraph: indica si la network se calcula como una matriz de correlaciones normal “cor”, como una matriz de correlaciones parciales “pcor” o como una matriz donde se busca una estimación dispersa óptima de la matriz de correlaciones parciales “glasso”
layout: tipo de figura “spring” o “circle”
Treshold: donde se activa el umbral o no para evitar el error de tipo 1
tuning: es un parametro que ayuda a obtener la estimación dispersa más optima, solo puede usarse con glasso, a mayor tunning más dispersión
Ahora bien, interpretar esta network en función de donde están colocados los síntomas (arriba, abajo, izquierda o derecha) es irrelevante, aquí solo podemos ver que variables están conectadas con otras de forma significativa teniendo en cuenta todas las variables de la network en el analisis, en caso de que haya alguna variable de interes especial, como en nuestro caso el suicidal event, podemos utilizar el comando flow, donde si tendrá importancia la situación de cada variable:
Net2 <- flow(Net1, "SE", palette = "colorblind", theme = "Hollywood", esize = 10, groups = groups, vsize = 7.5, nodeNames = NamesNet, legend = T, legend.cex = 0.375, layoutScale = c(1,1), layoutOffset= c(0,0), GLratio = 1.25, title = "Network flow")
En este caso el orden si tiene importancia, porque podemos ver que variable se relacionan en un primer orden con nuestra variable dependiente, cuales son de 2º orden ya que predicen las variables que predicen a la dependiente y así sucesivamente
3.1 Analisis de centralidad
Otro método para ver la importancia de un nodo en el network son los analisis de centralidad, hay distintos índices de centralidad, aquí indico los más comunes:
Strength: La suma de todas las relaciones positivas (todas las edges positivas)
Betweenness: Calcula que nodo muestra el camino más corto si se empieza desde cualquier otro nodo
Closeness: La suma de todas las relaciones directas e indirectas con un mismo nodo
Expected influence: La suma en valor absoluto de todas las relaciones tanto positivas como negativas directas
El problema es que muchos de ellos son muy criticados o se interpretan de forma erronea, aquí os dejo el cálculo:
centralityPlot <- centralityPlot(Net1, scale = "z-scores", include = c("Strength", "Betweenness", "Closeness", "ExpectedInfluence"), decreasing = F, theme_bw = T)## Note: z-scores are shown on x-axis rather than raw centrality indices.
centralityTable(Net2)## graph type node measure value
## 1 graph 1 NA Psy Betweenness 1.077692009
## 2 graph 1 NA Sad Betweenness 0.584845054
## 3 graph 1 NA Pes Betweenness 2.556232876
## 4 graph 1 NA Fai Betweenness 0.091998098
## 5 graph 1 NA Ple Betweenness -0.597987639
## 6 graph 1 NA Cul Betweenness 0.584845054
## 7 graph 1 NA Hat Betweenness 0.289136881
## 8 graph 1 NA Sui Betweenness -0.992265204
## 9 graph 1 NA Rel Betweenness -0.105140684
## 10 graph 1 NA Ind Betweenness 0.486275663
## 11 graph 1 NA Bod Betweenness -0.992265204
## 12 graph 1 NA Ene Betweenness -0.992265204
## 13 graph 1 NA Fat Betweenness -0.992265204
## 14 graph 1 NA App Betweenness -0.992265204
## 15 graph 1 NA SE Betweenness -0.006571293
## 16 graph 1 NA Psy Closeness 0.233423422
## 17 graph 1 NA Sad Closeness 1.003389319
## 18 graph 1 NA Pes Closeness 2.031522931
## 19 graph 1 NA Fai Closeness -0.002317995
## 20 graph 1 NA Ple Closeness 0.094667318
## 21 graph 1 NA Cul Closeness 0.879618740
## 22 graph 1 NA Hat Closeness 0.310742843
## 23 graph 1 NA Sui Closeness 0.117871245
## 24 graph 1 NA Rel Closeness 0.796107534
## 25 graph 1 NA Ind Closeness -0.077888689
## 26 graph 1 NA Bod Closeness -0.955500331
## 27 graph 1 NA Ene Closeness -1.315897572
## 28 graph 1 NA Fat Closeness -1.624307188
## 29 graph 1 NA App Closeness -1.410361463
## 30 graph 1 NA SE Closeness -0.081070113
## 31 graph 1 NA Psy Strength 0.847350020
## 32 graph 1 NA Sad Strength -0.142419721
## 33 graph 1 NA Pes Strength 1.855958032
## 34 graph 1 NA Fai Strength 0.146797165
## 35 graph 1 NA Ple Strength -0.203546311
## 36 graph 1 NA Cul Strength 1.063582769
## 37 graph 1 NA Hat Strength 0.951320606
## 38 graph 1 NA Sui Strength -0.691267663
## 39 graph 1 NA Rel Strength 0.190506886
## 40 graph 1 NA Ind Strength -0.056565227
## 41 graph 1 NA Bod Strength -1.600271284
## 42 graph 1 NA Ene Strength 0.451025610
## 43 graph 1 NA Fat Strength -0.148015568
## 44 graph 1 NA App Strength -1.989772032
## 45 graph 1 NA SE Strength -0.674683281
## 46 graph 1 NA Psy ExpectedInfluence 1.039998162
## 47 graph 1 NA Sad ExpectedInfluence 0.049600422
## 48 graph 1 NA Pes ExpectedInfluence 2.049246127
## 49 graph 1 NA Fai ExpectedInfluence 0.339000813
## 50 graph 1 NA Ple ExpectedInfluence -0.011564953
## 51 graph 1 NA Cul ExpectedInfluence -0.184460689
## 52 graph 1 NA Hat ExpectedInfluence 1.144034717
## 53 graph 1 NA Sui ExpectedInfluence -0.499595760
## 54 graph 1 NA Rel ExpectedInfluence -1.058090531
## 55 graph 1 NA Ind ExpectedInfluence 0.135509389
## 56 graph 1 NA Bod ExpectedInfluence -1.409176136
## 57 graph 1 NA Ene ExpectedInfluence 0.643422288
## 58 graph 1 NA Fat ExpectedInfluence 0.044001023
## 59 graph 1 NA App ExpectedInfluence -1.798924019
## 60 graph 1 NA SE ExpectedInfluence -0.4830008553.2 Estabilidad de la network
En este punto ya tenemos una network calculada con sus índices de centralidad, pero a continuación debemos testar que está network es valida, o mejor dicho, estable. ¿Qué quiere decir esto? Que aunque repitamos una y otra vez los mismos análisis probablemente obtendremos una network las mismas conexiones o al menos muy parecidas
Esto lo haremos con el paquete “bootnet”, un paquete creado recientemente por Eiko Fried y Satcha Epskmanp que nos permiten testar que nuestras networks son en realidad validas, importante, poned los mismos parametros que tenía vuestra network en cuanto a threshold, tunning etc.
boot1a <- bootnet(BBDDNet1, default ="glasso", threshold = T, tuning = 0.5, lambda.min.ratio = 0.01, nBoots = 500, nCores = 8)
plot(boot1a, labels = F, order = "sample") 
plot(boot1a, plot = "interval", split0 = TRUE, order="sample", labels = F)
Aquí tenemos a mano izquierda cada una de las edges, si hubiesemos indicado labels T en los plot, veriamos de que conexiones concretas se tratan, lo que nos tiene que preocupar es si hay muchas conexiones que incluyan el 0 en su intervalo de confianza, en ese caso la network no sería estable, en caso de que haya pocas, indica que precisamente esas conexiones son poco estables, si es alguna de las conexiones de interés deberiamos interpretarla con cautela
También vamos a testar si los índices de centralidad se mantienen estables, esto lo haremos con el mismo paquete
boot1b <- bootnet(BBDDNet1, default ="glasso", tuning = 0.5, threshold = T, lambda.min.ratio = 0.01, statistics = c("Strength", "Closeness", "Betweenness", "ExpectedInfluence"), nBoots = 1000, type = "case", nCores = 8)
plot(boot1b, statistics = c("Strength", "Betweenness", "Closeness", "ExpectedInfluence"))
corStability(boot1b)## === Correlation Stability Analysis ===
##
## Sampling levels tested:
## nPerson Drop% n
## 1 90 75.1 107
## 2 118 67.3 114
## 3 146 59.6 95
## 4 174 51.8 96
## 5 203 43.8 116
## 6 231 36.0 96
## 7 259 28.3 81
## 8 287 20.5 96
## 9 315 12.7 90
## 10 343 5.0 109
##
## Maximum drop proportions to retain correlation of 0.7 in at least 95% of the samples:
##
## betweenness: 0
## - For more accuracy, run bootnet(..., caseMin = 0, caseMax = 0.05)
##
## closeness: 0
## - For more accuracy, run bootnet(..., caseMin = 0, caseMax = 0.05)
##
## expectedInfluence: 0.283
## - For more accuracy, run bootnet(..., caseMin = 0.205, caseMax = 0.36)
##
## strength: 0.283
## - For more accuracy, run bootnet(..., caseMin = 0.205, caseMax = 0.36)
##
## Accuracy can also be increased by increasing both 'nBoots' and 'caseN'.En caso de que la centralidad el coeficiente CS no debe ser inferior a 0,25, y preferiblemente superior a 0,5, en caso contrario tenemos muy poca estabilidad de nuestros indicadores de centralidad
3.3 Comparación de dos networks
Otra cosa que podemos hacer con las networks con datos transversales es comparar networks de dos grupos diferentes, un buen ejemplo de este tipo de aproximación la encontramos en el artículo de Koenders et al. 2015
Para ello vamos a usar el paquete NCT (Network Comparison Test), primero dividiremos a nuestra muestra en pacientes con suicidal event y pacientes sin suicidal event, y cargaremos los datos, vamos a nombrar las variables y poner otra vez los nombres:
BBDD_SE <- read_excel("Psych_Posp_SE.xlsx")
BBDD_NSE <- read_excel("Psych_Posp_NSE.xlsx")
BBDDNet_SE <- BBDD_SE[, c("Psych_Pain", "BDI_1", "BDI_2", "BDI_3", "BDI_4", "BDI_5", "BDI_6",
"BDI_7", "BDI_8", "BDI_9", "BDI_10", "BDI_11", "BDI_12","BDI_13")]
BBDDNet_NSE <- BBDD_NSE[, c("Psych_Pain", "BDI_1", "BDI_2", "BDI_3", "BDI_4", "BDI_5", "BDI_6",
"BDI_7", "BDI_8", "BDI_9", "BDI_10", "BDI_11", "BDI_12","BDI_13")]
names(BBDDNet_SE) <-c("Psy", "Sad", "Pes", "Fai", "Ple", "Cul","Hat", "Sui", "Rel", "Ind", "Bod", "Ene","Fat", "App")
names(BBDDNet_NSE) <-c("Psy", "Sad", "Pes", "Fai", "Ple", "Cul","Hat", "Sui", "Rel", "Ind", "Bod", "Ene","Fat", "App")
BBDDNet_SE <- na.omit(BBDDNet_SE)
BBDDNet_NSE <- na.omit(BBDDNet_NSE)
NamesNet2 <- scan("NamesNet2.txt",what = "character", sep = "\n")Ahora vamos a calcular las dos networks
NetSE <- qgraph(cor_auto(BBDDNet_SE), graph = "glasso", sampleSize = nrow(BBDDNet_SE), layout = "circle", threshold = F, tuning = 0.5, cut = 0, palette = "colorblind", theme = "Hollywood", esize = 10, vsize = 7.5, nodeNames = NamesNet2, legend = T, legend.cex = 0.375, layoutScale = c(1,1), layoutOffset= c(0,0), GLratio = 1.25, title = "Network SE")## Variables detected as ordinal: Sad; Pes; Fai; Ple; Cul; Hat; Sui; Rel; Ind; Bod; Ene; Fat; App
NetNSE <- qgraph(cor_auto(BBDDNet_NSE), graph = "glasso", sampleSize = nrow(BBDDNet_NSE), layout = "circle", threshold = F, tuning = 0.5, cut = 0, palette = "colorblind", theme = "Hollywood", esize = 10, vsize = 7.5, nodeNames = NamesNet2, legend = T, legend.cex = 0.375, layoutScale = c(1,1), layoutOffset= c(0,0), GLratio = 1.25, title = "Network NSE")## Variables detected as ordinal: Sad; Pes; Fai; Ple; Cul; Hat; Sui; Rel; Ind; Bod; Ene; Fat; App
Con las nuevas network vamos a usar el paquete NCT, este paquete va a comparar si las dos networks son diferentes entres sí en varios aspectos, en primer lugar si alguna tiene más relaciones que otra de forma significativa, en segundo lugar si alguna tiene alguna relación concreta distinta de forma significativa y en tercer lugar vamos a ver si hay una centralidad distinta en alguno de los nodos de forma signifiactiva:
NetworkSE <-estimateNetwork(BBDDNet_SE,default = "EBICglasso" , threshold = F, corMethod = "cor_auto")
NetworkNSE <-estimateNetwork(BBDDNet_NSE,default = "EBICglasso" , threshold = F, corMethod = "cor_auto")
NCT(NetworkSE,NetworkNSE, gamma= 0.5, it=30, test.edges = T, edges = "all", test.centrality = T, centrality=c("strength","expectedInfluence"),nodes="all", p.adjust.methods = "fdr")##
|
| | 0%
|
|== | 3%
|
|===== | 7%
|
|======= | 10%
|
|========= | 13%
|
|============ | 17%
|
|============== | 20%
|
|================ | 23%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 33%
|
|========================== | 37%
|
|============================ | 40%
|
|============================== | 43%
|
|================================= | 47%
|
|=================================== | 50%
|
|===================================== | 53%
|
|======================================== | 57%
|
|========================================== | 60%
|
|============================================ | 63%
|
|=============================================== | 67%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|====================================================== | 77%
|
|======================================================== | 80%
|
|========================================================== | 83%
|
|============================================================= | 87%
|
|=============================================================== | 90%
|
|================================================================= | 93%
|
|==================================================================== | 97%
|
|======================================================================| 100%##
## NETWORK INVARIANCE TEST
## Test statistic M: 0.3356456
## p-value 0.2
##
## GLOBAL STRENGTH INVARIANCE TEST
## Global strength per group: 6.318459 6.127752
## Test statistic S: 0.190707
## p-value 0.8333333
##
## EDGE INVARIANCE TEST
##
## Var1 Var2 p-value
## 15 Psy Sad 1.000
## 29 Psy Pes 1.000
## 30 Sad Pes 1.000
## 43 Psy Fai 0.867
## 44 Sad Fai 1.000
## 45 Pes Fai 1.000
## 57 Psy Ple 1.000
## 58 Sad Ple 1.000
## 59 Pes Ple 1.000
## 60 Fai Ple 1.000
## 71 Psy Cul 1.000
## 72 Sad Cul 1.000
## 73 Pes Cul 1.000
## 74 Fai Cul 1.000
## 75 Ple Cul 1.000
## 85 Psy Hat 1.000
## 86 Sad Hat 1.000
## 87 Pes Hat 0.867
## 88 Fai Hat 1.000
## 89 Ple Hat 1.000
## 90 Cul Hat 1.000
## 99 Psy Sui 1.000
## 100 Sad Sui 1.000
## 101 Pes Sui 1.000
## 102 Fai Sui 1.000
## 103 Ple Sui 1.000
## 104 Cul Sui 1.000
## 105 Hat Sui 1.000
## 113 Psy Rel 1.000
## 114 Sad Rel 1.000
## 115 Pes Rel 1.000
## 116 Fai Rel 1.000
## 117 Ple Rel 1.000
## 118 Cul Rel 1.000
## 119 Hat Rel 1.000
## 120 Sui Rel 1.000
## 127 Psy Ind 1.000
## 128 Sad Ind 1.000
## 129 Pes Ind 1.000
## 130 Fai Ind 1.000
## 131 Ple Ind 1.000
## 132 Cul Ind 1.000
## 133 Hat Ind 1.000
## 134 Sui Ind 1.000
## 135 Rel Ind 1.000
## 141 Psy Bod 1.000
## 142 Sad Bod 1.000
## 143 Pes Bod 1.000
## 144 Fai Bod 1.000
## 145 Ple Bod 1.000
## 146 Cul Bod 1.000
## 147 Hat Bod 1.000
## 148 Sui Bod 1.000
## 149 Rel Bod 1.000
## 150 Ind Bod 0.758
## 155 Psy Ene 1.000
## 156 Sad Ene 1.000
## 157 Pes Ene 0.867
## 158 Fai Ene 1.000
## 159 Ple Ene 1.000
## 160 Cul Ene 0.758
## 161 Hat Ene 1.000
## 162 Sui Ene 0.827
## 163 Rel Ene 0.758
## 164 Ind Ene 0.758
## 165 Bod Ene 1.000
## 169 Psy Fat 1.000
## 170 Sad Fat 0.758
## 171 Pes Fat 1.000
## 172 Fai Fat 0.000
## 173 Ple Fat 1.000
## 174 Cul Fat 1.000
## 175 Hat Fat 1.000
## 176 Sui Fat 1.000
## 177 Rel Fat 1.000
## 178 Ind Fat 1.000
## 179 Bod Fat 1.000
## 180 Ene Fat 0.758
## 183 Psy App 0.827
## 184 Sad App 0.827
## 185 Pes App 1.000
## 186 Fai App 1.000
## 187 Ple App 1.000
## 188 Cul App 1.000
## 189 Hat App 1.000
## 190 Sui App 1.000
## 191 Rel App 0.000
## 192 Ind App 1.000
## 193 Bod App 1.000
## 194 Ene App 1.000
## 195 Fat App 1.000
##
## CENTRALITY INVARIANCE TEST
##
## strength expectedInfluence
## Psy 0.3111111 0.3111111
## Sad 1.0000000 1.0000000
## Pes 1.0000000 1.0000000
## Fai 1.0000000 1.0000000
## Ple 1.0000000 1.0000000
## Cul 1.0000000 1.0000000
## Hat 1.0000000 1.0000000
## Sui 1.0000000 1.0000000
## Rel 0.9333333 0.0000000
## Ind 1.0000000 1.0000000
## Bod 1.0000000 1.0000000
## Ene 1.0000000 1.0000000
## Fat 1.0000000 1.0000000
## App 1.0000000 1.0000000Los resultados muestran que en general las dos redes no son significantivamente diferentes, sin embargo las relaciones sad-app y rel-app si son diferentes
Calcularemos también el gráfico de centralidad para ver diferencias, aunque como hemos visto solo hay una diferencia significativa en expected influence el nodo Rel:
centralityPlot(list(SE = NetworkSE, NSE = NetworkNSE), scale = "z-scores", include = "all")
Aunque no lo vamos a hacer, deberíamos re
3.4 Inclusión de variables cualitativas
Para finalizar, hemos calculado networks con un paquete que nos permitía utilizar variables ordinales o cuantitativas, el qgraph, pero, ¿qué pasa cuando incluimos variables cualitativas?¿el SE no es ya de por sí una variable categórica con 2 categorías? Lo ideal sería aplicar otro método que tenga en cuenta la distribución de las variables, sobretodo si vamos a incluir varias variables categoricas en el modelo, en este caso vamos a tener que usar el paquete mgm, lo que permite combinar en una sola red variables con gausianas, ordinales y cualitativas
Primero como siempre cargaremos los datos, y esta vez vamos a incluir más variables, como algunas comorbilidades psiquiatricas, seguiremos el mismo protocolo, indicaremos los nombres nuevos, eliminaremos los NAs etc.
BBDDNet2 <- BBDD[, c("Psych_Pain", "BDI_1", "BDI_2", "BDI_3", "BDI_4", "BDI_5", "BDI_6",
"BDI_7", "BDI_8", "BDI_9", "BDI_10", "BDI_11", "BDI_12","BDI_13", "C_Al", "C_Dr", "C_anx", "C_TCA", "SB_0_12")]
names(BBDDNet2) <-c("Psy", "Sad", "Pes", "Fai", "Ple", "Cul","Hat", "Sui", "Rel", "Ind", "Bod", "Ene","Fat", "App", "Alc", "Dru", "Anx", "ED", "SE")
BBDDNet2 <- na.omit(BBDDNet2)
NamesNet3 <- scan("NamesNet3.txt",what = "character", sep = "\n")
groups2 <- structure(list(BDI = c(2,3,4,5,6,7,8,9,10,11,12,13,14),
EVA = c(1),
SE = c(19),
Disorder = c(15,16,17,18)),
Names = c("BDI",
"EVA", "SE", "Disorder"))Ahora estamos listos para aplicar el mgm, en orden indicaremos cada variable que tipo de escala tiene en el comando type, g = gausianas, p = poison (ordinales), c = categoricas, en levels indicaremos el nº de categorias para cada variable, si son gausianas o ordinales le pondremos 1, luego usaremos qgraph para poner que nos extraiga la network finalmente:
Network_mgm <- mgm(data = BBDDNet2,
type = c("g","p", "p", "p", "p", "p", "p", "p", "p", "p","p","p", "p", "p", "c", "c", "c", "c", "c"),
levels = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2),
k = 2,
lambdaSel = "EBIC", lambdaGam = 0.075,
ruleReg = "AND")##
|
| | 0%
|
|---- | 5%
|
|------- | 11%
|
|----------- | 16%
|
|--------------- | 21%
|
|------------------ | 26%
|
|---------------------- | 32%
|
|-------------------------- | 37%
|
|----------------------------- | 42%
|
|--------------------------------- | 47%
|
|------------------------------------- | 53%
|
|----------------------------------------- | 58%
|
|-------------------------------------------- | 63%
|
|------------------------------------------------ | 68%
|
|---------------------------------------------------- | 74%
|
|------------------------------------------------------- | 79%
|
|----------------------------------------------------------- | 84%
|
|--------------------------------------------------------------- | 89%
|
|------------------------------------------------------------------ | 95%
|
|----------------------------------------------------------------------| 100%
## Note that the sign of parameter estimates is stored separately; see ?mgmNet_mgm <- qgraph(Network_mgm$pairwise$wadj,
layout = 'spring', repulsion = 1.3,
edge.color = Network_mgm$pairwise$edgecolor, labels = colnames(BBDDNet2),
groups = groups2,palette = "colorblind", theme = "Hollywood", legend.cex = 0.375, layoutScale = c(1,1), layoutOffset= c(0,0), GLratio = 1.25, title = "Network mgm", nodeNames = NamesNet3, legend = T, sampleSize = nrow(Network_mgm))
4. Bibliografía y tutoriales recomendados
Curso de network analysis en abierto
Blog de Network analysis de Eiko Fried
Blog Network analysis y tutoriales
Epskamp, S. & Fried, Eiko I. (2018). A Tutorial on Regularized Partial Correlation Networks. Psychological Methods
Epskamp, S., Borsboom, D. & Fried, E.I. Estimating Psychological Networks and their Accuracy: A Tutorial Paper. (2017). Behavior Research Methods. doi:10.3758/s13428-017-0862-1
Jones, P. J., Mair, P., & McNally, R. (2018). Visualizing Psychological Networks: A Tutorial in R. Frontiers in Psychology, 9, 1742. https://doi.org/10.3389/fpsyg.2018.01742