Skip to content

Commit 6b87d09

Browse files
AI4HEALTH-LAB-THUAI4HEALTH-LAB-THU
authored
Update MainAnalysis.R
1 parent e24f428 commit 6b87d09

1 file changed

Lines changed: 95 additions & 43 deletions

File tree

‎Statistical_Analysis/MainAnalysis.R‎

Lines changed: 95 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -156,14 +156,14 @@ for(i in 1:length(disease)) {
156156
dat_cox <- subset(dat_age, time > 0)
157157
}
158158

159-
folds <- createFolds(dat_cox$event, k = 3)
159+
folds <- createFolds(dat_cox$event, k = 10)
160160
for (i in 1:length(var_ls)) {
161161
var <- var_ls[i]
162162
c_index_values <- c()
163163
c_index_lower_ls <- c()
164164
c_index_upper_ls <- c()
165165

166-
for(j in 1:3) {
166+
for(j in 1:10) {
167167
# train and test set
168168
test_indices <- folds[[j]]
169169
train_data <- dat_cox[-test_indices, ]
@@ -184,17 +184,31 @@ for(i in 1:length(disease)) {
184184
surv.time = test_data$time,
185185
surv.event = test_data$event)
186186
c_index <- concordance_result$c.index
187-
c_index_lower <- concordance_result$lower
188-
c_index_upper <- concordance_result$upper
187+
# c_index_lower <- concordance_result$lower
188+
# c_index_upper <- concordance_result$upper
189189
# save the c-index
190190
c_index_values <- c(c_index_values, c_index)
191-
c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
192-
c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
191+
# c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
192+
# c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
193193
print(paste0(item, " ------------ ", var, " ------------ fold ", j))
194194
}
195-
mean_c_index <- round(mean(c_index_values), digits = 3)
196-
mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
197-
mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
195+
# 计算均值和标准误（SE）
196+
mean_c_index <- mean(c_index_values)
197+
n_folds <- length(c_index_values)
198+
se_c_index <- sd(c_index_values) / sqrt(n_folds)
199+
200+
# 使用t分布计算置信区间（自由度为n_folds-1）
201+
t_value <- qt(0.975, df = n_folds - 1)
202+
mean_c_index_lower <- mean_c_index - t_value * se_c_index
203+
mean_c_index_upper <- mean_c_index + t_value * se_c_index
204+
205+
mean_c_index <- round(mean_c_index, digits = 3)
206+
mean_c_index_lower <- round(mean_c_index_lower, digits = 3)
207+
mean_c_index_upper <- round(mean_c_index_upper, digits = 3)
208+
209+
# mean_c_index <- round(mean(c_index_values), digits = 3)
210+
# mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
211+
# mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
198212

199213
var_mean_c_index <- c(var_mean_c_index, mean_c_index)
200214
var_mean_c_index_lower <- c(var_mean_c_index_lower, mean_c_index_lower)
@@ -624,20 +638,17 @@ for(i in 1:length(disease)) {
624638
# cause-specific cox
625639
dat_cox <- dat_merge_analysis
626640

627-
## data is small, use all data; results are similar
628-
# folds <- createFolds(dat_cox$event, k = 3)
641+
folds <- createFolds(dat_cox$event, k = 10)
629642
for (i in 1:length(var_ls)) {
630643
var <- var_ls[i]
631644
c_index_values <- c()
632645
c_index_lower_ls <- c()
633646
c_index_upper_ls <- c()
634647

635-
for(j in 1:1) {
636-
# test_indices <- folds[[j]]
637-
# train_data <- dat_cox[-test_indices, ]
638-
# test_data <- dat_cox[test_indices, ]
639-
train_data <- dat_cox
640-
test_data <- dat_cox
648+
for(j in 1:10) {
649+
test_indices <- folds[[j]]
650+
train_data <- dat_cox[-test_indices, ]
651+
test_data <- dat_cox[test_indices, ]
641652

642653
# build cause-specific cox model
643654
formula_covariates <- paste0("survobj ~ ", var)
@@ -654,18 +665,32 @@ for(i in 1:length(disease)) {
654665
surv.time = test_data$time,
655666
surv.event = test_data$event)
656667
c_index <- concordance_result$c.index
657-
c_index_lower <- concordance_result$lower
658-
c_index_upper <- concordance_result$upper
668+
# c_index_lower <- concordance_result$lower
669+
# c_index_upper <- concordance_result$upper
659670

660671
# save c-index
661672
c_index_values <- c(c_index_values, c_index)
662-
c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
663-
c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
673+
# c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
674+
# c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
664675
print(paste0(item, " ------------ ", var, " ------------ fold ", j))
665676
}
666-
mean_c_index <- round(mean(c_index_values), digits = 3)
667-
mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
668-
mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
677+
# 计算均值和标准误（SE）
678+
mean_c_index <- mean(c_index_values)
679+
n_folds <- length(c_index_values)
680+
se_c_index <- sd(c_index_values) / sqrt(n_folds)
681+
682+
# 使用t分布计算置信区间（自由度为n_folds-1）
683+
t_value <- qt(0.975, df = n_folds - 1)
684+
mean_c_index_lower <- mean_c_index - t_value * se_c_index
685+
mean_c_index_upper <- mean_c_index + t_value * se_c_index
686+
687+
mean_c_index <- round(mean_c_index, digits = 3)
688+
mean_c_index_lower <- round(mean_c_index_lower, digits = 3)
689+
mean_c_index_upper <- round(mean_c_index_upper, digits = 3)
690+
691+
# mean_c_index <- round(mean(c_index_values), digits = 3)
692+
# mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
693+
# mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
669694

670695
var_mean_c_index <- c(var_mean_c_index, mean_c_index)
671696
var_mean_c_index_lower <- c(var_mean_c_index_lower, mean_c_index_lower)
@@ -1233,14 +1258,14 @@ for(i in 1:length(disease)) {
12331258
dat_cox <- subset(dat_age, time > 0)
12341259
}
12351260

1236-
folds <- createFolds(dat_cox$event, k = 3)
1261+
folds <- createFolds(dat_cox$event, k = 10)
12371262
for (i in 1:length(var_ls)) {
12381263
var <- var_ls[i]
12391264
c_index_values <- c()
12401265
c_index_lower_ls <- c()
12411266
c_index_upper_ls <- c()
12421267

1243-
for(j in 1:3) {
1268+
for(j in 1:10) {
12441269
test_indices <- folds[[j]]
12451270
train_data <- dat_cox[-test_indices, ]
12461271
test_data <- dat_cox[test_indices, ]
@@ -1260,17 +1285,31 @@ for(i in 1:length(disease)) {
12601285
surv.time = test_data$time,
12611286
surv.event = test_data$event)
12621287
c_index <- concordance_result$c.index
1263-
c_index_lower <- concordance_result$lower
1264-
c_index_upper <- concordance_result$upper
1288+
# c_index_lower <- concordance_result$lower
1289+
# c_index_upper <- concordance_result$upper
12651290
# 存储c-index
12661291
c_index_values <- c(c_index_values, c_index)
1267-
c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
1268-
c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
1292+
# c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
1293+
# c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
12691294
print(paste0(item, " ------------ ", var, " ------------ fold ", j))
12701295
}
1271-
mean_c_index <- round(mean(c_index_values), digits = 3)
1272-
mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
1273-
mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
1296+
# 计算均值和标准误（SE）
1297+
mean_c_index <- mean(c_index_values)
1298+
n_folds <- length(c_index_values)
1299+
se_c_index <- sd(c_index_values) / sqrt(n_folds)
1300+
1301+
# 使用t分布计算置信区间（自由度为n_folds-1）
1302+
t_value <- qt(0.975, df = n_folds - 1)
1303+
mean_c_index_lower <- mean_c_index - t_value * se_c_index
1304+
mean_c_index_upper <- mean_c_index + t_value * se_c_index
1305+
1306+
mean_c_index <- round(mean_c_index, digits = 3)
1307+
mean_c_index_lower <- round(mean_c_index_lower, digits = 3)
1308+
mean_c_index_upper <- round(mean_c_index_upper, digits = 3)
1309+
1310+
# mean_c_index <- round(mean(c_index_values), digits = 3)
1311+
# mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
1312+
# mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
12741313

12751314
var_mean_c_index <- c(var_mean_c_index, mean_c_index)
12761315
var_mean_c_index_lower <- c(var_mean_c_index_lower, mean_c_index_lower)
@@ -2655,14 +2694,14 @@ for(i in 1:length(disease)) {
26552694
dat_cox <- subset(dat_age, time > 0)
26562695
}
26572696

2658-
folds <- createFolds(dat_cox$event, k = 3)
2697+
folds <- createFolds(dat_cox$event, k = 10)
26592698
for (i in 1:length(var_ls)) {
26602699
var <- var_ls[i]
26612700
c_index_values <- c()
26622701
c_index_lower_ls <- c()
26632702
c_index_upper_ls <- c()
26642703

2665-
for(j in 1:3) {
2704+
for(j in 1:10) {
26662705
# train and test split
26672706
test_indices <- folds[[j]]
26682707
train_data <- dat_cox[-test_indices, ]
@@ -2683,17 +2722,31 @@ for(i in 1:length(disease)) {
26832722
surv.time = test_data$time,
26842723
surv.event = test_data$event)
26852724
c_index <- concordance_result$c.index
2686-
c_index_lower <- concordance_result$lower
2687-
c_index_upper <- concordance_result$upper
2725+
# c_index_lower <- concordance_result$lower
2726+
# c_index_upper <- concordance_result$upper
26882727
# 存储c-index
26892728
c_index_values <- c(c_index_values, c_index)
2690-
c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
2691-
c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
2729+
# c_index_lower_ls <- c(c_index_lower_ls, c_index_lower)
2730+
# c_index_upper_ls <- c(c_index_upper_ls, c_index_upper)
26922731
print(paste0(item, " ------------ ", var, " ------------ fold ", j))
26932732
}
2694-
mean_c_index <- round(mean(c_index_values), digits = 3)
2695-
mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
2696-
mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
2733+
# 计算均值和标准误（SE）
2734+
mean_c_index <- mean(c_index_values)
2735+
n_folds <- length(c_index_values)
2736+
se_c_index <- sd(c_index_values) / sqrt(n_folds)
2737+
2738+
# 使用t分布计算置信区间（自由度为n_folds-1）
2739+
t_value <- qt(0.975, df = n_folds - 1)
2740+
mean_c_index_lower <- mean_c_index - t_value * se_c_index
2741+
mean_c_index_upper <- mean_c_index + t_value * se_c_index
2742+
2743+
mean_c_index <- round(mean_c_index, digits = 3)
2744+
mean_c_index_lower <- round(mean_c_index_lower, digits = 3)
2745+
mean_c_index_upper <- round(mean_c_index_upper, digits = 3)
2746+
2747+
# mean_c_index <- round(mean(c_index_values), digits = 3)
2748+
# mean_c_index_lower <- round(mean(c_index_lower_ls), digits = 3)
2749+
# mean_c_index_upper <- round(mean(c_index_upper_ls), digits = 3)
26972750

26982751
var_mean_c_index <- c(var_mean_c_index, mean_c_index)
26992752
var_mean_c_index_lower <- c(var_mean_c_index_lower, mean_c_index_lower)
@@ -4215,4 +4268,3 @@ p <- ggplot(res_telomere, aes(x = beta, y = var_name)) +
42154268
ggsave("fig6-f2.pdf", p, width = 12, height = 5)
42164269

42174270

4218-

0 commit comments

Comments
 (0)