Medical_code/paper_code/figure.r

# 基线表
# install.packages("compareGroups", repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
# install.packages("rmarkdown", repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")

library(compareGroups)
library(dplyr)

data <- read.csv("paper_data.csv")

#性别
data$Rgender <- factor(data$rgender, levels = c(1, 2), labels = c("Male", "Female"))

#年齡
data$Age <- cut(data$age,
                      breaks = c(45, 55, 65, Inf),
                      labels = c("45-54", "55-64",">=65"),
                      right = FALSE)
#婚姻
data$Marital_Status <- factor(data$marital_status, levels = c(1, 0), labels = c("married or partnered", "other marital status"))

#教育
data$Education <- factor(data$education, levels = c(0, 1, 2), labels = c("below high school", "high school", "college or above"))

#运动情况
data$Physical_Activity <- factor(data$Physical_activity, levels = c(0, 1, 2), labels = c("inactive", "moderate", "vigorous"))

#心理得分
data$Depressive_Symptoms <- cut(data$Psychiatric_score, breaks = c(0, 10, Inf),labels = c("No", "Yes"), right = FALSE)

#BMI
# 使用 cut() 创建因子类型的变量
data$BMI <- cut(data$BMI, breaks = c(0, 18.5, 24, Inf), labels = c("underweight", "normal", "overweight"), right = FALSE)
# 使用 dplyr::recode 直接替换因子水平
data$BMI <- recode(data$BMI, "underweight" = 1L, "normal" = 0L, "overweight" = 2L)
# 将 data$variable 转换回因子并保留原 levels
data$BMI <- factor(data$BMI, levels = c(0, 1, 2), labels = c("normal", "underweight", "overweight"))

#ADL
# 使用 cut() 创建因子类型的变量
data$ADL <- cut(data$ADL, breaks = c(0, 1, 3, Inf), labels = c("No impairment", "Mild impairment", "Severe impairment"), right = FALSE)

#抽烟
data$Smoke <- factor(data$Smoke, levels = c(0, 1), labels = c("No", "Yes"))

#喝酒
data$Drink <- factor(data$Drink, levels = c(0, 1), labels = c("No", "Yes"))

#结局
data$State <- factor(data$state, levels = c(1, 2, 3, 4), labels = c("No chronic disease", "One chronic disease", "Comorbidity", "Death"))

# 定义分层变量和汇总变量
vars <- c("State", "Age", "Rgender", "Marital_Status", "Education", "BMI","ADL","Smoke","Drink","Physical_Activity", "Depressive_Symptoms", "last_year_NO2", "last_year_O3")

formula = State ~ Age+Rgender+Marital_Status+Education+BMI+ADL+Smoke+Drink+Physical_Activity+Depressive_Symptoms+last_year_NO2+last_year_O3
# 创建基线表
table = descrTable(formula = formula, data = data, show.p.overall = FALSE)

export2word(table, file = "基线表.docx")



# 森林图
library(dplyr)
library(ggforce)
library(ggforestplot)
df_compare_traits_groups <- read.csv("output_new.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)

df_logodds <-
  df_compare_traits_groups %>%
  dplyr::arrange(num) %>%
  # Set the study variable to a factor to preserve order of appearance
  # Set class to factor to set order of display.
  dplyr::mutate(
    trait = factor(
      trait,
      levels = c("No chronic disease to One chronic disease", "No chronic disease to Comorbidity",
     "One chronic disease to Comorbidity", "One chronic disease to No chronic disease", "Comorbidity to One chronic disease",
     "Comorbidity to No chronic disease", "No chronic disease to Death", "One chronic disease to Death", "Comorbidity to Death")
    )
  )

# Draw a forestplot of cross-sectional, linear associations.
p=forestplot(
  df = df_logodds,
  estimate = beta,
  # pvalue = pvalue,
  # psignif = 0,
  logodds = TRUE,
  title = "The impact results of covariates",
  xlab = "Effect Sizie with 95% CI",
  colour = trait,
  shape=trait,
  xlim = c(0.1, 6),
  # You can explicitly define x-tick breaks
  xtickbreaks = c(0.1, 0.5, 0.8, 1.0,1.2, 1.5, 2.0, 3.0,4.0,5.0)
) +
    ggplot2::scale_shape_manual(
    values = c(6L, 6L, 6L, 2L, 2L,2L,21L,21L,21L),
    labels = c("No chronic disease to One chronic disease", "No chronic disease to Comorbidity",
     "One chronic disease to Comorbidity", "One chronic disease to No chronic disease", "Comorbidity to One chronic disease",
     "Comorbidity to No chronic disease", "No chronic disease to Death", "One chronic disease to Death", "Comorbidity to Death")
  )+
  # 修改字体大小
  ggplot2::theme(
    plot.title = ggplot2::element_text(size = 32, face = "bold"),  # 标题字体
    axis.title = ggplot2::element_text(size = 28),  # 轴标题字体
    axis.text = ggplot2::element_text(size = 28),  # 轴刻度字体
    strip.text = ggplot2::element_text(size = 30),  # facet标签字体
    legend.text = ggplot2::element_text(size = 24),  # 图例字体
    legend.title = element_blank(),  # 移除图例标题
    legend.position = "top",  # 默认位置，右边
  )+
  # 使用 facet_wrap 进行分组显示，按 trait 进行分组
  ggplot2::facet_wrap(~group, scales = "free_y", ncol = 1)

ggsave("forestplot_output.png", plot = p, width = 35, height = 30, dpi = 300, bg = "white")


# library(ggplot2)
# library(dplyr)

# Assuming df_logodds is your data
# df_logodds <- df_compare_traits_groups %>%
#   dplyr::arrange(num) %>%
#   dplyr::mutate(
#     trait = factor(
#       trait,
#       levels = c("No chronic disease to One chronic disease", "No chronic disease to Comorbidity",
#                  "One chronic disease to Comorbidity", "One chronic disease to No chronic disease", 
#                  "Comorbidity to One chronic disease", "Comorbidity to No chronic disease", 
#                  "No chronic disease to Death", "One chronic disease to Death", "Comorbidity to Death")
#     )
#   )

# # Draw the forestplot using ggplot2
# p <- ggplot(df_logodds, aes(x = beta, y = trait, color = trait, shape = trait)) +
#   # Add points (the estimates)
#   geom_point(size = 5) +
#   # Add error bars (95% CI)
#   geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
#   # Customizing the shape of points
#   scale_shape_manual(
#     values = c(6, 6, 6, 2, 2, 2, 21, 21, 21),
#     labels = c("No chronic disease to One chronic disease", "No chronic disease to Comorbidity",
#                "One chronic disease to Comorbidity", "One chronic disease to No chronic disease", 
#                "Comorbidity to One chronic disease", "Comorbidity to No chronic disease", 
#                "No chronic disease to Death", "One chronic disease to Death", "Comorbidity to Death")
#   ) +
#   # Set the color palette for different traits
#   scale_color_manual(values = c("red", "blue", "green", "purple", "orange", "brown", "pink", "yellow", "gray")) +
#   # Add title and axis labels
#   labs(
#     title = "The impact results of covariates",
#     x = "Effect Size with 95% CI",
#     y = NULL  # No y-axis title, it's already covered by trait names
#   ) +
#   # Adjust the x-axis limits to control the effect size range
#   xlim(0.1, 6) +
#   # Customize text size and legend position
#   theme(
#     plot.title = element_text(size = 32, face = "bold"),  # Title font size and style
#     axis.title = element_text(size = 28),  # Axis title font size
#     axis.text = element_text(size = 28),  # Axis ticks font size
#     strip.text = element_text(size = 30),  # Facet label font size (if any)
#     legend.text = element_text(size = 24),  # Legend font size
#     legend.title = element_blank(),  # Remove legend title
#     legend.position = "top",  # Position of the legend
#     legend.key.size = unit(1.5, "lines"),  # Adjust the size of the legend keys (points)
#     panel.grid.major = element_line(size = 0.5, linetype = "dashed", color = "gray90"),
#     panel.grid.minor = element_blank(),  # Remove minor gridlines
#     panel.background = element_blank(),  # White background for panel
#     plot.background = element_blank()  # White background for plot
#   ) +
#   # Customize x-axis breaks
#   scale_x_continuous(
#     breaks = c(0.1, 0.5, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0, 4.0, 5.0)
#   )

# # Print the plot
# print(p)