中医体质是指人体生命过程中，在先天禀赋和后天获得基础上所形成的形态结构、生理功能和心理状态方面综合的、相对稳定的固有特质。中医对于体质的论述起源可追溯至西汉时期的《黄帝内经》[1],其在《素问·四气调神大论》中以四时阴阳五行为中心来演绎，指明体质与四时节气、五运六气间存在密切关系，四时气候可影响人体体质而致病[2]。近代以来，以王琦教授为代表的研究团队对中医体质研究颇深，不仅提出奠定中医体质研究出发点和理论背景的9种体质学说[1],还表明体质具有可调性，可通过服用适宜的药食有效纠正偏颇体质，达到预防或延缓疾病的效果[3]。不仅我国中医自古以来研究体质，国外在历史上也对体质内涵有过研究，如《希波克拉底文集》中提及的四体液说、体型说、血型说、内分泌说和高级神经类型说等体质分类学说[4]。由此可见，无论是东方还是西方，自古以来人类对于体质变化的内涵探索就一直存在。

近些年以来，随着国家“健康中国2030”理念的提出，体病相关、体质健康管理等有关中医体质的研究如火如荼，面对此情形，为获得中医体质的长期走势，更好地帮助个人调理体质实现健康管理及中医“治未病”疾病预防，需要利用基于数学和统计理论的时间序列分析模型，分析和探索中医体质随节气变化数据，将其模型化进行拟合预测。利用模型构建预测未来中医体质类型变化，在一定程度上为个人体质饮食调理提供辅助参考，对国家“治未病”理念统筹发展具有一定实用意义。因此，采用时间序列分析方法，构建中医体质ARIMA预测模型，探索中医体质预测研究的新思路。

1、方法

1.1 数据来源

借助中医人工体质辨识量表，面向川西地区18～65岁健康或亚健康人群，采集该人群的中医体质辨识数据共1 920份，从中分别筛选经过预处理的平和质、气虚质、阳虚质、阴虚质、痰湿质、湿热质、血瘀质、气郁质、特禀质9种体质[1],在2020年6月夏至-2021年6月芒种二十四节气内共计1 574例中医体质辨识数据形成有效分析样本，如表1所示。

表1 二十四节气内体质类型变化数据分布 [n(%)]

1.2 数据量化

利用数字1～9分别量化平和质、气虚质、阳虚质、阴虚质、痰湿质、湿热质、血瘀质、气郁质、特禀质[5];根据采集日期完成节气的量化。

1.3 数据分析

1.3.1 中医体质预测与时间序列分析ARIMA模型

为更直观地观察体质类型随节气变化的数据特征，分别以节气和对应体质类型量化数值为横、纵坐标作时序图，详见图1。该图表明体质类型在2020年夏至节气-2021年芒种节气存在较为明显的季节性波动趋势。

由图1可以看出中医体质辨识类型随节气变化具有短期突变和长期趋变的特性，针对以上两种特性，本实验选择构建适用于该特性的ARIMA时间序列预测模型。ARIMA模型是20世纪70年代初由博客思(Box)和詹金斯(Jenkins)提出的一种著名的时间序列预测方法，又称差分整合移动平均自回归模型，基本思想是将预测对象随时间推移而形成的数据序列视为一个随机序列，若为非平稳序列，则将其处理转化为平稳序列，再做因变量与其滞后值和白噪声项当前值的回归，以此建立预测模型。该模型记为ARIMA(p,d,q),由ARMA(p,q)模型经过d阶差分序列平稳化转变之后得到，ARMA模型根据回归中所含部分的不同，包括自回归过程(AR)、移动平均过程(MA)、自回归移动平均过程(ARMA)。其中p为AR项参数，d为差分阶数，实现原始序列的平稳化转变，q为MA项参数[6,7]。

图1 2020年夏至-2021年芒种24节气内体质类型变化时序图   

在本实验中，体质类型变量的条件分布通常具有异方差和偏态性，为避免伪回归，缓解异方差，在不改变体质类型序列的性质及相关性前提下，通常会将原序列取自然对数，降低波动性，向平稳序列转化[8]。由表1可知，二十四节气内第一高发体质为量化值1的平和质，其在ARIMA建模过程中，经过取对数差分计算之后得到0,不具备继续研究的统计学意义，故此以下实证分析皆基于各节气内第二高发体质类型进行建模及预测，如表2所示。

表2 二十四节气与各节气第二高发体质类型对照表

首先构建ARMA模型确定p、q值，由于ARMA模型为自回归过程AR模型(即将时间序列看成是前期值和随机项的线性函数)和移动平均过程MA模型(即将时间序列看成是当期和前期的随机误差项的线性函数)的二者结合，所以设表2中的第二高发体质平稳化后的时间序列为yt,用以描述该平稳序列前期值、当期和前期的随机误差项二者关系的线性函数，包含p个自回归项和q个移动平均项，εt为第二高发体质类型序列的残差值(即残差序列),a为提高预测模型拟合度的常数，φp≠0,μq≠0,φi,μi为模型参变量，其中φi又称自回归系数，μi又称移动平均系数。所以首先根据式(1)确定p、q值：

yt=φ0+φ1yt-1+φ2yt-2+……+φpyt-p+a+εt-μ1εt-1-μ2εt-2-……-μqεt-q (1)

其次，利用体质类型变量的本期值与滞后值相减运算得出差分阶数，即确定d值(公式2),其中Δ为一阶差分算子，即d=1。

Δdyt=Δd-1yt-Δd-1yt-1 (2)

在公式1、2的基础上，获取合适的p,d,q值，建立起ARIMA(p,d,q)模型。

为保证ARIMA模型通过适应性检验，即εt通过白噪声检验，所以要满足εt的期望值为0,方差为σ2,对于任意的s≠t,εt和εs不相关，即E(εtεs)=0,详见式(3)。

E(εt)=0,Var(εt)=σ2,E(εtεs)=0,∀s≠t (3)

1.3.2 ARIMA模型建模步骤

ARIMA模型在考虑季节、趋势等不规则变动因素后对时间序列的分析预测基本流程如图2所示。

图2 ARIMA模型对时间序列分析预测流程图   

根据图2可知，ARIMA模型建模步骤如下：

(1)序列变换与数据平稳性检验。

可通过时序图观察及ADF单位根检验(即对第二高发体质类型序列的平稳性进行检验，若序列存在单位根，则该序列为非平稳序列，而ADF是一种常用的单位根检验方法，它的原假设为第二高发体质类型序列具有单位根，即非平稳，对于一个平稳体质类型序列数据，就需在给定的置信水平上显著拒绝原假设进行判断，若表2中体质类型序列呈现非平稳性，一般可以通过公式(2)中差分方法来消除单位根，实现第二高发体质类型序列的平稳化转变，只是需确定该序列平稳化转变过程所需的差分阶数d(若第二高发体质序列为平稳序列，即通过ADF检验，则跳过该步骤，进入第二步)。

(2)模型识别与定阶。

第一步差分处理完成后，通过绘制差分后平稳的体质类型序列的自相关图(ACF,描述两个不同节气体质之间的相关性，用以度量历史节气的体质类型对当下节气体质类型的影响)和偏自相关图(PACF,鉴于各节气体质间存在相互影响，所以在ACF基础上，剔除中间节气体质带来的干扰，单独研究这两个不同节气体质之间的相关程度),初步识别模型参数q和p。若该平稳序列的自相关图呈现拖尾(图像自q阶后，缓慢收敛至0附近)、偏自相关图呈现截尾(图像自p阶后，快速收敛至0附近)时，建立AR模型；若平稳序列偏自相关图呈现拖尾、自相关图呈现截尾时，则建立MA模型；若平稳序列的自相关和偏自相关图均呈现拖尾，再结合第1)步得出的差分阶数d,代入公式(1)、(2),则建立起ARIMA模型[6]。

(3)参数估计。

根据前两步得出的(p,d,q)参数值，计算模型显著性结果，AIC值和BIC值，调整参数p,q,同理依次计算上述三项结果，将参数调整前后的三项计算结果进行比较。

(4)模型预测及拟合。

依据第三步得出多种模型参数组合，绘制其模型预测拟合图，二者结合进行比较，依据AIC、BIC准则(数值越小，模型越优)和图中体质预测模型拟合程度的情形，确定预测体质类型的最佳模型参数(p,d,q)。

(5)适应性检验。

以公式(3)为基础，通过对平稳体质类型序列的残差值绘制自相关图(ACF)和偏自相关图(PACF),判断模型残差序列是否通过白噪声检验(即平稳体质类型序列的各残差值之间相互独立，没有任何相关性，为白噪声序列，又称纯随机序列，不存在分析价值)。

(6)误差分析及拟合优度判断。

检验该预测模型的体质预测值与体质拟合值二者之间的误差大小和体质预测模型的拟合性，至此，ARIMA体质预测模型建模完成[9]。

2、结果

2.1 平稳性检验

根据表2,分别以节气和体质类型量化值为横、纵坐标，作第二高发体质类型随节气变化图，如图3所示。图3可见随节气变化的第二高发体质类型序列(以下称为“second_constitution序列”)属于非平稳序列。

图3 第二高发体质类型随节气变化图   

本文使用ADF单位根检验法对second_constitution序列平稳性做进一步判断，该方法通过在second_constitution序列随节气变化的回归方程右边加上序列的滞后差分项ut来控制高阶序列相关。针对图3中存在明显范围波动的second_constitution序列，选择含有常数项和时间趋势项的ADF检验公式：

Δyt = γyt-1 + a + δt +∑pi=1βiΔyt-i + ut (4)

公式(4)中，a为常数，δt为时间趋势项[9]。ADF检验结果如表3所示，表中t统计量值为-0.899,即落在10%的置信区间之外，second_constitution序列以较大的P值(78.82%)接受原假设，由此可以判断该序列为非平稳序列。

表3 second_constitution序列ADF检验

对于非平稳的second_constitution序列，本文采用取对数后进行差分处理的方法将其转化为平稳序列，其间通过取两次自然对数降低序列波动性，缓解异方差情况，通过逐项相减差分消除序列中的线性趋势。运用公式(2),先令d值取1,计算一阶差分结果，并对一阶差分结果进行ADF平稳性检验，效果若不理想则令d=2进行二阶差分，以此类推，观察动态结果及ADF计算结果，直至second_constitution序列实现平稳化转变。

图4 second_constitution序列差分处理前后对比图   

图4是second_constitution序列数据的一阶、二阶差分处理的前后对比图，图中二阶差分后的序列(以下称为“second_constitution_d2序列”)在数值0的上下波动较为稳定，符合平稳序列的特点，直观上可以将该second_constitution_d2序列初步判断为平稳序列。对该序列做进一步ADF检验，结果如表4所示，得到t统计量值为-22.185,小于-3.452,落在1%的置信区间内，其对应的概率小于0.05,因此可以将second_constitution_d2序列判断为平稳序列，即确定出差分次数d的值为2。

表4 second_constitution_d2序列ADF检验

2.2 ARIMA模型识别与定阶

通过ADF检验过后的second_constitution_d2序列数据已符合ARIMA预测模型对数据的平稳性要求，之后通过绘制second_constitution_d2序列的自相关图(ACF)和偏自相关图(PACF)初步识别q、p的值(如图5所示)[6,9],用于描述second_constitution_d2序列当下节气体质与过去某一节气体质二者之间的相关性，并以当下节气体质为基础，指示预测未来某一节气的体质。图5中ACF自第一阶后逐渐衰减趋于0,呈拖尾状，即q值为1,PACF自第二阶后同样逐渐衰减趋于0,呈拖尾状，即p值为2,可初步判second_constitution_d2序列适合ARIMA(2,2,1)模型。

图5 second_constitution_d2序列的ACF图和PACF图  

2.3 ARIMA模型参数估计及预测拟合

重复拟合模型中参数p、q的可能取值，并计算相应参数对应的显著性Wald、AIC值、BIC值来初步判断模型的最佳参数，其中AIC值依据赤池信息准则(AIC)计算，BIC值则依据贝叶斯信息准则(BIC)进行计算[9,10]9-10],计算结果见表5。

表5 second_constitution_d2序列的显著性、AIC值、BIC值

观察表5中的计算结果，其中参数(2,2,1)对应的显著性Wald.值最大，参数(1,2,1)对应的AIC值、BIC值最小。依据AIC、BIC准则，初步判断ARIMA(1,2,1)模型更有利于对未来体质进行预测，但为进一步比较2个参数模型的实际拟合效果，对参数(2,2,1)、(1,2,1)对应的体质预测模型分别作模型拟合图，详见图6。

综合表5计算结果和图6模型预测拟合图结果，判定当p,q的取值均为1时，能够较好地拟合各节气时间的体质类型，与最初判断结果一致，所以将体质预测模型的最佳参数确定为ARIMA(1,2,1)。

2.4 模型适应性检验

确定模型参数之后，通过绘制模型残差序列的ACF图和PACF图来检验模型的适应性，检验结果(见图7)表明，模型残差序列不存在明显的自相关性，可判定其为白噪声序列[9],即平稳体质类型序列的各残差值之间相互独立，不存在任何相关性，属于纯随机序列，不再具有研究价值。所以ARIMA(1,2,1)模型的构建对second_constitution序列数据是适应的，即构建ARIMA(1,2,1)模型对预测未来体质是合适的。

图6 ARIMA模型各参数对应的模型拟合图   

图7 残差序列ACF、PACF检验结果   

2.5 模型误差分析及拟合优度判断

基于图7结果，分别对ARIMA模型(2,2,1)、(1,2,1)对应的体质类型预测结果进行误差分析，选用均方根误差(RMSE)、最大绝对误差(MaxAE)、平均绝对误差(MAE)和R-Squared4项指标评估模型预测拟合效果，前3项评估预测值误差，第4项评估拟合程度，计算结果如表6所示。

表6 模型各参数预测误差、拟合程度分析

由表6计算值可以得出，参数(1,2,1)模型的均方根误差(RMSE,度量体质类型序列真实值与其模型预测水平的相差程度，对离群体质数据较为敏感，可突出影响较大的误差值)、最大绝对误差(MaxAE,度量体质预测模型的最大误差，适用于想象体质预测的最坏结果)和平均绝对误差(MAE,同样衡量体质类型序列真实值与预测值的偏离程度，但不易受极端值影响，可准确反映实际体质预测误差的大小)三项误差分析结果均低于参数(2,2,1)模型的计算结果，故此得出参数(1,2,1)模型的预测值与拟合值之间的误差相对较小，预测的平均绝对误差仅为0.002,加之最后一项模型的R-Squared计算(衡量模型体质预测值与真实体质记录值的拟合程度，判断模型质量好坏),参数(1,2,1)优于参数(2,2,1),所以参数(1,2,1)模型的拟合程度更高。由此可见，ARIMA(1,2,1)模型能够较为准确地预测出体质的变化趋势，可为体质类型的未来预测提供依据。

3、结论与讨论

不同个体体质的存在和演变时间的不间断性体现了体质的连续性，而体质的特征在整个生命过程中，具有循着某种体质类型固有的发展演变规律缓慢演化的趋势，这就使得体质具有可预测性。中医体质通过中医体质量表(CCMQ)测量得出，能够较为科学地体现人群的生理、心理特征及发病倾向，具有低偏倚风险和良好的适用性，故可将中医体质作为预测变量[11]。本文从节气入手，以节气为时间序列进行中医体质预测研究，经过模型识别定阶、参数优选、适应性检验、误差分析等多步骤检验，最终构建ARIMA(1,2,1)预测模型的拟合图显示预测数据与真实记录数据较为接近，模型预测精度较高，拟合优度良好，可以较好地反映体质类型序列随节气的变化趋势。同时在一定程度上有望基于历史节气的中医体质辨识数据对短期未来体质类型变化进行动态追踪及预测，从而为个人体质健康提供饮食调理的科学依据。此外，体质是疾病的土壤，是预测疾病发生发展的重要因素，在CNKI上以“中医体质”“预测模型”为关键词进行检索[11],发现目前关于预测模型的研究多是与某种疾病相结合，进行体病相关预测研究，而对于中医体质本身进行建模预测分析的研究较少。加之鉴于《素问·四气调神大论》中表明四时气候可影响人体体质而致病[2],文章对四时节气高发体质类型进行统计分析后得出气虚质常发于秋冬之交(寒露、霜降、立冬、小雪)、湿热质常发于偏春季(大雪、大寒、雨水、春分、清明、立夏、小满)等结果。由此发现，在季节更迭之时，各体质类型之间有较高的倾向性发生传变。这也方便人们对照学习《黄帝内经》中的摄生法则，适应四时气候变化，为做好疾病预防提供了可能[2]。

所以本文将中医体质作为预测变量，构建ARIMA预测模型，是为预测中医体质未来的走向和趋势，提前把握各节气之间体质类型转变的可能趋向，以此提前做好饮食等方面的改善，调理体质，预防因四时气候变化而导致疾病的发生，实现“治未病”的效果。但由于实验中采用的样本统计量较少、数据采集地域范围较小，存在一定局限性，故计算得出的预测模型仅限于为未来开展类似研究提供些许参考及思路。后续或可针对目前已建立的模型采集不同节气、不同地域的大样本、高质量的中医体质辨识结果来进行外部验证，进一步优化、完善预测模型。总之，ARIMA(1,2,1)预测模型的构建为个人养生制定行之有效的饮食调理方案提供科学依据，通过调理偏颇体质使个体实现未病先防、欲病早治、既病防变[11]，对国家在“健康中国2030”战略下“治未病”理念的统筹发展具有一定实用意义[12]。

参考文献:

[1]王琦.中医体质学[M].北京:中国医药科技出版社,1995:5-9,35.

[2]李宝锋.中医人文关怀在老年桡骨骨折患者诊治中的效果分析[J].中医临床研究,2018,10(33):97-98.

[3]王琦.论中医体质研究的3个关键问题(下)[J].中医杂志,2006,47(5):329-332.

[4]马晓峰.中医体质学术发展史及中西医学体质学说比较研究[D].北京:北京中医药大学,2008.

[5]罗悦,林冰,温川飙.基于神经网络的中医体质与体检指标关联模型算法研究[J].时珍国医国药,2018,29(3):763-766.

[6]丁海峰,李立清.基于ARIMA模型的我国长三角地区卫生总费用时间序列预测分析研究[J].中国医疗管理科学,2022,12(2):4-10.

[7]董大勇,刘珂言.基于ARIMA模型的工业品出厂价格指数预测与分析[J].统计与决策,2016,32(1):179-181.

[8]王强.山西省FDI对就业效应的实证分析[J].商情,2019,31(13):86.

[9]葛娜,孙连英,赵平,等.基于ARIMA时间序列模型的销售量预测分析[J].北京联合大学学报,2018,32(4):27-33.

[10]李雅晴,谢平,桑燕芳,等.水文序列相依变异识别的RIC定阶准则——以自回归模型为例[J].水利学报,2019,50(6):721-731.

[11]吴涛涛,王琦,李岩,等.基于体病相关探讨临床预测模型的构建[J].中医杂志,2022,63(16):1507-1510,1531.

[12]方旖旎,王琦,张国辉,等.中医体质学在“治未病”中的应用研究[J].中医杂志,2020,61(7):581-585.

基金资助:国家自然科学基金(81904324);成都中医药大学大学生创新创业训练计划省级项目(S202110633091);

文章来源:张新格,罗悦.基于ARIMA时间序列模型的中医体质预测分析[J].亚太传统医药,2024,20(04):156-162.