合理进行均值比较——具有组间异质性的增长重复测量模型.pdf

四川精神卫生 2023 年第 36 卷第 S1期http：/合理进行均值比较具有组间异质性的增长重复测量模型胡纯严1，胡良平1，2*（1.军事科学院研究生院，北京 100850；2.世界中医药学会联合会临床科研统计学专业委员会，北京 100029*通信作者：胡良平，E-mail：）【摘要】本文目的是介绍与具有组间异质性的增长重复测量模型有关的5个基本概念、计算方法、一个临床试验实例及其SAS实现。基本概念包括重复试验与重复测量、方差协方差结构、均衡性与异质性、先验分布与后验分布、混合模型与多水平模型。计算方法涉及抽样方法和先验分布。临床试验实例涉及具有异质性的两组增长重复测量数据。本文呈现了采用SAS处理实例中增长重复测量资料的全过程，包括呈现两组儿童口腔内一特定指标（简称距离）随年龄变化的增长图、构建两组儿童线性重复测量模型、构建一个既考虑到异质性又能估计性别效应的模型、直接获得与性别有关的截距和斜率的后验分布，以及基于后验估计值直接推断概率。对输出结果作出了解释，并讨论了改用多水平模型所产生的效果。【关键词】重复测量；方差协方差结构；先验分布；后验分布；抽样方法中图分类号：R195.1 文献标识码：A doi：10.11886/scjsws20230201004 Reasonably carry out mean value comparison：repeated measurement model of the growth with inter-group heterogeneityHu Chunyan1，Hu Liangping1，2*（1.Graduate School，Academy of Military Sciences PLA China，Beijing 100850，China；2.Specialty Committee of Clinical Scientific Research Statistics of World Federation of Chinese Medicine Societies，Beijing 100029，China*Corresponding author:Hu Liangping，E-mail：）【Abstract】The purpose of this article was to introduce 5 basic concepts，calculation methods，a clinical trial example and its SAS implementation related to the repeated measurement model of growth with intergroup heterogeneity.The basic concepts included the repeated experiments versus repeated measurements，variance-covariance structures，equilibrium versus heterogeneity，prior versus posterior distributions，mixed versus multilevel models.The computational methods involved the sampling methods and the prior distributions.The clinical trial example involved two groups of repeated measurement of growth with heterogeneity.This paper gave the whole process of using SAS software to analyze the growth repeated measurement data in the example，including presenting the growth graph of a specific index（short for distance）in the mouth of the two groups of children with age，constructing the linear repeated measurement model of the two groups of children，constructing a model that took into account heterogeneity while estimating gender effects，directly obtaining posterior distributions for gender-dependent intercepts and slopes，and directly inferring probabilities based on the posterior estimates.The output results were explained and the effects of switching to a multilevel model were discussed.【Keywords】Repeated measurement；Variance-covariance structure；Prior distribution；Posterior distribution；Sampling method重复测量设计定量资料在临床研究中出现的频率非常高1-2，其中，对于具有组间异质性的增长重复测量资料，统计分析方法比较复杂。本文扼要介绍了基于贝叶斯理论和哈密顿蒙特卡罗求解混合模型中参数的两个主要内容3-4，针对一个临床试验实例，介绍了如何以图形的方式直观呈现增长重复测量资料具有组间异质性的情况；构建了两种模型，即两组儿童线性重复测量模型和一个既考虑到异质性又能估计性别效应的模型；还介绍了两种方法，即直接获得与性别有关的截距和斜率的后验分布以及基于后验估计值直接推断概率。1 基本概念 1.1重复试验与重复测量所谓重复试验，就是在特定的试验条件下，对一组具有同质性的n个个体进行试验，得到某定量指标n个试验结果之间是彼此独立的。所谓重复测182023 年第 36 卷第 S1期四川精神卫生http：/量，就是在某些因素取m个不同水平的条件下对众多个体中的每一个进行m次观测或试验，得到某定量指标彼此有关联性的m个试验结果，它们之间联系的密切程度取决于多个方面，包括与重复测量有关的因素的个数、水平数、相邻观测的时间间隔以及定量观测指标的性质等。一般来说，相邻观测结果之间的联系较密切，相隔越远的观测结果之间的联系较疏远。1.2方差协方差结构在多个重复测量的时间点上观测到一个定量指标的取值，计算每个时间点上定量指标的方差，再计算出任何两个时间点上定量指标之间的协方差，将这些方差和协方差按时间顺序排列成一个方阵，就是方差协方差矩阵。假设只有4个时间点，则最简单的方差协方差结构TYPE=VC和最复杂的方差协方差结构TYPE=UN的形式分别见式（1）、式（2）5。2000020000002002（1）21213141212232423141324223434324（2）1.3均衡性与异质性所谓均衡性，就是两组或多组受试对象在多项观测指标上的分布相同或差异无统计学意义。所谓异质性，就是两组或多组受试对象在多个重复测量的时间点上定量指标的方差协方差结构不同或差异有统计学意义。1.4先验分布与后验分布所谓先验分布，就是事先知道或假设某随机变量（在贝叶斯统计理论中，认为总体参数也是随机变量）服从某种特定的分布。所谓后验分布，就是基于已知样本和先验分布等信息，代入贝叶斯公式计算，得到所研究问题中关于某特定参数的条件分布。后验分布集中了总体、样本和先验三种信息中关于该参数的一切信息，同时，它又排除了一切与该参数无关的信息6。1.5混合模型与多水平模型所谓混合模型，就是既包含固定效应又包含随机效应的回归模型。其中，固定效应指某自变量或因素的各水平对结果变量的影响是相对确定的，在各水平组之间的变异度无统计学意义；而随机效应指某自变量或因素的各水平对结果变量的影响不是相对确定的，在各水平组之间的变异度有统计学意义。固定效应和随机效应通常都是针对自变量或因素而言的。多水平模型也称为随机系数模型，通常是把回归模型中的截距和某自变量之前的回归系数（或称斜率）视为“随机变量”，通过对资料中某些因素进行分层，从而计算出截距和斜率的“固定效应与随机效应”的数量大小。如果用于分层的因素只有一个，此时的回归模型就叫做二水平模型（分层因素为二水平，个体为一水平）；以此类推，还有三水平、四水平模型7。2 计算方法 SAS/STAT的proc bglimm过程的计算方法主要包括抽样方法和先验分布两部分内容，现扼要介绍如下。2.1抽样方法2.1.1共轭抽样共轭先验分布是一种先验分布，其中先验分布和后验分布在同一个分布族中。例如，如果通过使用方差为2的正态似然对独立同分布的随机变量yi建模，见式（3）；关于的正态先验，见式（4）。yinormal(xi，2)（3）normal(0，0)（4）则的后验分布也是正态分布，它的协方差和平均值分别见式（5）、式（6）。=-10+12i=1Nxixi（5）=(-10+12i=1Nxixi)-1(-100+12i=1Nxiyi)（6）Proc bglimm 过程在正态和多变量正态的情况下使用共轭抽样算法，见表2。使用时，对固定效应参数以及协方差G和R进行抽样。表2Proc bglimm过程中的共轭抽样算法Table 2Conjugate sampling in proc bglimmFamilyNormalNormalMultivariate normalMultivariate normalParameterVariance 2Covariance PriorMultivariate normalInverse gammaMultivariate normalInverse wishart19四川精神卫生 2023 年第 36 卷第 S1期http：/2.1.2Gamerman抽样算法Gamerman 抽样算法以发明家 Dani Gamerman的名字命名，是梅特罗波利斯（Metropolis）算法的一个特例。在Metropolis算法中，建议分布源自迭代加权最小二乘算法的一次迭代8。加权最小二乘算法在迭代循环中运行，对于每次迭代，观测值的一组权重用于最小二乘拟合。通过对当前残差应用权重函数来构造权重，使用当前迭代的参数值来形成建议分布，从中生成建议的随机值。Gamerman算法通过对变换后的观测值迭代应用加权最小二乘算法，得到一般线性模型中的最大似然估计量和渐近方差。变换后的响应变量定义如下，见式（7）。y i()=i+(yi-i)g(i)（7）对应的权重定义式见式（8）。w-1i()=b(i)g(i)2（8）因篇幅所限，Gamerman算法的具体过程从略。2.1.3哈密顿蒙特卡罗抽样算法哈密顿蒙特卡罗（HMC）算法也称为混合蒙特卡罗算法，是Metropolis算法的一个版本，该算法使用梯度信息和辅助动量变量从后验分布中提取样本9。HMC算法仅适用于连续参数。HMC将目标密度函数转换为势能函数，并为每个模型参数添加一个辅助动量变量r。由此产生的联合密度具有以下形式，见式（9）。p(，r)p()exp(-12rr)（9）式（9）中，p（）是参数的后验值。HMC从（，r）的联合空间中提取，丢弃r，保留作为p（）的样本。该算法利用哈密顿动力学的思想来保持物理系统的总能量，其中，是势能函数的一部分，r 是动能（速度）的一部分。随着速度的变化，势能也随之变化，从而导致参数空间中的运动。2.2先验分布GLMM模型包含各种类型的参数：模型语句中指定的固定效应系数、随机语句中规定的随机效应系数以及协方差矩阵的参数，包括G侧和R侧协方差矩阵。2.2.1固定效应系数的先验在混合效应模型中，固定效应系数的平坦先验被认为是无信息性的。平坦先验为参数()1的所有可能值分配相等的可能性。这是proc bglimm过程中固定效应系数的默认先验。除了平坦先验外，方差非常大的正态先验也被认为是无信息性或弱信息性的。另一方面，可以使用信息先验使方差在正态分布中变小。如果想获得正态先验的特定平均值和协方差，可通过提供包含正态先验的平均值和协方差信息的SAS数据集来实现。2.2.2随机效应系数的先验假设随机效应系数具有正态先验分布，见式（10）。Normal(0，G)（10）式（10）中，=(1，2，n)T和G是块对角矩阵，每个块是Gi。对于特定对象i，随机效应系数i正态分布，平均值为零，协方差矩阵为Gi。当尺寸相同时，任何i和j的Gi=Gj。可以在random语句中使用type=选项来指定G矩阵的类型。因篇幅所限，有关“G 侧协方差的先验”“尺度参数的先验”和“R侧协方差的先验”的内容从略。3 实例与SAS实现 3.1问题与数据结构3.1.1一个临床试验问题及数据【例1】本实例说明了如何使用两组增长重复测量数据对协方差结构中的异质性进行建模。在正畸治疗中，从垂体中心到牙齿附近的翼上颌裂的距离（Distance，单位为mm）的变化非常重要。一项牙科研究纳入27名（11名女生和16名男生）儿童为研究对象，他们在8、10、12和14岁时被重复测量前述体积的“距离”指标的数值10。试验安排和数据见表1。表127名儿童在8、10、12和14岁时被重复测量“距离”的结果Table 1Results of repeated measurement of distance of 27 children at the age of 8，10，12 and 14Person123GenderFFFy121.021.020.5y220.021.524.0y321.524.024.5y423.025.526.0注：Person是27名儿童的编号；Gender代表性别；y1y4分别代表每名儿童在8、10、12和14岁时被重复测量所得到的“距离（Distance）”202023 年第 36 卷第 S1期四川精神卫生http：/3.1.2对数据结构的分析本例资料是具有一个重复测量因素（即年龄）的两因素（即性别与年龄）重复测量设计定量资料。27名儿童在 4个时间点（8、10、12和 14岁）被重复测量 Distance 值，试验分组因素为儿童的性别（Gender），重复测量因素为年龄（age，可以转变为4个时间点，用time表示），定量观测指标为Distance。3.1.3研究者希望研究的问题垂体中心到翼上颌裂的距离是否随年龄而变化？变化模式如何？男生和女生的变化模式是否不同？同一受试者不同年龄段的重复测量结果之间的相关性如何？男生和女生的相关性是否相似？3.1.4创建SAS数据集设所需要的SAS程序如下：data pr；input Person Gender$y1 y2 y3 y4；Distance=y1；Age=8；Time=1；output；Distance=y2；Age=10；Time=2；output；Distance=y3；Age=12；Time=3；output；Distance=y4；Age=14；Time=4；output；drop y1-y4；datalines；1 F 21.0 20.0 21.5 23.02 F 21.0 21.5 24.0 25.53 F 20.5 24.0 24.5 26.025 M 22.5 25.5 25.5 26.026 M 23.0 24.5 26.0 30.027 M 22.0 21.5 23.5 25.0；run；【变量说明】Person代表儿童编号；Gender代表性别；y1、y2、y3、y4 分别代表每位儿童在 8、10、12和14岁时被重复测量所得到的“距离（Distance）”；Distance代表从垂体中心到牙齿附近的翼上颌裂的距离；age代表年龄。3.2用SAS实现统计分析3.2.1呈现两组儿童随年龄变化的距离增长图设所需要的SAS程序如下：proc sgpanel data=pr noautolegend；panelby Gender/columns=2 sparse；scatter x=Age y=Distance/group=Personmarkerattrs=（color=black symbol=circlefilled）；series x=Age y=Distance/group=Personlineattrs=GraphDataDefault（color=blue pattern=1）；label age=Age in Years Distance=Distance（mm）；run；【SAS主要输出结果及解释】女生和男生随年龄变化的距离增长情况分别见图1、图2。图1显示了明显的迹象（对于起初具有较小距离的儿童来说，在后来的几个时间点上，他们的距离增加幅度均较小）和近似线性的增长。然而，男生和女生的增长率明显不同。3.2.2呈现两组儿童线性重复测量模型中协方差矩阵的异质性所需要的SAS程序非常复杂，此处从略。下面给出“呈现两组儿童线性重复测量模型中协方差矩阵的异质性”的输出结果，见图3。由图3可知，独立分析解释了数据中的异质性（不同的群体具有不同的协方差结构），但模型无法呈现性别效应。图1女生随年龄变化的距离增长图Figure 1Distance growth with age for girls图2男生随年龄变化的距离增长图Figure 2Distance growth with age for boys21四川精神卫生 2023 年第 36 卷第 S1期http：/3.2.3构建两组儿童线性重复测量模型设所需要的SAS程序如下：proc bglimm data=pr seed=475193 outpost=pr_out DIC；by Gender；class Person Time；model Distance=Age；random Time/subject=Person type=un r；run；【SAS程序说明】如果将random语句换成repeated语句，拟合效果会变差。【SAS 主要输出结果及解释】女生组模型中各参数的后验汇总和区间的计算结果见表 2。输出结果为女生组线性重复测量模型的参数估计结果（DIC=114.852），从各行最后两个数据可看出，所有参数与0之间的差异均有统计学意义。男生组模型中各参数的后验汇总和区间的计算结果见表3。输出结果为男生组线性重复测量模型的参数估计结果（DIC=249.984），从各行最后两个数据可看出，截距、年龄和尺度参数与0之间的差异均有统计学意义；而协方差参数中，仅4个时间点上的方差有统计学意义，任何两个时间点之间的协方差均无统计学意义。3.2.4构建一个既考虑到异质性又能估计性别效应的模型设所需要的SAS程序如下：proc bglimm data=pr seed=475193 outpost=pr_out；class Person Gender Time；model Distance=Age|Gender；repeated Time/type=un subject=Person group=Gender r；run；【SAS 程序说明】model 语句使用年龄、性别和年龄*性别（即交互作用）3 项。Group=Gender（在repeated 语句中）为每个性别建立单独的R侧协方差矩阵。R选项打印估计的协方差矩阵。【SAS主要输出结果及解释】两组模型中各参数的后验汇总和区间的计算结果见表4。男生截距的后验均值估计值为15.907，女生截距的后验均值估计值为15.907+1.477=17.384。同样，男生的斜率估计值为0.825，女生的斜率估计为0.8250.344=0.481。女生的起点大于男生的起点，但女生的增长率约为男生增长率的一半。图3基于独立分析的关于R协方差矩阵的后验密度比较Figure 3Comparison of posterior density of R-covariance matrix based on independent analysis表2女生组模型中后验汇总和区间的计算结果Table 2Calculation results of posterior summaries and intervals in the girl group model参数InterceptAgeScaleRandom UN（1，1）数目5 0005 0005 0005 000均值17.3950.4790.4963.928标准差0.9460.0840.2041.91795%HPD15.5580.3090.1791.28619.2780.6420.8877.811表3男生组模型中后验汇总和区间的计算结果Table 3Calculation results of posterior summaries and intervals in the boy group model参数InterceptAgeScaleRandom UN（1，1）数目5 0005 0005 0005 000均值16.2410.7941.8533.714标准差1.2450.1070.7441.88195%HPD13.6270.5630.5220.68818.5580.9853.2987.246表4两组模型中后验汇总和区间的计算结果Table 4Calculation results of posterior summaries and intervals in the two group model参数InterceptAgeGender FGender MAge*Gender FAge*Gender MResidual UN（1，1）Residual UN（4，4）组Gender FGender M数目5 0005 0005 00005 00005 0005 000均值15.9070.8251.477-0.344-4.2324.261标准差1.1410.0961.411-0.122-1.8221.55295%HPD区间13.7920.637-1.165-0.589-1.6451.99518.2701.0144.371-0.101-7.8027.461222023 年第 36 卷第 S1期四川精神卫生http：/3.2.5直接获得与性别有关的截距和斜率的后验分布设所需要的SAS程序如下：proc bglimm data=pr seed=475193 outpost=pr_out；class Person Gender Time；model Distance=Age|Gender；repeated Time/type=un subject=Person group=Gender r；estimate Boys Intercept Int 1 Gender 1 0；estimate Girls Intercept Int 1 Gender 0 1；estimate Boys Slope Age 1 Age*Gender 0 1；estimate Girls Slope Age 1 Age*Gender 1 0；run；【SAS主要输出结果及解释】在过程步程序中使用estimate语句直接获得与性别有关的截距和斜率的后验分布的计算结果，见表5。由表5可知，直接获得了女生与男生的截距、斜率的后验均值及其95%HPD的估计值。3.2.6基于后验估计值直接推断男生斜率大于女生斜率的概率设所需要的SAS程序如下：data prob；set pr_out；pDiff=boys_slope-girls_slope；prob=（pDiff 0）；keep pDiff prob；run；%sumint（data=prob，var=pDiff prob）【SAS主要输出结果及解释】基于后验估计值直接推断男生斜率大于女生斜率的概率的计算结果，见表6。男生的增长率大于女生的增长率。预期斜率差为0.344。男生的生长斜率大于女生的概率超过99%（约为99.7%）。3.3结论本例两组儿童Distance的重复测量增长资料之间具有非齐性，分别在不使用和使用 estimate语句的条件下，按性别分组对两组儿童的线性重复测量模型建模，可间接或直接获得男生和女生截距和斜率后验均值的估计值，即男生截距的后验均值估计值 为 15.907，女 生 截 距 的 后 验 均 值 估 计 值 为17.384。同样，男生的斜率估计值为0.825，女生的斜率估计值为0.481。因此，女生的起点大于男生的起点，但女生的增长率约为男生增长率的一半。并且，男生的生长斜率大于女生的概率超过99%，它们之间的预期斜率差为0.344。4 讨论与小结 4.1讨论由于本例男生与女生的资料之间不具有齐性，故需要按性别分组（即使用by语句）进行建模。若不按性别分组计算，分别构建以下3个模型，得到的结果和结论如下。模型1（一般混合模型5，11）：在“第3.2.3节”的SAS程序中不按性别分组（即删除by语句），可得到如下主要输出结果。pD（有效参数个数）=11.718，DIC=449.986，截距=17.471，95%HPD为（15.74819.182）。模型2（随机截距模型5，7）：在模型1中SAS程序的基础上，在“run语句”之前增加一个语句“random int/subject=Person type=un r；”，可得到如下主要输出结果。pD（有效参数个数）=16.162，DIC=382.912，截距=17.543，95%HPD为（15.85319.486）；截距的随机方差为3.298（有统计学意义）。模型3（随机截距-随机斜率模型5，7）：在模型1 中 SAS 程序的基础上，在“run 语句”之前增加一个语句“random int Age/subject=Person type=un r；”，可得到如下主要输出结果。pD（有效参数个数）=15.980，DIC=372.047，截距=16.470，95%HPD 为（14.78118.203）；截距的随机方差为 2.226（有统计学意义），斜率的随机方差为0.233（有统计学意义），截距与斜率的随机协方差为-0.153（无统计表5直接获得与性别有关的截距和斜率的后验分布的计算结果Table 5Calculation results of posterior distribution of intercept and slope directly related to gender标签Girls InterceptBoys InterceptGirls SlopeBoys Slope均值17.38415.9070.4810.825标准差0.8261.1410.0740.09695%HPD15.70813.7920.3410.63718.99618.2700.6301.014表6基于后验估计值直接推断男生斜率大于女生斜率的概率Table 6Calculation results of the probability of directly inferring the slope of boys is greater than that of girls based on the posterior estimates参数pDiffprob数目5 0005 000均值0.3440.997标准差0.1220.05395%HPD0.1011.0000.5891.00023四川精神卫生 2023 年第 36 卷第 S1期http：/学意义）。比较上述3个模型可知，模型3优于模型1和模型 2。已 知 模 型 3 的 截 距=16.470，95%HPD 为（14.78118.203）；而基于“第3.2.3节”计算的女生与男生的截距的平均值为（17.395+16.241）/2=16.818，该数值与未按性别分组但考虑了随机截距-随机斜率模型计算的结果16.470非常接近。由此可知，当所分析的资料在某重要影响因素（例如本例中的性别）的不同水平之间存在异质性时，有两种解决方案：一是按重要影响因素分层建模计算；二是直接构建多水平模型。从“模型误差项的自由度越大、模型的稳定性越好”这一统计学常识来判断，在前述两种方案中，以第二种方案为佳。4.2小结本文介绍了与重复测量及其统计模型有关的5个基本概念，介绍了基于贝叶斯理论和马尔科夫链蒙特卡罗方法求解统计模型中参数的两个重要内容；针对一个具有组间异质性的增长资料，采用3种不同方法构建统计模型，实现了间接和直接估计模型中固定效应参数和随机效应参数的目的，并对输出结果作出了解释，给出了明确的专业结论。参考文献1王芳，马元业，杨冬冰，等.舍曲林联合小剂量奥氮平对抑郁焦虑共病患者的疗效J.四川精神卫生，2021，34（2）：135-138.Wang F，Ma YY，Yang DB，et al.Efficacy of sertraline combined with low-dose olanzapine in patients with depression and anxiety comorbidityJ.Sichuan Mental Health，2021，34（2）：135-138.2任俊赏，蒲澜.高频重复经颅磁刺激联合舍曲林治疗儿童青少年抑郁发作的临床对照研究 J.四川精神卫生，2022，35（2）：153-156.Ren JS，Pu L.Clinical controlled study of high-frequency repeated transcranial magnetic stimulation combined with sertraline in the treatment of depression in children and adolescents J.Sichuan Mental Health，2022，35（2）：153-156.3刘金山，夏强.基于 MCMC 算法的贝叶斯统计方法 M.北京：科学出版社，2016：4-117.Liu JS，Xia Q.Bayesian statistical method based on MCMC algorithm M.Beijing：Science Press，2016：4-117.4康崇禄.蒙特卡罗方法理论和应用 M.北京：科学出版社，2015：86-149.Kang CL.Monte Carlo method theory and applicationM.Beijing：Science Press，2015：86-149.5SAS Institute Inc.SAS/STAT15.1 users guide M.Cary，NC：SAS Institute Inc，2018：1205-1306，6533-6728.6茆诗松.统计手册 M.北京：科学出版社，2003：351-370.Mao SS.Statistical manualM.Beijing：Science Press，2003：351-370.7Wang J，Xie H，Fisher JH.Multilevel models：application using SAS M.Beijing：Higher Education Press，2009：1-112.8Gamerman D.Sampling from the posterior distribution in generalized linear mixed models J.Stat Comput，1997，7（1）：57-68.9Neal RM.MCMC using Hamiltonian dynamics M.Boca Raton，FL：CRC Press，2011：113-161.10 Pothoff RF，Roy SN.A generalized multivariate analysis of variance model useful especially for growth curve problemsJ.Biometrika，1964，51（3-4）：313-326.11 Littell RC，Milliken GA，Stroup WW，et al.SAS system for mixed models M.Cary，NC：SAS Institute Inc，1996：31-170，253-302.（收稿日期：2023-02-01）（本文编辑：陈霞）24