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

library(msm)
library(survival)

# View(cav)

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

# data$age_group <- cut(data$age,
#                       breaks = c(45, 55, 65, Inf),
#                       labels = c("45-54", "55-64",">=65"),
#                       right = FALSE)

data$gender_ <- factor(data$rgender, levels = c(1, 2), labels = c("男", "女"))

summary(data$age_group)
# View(data[, c("age", "age_group")])
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 = ~ gender_, 
                 death = 4,
                 method = "BFGS", control = list(fnscale = 4000, maxit = 10000)
                 )

# 计算状态转移概率矩阵
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)

# 计算均衡状态概率

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

rm(list = ls())