SQEIR模型

使用R实现SQEIR模型,模拟流行病传播过程


日期

SQEIR模型

1 概述

SQEIR模型是一种用于描述传染病在人群中传播的模型。通过将人群划分为不同的健康状态仓室,模拟疾病传播的动态过程。SQEIR模型常用于传染病的流行趋势预测、干预效果评估及疫苗接种策略分析等方面。

1.1 模型假设

1.在疾病观察期间,该地区的总人口保持不变;

2.每个人被感染的可能性相同,仅考虑人际传播,传播概率均等;

3.隔离仓室Q和恢复仓室R的变化假定通过自愈或治疗实现,病例不会在短时间内转化为易感病例。

4.被隔离病例包括部分接触者密切接触者和感染者,被隔离病例不会传染给他人;

5.潜伏病例在潜伏期后被检测为确诊病例;

6.潜伏和确诊病例都具有传染性,感染率相同;

7.确诊病例得到治疗,最终治愈或死亡。

8.假设感染者通过直接接触或者空气传播将病毒传播给易感者。

9.假定恢复者一旦获得免疫就会持续维持,不考虑免疫衰退。

1.2 模型描述

SQEIR模型是一个常用的传染病动力学模型,常用于描述和模拟疾病传播的过程。 它扩展了传统的SIR模型(易感-感染-恢复模型),引入了更多的状态变量,能够更精确地刻画疫情的传播过程。 它通常包括:易感者、隔离者、潜伏者、感染者和康复者。 SQEIR模型可以广泛应用于多种传染病的传播分析,包括但不限于流感、COVID-19、SARS等。它的优势在于可以更细致地描述疫情传播的不同阶段,尤其是在对高风险人群(如被隔离者)进行控制时,模型能够帮助预测不同干预措施对疫情发展的影响。

1.3 模型仓室及参数含义

参数 含义
β 每单位时间的有效传播率
ω 当易感个体被感染时,首先进入潜伏期或潜伏状态,从潜伏个体到感染个体的转化率
q 隔离率(暴露个体转为隔离状态的比例)率
p 表现出症状并发展为确诊患者的感染者的比例
γ_1 感染者康复率
γ_2 隔离感染者康复率
μ 感染者的隔离率
$alpha^E$ 单位时间内暴露者被隔离的速率
$alpha^I$ 单位时间内感染者被隔离的速率

SQEIR 模型将总人口分为以下几类:

易感 (S):未感染但有机会感染疾病的人;

隔离暴露病例 (Q_E):暴露者中未被确诊且被隔离的人;

隔离感染病例 (Q_I)):感染者中被隔离的病例;

暴露 (E):指处于潜伏期并在潜伏期后将被感染的个体;

感染者 (I):指具有传染性并出现症状的个体;

康复 (R):指在感染者中已康复的个体;

2 流程图及方程

$$ \frac{dS}{dt} = -\beta S I $$

$$ \frac{dE}{dt} = \beta S I - \alpha^E qE - (1-\alpha^E)(1-q)p\omega E $$

$$ \frac{dQ_E}{dt} = \alpha^E qE - p\omega Q_E $$

$$ \frac{dI}{dt} = (1-\alpha^E)(1-q)p\omega E - (1-\mu)\gamma_1 I - \alpha^I \mu I $$

$$ \frac{dQ_I}{dt} = p\omega Q_E + \alpha^I \mu I - \gamma_2 Q_I $$

$$ \frac{dR}{dt} = (1-\mu)\gamma_1 I + \gamma_2 Q_I $$

3 适用疾病流行条件

SQEIR模型适用于以下疾病的流行研究:

1.传染性较强的疾病:如流感、COVID-19等疾病,能够通过人与人之间的接触传播。

2.潜伏期和感染期可分, 且需要隔离的疾病:并且适用于隔离环境封闭性较强的场景,例如集中隔离或医疗隔离。

4 代码(R)Running Code

4.1 代码

library(deSolve)
library(ggplot2)
library(tidyr)
library(plotly)
# 定义初始参数
parameters <- list(
  beta = 0.6,       # 感染率
  omega = 0.2,     # 潜伏期的逆
  q = 0.6,          # 暴露者隔离率
  p = 0.5,          # 隔离成功比例
  gamma1 = 0.1,     # 感染者康复率
  gamma2 = 0.05,     # 隔离感染者康复率
  mu = 0.1,         # 感染者的隔离率
  alphaE = 0.7,     # 暴露者隔离速率
  alphaI = 0.8      # 感染者隔离速率
)



# 定义初始状态
init_state <- c(
  S = 9800,         # 易感个体数量
  E = 150,          # 暴露个体数量
  QE = 100,         # 隔离暴露个体数量
  I = 50,          # 感染者数量
  QI = 30,          # 隔离感染者数量
  R = 20             # 康复者数量
)


# 定义时间范围
time <- seq(0, 60, by = 1)  

# 定义微分方程
sqeir_model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    # 微分方程
    dS <- -beta * S * I
    dE <- beta * S * I - alphaE * q * E - (1 - alphaE) * (1 - q) * p * omega * E
    dQE <- alphaE * q * E - p * omega * QE
    dI <- (1 - alphaE) * (1 - q) * p * omega * E - (1 - mu) * gamma1 * I - alphaI * mu * I
    dQI <- p * omega * QE + alphaI * mu * I - gamma2 * QI
    dR <- (1 - mu) * gamma1 * I + gamma2 * QI
    
    # 返回结果
    list(c(dS, dE, dQE, dI, dQI, dR))
  })
}

# 求解微分方程
output <- ode(
  y = init_state, 
  times = time, 
  func = sqeir_model, 
  parms = parameters
)

# 将结果转换为数据框
output_df <- as.data.frame(output)

# 查看结果
head(output_df)


    
# 绘制结果

output_long <- tidyr::pivot_longer(output_df, cols = -time, names_to = "Compartment", values_to = "Count")

ggplot(output_long, aes(x = time, y = Count, color = Compartment)) +
  geom_line(size = 1.2) +
  labs(title = "SQEIR 模型动态变化 ", x = "时间", y = "人数", color = "仓室") +
  theme_minimal() +
  scale_color_manual(values = c(
    "S" = "#A3C1DA",  # 易感个体(浅蓝色)
    "E" = "#00FFFF",  # 暴露个体(浅青色)
    "QE" = "#F1D0A9", # 隔离暴露个体(浅橙色)
    "I" = "#F9B9B6",  # 感染者(浅红色)
    "QI" = "#D2A1C7", # 隔离感染者(浅紫色)
    "R" = "#B0E0A4"   #康复者(浅绿色)
  )) +
  theme(
    legend.position = "top", 
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 16),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12)
  ) 




### 敏感性分析
library(plotly)
# 定义敏感性分析的参数及其取值范围
sensitivity_parameters <- list(
  q = c(0.4, 0.6, 0.8),    # 暴露者隔离率
  p = c(0.4, 0.5, 0.6),    # 隔离成功比例
  mu = c(0.1, 0.2, 0.3)    # 感染者的隔离率
)

# 定义敏感性分析函数(单个参数)
run_sensitivity <- function(param_name, param_values, base_parameters, init_state, time) {
  results_list <- list()
  
  for (val in param_values) {
    # 更新参数
    current_parameters <- base_parameters
    current_parameters[[param_name]] <- val
    
    # 求解微分方程
    output <- ode(
      y = init_state, 
      times = time, 
      func = sqeir_model, 
      parms = current_parameters
    )
    
    # 转换为数据框并添加参数值信息
    output_df <- as.data.frame(output)
    output_df$Parameter <- param_name
    output_df$Value <- val
    results_list[[length(results_list) + 1]] <- output_df
  }
  
  # 合并所有结果
  all_results <- do.call(rbind, results_list)
  return(all_results)
}

# 初始化一个空的数据框来存储所有灵敏度分析结果
all_sensitivity_results <- data.frame()

# 对每个参数进行灵敏度分析
for (param in names(sensitivity_parameters)) {
  param_values <- sensitivity_parameters[[param]]
  sensitivity_output <- run_sensitivity(param, param_values, parameters, init_state, time)
  all_sensitivity_results <- rbind(all_sensitivity_results, sensitivity_output)
}

# 将结果转换为长格式以便绘图
all_sensitivity_long <- pivot_longer(
  all_sensitivity_results, 
  cols = c("S", "E", "QE", "I", "QI", "R"), 
  names_to = "Compartment", 
  values_to = "Count"
)

# 可视化灵敏度分析结果
p_sensitivity <- ggplot(all_sensitivity_long, aes(x = time, y = Count, color = Compartment)) +
  geom_line(size = 1) +
  facet_grid(Parameter ~ Value, scales = "free_y") +
  labs(
    title = "SQEIR 模型灵敏度分析",
    x = "时间 (天)",
    y = "人数",
    color = "类别"
  ) +
  theme_minimal() +
  scale_color_manual(values = c(
    "S" = "#A3C1DA",  # 易感个体(浅蓝色)
    "E" = "#00FFFF",  # 暴露个体(浅青色)
    "QE" = "#F1D0A9", # 隔离暴露个体(浅橙色)
    "I" = "#F9B9B6", # 感染者(浅红色)
    "QI" = "#D2A1C7",# 隔离感染者(浅紫色)
    "R" = "#B0E0A4"   # 康复者(浅绿色)
  )) +
  theme(
    legend.position = "bottom", 
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 16),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    strip.text = element_text(size = 12)
  )

# 显示灵敏度分析图
print(p_sensitivity)

4.2结果分析

S(易感个体)数量随着时间快速下降,从初始的接近 10,000 人迅速减少到几乎为零。随着疾病的传播,易感人群逐步被感染并转变为暴露个体(E)或进入其他状态

E(暴露个体)数量在开始时快速上升,达到一个峰值后下降,最终趋于零。随着易感个体接触感染者,感染人数增多,暴露人群增加。潜伏期结束或隔离措施的实施,暴露个体逐步转变为感染者或隔离暴露个体。

QE(隔离暴露个体)的数量先快速上升,达到峰值后逐渐减少,最终接近于零。

I(感染者)的数量在开始时快速上升,达到峰值后逐步下降,最终趋于零。

QI(隔离感染者)的数量先逐步增加,达到峰值后下降,最终接近于零。

R(康复者)的数量随着时间不断增加,最终占据总人口的绝大多数。

1.初期阶段(0-10天)):暴露者和感染者数量迅速增加,隔离暴露者和隔离感染者数量也在逐步积累。 疫情处于快速扩散阶段。

2.疫情高峰期(10-30 天):感染者数量达到峰值并开始下降,同时暴露者数量迅速减少,表明疫情逐渐被控制。 隔离暴露者和隔离感染者数量达到峰值,说明隔离措施的作用明显。

3.疫情后期(30-60 天):易感者数量趋于稳定,感染者和隔离感染者几乎消失,康复者数量接近总人口。 疫情基本结束。

隔离暴露者和隔离感染者数量的峰值表明隔离策略对疫情控制的影响显著。 提高隔离率可能进一步减少感染者的峰值,同时提高康复率和隔离速率可以更快控制感染。

4.3敏感性结果分析

更高的 q 值(隔离更多的暴露个体)显著减少了暴露者向感染者的转化,暴露者 E 的峰值显著下降,感染者 I 的峰值略有下降,从而降低了疫情的传播速度和峰值。

提高隔离成功比例(p)能有效控制社区传播,减少未被隔离的感染者,从而降低疫情高峰。

更高的 mu 值意味着更多的感染者被及时隔离,提高隔离率对抑制疫情传播有显著作用。有效减少了社区传播和感染者的峰值。

整体来看,优先提高感染者隔离率 mu(增加感染者的隔离措施),如加速隔离转移、提升隔离效率,对控制疫情峰值最为关键。加强接触者追踪和隔离政策,能够有效减少潜伏期内的传播.确保隔离环境的管理和执行,提升隔离的实际效果,进一步降低暴露者对感染者的转化。

需改进之处

1.如果有实际数据,可以将实际值与模拟值进行对比,评估模型的拟合度。

2.可模拟不同隔离策略或疫苗接种策略对疫情的影响,进一步提升模型的实用性。

王嘉鹤
励志拯救塞尔达但只会上树的呀哈哈大姐头