一种三角函数建模季节性变化的方法——以多发性硬化症复发数据为例

摘要

将图示分析与三角回归相结合，是研究复杂周期性现象（如多发性硬化症（MS）复发发作时间）的一种强有力方法。该方法能够无偏倚地表征复发发作的季节性趋势，从而对季节变化、紫外线辐射（UVR）和纬度的影响得出新的推论。

摘要

本报告介绍了一种基于 Stata 的新型方法，应用三角回归模型分析了来自20个国家、横跨两个半球的46个临床中心长达55年的多发性硬化症复发数据。该方法成功的关键在于战略性地结合图示分析，以指导并验证统计回归建模过程。初步的图示分析对于建立关于季节性和纬度因素对复发概率影响的存在性及其结构形式的合理假设至关重要，并用于检验所构建模型的表现效果。随后，必须采用三角回归来量化这些关系，校正重要混杂因素，并评估这些关联的可信程度。将绘图技术与回归建模同步结合，使模型得以系统性优化，直至达到最佳拟合收敛，从而能够对季节和纬度在预测多发性硬化症复发时间中的独立影响作出新的推断。这些方法有望应用于其他已知或怀疑随季节和/或地理区域系统性变化的复杂疾病及流行病学现象。

引言

多发性硬化症（MS）最常见的类型是复发缓解型多发性硬化症（RRMS）。RRMS 的特征是神经功能出现 episodic（发作性）恶化，随后部分或完全恢复。全球范围内，无论南北半球，MS 的发病率和患病率均随着距离赤道越远而升高^1-3。然而，RRMS 患者复发事件的发生频率是否也随纬度变化，以及这种关联是否存在潜在的季节性变异，目前尚不明确。迄今为止，关于复发时间季节性的研究局限于单一临床中心，导致对复发时间季节趋势的推论仅限于单一地理区域，因而无法探讨更广泛的纬度影响^4-14。这些研究还受到样本量小和复发数据稀疏的进一步限制。2000 年一项针对欧洲、美国和加拿大临床中心十项研究的荟萃分析显示，每项研究均至少纳入30例报告复发发作季节的病例，结果揭示了复发发作时间存在明显的季节性趋势，复发在春季达到高峰，冬季则处于低谷⁴。随后在日本¹⁵和西班牙¹⁶开展的较小规模研究也观察到了类似的年度周期性趋势。然而，一项类似的美国研究未能证实这一模式¹⁷。迄今为止，这些研究和观察均局限于北半球。MSBase 研究组最近分析了一个涵盖南北半球的全球性 MS 复发大数据集，旨在探讨复发发作时间的季节性趋势，以及纬度对复发概率峰值与季节性紫外线辐射（UVR）低谷之间关系的影响¹⁸。该方法的核心是应用三角回归分析，以可视化并评估复发发作时间及 UVR 分布的趋势。

本研究的总体目标是检验以下假设：在北半球和南半球，多发性硬化症（MS）复发发作时间的时序变化均随季节呈现可预测的规律，且这种季节性受到纬度的影响。采用三角函数建模方法来研究这些问题的理由在于，该方法具有灵活性，可用于刻画已知或疑似呈现离散、可预测且一致形态或模式的二维或三维现象，例如在具有季节性的生物学或流行病学现象中常见的年度峰谷周期。^19-22 传统时间序列分析方法（包括傅里叶分析）的一个缺点是其通常假定时间序列由随机过程所主导。^21,23,24 相比之下，将三角函数引入回归类模型的优势在于，既能促进对周期性数据中规律性和系统性结构的探索，又能利用回归模型的结构来分析其他相关因素，或对季节性的混杂因素进行调整。

三角回归法此前已在医学流行病学文献中得到广泛应用，用于探讨多种主题的时间规律，包括传染病暴发的检测、昼夜节律在自主神经系统功能障碍到早产胎盘早剥中的作用，以及先天性畸形的季节性相关因素和急诊就诊时间的分布等。^25-32 这类模型通常需要比传统时间序列分析更大的样本量，因此本研究是首次将其应用于全球多发性硬化症（MS）复发发作的数据集。本文所述的三角回归法是一种适用于研究任何已知或疑似随时间系统性周期变化现象的有效工具。此类建模不仅有助于描述和可视化这些模式，还可进一步帮助研究者探索这些趋势的潜在驱动因素及相关变量。

关于此处呈现的多发性硬化症（MS）复发发作的具体实例，使用散点图和残差图来可视化并评估所假设的三角函数模型形式与数据的拟合程度，是确定以下两点的关键步骤：1）观察到的数据是否提供了充分的证据支持复发发作时间存在季节性或其他时间趋势的假设；2）定义特定三角函数模型的正弦和余弦函数的频率与排列方式是否足够充分，从而允许使用该模型进行后续的推断与预测。回归建模还能够控制可能影响任何观察到的季节性或纬度效应的重要混杂因素，例如患者个体层面的复发倾向，尤其是那些本身随时间变化的因素，如复发前暴露于疾病修饰药物（DMD）治疗的持续时间。分离出MS复发发作时间的独立地理和时间预测因子及相关因素，有望指导对复发事件生物学机制的研究，进而为开发预防或延缓疾病加重的未来治疗干预措施提供依据。

MSBase 登记库

参与本次分析复发数据的多发性硬化症（MS）患者来自国际MSBase登记系统。该登记系统建立于2004年，通过基于互联网、由医生拥有并运营的系统，长期收集来自同意参与的MS门诊患者的流行病学、疾病活动性、临床检查以及相关检测特征和指标数据³³。各成员单位遵循统一的方案，规定了必须在约定的定期时间间隔内上传的最小数据集，以确保复发事件等结局数据能够被持续且前瞻性地汇总。复发起始日期被列为最小数据集中的强制性变量。此外，通常还会收集与这些复发事件相关的临床数据，包括皮质类固醇治疗情况以及受累的功能系统。采用统一的iMed数据录入系统进一步确保了各中心在数据收集和报告方面的一致性。本项目在每个参与中心均获得了人类研究伦理委员会的批准或豁免。所有纳入分析的患者均须根据当地法律签署知情同意书。

纳入标准

共纳入9811名患者，贡献了32,762次复发事件用于分析。截至该时间点，共有至少注册20例患者的临床多发性硬化中心同意参与，并将患者数据上传和追踪至登记系统中^st 2013年12月（数据汇总日期）之前符合条件的患者被纳入分析。为确保分析中包含的所有复发事件均为前瞻性观察所得，仅将首次记录患者残疾评估（采用Kurtzke扩展残疾状态量表（EDSS））之后发生的复发起始时间纳入分析。所有为分析提供复发数据的患者均符合多发性硬化症（MS）的正式诊断标准。^34,35

结局指标

本研究考虑了两个主要结局：1）在地理位置、半球和/或全球范围内，多发性硬化复发的起始概率是否存在时间上的变化；2）纬度与季节性紫外线辐射（UVR）低谷时间同随后的复发概率峰值日期之间的时间间隔（以月为单位）是否存在关联。MSBase研究组提出假设：由于远离赤道的地区绝对维生素D水平较低，且这些较远地区特定地点的人群维生素D水平在冬至后可能更早达到季节性最低点，因此低维生素D水平对多发性硬化复发概率升高的影响也应呈现出相应的时间和纬度分布模式。

复发定义及日期

复发被定义为出现新的症状或原有症状加重，且持续至少24小时，在无伴随疾病或发热的情况下，并与前一次发作相隔至少30天。该定义此前已在MSBase复发表型分析中应用。 ³⁶ 每位符合条件的患者的随访期，即可能观察到复发事件的时间段，定义为从首次EDSS评估日期起，至数据提取与汇总日期前该患者在登记系统中记录的最近一次EDSS评估日期为止的期间。当某一特定月份的复发确切起始日期无法获得或无法确定时，临床机构采用该月的1号作为复发日期^st 或 15^th 月份中的日期作为默认日期。在本报告分析的 32,762 例复发病例中，有 7913 例（24.2%）和 4594 例（14.0%）分别记录在 1^st 和 15^th 分别为该月各日中显著更高的比例，远高于该月其他各日记录到的0.8%至5.6%的范围。为校正此情况，将记录在1日或^st 15 中的^th 月份中的日期被随机分配至这两个默认日期前后15天区间内的某一天。通过敏感性分析验证了该方法的内部有效性，结果表明，在默认日期随机化条件下所建模型对复发高峰日期的估计值，与使用原始报告日期的模型或完全排除默认日期的模型相比，差异无统计学意义。

方案

注意：所述的每个步骤均对应所提供代码文件中具有相同编号的 Stata 代码部分。以下方案中的 Stata 命令名称已用斜体表示。

1. 准备并绘制观测到的复发起始数据

1. 通过单击打开一个 do 文件 "新建 Do 文件编辑器" 按钮并使用 生成 计算每个建模地理层级（地点、半球和全球）在十二个日历月中的复发起始次数的命令。通过单击执行该操作命令 "执行（操作）" do-file 中的 do-file 操作按钮
2. 使用 swilk 或sktest 使用Shapiro-Wilk检验或修正的Jarque-Bera检验，分别对汇总的复发计数数据或个体患者的未汇总复发数据进行正态性检验。选择代码并点击 "执行（操作）".
   注意：当存在显著偏态时，应对复发计数变量进行自然对数转换，并随后使用Shapiro-Wilk检验或修正的Jarque-Bera检验重新检验对数转换后变量的近似正态性。^38,39
3. 使用 生成 创建新变量的命令 "北半球月份" 为南半球日历月份增加+6的偏移量，以便在同一水平轴上按季节绘制北半球和南半球的复发情况。选择代码并点击 "执行（操作）".
   1. 绘制散点图，以每个半球的观察到的月度复发起始情况为数据点，y 轴表示复发频率，x 轴表示日历月份 双向散点图 命令。对每个位置重复此操作。通过查看图形查看器中自动打开的屏幕上的各个图表，观察一年中复发 onset 的峰值与低谷模式。
4. 使用 雷达 绘制按日历月份分布的复发频率雷达图的命令，每个雷达轴代表一个月份，按顺时针方向排列。选择代码并点击 "执行（操作）".
   1. 对所有位点重复此操作。通过查看图形查看器中自动弹出的每个图表，观察一年中复发 onset 的峰谷分布模式。
5. 运行 海洋 用于对观察到的复发数据应用Edward季节性检验的命令。^40-42 对所有地理层级重复上述步骤。

2. 模型构建与选择

1. 使用 generate 命令指定将在回归中使用的年度周期正弦和余弦三角函数。选择代码并点击 "执行 (do)"。
2. 使用 regress 命令指定基础模型的形式，以复发计数作为因变量结果，以步骤 2.1 中计算的正弦和余弦项作为主要解释变量。
   1. 将地点特异性的 UVR³⁷ 添加到基础模型中，作为额外的调整协变量，并使用分析权重 aweight 选项，根据每个地点贡献的患者数量对模型进行加权。选择代码并点击 "执行 (do)"。
      注：在自动弹出的结果窗口中记录模型的决定系数（R²）和残差误差。紫外线辐射：分析中包含的所有单独地点，其1979年至2004年（含）期间每个月的每日平均红斑加权环境UVR数据，均来源于美国国家航空航天局地球探测卫星总臭氧绘图光谱仪（NASA Earth Probe Total Ozone Mapping Spectrometer）。³⁷
3. 使用 predict 命令存储模型预测的月度 log(复发) 值。通过使用 generate 命令对 log(复发) 项取指数，将对数复发数转换回整数复发计数。选择代码并点击 "执行 (do)"。对所有地点重复此操作。
4. 使用 twoway scatter 命令，将步骤 2.3 中指数化后的预测月度复发估计值叠加到观测到的月度复发数据上。选择代码并点击 "执行 (do)"。
   1. 对所有地点重复此操作。在图形查看器中查看每个图表。
5. 使用 regress 命令，在步骤 2.2 中指定的模型基础上增加一对额外的谐波正弦/余弦项以扩展模型。选择代码并点击 "执行 (do)"。
   注：记录残差误差和决定系数。按照步骤 2.3 的方法保存并转换模型估计值，并按照步骤 2.4 的方法将模型估计值绘制在观测数据之上。对所有地点重复此操作。
6. 使用 regress 命令，在步骤 2.2 中指定的模型基础上进一步增加两对额外的谐波正弦/余弦项以扩展模型。选择代码并点击 "执行 (do)"。
   注：记录残差和决定系数。使用似然比检验将此模型与基础模型直接比较。使用 estat ic 事后估计命令生成赤池信息准则（Akaike Information Criterion）和贝叶斯信息准则（Bayesian Information Criterion）。按照步骤 2.3 的方法保存并转换模型估计值，并按照步骤 2.4 的方法将模型估计值绘制在观测数据之上。对所有地理层级重复此操作。

3. 估算复发高峰概率

1. 使用估计量的非线性组合函数（nlcom）以计算相位偏移的点估计值和95%置信区间，使用从步骤2.1至2.6中确定的最佳拟合模型。选择代码并单击 "执行（操作）".
   1. 将这些点估计值及相关的置信区间转换为表示复发频率峰值对应日历日期的数字（T_最大值）以及缓解期的复发频率（T_分钟），其中 1=1^st 一月和365=31^st 十二月和T_最大值 = 相移 + (365/4) 和 T_分钟 = 相位偏移 + ((365/4)×3)。对所有地理层级重复此操作。匹配 T_最大值 和 T_分钟 通过 Excel 查找文件转换为日历日期。
2. 使用 生成 计算峰谷差值（T）的命令_最大值 减T_分钟）针对每个地点，以每站点 100 名患者为标准进行标准化。采用 Wilcoxon 秩和检验，根据纬度范围比较标准化的峰谷差值。选择代码并点击 "执行（进行）".

4. 紫外辐射数据建模

1. 运行 使用 加载UVR数据的命令。利用该数据计算每个地点的月平均UVR。 egen 命令。选择代码并单击 "执行（进行）".
2. 绘制每个地点日历月份（x轴）与月均紫外线辐射（UVR）（y轴）的散点图，使用 双向散点图 功能。在图形查看器中查看每个图表，该查看器会自动打开并显示在屏幕上。
3. 对UVR数据重复步骤1.2，并使用 回归 指定一个基于地点层级的年际紫外线辐射（UVR）趋势的基础模型，其中将月度UVR设为因变量，并将步骤2.1中指定的正弦和余弦三角函数作为解释变量纳入模型。
4. 对UVR模型重复步骤2.4至2.6，且仅限于特定位置模型。这包括重新运行 双向散点图 在观测数据上叠加预测估计值的命令，并使用 回归 运行扩展谐波模型替代方案的命令。
5. 使用在步骤4.2至4.4中确定的最优的地点特异性月均紫外线辐射（UVR）模型，利用 生成 使用步骤3.1中指定的双角公式再次计算UVR的相位偏移点估计值及其相应的95%置信区间。计算T_分钟 （使用步骤 3.1 中所述公式计算每个地点的最低紫外线辐射值日期）。选择代码并点击 "执行（操作）".

5. 紫外线辐射从谷值到复发峰值滞后期的建模

1. 使用 merge 命令，将步骤 4.5 中模型估计的各地区季节性紫外线辐射（UVR）谷值日期与步骤 3.1 中的复发高峰日期进行合并。使用 generate 命令计算从 UVR 谷值日期到后续复发高峰日期之间的时间间隔（以月为单位）。选中代码并点击 "Execute (do)"。
2. 使用 sktest 命令，通过 Shapiro-Wilk 检验对 UVR 谷值至复发高峰的时间滞后变量进行正态性偏离显著性检验。选中代码并点击 "Execute (do)"。
3. 使用 merge 命令将各地区的纬度数据添加至数据集中。使用 abs(x) 函数将相对纬度转换为绝对纬度。选中代码并点击 "Execute (do)"。
4. 使用 regress 命令，通过拟合线性与二次回归模型并比较残差，检验时间滞后与绝对纬度之间关系的线性性。选中代码并点击 "Execute (do)"。
5. 使用 regress 命令，设定一个线性均值回归模型，以 UVR 谷值至复发高峰的时间滞后（月）为因变量，以每 10 度为单位的绝对纬度为预测变量。使用 aweights 回归选项，根据各地区贡献的患者数量对模型进行加权。选中代码并点击 "Execute (do)"。
6. 使用 twoway scatter 命令，将绝对纬度绘制在 y 轴，UVR 谷值至复发高峰的时间滞后（月）绘制在 x 轴。使用 lfit 图形选项叠加最佳拟合线，并使用 aweight 分析权重选项可视化各地区患者数量的相对权重。选中代码并点击 "Execute (do)"。

6. 患者层面复发倾向的敏感性分析

1. 使用 mepossion 命令指定一个混合效应泊松回归模型，其中月复发计数为因变量，第 2.1 步中定义的正弦和余弦三角函数再次作为固定变量纳入模型，基线 EDSS、MS 发病年龄以及既往暴露于 MS 特异性疾病修饰治疗被作为潜在混杂因素纳入，同时将唯一的患者标识符指定为随机效应。选择代码并点击 "Execute (do)"。
2. 重复第 2.4 至 2.6 步，以确定最优拟合的泊松模型。此过程包括重新运行 twoway scatter 命令，将预测估计值叠加到观测数据上，并使用 regress 命令运行扩展的谐波模型变体。
3. 使用估计量的非线性组合函数（nlcom）计算相位偏移的点估计值及其 95% 置信区间，并计算复发频率峰值的日期。将结果与主要分析进行比较。
4. 使用 generate 命令，按照第 5.1 步所述方法，基于第 6.3 步中得出的个体患者水平泊松模型对复发峰值日期的估计值，重新计算每个地区从紫外线辐射谷值到复发峰值之间的滞后月数。选择代码并点击 "Execute (do)"。
5. 使用 regress 命令，按照第 5.5 步所述方法，将绝对纬度重新建模为滞后的预测因子，并将结果与主要分析进行比较。选择代码并点击 "Execute (do)"。

结果

对来自20个国家46个临床中心的32,762例复发事件应用三角回归分析，为多发性硬化症（MS）复发发作时间在两个半球均呈现周期性和季节性这一观察结果提供了有力的统计学依据，并表明季节性紫外线辐射（UVR）低谷与随后的复发高峰之间的时间间隔与纬度相关。实现这一结论的关键在于依赖图示分析，以指导模型构建、评估和优化所必需的迭代过程。

按日历月份绘制的复发频率散点图分析显示，在所有地理层级上均存在年度周期，表现为春季高峰和秋季低谷。全球复发数据的雷达图进一步证实，北半球的复发高峰出现在5月（图1A）。当将南半球的复发起始数据与北半球数据合并后，春季高峰仍然存在（图1B），其中南半球地区在11月出现高峰。两个半球均记录到秋季低谷，北半球和南半球复发频率最低的月份分别为11月和5月。Edward检验进一步证实，复发起始时间显著偏离均匀、非季节性分布。综上所述，这些结果表明，在全部三个地理层级上观察到的多发性硬化（MS）复发起始时间的周期性变化，最符合一个单一的年度周期，即包含一个高峰和一个低谷，两者相隔约六个月。因此，选择包...

讨论

本文所述方案详细介绍了一种基于系统回归的技术，该技术通过可视化图表分析，对全球多发性硬化症（MS）复发发作数据进行研究。本方法以对来自南北两个半球共20个国家的复发数据进行相对简单的描述性分析为起点，使研究者能够探索关于MS复发发作时间规律的理论，并通过三角函数模型对这些理论进行正式检验。通过逐步绘制全球复发发作数据，并系统地绘制和评估观测数据的候选几何拟合曲线，研究发现紫外线辐射（UVR）季节性低谷与随后MS复发发作概率峰值之间存在依赖于纬度的相关性，这一发现在此前的MS流行病学研究中从未有过报道。此外，通过将趋势可视化与正式的统计建模相结合，该分析不仅证实了先前的荟萃分析结果，即季节性是北半球MS复发发作时间的影响因素，而且首次将这一观察结果扩展至南半球。

三角函数回归建模是一种灵活的工具，可用于正式探索具有周期性、时间或季节依赖性的周期现象，能够对符合几何形状趋势的数据进行统计表征，例如本报告中所分析的复发起始时间与紫外线辐射（UVR）数据均呈现出的年度周期性正弦曲线。然而，鉴于诸如复发起始时间等复杂、多因素流行病学趋势现象可能呈现多种形状和结构，有必要对原始数据以及观测数据...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Stata SE 版本 13	StataCorp, College Station, Texas	版本 13	用于分析的统计分析软件
Microsoft Excel 2010	Microsoft	2010	用于日历日期查询的电子表格程序