题西林壁

宋·苏轼

横看成岭侧成峰，远近高低各不同。

不识庐山真面目，只缘身在此山中。

正式开讲之前，先放定场诗一首。

之所以要引用这首诗, 是因为它和我们今天的话题，虽然不能说一模一样，但可以说多少沾边。在临床试验的统计分析中，也会常常遇到这样的问题，“横看成岭侧成峰”，换个角度去观察和思考，结论也会截然不同。那么有没有跳出三界外，不在五行中，而从云端去俯瞰全局的方法呢？某种意义上来讲，生存分析中的MaxCombo方法，可说虽不中，亦不远矣。

在肿瘤试验，尤其是确证性试验中，受试者的生存获益是最受关注的。换言之，两种药的药效比较，通常就是两条生存曲线（即生存率随时间变化的曲线）的比较。如下图所示，在时间零点的时候，两组受试者的生存率都是100%，而后随着时间的推移，两组受试者的生存率都在逐渐的下降，而试验药（A药）组的生存率始终高于对照药组，领先的“比例“几乎固定（术语叫做“符合等比例风险模型”）。这是我们最想看到的结果，说明A药可以在整个试验观察期内，给患者带来毫无争议的稳定的生存获益。

这里的A药由于表现太过优异，有时候难免会成为别人家的孩子。

Anyway，如果自己家孩子的实际表现没有那么的尽如人意，而是这样的（图B）：

或者这样的（图C）：

甚至这样的（图D）：

其中，图B是early/Diminishing effect, 也就是说疗效随时间逐渐消失，好比小学成绩领先，但初中被别人迎头赶上；图C是late/delayed effect, 即延迟出现的疗效，好比前期成绩平平无奇，后面突然发力，一举考上985；图D为效应交叉，包括一些更为复杂的情况。

那么，如果遇到像图B、C、D这样的情况（术语为“非等比例风险”），应该怎么办？实际上也有多种方法可以处理，而MaxCombo法就是最常用的“抢救”A药的方法之一。常规的生存分析采用log-rank检验来做组间比较，log-rank检验给所有的死亡事件相同的权重；而MaxCombo法则尝试给不同时刻发生的死亡事件以不同的权重，看什么样的权重可以最好的体现出两组之间的差异。

这个思路类似于，现在要比较两个孩子从小学到高中的总学习成绩，在不能保证其中一个孩子从头一直领先到尾的情况下，我们分别尝试给小学阶段更高的权重，给初中阶段更高的权重，给高中阶段更高的权重，以及给所有阶段相同的权重，看看哪一种权重的给法最能体现出两个孩子之间学习情况的差异。这样一来，无论我们的孩子是小学阶段格外优秀还是高中阶段格外优秀（假设其他阶段也没有落后太多），这个优势都能被很好的发掘出来。

如果出现图D那种大交叉的情况，最后鹿死谁手还未可知，可能要具体情况具体分析。这里面其实还涉及到多重检验的问题，不过暂且按下不表。由于这种方法更倾向于想方设法去“发现”差异，一般被认为是不够保守的，所以通常只能做敏感性分析，而不能做主要分析。

综上所述，MaxCombo法其实是一组权重赋值方法不同的weighted log-rank检验，而每个事件所得到的权重，与它们发生的时间点是相关的。计算权重的时候，通常采用Fleming-Harrington （简称FH）公式，此公式的具体形式如下：

其中，是t时刻（也就是事件发生的时刻）的生存率，对于早期事件来说，这个生存率会比较高（接近于1）；而对于晚期事件来说，这个生存率会比较低。p和q是我们可以任意指定的参数，只要非负即可。因此，我们可以通过调节p和q，达到给发生在不同时刻的事件以不同权重的目的。理论上可以做任意多组尝试，试的多了，可能总有一款适合您。但是一般来讲，针对以下4种具有代表性的情况进行尝试就足够了：

1. G (0,0)：p=0，q=0;

2. G (1,0)：p=1，q=0;

3. G (0,1)：p=0，q=1;

4. G (1,1)：p=1，q=1;

在这4种情况下，权重随事件发生的时间的变化大致如下：

其中，G (0,0) 相当于普通的log-rank检验，给所有事件相同的权重；G (1,0) 给早期事件（相当于小学阶段）更多的权重；G (0.1) 给晚期事件（相当于高中阶段）更多的权重；G (1.1) 给中期事件（相当于初中阶段）更多的权重。

分别采用以上四组不同的（p,q）参数设置，计算每个事件的权重，并进行weighted log-rank检验，得到相应的weighted log-rank统计量及其方差。为了避免命名上的混乱，我们先把后面会出现的各位“演员“朋友们列个表：

注意：虽然我们一般做log-rank卡方检验，但实际上log-rank也是有z值（也就是标准正态分布）形式的。z值可以通过简单的计算得出，即z=统计量/sqrt(方差)。如

的计算方法可以此类推。

Z1，z2，z3，z4 是我们的联合主演，因为它们共同构成了最终的MaxCombo统计量。MaxCombo顾名思义，就是从4个检验构成的combination中取Maximum，也就是取结果最好、最显著的，换句话说，z的绝对值最大的。所以最终的MaxCombo统计量（z统计量）可以表示为：

好，统计量有了，下一步就是P值的计算。Z1，z2，z3和z4单独看来，都可认为服从方差为1的正态分布，但由于它们都是用同一套数据同一个方法计算出来的，仅仅权重不同而已，它们相互之间不独立，而具有不可忽视的相关性。因此，我们需要把 Z = {Z1，z2，z3，z4}作为一个向量来整体考虑，这个向量服从多元正态分析。在零假设（即4个数均值都为零）的情况下，

为了避免混淆，对于零假设下Z1，z2，z3和z4，我们分别Z01，z02，z03和z04来表示。上式中，ρ12是Z1和z2之间的相关系数（等同于WLR(0,0)和WLR(1,0)之间的相关系数），ρ13是Z1和z3之间的相关系数，以此类推。最终的MaxCombo检验P值可根据上述多元正态分布计算而得，也就是在均值为0的这种假设条件下Z01，z02，z03和z04这四个数中，至少有一个数的绝对值大于由真实数据算出的zmax的概率。注意这种基于多元正态分布的方法其实已经自动的处理了多重检验的问题，我们不需要再为之烦恼了。

至此，P值几乎呼之欲出，唯一欠的东风就是，上面的那些相关系数ρ如何确定。这就需要4个WLR统计量的方差-协方差矩阵中，16位金牌配角的大力支持：

我们知道，对于两个变量x1和x2，相关系数：

所以只要把上面这个矩阵计算出来，带入相关系数公式，就万事大吉了。

好，理论解说到此结束，下面是如何用SAS实现这个方法。

我是如何实现的

MaxCombo这个算法其实可以很方便的用R来实现，比如调用nph这个包，里面有现成的函数。而对于喜欢用SAS来铁锅炖一切的朋友来说，编程实现也是不难的。笔者主要参考了这个程序：https://github.com/dreaknezevic/combo-wlr/blob/master/combo_wlr.sas，为了增加可读性，在它的基础上做了一些简化和整理，感兴趣的朋友可以去看原文。

要想用SAS编程实现，大致需要做这么几件事情：

1） 采用上面4种不同的权重，分别计算weighted log-rank统计量WLR(0,0)、WLR(1,0)、WLR(0,1)和WLR(1,1)，以及它们各自的方差，并用z=统计量/sqrt(方差)的公式，将WLR统计量转化成z值，得到Z1，z2，z3和z4。

2） 按照zmax=max（∣z1∣,∣ z2∣,∣z3∣,∣z4∣），算出MaxCombo检验的统计量zmax

3） 计算4个WLR统计量之间的协方差

4） 构造4个WLR统计量的方差-协方差矩阵，并把它转化为相关系数矩阵

5） 利用零假设下的多元正态分布，求最终的MaxCombo检验P值

这五步里面，第一步很简单，可以用proc lifetest直接实现；第二步更简单，直接计算即可；第四步和第五步看上去也是常规操作，感觉难度应该在第二步。要推算出这4个WLR统计量之间的协方差，对我辈学渣来讲，可能是有挑战的，但是好消息是，我们可以站在巨人的肩上。根据Karrison（2016），

也就是说，WLR(0,0)和WLR(1,1)之间的协方差，就等于WLR(0.5, 0.5)的方差。而WLR(0.5, 0.5)的方差是可以直接用proc lifetest来计算的。怎么样，惊不惊喜，意不意外？到此是不是有种所有问题都已迎刃而解，一步踏入逍遥天境的感觉？其实还是高兴的有点早，不过咱们先看这几步相应的代码吧。

先来做weighted log-rank检验，并计算方差-协方差矩阵：

在输出的all_cov数据集中，几个关键的变量是这样子的：

检验：FH公式的具体参数设置 (p,q)

Fleming：WLR统计量，WLR(0,0), WLR(1,0), ……

Variance: WLR统计量的方差

i,j: 循环数

z：weighted log-rank检验z值

abs_z: z值的绝对值

注意，我们其实只关心 (p=0, q=0), (p=0, q=1), (p=1, q=0)和(p=1,q=1)对应的z值（即红框框里的，同时也是满足i=j条件的）。对于其他12条记录，我们只使用variance这个变量（相当于这4个感兴趣的WLR统计量两两之间的协方差），而这12条记录里面的WLR统计量和z值，是我们所不关心，也会弃之不用的。

计算方差-协方差矩阵时，由于对称性，最适合的方法是只计算上三角矩阵，然后再折叠下来。这里为了编程简便，直接把16个数字都算出来了，偷懒的同时浪费了一点点算力：）

下面计算zmax统计量，并确定zmax统计量是由哪一组(p,q)组合产生的：

将相关变量捞出乘盘，哦不，应该说塞入宏变量待用。

到这里，准备工作已经差不多了，下面该开始计算P值了，P值应该怎么算呢？

众所周知，正态分布是没有累积概率公式的，更不用说多元正态分布了。不过，在SAS中，一次元和二次元正态分布是可以通过调用内置函数来计算累积概率的。但是，我们抱歉的通知，现在已经是四次元了，所以P值的计算只能用数值积分的方法来解决。

参考程序采用了一个简单粗暴的办法，就是从零假设的多元正态分布中大量的抽样，然后直接数有多少样本是满足条件的，这类样本在总样本中所占的比例，就是P值。[话说笔者也只会这一招，看别人也在用这招，我就放心了]。大家也许会担心如果抽样次数多了，运行起来会不会很慢。笔者实测结果还好，至少在抽样100000次的情况下，运行速度是可以接受的，延迟不是很明显。

抽样得到的矩阵是这个样子的：

每一行都是零假设下的一个独立的样本向量，COL1-COL4相当于上面提到的z01-z04 （取了绝对值以后）。下面统计在这100000个样本中，有多少个样本是满足条件的（即z01, z02, z03 和z04 这四个数的绝对值中，至少有一个大于zmax ），并由此计算P值。

至此P值得到，基本完工。参考程序中还做了另外一件事件，就是处理计算出来的相关系数矩阵非半正定的这种特殊情况。相关系数矩阵应该满足两个条件：（1）半正定（2）对角线为1。先对相关系数矩阵做特征分解，分解成的形式，其中Q为特征向量组成的矩阵，P为由特征值组成的对角矩阵。如果有任何特征值<0.0001, 则将其赋值为0.0001。这样把改造后的P乘回去之后，可能对角线为1的条件又不满足了。于是强令对角线为1，然后再检查特征值是否需要调整，如此循环往复，直至收敛，目的是在满足半正定和对角线为1这两个条件的情况下，找到与初始的相关系数矩阵最接近的一个矩阵。由于笔者在实际的模拟中，并没有遇到因为相关系数矩阵非半正定而算不下去的情况，所以这一趴就先忽略了。

最后，与R的结果做对比，检查计算是否正确。

注意event这个参数，1=事件，0=删失，与我们常规的编程习惯是相反的，所以需要调整一下数据。随手产生一套试验组中位生存9个月，对照组中位生存7个月，100人每组的数据，SAS算出来的Z1，z2，z3，z4以及相关系数矩阵与R完全一致，在抽样100000次的情况下，SAS算来的P值是0.0893, R算出来的P值是0.08912, 还是非常接近的。因为R在做数值积分的过程中也用到了某些近似，可以说SAS的这个抽样方法反而是可以把P值算到任何精确度的方法，只要抽样的次数足够多。

关于生存分析中的MaxCombo方法浅析及SAS代码示例分享就到这儿，各位看官再见，敬请期待下一期。