网络分析教程

1 基本的网络分析原理与代码

  本部分是网络分析的基础部分，入门和基本概念等。基于的资料和代码参考有几下几个。

1. 网络分析（1）：横断面网络分析及代码实例；给出了公开的代码和数据库，便于复现和理解；参考文献

2. 网络分析（2）：桥梁分析、网络比较分析及代码实例与纵向网络分析；作为第1部分的拓展，也有相应的代码和教程复现；这个公众号关于网络分析的介绍就是这两个文章；参考文献

3. 真的有写的这么详细的网络分析教程嘛？；原理的讲解非常清楚，但是数据和代码复现要收费，作为补充材料；

4. 你见过有调节的网络分析吗？保姆级教程！！！，3的拓展，原理清楚，无法代码复现。补充一个有调节的网络分析——另一种方法；包括论文全文和数据代码复现材料：Effects of an MPPI on mental health文件夹。

---

1.1 基本复现例子一

  首先基于1的代码复现进行操作。网络分析（Network Analysis），原本较常见于社会学领域，名为“社会网络分析”（Social Network Analysis），又被称作“结构分析方法”。起初，是由社会学家根据图论等方法发展起来的定量分析方法。近些年来，社会网络分析方法经过心理学科研工作者的调整和引入，在心理学科研领域占据了越来越重要的位置。本文将从心理学科研领域视角介绍该方法，而非社会学领域的视角，因此不会过多讨论社会学视角中的“社会网络”等成分

rm(list = ls()) #痕迹清除
#安装必要的包
library(bootnet) # Bootstrap Methods for Various Network Estimation Routines
library(qgraph) # Graph Plotting Methods, Psychometric Data Visualization and Graphical Model Estimation
data('attitude') #自带的数据库
dataw22=as.data.frame(attitude[,c(1:7)]) #选择你的待分析数据
###把数据变成矩阵形式
data222.cor=cor_auto(dataw22,detectOrdinal =FALSE) #得到一个偏相关系数矩阵
###命名节点和变量
names=c("RA","CO","PR","LE","RS","CR","AD")
longnames = c("Rating","Complaints","Privileges","Learning","Raises","Critical","Advance")

### A. Figure 1 留意长短标签的用法
graph_ptsd.g=qgraph(data222.cor, labels=names, graph="pcor", layout="spring", vsize=6, cut=0, maximum=.99, sampleSize = nrow(dataw22),      border.width=1.5, border.color="black", minimum=.1,                  nodeNames =longnames,legend.cex=.4)
# dev.off()

### B. Figure 2 
centralityPlot(graph_ptsd.g,include = c('Strength','Closeness','Betweenness',"ExpectedInfluence"))
# dev.off()
centralityTable(graph_ptsd.g) #the details of centrality

### 稳健性检验
network1=estimateNetwork(dataw22, default = "pcor") #连续变量
#network1=estimateNetwork(dataw22, default = "EBICglasso") #顺序变量

boot1 = bootnet(network1, nBoots = 100,  nCores = 4) # 边估计检验
boot3 = bootnet(network1, nBoots = 100, type = 'case', nCores = 4) # 中心性估计检验

plot(boot1, labels = F, order = "sample" ) 
# dev.off()
plot(boot3)
#dev.off()
corStability(boot3) #输出CS系数，必须针对type = 'person'才行

1.2 分组复现例子二

  在一些情况中，研究者可能想要将比较不同组变量间的关系，此时就可以将所纳入的变量进行分组（如同一心理疾病的不同症状组），以此更直观地比较不同组变量间的联系。为了在定量基础上更好地探究不同组变量间的关联大小，常使用的一种分析方式是桥梁分析。

  桥梁分析由Jones等人在中心性指标的基础上提出，包括桥强度（bridge strength）、桥中介度（bridge betweenness）、桥亲近度（bridge closeness）和桥预期影响度（bridge expected influence）(Jones et al., 2021)。各指标定义类似与各中心性指标，不同之处在于，中心性指标衡量的是某节点与所有其他节点的关系，而桥梁中心性是基于节点分组，分析某一组（A组）内某一节点N1与另一组（B组）内所有节点的关系，而不计算N1与同组内其他节点的关联。因此，某一节点的桥中心性越高（如N1），表示它可能是在网络中起到连接它所在组与其他组关联的作用，即A组的其他节点通过N1与其他组的节点产生联系，这样的节点被称作“桥梁节点”。类似地，在桥梁分析中，桥中心性和桥EI表现出更好的稳定（Jones et al., 2021）。在心理领域中，桥梁分析可以用来分析两组心理学变量间的关联。以抑郁症状和焦虑症状的共患病研究为例，桥梁分析可以标记出能够解释两者共患病的症状（Kaiser et al., 2021），为理解两病共患病的潜在机理起到了提示作用的同时，也能起到揭示潜在干预方式的作用。

  研究者们可能会好奇不同网络结构是否有明显差异，以此探究某些变量在两个不同群体间的差异。限于篇幅，此处不针对纵向数据的网络分析进行过多的阐释和举例，感兴趣的读者也可以自行阅读这一参考文献（Borsboom et al., 2021），获得更全面的介绍。代码和例题可以见下方。

rm(list = ls()) #痕迹清除
#安装必要的包
library(qgraph) # Graph Plotting Methods, Psychometric Data Visualization and Graphical Model Estimation
library(NetworkComparisonTest) # Statistical Comparison of Two Networks Based on Several Invariance Measures
library(networktools) # Tools for Identifying Important Nodes in Networks
data('attitude')

#分组比较，首先创造两个组；obs做成两个组;变量也分了两个组
dataw22=as.data.frame(attitude[c(1:29),]) #选择你的待分析数据
dataw33=as.data.frame(attitude[c(2:30),]) #选择你的待分析数据

###把数据变成矩阵形式
data222.cor=cor_auto(dataw22,detectOrdinal =FALSE) 
data333.cor=cor_auto(dataw33,detectOrdinal =FALSE)

###命名节点/变量
names=c("RA","CO","PR","LE","RS","CR","AD")
longnames = c("Rating","Complaints","Privileges","Learning","Raises","Critical","Advance")
gr=list("TYPEA"=c(1:3), "TYEPB"=c(4:7))  #变量做成两个组

### A. Figure 1##
graph_ptsd.g = qgraph(
  data222.cor,
  labels = names,
  graph = "pcor",
  layout = "spring",
  vsize = 6,
  cut = 0,
  maximum = .99,
  sampleSize = nrow(dataw22),
  border.width = 1.5,
  border.color = "black",
  minimum = .1,
  groups = gr,
  color = c('#a8e6cf', '#dcedc1'),
  nodeNames = longnames,
  legend.cex = .4
)
#dev.off()
### B. Figure 2##
graph_ptsd.g <- qgraph(
  data333.cor,
  labels = names,
  graph = "pcor",
  layout = "spring",
  vsize = 6,
  cut = 0,
  maximum = .99,
  sampleSize = nrow(dataw33),
  border.width = 1.5,
  border.color = "black",
  minimum = .1,
  groups = gr,
  color = c('#a8e6cf', '#dcedc1'),
  nodeNames = longnames,
  legend.cex = .4
)

#dev.off()

#桥梁分析部分
b2=data333.cor
colnames(b2)=c("RA","CO","PR","LE","RS","CR","AD")
a2=bridge(b2,communities = list("TYPEA"=c(1:3),"TYEPB"=c(4:7))) #计算桥梁中心性;这里的意思是b2的这个亚组（obs）区分；其中的每个变量和本组其他变量的连接重要性；不是两个网络，而是同一个网络中变量分组的所谓桥梁性质

plot(a2,include = c("Bridge Strength","Bridge Expected Influence (1-step)"),zscore = T,order="value")
#桥梁中心性图
#dev.off()

#网络比较分析；这里就是两个网络图像不像的问题了；显著就代表存在显著差异的意思，两个图的中心性和边的问题
colnames(dataw22)=c("RA","CO","PR","LE","RS","CR","AD")
colnames(dataw33)=c("RA","CO","PR","LE","RS","CR","AD")

com1 = NCT(
  dataw22,
  dataw33,
  0.5,
  it = 1000,
  binary.data = FALSE,
  paired = FALSE,
  weighted = TRUE,
  AND = FALSE,
  abs = FALSE,
  test.edges = TRUE,
  edges = "all",
  test.centrality = TRUE,
  progressbar = TRUE,
  make.positive.definite = TRUE,
  verbose = TRUE
)
com1

  以上两个例题作者认为基本的应付简单论文已经够用了。配有不算太详细的视频讲解备查！

1.3 原理补充解读

  网络分析（Network analysis）正是解决此类问题的优势方法，近年来正在不断受到更广泛的关注，并且正在逐步成为代替潜变量建模的优势方法。举例而言，在使用GAD-7测量焦虑症状时，将各个测量条目视为一个焦虑症状放入网络分析进行关联呈现，最后可以解析出最核心的焦虑症状，解析各个症状之间的复杂关联。网络分析以图论作为理论基础，用可视化方式呈现出统计结果，可以更加直观地了解各变量之间的关系。

  网络结构图包括不同的类型，主要是涉及选择不同的矩阵类型。一般来说，包括两种类型:相关矩阵网络(Correlation networks)和偏相关网络(Partial correlation networks)。前面的例子给定的就是偏相关的矩阵作图和比较计算。

  具体的语法不复杂，查阅相应的代码帮助文件即可。除了qgraph包之外，还可以使用“bootnet”包的“estimateNetwork”函数进行网络分析估计

Network<- estimateNetwork(data, default = "cor")#cor是矩阵类型为相关矩阵，pcor偏向关矩阵，glasso为自适应Lasso算法矩阵

2 调节的网络分析

2.1 modnets包实现过程

  基于前面4，研究大部分仅探讨了节点之间的两两关系，但对于更复杂机制的考察很少，例如调节作用。并且，在考察调节作用的网络分析研究中，相当一部分仅仅是简单的将整个样本分为几个子样本（例如分为男性和女性）分别进行网络分析。因此，为了更准确的估计调节变量对节点网络的作用，有研究者开发了有调节的网络分析技术。

  众所周知，采用GGM方法估计无向网络模型是构建偏相关矩阵。然而，当模型中包含调节因子、协变量或任何高阶交互项时，这种方法就不合适了。因此，可以使用另一种称为节点智能回归的方法(Epskamp et al., 2018)，该方法基于一种称为邻域选择(neighborhood selection)的结构学习图论方法，并需要通过一系列单变量回归模型来估计网络结构。也就是说，为网络中为每个节点构建一个单独的回归模型，然后跨模型聚合系数以形成最终的网络结构。

rm(list = ls()) #痕迹清除
library(modnets) # Modeling Moderated Networks
skim(ggmDat)

fit1=fitNetwork(ggmDat)
plotNet(fit1) #plot出错，缺乏对象
fit1$adjMat
fit2=estimateNetwork(ggmDat,default = 'ggmModSelect')
plot(fit2)
fit2$graph
#可以看到，两种方法所得到的网络是相似的

fit3=fitNetwork(ggmDat,'M') #其中’M’为调节变量的变量名
plotNet(fit3)
#由于我们使用的是默认设置，该网络采用的是“OR”规则
fit3.1=fitNetwork(ggmDat,'M',rule = 'AND') #AND规则
plotNet(fit3.1)

#将调节变量在图中显示出来,同时显示边的系数。
plotNet(fit3,mnet = TRUE,elabs = T)

#加上将每个节点相关的预测误差在节点周围绘制为饼状图（不能和mnet = TRUE同时使用）
plotNet(fit3,predict = TRUE,elabs = T,con='R2') #饼子代表R2

#估计两个节点之间单向预测关系的调节边际效应图
library(arm) # Data Analysis Using Regression and Multilevel/HierarchicalModels
condPlot(fit3,from='V4',to='V5')

#估计节点的中心性指标
centPlot(fit3,scale = c("z-scores"))
centTable(fit3,scale = T)

#边的稳定性/准确性估计和差异检验
boot1=modnets::bootNet(ggmDat,'M',nboot=100)
plotBoot(boot1) #修订了推文命令左图代表交互作用的边，图右代表成对节点的边

boot2=modnets::bootNet(ggmDat,'M',nboot=100,caseDrop = T)
plotBoot(boot2) #进一步通过caseDrop bootstrap方法估计边的稳定性

plotBoot(boot1,difference = T) #检验边之间是否存在差异

#节点的稳定性/准确性估计

plotBoot(boot1,type='strength')
plotBoot(boot2,type='strength') #进一步通过caseDrop bootstrap方法估计节点的稳定性
plotBoot(boot1,difference = T,type='strength') #检验节点中心性之间是否存在差异

2.2 mgm包实现过程

  以2022年在Applied Psychology: Health and Well-Being上发表的文章《Changes in the network structure of mental health after a multicomponent positive psychology intervention in adolescents: A moderated network analysis》为例(Tejada‐Gallardo et al., 2022)，详细的代码和数据已经下载下载；基于下载的代码来进行复现和理解。

  本研究的结果分为三个数据集DataPre、DataPost和DataFu，分别代表三个时间点所测得的数据。每个数据集包含7个变量。

rm(list = ls()) #痕迹清除
## Changes in the structure of mental health after a multicomponent positive psychology 
##intervention in adolescents: A longitudinal network analysis

#setwd("~/...")
library("foreign") # Read Data Stored by 'Minitab', 'S', 'SAS', 'SPSS', 'Stata','Systat', 'Weka', 'dBase', ...
library(mgm) # Estimating Time-Varying k-Order Mixed Graphical Models
library(psych) # Procedures for Psychological, Psychometric, and Personality Research
library(qgraph) # Graph Plotting Methods, Psychometric Data Visualization and Graphical Model Estimation
library(bootnet) # Bootstrap Methods for Various Network Estimation Routines


##Load Data##
DataPre <-read.spss('/Users/mengtaogao/Library/CloudStorage/OneDrive-共享的库-onedrive/outsyn/station/Effects of an MPPI on mental health/Pre_MNM.sav', to.data.frame=T)
DataPost <-read.spss('/Users/mengtaogao/Library/CloudStorage/OneDrive-共享的库-onedrive/outsyn/station/Effects of an MPPI on mental health/Post_MNM.sav', to.data.frame=T)
DataFu <-read.spss('/Users/mengtaogao/Library/CloudStorage/OneDrive-共享的库-onedrive/outsyn/station/Effects of an MPPI on mental health/Fu_MNM.sav', to.data.frame=T)

# class(DataPre) #检查一下数据框架包括view
###Network models##
Net_Pre<-estimateNetwork(DataPre, default="ggmModSelect", stepwise=TRUE)
Net_Post<-estimateNetwork(DataPost, default="ggmModSelect", stepwise=TRUE)
Net_Fu<-estimateNetwork(DataFu, default="ggmModSelect", stepwise=TRUE)

plot(Net_Pre)
plot(Net_Post)
plot(Net_Fu)
#从左到右分别是三个时间点横断网络，可以发现调节变量与任何节点均无显著相关
#将分组作为调节变量，估计MNMs;常规操作运行正常暂时略过
# centralityPlot(Net_Fu, include=c("Strength", "ExpectedInfluence"), scale=c("z-score"))

#weight-edge accuracy
# boot1<-bootnet(Net_Fu, nBoots=100, type="parametric", nCores=2)
# plot(boot1, order = "sample")
# plot(boot1)
# dev.off()

##centrality stability##
# boot2 <- bootnet(Net_Pre, nBoots = 100, type = "case", nCores = 2)
# plot(boot2)
# corStability(boot2)

#Pre-test#
set.seed(30)
mgm_objPre <- mgm(
  data = as.matrix(DataPre),
  type = c("c", "g", "g", "g", "g", "g", "g"),
  level = c(2, 1, 1, 1, 1, 1, 1),
  lambdaSel = "CV",
  lambdaGam = .25,
  ruleReg = "AND",
  moderators = 1,
  scale = TRUE
)

mgm_objPre


#Pre-test#
l_mgm_cond <- list()
for(g in 1:2) l_mgm_cond[[g]] <- condition(object = mgm_objPre,
                                           values = list("1" = g))

v_max <- rep(NA, 2)
for(g in 1:2) v_max[g] <- max(l_mgm_cond[[g]]$pairwise$wadj)
par(mfrow=c(1, 2))
for(g in 1:2) {
  qgraph(input = l_mgm_cond[[g]]$pairwise$wadj, #plot from web tutorial
         edge.color = l_mgm_cond[[g]]$pairwise$edgecolor_cb,
         lty = l_mgm_cond[[g]]$pairwise$edge_lty,
         layout = "circle", mar = c(2, 3, 5, 3),
         title = "Time 1", labels = colnames(DataPre),
         maximum = max(v_max), vsize = 10, esize = 23,
         edge.labels  = TRUE, edge.label.cex = 3)
  mtext(text = paste0("Group", g), line = 2.5,
  )
}
# dev.off()

#Post-test#
set.seed(30)
mgm_objPost <- mgm(data = as.matrix(DataPost),
                   type = c("c", "g", "g", "g", "g", "g", "g"),
                   level = c(2, 1, 1, 1, 1, 1, 1),
                   lambdaSel = "CV",
                   lambdaGam =.25,
                   ruleReg = "AND",
                   moderators = 1,
                   scale = TRUE)

# Assess Stability，比较耗时
res_obj <- resample(object = mgm_objPost,
                    data = as.matrix(DataPost),
                    nB = 50,
                    pbar = TRUE)

plotRes(res_obj)

  稳定性左图显示两个节点之间的成对效应，右图显示调节效应。圆圈中的值表示估计给定参数为非零的bootstrap样本的比例，黑色线条的两段表示5%和95%bootstra区间。可以发现，在不同时间下，成对和调节效应都是相当稳定的。调节效应在时间点2时明显更强，因为有几个值远离零。

  感觉这个部分的原理还不是太清楚。

3 时间网络分析

  纵向数据可以分为两种类型：发展过程数据和稳定过程数据。这两种数据是通过不同的研究设计和测量方法得到的，两者的主要差别可以见下表。针对2~3次追踪的面板数据，目前常用交叉滞后网络分析。而针对5次以上追踪的面板数据和密集纵向数据则常用时间网络分析。

  主要参考推文：时间网络分析(Temporal Network Analysis)入门教程。在这里我们使用《A Temporal Network Approach to Paranoia: A Pilot Study》中的数据为例(Contreras et al., 2020)。时间网络分析方法建立于“向量自回归”(Vector autoregression, VAR)模型之上。简单来说，向量自回归模型中的“向量”指的是多元特征，是单变量自回归模型的延伸。而向量自回归模型由一组回归方程组成，其中所有变量都被视为内生变量。因此，它们既可以作为自变量，也可以作为因变量。所以，向量自回归分析可以在没有关于变量之间关联方向的先验假设的情况下进行。通过这种方法，我们就可以知道多个变量之间的跨时间动态关系。

  在R中，我们使用mlVAR包实现时间网络分析。在R中，我们使用mlVAR包实现时间网络分析。

rm(list = ls()) #痕迹清除
# Load packages:
#---------------                         
library("mlVAR")
library("graphicalVAR")
library("qgraph")
library(foreign)      
library(qgraph)
library(networktools) 
library(corpcor)      
library(tseries)      
library(TSA)          
library(dplyr)

  dataset<- read.spss('/Users/mengtaogao/Library/CloudStorage/OneDrive-共享的库-onedrive/outsyn/station/A Temporal Network Approach to Paranoia/Temporal network approach to paranoia_data.sav',to.data.frame=TRUE, use.value.labels = F)
  baseline<-dataset[,c(1:4,10:13)]
  
  #Preparation:
#------------
  names(baseline) <- c("subject","day","beep","sad","others","EA","SE","paranoia")
  vars <- c ("sad","others","EA","SE","paranoia")
  idvar <- "subject"
  dayvar<-"day"
  beepvar<-"beep"
  
  #Estimation: two-step multilevel Vector Autoregressive (mlVAR):
#--------------------------------------------------------------
baseline_mlVAR <- mlVAR(baseline, vars, idvar, lags = 1, estimator = "lmer",temporal = "correlated", contemporaneous = "correlated")
  
  #Visualization:
#--------------
  
  #Separated plot: 
contemporaneous<-plot(baseline_mlVAR, "contemporaneous", layout = "circle", nonsig = "hide", theme = "colorblind", title = "Contemporaneous",
labels = vars, vsize = 13, rule = "and", asize = 10, mar = c(5,5,5,5))

temporal<-plot(baseline_mlVAR, "temporal", layout = "circle", nonsig = "hide", theme = "colorblind", title = "Temporal",labels = vars, vsize = 13, rule = "and", asize = 5, mar = c(5,5,5,5))

between_subjects<-plot(baseline_mlVAR, "between", layout = "circle", nonsig = "hide",  theme = "colorblind", title = "Between-subjects",labels = vars, vsize = 13, rule = "and", asize = 10, mar = c(5,5,5,5))

#节点中心性估计
centralityPlot(contemporaneous,include = c('Strength','Closeness','Betweenness',"ExpectedInfluence"))
centrality(contemporaneous)
centralityPlot(temporal,include = c('Strength','Closeness','Betweenness',"ExpectedInfluence"))
centrality(temporal)
centralityPlot(between_subjects,include = c('Strength','Closeness','Betweenness',"ExpectedInfluence"))
centrality(between_subjects)

• 时间网络（Temporal network）。时间网络中边表示在控制了前一个测量点上的所有其他变量后，前一个测量点对下一个测量点的节点的预测能力。例如图中我们可以发现，悲伤、与他人的亲密度、体验回避、自尊、偏执信念都存在自回归效应。与他人的亲密度消极地预测了下一个时间点的偏执，悲伤积极地预测了体验回避，而体验回避反过来又积极地预测了下一个时间点的自尊;其实就是within估计类似。

• 同时网络（Contemporaneous network）。同时网络中的边表示在控制了时间关联之后，同一时间框架内变量之间的关联。可以发现有几个联系非常突出，偏执信念和悲伤呈正相关，这表明在同一时间段内，偏执信念得分越高，悲伤水平越高。同样，悲伤与自尊和与他人的亲密度都是负相关。换句话说，悲伤程度越高，自尊水平越低，同时与他人的亲密感也越低。此外，自尊和与他人的亲密度是正相关的。最后，体验回避与任何其他变量无关;

• 个体间网络（Between-subject network）。个体间网络模型边表示，在考虑到网络中的其余变量后，个体内部平均水平之间的相关性。可以发现，与他人的亲密度的平均水平与偏执信念的平均水平呈负相关，与自尊的平均水平呈正相关。悲伤的平均水平与自尊的平均水平呈负相关，与体验回避的平均水平呈正相关。其实就是Between估计量；

4 针对面板数据的时间网络分析(Temporal Network Analysis)

  介绍针对面板数据的时间网络分析(Temporal Network Analysis)。需要注意，这种纵向网络分析所需要的面板数据至少要有三个wave。并且，面板数据的时间网络分析与密集追踪数据的时间网络分析相似，同样区分了个体内和个体间效应，并可以估计时间网络、同时网络和个体间网络。

  命令主要是：mlGVAR {modnets}:Fit GVAR models with multilevel data Fits a graphical vector autoregressive model to data containing multiple time points measured for multiple individuals.

  纵向网络分析主要通过图形向量自回归模型(the graphical vector-autoregression model, GVAR)实现(Epskamp et al., 2018; Wild et al., 2010)。1阶滞后(lag-1) GVAR模型是由单主体时间序列的高斯图形模型(the Gaussian graphical model, GGM)推广而来。GGM模型大家都很熟悉，在横断网络分析中被用来估计部分相关，形成无向网络模型。而在(lag-1) GVAR模型中，时间依赖关系通过对先前测量情况的回归来建模，这导致回归系数矩阵也可用于绘制有向网络模型(Bringmann et al., 2013)，通常被称为时间网络。剩余的方差和协方差(即控制前一次测量场合后的协方差结构)可以建模为一个GGM模型，也称为同时网络。当多个主体的时间序列可用时，可以形成关于个体间效应(稳定均值之间的关系)的第三个GGM模型，也称为个体间网络。具体统计原理大家可以参考(Epskamp, 2020)。

  简单通过例题数据来进行演示。按照这里的面板逻辑，仅仅fixed和between是有用的，时序上的ar过程其实没有太大意义。

rm(list = ls()) #痕迹清除
library(modnets)
head(mlgvarDat)
fit1 <- mlGVAR(mlgvarDat, 'M',idvar = "ID") #Estimating fixed networks
(p1=plotNet(fit1,mnet = TRUE,elabs = T))
fit2 <- mlGVAR(mlgvarDat, 'M', bm = TRUE) # Fit the same moderator in the between-subjects network
(p2=plotNet(fit2,mnet = TRUE,elabs = T))
fit3=mlGVAR(mlgvarDat, 'M', selectFUN = 'varSelect')
(p3=plotNet(fit3,mnet = TRUE,elabs = T))

#极为耗时的一些计算过程
(x <- lmerVAR(mlgvarDat, 'M', temporal = "fixed", contemp = "orthogonal"))
fit1 <- lmerVAR(mlgvarDat, temporal = "fixed", contemp = "orthogonal")
fit2 <- lmerVAR(mlgvarDat, temporal = "orthogonal", contemp = "orthogonal")

compareVAR(fit1, fit2)

  目前，已有一些研究使用了基于面板数据的纵向网络分析，例如有研究考察了有问题的网络使用、心理困扰和饮食失调症状之间的关系(Lin et al., 2023)。但这种方法相比于横断网络分析，在评估网络和节点的中心性、稳定性、可重复性等方面存在不足，有待研究者们的进一步开发。2024年五月份的意见：基于面板数据的纵向网络分析-时间网络分析(Temporal Network Analysis)入门教程（2）。

5 探索性图论分析

  根据推文：探索性图论分析；这里原理进行了讲解，但是数据没有提供，直接采用?EGAnet的例题自带数据进行处理就足够了。

  探索性图论分析 (Exploratory Graph Analysis) 是近年来兴起的一种维度检测方法。这种方法有别于之前的探索性因子分析 (Exploratory Factor Analysis)，它将变量构建为多元正态分布网络，并将显性变量视为由边线（Edges）连接的节点（Nodes）所构成的网络，而因子模型中的潜在变量则使这些节点相互聚集。

  通过将高斯图模型（Gaussian graphical model，GGM）与加权网络的聚类算法相结合，EGA 可以对心理构念中条目维度进行评估。模拟研究显示，与作为评估维度“金标准”的因素分析法相比，EGA 在准确检测心理量表的正确维度方面表现更为出色(Golino et al., 2020)。EGA 以其以下优点而备受关注：（1）可视化能力：可以通过使用不同颜色编码的网络图直观地判断条目维度的归属；（2）自动估计：无需考虑采用何种旋转方法；（3）稳定性检验：可通过 bootstrap 方法进行维度稳定性检验。因此，Cosemans等人(2022)将EGA作为估计连续数据维度的首选方法。

rm(list = ls()) #痕迹清除
library(EGAnet) # Exploratory Graph Analysis – a Framework for Estimating theNumber of Dimensions in Multivariate Data using Network Psychometrics

# Get data
data = psych::bfi[,1:25]

## Not run: # Compute LCT
## Factor model;默认it=100
LCT(data,seed=123)
#结果表明，随机抽取的子样本总，有多少适合网络分析

#Bootstrap Exploratory Graph Analysis
# Load data
wmt = wmt2[,7:24]

## Not run: 
# Standard EGA parametric example
boot.wmt = bootEGA(
  data = wmt, iter = 500,
  type = "parametric", ncores = 2,seed=123
) #项目的稳定性应该大于0.65
plot(boot.wmt$EGA)

boot.wmt$summary.table #得到的因子中位数与其他统计指标
boot.wmt$frequency  #得到的因子数量估计分布
dimensionStability(boot.wmt) #维数的稳定性应该大于0.75

net.loads(boot.wmt$EGA) #载荷为0.15表示非常小

  结合item stability和network loading对各个项目进行检查，当itme stability小于0.65，network loading也较小时，则考虑删除该项目。值得注意的是，项目的删除与否需要考虑多方面的原因，例如理论层面上，该项目是否反映了核心变量；实践层面上，该项目是否是主要关注的变量。因此，项目的删除与否需要结合文献具体分析。

6 网络分析中控制协变量

  主要根据推文：如何在网络分析中控制协变量？和推文上的两篇论文来进行实操和理解。目前主要有两种方式控制协变量：一种直接把每个节点回归掉协变量后取残差，另一种将协变量也当做节点纳入网络模型。

  两种方法各有利弊，第一种的问题在于节点是基于总样本的残差，在后续bootstrap时会引入误差（另外，有时候，使用这种方法与不控制协变量相比，在进行网络中心性估计时，稳定性和准确性会相对差一些）；第二种的问题在于协变量的引入会造成网络的结构变化。采用的时候要说清楚利弊何在。其实一般来说，第二种方法实施起来似乎简单。从实际运用的角度，采用第二种方法更为便利一些。

  示例数据来自于networktools中的depression数据集，这个模拟数据集包含1000个人中9种重度抑郁症症状（PHQ-9）的严重程度评级。症状评级假设是在100点滑动量表上自我报告的。

rm(list = ls()) #痕迹清除
#3.1载入包与确定代码路径（载入出现错误的话需要先去安装包）----
library(networktools) # Tools for Identifying Important Nodes in Networks
library(bootnet) # Bootstrap Methods for Various Network Estimation Routines
library(mgm) # Estimating Time-Varying k-Order Mixed Graphical Models

#3.2 数据和网络估计准备------
data = networktools::depression #用depression数据作为本教程的示例数据;指明究竟哪个包
#生成协变量并组合成新数据集

gender = as.numeric(sample(x = c(1, 2), size = nrow(data), replace = TRUE))
# 随机生成1000个性别数据
age = sample(x = 20:30, size = nrow(data), replace = TRUE) #随机生成1000个年龄数据
usedata = data.frame(data, age,gender) # 将生成的数据合并成一个数据框

# 3.3 协变量网络估计
##3.3.1 第一种方法
###回归掉协变量

depression_mat=as.data.frame(matrix(nrow=1000,ncol=9)) #做出一个空的数据库来存储残差
colnames(depression_mat)=colnames(usedata)[1:9] #空数据库变量命名
#这个函数的意义是将每一个变量进行回归，然后将残差存储在depression_mat中
for (i in 1:(ncol(depression_mat))) { #留意一下循环的操作
  formula <- paste(colnames(usedata)[i], '~age+gender')
  lm1 <- glm(formula = formula,family = "gaussian", data = usedata) 
  y_res <- lm1$residuals
  depression_mat[, i] <- y_res
}

head(depression_mat) #检查归结残差的数据

network = as.matrix(depression_mat)
p=ncol(network)
fit_obj = mgm(data = network, type = rep('g', p), level = rep(1, p), lambdaSel = 'CV', ruleReg = 'OR', pbar = TRUE)
pred_obj = predict(object = fit_obj, data = network, errorCon = 'R2')
pred_obj$error #可预测性计算结果，其实就是得到全部的R2而已

#估计网络
Network = estimateNetwork(data=network, default = "EBICglasso", tuning = 0.5, corMethod="cor",corArgs=list(method="spearman",use="pairwise.complete.obs"))
groups=list("PHQ-9 symptoms"=c(1:9))
colnames(usedata)=c("D1","D2","D3","D4","D5","D6","D7","D8","D9")
items =list("Angedonia","Sadmood","Sleep","Fatigue","Appeitite","Worthless","Concentration","Montor","Death")

#网络图
plot_network <- plot(
  Network,
  layout = "spring",
  groups = groups,
  label.cex = 1,
  #调整节点中标签的字体大小
  label.color = 'black',
  # 调整节点中标签的字体颜色
  negDashed = T,
  #负相关的边为虚线
  legend = T,
  nodeNames = items,
  legend.cex = 0.55,
  #图例的字体大小
  legend.mode = 'style1',
  #图例的风格，默认是'style1'
  pie = pred_obj$error[, 2]
  #节点的环形图显示的是R2
)

#3.3.2 第二种方法
groups=list("PHQ-9 symptoms"=c(1:9),"Covariates"=c(10:11)) #协变量一组
colnames(usedata) =c("D1","D2","D3","D4","D5","D6","D7","D8","D9","C1","C2")
items = list("Angedonia","Sad mood","Sleep","Fatigue","Appeitite","Worthless","Concentration","Montor","Death","age","gender")
network = as.matrix(usedata) #整个运行的网络框架11个变量，包含两个协变量
p=ncol(network)

fit_obj = mgm(
  data = network,
  type = rep('g', p),
  #这里将性别也作为连续变量进行处理
  level = rep(1, p),
  lambdaSel = 'CV',
  ruleReg = 'OR',
  pbar = TRUE
)
pred_obj = predict(object = fit_obj, data = network, errorCon = 'R2')
pred_obj$error#可预测性计算结果

#估计网络
Network = estimateNetwork(
  data = network,
  default = "EBICglasso",
  tuning = 0.5,
  corMethod = "cor",
  corArgs = list(method = "spearman", use = "pairwise.complete.obs")
)

#网络图
plot_network = plot(
  Network,
  layout = "spring",
  groups = groups,
  label.cex = 1,
  #调整节点中标签的字体大小
  label.color = 'black',
  # 调整节点中标签的字体颜色
  negDashed = T,
  #负相关的边为虚线
  legend = T,
  nodeNames = items,
  #图例中的节点标签名
  legend.cex = 0.55,
  #图例的字体大小
  legend.mode = 'style1',
  #图例的风格，默认是'style1'
  pie = pred_obj$error[, 2]
)

7 贝叶斯网络与因果推断之类

  例子和演示的第一篇推文：一文搞懂贝叶斯网络模型及其R实现；贝叶斯网络（Bayesian Network, BN）是一种概率图模型，它通过有向无环图（Directed Acyclic Graph, DAG）表示一组随机变量及其相互之间的条件依赖性。在处理不确定性和复杂性问题时，贝叶斯网络提供了一种直观且强大的工具。

7.1 贝叶斯网络的构造

  贝叶斯网络由两部分组成：网络结构和条件概率表（Conditional Probability Table, CPT）。

• 网络结构：网络结构是指用DAG表示的变量间的拓扑关系。在DAG中，节点代表随机变量，有向边（从A指向B）表示变量A是变量B的父节点，即A对B有直接影响。根据DAG，我们可以得出变量间的条件独立性，即在给定父节点的情况下，一个节点的条件概率分布仅依赖于其父节点。

• 条件概率表（CPT）：对于网络中的每个节点，我们需要定义一个CPT，它描述了该节点在其父节点取特定值时的概率分布。CPT是局部的，因为它只依赖于父节点的状态。

7.2 贝叶斯网络模型的学习分为两部分

• 结构学习（Structure Learning）： 结构学习是指从数据中学习变量之间连接的拓扑结构。这个过程通常涉及以下步骤：

  • 变量选择：选择需要建模的变量集合。

  • 条件独立性检验：利用条件独立性原则来确定哪些变量之间需要建立连接。条件独立性是指在给定某些变量的值的情况下，其他变量的值是相互独立的。

  • 评分和搜索：使用评分函数（如贝叶斯信息准则BIC、阿卡伊克信息准则AIC等）来评估不同结构的好坏，并结合搜索算法（如贪婪等价搜索、模拟退火等）来找到最优的结构。

• 参数学习（Parameter Learning）： 参数学习是指在给定结构的情况下，估计条件概率表（CPTs）中的参数。这个过程通常包括：

  • 条件概率分布的确定：对于每个节点，确定其在父节点取不同值的情况下的条件概率分布。

  • 参数估计：使用最大似然估计（MLE）、贝叶斯估计等方法，根据数据来估计这些条件概率。贝叶斯网络的构造和参数学习是一个相互关联的过程。

结构学习提供了变量之间关系的框架，而参数学习则填充了这个框架中的具体内容。通过这两个步骤，我们可以得到一个能够对数据进行有效建模和推理的贝叶斯网络。

  介绍四类主要的网络结构学习方法：基于约束、基于分数、配对和混合方法，以及他们的优缺点和比较网络结构学习方法的一些基本原则。

  bnlearn是R中一个强大的包，用于学习、推理、模拟和分析贝叶斯网络。它提供了一系列工具来处理贝叶斯网络的不同方面，包括网络结构的学习、参数的估计、条件概率的查询、数据集的模拟，以及多种网络的图形表示方法。以下是如何使用bnlearn包构建贝叶斯网络的详细介绍。 本次使用的数据集是：BreastCancer数据集是著名的机器学习和数据挖掘领域的标准数据集之一，通常用于分类和预测任务。这个数据集包含了来自乳腺癌患者的各种特征，旨在通过这些特征预测肿瘤是良性还是恶性。

rm(list = ls()) #痕迹清除
# 加载bnlearn包
library(bnlearn) # Bayesian Network Structure Learning, Parameter Learning and Inference
# 导入BreastCancer数据集
library(mlbench) # Machine Learning Benchmark Problems

data(BreastCancer)
# 查看变量类型
str(BreastCancer)
# 查看摘要统计信息
skim(BreastCancer)

# 将数据集转换为离散变量
# 在构建贝叶斯网络模型之前，我们需要将连续变量转换为离散变量。可以使用discretize()函数将连续变量转换为离散变量。以下代码将数据集中的连续变量转换为离散变量：
# library(bnlearn)
data.bc <- discretize(BreastCancer[,2:ncol(BreastCancer)], method = "interval", breaks = "frequency")

# 数据预处理# 使用数据集
data.bc <- BreastCancer[,2:ncol(BreastCancer)]# 使用VIM::aggr()查看函数的缺失值
library(VIM) # Visualization and Imputation of Missing Values
aggr(data.bc) #缺失值的情况的可视化
# 使用随机森林填补缺失值
library(missForest) # Nonparametric Missing Value Imputation using Random Forest
data.bc_new <- missForest(data.bc,ntree = 100)
data.bc_new<- data.bc_new$ximp


# ① 使用希尔-克莱姆算法学习网络结构：
bn.structure <- hc(data.bc_new, score = "bic")
plot(bn.structure) # 绘制学习到的结构

# ② 使用基于约束的PC算法：
bn.structure <- pc.stable(data.bc_new,test="mi")
plot(bn.structure) # 绘制学习到的结构

# 1. 基于分数的方法# 贝叶斯网络结构的禁忌搜索（Tabu Search）
bn.tabu <- tabu(data.bc_new)
plot(bn.tabu) # 绘制学习到的结构

# 2. 基于约束的方法# Grow-Shrink算法（GS）
bn.gs <- gs(data.bc_new)
plot(bn.gs) # 绘制学习到的结构

# 3. 配对方法# 最大权重生成树算法（Chow-Liu）
bn.chowliu <- chow.liu(data.bc_new)
plot(bn.chowliu) # 绘制学习到的结构

#4. 混合方法#Max-Min Hill-Climbing算法（MMHC）
bn.mmhc <- mmhc(data.bc_new)
plot(bn.mmhc ) # 绘制学习到的结构

# 参数学习# 结构学习完成后，下一步是参数学习。参数学习是指学习网络中每个变量的条件概率分布（对于离散变量）或条件概率密度（对于连续变量）。
bn.parameters <- bn.fit(bn.structure, data.bc_new) #结构选择基于约束
bn.parameters

# 绘制网络
# graphviz.plot(bn.structure) #结构选择基于约束

# 计算在给定某些证据下肿瘤是良性的概率# 假设我们知道 Bare.nuclei = "1" 和 Marg.adhesion = "1"
prob_benign_given_evidence <- cpquery(fitted = bn.parameters, event = (Class == "benign"),evidence = ( Marg.adhesion == "1"))

prob_benign_given_evidence

# 生成给定证据下Class变量的条件概率分布
distribution_class_given_evidence <- cpdist(
  fitted = bn.parameters,
  nodes = "Class",
  evidence = (Bare.nuclei == "1")
)
tabyl(distribution_class_given_evidence$Class)
# 比较模型,这将输出一个列表，其中包括两个网络共有的边、只在第一个网络中出现的边、以及只在第二个网络中出现的边。
compare(bn.tabu, bn.gs)

# 对Tabu搜索算法进行引导稳定性分析
bs.tabu <- boot.strength(data = data.bc_new, algorithm = "tabu", R = 100)
bs.tabu

#Rgraphviz没有安装，无法可视化

# 条件独立性测试# 条件独立性测试是评估变量间是否独立的一种方法，对于基于约束的网络学习方法尤其重要。bnlearn提供了ci.test函数来执行条件独立性测试。以下是如何在变量之间进行条件独立性测试的示例：
ci.test("Marg.adhesion", "Bare.nuclei", data = data.bc_new)

7.3 横断网络中的因果推断

  基于推文：贝叶斯网络分析入门-横断网络中的因果推断下载了数据。以此为例子。doi顺手把参考文献添加。

rm(list = ls()) #痕迹清除
# A tutorial on Bayesian networks for psychopathology researchers
# Briganti G, Scutari M, McNally RJ

library(stats)
library(qgraph)
library(readr)
library(bootnet)
library(dplyr)
library(corpcor)
library(bnlearn)
library(psych)
library (ggplot2)

data <- read_delim('/Users/mengtaogao/Library/CloudStorage/OneDrive-共享的库-onedrive/outsyn/station/A tutorial on Bayesian Networks/data.csv', ";", escape_double = FALSE, trim_ws = TRUE)

names <- c("DepMood", 
           "Sleep", 
           "Weight", 
           "Fatigue", 
           "Irritable", 
           "SuicideId", 
           "Anhedonia")

longnames <- c("Depressive Mood", 
               "Sleep Problems", 
               "Weight", 
               "Fatigue", 
               "Irritable", 
               "Suicidal Ideation", 
               "Anhedonia")

data <- as.data.frame(matrix(as.numeric(as.matrix(data)), 
                             ncol=ncol(data), 
                             byrow=TRUE))  #set data as numerics
# PC algorithm

BNpc<-pc.stable(data)
qgraph(BNpc, labels=names, 
       nodeNames=longnames, 
       legend.cex=.5, 
       vsize=6, 
       layout="circle")

# Boot stability for the DAG
# Strength: connection strength, e.g. 0.85 > 
#              connection appears in 86% of the fitted networks.
# Direction: probability of the direction
#     e.g. 0.57 means that in 57% of the fitted networks the connection goes in 
#                the direction depicted in the graph.

BST <- boot.strength(data, 
                     R = 1000, 
                     algorithm = "pc.stable") 
head(BST)
qgraph(BST) # visualize output with qgraph
BST[BST$strength > 0.85 & BST$direction > 0.5, ]

avgnet1 <- averaged.network(BST, threshold = 0.85)
astr1 <- arc.strength(avgnet1, data, "bic-g")  # compute edge strengths
qgraph(
  astr1,
  labels = names,
  nodeNames = longnames,
  legend.cex = .5,
  vsize = 6,
  layout = "circle"
)

#公众号上的演示不完整

# Hill Climbing algorithm

# Boot stability for the DAG
# Strength: connection strength, e.g. 0.85 > 
#              connection appears in 86% of the fitted networks.
# Direction: probability of the direction
#     e.g. 0.57 means that in 57% of the fitted networks the connection goes in 
#                the direction depicted in the graph.

BST <- boot.strength(data, R = 1000, algorithm = "hc", debug = TRUE)
head(BST)
BST[BST$strength > 0.85 & BST$direction > 0.5, ]
qgraph(BST)

avgnet1 <- averaged.network(BST, threshold = 0.85)
avgnet1

bnlearn::score(avgnet1, data = data)

astr1 <- arc.strength(avgnet1, data, "bic-g")   ## compute edge strengths

qgraph(astr1, 
       labels=names, 
       nodeNames=longnames, 
       legend.cex=.5, 
       vsize=6, 
       layout="circle")


# Equivalence classes

BNpc<-pc.stable(data)
BNcp <- cpdag(BNpc)
qgraph(
  BNcp,
  labels = names,
  nodeNames = longnames,
  legend.cex = .5,
  vsize = 6,
  layout = "circle"
)

vstructs(BNcp)

BST <- boot.strength(
  data,
  R = 1000,
  algorithm = "pc.stable",
  debug = TRUE,
  cpdag = TRUE
)


head(BST)
BST[BST$strength > 0.85 & BST$direction > 0.5, ]
qgraph(BST)

avgnet1 <- averaged.network(BST, threshold = 0.85)
avgnet1

bnlearn::score(avgnet1, data = data)

astr1 <- arc.strength(avgnet1, data, "bic-g")   ## compute edge strengths

qgraph(
  astr1,
  labels = names,
  nodeNames = longnames,
  legend.cex = .5,
  vsize = 6,
  layout = "circle"
)

# Grow shrink
BNgs <- gs(data)
qgraph(
  BNgs,
  labels = names,
  nodeNames = longnames,
  legend.cex = .5,
  vsize = 6,
  layout = "circle"
)

BST <- boot.strength(data, R = 1000, 
                     algorithm = "gs", 
                     debug = TRUE)  
head(BST)
BST[BST$strength > 0.85 & BST$direction > 0.5, ]
qgraph(BST)

avgnet1 <- averaged.network(BST, threshold = 0.85)
avgnet1

bnlearn::score(avgnet1, data = data)

astr1 <- arc.strength(avgnet1, data, "bic-g")   ## compute edge strengths

qgraph(
  astr1,
  labels = names,
  nodeNames = longnames,
  legend.cex = .5,
  vsize = 6,
  layout = "circle"
)

# Restricted Maximization hybrid algorithm

BNrs <- rsmax2(data)
BST <- boot.strength(data, 
                     R = 1000, 
                     algorithm = "rsmax2", 
                     debug = TRUE)  
head(BST)
BST[BST$strength > 0.85 & BST$direction > 0.5, ]

avgnet1 <- averaged.network(BST, threshold = 0.85)
avgnet1

bnlearn::score(avgnet1, data = data)

astr1 <- arc.strength(avgnet1, data, "bic-g")   ## compute edge strengths

qgraph(
  astr1,
  labels = names,
  nodeNames = longnames,
  legend.cex = .5,
  vsize = 6,
  layout = "circle"
)

## GGM for comparison
n1 <- estimateNetwork(data, default="EBICglasso", threshold=TRUE)
plot(
  n1,
  layout = "circle",
  labels = names,
  nodeNames = longnames,
  theme = "colorblind",
  legend.cex = 0.5,
  vsize = 6
)

  原理与工具的一些补充。神器dagitty；已经安装了dagitty包。基于推文：（因果）贝叶斯网络与有向无环图（DAG）在因果推断中的应用。原理以这个为基础：贝叶斯网络在精神病理学研究中的应用。

• 贝叶斯网络的优势：

  • a.可视化依赖关系：通过图结构直观展示变量之间的依赖性和因果方向，有助于理解复杂交互机制。

  • b.处理复杂交互关系：适合多变量、多层次交互分析，尤其是涉及心理学和医学中的复杂系统。

  • c.缺失数据处理能力：贝叶斯网络能够在数据缺失的情况下进行推断，不需要完整数据集即可进行分析。

  • d.因果推断：不仅可以分析相关性，还能明确变量之间的因果关系，有助于决策和干预方案的设计。

  • e.灵活性和扩展性：支持动态贝叶斯网络（DBN）分析时间序列数据，进一步研究变量随时间变化的动态影响。

• 贝叶斯网络与横断网络分析的区别

  • 横断网络分析：分析变量间的相关性或联系强度；使用无向图表示变量间关系，不能区分方向性或因果性；适合描述变量之间的静态关联结构。

  • 贝叶斯网络分析：强调因果推断和方向依赖关系；使用有向无环图（DAG）明确因果路径，能够推断变量之间的因果链条；适合探索复杂交互和预测分析，尤其在心理学和医学领域。

  贝叶斯网络是分析复杂多变量关系的强大工具，尤其适用于心理学和医学领域的因果推断。它不仅可以揭示变量之间的条件依赖和方向性，还能有效处理缺失数据，并对变量之间的交互作用进行可视化和量化分析。

  此外，贝叶斯网络具有灵活性和扩展性，支持动态贝叶斯网络（DBN）处理时间序列数据，进一步探索随时间变化的动态过程。因此，它在精神病理学等领域具有重要的应用价值，为诊断、预测和干预策略提供了科学依据和决策支持。

8 传统意义的社会网络分析（基于obs的处理）

  前面的网络分析都是基于列的处理；传统的SNA是基于行的计算与可视化。这个部分就是补充一些关于行的计算。