修改四状态马尔可夫链模拟，记录首次进入状态4的患者

我在R中构建了一个四状态马尔可夫链，可视化代码如下：

grViz("
  digraph {
    rankdir=LR;
    node [shape = circle]
    1 [label = '1', style = filled, fillcolor = green]
    2 [label = '2']
    3 [label = '3']
    4 [label = '4', style = filled, fillcolor = red]
    
    1 -> 1 [label = '1/3']
    1 -> 2 [label = '1/3']
    1 -> 3 [label = '1/3']
    2 -> 1 [label = '1/2']
    2 -> 2 [label = '1/2']
    3 -> 3 [label = '1/2']
    3 -> 4 [label = '1/2']
    4 -> 4 [label = '1']
  }
")

已完成包含100次迭代、100名初始状态为1的患者、重复1000次的模拟，代码如下：

library(reshape2)
library(ggplot2)
library(dplyr)
library(tidyr)

# 转移矩阵
P <- matrix(c(1/3, 1/3, 1/3, 0,
              1/2, 1/2, 0,   0,
              0,   0,   1/2, 1/2,
              0,   0,   0,   1),
            nrow = 4, ncol = 4, byrow = TRUE)

n_iterations <- 100
n_people <- 100
n_simulations <- 100

run_simulation <- function(sim_number) {
    state_matrix <- matrix(0, nrow = n_iterations, ncol = n_people)
    state_matrix[1, ] <- 1  # 所有患者初始状态为1
    
    for (i in 2:n_iterations) {
        for (person in 1:n_people) {
            current_state <- state_matrix[i-1, person]
            state_matrix[i, person] <- sample(1:4, 1, prob = P[current_state, ])
        }
    }
    
    state_counts <- apply(state_matrix, 1, function(x) table(factor(x, levels = 1:4)))
    percentages <- as.data.frame(t(state_counts) / n_people * 100)  # 转换为百分比（0-100）
    
    # 生成患者汇总表
    patient_summary <- data.frame(
        patient = 1:n_people,
        simulation = sim_number,
        state1 = rowSums(state_matrix == 1),
        state2 = rowSums(state_matrix == 2),
        state3 = rowSums(state_matrix == 3),
        state4 = rowSums(state_matrix == 4)
    )
    
    list(percentages = percentages, patient_summary = patient_summary)
}

# 运行模拟
set.seed(123)  # 保证结果可复现
all_simulations <- vector("list", n_simulations)
all_patient_summaries <- vector("list", n_simulations)

for (sim in 1:n_simulations) {
    simulation_result <- run_simulation(sim)
    all_simulations[[sim]] <- simulation_result$percentages
    all_patient_summaries[[sim]] <- simulation_result$patient_summary
    if (sim %% 10 == 0) {  # 每完成10次模拟打印进度
        cat(sprintf("已完成第%d次模拟，共%d次\n", sim, n_simulations))
    }
}

# 合并所有患者汇总表
final_patient_summary <- do.call(rbind, all_patient_summaries)

# 调整列顺序，将patient和simulation放在前面
final_patient_summary <- final_patient_summary[, c("patient", "simulation", "state1", "state2", "state3", "state4")]

# 计算所有模拟的平均百分比
average_percentages <- Reduce(`+`, all_simulations) / n_simulations
average_percentages$iteration <- 1:n_iterations

# 计算置信区间的函数
calculate_ci <- function(x) {
    quantile(x, probs = c(0.05, 0.95))
}

# 生成置信区间数据
ci_data <- lapply(1:n_iterations, function(i) {
    data.frame(
        iteration = i,
        state = as.character(1:4),
        lower = sapply(1:4, function(j) calculate_ci(sapply(all_simulations, function(sim) sim[i, j]))[1]),
        upper = sapply(1:4, function(j) calculate_ci(sapply(all_simulations, function(sim) sim[i, j]))[2])
    )
})
ci_data <- do.call(rbind, ci_data)

# 转换为长格式
average_percentages_long <- average_percentages %>%
    pivot_longer(cols = as.character(1:4),
                 names_to = "state",
                 values_to = "percentage")

plot_data <- merge(average_percentages_long, ci_data, by = c("iteration", "state"))

ggplot(plot_data, aes(x = iteration, y = percentage, color = state, fill = state)) +
    geom_line() +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
    scale_color_manual(values = c("1" = "blue", "2" = "red", "3" = "green", "4" = "purple"),
                       labels = paste("状态", 1:4)) +
    scale_fill_manual(values = c("1" = "blue", "2" = "red", "3" = "green", "4" = "purple"),
                      labels = paste("状态", 1:4)) +
    labs(x = "迭代次数", y = "各状态下人群平均百分比", color = "状态", fill = "状态") +
    ggtitle("1000次模拟中各状态人群占比随时间变化的平均值") +
    theme_minimal() +
    scale_y_continuous(limits = c(0, 100))

print(head(final_patient_summary))

同时已完成模拟结果的可视化。

需要修改上述模拟代码，在每次迭代中记录以下内容：

• 首次进入状态4的患者（标记为state_4_new）
• 已处于状态4的患者（标记为state_4_existing）

最终生成如下格式的记录：

simulation iteration patient state
          1        1      87     ...
          1        1      88     ...
          ...      ...      ...   ...

其中state字段的取值为：state1、state2、state3、state_4_new、state_4_existing。

内容的提问来源于stack exchange，提问作者farrow90