Medical_code/paper_code/code.R

# install.packages("msm", repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")

library(msm)
library(survival)
library(dplyr)

# View(cav)

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


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

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


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

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

#心理得分
data$psychiatric_group <- cut(data$Psychiatric_score, breaks = c(0, 10, Inf),labels = c("无抑郁", "有抑郁"), right = FALSE)

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

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

summary(data)
View(data[, c("ADL", "variable")])
View(data$age_group)
# 计算状态转移频数表
freq_table <- statetable.msm(state, ID, data = data)
print(freq_table)

# 初始化转移速率矩阵
qmatrix_init <- matrix(c(-0.5, 0.25, 0.15, 0.1,
                          0.1, -0.3, 0.1, 0.1,
                          0.3, 0.1, -0.5, 0.1,
                          0,   0,    0,   0), 
                        nrow = 4, byrow = TRUE)

# 创建初始模型
crude_init <- crudeinits.msm(state ~ wave, subject = ID, data = data, qmatrix = qmatrix_init)

View(crude_init)

# 进行多状态模型分析
msm_model <- msm(state ~ wave, subject = ID, data = data,
                 qmatrix = crude_init,
                 covariates = ~ rgender_group+age_group+marital_group+education_group+activity_group+psychiatric_group+BMI_group+ADL_group+Smoke+Drink+last_year_O3+last_year_pm1+last_year_NO2+last_year_NH4+last_year_nl,
                 death = 4,
                 method = "BFGS", control = list(fnscale = 5000, maxit = 10000)
                 )


# 获取模型结果并转换为字符串
model_summary <- capture.output(summary(msm_model))
# 创建一个空的结果文件
result_file <- "msm_model_results_test.txt"
write("", file = result_file)  # 清空文件内容
# 写入文件，附加协变量名称
cat("Results for covariates:", file = result_file, append = TRUE)
cat(model_summary, file = result_file, sep = "\n", append = TRUE)


# 查看模型的详细结果
summary(msm_model)

# 计算状态转移概率矩阵
prob_matrix <- pmatrix.msm(msm_model, t = 5)  # t = 1 代表随访之间的间隔时间
print(prob_matrix)

# 输出拟合模型的速率矩阵
q_matrix <- qmatrix.msm(msm_model)
print(q_matrix)

# 提取转移强度
transition_intensity <- msm_model$qmatrix
print(transition_intensity)

# 计算在每个状态中的平均逗留时间
so_journ <- sojourn.msm(msm_model)
print(so_journ)

# 计算均衡状态概率

rm(list = ls())


# 定义协变量列表
covariates_list <- list(~last_year_NO3, ~before_last_NO3, ~last_year_NH4, ~before_last_NH4, ~last_year_OM, ~before_last_OM, ~last_year_BC, ~before_last_BC, ~last_year_nl, ~before_last_nl)

# 创建一个空的结果文件
result_file <- "msm_model_results.txt"
write("", file = result_file)  # 清空文件内容

# 循环计算不同协变量的模型
for (cov in covariates_list) {
  
  # 运行msm模型
  msm_model <- msm(state ~ wave, subject = ID, data = data,
                   qmatrix = crude_init,
                   covariates = cov,
                   death = 4,
                   method = "BFGS", control = list(fnscale = 4000, maxit = 10000)
                   )
  
  # 获取模型结果并转换为字符串
  model_summary <- capture.output(summary(msm_model))
  
  # 写入文件，附加协变量名称
  cat("Results for covariates:", deparse(cov), "\n", file = result_file, append = TRUE)
  cat(model_summary, file = result_file, sep = "\n", append = TRUE)
  cat("\n\n", file = result_file, append = TRUE)  # 空行分隔每个协变量结果
}