分年龄组SEIAVR模型

使用MATLAB实现分年龄组SEIAVR模型β半自动化拟合，并模拟流行病传播过程

林宇豪

日期

SEIAVR模型

1 概述

1.1 模型假设

1.2 模型及分年龄组拟合描述

SEIAVR分年龄组模型既考虑了疫苗的作用，同时也考虑了无症状患者、暴露者在传染病流行过程中的作用，更贴合实际流行过程中的情况。但不同年龄组之间的传播能力难以量化，因而需要拟合传播能力($ \beta $)矩阵，这也是分年龄组模型的关键点。本人开发了一个半自动分年龄组($ \beta$ )矩阵拟合机，基于单次拟合的代码（首先得能够拟合），主要实现给出不同初始的( $\beta $)值，按照指定的次数进行多次拟合，之后可根据结果，选择误差最小的一组初始值检查效果。

1.3 模型仓室及参数含义

表 1: 仓室及其含义

仓室 含义
( $S$ ) 易感者人数
( $V_s$ ) 疫苗接种者（仍可能被感染）
( $V_r$ ) 疫苗生效者
( $E$ ) 暴露者人数（已感染但尚未传染他人）
( $ I$ ) 感染者人数（有症状且具有传染性）
( $A $ ) 无症状感染者人数（无症状但具有传染性）
( $ R$ ) 康复者人数
( $ N$ ) 总人口数 ($ N = S + E + I + A + R + V_s + V_r$)

表 2: 参数及其含义

参数 含义
( $\beta$ ) 传播率，不同年龄组人群间的传播能力（矩阵）
( $ \gamma $ ) 感染者的康复率
( $\gamma_1$ ) 无症状感染者的康复率
( $\omega$ ) 暴露者转化为感染者的速率
( $\omega_1$ ) 暴露者转化为无症状感染者的速率
( $p$ ) 暴露者转化为无症状感染者的比例
( $k$ ) 无症状感染者的相对传播能力系数
( $\phi$ ) 疫苗接种速率（向量）
( $f$ ) 抗体生成系数
( $q$ ) 接种疫苗但未产生抗体的比例（向量）
( $br$ ) 出生率
( $dr$ ) 死亡率

2 模型流程图及方程

2.1 模型流程图

以下是SEIAVR模型的流程图，展示了易感者、接种疫苗者、感染者和恢复者之间的转化过程：

SEIAVR模型流程图

图 1: 分年龄组SEIAVR模型流程图

2.2 模型方程

SEIAVR 模型由以下常微分方程组描述：

$$ \begin{cases} \frac{d\mathbf{S}}{dt} = br \cdot \mathbf{N} - dr \cdot \mathbf{S} - \boldsymbol{\phi} \odot \mathbf{S} -\boldsymbol{\beta}^\top \cdot \frac{\mathbf{I} + k \cdot \mathbf{A}}{\mathbf{N}} \odot \mathbf{S} \\
\frac{d\mathbf{E}}{dt} = \boldsymbol{\beta}^\top \cdot \frac{\mathbf{I} + k \cdot \mathbf{A}}{\mathbf{N}} \odot \mathbf{S} + q \odot \boldsymbol{\beta}^\top \cdot \frac{\mathbf{I} + k \cdot \mathbf{A}}{\mathbf{N}} \odot \mathbf{V_s} - p \cdot \omega_1 \cdot \mathbf{E} - (1 - p) \cdot \omega \cdot \mathbf{E} - dr \cdot \mathbf{E} \\
\frac{d\mathbf{I}}{dt} = (1 - p) \cdot \omega \cdot \mathbf{E} - dr \cdot \mathbf{I} - \gamma \cdot \mathbf{I} \\
\frac{d\mathbf{A}}{dt} = p \cdot \omega_1 \cdot \mathbf{E} - dr \cdot \mathbf{A} - \gamma_1 \cdot \mathbf{A} \\
\frac{d\mathbf{R}}{dt} = \gamma_1 \cdot \mathbf{A} + \gamma \cdot \mathbf{I} - dr \cdot \mathbf{R} \\
\frac{d\mathbf{V_s}}{dt} = \boldsymbol{\phi} \odot \mathbf{S} - dr \cdot \mathbf{V_s} - f \cdot \mathbf{V_s} - q \odot \boldsymbol{\beta}^\top \cdot \frac{\mathbf{I} + k \cdot \mathbf{A}}{\mathbf{N}} \odot \mathbf{V_s} \\
\frac{d\mathbf{V_r}}{dt} = f \cdot \mathbf{V_s} - dr \cdot \mathbf{V_r} \end{cases}$$

3 适用疾病流行条件

4 代码（Matlab）

4.1 条件设定

4.2 代码

%%
% 定义循环参数
n = 6;
aValues = linspace(0.01, 1.3, 51);  

% 预分配结果存储
results = zeros(length(aValues), 2);  % 第一列存储 a，第二列存储 finalErrorValue

% 开始循环
for i = 1:length(aValues)
    % 设置当前 a 值
    clear finalErrorValue ;
    a = aValues(i);

 % 在主脚本中
    aInitial = a * ones(n);  % 设置 a
    assignin('base', 'aInitial', aInitial);  % 将其放入 base 工作区
    run('plotFitting.m');  
    results(i, :) = [a, finalErrorValue];

    % 显示当前的 a 和最终误差值
    fprintf('a = %.2f, Final Error Value = %.4f\n', a, finalErrorValue);
end

% 显示所有结果
disp('All results (a and Final Error Value):');
disp(results);


% close all; clear; clc;
close all;

% 设置模型参数
n = 6;  % 年龄组数量
load('VaccinationTable.mat') % 加载疫苗接种数据
load('processedTable.mat') % 加载疫情数据

% 模拟时间范围
params.tSpan = [0, 58] + 0.5; % 从第0.5天到第58.5天，共59天
params.k  = 0.50; % 无症状患者传播能力系数
params.omega  = 0.53; % 有症状患者潜伏期结束速率（单位：1/天）
params.omega1  = 0.83; % 无症状患者潜伏期结束速率（单位：1/天）
params.p = 0.333; % 无症状感染者比例
params.gamma  = 0.31; % 有症状感染者传染期结束速率
params.gamma1  = 0.25; % 无症状感染者传染期结束速率
params.f = 0.1; % 疫苗接种后产生抗体的速率
params.q = [0.42 0.45 0.45 0.41 0.42 0.42]'; % 每个年龄组接种疫苗但未产生抗体的比例
params.dr = .00651/365; % 自然死亡率（按日计算）
params.br = .00567/365; % 自然出生率（按日计算）
params.n = 6; % 年龄组数量

% 振幅初始值猜测
aValues = linspace(0.01, 1.3, 101); 
aInitial = aValues(3) * ones(n);

% 设置初始猜测值
b =  0 * ones(n); % 初始化参数b
c = [aInitial(:); b(:)]; % 将aInitial和b合并为一个向量

% 提取疫情和疫苗接种数据
E0 = T{1125, 2:7}'; % 初始潜伏者数量（从疫情表中提取）
I0 = E0 / (1 - params.p) / params.omega; % 计算初始有症状感染者数量
A0 = I0 * params.p; % 初始无症状感染者数量
R0 = [0, 0, 0, 0, 0, 0]'; % 初始康复者数量
pre27VacNum = sum(VT{948:1115 - 1 - 1/params.f, 2:7})'; % 疫情流行前27天的疫苗接种数量
pre10VacNum = sum(VT{1115 - 1/params.f:1114, 2:7})'; % 疫情流行前10天的疫苗接种数量

% 人口数据（按年龄组划分）
population = [103.26, 129.98, 72.28, 1413.68, 397.8, 67.3]' * 1e4; % 各年龄组人口数量
totalVacNum = sum(VT{948:1173, 2:7})'; % 总疫苗接种数量
totalVacRate = totalVacNum ./ population; % 各年龄组的疫苗接种率
popuVacRate = sum(totalVacNum) / sum(population); % 总人口疫苗接种率
epiVacNum = totalVacNum - pre27VacNum - pre10VacNum; % 流行季节接种的疫苗数量
day = days(VT.Date(1173) - VT.Date(1115)) + 1; % 流行季节天数
params.fai = epiVacNum / day ./ population; % 流行季节每日疫苗接种率

% 计算初始状态
Vs0 = pre27VacNum .* params.q + pre10VacNum; % 接种但未产生抗体者数量
Vr0 = pre27VacNum .* (1 - params.q); % 接种并产生抗体者数量
S0 = population - E0 - I0 - A0 - Vs0 - Vr0; % 易感者数量
params.initialStates = [S0; E0; I0; A0; R0; Vs0; Vr0]; % 初始状态向量

%% 拟合β
observedIncidence = T{1125:1183,2:7};% 观察到的发病数据
objectiveFunction = @(unknowns) computeParameterLikelihood2(reshape(unknowns(1:n^2),n,n), reshape(unknowns(1+n^2:end),n,n),params, observedIncidence);

% options = optimoptions('fmincon', 'MaxFunctionEvaluations', 1e5, 'MaxIterations', ...
    % 1e3, 'PlotFcn', 'optimplotfval'); 
options = optimoptions('fmincon', 'MaxFunctionEvaluations', 1e5, 'MaxIterations', 1e3); 
[d, finalErrorValue] = fmincon(objectiveFunction, c(:), [], [], [], [], ...
    [0*ones(n^2,1), -2*pi*ones(n^2,1)], [1.3*ones(n^2,1); 2*pi*ones(n^2,1)], [], options);
disp(finalErrorValue)% 显示最终误差值
clear options;

aEstimated = reshape(d(1:n^2), n, n);% 振幅矩阵a
bEstimated = reshape(d(1+n^2:end), n, n);% 相位矩阵b

% 在初始的Beta下模拟得到模型预测值
dxdt = @(t,x) computeDerivative2(aEstimated, bEstimated, params, t, x);
% [tt, xPredicted] = ode23s(dxdt, params.tSpan, params.initialStates); 
[tt, xPredicted] = ode45(dxdt, params.tSpan, params.initialStates);% 使用ODE45求解微分方程
yPredicted = (1 - params.p) * params.omega * xPredicted(:, n+1:2*n); % 预测的每日新增病例数


%%% 取消下面注释内容后查看拟合的状态
% figure; heatmap(aEstimated)
% figure; heatmap(bEstimated)
% fig = figure;
% fig.Position = [248.2,291.4,1374.4,584];
% tiledlayout(2,1)
% ax1 = nexttile; plot(T{1125:1183,2:7});
% legend('age1','age2','age3','age4','age5','age6')
% ax2 = nexttile; plot(tt, yPredicted);
% legend('age1','age2','age3','age4','age5','age6')
% ax2.YLim = ax1.YLim;


%% 画出拟合效果图

% 年龄组名称
ageGroups = {'age0-5', 'age6-12', 'age13-18', 'age19-59', 'age60-79', 'age80+'};
% 实际数据对应的时间范围
timeActual = datetime(2023,02,16) + caldays(0:(1183-1125));

% 设置颜色（每个年龄组对应一个颜色）
colors = ["#f7df87", "#eca680", "#7ac7e2", "#e3716e", "#f3deb7", "#54beaa"];

% 创建绘图窗口
fig = figure;
fig.Position = [248.2, 291.4, 1374.4, 584]; % 设置窗口大小和位置
tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 创建1行2列的布局

% 第一个子图：按年龄组绘制实际数据与拟合数据

% 子图1：绘制实际数据
nexttile;
hActual = plot(timeActual, T{1125:1183, 2:7}); % 绘制实际观察数据（每列对应一个年龄组）
hold on;

% 设置实际数据的线条宽度和颜色
for i = 1:length(hActual)
    set(hActual(i), 'LineWidth', 1.5, 'Color', colors(i)); % 为每个年龄组分配颜色
end

% 设置X轴时间范围
xlim([datetime(2023,02,16) datetime(2023,04,15)]);
datetick('x', 'yyyy-mm-dd', 'keeplimits'); % 格式化日期
xlabel('Date'); % 设置X轴标签

% 设置Y轴标签和范围
ylabel('New Cases'); % 设置Y轴标签
ylim([0, max(ylim) * 1.1]); 
yticks1 = get(gca, 'YTick'); % 获取当前Y轴的刻度
yticklabels1 = strcat(string(yticks1 / 1000), 'k'); % 将刻度转化为千人单位
set(gca, 'YTickLabel', yticklabels1); % 应用千人单位刻度标签

% 子图1：绘制拟合数据
timePredicted = datetime(2023,02,16) + days(tt - tt(1)); % 将预测时间转化为实际时间
hFitted = plot(timePredicted, yPredicted, '--'); % 绘制拟合的每日新增病例数

% 设置拟合数据的线条宽度和颜色
for i = 1:length(hActual)
    set(hFitted(i), 'Color', colors(i), 'LineWidth', 1.5); % 为每个年龄组分配颜色
end

% 添加图例
legendLabels = [strcat(ageGroups, ' (Actual)'), strcat(ageGroups, ' (Fitted)')];
legend([hActual; hFitted], legendLabels, 'Location', 'northeast');

% 添加标注
text(datetime(2023,02,17), max(ylim) * 0.95, 'A', 'FontSize', 14, 'FontWeight', 'bold'); % 标注"A"

% 计算并显示每个年龄组的R²值

% 初始化R²存储数组
R2_values = zeros(1, 6); % 用于存储每个年龄组的R²值

% 遍历每个年龄组，计算R²值
for i = 1:6
    % 将拟合值插值到与实际数据对齐的时间点
    alignedFitted = interp1(timePredicted, yPredicted(:, i), timeActual, 'linear', 'extrap');

    % 提取实际值
    actualValues = T{1125:1183, i + 1}; % 提取实际数据（1125到1183行是对应的时间段）
    ResidualSumSquares = sum((actualValues - alignedFitted').^2); % 残差平方和
    TotalSumSquares = sum((actualValues - mean(actualValues)).^2); % 总平方和

    % 计算R²值（避免除以0）
    if TotalSumSquares ~= 0
        R2_values(i) = 1 - ResidualSumSquares / TotalSumSquares;
    else
        R2_values(i) = 0; % 如果总平方和为0，直接将R²设为0
    end
end

% 将每个年龄组的R²值添加到图中
text(datetime(2023,04,06), max(T{1125:1183, 2}) * 0.87, ...
    sprintf('R^2 = %.2f', R2_values(1)), ...
    'Color', colors(1), 'FontSize', 12, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'center');

text(datetime(2023,04,06), max(T{1125:1183, 3}) * 0.85, ...
    sprintf('R^2 = %.2f', R2_values(2)), ...
    'Color', colors(2), 'FontSize', 12, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'center');

text(datetime(2023,04,06), max(T{1125:1183, 5}) * 0.60, ...
    sprintf('R^2 = %.2f', R2_values(3)), ...
    'Color', colors(3), 'FontSize', 12, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'center');

text(datetime(2023,04,06), max(T{1125:1183, 5}) * 0.65, ...
    sprintf('R^2 = %.2f', R2_values(4)), ...
    'Color', colors(4), 'FontSize', 12, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'center');

text(datetime(2023,04,06), max(T{1125:1183, 5}) * 0.45, ...
    sprintf('R^2 = %.2f', R2_values(5)), ...
    'Color', colors(5), 'FontSize', 12, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'center');

text(datetime(2023,04,06), max(T{1125:1183, 5}) * 0.4, ...
    sprintf('R^2 = %.2f', R2_values(6)), ...
    'Color', colors(6), 'FontSize', 12, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'center');

% 第二个子图：总体实际数据与拟合数据

% 子图2：绘制总体数据
totalT = sum(T{1125:1183, 2:7}, 2); % 实际的总新增病例数（所有年龄组总和）
totalyPredicted = sum(yPredicted, 2); % 拟合的总新增病例数

nexttile;
hTotalActual = plot(timeActual, totalT, 'Color', [0.349, 0.612, 0.749], 'LineWidth', 1.5); % 绘制实际数据
hold on;
hTotalFitted = plot(timePredicted, totalyPredicted, 'Color', [0.761, 0.341, 0.349], 'LineStyle', '--', 'LineWidth', 1.5); % 绘制拟合数据

% 设置X轴
xlim([datetime(2023,02,16) datetime(2023,04,15)]);
datetick('x', 'yyyy-mm-dd', 'keeplimits');
xlabel('Date'); % 设置X轴标签

% 设置Y轴单位为千人（k）
yticks2 = get(gca, 'YTick');
yticklabels2 = strcat(string(yticks2 / 1000), 'k');
set(gca, 'YTickLabel', yticklabels2);

% 添加图例
legend([hTotalActual, hTotalFitted], {'Actual Daily Cases', 'Fitted Daily Cases'}, 'Location', 'northeast');

% 添加标注
text(datetime(2023,02,17), max(ylim) * 0.95, 'B', 'FontSize', 14, 'FontWeight', 'bold'); % 标注"B"

% 计算总体R²
AlignedYPredicted = interp1(timePredicted, yPredicted, timeActual, 'linear'); % 拟合值插值到实际时间点
TotalActualCases = sum(T{1125:1183, 2:7}, 2); % 总实际病例数
TotalPredictedCases = sum(AlignedYPredicted, 2); % 总拟合病例数
ResidualSumSquares = sum((TotalActualCases - TotalPredictedCases).^2); % 残差平方和
TotalSumSquares = sum((TotalActualCases - mean(TotalActualCases)).^2); % 总平方和
TotalRSquared = 1 - ResidualSumSquares / TotalSumSquares; % 总R²值

% 显示总体R²值
text(datetime(2023,04,05), max(ylim) * 0.85, sprintf('R^2 = %.2f', TotalRSquared), ...
    'FontSize', 12, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');

% 导出图像为高分辨率PNG文件
exportgraphics(fig, '6个年龄组拟合效果图.png', 'Resolution', 300);

%% 计算微分方程
function dxdt = computeDerivative2(aInitial, b, params, t, x)

    n = numel(x) / 7; % number of age groups
    idx = 1:n;

    S = x(0*n + idx); % 易感者
    E = x(1*n + idx); % 潜伏者
    I = x(2*n + idx); % 感染者
    A = x(3*n + idx); % 无症状感染者
    R = x(4*n + idx); % 康复者
    Vs = x(5*n + idx); % 接种未抗体者
    Vr = x(6*n + idx); % 接种抗体者

    N = S + E + I + A + R + Vs + Vr; % 总人口
    Betat = aInitial .* max(0, sin(2 * pi / 68* t + b));  % β矩阵添加sin函数
    
    % 计算状态变化率
    interactionTerm = Betat' * (I + params.k * A)./N;  % 交互项
    dSdt = params.br * N - params.dr * S - params.fai .* S - interactionTerm .* S ;
    dEdt = interactionTerm .* S + params.q .* interactionTerm .* Vs - params.p * params.omega1 * E ...
           - (1 - params.p) * params.omega * E - params.dr * E;
    dIdt = (1 - params.p) * params.omega * E - params.dr * I - params.gamma * I;
    dAdt = params.p * params.omega1 * E - params.dr * A - params.gamma1 * A;
    dRdt = params.gamma1 * A + params.gamma * I - params.dr * R;
    dVsdt = params.fai .* S - params.dr * Vs - params.f * Vs - params.q .* interactionTerm .* Vs;
    dVrdt = params.f * Vs - params.dr * Vr;
    
    dxdt = [dSdt; dEdt; dIdt; dAdt; dRdt; dVsdt; dVrdt];
end

%% 计算残差
function residue = computeParameterLikelihood2(aInitial, b, params, observedIncidence)

n = size(aInitial, 1);
idx = 1:n;
% 微分方程求解
dxdt = @(t,x) computeDerivative2(aInitial, b, params, t, x);
[tt, xPredicted] = ode45(dxdt, params.tSpan, params.initialStates);

% 插值对齐预测值
yPredicted = interp1(tt, (1 - params.p) * params.omega * xPredicted(:, 1*n + idx), ...
    ((0.5:params.tSpan(2))), "linear", 0); % omega * E

% 计算残差
yObserved = observedIncidence;
residue = norm(yObserved - yPredicted, 'fro');
end

4.3 结果

运行以上代码后，可以观察到5个年龄组的拟合效果都不错，但第6个年龄组因为病例数较少的原因R方较低，但总体（六个年龄组合并）的每日新发病例数拟合较高。 以下是模拟结果的可视化图表以及拟合的参数展示：

SEIAVR 模型拟合结果

图 2: 6个年龄组拟合效果图

SEIAVR 模型拟合参数

图 3: 拟合的参数值

5 扩展（保证仓室不变）

5.1 更换年龄组

如果给出的数据是3个年龄组的数据，或者数据在多个年龄组中某一个年龄组数据很小,年龄组也可以换成3个年龄组进行拟合，或者其它数量的年龄组，只要数据支持。

SEIAVR 模型拟合结果

图 3: 3个年龄组拟合效果图

5.2 模拟疫苗接种策略

通过拟合出来的( $\beta$ )可以很好的代表现实的传染病流行情况，通过修改模型中疫苗接种速率，或者调整初始( $V_r$ )仓室人口数量等情景。

$$\phi = \phi * 1.1 , \quad V_r = V_r + 0.4 \cdot N$$

注：以上扩展均在仓室设置 ( S, E, V_s, V_r, I, A, R ) 不变的前提下进行。
林宇豪
点外卖要用券，健身从入门到菜鸟，耳机仙人