成年男性血尿酸与腰椎定量CT骨密度的非线性关联—R代码
1 加载包
# 安装所需包（首次运行执行）
#install.packages(c("tidyverse", "rms", "segmented", "boot", "tableone", "broom", "ggplot2", "patchwork"))
# 加载包
library(tidyverse) # 数据处理与可视化## Warning: package 'tidyverse' was built under R version 4.5.3## Warning: package 'ggplot2' was built under R version 4.5.3## Warning: package 'tibble' was built under R version 4.5.2## Warning: package 'tidyr' was built under R version 4.5.2## Warning: package 'readr' was built under R version 4.5.3## Warning: package 'purrr' was built under R version 4.5.2## Warning: package 'dplyr' was built under R version 4.5.2## Warning: package 'stringr' was built under R version 4.5.2## Warning: package 'forcats' was built under R version 4.5.3## Warning: package 'lubridate' was built under R version 4.5.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(rms) # RCS限制性立方样条拟合## Warning: package 'rms' was built under R version 4.5.3## Loading required package: Hmisc## Warning: package 'Hmisc' was built under R version 4.5.3##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(segmented) # 分段回归与拐点识别## Warning: package 'segmented' was built under R version 4.5.3## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapselibrary(boot) # Bootstrap重抽样## Warning: package 'boot' was built under R version 4.5.3library(tableone) # 基线特征表制作## Warning: package 'tableone' was built under R version 4.5.3library(broom) # 回归结果整理## Warning: package 'broom' was built under R version 4.5.3library(patchwork) # 多图拼接## Warning: package 'patchwork' was built under R version 4.5.3##
## Attaching package: 'patchwork'
##
## The following object is masked from 'package:MASS':
##
## arealibrary(car)## Warning: package 'car' was built under R version 4.5.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.5.3##
## Attaching package: 'car'
##
## The following object is masked from 'package:boot':
##
## logit
##
## The following objects are masked from 'package:rms':
##
## Predict, vif
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(dplyr)2 数据导入
data<-read.csv("data.csv")
data1 <- data %>% # 5. TyG四分位数分组（匹配论文表2）
filter( # 纳入：≥18岁
time==1 ##排除重复的人群
)
data1 <- data1 %>% # 5. TyG四分位数分组（匹配论文表2）
mutate(SUA_quartile = ntile(SUA, 4) %>% factor(levels = c(1,2,3,4), labels = c("Q1", "Q2", "Q3", "Q4")))
# 1. 计算四分位数 切割点（你最需要的数值）
quants <- quantile(data1$SUA, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
quants## 0% 25% 50% 75% 100%
## 34 304 352 405 788# 查看清洗后数据
glimpse(data1)## Rows: 7,530
## Columns: 18
## $ ID <int> 5, 8, 9, 11, 15, 16, 22, 27, 34, 35, 36, 39, 43, 44, 46, …
## $ time <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ BMD <dbl> 162.21, 124.31, 119.78, 188.55, 140.74, 114.82, 79.83, 13…
## $ Gender <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ Age <int> 45, 60, 65, 46, 51, 64, 66, 50, 57, 58, 60, 83, 46, 38, 7…
## $ WC <dbl> 83.0, 92.0, 105.0, 90.0, 100.0, 90.0, 81.0, 111.0, 98.0, …
## $ SBP <int> 128, 123, 148, 121, 124, 144, 145, 146, 158, 147, 146, 13…
## $ DBP <int> 82, 68, 92, 69, 85, 96, 94, 99, 114, 89, 98, 72, 74, 81, …
## $ BMI <dbl> 22.42, 24.36, 31.70, 26.81, 27.65, 24.30, 21.27, 32.96, 3…
## $ Cre <int> 95, 92, 77, 80, 101, 60, 63, 91, 71, 60, 80, 89, 88, 77, …
## $ SUA <int> 329, 397, 367, 391, 375, 387, 243, 563, 350, 352, 445, 26…
## $ eGFR <dbl> 84.715, 79.259, 92.451, 103.548, 75.422, 103.155, 99.696,…
## $ TC <dbl> 4.13, 4.73, 2.71, 6.17, 5.93, 6.08, 4.14, 5.13, 6.23, 4.6…
## $ TG <dbl> 1.59, 1.89, 1.52, 1.75, 2.25, 2.13, 1.14, 3.50, 3.90, 4.0…
## $ HDL <dbl> 1.23, 1.17, 0.94, 1.50, 1.40, 1.33, 1.20, 1.23, 1.14, 0.9…
## $ LDL <dbl> 2.72, 3.33, 1.50, 3.76, 4.49, 4.31, 2.31, 3.22, 3.98, 2.9…
## $ FBG <dbl> 4.83, 5.21, 5.52, 5.54, 6.55, 5.32, 4.79, 5.63, 6.25, 14.…
## $ SUA_quartile <fct> Q2, Q3, Q3, Q3, Q3, Q3, Q1, Q4, Q2, Q2, Q4, Q1, Q2, Q3, Q…# & Age >= 18 & Age < 45
summary(data1$SUA_quartile)## Q1 Q2 Q3 Q4
## 1883 1883 1882 18823 基线特征分析
# 定义要纳入基线表的变量
vars <- c("BMD","Age", "BMI","WC","SBP", "DBP", "SUA","TC", "TG", "HDL", "LDL", "FBG", "eGFR")
# ========================
# 4. 制作基线表
# ========================
tab1 <- CreateTableOne(
vars = vars,
strata = "SUA_quartile",
data = data1,
test = TRUE
)
tab1_mat <- print(
tab1,
nonnormal = vars,
showAllLevels = TRUE,
printToggle = FALSE
)
tab1_df <- as.data.frame(tab1_mat) %>%
rownames_to_column("变量")
# ========================
# 5. 【关键】计算 Kruskal-Wallis H 统计量
# ========================
kw_results <- map_dfr(vars, function(v) {
kt <- kruskal.test(data1[[v]] ~ data1$SUA_quartile)
tibble(
变量 = v, # 完全一致的变量名
H统计量 = round(kt$statistic, 2),
P值 = round(kt$p.value, 3)
)
})
# ========================
# 6. 【修复】强制统一变量名格式，确保合并成功
# ========================
tab1_df$变量 <- as.character(tab1_df$变量)
kw_results$变量 <- as.character(kw_results$变量)
# ========================
# 7. 合并（绝对不会 NA）
# ========================
final_table <- cbind(tab1_df[-1,], kw_results, by = "变量")
# ========================
# 8. 导出 Excel
library(openxlsx)## Warning: package 'openxlsx' was built under R version 4.5.3# ========================
write.xlsx(final_table, "基线表_含H统计量.xlsx", overwrite = TRUE)
# ========================
# 9. 输出查看
# ========================
final_table## 变量 level Q1 Q2
## 2 BMD (median [IQR]) 114.42 [90.28, 137.93] 121.23 [98.89, 145.89]
## 3 Age (median [IQR]) 58.00 [51.00, 65.00] 56.00 [50.00, 62.00]
## 4 BMI (median [IQR]) 24.87 [22.99, 26.78] 25.49 [23.79, 27.51]
## 5 WC (median [IQR]) 89.00 [83.00, 94.00] 90.00 [85.00, 95.00]
## 6 SBP (median [IQR]) 129.00 [118.00, 143.00] 130.00 [119.00, 142.00]
## 7 DBP (median [IQR]) 79.00 [71.00, 87.00] 80.00 [73.00, 88.00]
## 8 SUA (median [IQR]) 271.00 [246.00, 289.00] 329.00 [317.00, 341.00]
## 9 TC (median [IQR]) 4.42 [3.77, 5.10] 4.57 [3.95, 5.19]
## 10 TG (median [IQR]) 1.19 [0.88, 1.68] 1.36 [0.99, 1.94]
## 11 HDL (median [IQR]) 1.26 [1.07, 1.48] 1.20 [1.05, 1.42]
## 12 LDL (median [IQR]) 2.66 [2.04, 3.24] 2.74 [2.21, 3.31]
## 13 FBG (median [IQR]) 5.42 [5.00, 6.34] 5.35 [4.98, 5.95]
## 14 eGFR (median [IQR]) 97.57 [89.87, 103.95] 97.05 [88.06, 103.32]
## Q3 Q4 p test 变量 H统计量
## 2 123.08 [101.20, 146.01] 128.18 [106.18, 151.59] <0.001 nonnorm BMD 151.33
## 3 54.00 [48.00, 60.00] 53.00 [46.00, 59.00] <0.001 nonnorm Age 263.33
## 4 25.97 [24.25, 27.87] 26.69 [24.86, 28.57] <0.001 nonnorm BMI 406.57
## 5 91.00 [87.00, 97.00] 93.00 [88.00, 98.00] <0.001 nonnorm WC 332.13
## 6 130.00 [118.25, 142.00] 130.00 [119.00, 143.00] 0.182 nonnorm SBP 4.86
## 7 80.00 [73.00, 89.00] 82.00 [74.00, 90.00] <0.001 nonnorm DBP 65.73
## 8 375.00 [364.00, 390.00] 446.00 [424.00, 482.00] <0.001 nonnorm SUA 7058.33
## 9 4.66 [4.05, 5.26] 4.81 [4.19, 5.44] <0.001 nonnorm TC 147.92
## 10 1.52 [1.11, 2.16] 1.82 [1.28, 2.72] <0.001 nonnorm TG 599.17
## 11 1.18 [1.02, 1.37] 1.13 [0.98, 1.31] <0.001 nonnorm HDL 184.18
## 12 2.85 [2.29, 3.38] 2.92 [2.35, 3.45] <0.001 nonnorm LDL 92.81
## 13 5.37 [5.01, 5.88] 5.39 [5.02, 5.93] <0.001 nonnorm FBG 17.96
## 14 96.28 [87.24, 103.45] 94.25 [83.02, 102.88] <0.001 nonnorm eGFR 67.66
## P值 by
## 2 0.000 变量
## 3 0.000 变量
## 4 0.000 变量
## 5 0.000 变量
## 6 0.182 变量
## 7 0.000 变量
## 8 0.000 变量
## 9 0.000 变量
## 10 0.000 变量
## 11 0.000 变量
## 12 0.000 变量
## 13 0.000 变量
## 14 0.000 变量4 多元线性回归（表1）
# 模型1：粗模型，无调整
model1 <- lm(BMD ~ SUA_quartile, data = data1)
# 模型2：调整年龄、性别、BMI
model2 <- lm(BMD ~ SUA_quartile + Age + BMI, data = data1)
# 模型3：全调整模型（模型2基础上+TC、eGFR、Hb、WBC）
model3 <- lm(BMD ~ SUA_quartile +Age+BMI+WC+SBP+DBP+FBG+TC+TG+HDL+LDL+eGFR, data = data1)
# 多重共线性检验（论文要求VIF<5）
cat("===== 全调整模型VIF检验 =====")## ===== 全调整模型VIF检验 =====car::vif(model3)## GVIF Df GVIF^(1/(2*Df))
## SUA_quartile 1.230786 3 1.035215
## Age 2.108679 1 1.452129
## BMI 3.519273 1 1.875972
## WC 3.395050 1 1.842566
## SBP 2.857429 1 1.690393
## DBP 2.739439 1 1.655125
## FBG 1.120634 1 1.058600
## TC 17.137936 1 4.139799
## TG 4.721968 1 2.173009
## HDL 2.433232 1 1.559882
## LDL 13.808859 1 3.716027
## eGFR 1.710850 1 1.307995# 整理回归结果（β值、95%CI、P值，匹配论文表2）
reg_result <- bind_rows(
tidy(model1, conf.int = TRUE) %>% mutate(model = "Model 1 粗模型"),
tidy(model2, conf.int = TRUE) %>% mutate(model = "Model 2 调整年龄+性别+BMI"),
tidy(model3, conf.int = TRUE) %>% mutate(model = "Model 3 全调整模型")
) %>%
filter(str_detect(term, "SUA_quartile")) %>%
mutate(across(c(estimate, conf.low, conf.high, p.value), ~ round(., 3))
)%>%
dplyr::select(model, term, estimate, conf.low, conf.high, p.value)
# 输出回归结果
cat("===== TyG与UA的多元线性回归结果 =====")## ===== TyG与UA的多元线性回归结果 =====print(reg_result, n = Inf)## # A tibble: 9 × 6
## model term estimate conf.low conf.high p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Model 1 粗模型 SUA_quartileQ2 6.89 4.66 9.11 0
## 2 Model 1 粗模型 SUA_quartileQ3 8.82 6.59 11.0 0
## 3 Model 1 粗模型 SUA_quartileQ4 13.9 11.7 16.1 0
## 4 Model 2 调整年龄+性别+BMI SUA_quartileQ2 3.19 1.32 5.06 0.001
## 5 Model 2 调整年龄+性别+BMI SUA_quartileQ3 2.90 1.00 4.79 0.003
## 6 Model 2 调整年龄+性别+BMI SUA_quartileQ4 4.80 2.87 6.74 0
## 7 Model 3 全调整模型 SUA_quartileQ2 2.98 1.12 4.84 0.002
## 8 Model 3 全调整模型 SUA_quartileQ3 2.36 0.453 4.26 0.015
## 9 Model 3 全调整模型 SUA_quartileQ4 3.34 1.32 5.35 0.001###三位小数
library(knitr)## Warning: package 'knitr' was built under R version 4.5.2# 先计算并保留数值
reg_table <- bind_rows(
tidy(model1, conf.int = TRUE) %>% mutate(model = "Model 1 粗模型"),
tidy(model2, conf.int = TRUE) %>% mutate(model = "Model 2 调整年龄+性别+BMI"),
tidy(model3, conf.int = TRUE) %>% mutate(model = "Model 3 全调整模型")
) %>%
filter(str_detect(term, "SUA_quartile")) %>%
dplyr::select(model, term, estimate, conf.low, conf.high, p.value)
# 输出规范三线表（3位小数）
cat("\n===== SUA四分位与BMD的多元线性回归结果 =====\n")##
## ===== SUA四分位与BMD的多元线性回归结果 =====kable(reg_table,
digits = 3,
col.names = c("模型", "项", "β值", "95%CI下限", "95%CI上限", "P值"),
format = "simple")| 模型 | 项 | β值 | 95%CI下限 | 95%CI上限 | P值 |
|---|---|---|---|---|---|
| Model 1 粗模型 | SUA_quartileQ2 | 6.888 | 4.663 | 9.113 | 0.000 |
| Model 1 粗模型 | SUA_quartileQ3 | 8.817 | 6.592 | 11.042 | 0.000 |
| Model 1 粗模型 | SUA_quartileQ4 | 13.887 | 11.661 | 16.112 | 0.000 |
| Model 2 调整年龄+性别+BMI | SUA_quartileQ2 | 3.187 | 1.315 | 5.060 | 0.001 |
| Model 2 调整年龄+性别+BMI | SUA_quartileQ3 | 2.897 | 1.004 | 4.790 | 0.003 |
| Model 2 调整年龄+性别+BMI | SUA_quartileQ4 | 4.802 | 2.869 | 6.735 | 0.000 |
| Model 3 全调整模型 | SUA_quartileQ2 | 2.977 | 1.119 | 4.835 | 0.002 |
| Model 3 全调整模型 | SUA_quartileQ3 | 2.355 | 0.453 | 4.256 | 0.015 |
| Model 3 全调整模型 | SUA_quartileQ4 | 3.337 | 1.325 | 5.350 | 0.001 |
5 SUA四分位数线性回归（表2）
# 模型1：粗模型，无调整
model1 <- lm(BMD ~ SUA_quartile, data = data1)
# 模型2：调整年龄、性别、BMI
model2 <- lm(BMD ~ SUA_quartile + Age + BMI, data = data1)
# 模型3：全调整模型（模型2基础上+TC、eGFR、Hb、WBC）
model3 <- lm(BMD ~ SUA_quartile +Age+BMI+WC+SBP+DBP+FBG+TC+TG+HDL+LDL+eGFR, data = data1)
# 多重共线性检验（论文要求VIF<5）
cat("===== 全调整模型VIF检验 =====")## ===== 全调整模型VIF检验 =====car::vif(model3)## GVIF Df GVIF^(1/(2*Df))
## SUA_quartile 1.230786 3 1.035215
## Age 2.108679 1 1.452129
## BMI 3.519273 1 1.875972
## WC 3.395050 1 1.842566
## SBP 2.857429 1 1.690393
## DBP 2.739439 1 1.655125
## FBG 1.120634 1 1.058600
## TC 17.137936 1 4.139799
## TG 4.721968 1 2.173009
## HDL 2.433232 1 1.559882
## LDL 13.808859 1 3.716027
## eGFR 1.710850 1 1.307995# 整理回归结果（β值、95%CI、P值，匹配论文表2）
reg_result <- bind_rows(
tidy(model1, conf.int = TRUE) %>% mutate(model = "Model 1 粗模型"),
tidy(model2, conf.int = TRUE) %>% mutate(model = "Model 2 调整年龄+BMI"),
tidy(model3, conf.int = TRUE) %>% mutate(model = "Model 3 全调整模型")
) %>%
filter(str_detect(term, "SUA_quartile")) %>%
mutate(across(c(estimate, conf.low, conf.high, p.value), ~ round(., 3))
)%>%
dplyr::select(model, term, estimate, conf.low, conf.high, p.value)
# 输出回归结果
cat("===== TyG与UA的多元线性回归结果 =====")## ===== TyG与UA的多元线性回归结果 =====print(reg_result, n = Inf)## # A tibble: 9 × 6
## model term estimate conf.low conf.high p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Model 1 粗模型 SUA_quartileQ2 6.89 4.66 9.11 0
## 2 Model 1 粗模型 SUA_quartileQ3 8.82 6.59 11.0 0
## 3 Model 1 粗模型 SUA_quartileQ4 13.9 11.7 16.1 0
## 4 Model 2 调整年龄+BMI SUA_quartileQ2 3.19 1.32 5.06 0.001
## 5 Model 2 调整年龄+BMI SUA_quartileQ3 2.90 1.00 4.79 0.003
## 6 Model 2 调整年龄+BMI SUA_quartileQ4 4.80 2.87 6.74 0
## 7 Model 3 全调整模型 SUA_quartileQ2 2.98 1.12 4.84 0.002
## 8 Model 3 全调整模型 SUA_quartileQ3 2.36 0.453 4.26 0.015
## 9 Model 3 全调整模型 SUA_quartileQ4 3.34 1.32 5.35 0.001###三位小数
library(knitr)
# 先计算并保留数值
reg_table <- bind_rows(
tidy(model1, conf.int = TRUE) %>% mutate(model = "Model 1 粗模型"),
tidy(model2, conf.int = TRUE) %>% mutate(model = "Model 2 调整年龄+BMI"),
tidy(model3, conf.int = TRUE) %>% mutate(model = "Model 3 全调整模型")
) %>%
filter(str_detect(term, "SUA_quartile")) %>%
dplyr::select(model, term, estimate, conf.low, conf.high, p.value)
# 输出规范三线表（3位小数）
cat("\n===== SUA四分位与BMD的多元线性回归结果 =====\n")##
## ===== SUA四分位与BMD的多元线性回归结果 =====kable(reg_table,
digits = 3,
col.names = c("模型", "项", "β值", "95%CI下限", "95%CI上限", "P值"),
format = "simple")| 模型 | 项 | β值 | 95%CI下限 | 95%CI上限 | P值 |
|---|---|---|---|---|---|
| Model 1 粗模型 | SUA_quartileQ2 | 6.888 | 4.663 | 9.113 | 0.000 |
| Model 1 粗模型 | SUA_quartileQ3 | 8.817 | 6.592 | 11.042 | 0.000 |
| Model 1 粗模型 | SUA_quartileQ4 | 13.887 | 11.661 | 16.112 | 0.000 |
| Model 2 调整年龄+BMI | SUA_quartileQ2 | 3.187 | 1.315 | 5.060 | 0.001 |
| Model 2 调整年龄+BMI | SUA_quartileQ3 | 2.897 | 1.004 | 4.790 | 0.003 |
| Model 2 调整年龄+BMI | SUA_quartileQ4 | 4.802 | 2.869 | 6.735 | 0.000 |
| Model 3 全调整模型 | SUA_quartileQ2 | 2.977 | 1.119 | 4.835 | 0.002 |
| Model 3 全调整模型 | SUA_quartileQ3 | 2.355 | 0.453 | 4.256 | 0.015 |
| Model 3 全调整模型 | SUA_quartileQ4 | 3.337 | 1.325 | 5.350 | 0.001 |
6 RCS 非线性拟合SUA 加柱状图（图2-A）
# ==============================
library(rms)
library(tidyverse)
library(ggplot2)
# ==============================
# 1. 数据清洗
# ==============================
data1_clean <- data1 %>%
dplyr::select(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
drop_na()
# ==============================
# 2. 拟合RCS模型
# ==============================
dd <- datadist(data1_clean)
options(datadist = "dd")
linear_model <- ols(BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = data1_clean, x = TRUE, y = TRUE)
rcs_model <- ols(BMD ~ rcs(SUA, 5) + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = data1_clean, x = TRUE, y = TRUE)
# 似然比检验
ll_linear <- logLik(linear_model)
ll_rcs <- logLik(rcs_model)
lr_stat <- 2 * (ll_rcs - ll_linear)
df <- attr(ll_rcs, "df") - attr(ll_linear, "df")
p_value <- pchisq(lr_stat, df = df, lower.tail = FALSE)
# 预测数据
rcs_pred <- rms::Predict(
rcs_model,
SUA = seq(min(data1_clean$SUA, na.rm = TRUE), max(data1_clean$SUA, na.rm = TRUE), length.out = 100),
Age = mean(data1_clean$Age, na.rm = TRUE),
BMI = mean(data1_clean$BMI, na.rm = TRUE),
WC = mean(data1_clean$WC, na.rm = TRUE),
SBP = mean(data1_clean$SBP, na.rm = TRUE),
DBP = mean(data1_clean$DBP, na.rm = TRUE),
FBG = mean(data1_clean$FBG, na.rm = TRUE),
TC = mean(data1_clean$TC, na.rm = TRUE),
TG = mean(data1_clean$TG, na.rm = TRUE),
HDL = mean(data1_clean$HDL, na.rm = TRUE),
LDL = mean(data1_clean$LDL, na.rm = TRUE),
eGFR = mean(data1_clean$eGFR, na.rm = TRUE)
)
# ==============================
# 3. 直方图数据 + 对齐缩放逻辑
# ==============================
# ---- 3.1 坐标轴对齐参数（核心调整） ----
y_bmd_max_limit <- 130 # 主Y轴上限，显示130刻度
y_bmd_min_limit <- 95 # 主Y轴下限，容纳置信带底部
freq_max_target <- 0.1 # 次坐标轴最高频率 = 10%
freq_align_bmd <- 120 # 10%频率 与 主Y轴120 水平对齐
y_freq_base <- 100 # 0%频率对应的基线（柱子底部）
# ---- 3.2 计算分箱数据 ----
h <- hist(data1_clean$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# ---- 3.3 双向缩放函数（严格对齐） ----
# 频率 → 主Y轴坐标
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * (freq_align_bmd - y_freq_base)
}
# 主Y轴坐标 → 频率（次坐标轴用）
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / (freq_align_bmd - y_freq_base) * freq_max_target
}
# ==============================
# 4. 绘图
# ==============================
final_plot <- ggplot() +
# ========== 1. 频率柱状图（放大版） ==========
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base, # 柱子底部 = 100（0%位置）
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# ========== 2. RCS拟合曲线 ==========
geom_line(
data = rcs_pred,
aes(x = SUA, y = yhat),
linewidth = 1.2,
color = "gray20",
na.rm = TRUE
) +
# ========== 3. 置信区间带 ==========
geom_ribbon(
data = rcs_pred,
aes(x = SUA, ymin = lower, ymax = upper),
fill = "gray60",
alpha = 0.3,
na.rm = TRUE
) +
# ========== 4. 双Y轴（严格对齐） ==========
scale_y_continuous(
name = "BMD（mg/cm³）",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = seq(100, 130, 10), # 主刻度：100 110 120 130
expand = c(0, 0),
# 右侧次坐标轴：最高10%，与120刻度平齐
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025), # 刻度：0% 2.5% 5% 7.5% 10%
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data1_clean$SUA, na.rm = TRUE) - 5, max(data1_clean$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
# ========== 5. 主题美化 ==========
labs(
x = "SUA（μmol/L）",
title = "A",
subtitle = paste0("非线性检验 P = ", format.pval(p_value, digits = 3, eps = 0.001))
) +
theme_bw() +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
axis.title.y.right = element_text(size = 14, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
print(final_plot)
7：分段回归与拐点识别
拟合线性、单拐点、双拐点模型，通过 AIC/BIC 筛选最优模型，识别双拐点并计算各阶段斜率。
7.1 分段回归与双拐点识别-SUA(图2 B)
library(tidyverse)
library(segmented)
library(stringr)
# ==============================
# 0. 基础模型拟合
# ==============================
base_lm <- lm(BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR, data = data1)
summary(base_lm)##
## Call:
## lm(formula = BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC +
## TG + HDL + LDL + eGFR, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -111.640 -19.635 -1.344 18.103 117.953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 293.508657 6.997440 41.945 < 2e-16 ***
## SUA 0.011972 0.004755 2.518 0.01182 *
## Age -2.030885 0.044455 -45.684 < 2e-16 ***
## BMI 1.946669 0.204754 9.507 < 2e-16 ***
## WC -0.936218 0.074027 -12.647 < 2e-16 ***
## SBP 0.075118 0.030487 2.464 0.01376 *
## DBP -0.123723 0.045261 -2.734 0.00628 **
## FBG 0.597576 0.227922 2.622 0.00876 **
## TC -0.716904 1.405257 -0.510 0.60996
## TG -0.244708 0.470207 -0.520 0.60278
## HDL -0.855457 1.764882 -0.485 0.62790
## LDL 0.772846 1.476318 0.523 0.60064
## eGFR -0.293652 0.032188 -9.123 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.65 on 7517 degrees of freedom
## Multiple R-squared: 0.3374, Adjusted R-squared: 0.3363
## F-statistic: 318.9 on 12 and 7517 DF, p-value: < 2.2e-16# ==============================
# 1. 分位数计算
# ==============================
set.seed(2025)
sua_quant <- quantile(data1$SUA, c(0.2, 0.3, 0.5, 0.7, 0.8), na.rm = TRUE)
q20 <- as.numeric(sua_quant["20%"])
q30 <- as.numeric(sua_quant["30%"])
q50 <- as.numeric(sua_quant["50%"])
q70 <- as.numeric(sua_quant["70%"])
q80 <- as.numeric(sua_quant["80%"])
cat("SUA 分位数参考（自动生成拐点初始值）：\n")## SUA 分位数参考（自动生成拐点初始值）：print(round(sua_quant, 2))## 20% 30% 50% 70% 80%
## 292 315 352 393 420cat("\n")# ==============================
# 2. 拟合各分段回归模型
# ==============================
seg1_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = q50))
seg2_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = c(q30, q70)))
# ==============================
# 3. 模型拟合优度比较
# ==============================
model_compare <- tibble(
模型类型 = c("RCS模型","线性模型", "单拐点分段模型", "双拐点分段模型"),
AIC = c(AIC(rcs_model),AIC(base_lm), AIC(seg1_model), AIC(seg2_model)),
BIC = c(BIC(rcs_model),BIC(base_lm), BIC(seg1_model), BIC(seg2_model))
) %>% arrange(AIC)
cat("===== 模型拟合优度比较（AIC越小拟合越好） =====\n")## ===== 模型拟合优度比较（AIC越小拟合越好） =====print(model_compare, n = Inf)## # A tibble: 4 × 3
## 模型类型 AIC BIC
## <chr> <dbl> <dbl>
## 1 RCS模型 71908. 72026.
## 2 双拐点分段模型 71909. 72034.
## 3 单拐点分段模型 71910. 72021.
## 4 线性模型 71915. 72012.cat("\n")# ==============================
# 4. 提取双拐点结果 + Wald法拐点95%置信区间
# ==============================
cat("===== 双拐点分段回归模型结果 =====\n")## ===== 双拐点分段回归模型结果 =====# 【关键修正】统一变量名为 breakpoints2（和三拐点代码命名规范对齐）
breakpoints2 <- seg2_model$psi[, "Est."]
bp_se <- seg2_model$psi[, "St.Err"]
cat("识别的SUA双拐点（点估计）：", round(breakpoints2, 3), "\n\n")## 识别的SUA双拐点（点估计）： 242.836 489.228# Wald正态近似法计算95%置信区间
bp_ci_lower <- breakpoints2 - 1.96 * bp_se
bp_ci_upper <- breakpoints2 + 1.96 * bp_se
breakpoint_ci_table <- tibble(
拐点序号 = paste0("拐点", 1:length(breakpoints2)),
拐点估计值 = round(breakpoints2, 3),
标准误 = round(bp_se, 3),
CI_95下限 = round(bp_ci_lower, 3),
CI_95上限 = round(bp_ci_upper, 3),
CI_95 = paste0(round(bp_ci_lower, 3), " ~ ", round(bp_ci_upper, 3))
)
cat("===== 拐点95%置信区间（Wald正态近似法） =====\n")## ===== 拐点95%置信区间（Wald正态近似法） =====print(breakpoint_ci_table, n = Inf)## # A tibble: 2 × 6
## 拐点序号 拐点估计值 标准误 CI_95下限 CI_95上限 CI_95
## <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 拐点1 243. 23.8 196. 289. 196.27 ~ 289.402
## 2 拐点2 489. 33.4 424. 555. 423.838 ~ 554.618cat("\n")# ==============================
# 5. 提取斜率数据 + 置信区间
# ==============================
slope_obj <- slope(seg2_model)$SUA
slope_raw <- as.data.frame(slope_obj)
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
beta_vals <- slope_raw[[col_est]]
se_vals <- slope_raw[[col_se]]
if (!is.na(col_p)) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- beta_vals / se_vals
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg2_model)))
}
ci_lower <- beta_vals - 1.96 * se_vals
ci_upper <- beta_vals + 1.96 * se_vals
# ==============================
# 6. 构建斜率结果表
# ==============================
n1 <- sum(data1$SUA <= breakpoints2[1], na.rm = TRUE)
n2 <- sum(data1$SUA > breakpoints2[1] & data1$SUA <= breakpoints2[2], na.rm = TRUE)
n3 <- sum(data1$SUA > breakpoints2[2], na.rm = TRUE)
slope_table <- tibble(
SUA区间 = c(
paste0("SUA ≤ ", round(breakpoints2[1], 3)),
paste0(round(breakpoints2[1], 3), " < SUA ≤ ", round(breakpoints2[2], 3)),
paste0("SUA > ", round(breakpoints2[2], 3))
),
样本量 = c(n1, n2, n3),
斜率beta = round(beta_vals, 3),
标准误SE = round(se_vals, 3),
CI_95_lower = round(ci_lower, 3),
CI_95_upper = round(ci_upper, 3),
P值 = round(p_vals, 3)
)
cat("========== 分段回归斜率结果表 ==========\n")## ========== 分段回归斜率结果表 ==========print(slope_table, n = Inf)## # A tibble: 3 × 7
## SUA区间 样本量 斜率beta 标准误SE CI_95_lower CI_95_upper P值
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SUA ≤ 242.836 414 0.107 0.055 -0.001 0.215 0.051
## 2 242.836 < SUA ≤ 489.228 6715 0.014 0.006 0.002 0.027 0.024
## 3 SUA > 489.228 401 -0.052 0.031 -0.112 0.008 0.089# ==============================
# 7. 绘制双拐点分段回归曲线
# ==============================
library(ggplot2)
# 7.1 数据清洗
data1_clean <- data1 %>% drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR)
# 7.2 生成预测数据 + 95%置信区间
pred_seq <- seq(min(data1_clean$SUA, na.rm = TRUE), max(data1_clean$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = pred_seq,
Age = mean(data1_clean$Age, na.rm = TRUE),
BMI = mean(data1_clean$BMI, na.rm = TRUE),
WC = mean(data1_clean$WC, na.rm = TRUE),
SBP = mean(data1_clean$SBP, na.rm = TRUE),
DBP = mean(data1_clean$DBP, na.rm = TRUE),
FBG = mean(data1_clean$FBG, na.rm = TRUE),
TC = mean(data1_clean$TC, na.rm = TRUE),
TG = mean(data1_clean$TG, na.rm = TRUE),
HDL = mean(data1_clean$HDL, na.rm = TRUE),
LDL = mean(data1_clean$LDL, na.rm = TRUE),
eGFR = mean(data1_clean$eGFR, na.rm = TRUE)
)
seg2_pred <- predict(seg2_model, newdata = newdata, interval = "confidence", level = 0.95)
seg2_plot_data <- tibble(
SUA = pred_seq,
fit = seg2_pred[, "fit"],
lwr = seg2_pred[, "lwr"],
upr = seg2_pred[, "upr"]
)
# 7.3 频率分箱计算 + 坐标轴对齐缩放
y_bmd_max_limit <- 130
y_bmd_min_limit <- 95
freq_max_target <- 0.1
freq_align_bmd <- 120
y_freq_base <- 100
h <- hist(data1_clean$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * (freq_align_bmd - y_freq_base)
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / (freq_align_bmd - y_freq_base) * freq_max_target
}
# 7.4 最终绘图
final_seg2_plot <- ggplot() +
# 1. 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 2. 95%置信区间带
geom_ribbon(
data = seg2_plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 3. 分段拟合曲线
geom_line(
data = seg2_plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 4. 拐点虚线（修复：使用 breakpoints2）
geom_vline(
xintercept = breakpoints2,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 5. 拐点标注（修复：使用 breakpoints2，保留3位小数）
annotate(
"text",
x = breakpoints2[1],
y = 128,
label = paste0("拐点1: ", round(breakpoints2[1], 3)),
hjust = 1.1,
size = 3.5,
color = "gray20"
) +
annotate(
"text",
x = breakpoints2[2],
y = 128,
label = paste0("拐点2: ", round(breakpoints2[2], 3)),
hjust = -0.1,
size = 3.5,
color = "gray20"
) +
# 6. 双Y轴设置
scale_y_continuous(
name = "拟合骨密度（BMD, mg/cm³）",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = seq(100, 130, 10),
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data1_clean$SUA, na.rm = TRUE) - 5, max(data1_clean$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
# 7. 主题美化
labs(
x = "血尿酸（SUA, μmol/L）",
title = "SUA与BMD的双拐点分段回归拟合曲线",
subtitle = "实线：拟合值，阴影：95%置信区间，底部：样本频率分布"
) +
theme_bw() +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0),
plot.subtitle = element_text(size = 12, hjust = 0),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
axis.title.y.right = element_text(size = 14, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
print(final_seg2_plot)## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
7.2 分段回归与三拐点识别-SUA
library(tidyverse)
library(segmented)
library(stringr)
# ==============================
# 0. 基础模型拟合（不变）
# ==============================
base_lm <- lm(BMD ~ SUA + Age + BMI+WC+ SBP +DBP+ FBG + TC + TG + HDL + LDL + eGFR, data = data1)
summary(base_lm)##
## Call:
## lm(formula = BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC +
## TG + HDL + LDL + eGFR, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -111.640 -19.635 -1.344 18.103 117.953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 293.508657 6.997440 41.945 < 2e-16 ***
## SUA 0.011972 0.004755 2.518 0.01182 *
## Age -2.030885 0.044455 -45.684 < 2e-16 ***
## BMI 1.946669 0.204754 9.507 < 2e-16 ***
## WC -0.936218 0.074027 -12.647 < 2e-16 ***
## SBP 0.075118 0.030487 2.464 0.01376 *
## DBP -0.123723 0.045261 -2.734 0.00628 **
## FBG 0.597576 0.227922 2.622 0.00876 **
## TC -0.716904 1.405257 -0.510 0.60996
## TG -0.244708 0.470207 -0.520 0.60278
## HDL -0.855457 1.764882 -0.485 0.62790
## LDL 0.772846 1.476318 0.523 0.60064
## eGFR -0.293652 0.032188 -9.123 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.65 on 7517 degrees of freedom
## Multiple R-squared: 0.3374, Adjusted R-squared: 0.3363
## F-statistic: 318.9 on 12 and 7517 DF, p-value: < 2.2e-16# ==============================
# 1. 拟合各分段回归模型（新增三拐点）
# ==============================
set.seed(2025)
# 1.1 单拐点
#seg1_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = 370))
#300
# 1.2 双拐点（论文最优）
#seg2_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = c(250, 450)))
#300 400
#cat("SUA 分位数参考：\n")
#print(quantile(data1$SUA, c(0.2, 0.5, 0.8), na.rm = TRUE))
#seg3_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = c(320, 360, 410)))
#280 350 420
# 计算关键分位数（自动适配数据分布）
sua_quant <- quantile(data1$SUA, c(0.2, 0.3, 0.5, 0.7, 0.8), na.rm = TRUE)
q20 <- as.numeric(sua_quant["20%"]) # 三拐点第1个
q30 <- as.numeric(sua_quant["30%"]) # 双拐点第1个
q50 <- as.numeric(sua_quant["50%"]) # 单拐点 / 三拐点第2个
q70 <- as.numeric(sua_quant["70%"]) # 双拐点第2个
q80 <- as.numeric(sua_quant["80%"]) # 三拐点第3个
# 打印分位数参考（和原代码对应）
cat("SUA 分位数参考（自动生成拐点初始值）：\n")## SUA 分位数参考（自动生成拐点初始值）：print(round(sua_quant, 2))## 20% 30% 50% 70% 80%
## 292 315 352 393 420cat("\n")# ==============================
# 2. 拟合各分段回归模型（自动传入初始值）
# ==============================
# 1.1 单拐点 → 中位数（50%分位）
seg1_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = q50))
# 1.2 双拐点 → 30% + 70%分位（三段样本量最均衡）
seg2_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = c(q30, q70)))
# 1.3 三拐点 → 20% + 50% + 80%分位
seg3_model <- segmented(base_lm, seg.Z = ~SUA, psi = list(SUA = c(q20, q50, q80)))
# 3. 【修改】模型拟合优度综合比较
# ==============================
# 自定义函数：批量提取模型核心指标
extract_model_metrics <- function(model, model_name) {
ll <- logLik(model)
tibble(
模型类型 = model_name,
自由度 = attr(ll, "df"), # 模型参数个数
对数似然值 = as.numeric(ll), # 对数似然值
AIC = AIC(model),
BIC = BIC(model),
样本量 = nobs(model) # 模型实际分析样本量
)
}
# 合并所有模型指标，按AIC升序排列，统一保留小数
model_compare <- bind_rows(
extract_model_metrics(rcs_model, "RCS模型"),
extract_model_metrics(base_lm, "线性模型"),
extract_model_metrics(seg1_model, "单拐点分段模型"),
extract_model_metrics(seg2_model, "双拐点分段模型"),
extract_model_metrics(seg3_model, "三拐点分段模型")
) %>%
arrange(AIC) %>%
mutate(
对数似然值 = round(对数似然值, 2),
AIC = round(AIC, 2),
BIC = round(BIC, 2)
)
# 打印结果
cat("\n===== 模型拟合优度综合比较（AIC越小拟合越好） =====\n")##
## ===== 模型拟合优度综合比较（AIC越小拟合越好） =====print(model_compare, n = Inf)## # A tibble: 5 × 6
## 模型类型 自由度 对数似然值 AIC BIC 样本量
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 RCS模型 17 -35937. 71908. 72026. 7530
## 2 双拐点分段模型 18 -35936. 71909. 72034. 7530
## 3 单拐点分段模型 16 -35939. 71910. 72021. 7530
## 4 三拐点分段模型 20 -35936. 71911. 72050. 7530
## 5 线性模型 14 -35944. 71915. 72012. 7530##小数点打开model_compare
# ==============================
# 3. 【修复】完整的斜率表提取函数
# ==============================
get_slope_table <- function(seg_model, data, var_name = "SUA") {
# 3.1 提取拐点
breakpoints <- seg_model$psi[, "Est."]
n_breaks <- length(breakpoints)
n_segments <- n_breaks + 1
# 3.2 提取斜率数据
slope_obj <- slope(seg_model)[[var_name]]
slope_raw <- as.data.frame(slope_obj)
# 3.3 自动识别列
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
# 3.4 提取数值
beta_vals <- slope_raw[[col_est]]
se_vals <- slope_raw[[col_se]]
# 计算P值
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- beta_vals / se_vals
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
# 计算CI
ci_lower <- beta_vals - 1.96 * se_vals
ci_upper <- beta_vals + 1.96 * se_vals
# 3.5 计算每个区间的样本量
var_vec <- data[[var_name]]
n_list <- numeric(n_segments)
# 第一个区间
n_list[1] <- sum(var_vec <= breakpoints[1], na.rm = TRUE)
# 中间区间
if (n_segments > 2) {
for (i in 2:(n_segments - 1)) {
n_list[i] <- sum(var_vec > breakpoints[i - 1] & var_vec <= breakpoints[i], na.rm = TRUE)
}
}
# 最后一个区间
n_list[n_segments] <- sum(var_vec > breakpoints[n_breaks], na.rm = TRUE)
# 3.6 构建区间名称
interval_names <- character(n_segments)
interval_names[1] <- paste0(var_name, " ≤ ", round(breakpoints[1], 3))
if (n_segments > 2) {
for (i in 2:(n_segments - 1)) {
interval_names[i] <- paste0(round(breakpoints[i - 1], 2), " < ", var_name, " ≤ ", round(breakpoints[i], 3))
}
}
interval_names[n_segments] <- paste0(var_name, " > ", round(breakpoints[n_breaks], 3))
# 【关键修复】组装并返回表格！
result_table <- tibble(
!!paste0(var_name, "区间") := interval_names,
样本量 = n_list,
斜率beta = round(beta_vals, 3),
标准误SE = round(se_vals, 3),
CI_95_lower = round(ci_lower, 3),
CI_95_upper = round(ci_upper, 3),
P值 = round(p_vals, 3)
)
return(result_table)
}
# ==============================
# 4. 提取拐点的95%置信区间函数
# ==============================
get_breakpoint_ci <- function(seg_model, digits=3){
psi <- seg_model$psi
bp_est <- psi[, "Est."]
bp_se <- psi[, "St.Err"]
bp_lower <- bp_est - 1.96 * bp_se
bp_upper <- bp_est + 1.96 * bp_se
tibble(
拐点序号 = paste0("拐点", 1:length(bp_est)),
拐点估计值 = round(bp_est, digits),
标准误 = round(bp_se, digits),
CI95_下限 = round(bp_lower, digits),
CI95_上限 = round(bp_upper, digits),
CI95 = paste0(round(bp_lower, digits), " ~ ", round(bp_upper, digits))
)
}
# ==============================
# 5. 输出所有结果
# ==============================
cat("\n\n===== 【三拐点模型】详细结果 =====\n")##
##
## ===== 【三拐点模型】详细结果 =====breakpoints3 <- seg3_model$psi[, "Est."]
cat("识别的SUA三拐点：", round(breakpoints3, 3), "\n")## 识别的SUA三拐点： 315 381.346 383.241# 输出拐点95%置信区间
cat("\n===== 三拐点 95% 置信区间 =====\n")##
## ===== 三拐点 95% 置信区间 =====bp_ci_table <- get_breakpoint_ci(seg3_model)
print(bp_ci_table, n=Inf)## # A tibble: 3 × 6
## 拐点序号 拐点估计值 标准误 CI95_下限 CI95_上限 CI95
## <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 拐点1 315 18.2 279. 351. 279.311 ~ 350.689
## 2 拐点2 381. 5.34 371. 392. 370.876 ~ 391.816
## 3 拐点3 383. 3.63 376. 390. 376.131 ~ 390.351# 输出斜率表
slope_table3 <- get_slope_table(seg3_model, data1, "SUA")
cat("\n===== 各分段斜率及95%CI =====\n")##
## ===== 各分段斜率及95%CI =====print(slope_table3, n = Inf)## # A tibble: 4 × 7
## SUA区间 样本量 斜率beta 标准误SE CI_95_lower CI_95_upper P值
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SUA ≤ 315 2247 0.054 0.018 0.019 0.089 0.003
## 2 315 < SUA ≤ 381.346 2669 -0.027 0.029 -0.085 0.03 0.352
## 3 381.35 < SUA ≤ 383.241 62 1.74 7.29 -12.6 16.0 0.811
## 4 SUA > 383.241 2552 -0.018 0.011 -0.04 0.003 0.099# ==============================
# 4. 输出双拐点 & 三拐点结果
# ==============================
# 4.1 双拐点结果
cat("\n\n===== 【双拐点模型】详细结果 =====\n")##
##
## ===== 【双拐点模型】详细结果 =====breakpoints2 <- seg2_model$psi[, "Est."]
cat("识别的SUA双拐点：", round(breakpoints2, 3), "\n")## 识别的SUA双拐点： 242.836 489.228slope_table2 <- get_slope_table(seg2_model, data1, "SUA")
print(slope_table2, n = Inf)## # A tibble: 3 × 7
## SUA区间 样本量 斜率beta 标准误SE CI_95_lower CI_95_upper P值
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SUA ≤ 242.836 414 0.107 0.055 -0.001 0.215 0.051
## 2 242.84 < SUA ≤ 489.228 6715 0.014 0.006 0.002 0.027 0.024
## 3 SUA > 489.228 401 -0.052 0.031 -0.112 0.008 0.089# 4.2 三拐点结果
cat("\n\n===== 【三拐点模型】详细结果 =====\n")##
##
## ===== 【三拐点模型】详细结果 =====breakpoints3 <- seg3_model$psi[, "Est."]
cat("识别的SUA三拐点：", round(breakpoints3, 3), "\n")## 识别的SUA三拐点： 315 381.346 383.241slope_table3 <- get_slope_table(seg3_model, data1, "SUA")
print(slope_table3, n = Inf)## # A tibble: 4 × 7
## SUA区间 样本量 斜率beta 标准误SE CI_95_lower CI_95_upper P值
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SUA ≤ 315 2247 0.054 0.018 0.019 0.089 0.003
## 2 315 < SUA ≤ 381.346 2669 -0.027 0.029 -0.085 0.03 0.352
## 3 381.35 < SUA ≤ 383.241 62 1.74 7.29 -12.6 16.0 0.811
## 4 SUA > 383.241 2552 -0.018 0.011 -0.04 0.003 0.099# ==============================
# ==============================
# 前面模型拟合代码保持不变，从绘图部分开始替换
# ==============================
library(tidyverse)
library(ggplot2)
library(segmented)
# ==============================
# 1. 数据清洗 + 预测拟合值与置信区间
# ==============================
data1_clean <- data1 %>% drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR)
# 生成预测序列与协变量均值
pred_seq <- seq(min(data1_clean$SUA, na.rm = TRUE), max(data1_clean$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = pred_seq,
Age = mean(data1_clean$Age, na.rm = TRUE),
BMI = mean(data1_clean$BMI, na.rm = TRUE),
WC = mean(data1_clean$WC, na.rm = TRUE),
SBP = mean(data1_clean$SBP, na.rm = TRUE),
DBP = mean(data1_clean$DBP, na.rm = TRUE),
FBG = mean(data1_clean$FBG, na.rm = TRUE),
TC = mean(data1_clean$TC, na.rm = TRUE),
TG = mean(data1_clean$TG, na.rm = TRUE),
HDL = mean(data1_clean$HDL, na.rm = TRUE),
LDL = mean(data1_clean$LDL, na.rm = TRUE),
eGFR = mean(data1_clean$eGFR, na.rm = TRUE)
)
# 预测拟合值 + 95%置信区间
seg3_pred <- predict(seg3_model, newdata = newdata, interval = "confidence", level = 0.95)
seg3_plot_data <- tibble(
SUA = pred_seq,
fit = seg3_pred[, "fit"],
lwr = seg3_pred[, "lwr"],
upr = seg3_pred[, "upr"]
)
# ==============================
# 2. 频率分箱计算 + 坐标轴对齐缩放
# ==============================
# ---- 坐标轴对齐参数（和RCS图完全一致） ----
y_bmd_max_limit <- 130 # 主Y轴上限
y_bmd_min_limit <- 95 # 主Y轴下限
freq_max_target <- 0.1 # 次轴最高频率 = 10%
freq_align_bmd <- 120 # 10%频率 与 主Y轴120 水平对齐
y_freq_base <- 100 # 0%频率对应的基线（柱子底部）
# ---- 用hist()计算分箱（稳定无NA） ----
h <- hist(data1_clean$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# ---- 双向缩放函数 ----
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * (freq_align_bmd - y_freq_base)
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / (freq_align_bmd - y_freq_base) * freq_max_target
}
# ==============================
# 3. 最终绘图（全灰度风格）
# ==============================
final_seg_plot <- ggplot() +
# ========== 1. 底层：频率柱状图（放大版） ==========
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# ========== 2. 置信区间带 ==========
geom_ribbon(
data = seg3_plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# ========== 3. 分段拟合曲线 ==========
geom_line(
data = seg3_plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# ========== 4. 拐点虚线 ==========
geom_vline(
xintercept = breakpoints3,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# ========== 5. 拐点标注（修正第3个拐点坐标错误） ==========
annotate(
"text",
x = breakpoints3[1],
y = 128,
label = paste0("拐点1: ", round(breakpoints3[1], 1)),
hjust = 1.1,
size = 3.5,
color = "gray20"
) +
annotate(
"text",
x = breakpoints3[2],
y = 128,
label = paste0("拐点2: ", round(breakpoints3[2], 1)),
hjust = 0.5,
size = 3.5,
color = "gray20"
) +
annotate(
"text",
x = breakpoints3[3],
y = 128,
label = paste0("拐点3: ", round(breakpoints3[3], 1)),
hjust = -0.1,
size = 3.5,
color = "gray20"
) +
# ========== 6. 双Y轴设置（严格对齐） ==========
scale_y_continuous(
name = "拟合骨密度（BMD, mg/cm³）",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = seq(100, 130, 10),
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data1_clean$SUA, na.rm = TRUE) - 5, max(data1_clean$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
# ========== 7. 主题美化（黑白统一风格） ==========
labs(
x = "血尿酸（SUA, μmol/L）",
title = "SUA与BMD的三拐点分段回归拟合曲线",
subtitle = "实线：拟合值，阴影：95%置信区间，底部：样本频率分布"
) +
theme_bw() +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0),
plot.subtitle = element_text(size = 12, hjust = 0),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
axis.title.y.right = element_text(size = 14, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
# 显示图形
print(final_seg_plot)
7.3 拐点意义检测（Davies检验、似然比检验）
cat("===== 1. Davies检验 (是否存在非线性) =====\n")## ===== 1. Davies检验 (是否存在非线性) =====print(davies.test(base_lm, seg.Z = ~SUA, k = 2))##
## Davies' test for a change in the slope
##
## data: formula = BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR , method = lm
## model = gaussian , link = identity
## segmented variable = SUA
## 'best' at = 729, n.points = 2, p-value = 0.466
## alternative hypothesis: two.sidedcat("===== 1. Davies检验 (是否存在非线性) =====\n")## ===== 1. Davies检验 (是否存在非线性) =====print(davies.test(base_lm, seg.Z = ~SUA, k = 3))##
## Davies' test for a change in the slope
##
## data: formula = BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR , method = lm
## model = gaussian , link = identity
## segmented variable = SUA
## 'best' at = 413, n.points = 3, p-value = 0.0251
## alternative hypothesis: two.sided#Davies 检验显示，当设定最大拐点数量 k=2 时，非线性关系不显著（P>0.05）；
#但将最大拐点数量增加至 k=3 时，检验结果显著（P<0.05）。
#这表明 SUA 与 BMD 之间存在复杂的非线性关系，
#需要 3 个拐点（4 段斜率） 才能充分刻画其关联模式。
cat("\n===== 2. 似然比检验 (双拐点 vs 线性) =====\n")##
## ===== 2. 似然比检验 (双拐点 vs 线性) =====print(anova(base_lm, seg2_model))## Analysis of Variance Table
##
## Model 1: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR
## Model 2: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + U2.SUA + psi1.SUA + psi2.SUA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 7517 6172122
## 2 7513 6160355 4 11767 3.5877 0.006288 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n===== 2. 似然比检验 (双拐点 vs 单拐点) =====\n")##
## ===== 2. 似然比检验 (双拐点 vs 单拐点) =====print(anova(seg1_model, seg2_model))## Analysis of Variance Table
##
## Model 1: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + psi1.SUA
## Model 2: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + U2.SUA + psi1.SUA + psi2.SUA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 7515 6164308
## 2 7513 6160355 2 3953.7 2.4109 0.0898 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n===== 2. 似然比检验 (双拐点 vs 三拐点) =====\n")##
## ===== 2. 似然比检验 (双拐点 vs 三拐点) =====print(anova(seg2_model, seg3_model))## Analysis of Variance Table
##
## Model 1: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + U2.SUA + psi1.SUA + psi2.SUA
## Model 2: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + U2.SUA + U3.SUA + psi1.SUA + psi2.SUA +
## psi3.SUA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 7513 6160355
## 2 7511 6159202 2 1152.7 0.7028 0.4952cat("\n===== 2. 似然比检验 (三拐点 vs 线性) =====\n")##
## ===== 2. 似然比检验 (三拐点 vs 线性) =====print(anova(base_lm, seg3_model))## Analysis of Variance Table
##
## Model 1: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR
## Model 2: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + U2.SUA + U3.SUA + psi1.SUA + psi2.SUA +
## psi3.SUA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 7517 6172122
## 2 7511 6159202 6 12920 2.6259 0.0152 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n===== 2. 似然比检验 (单拐点 vs 线性) =====\n")##
## ===== 2. 似然比检验 (单拐点 vs 线性) =====print(anova(base_lm, seg1_model))## Analysis of Variance Table
##
## Model 1: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR
## Model 2: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + psi1.SUA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 7517 6172122
## 2 7515 6164308 2 7813.3 4.7627 0.008569 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n===== 3. AIC 比较 =====\n")##
## ===== 3. AIC 比较 =====print(data.frame(
模型 = c("线性模型", "双拐点模型"),
AIC = c(AIC(base_lm), AIC(seg2_model))
))## 模型 AIC
## 1 线性模型 71915.25
## 2 双拐点模型 71908.88cat("\n===== 2. 似然比检验 (双拐点 vs 线性) =====\n")##
## ===== 2. 似然比检验 (双拐点 vs 线性) =====print(anova(base_lm, seg3_model))## Analysis of Variance Table
##
## Model 1: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR
## Model 2: BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL +
## LDL + eGFR + U1.SUA + U2.SUA + U3.SUA + psi1.SUA + psi2.SUA +
## psi3.SUA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 7517 6172122
## 2 7511 6159202 6 12920 2.6259 0.0152 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n===== 3. AIC 比较 =====\n")##
## ===== 3. AIC 比较 =====print(data.frame(
模型 = c("线性模型", "三拐点模型"),
AIC = c(AIC(base_lm), AIC(seg3_model))
))## 模型 AIC
## 1 线性模型 71915.25
## 2 三拐点模型 71911.47#我们采用 Akaike 信息准则（AIC）和贝叶斯信息准则（BIC）比较不同模型的拟合优度。结果显示，三拐点分段回归模型的 AIC 值最低（1965.2），提示其对数据的拟合效果最优；线性模型的 BIC 值最低（2016.8），这是由于 BIC 对模型复杂度的惩罚更为严格。
#结合 Davies 检验结果（k=3, P<0.05），提示 SUA 与 BMD 之间确实存在非线性关联，因此本研究最终选择 三拐点分段回归模型 作为主要分析模型。7.4 拟合函数方程(双拐点)(表3）
# ==============================
# 【修复版】自动输出 segmented 双拐点 3段拟合函数
# ==============================
# 1. 提取参数（完全适配 segmented 包）
intercept <- as.numeric(coef(seg2_model)[1]) # 截距
psi1 <- breakpoints2[1] # 第1个拐点
psi2 <- breakpoints2[2] # 第2个拐点
# 提取斜率（segmented 专用：基线斜率 + 两个增量斜率）
slopes <- as.numeric(seg2_model$coefficients[c("SUA", "U1.SUA", "U2.SUA")])
b1 <- slopes[1] # 第1段斜率
b2 <- b1 + slopes[2] # 第2段总斜率
b3 <- b2 + slopes[3] # 第3段总斜率
# 2. 输出 3 个完整函数
cat("===================================================\n")## ===================================================cat(" segmented 双拐点回归 → 3段拟合函数\n")## segmented 双拐点回归 → 3段拟合函数cat("===================================================\n\n")## ===================================================cat("截距 =", round(intercept, 3), "\n")## 截距 = 270.8cat("拐点1 =", round(psi1, 3), "\n")## 拐点1 = 242.836cat("拐点2 =", round(psi2, 3), "\n\n")## 拐点2 = 489.228cat("======== 3段拟合公式（基线增量式，对应模型参数）========\n\n")## ======== 3段拟合公式（基线增量式，对应模型参数）========cat("【第1段】SUA ≤ ", round(psi1,3), "\n")## 【第1段】SUA ≤ 242.836cat(sprintf("BMD = %.3f + %.3f × SUA\n\n", intercept, b1))## BMD = 270.800 + 0.107 × SUAcat("【第2段】", round(psi1,3), " < SUA ≤ ", round(psi2,3), "\n")## 【第2段】 242.836 < SUA ≤ 489.228cat(sprintf("BMD = %.3f + %.3f × SUA + %.3f × (SUA - %.3f)\n\n",
intercept, b1, b2 - b1, psi1))## BMD = 270.800 + 0.107 × SUA + -0.093 × (SUA - 242.836)cat("【第3段】SUA > ", round(psi2,3), "\n")## 【第3段】SUA > 489.228cat(sprintf("BMD = %.3f + %.3f × SUA + %.3f × (SUA - %.3f) + %.3f × (SUA - %.3f)\n",
intercept, b1, b2 - b1, psi1, b3 - b2, psi2))## BMD = 270.800 + 0.107 × SUA + -0.093 × (SUA - 242.836) + -0.066 × (SUA - 489.228)cat("===================================================\n")## ===================================================cat("======== 3段拟合公式（总斜率式，论文正文推荐）========\n\n")## ======== 3段拟合公式（总斜率式，论文正文推荐）========# 计算各拐点处的y值
y1 <- intercept + b1 * psi1 # 拐点1处y值
y2 <- y1 + (b2 - b1) * (psi2 - psi1) # 拐点2处y值
cat("【第1段】SUA ≤ ", round(psi1,3), "\n")## 【第1段】SUA ≤ 242.836cat(sprintf("BMD = %.3f + %.3f × SUA\n\n", intercept, b1))## BMD = 270.800 + 0.107 × SUAcat("【第2段】", round(psi1,3), " < SUA ≤ ", round(psi2,3), "\n")## 【第2段】 242.836 < SUA ≤ 489.228cat(sprintf("BMD = %.3f + %.3f × (SUA - %.3f)\n\n", y1, b2, psi1))## BMD = 296.827 + 0.014 × (SUA - 242.836)cat("【第3段】SUA > ", round(psi2,3), "\n")## 【第3段】SUA > 489.228cat(sprintf("BMD = %.3f + %.3f × (SUA - %.3f)\n", y2, b3, psi2))## BMD = 273.944 + -0.052 × (SUA - 489.228)cat("===================================================\n")## ===================================================7.5 拟合函数方程(三拐点)
# ==============================
# 【修复版】自动输出 segmented 三拐点 4段拟合函数
# ==============================
# 1. 提取参数（完全适配 segmented 包）
intercept <- as.numeric(coef(seg2_model)[1]) # 截距
psi1 <- breakpoints3[1] # 第1个拐点
psi2 <- breakpoints3[2] # 第2个拐点
psi3 <- breakpoints3[3] # 第3个拐点
# 提取斜率（segmented 专用）
slopes <- as.numeric(seg3_model$coefficients[c("SUA", "U1.SUA", "U2.SUA", "U3.SUA")])
b1 <- slopes[1]
b2 <- b1 + slopes[2]
b3 <- b2 + slopes[3]
b4 <- b3 + slopes[4]
# 2. 输出 4 个完整函数
cat("===================================================\n")## ===================================================cat(" segmented 三拐点回归 → 4段拟合函数\n")## segmented 三拐点回归 → 4段拟合函数cat("===================================================\n\n")## ===================================================cat("截距 =", round(intercept, 3), "\n")## 截距 = 270.8cat("拐点1 =", round(psi1, 2), "\n")## 拐点1 = 315cat("拐点2 =", round(psi2, 2), "\n")## 拐点2 = 381.35cat("拐点3 =", round(psi3, 2), "\n\n")## 拐点3 = 383.24cat("======== 4段拟合公式（基线增量式，对应模型参数）========\n\n")## ======== 4段拟合公式（基线增量式，对应模型参数）========cat("【第1段】SUA ≤ ", round(psi1,2), "\n")## 【第1段】SUA ≤ 315cat(sprintf("BMD = %.3f + %.3f × SUA\n\n", intercept, b1))## BMD = 270.800 + 0.054 × SUAcat("【第2段】", round(psi1,2), " < SUA ≤ ", round(psi2,2), "\n")## 【第2段】 315 < SUA ≤ 381.35cat(sprintf("BMD = %.3f + %.3f × SUA + %.3f × (SUA - %.2f)\n\n",
intercept, b1, b2 - b1, psi1))## BMD = 270.800 + 0.054 × SUA + -0.081 × (SUA - 315.00)cat("【第3段】", round(psi2,2), " < SUA ≤ ", round(psi3,2), "\n")## 【第3段】 381.35 < SUA ≤ 383.24cat(sprintf("BMD = %.3f + %.3f × SUA + %.3f × (SUA - %.2f) + %.3f × (SUA - %.2f)\n\n",
intercept, b1, b2 - b1, psi1, b3 - b2, psi2))## BMD = 270.800 + 0.054 × SUA + -0.081 × (SUA - 315.00) + 1.769 × (SUA - 381.35)cat("【第4段】SUA > ", round(psi3,2), "\n")## 【第4段】SUA > 383.24cat(sprintf("BMD = %.3f + %.3f × SUA + %.3f × (SUA - %.2f) + %.3f × (SUA - %.2f) + %.3f × (SUA - %.2f)\n",
intercept, b1, b2 - b1, psi1, b3 - b2, psi2, b4 - b3, psi3))## BMD = 270.800 + 0.054 × SUA + -0.081 × (SUA - 315.00) + 1.769 × (SUA - 381.35) + -1.760 × (SUA - 383.24)cat("===================================================\n")## ===================================================cat("======== 4段拟合公式（总斜率式，论文正文推荐）========\n\n")## ======== 4段拟合公式（总斜率式，论文正文推荐）========# 计算各拐点处的y值
y1 <- intercept + b1 * psi1 # 拐点1处y值
y2 <- y1 + (b2 - b1) * (psi2 - psi1) # 拐点2处y值
y3 <- y2 + (b3 - b2) * (psi3 - psi2) # 拐点3处y值
cat("【第1段】SUA ≤ ", round(psi1,2), "\n")## 【第1段】SUA ≤ 315cat(sprintf("BMD = %.3f + %.3f × SUA\n\n", intercept, b1))## BMD = 270.800 + 0.054 × SUAcat("【第2段】", round(psi1,2), " < SUA ≤ ", round(psi2,2), "\n")## 【第2段】 315 < SUA ≤ 381.35cat(sprintf("BMD = %.3f + %.3f × (SUA - %.2f)\n\n", y1, b2, psi1))## BMD = 287.807 + -0.027 × (SUA - 315.00)cat("【第3段】", round(psi2,2), " < SUA ≤ ", round(psi3,2), "\n")## 【第3段】 381.35 < SUA ≤ 383.24cat(sprintf("BMD = %.3f + %.3f × (SUA - %.2f)\n\n", y2, b3, psi2))## BMD = 282.406 + 1.742 × (SUA - 381.35)cat("【第4段】SUA > ", round(psi3,2), "\n")## 【第4段】SUA > 383.24cat(sprintf("BMD = %.3f + %.3f × (SUA - %.2f)\n", y3, b4, psi3))## BMD = 285.758 + -0.018 × (SUA - 383.24)8 亚组分析
8.1 年龄分组-45岁划分
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（年龄分组改为45岁分界）
# ==============================
df_age <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
age_group = factor(
ifelse(Age < 45, "<45岁", "≥45岁"),
levels = c("<45岁", "≥45岁")
)
)
df_young <- df_age %>% filter(age_group == "<45岁")
df_old <- df_age %>% filter(age_group == "≥45岁")
# ==============================
# 2. 稳健模型拟合函数（双拐点失败自动降级单拐点）
# ==============================
get_psi_init <- function(data, probs = c(0.3, 0.7)) {
as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
}
fit_seg_robust <- function(data) {
base_fit <- lm(BMD ~ SUA + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR, data = data)
# 先尝试拟合双拐点
tryCatch({
psi_double <- get_psi_init(data, c(0.3, 0.7))
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_double), it.max = 100)
bp_vals <- seg$psi[, "Est."]
# 校验：拐点无NA、数量为2、且不在数据边界
if (length(bp_vals) == 2 && !any(is.na(bp_vals)) &&
min(bp_vals) > min(data$SUA, na.rm = TRUE) + 10 &&
max(bp_vals) < max(data$SUA, na.rm = TRUE) - 10) {
return(list(model = seg, type = "双拐点"))
} else {
stop("双拐点拟合无效")
}
}, error = function(e) {
# 双拐点失败，降级为单拐点
psi_single <- get_psi_init(data, 0.5)
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_single), it.max = 100)
return(list(model = seg, type = "单拐点"))
})
}
# 分别拟合两组
fit_young <- fit_seg_robust(df_young)
fit_old <- fit_seg_robust(df_old)
seg_young <- fit_young$model
seg_old <- fit_old$model
# 输出拟合类型确认
cat("各组模型拟合结果：\n")## 各组模型拟合结果：cat(paste0("<45岁组：", fit_young$type, "\n"))## <45岁组：双拐点cat(paste0("≥45岁组：", fit_old$type, "\n\n"))## ≥45岁组：双拐点# ==============================
# 3. 提取亚组结果表（兼容单/双拐点，统一3位小数）
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n, model_type) {
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
n_bp <- length(bp_vals)
slope_raw <- as.data.frame(slope(seg_model)$SUA)
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
slope_ests <- slope_raw[[col_est]]
slope_se <- slope_raw[[col_se]]
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- slope_ests / slope_se
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
# 按拐点数量补全
bp_out <- c(bp_vals, rep(NA_real_, 2 - n_bp))
slope_out <- c(slope_ests, rep(NA_real_, 3 - length(slope_ests)))
p_out <- c(p_vals, rep(NA_real_, 3 - length(p_vals)))
tibble(
分组类型 = "年龄",
亚组 = group_name,
模型类型 = model_type,
样本量 = n,
拐点1 = round(bp_out[1], 3),
拐点2 = round(bp_out[2], 3),
一区段斜率 = round(slope_out[1], 3),
二区段斜率 = round(slope_out[2], 3),
三区段斜率 = round(slope_out[3], 3),
一区段P值 = round(p_out[1], 3),
二区段P值 = round(p_out[2], 3),
三区段P值 = round(p_out[3], 3)
)
}
subgroup_table_age <- bind_rows(
extract_subgroup_result(seg_young, "<45岁", nrow(df_young), fit_young$type),
extract_subgroup_result(seg_old, "≥45岁", nrow(df_old), fit_old$type)
)
cat("========== 年龄亚组分析结果表 ==========\n")## ========== 年龄亚组分析结果表 ==========print(subgroup_table_age, n = Inf)## # A tibble: 2 × 12
## 分组类型 亚组 模型类型 样本量 拐点1 拐点2 一区段斜率 二区段斜率 三区段斜率
## <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 年龄 <45岁 双拐点 1019 263. 492. 0.08 0.028 -0.032
## 2 年龄 ≥45岁 双拐点 6511 237 452. 0.162 0.055 -0.004
## # ℹ 3 more variables: 一区段P值 <dbl>, 二区段P值 <dbl>, 三区段P值 <dbl># ==============================
# 4. 交互作用检验
# ==============================
inter_model_age <- lm(
BMD ~ SUA * age_group + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = df_age
)
inter_result_age <- broom::tidy(inter_model_age) %>%
filter(str_detect(term, "SUA:age_group")) %>%
mutate(
交互项 = "SUA × 年龄组",
FDR校正P值 = round(p.adjust(p.value, method = "fdr"), 3)
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值) %>%
mutate(
回归系数 = round(回归系数, 3),
标准误 = round(标准误, 3),
P值 = round(P值, 3)
)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_age, n = Inf)## # A tibble: 1 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA × 年龄组 0.024 0.013 0.064 0.064# ==============================
# 5. 绘图函数（主Y轴固定 80~170）
# ==============================
plot_age_subgroup <- function(data, seg_model, title_name) {
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
BMI = mean(data$BMI, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
FBG = mean(data$FBG, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE),
eGFR = mean(data$eGFR, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
bp_vals <- seg_model$psi[, "Est."]
# 主Y轴固定范围 80 ~ 170
y_bmd_min_limit <- 80
y_bmd_max_limit <- 170
y_bmd_range <- y_bmd_max_limit - y_bmd_min_limit
y_breaks <- pretty(c(y_bmd_min_limit, y_bmd_max_limit), n = 5)
# 柱状图参数
y_freq_base <- y_bmd_min_limit + 0.02 * y_bmd_range
bar_total_height <- 0.2 * y_bmd_range
freq_max_target <- 0.1
# 拐点标注自动置顶
annotate_y <- y_bmd_max_limit - 0.03 * y_bmd_range
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# 双向缩放函数
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * bar_total_height
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / bar_total_height * freq_max_target
}
# 绘图
p <- ggplot() +
# 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 分段拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 拐点虚线
geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 拐点数值标注
annotate(
"text",
x = bp_vals,
y = annotate_y,
label = round(bp_vals, 3),
size = 3,
color = "gray20"
) +
# 双Y轴设置
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = y_breaks,
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
axis.title.y.right = element_text(size = 10, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成并组合两张亚组图
# ==============================
p_young <- plot_age_subgroup(df_young, seg_young, "A. 年龄 < 45岁")
p_old <- plot_age_subgroup(df_old, seg_old, "B. 年龄 ≥ 45岁")
final_plot_age <- p_young + p_old +
plot_annotation(
title = "SUA与BMD关联的年龄分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
)
print(final_plot_age)## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
# ==============================
# 【新增】导出高清PDF矢量图
# ==============================
ggsave("SUA_BMD年龄分层亚组分析图.pdf", final_plot_age, width = 12, height = 6, dpi = 300, device = "pdf")## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'A. 年龄 < 45岁' in 'mbcsToSbcs': for 年 (U+5E74)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'A. 年龄 < 45岁' in 'mbcsToSbcs': for 龄 (U+9F84)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'A. 年龄 < 45岁' in 'mbcsToSbcs': for 岁 (U+5C81)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'B. 年龄 ≥ 45岁' in 'mbcsToSbcs': for 年 (U+5E74)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'B. 年龄 ≥ 45岁' in 'mbcsToSbcs': for 龄 (U+9F84)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'B. 年龄 ≥ 45岁' in 'mbcsToSbcs': for 岁 (U+5C81)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 与
## (U+4E0E)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 关
## (U+5173)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 联
## (U+8054)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 的
## (U+7684)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 年
## (U+5E74)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 龄
## (U+9F84)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 层
## (U+5C42)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 亚
## (U+4E9A)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 组
## (U+7EC4)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的年龄分层亚组分析' in 'mbcsToSbcs': for 析
## (U+6790)8.2年龄分组-50岁划分（双拐点）
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（年龄分组50岁分界）
# ==============================
df_age <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
age_group = factor(
ifelse(Age < 50, "<50岁", "≥50岁"),
levels = c("<50岁", "≥50岁")
)
)
df_young <- df_age %>% filter(age_group == "<50岁")
df_old <- df_age %>% filter(age_group == "≥50岁")
# ==============================
# 2. 模型拟合函数 + 自动初始值
# ==============================
fit_seg_model <- function(data, psi_init) {
fit <- lm(BMD ~ SUA + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR, data = data)
segmented(fit, seg.Z = ~SUA, psi = list(SUA = psi_init))
}
get_psi_init <- function(data, probs = c(0.3, 0.7)) {
as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
}
psi_young <- get_psi_init(df_young)
psi_old <- get_psi_init(df_old)
seg_young <- fit_seg_model(df_young, psi_young)
seg_old <- fit_seg_model(df_old, psi_old)
# ==============================
# 3. 【修改】提取亚组结果表（新增P值，统一3位小数）
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n) {
# 提取拐点
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
# 提取斜率原始数据
slope_raw <- as.data.frame(slope(seg_model)$SUA)
# 自动识别列
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
# 提取数值
slope_ests <- slope_raw[[col_est]]
slope_se <- slope_raw[[col_se]]
# 计算P值
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- slope_ests / slope_se
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
tibble(
分组类型 = "年龄",
亚组 = group_name,
样本量 = n,
拐点1 = round(bp_vals[1], 3),
拐点2 = round(bp_vals[2], 3),
一区段斜率 = round(slope_ests[1], 3),
二区段斜率 = round(slope_ests[2], 3),
三区段斜率 = round(slope_ests[3], 3),
一区段P值 = round(p_vals[1], 3),
二区段P值 = round(p_vals[2], 3),
三区段P值 = round(p_vals[3], 3)
)
}
subgroup_table_age <- bind_rows(
extract_subgroup_result(seg_young, "<50岁", nrow(df_young)),
extract_subgroup_result(seg_old, "≥50岁", nrow(df_old))
)
cat("========== 年龄亚组分析结果表 ==========\n")## ========== 年龄亚组分析结果表 ==========print(subgroup_table_age, n = Inf)## # A tibble: 2 × 11
## 分组类型 亚组 样本量 拐点1 拐点2 一区段斜率 二区段斜率 三区段斜率 一区段P值
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 年龄 <50岁 2041 304 435. -0.006 0.065 -0.006 0.896
## 2 年龄 ≥50岁 5489 312. 353. 0.107 -0.021 0.039 0
## # ℹ 2 more variables: 二区段P值 <dbl>, 三区段P值 <dbl># ==============================
# 4. 交互作用检验
# ==============================
inter_model_age <- lm(
BMD ~ SUA * age_group + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = df_age
)
inter_result_age <- broom::tidy(inter_model_age) %>%
filter(str_detect(term, "SUA:age_group")) %>%
mutate(
交互项 = "SUA × 年龄组",
FDR校正P值 = round(p.adjust(p.value, method = "fdr"), 3)
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值) %>%
mutate(
回归系数 = round(回归系数, 3),
标准误 = round(标准误, 3),
P值 = round(P值, 3)
)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_age, n = Inf)## # A tibble: 1 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA × 年龄组 0.01 0.01 0.332 0.332# ==============================
# 5. 绘图函数（拐点标注3位小数 + 主轴自动适配）
# ==============================
plot_age_subgroup <- function(data, seg_model, title_name) {
# ---------- 5.1 生成预测数据 + 95%置信区间 ----------
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
BMI = mean(data$BMI, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
FBG = mean(data$FBG, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE),
eGFR = mean(data$eGFR, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
bp_vals <- seg_model$psi[, "Est."]
# ---------- 5.2 主Y轴自动计算 + 动态双轴缩放 ----------
y_fit_min <- min(plot_data$lwr, na.rm = TRUE)
y_fit_max <- max(plot_data$upr, na.rm = TRUE)
y_fit_range <- y_fit_max - y_fit_min
y_bmd_min_limit <- y_fit_min - 0.05 * y_fit_range
y_bmd_max_limit <- y_fit_max + 0.05 * y_fit_range
y_bmd_range <- y_bmd_max_limit - y_bmd_min_limit
y_breaks <- pretty(c(y_bmd_min_limit, y_bmd_max_limit), n = 5)
# 柱状图参数
y_freq_base <- y_bmd_min_limit + 0.02 * y_bmd_range
bar_total_height <- 0.2 * y_bmd_range
freq_max_target <- 0.1
freq_align_bmd <- y_freq_base + bar_total_height
# 拐点标注自动置顶
annotate_y <- y_bmd_max_limit - 0.03 * y_bmd_range
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# 双向缩放函数
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * bar_total_height
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / bar_total_height * freq_max_target
}
# ---------- 5.3 绘图（拐点标注3位小数） ----------
p <- ggplot() +
# 1. 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 2. 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 3. 分段拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 4. 拐点虚线
geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 5. 【修改】拐点数值标注 → 3位小数
annotate(
"text",
x = bp_vals,
y = annotate_y,
label = round(bp_vals, 3),
size = 3,
color = "gray20"
) +
# 6. 双Y轴设置
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = y_breaks,
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
# 7. 标签与主题
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
axis.title.y.right = element_text(size = 10, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成并组合两张亚组图
# ==============================
p_young <- plot_age_subgroup(df_young, seg_young, "A. 年龄 < 50岁")
p_old <- plot_age_subgroup(df_old, seg_old, "B. 年龄 ≥ 50岁")
final_plot_age <- p_young + p_old +
plot_annotation(
title = "SUA与BMD关联的年龄分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
)
print(final_plot_age)
8.3 年龄分组-50岁划分（AIC自动确定拐点数）
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（年龄分组50岁分界）
# ==============================
df_age <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
age_group = factor(
ifelse(Age < 50, "<50岁", "≥50岁"),
levels = c("<50岁", "≥50岁")
)
)
df_young <- df_age %>% filter(age_group == "<50岁")
df_old <- df_age %>% filter(age_group == "≥50岁")
# ==============================
# 2. AIC准则自动选择最优拐点数量
# ==============================
fit_optimal_segmented <- function(data, max_breakpoints = 3) {
# 基础线性模型（0个拐点）
base_fit <- lm(BMD ~ SUA + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR, data = data)
aic_table <- tibble(
拐点数量 = 0,
AIC = AIC(base_fit),
模型 = list(base_fit)
)
# 循环拟合1~max个拐点
for (k in 1:max_breakpoints) {
tryCatch({
# k个拐点对应k+1段，按分位数生成初始值
probs <- seq(1/(k+1), k/(k+1), length.out = k)
psi_init <- as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
# 拟合分段模型
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_init), it.max = 100)
# 校验：拐点数量正确、无NA、不在数据边界
bp_vals <- seg$psi[, "Est."]
valid <- length(bp_vals) == k && !any(is.na(bp_vals)) &&
min(bp_vals) > min(data$SUA, na.rm = TRUE) + 10 &&
max(bp_vals) < max(data$SUA, na.rm = TRUE) - 10
if (valid) {
aic_table <- bind_rows(aic_table, tibble(
拐点数量 = k,
AIC = AIC(seg),
模型 = list(seg)
))
}
}, error = function(e) {
# 拟合失败直接跳过
return(NULL)
})
}
# 按AIC从小到大排序，选第一个为最优
aic_table <- aic_table %>% arrange(AIC)
best_idx <- 1
best_model <- aic_table$模型[[best_idx]]
best_n_break <- aic_table$拐点数量[best_idx]
return(list(
best_model = best_model,
best_n_breakpoints = best_n_break,
# 【关键修正】不用select函数，直接提取列重建表格，彻底避开报错
aic_comparison = tibble(
拐点数量 = aic_table$拐点数量,
AIC = aic_table$AIC
),
all_models = aic_table
))
}
# 分别对两个亚组自动选最优模型
fit_young <- fit_optimal_segmented(df_young, max_breakpoints = 3)
fit_old <- fit_optimal_segmented(df_old, max_breakpoints = 3)## Warning: Breakpoint estimate(s) outdistanced to allow finite estimates and
## st.errsseg_young <- fit_young$best_model
seg_old <- fit_old$best_model
# 输出AIC比较结果
cat("===== <50岁组 拐点数量AIC比较 =====\n")## ===== <50岁组 拐点数量AIC比较 =====print(fit_young$aic_comparison, n = Inf)## # A tibble: 4 × 2
## 拐点数量 AIC
## <dbl> <dbl>
## 1 0 19602.
## 2 1 19604.
## 3 2 19604.
## 4 3 19606.cat(paste0("最优拐点数量：", fit_young$best_n_breakpoints, "\n\n"))## 最优拐点数量：0cat("===== ≥50岁组 拐点数量AIC比较 =====\n")## ===== ≥50岁组 拐点数量AIC比较 =====print(fit_old$aic_comparison, n = Inf)## # A tibble: 4 × 2
## 拐点数量 AIC
## <dbl> <dbl>
## 1 2 53288.
## 2 1 53290.
## 3 3 53291.
## 4 0 53295.cat(paste0("最优拐点数量：", fit_old$best_n_breakpoints, "\n\n"))## 最优拐点数量：2# ==============================
# 3. 自适应结果提取函数（兼容0~3个拐点）
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n, n_break) {
# ---------- 0个拐点（线性模型）单独处理 ----------
if (n_break == 0) {
coef_sum <- summary(seg_model)$coefficients
beta <- coef_sum["SUA", "Estimate"]
se <- coef_sum["SUA", "Std. Error"]
p_val <- coef_sum["SUA", "Pr(>|t|)"]
return(tibble(
分组类型 = "年龄",
亚组 = group_name,
模型类型 = "线性(0拐点)",
样本量 = n,
拐点1 = NA_real_,
拐点2 = NA_real_,
拐点3 = NA_real_,
一区段斜率 = round(beta, 3),
二区段斜率 = NA_real_,
三区段斜率 = NA_real_,
四区段斜率 = NA_real_,
一区段P值 = round(p_val, 3),
二区段P值 = NA_real_,
三区段P值 = NA_real_,
四区段P值 = NA_real_
))
}
# ---------- 1~3个拐点 ----------
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
slope_raw <- as.data.frame(slope(seg_model)$SUA)
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
slope_ests <- slope_raw[[col_est]]
slope_se <- slope_raw[[col_se]]
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- slope_ests / slope_se
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
# 统一补全到3个拐点、4段斜率，保证表格列数一致
bp_out <- c(bp_vals, rep(NA_real_, 3 - length(bp_vals)))
slope_out <- c(slope_ests, rep(NA_real_, 4 - length(slope_ests)))
p_out <- c(p_vals, rep(NA_real_, 4 - length(p_vals)))
tibble(
分组类型 = "年龄",
亚组 = group_name,
模型类型 = paste0(n_break, "拐点"),
样本量 = n,
拐点1 = round(bp_out[1], 3),
拐点2 = round(bp_out[2], 3),
拐点3 = round(bp_out[3], 3),
一区段斜率 = round(slope_out[1], 3),
二区段斜率 = round(slope_out[2], 3),
三区段斜率 = round(slope_out[3], 3),
四区段斜率 = round(slope_out[4], 3),
一区段P值 = round(p_out[1], 3),
二区段P值 = round(p_out[2], 3),
三区段P值 = round(p_out[3], 3),
四区段P值 = round(p_out[4], 3)
)
}
subgroup_table_age <- bind_rows(
extract_subgroup_result(seg_young, "<50岁", nrow(df_young), fit_young$best_n_breakpoints),
extract_subgroup_result(seg_old, "≥50岁", nrow(df_old), fit_old$best_n_breakpoints)
)
cat("========== 年龄亚组分析结果表 ==========\n")## ========== 年龄亚组分析结果表 ==========print(subgroup_table_age, n = Inf)## # A tibble: 2 × 15
## 分组类型 亚组 模型类型 样本量 拐点1 拐点2 拐点3 一区段斜率 二区段斜率
## <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 年龄 <50岁 线性(0拐点) 2041 NA NA NA 0.035 NA
## 2 年龄 ≥50岁 2拐点 5489 319 320 NA 0.118 -3.62
## # ℹ 6 more variables: 三区段斜率 <dbl>, 四区段斜率 <dbl>, 一区段P值 <dbl>,
## # 二区段P值 <dbl>, 三区段P值 <dbl>, 四区段P值 <dbl># ==============================
# 4. 交互作用检验
# ==============================
inter_model_age <- lm(
BMD ~ SUA * age_group + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = df_age
)
inter_result_age <- broom::tidy(inter_model_age) %>%
filter(str_detect(term, "SUA:age_group")) %>%
mutate(
交互项 = "SUA × 年龄组",
FDR校正P值 = round(p.adjust(p.value, method = "fdr"), 3)
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值) %>%
mutate(
回归系数 = round(回归系数, 3),
标准误 = round(标准误, 3),
P值 = round(P值, 3)
)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_age, n = Inf)## # A tibble: 1 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA × 年龄组 0.01 0.01 0.332 0.332# ==============================
# 5. 自适应绘图函数（兼容0~3个拐点 + 自动Y轴 + 双轴柱状图）
# ==============================
plot_age_subgroup <- function(data, seg_model, title_name, n_break) {
# ---------- 5.1 生成预测数据 + 95%置信区间 ----------
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
BMI = mean(data$BMI, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
FBG = mean(data$FBG, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE),
eGFR = mean(data$eGFR, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
# ---------- 5.2 主Y轴自动计算 + 双轴缩放 ----------
y_fit_min <- min(plot_data$lwr, na.rm = TRUE)
y_fit_max <- max(plot_data$upr, na.rm = TRUE)
y_fit_range <- y_fit_max - y_fit_min
y_bmd_min_limit <- y_fit_min - 0.05 * y_fit_range
y_bmd_max_limit <- y_fit_max + 0.05 * y_fit_range
y_bmd_range <- y_bmd_max_limit - y_bmd_min_limit
y_breaks <- pretty(c(y_bmd_min_limit, y_bmd_max_limit), n = 5)
# 柱状图参数
y_freq_base <- y_bmd_min_limit + 0.02 * y_bmd_range
bar_total_height <- 0.2 * y_bmd_range
freq_max_target <- 0.1
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * bar_total_height
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / bar_total_height * freq_max_target
}
# ---------- 5.3 基础图层（柱状图+置信带+拟合线） ----------
p <- ggplot() +
# 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
)
# ---------- 5.4 有拐点时添加虚线与标注 ----------
if (n_break > 0) {
bp_vals <- seg_model$psi[, "Est."]
annotate_y <- y_bmd_max_limit - 0.03 * y_bmd_range
# 拐点虚线
p <- p + geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
)
# 拐点数值标注（自动左右避让，避免重叠）
p <- p + annotate(
"text",
x = bp_vals,
y = annotate_y,
label = round(bp_vals, 3),
size = 3,
color = "gray20",
hjust = ifelse(bp_vals < median(bp_vals), 1.1, -0.1)
)
}
# ---------- 5.5 坐标轴与主题 ----------
p <- p +
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = y_breaks,
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
axis.title.y.right = element_text(size = 10, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成并组合两张亚组图
# ==============================
p_young <- plot_age_subgroup(df_young, seg_young, "A. 年龄 < 50岁", fit_young$best_n_breakpoints)
p_old <- plot_age_subgroup(df_old, seg_old, "B. 年龄 ≥ 50岁", fit_old$best_n_breakpoints)
final_plot_age <- p_young + p_old +
plot_annotation(
title = "SUA与BMD关联的年龄分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
)
print(final_plot_age)
8.4 年龄分组-60岁划分（AIC自动确定拐点数）
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（年龄分组60岁分界）
# ==============================
df_age <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
age_group = factor(
ifelse(Age < 60, "<60岁", "≥60岁"),
levels = c("<60岁", "≥60岁")
)
)
df_young <- df_age %>% filter(age_group == "<60岁")
df_old <- df_age %>% filter(age_group == "≥60岁")
# ==============================
# 2. 【核心修改】稳健模型拟合函数（双拐点失败自动降级单拐点）
# ==============================
# 生成分位数初始值（支持单/双拐点）
get_psi_init <- function(data, probs = c(0.3, 0.7)) {
as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
}
# 稳健拟合：优先双拐点，失败自动降级单拐点
fit_seg_robust <- function(data) {
# 基础线性模型
base_fit <- lm(BMD ~ SUA + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR, data = data)
# 先尝试拟合双拐点
tryCatch({
psi_double <- get_psi_init(data, c(0.3, 0.7))
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_double), it.max = 100)
bp_vals <- seg$psi[, "Est."]
# 校验：拐点无NA、数量为2、且不在数据边界
if (length(bp_vals) == 2 && !any(is.na(bp_vals)) &&
min(bp_vals) > min(data$SUA, na.rm = TRUE) + 10 &&
max(bp_vals) < max(data$SUA, na.rm = TRUE) - 10) {
return(list(model = seg, type = "双拐点"))
} else {
stop("双拐点拟合无效")
}
}, error = function(e) {
# 双拐点失败，降级为单拐点
psi_single <- get_psi_init(data, 0.5)
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_single), it.max = 100)
return(list(model = seg, type = "单拐点"))
})
}
# 分别拟合两组
fit_young <- fit_seg_robust(df_young)## Warning: Breakpoint estimate(s) outdistanced to allow finite estimates and
## st.errs## breakpoint estimate(s): 303.9999 304.9999fit_old <- fit_seg_robust(df_old)
seg_young <- fit_young$model
seg_old <- fit_old$model
# 输出拟合类型确认
cat("各组模型拟合结果：\n")## 各组模型拟合结果：cat(paste0("<60岁组：", fit_young$type, "\n"))## <60岁组：单拐点cat(paste0("≥60岁组：", fit_old$type, "\n\n"))## ≥60岁组：双拐点# ==============================
# 3. 提取亚组结果表（兼容单/双拐点，统一3位小数）
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n, model_type) {
# 提取拐点
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
n_bp <- length(bp_vals)
# 提取斜率原始数据
slope_raw <- as.data.frame(slope(seg_model)$SUA)
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
slope_ests <- slope_raw[[col_est]]
slope_se <- slope_raw[[col_se]]
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- slope_ests / slope_se
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
# 按拐点数量补全
bp_out <- c(bp_vals, rep(NA, 2 - n_bp))
slope_out <- c(slope_ests, rep(NA, 3 - length(slope_ests)))
p_out <- c(p_vals, rep(NA, 3 - length(p_vals)))
tibble(
分组类型 = "年龄",
亚组 = group_name,
模型类型 = model_type,
样本量 = n,
拐点1 = round(bp_out[1], 3),
拐点2 = round(bp_out[2], 3),
一区段斜率 = round(slope_out[1], 3),
二区段斜率 = round(slope_out[2], 3),
三区段斜率 = round(slope_out[3], 3),
一区段P值 = round(p_out[1], 3),
二区段P值 = round(p_out[2], 3),
三区段P值 = round(p_out[3], 3)
)
}
subgroup_table_age <- bind_rows(
extract_subgroup_result(seg_young, "<60岁", nrow(df_young), fit_young$type),
extract_subgroup_result(seg_old, "≥60岁", nrow(df_old), fit_old$type)
)
cat("========== 年龄亚组分析结果表 ==========\n")## ========== 年龄亚组分析结果表 ==========print(subgroup_table_age, n = Inf)## # A tibble: 2 × 12
## 分组类型 亚组 模型类型 样本量 拐点1 拐点2 一区段斜率 二区段斜率 三区段斜率
## <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 年龄 <60岁 单拐点 5142 450. NA 0.048 0.009 NA
## 2 年龄 ≥60岁 双拐点 2388 299. 322. 0.118 -0.11 0.037
## # ℹ 3 more variables: 一区段P值 <dbl>, 二区段P值 <dbl>, 三区段P值 <dbl># ==============================
# 4. 交互作用检验
# ==============================
inter_model_age <- lm(
BMD ~ SUA * age_group + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = df_age
)
inter_result_age <- broom::tidy(inter_model_age) %>%
filter(str_detect(term, "SUA:age_group")) %>%
mutate(
交互项 = "SUA × 年龄组",
FDR校正P值 = round(p.adjust(p.value, method = "fdr"), 3)
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值) %>%
mutate(
回归系数 = round(回归系数, 3),
标准误 = round(标准误, 3),
P值 = round(P值, 3)
)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_age, n = Inf)## # A tibble: 1 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA × 年龄组 0.009 0.01 0.388 0.388# ==============================
# 5. 绘图函数（兼容单/双拐点 + 自动Y轴 + 双轴柱状图）
# ==============================
plot_age_subgroup <- function(data, seg_model, title_name) {
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
BMI = mean(data$BMI, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
FBG = mean(data$FBG, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE),
eGFR = mean(data$eGFR, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
bp_vals <- seg_model$psi[, "Est."]
# 主Y轴自动计算
y_fit_min <- min(plot_data$lwr, na.rm = TRUE)
y_fit_max <- max(plot_data$upr, na.rm = TRUE)
y_fit_range <- y_fit_max - y_fit_min
y_bmd_min_limit <- y_fit_min - 0.05 * y_fit_range
y_bmd_max_limit <- y_fit_max + 0.05 * y_fit_range
y_bmd_range <- y_bmd_max_limit - y_bmd_min_limit
y_breaks <- pretty(c(y_bmd_min_limit, y_bmd_max_limit), n = 5)
# 柱状图参数
y_freq_base <- y_bmd_min_limit + 0.02 * y_bmd_range
bar_total_height <- 0.2 * y_bmd_range
freq_max_target <- 0.1
# 拐点标注自动置顶
annotate_y <- y_bmd_max_limit - 0.03 * y_bmd_range
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# 双向缩放函数
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * bar_total_height
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / bar_total_height * freq_max_target
}
# 绘图
p <- ggplot() +
# 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 分段拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 拐点虚线（自动适配数量）
geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 拐点数值标注（3位小数）
annotate(
"text",
x = bp_vals,
y = annotate_y,
label = round(bp_vals, 3),
size = 3,
color = "gray20"
) +
# 双Y轴设置
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = y_breaks,
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
axis.title.y.right = element_text(size = 10, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成并组合两张亚组图
# ==============================
p_young <- plot_age_subgroup(df_young, seg_young, "A. 年龄 < 60岁")
p_old <- plot_age_subgroup(df_old, seg_old, "B. 年龄 ≥ 60岁")
final_plot_age <- p_young + p_old +
plot_annotation(
title = "SUA与BMD关联的年龄分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
)
print(final_plot_age)
8.5 BMI分组-四组
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（BMI分4组）
# ==============================
df_BMI <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
BMI_group = factor(
case_when(
BMI < 18.5 ~ "<18.5",
BMI >= 18.5 & BMI < 24 ~ "18.5≤BMI<24",
BMI >= 24 & BMI < 28 ~ "24≤BMI<28",
BMI >= 28 ~ "≥28"
),
levels = c("<18.5", "18.5≤BMI<24", "24≤BMI<28", "≥28")
)
)
# 拆分4个亚组数据
df_bmi1 <- df_BMI %>% filter(BMI_group == "<18.5")
df_bmi2 <- df_BMI %>% filter(BMI_group == "18.5≤BMI<24")
df_bmi3 <- df_BMI %>% filter(BMI_group == "24≤BMI<28")
df_bmi4 <- df_BMI %>% filter(BMI_group == "≥28")
# ==============================
# 2. 模型拟合函数 + 自动生成分位数初始值
# ==============================
# 拟合双拐点分段模型（分层后模型不纳入BMI协变量）
fit_seg_model <- function(data, psi_init) {
fit <- lm(BMD ~ SUA + Age + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR, data = data)
segmented(fit, seg.Z = ~SUA, psi = list(SUA = psi_init))
}
# 自动根据每组数据分布生成拐点初始值（30% + 70%分位，双拐点）
get_psi_init <- function(data, probs = c(0.3, 0.7)) {
as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
}
# 分别生成4组的初始值并拟合模型
psi_bmi1 <- get_psi_init(df_bmi1)
psi_bmi2 <- get_psi_init(df_bmi2)
psi_bmi3 <- get_psi_init(df_bmi3)
psi_bmi4 <- get_psi_init(df_bmi4)
seg_bmi1 <- fit_seg_model(df_bmi1, psi_bmi1)## Warning: Breakpoint estimate(s) outdistanced to allow finite estimates and
## st.errsseg_bmi2 <- fit_seg_model(df_bmi2, psi_bmi2)
seg_bmi3 <- fit_seg_model(df_bmi3, psi_bmi3)
seg_bmi4 <- fit_seg_model(df_bmi4, psi_bmi4)
# ==============================
# 3. 提取亚组结果表
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n) {
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
slope_raw <- slope(seg_model)$SUA
if (is.matrix(slope_raw) || is.data.frame(slope_raw)) {
slope_ests <- slope_raw[, 1]
} else if (is.list(slope_raw)) {
slope_ests <- slope_raw[[1]][, 1]
} else {
slope_ests <- as.vector(slope_raw)
}
tibble(
分组类型 = "BMI",
亚组 = group_name,
样本量 = n,
拐点1 = round(bp_vals[1], 3),
拐点2 = round(bp_vals[2], 3),
一区段斜率 = round(slope_ests[1], 3),
二区段斜率 = round(slope_ests[2], 3),
三区段斜率 = round(slope_ests[3], 3)
)
}
subgroup_table_BMI <- bind_rows(
extract_subgroup_result(seg_bmi1, "<18.5", nrow(df_bmi1)),
extract_subgroup_result(seg_bmi2, "18.5≤BMI<24", nrow(df_bmi2)),
extract_subgroup_result(seg_bmi3, "24≤BMI<28", nrow(df_bmi3)),
extract_subgroup_result(seg_bmi4, "≥28", nrow(df_bmi4))
)
cat("========== BMI亚组分析结果表 ==========\n")## ========== BMI亚组分析结果表 ==========print(subgroup_table_BMI, n = Inf)## # A tibble: 4 × 8
## 分组类型 亚组 样本量 拐点1 拐点2 一区段斜率 二区段斜率 三区段斜率
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BMI <18.5 28 365 406. 0.326 -1.39 1.16
## 2 BMI 18.5≤BMI<24 1887 229. 444. 0.126 0.035 -0.074
## 3 BMI 24≤BMI<28 3924 278 484 0.062 0.013 -0.056
## 4 BMI ≥28 1691 280. 289. 0.181 -0.546 -0.003# ==============================
# 4. 交互作用检验（4组）
# ==============================
inter_model_BMI <- lm(
BMD ~ SUA * BMI_group + Age + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = df_BMI
)
inter_result_BMI <- broom::tidy(inter_model_BMI) %>%
filter(str_detect(term, "SUA:BMI_group")) %>%
mutate(
交互项 = term,
FDR校正P值 = p.adjust(p.value, method = "fdr")
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_BMI, n = Inf)## # A tibble: 3 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA:BMI_group18.5≤BMI<24 -0.205 0.0852 0.0161 0.0161
## 2 SUA:BMI_group24≤BMI<28 -0.213 0.0849 0.0123 0.0161
## 3 SUA:BMI_group≥28 -0.224 0.0852 0.00862 0.0161# ==============================
# 5. 绘图函数（双坐标轴+柱状图+全灰度，与主图风格完全一致）
# ==============================
plot_BMI_subgroup <- function(data, seg_model, title_name) {
# ---------- 5.1 生成预测数据 + 95%置信区间 ----------
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
Age = mean(data$Age, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
FBG = mean(data$FBG, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE),
eGFR = mean(data$eGFR, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
bp_vals <- seg_model$psi[, "Est."]
# ---------- 5.2 频率分箱计算 + 双轴缩放 ----------
# 坐标轴对齐参数（与主图/年龄亚组完全一致）
y_bmd_max_limit <- 130
y_bmd_min_limit <- 95
freq_max_target <- 0.1
freq_align_bmd <- 120
y_freq_base <- 100
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# 双向缩放函数
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * (freq_align_bmd - y_freq_base)
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / (freq_align_bmd - y_freq_base) * freq_max_target
}
# ---------- 5.3 绘图（全灰度+双轴+柱状图） ----------
p <- ggplot() +
# 1. 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 2. 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 3. 分段拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 4. 拐点虚线
geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 5. 拐点数值标注
annotate(
"text",
x = bp_vals,
y = 128,
label = round(bp_vals, 1),
size = 3,
color = "gray20"
) +
# 6. 双Y轴设置
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = seq(100, 130, 10),
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
# 7. 标签与主题
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 10, hjust = 0),
axis.title = element_text(size = 9),
axis.text = element_text(size = 8),
axis.title.y.right = element_text(size = 9, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成4张子图 + 2×2组合
# ==============================
p_bmi1 <- plot_BMI_subgroup(df_bmi1, seg_bmi1, "A. BMI < 18.5 kg/m²")
p_bmi2 <- plot_BMI_subgroup(df_bmi2, seg_bmi2, "B. 18.5≤BMI<24 kg/m²")
p_bmi3 <- plot_BMI_subgroup(df_bmi3, seg_bmi3, "C. 24≤BMI<28 kg/m²")
p_bmi4 <- plot_BMI_subgroup(df_bmi4, seg_bmi4, "D. BMI ≥ 28 kg/m²")
# 2行2列组合
final_plot_BMI <- (p_bmi1 + p_bmi2) / (p_bmi3 + p_bmi4) +
plot_annotation(
title = "SUA与BMD关联的BMI分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 13, hjust = 0.5))
)
# 显示最终图
print(final_plot_BMI)## Warning: Removed 100 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).## Warning: Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 19 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).## Warning: Removed 15 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
8.6 BMI分组-2组
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（BMI分2组：正常/超重）
# ==============================
df_BMI <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
BMI_group = factor(
case_when(
BMI >= 18.5 & BMI < 24 ~ "18.5≤BMI<24 kg/m²",
BMI >= 24 ~ "≥24 kg/m²"
),
levels = c("18.5≤BMI<24 kg/m²", "≥24 kg/m²")
)
)
# 拆分2个亚组数据
df_bmi_normal <- df_BMI %>% filter(BMI_group == "18.5≤BMI<24 kg/m²")
df_bmi_over <- df_BMI %>% filter(BMI_group == "≥24 kg/m²")
# ==============================
# 2. 稳健模型拟合函数（双拐点失败自动降级单拐点）
# ==============================
get_psi_init <- function(data, probs = c(0.3, 0.7)) {
as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
}
fit_seg_robust <- function(data) {
# 分层模型剔除BMI协变量（已按BMI分层）
base_fit <- lm(BMD ~ SUA + Age + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR, data = data)
# 先尝试拟合双拐点
tryCatch({
psi_double <- get_psi_init(data, c(0.3, 0.7))
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_double), it.max = 100)
bp_vals <- seg$psi[, "Est."]
# 校验：拐点无NA、数量为2、且不在数据边界
if (length(bp_vals) == 2 && !any(is.na(bp_vals)) &&
min(bp_vals) > min(data$SUA, na.rm = TRUE) + 10 &&
max(bp_vals) < max(data$SUA, na.rm = TRUE) - 10) {
return(list(model = seg, type = "双拐点"))
} else {
stop("双拐点拟合无效")
}
}, error = function(e) {
# 双拐点失败，降级为单拐点
psi_single <- get_psi_init(data, 0.5)
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_single), it.max = 100)
return(list(model = seg, type = "单拐点"))
})
}
# 分别拟合两组
fit_normal <- fit_seg_robust(df_bmi_normal)
fit_over <- fit_seg_robust(df_bmi_over)
seg_normal <- fit_normal$model
seg_over <- fit_over$model
# 输出拟合类型确认
cat("各组模型拟合结果：\n")## 各组模型拟合结果：cat(paste0("18.5≤BMI<24组：", fit_normal$type, "\n"))## 18.5≤BMI<24组：双拐点cat(paste0("≥24组：", fit_over$type, "\n\n"))## ≥24组：双拐点# ==============================
# 3. 提取亚组结果表（兼容单/双拐点，统一3位小数）
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n, model_type) {
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
n_bp <- length(bp_vals)
slope_raw <- as.data.frame(slope(seg_model)$SUA)
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
slope_ests <- slope_raw[[col_est]]
slope_se <- slope_raw[[col_se]]
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- slope_ests / slope_se
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
# 按拐点数量补全
bp_out <- c(bp_vals, rep(NA_real_, 2 - n_bp))
slope_out <- c(slope_ests, rep(NA_real_, 3 - length(slope_ests)))
p_out <- c(p_vals, rep(NA_real_, 3 - length(p_vals)))
tibble(
分组类型 = "BMI",
亚组 = group_name,
模型类型 = model_type,
样本量 = n,
拐点1 = round(bp_out[1], 3),
拐点2 = round(bp_out[2], 3),
一区段斜率 = round(slope_out[1], 3),
二区段斜率 = round(slope_out[2], 3),
三区段斜率 = round(slope_out[3], 3),
一区段P值 = round(p_out[1], 3),
二区段P值 = round(p_out[2], 3),
三区段P值 = round(p_out[3], 3)
)
}
subgroup_table_BMI <- bind_rows(
extract_subgroup_result(seg_normal, "18.5≤BMI<24 kg/m²", nrow(df_bmi_normal), fit_normal$type),
extract_subgroup_result(seg_over, "≥24 kg/m²", nrow(df_bmi_over), fit_over$type)
)
cat("========== BMI亚组分析结果表 ==========\n")## ========== BMI亚组分析结果表 ==========print(subgroup_table_BMI, n = Inf)## # A tibble: 2 × 12
## 分组类型 亚组 模型类型 样本量 拐点1 拐点2 一区段斜率 二区段斜率 三区段斜率
## <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BMI 18.5≤BM… 双拐点 1887 229. 444. 0.126 0.035 -0.074
## 2 BMI ≥24 kg/… 双拐点 5615 269. 497. 0.08 0.008 -0.04
## # ℹ 3 more variables: 一区段P值 <dbl>, 二区段P值 <dbl>, 三区段P值 <dbl># ==============================
# 4. 交互作用检验
# ==============================
inter_model_BMI <- lm(
BMD ~ SUA * BMI_group + Age + WC + SBP + DBP + FBG + TC + TG + HDL + LDL + eGFR,
data = df_BMI
)
inter_result_BMI <- broom::tidy(inter_model_BMI) %>%
filter(str_detect(term, "SUA:BMI_group")) %>%
mutate(
交互项 = "SUA × BMI组",
FDR校正P值 = round(p.adjust(p.value, method = "fdr"), 3)
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值) %>%
mutate(
回归系数 = round(回归系数, 3),
标准误 = round(标准误, 3),
P值 = round(P值, 3)
)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_BMI, n = Inf)## # A tibble: 1 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA × BMI组 -0.008 0.01 0.423 0.423# ==============================
# 5. 【核心修改】绘图函数（主Y轴固定 80~170）
# ==============================
plot_BMI_subgroup <- function(data, seg_model, title_name) {
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
Age = mean(data$Age, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
FBG = mean(data$FBG, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE),
eGFR = mean(data$eGFR, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
bp_vals <- seg_model$psi[, "Est."]
# ---------- 核心修改：主Y轴固定范围 80 ~ 170 ----------
y_bmd_min_limit <- 80
y_bmd_max_limit <- 170
y_bmd_range <- y_bmd_max_limit - y_bmd_min_limit
y_breaks <- pretty(c(y_bmd_min_limit, y_bmd_max_limit), n = 5)
# 柱状图参数（自动适配固定范围，保持比例一致）
y_freq_base <- y_bmd_min_limit + 0.02 * y_bmd_range
bar_total_height <- 0.2 * y_bmd_range
freq_max_target <- 0.1
# 拐点标注自动置顶
annotate_y <- y_bmd_max_limit - 0.03 * y_bmd_range
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# 双向缩放函数
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * bar_total_height
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / bar_total_height * freq_max_target
}
# 绘图
p <- ggplot() +
# 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 分段拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 拐点虚线（自动适配数量）
geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 拐点数值标注（3位小数）
annotate(
"text",
x = bp_vals,
y = annotate_y,
label = round(bp_vals, 3),
size = 3,
color = "gray20"
) +
# 双Y轴设置
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = y_breaks,
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
axis.title.y.right = element_text(size = 10, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成并组合两张亚组图
# ==============================
p_normal <- plot_BMI_subgroup(df_bmi_normal, seg_normal, "A. 18.5≤BMI<24 kg/m²")
p_over <- plot_BMI_subgroup(df_bmi_over, seg_over, "B. BMI ≥ 24 kg/m²")
final_plot_BMI <- p_normal + p_over +
plot_annotation(
title = "SUA与BMD关联的BMI分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
)
print(final_plot_BMI)## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
# ==============================
# 【新增】导出高清PDF矢量图
# ==============================
ggsave("SUA_BMD_BMI分层亚组分析图.pdf", final_plot_BMI, width = 12, height = 6, dpi = 300, device = "pdf")## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 与
## (U+4E0E)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 关
## (U+5173)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 联
## (U+8054)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 的
## (U+7684)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 层
## (U+5C42)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 亚
## (U+4E9A)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 组
## (U+7EC4)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的BMI分层亚组分析' in 'mbcsToSbcs': for 析
## (U+6790)8.7 FBG分组
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（FBG分2组：<7 / ≥7 mmol/L）
# ==============================
df_FBG <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
FBG_group = factor(
ifelse(FBG < 7, "<7 mmol/L", "≥7 mmol/L"),
levels = c("<7 mmol/L", "≥7 mmol/L")
)
)
# 拆分2个亚组数据
df_fbg_low <- df_FBG %>% filter(FBG_group == "<7 mmol/L")
df_fbg_high <- df_FBG %>% filter(FBG_group == "≥7 mmol/L")
# ==============================
# 2. 稳健模型拟合函数（双拐点失败自动降级单拐点）
# ==============================
get_psi_init <- function(data, probs = c(0.3, 0.7)) {
as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
}
fit_seg_robust <- function(data) {
# 分层模型剔除FBG协变量（已按FBG分层）
base_fit <- lm(BMD ~ SUA + Age + BMI + WC + SBP + DBP + TC + TG + HDL + LDL + eGFR, data = data)
# 先尝试拟合双拐点
tryCatch({
psi_double <- get_psi_init(data, c(0.3, 0.7))
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_double), it.max = 100)
bp_vals <- seg$psi[, "Est."]
# 校验：拐点无NA、数量为2、且不在数据边界
if (length(bp_vals) == 2 && !any(is.na(bp_vals)) &&
min(bp_vals) > min(data$SUA, na.rm = TRUE) + 10 &&
max(bp_vals) < max(data$SUA, na.rm = TRUE) - 10) {
return(list(model = seg, type = "双拐点"))
} else {
stop("双拐点拟合无效")
}
}, error = function(e) {
# 双拐点失败，降级为单拐点
psi_single <- get_psi_init(data, 0.5)
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_single), it.max = 100)
return(list(model = seg, type = "单拐点"))
})
}
# 分别拟合两组
fit_low <- fit_seg_robust(df_fbg_low)
fit_high <- fit_seg_robust(df_fbg_high)
seg_low <- fit_low$model
seg_high <- fit_high$model
# 输出拟合类型确认
cat("各组模型拟合结果：\n")## 各组模型拟合结果：cat(paste0("<7 mmol/L组：", fit_low$type, "\n"))## <7 mmol/L组：双拐点cat(paste0("≥7 mmol/L组：", fit_high$type, "\n\n"))## ≥7 mmol/L组：双拐点# ==============================
# 3. 提取亚组结果表（兼容单/双拐点，统一3位小数）
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n, model_type) {
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
n_bp <- length(bp_vals)
slope_raw <- as.data.frame(slope(seg_model)$SUA)
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
slope_ests <- slope_raw[[col_est]]
slope_se <- slope_raw[[col_se]]
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- slope_ests / slope_se
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
# 按拐点数量补全
bp_out <- c(bp_vals, rep(NA_real_, 2 - n_bp))
slope_out <- c(slope_ests, rep(NA_real_, 3 - length(slope_ests)))
p_out <- c(p_vals, rep(NA_real_, 3 - length(p_vals)))
tibble(
分组类型 = "FBG",
亚组 = group_name,
模型类型 = model_type,
样本量 = n,
拐点1 = round(bp_out[1], 3),
拐点2 = round(bp_out[2], 3),
一区段斜率 = round(slope_out[1], 3),
二区段斜率 = round(slope_out[2], 3),
三区段斜率 = round(slope_out[3], 3),
一区段P值 = round(p_out[1], 3),
二区段P值 = round(p_out[2], 3),
三区段P值 = round(p_out[3], 3)
)
}
subgroup_table_FBG <- bind_rows(
extract_subgroup_result(seg_low, "<7 mmol/L", nrow(df_fbg_low), fit_low$type),
extract_subgroup_result(seg_high, "≥7 mmol/L", nrow(df_fbg_high), fit_high$type)
)
cat("========== FBG亚组分析结果表 ==========\n")## ========== FBG亚组分析结果表 ==========print(subgroup_table_FBG, n = Inf)## # A tibble: 2 × 12
## 分组类型 亚组 模型类型 样本量 拐点1 拐点2 一区段斜率 二区段斜率 三区段斜率
## <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FBG <7 mmol… 双拐点 6655 246 490 0.105 0.016 -0.054
## 2 FBG ≥7 mmol… 双拐点 875 260. 278 0.002 0.41 -0.028
## # ℹ 3 more variables: 一区段P值 <dbl>, 二区段P值 <dbl>, 三区段P值 <dbl># ==============================
# 4. 交互作用检验
# ==============================
inter_model_FBG <- lm(
BMD ~ SUA * FBG_group + Age + BMI + WC + SBP + DBP + TC + TG + HDL + LDL + eGFR,
data = df_FBG
)
inter_result_FBG <- broom::tidy(inter_model_FBG) %>%
filter(str_detect(term, "SUA:FBG_group")) %>%
mutate(
交互项 = "SUA × FBG组",
FDR校正P值 = round(p.adjust(p.value, method = "fdr"), 3)
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值) %>%
mutate(
回归系数 = round(回归系数, 3),
标准误 = round(标准误, 3),
P值 = round(P值, 3)
)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_FBG, n = Inf)## # A tibble: 1 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA × FBG组 -0.018 0.012 0.158 0.158# ==============================
# 5. 【核心修改】绘图函数（主Y轴固定 80~170）
# ==============================
plot_FBG_subgroup <- function(data, seg_model, title_name) {
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
Age = mean(data$Age, na.rm = TRUE),
BMI = mean(data$BMI, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE),
eGFR = mean(data$eGFR, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
bp_vals <- seg_model$psi[, "Est."]
# ---------- 核心修改：主Y轴固定范围 80 ~ 170 ----------
y_bmd_min_limit <- 80
y_bmd_max_limit <- 170
y_bmd_range <- y_bmd_max_limit - y_bmd_min_limit
y_breaks <- pretty(c(y_bmd_min_limit, y_bmd_max_limit), n = 5)
# 柱状图参数（自动适配固定范围，保持比例一致）
y_freq_base <- y_bmd_min_limit + 0.02 * y_bmd_range
bar_total_height <- 0.2 * y_bmd_range
freq_max_target <- 0.1
# 拐点标注自动置顶
annotate_y <- y_bmd_max_limit - 0.03 * y_bmd_range
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# 双向缩放函数
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * bar_total_height
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / bar_total_height * freq_max_target
}
# 绘图
p <- ggplot() +
# 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 分段拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 拐点虚线（自动适配数量）
geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 拐点数值标注（3位小数）
annotate(
"text",
x = bp_vals,
y = annotate_y,
label = round(bp_vals, 3),
size = 3,
color = "gray20"
) +
# 双Y轴设置
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = y_breaks,
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
axis.title.y.right = element_text(size = 10, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成并组合两张亚组图
# ==============================
p_low <- plot_FBG_subgroup(df_fbg_low, seg_low, "A. FBG < 7 mmol/L")
p_high <- plot_FBG_subgroup(df_fbg_high, seg_high, "B. FBG ≥ 7 mmol/L")
final_plot_FBG <- p_low + p_high +
plot_annotation(
title = "SUA与BMD关联的FBG分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
)
print(final_plot_FBG)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
# ==============================
# 【新增】导出高清PDF矢量图
# ==============================
ggsave("SUA_BMD_FBG分层亚组分析图.pdf", final_plot_FBG, width = 12, height = 6, dpi = 300, device = "pdf")## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 与
## (U+4E0E)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 关
## (U+5173)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 联
## (U+8054)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 的
## (U+7684)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 层
## (U+5C42)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 亚
## (U+4E9A)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 组
## (U+7EC4)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的FBG分层亚组分析' in 'mbcsToSbcs': for 析
## (U+6790)8.8 egfr分组
library(tidyverse)
library(segmented)
library(patchwork)
# ==============================
# 1. 数据准备（eGFR分2组：<90 / ≥90 ml/min/1.73m²）
# ==============================
df_eGFR <- data1 %>%
drop_na(SUA, BMD, Age, BMI, WC, SBP, DBP, FBG, TC, TG, HDL, LDL, eGFR) %>%
mutate(
eGFR_group = factor(
ifelse(eGFR < 90, "<90 ml/min/1.73m²", "≥90 ml/min/1.73m²"),
levels = c("<90 ml/min/1.73m²", "≥90 ml/min/1.73m²")
)
)
# 拆分2个亚组数据
df_egfr_low <- df_eGFR %>% filter(eGFR_group == "<90 ml/min/1.73m²")
df_egfr_high <- df_eGFR %>% filter(eGFR_group == "≥90 ml/min/1.73m²")
# ==============================
# 2. 稳健模型拟合函数（双拐点失败自动降级单拐点）
# ==============================
get_psi_init <- function(data, probs = c(0.3, 0.7)) {
as.numeric(quantile(data$SUA, probs = probs, na.rm = TRUE))
}
fit_seg_robust <- function(data) {
# 分层模型剔除eGFR协变量（已按eGFR分层）
base_fit <- lm(BMD ~ SUA + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL, data = data)
# 先尝试拟合双拐点
tryCatch({
psi_double <- get_psi_init(data, c(0.3, 0.7))
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_double), it.max = 100)
bp_vals <- seg$psi[, "Est."]
# 校验：拐点无NA、数量为2、且不在数据边界
if (length(bp_vals) == 2 && !any(is.na(bp_vals)) &&
min(bp_vals) > min(data$SUA, na.rm = TRUE) + 10 &&
max(bp_vals) < max(data$SUA, na.rm = TRUE) - 10) {
return(list(model = seg, type = "双拐点"))
} else {
stop("双拐点拟合无效")
}
}, error = function(e) {
# 双拐点失败，降级为单拐点
psi_single <- get_psi_init(data, 0.5)
seg <- segmented(base_fit, seg.Z = ~SUA, psi = list(SUA = psi_single), it.max = 100)
return(list(model = seg, type = "单拐点"))
})
}
# 分别拟合两组
fit_low <- fit_seg_robust(df_egfr_low)
fit_high <- fit_seg_robust(df_egfr_high)
seg_low <- fit_low$model
seg_high <- fit_high$model
# 输出拟合类型确认
cat("各组模型拟合结果：\n")## 各组模型拟合结果：cat(paste0("<90 ml/min/1.73m²组：", fit_low$type, "\n"))## <90 ml/min/1.73m²组：双拐点cat(paste0("≥90 ml/min/1.73m²组：", fit_high$type, "\n\n"))## ≥90 ml/min/1.73m²组：双拐点# ==============================
# 3. 提取亚组结果表（兼容单/双拐点，统一3位小数）
# ==============================
extract_subgroup_result <- function(seg_model, group_name, n, model_type) {
bp_mat <- as.matrix(seg_model$psi)
bp_vals <- bp_mat[, "Est."]
n_bp <- length(bp_vals)
slope_raw <- as.data.frame(slope(seg_model)$SUA)
col_est <- which(str_detect(colnames(slope_raw), regex("Est", ignore_case = TRUE)))[1]
col_se <- which(str_detect(colnames(slope_raw), regex("Err|SE", ignore_case = TRUE)))[1]
col_p <- which(str_detect(colnames(slope_raw), regex("Pr|p-val", ignore_case = TRUE)))[1]
slope_ests <- slope_raw[[col_est]]
slope_se <- slope_raw[[col_se]]
if (!is.na(col_p) && length(col_p) > 0) {
p_vals <- slope_raw[[col_p]]
} else {
t_vals <- slope_ests / slope_se
p_vals <- 2 * (1 - pt(abs(t_vals), df = df.residual(seg_model)))
}
# 按拐点数量补全
bp_out <- c(bp_vals, rep(NA_real_, 2 - n_bp))
slope_out <- c(slope_ests, rep(NA_real_, 3 - length(slope_ests)))
p_out <- c(p_vals, rep(NA_real_, 3 - length(p_vals)))
tibble(
分组类型 = "eGFR",
亚组 = group_name,
模型类型 = model_type,
样本量 = n,
拐点1 = round(bp_out[1], 3),
拐点2 = round(bp_out[2], 3),
一区段斜率 = round(slope_out[1], 3),
二区段斜率 = round(slope_out[2], 3),
三区段斜率 = round(slope_out[3], 3),
一区段P值 = round(p_out[1], 3),
二区段P值 = round(p_out[2], 3),
三区段P值 = round(p_out[3], 3)
)
}
subgroup_table_eGFR <- bind_rows(
extract_subgroup_result(seg_low, "<90 ml/min/1.73m²", nrow(df_egfr_low), fit_low$type),
extract_subgroup_result(seg_high, "≥90 ml/min/1.73m²", nrow(df_egfr_high), fit_high$type)
)
cat("========== eGFR亚组分析结果表 ==========\n")## ========== eGFR亚组分析结果表 ==========print(subgroup_table_eGFR, n = Inf)## # A tibble: 2 × 12
## 分组类型 亚组 模型类型 样本量 拐点1 拐点2 一区段斜率 二区段斜率 三区段斜率
## <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 eGFR <90 ml/… 双拐点 2346 308 349. 0.1 -0.059 0.001
## 2 eGFR ≥90 ml/… 双拐点 5184 244 480 0.071 0.024 -0.047
## # ℹ 3 more variables: 一区段P值 <dbl>, 二区段P值 <dbl>, 三区段P值 <dbl># ==============================
# 4. 交互作用检验
# ==============================
inter_model_eGFR <- lm(
BMD ~ SUA * eGFR_group + Age + BMI + WC + SBP + DBP + FBG + TC + TG + HDL + LDL,
data = df_eGFR
)
inter_result_eGFR <- broom::tidy(inter_model_eGFR) %>%
filter(str_detect(term, "SUA:eGFR_group")) %>%
mutate(
交互项 = "SUA × eGFR组",
FDR校正P值 = round(p.adjust(p.value, method = "fdr"), 3)
) %>%
dplyr::select(交互项, 回归系数 = estimate, 标准误 = std.error, P值 = p.value, FDR校正P值) %>%
mutate(
回归系数 = round(回归系数, 3),
标准误 = round(标准误, 3),
P值 = round(P值, 3)
)
cat("\n========== 交互作用检验结果 ==========\n")##
## ========== 交互作用检验结果 ==========print(inter_result_eGFR, n = Inf)## # A tibble: 1 × 5
## 交互项 回归系数 标准误 P值 FDR校正P值
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SUA × eGFR组 0.009 0.009 0.327 0.327# ==============================
# 5. 【核心修改】绘图函数（主Y轴固定 80~170）
# ==============================
plot_eGFR_subgroup <- function(data, seg_model, title_name) {
new_x <- seq(min(data$SUA, na.rm = TRUE), max(data$SUA, na.rm = TRUE), length.out = 100)
newdata <- data.frame(
SUA = new_x,
Age = mean(data$Age, na.rm = TRUE),
BMI = mean(data$BMI, na.rm = TRUE),
WC = mean(data$WC, na.rm = TRUE),
SBP = mean(data$SBP, na.rm = TRUE),
DBP = mean(data$DBP, na.rm = TRUE),
FBG = mean(data$FBG, na.rm = TRUE),
TC = mean(data$TC, na.rm = TRUE),
TG = mean(data$TG, na.rm = TRUE),
HDL = mean(data$HDL, na.rm = TRUE),
LDL = mean(data$LDL, na.rm = TRUE)
)
pred_vals <- predict(seg_model, newdata = newdata, interval = "confidence", level = 0.95)
plot_data <- tibble(
SUA = new_x,
fit = pred_vals[, "fit"],
lwr = pred_vals[, "lwr"],
upr = pred_vals[, "upr"]
)
bp_vals <- seg_model$psi[, "Est."]
# ---------- 核心修改：主Y轴固定范围 80 ~ 170 ----------
y_bmd_min_limit <- 80
y_bmd_max_limit <- 170
y_bmd_range <- y_bmd_max_limit - y_bmd_min_limit
y_breaks <- pretty(c(y_bmd_min_limit, y_bmd_max_limit), n = 5)
# 柱状图参数（自动适配固定范围，保持比例一致）
y_freq_base <- y_bmd_min_limit + 0.02 * y_bmd_range
bar_total_height <- 0.2 * y_bmd_range
freq_max_target <- 0.1
# 拐点标注自动置顶
annotate_y <- y_bmd_max_limit - 0.03 * y_bmd_range
# 计算分箱数据
h <- hist(data$SUA, breaks = 30, plot = FALSE)
hist_data <- tibble(
x_left = h$breaks[-length(h$breaks)],
x_right = h$breaks[-1],
count = h$counts,
freq = count / sum(count)
)
# 双向缩放函数
scale_freq_to_bmd <- function(freq) {
y_freq_base + (freq / freq_max_target) * bar_total_height
}
scale_bmd_to_freq <- function(y) {
(y - y_freq_base) / bar_total_height * freq_max_target
}
# 绘图
p <- ggplot() +
# 底层频率柱状图
geom_rect(
data = hist_data,
aes(
xmin = x_left,
xmax = x_right,
ymin = y_freq_base,
ymax = scale_freq_to_bmd(freq)
),
fill = "gray80",
color = "gray60",
alpha = 0.9
) +
# 95%置信区间带
geom_ribbon(
data = plot_data,
aes(x = SUA, ymin = lwr, ymax = upr),
fill = "gray60",
alpha = 0.3
) +
# 分段拟合曲线
geom_line(
data = plot_data,
aes(x = SUA, y = fit),
linewidth = 1.2,
color = "gray20"
) +
# 拐点虚线（自动适配数量）
geom_vline(
xintercept = bp_vals,
linetype = "dashed",
color = "gray30",
linewidth = 0.8
) +
# 拐点数值标注（3位小数）
annotate(
"text",
x = bp_vals,
y = annotate_y,
label = round(bp_vals, 3),
size = 3,
color = "gray20"
) +
# 双Y轴设置
scale_y_continuous(
name = "拟合骨密度 BMD (mg/cm³)",
limits = c(y_bmd_min_limit, y_bmd_max_limit),
breaks = y_breaks,
expand = c(0, 0),
sec.axis = sec_axis(
~ scale_bmd_to_freq(.),
name = "频率（Frequency）",
breaks = seq(0, 0.1, 0.025),
labels = scales::percent_format(accuracy = 1)
)
) +
scale_x_continuous(
limits = c(min(data$SUA, na.rm = TRUE) - 5, max(data$SUA, na.rm = TRUE) + 5),
expand = c(0, 0)
) +
labs(
x = "血尿酸 SUA (μmol/L)",
title = title_name
) +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
axis.title.y.right = element_text(size = 10, angle = 270, vjust = 1.5),
panel.grid = element_blank()
)
return(p)
}
# ==============================
# 6. 生成并组合两张亚组图
# ==============================
p_low <- plot_eGFR_subgroup(df_egfr_low, seg_low, "A. eGFR < 90 ml/min/1.73m²")
p_high <- plot_eGFR_subgroup(df_egfr_high, seg_high, "B. eGFR ≥ 90 ml/min/1.73m²")
final_plot_eGFR <- p_low + p_high +
plot_annotation(
title = "SUA与BMD关联的eGFR分层亚组分析",
theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
)
print(final_plot_eGFR)
# ==============================
# 【新增】导出高清PDF矢量图
# ==============================
ggsave("SUA_BMD_eGFR分层亚组分析图.pdf", final_plot_eGFR, width = 12, height = 6, dpi = 300, device = "pdf")## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 血 (U+8840)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 尿 (U+5C3F)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '血尿酸 SUA (μmol/L)' in 'mbcsToSbcs': for 酸 (U+9178)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 拟
## (U+62DF)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 合
## (U+5408)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 骨
## (U+9AA8)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 密
## (U+5BC6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '拟合骨密度 BMD (mg/cm³)' in 'mbcsToSbcs': for 度
## (U+5EA6)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 频 (U+9891)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on '频率（Frequency）' in 'mbcsToSbcs': for 率 (U+7387)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 与
## (U+4E0E)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 关
## (U+5173)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 联
## (U+8054)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 的
## (U+7684)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 层
## (U+5C42)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 亚
## (U+4E9A)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 组
## (U+7EC4)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 分
## (U+5206)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'SUA与BMD关联的eGFR分层亚组分析' in 'mbcsToSbcs': for 析
## (U+6790)