library(glmnet)
library(qgraph)
library(lavaan)
library(psych)
library(bootnet)
setwd("D:/R/新数据")
R <- read.table(file = "交叉数据.txt", header = TRUE)
#View(R)（消除#以查看）
R$time <- as.factor(R$time)
R$RX <- as.factor(R$RX)
R$ID <- as.factor(R$ID)
R1 <- R[R$time == 1, ]
R2 <- R[R$time == 2, ]
nfac <- 1
GSE1 <- R1[,grep(colnames(R1), pattern = "gse")] #10个项目提取一个主成分
GSE2 <- R2[,grep(colnames(R2), pattern = "gse")]
pcmodel1 <- principal(GSE1, nfactors = nfac, scores = TRUE, rotate = "varimax")
pcmodel2 <- principal(GSE2, nfactors = nfac, scores = TRUE, rotate = "varimax")
plot(pcmodel1$values,type = "b")
lines(c(0, 50), c(1, 1))
plot(pcmodel2$values,type = "b")
lines(c(0, 50), c(1, 1))
r1 <- R1[R1$time == 1, c(grep(names(R1), pattern = "cesd"),
grep(names(R1), pattern = "RX"),
grep(names(R1), pattern = "GSE"))]
r2 <- R2[R2$time == 2, c(grep(names(R2), pattern = "cesd"),
grep(names(R2), pattern = "RX"),
grep(names(R2), pattern = "GSE"))]
colnames(r1) <- c("cesd1.1","cesd1.2","cesd1.3","cesd1.4","cesd1.5","cesd1.6","cesd1.7","cesd1.8","cesd1.9","cesd1.10","cesd1.11",
"cesd1.12","cesd1.13","cesd1.14","cesd1.15","cesd1.16","cesd1.17","cesd1.18","cesd1.19","cesd1.20","RX1","GSE1")
colnames(r2) <- c("cesd2.1","cesd2.2","cesd2.3","cesd2.4","cesd2.5","cesd2.6","cesd2.7","cesd2.8","cesd2.9","cesd2.10","cesd2.11",
"cesd2.12","cesd2.13","cesd2.14","cesd2.15","cesd2.16","cesd2.17","cesd2.18","cesd2.19","cesd2.20","RX1","GSE2")
cols1 <- c(grep(names(r1), pattern = "RX1"), grep(names(r1), pattern = "cesdtot"))
cols2 <- c(grep(names(r2), pattern = "RX1"), grep(names(r2), pattern = "cesdtot"))
r1_nct <- r1[as.character(R1$ID) %in% as.character(R2$ID), -cols1]
r2_nct <- r2[, -cols2]
r1_nct <- r1_nct[rowSums(is.na(r1_nct)) == 0,]
r2_nct <- r2_nct[rowSums(is.na(r2_nct)) == 0,]
data <- data.frame(r1_nct,r2_nct)
#View(data)（消除#以查看）
# 在生成r1和r2时保留ID列
r1 <- R1[R1$time == 1, c(
grep(names(R1), pattern = "cesd"),
grep(names(R1), pattern = "RX"),
grep(names(R1), pattern = "GSE"),
which(names(R1) == "ID")
)]
r2 <- R2[R2$time == 2, c(
grep(names(R2), pattern = "cesd"),
grep(names(R2), pattern = "RX"),
grep(names(R2), pattern = "GSE"),
which(names(R2) == "ID")
)]
# 调整列名（确保与列顺序对应）
colnames(r1) <- c(paste0("cesd1.", 1:20), "RX1", "GSE1", "ID")
colnames(r2) <- c(paste0("cesd2.", 1:20), "RX1", "GSE2", "ID")
# 定义需要排除的列（排除RX1和CESDtot）
cols1 <- c(grep("RX1", names(r1)), grep("CESDtot", names(r1)))
cols2 <- c(grep("RX1", names(r2)), grep("CESDtot", names(r2)))
# 筛选共有的ID
common_ids <- intersect(R1$ID, R2$ID)
r1_nct <- r1[r1$ID %in% common_ids, -cols1]
r2_nct <- r2[r2$ID %in% common_ids, -cols2]
# 处理缺失值（严格删除含NA的行）
r1_nct <- r1_nct[complete.cases(r1_nct), ]
r2_nct <- r2_nct[complete.cases(r2_nct), ]
# 再次筛选共有的ID（确保处理缺失值后一致）
final_common_ids <- intersect(r1_nct$ID, r2_nct$ID)
r1_nct <- r1_nct[r1_nct$ID %in% final_common_ids, ]
r2_nct <- r2_nct[r2_nct$ID %in% final_common_ids, ]
# 按ID排序后合并
r1_nct <- r1_nct[order(r1_nct$ID), ]
r2_nct <- r2_nct[order(r2_nct$ID), ]
# 检查行数是否一致
stopifnot(nrow(r1_nct) == nrow(r2_nct))
# 合并数据框
data <- data.frame(r1_nct, r2_nct)
#View(data)（消除#以查看）
k <- 21 #k is the number of variables at each time point
adjMat <- matrix(0, k, k) #set up empty matrix of coefficients
rsquarelist <- rep(0, k)
for (i in 1:k){
set.seed(100)
#i <- 2
lassoreg <- cv.glmnet(as.matrix(data[,1:k]), data[,(k+i)],
family = "gaussian", alpha = 1, standardize=FALSE)
lambda <- lassoreg$lambda.min
adjMat[1:k,i] <- coef(lassoreg, s = lambda, exact = FALSE)[2:(k+1)]
}
cor.data <- cor_auto(data)
groups <- c(rep("cesd",20),rep("GSE",1))
L <- qgraph(cor.data, graph = "glasso", layout = "spring", sampleSize = nrow(data), DoNotPlot = FALSE)$layout
########################第二种数据导入
setwd("D:/R/新数据")
R1 <- read.table(file = "数据一.txt", header = TRUE)
R2 <- read.table(file = "数据二.txt", header = TRUE)
#View(R1)（消除#以查看）
#View(R2)（消除#以查看）
R1_nct <- R1[, c(-2,-34)]##去除“time”列和“RX”列
#View(R1_nct)（消除#以查看）
R2_nct <- R2[, c(-2,-34)]
#View(R2_nct)（消除#以查看）
r <- merge(R1_nct,R2_nct,by="ID")#合并数据
#View(r)（消除#以查看）
is.na(r) == 0
r[r == NA] <- NA
nrow(r)
r <- r[complete.cases(r),]
nrow(r)
r <- r[, -1]
#View(r)（消除#以查看）
k <- 21 #k is the number of variables at each time point
adjMat <- matrix(0, k, k) #set up empty matrix of coefficients
rsquarelist <- rep(0, k)
# 生成邻接矩阵 adjMat 的循环部分
for (i in 1:k){
set.seed(100)
lassoreg <- cv.glmnet(
as.matrix(r[, 1:k]),
r[, (k + i)],
family = "gaussian",
alpha = 1,
standardize = TRUE
)
lambda <- lassoreg$lambda.min
adjMat[1:k, i] <- coef(lassoreg, s = lambda, exact = FALSE)[2:(k + 1)]
}
# 消除自回归路径：将对角线元素设为0
diag(adjMat) <- 0 # 添加这一行
# 修改后的绘图部分
groups <- c(rep("cesd", 20), rep("GSE", 1))
labels <- c("cesd1", "cesd2","cesd3","cesd4","cesd5","cesd6","cesd7","cesd8","cesd9","cesd10","cesd11",
"cesd12","cesd13","cesd14","cesd15","cesd16","cesd17","cesd18","cesd19", "cesd20", "GSE") # 确保标签完整
# 定义颜色参数（长度需与组数一致）
group_colors <- c("#E69F00", "#56B4E9") # 两个颜色对应两个组
# 绘图
layout(mat = matrix(c(1, 2), 1, 2), widths = c(2, 1))
qgraph(
adjMat,
groups = groups,
labels = labels,
legend = TRUE,
color = group_colors # 使用正确的参数名 color
)
plot.new()
adjMat2 <- adjMat
diag(adjMat2) <- 0 #set the AR paths to 0 for plotting purposes:
#######可预测性和影响
library(lavaan)
pred0 <- rep(0, k) #set up empty vectors to hold predictability coefficients
pred1 <- rep(0, k)
pred2 <- rep(0, k)
for (i in 1:k){ #loop through each variable
name.i <- names(r)[k+1] #what is the name of the DV
group.i <- groups[i] #what construct does the DV belong to
include.mod0 <- names(r)[1:k] #include self at previous time point
include.mod1 <- names(r)[1:k][-i] #exclude self at previous time point
include.mod2 <- names(r)[1:k][groups != group.i] #include predictors belonging to same construct
betas.mod0 <- adjMat[1:k,i] #fix regression coefficients to the ones obtained by glmnet
betas.mod1 <- adjMat[1:k,i][-i]
betas.mod2 <- adjMat[1:k,i][groups != group.i]
mod0 <- paste(name.i, " ~ ", paste(paste(betas.mod0,'*', include.mod0, sep = ""), collapse = "+"))
mod1 <- paste(name.i, " ~ ", paste(paste(betas.mod1,'*', include.mod1, sep = ""), collapse = "+"))
mod2 <- paste(name.i, " ~ ", paste(paste(betas.mod2,'*', include.mod2, sep = ""), collapse = "+"))
fit0 <- sem(mod0, sample.cov = cor(r), sample.nobs = nrow(r)) #fit model with subset of predictors to get r-square
pred0[i] <- inspect(fit0, "rsquare")
fit1 <- sem(mod1, sample.cov = cor(r), sample.nobs = nrow(r)) #fit model with subset of predictors to get r-square
pred1[i] <- inspect(fit1, "rsquare")
fit2 <- sem(mod2, sample.cov = cor(r), sample.nobs = nrow(r)) #fit model with subset of predictors to get r-square
pred2[i] <- inspect(fit2, "rsquare")
}
pred0
pred1
pred2
#nodewise influence: sum squared outgoing regression paths/node
infl0 <- apply(adjMat^2, 1, sum) #all outgoing nodes including AR
infl1 <- apply(adjMat2^2, 1, sum) #all outgoing nodes excluding AR
adjMat3 <- adjMat2 #adjMat3 sets paths within the same construct to 0
adjMat3[groups[row(adjMat3)] == groups[col(adjMat3)]] <- 0
infl2 <- apply(adjMat3^2, 1, sum) #rows = T1 predictors; columns = T2 DVs
infl0
infl1
infl2
# 定义 brieflegend
brieflegend <- c("cesd1", "cesd2","cesd3","cesd4","cesd5","cesd6","cesd7","cesd8","cesd9","cesd10","cesd11",
"cesd12","cesd13","cesd14","cesd15","cesd16","cesd17","cesd18","cesd19", "cesd20", "GSE") # 定义图示
# 取消注释以下两行启用图形设备
# jpeg("PredictabilityInfluencePlot.jpg", width=8, height=6, res=800, units="in")
# 调整布局和边距
layout(mat = matrix(c(1, 2, 3, 3), 2, 2), widths = c(1.5, 1))
par(mar = c(2, 6, 0, 0), las = 1, mgp = c(4, 1, 0)) # 关键修改点
# 绘制柱状图（手动控制 Y 轴标签位置）
barplot(rev(pred1), horiz = TRUE, names.arg = rev(labels),
col = rainbow(k), ylab = "")
title(ylab = "predictability", line = 4.5) # 关键修改点：手动定位标签
barplot(rev(infl1), horiz = TRUE, names.arg = rev(labels),
col = rainbow(k), ylab = "")
title(ylab = "influence", line = 4.5) # 关键修改点：手动定位标签
# 添加图例
plot.new()
legend(x = "center", inset = c(0, 0), bty = "n", fill=rev(rainbow(k)),
legend = brieflegend, cex = 1.2)
dev.off() #(uncomment this code to save plot as jpeg)
# 出预期影响: 节点作为预测变量时对其他所有变量的总影响
outgoing_IE <- rowSums(abs(adjMat2))
# 入预期影响: 节点作为结果变量时受到其他所有变量的总影响
incoming_IE <- colSums(abs(adjMat2))
# 创建数据框
influence_data <- data.frame(
node = labels,
outgoing = outgoing_IE,
incoming = incoming_IE
)
# 设置绘图布局
layout(matrix(c(1, 2), nrow = 1))
par(mar = c(5, 10, 4, 2)) # 调整边距：下、左、上、右
# 1. 出预期影响柱状图
barplot(rev(influence_data$outgoing),
horiz = TRUE,
names.arg = rev(influence_data$node),
las = 1, # 水平标签
col = "#4E84C4", # 蓝色系
border = NA,
main = "出预期影响 (Outgoing Expected Influence)",
xlab = "影响强度",
cex.names = 0.8,
xlim = c(0, max(outgoing_IE) * 1.1))
# 2. 入预期影响柱状图
barplot(rev(influence_data$incoming),
horiz = TRUE,
names.arg = rev(influence_data$node),
las = 1, # 水平标签
col = "#D16103", # 橙色系
border = NA,
main = "入预期影响 (Incoming Expected Influence)",
xlab = "影响强度",
cex.names = 0.8,
xlim = c(0, max(incoming_IE) * 1.1))
# 计算预期影响指标
outgoing_IE <- rowSums(abs(adjMat2)) # 出预期影响
incoming_IE <- colSums(abs(adjMat2)) # 入预期影响
# 创建数据框
influence_data <- data.frame(
node = labels,
outgoing = outgoing_IE,
incoming = incoming_IE
)
# 按出预期影响降序排序
influence_data <- influence_data[order(influence_data$outgoing, decreasing = TRUE), ]
# 输出到控制台
cat("\n=== 节点预期影响分析 ===\n")
print(influence_data, row.names = FALSE)
# 输出详细报告
cat("\n=== 详细报告 ===\n")
cat(sprintf("总节点数: %d\n", nrow(influence_data)))
cat(sprintf("平均出预期影响: %.3f\n", mean(influence_data$outgoing)))
cat(sprintf("平均入预期影响: %.3f\n", mean(influence_data$incoming)))
cat(sprintf("最大出预期影响: %.3f (节点: %s)\n",
max(influence_data$outgoing),
influence_data$node[which.max(influence_data$outgoing)]))
cat(sprintf("最大入预期影响: %.3f (节点: %s)\n",
max(influence_data$incoming),
influence_data$node[which.max(influence_data$incoming)]))
cat("\n节点影响度排名:\n")
for(i in 1:nrow(influence_data)) {
cat(sprintf("%d. %s: 出影响=%.3f, 入影响=%.3f\n",
i,
influence_data$node[i],
influence_data$outgoing[i],
influence_data$incoming[i]))
}
