covariance matrix r语言_NONMEM中$COVARIANCE失败时的救济方法，谈一下参数uncertainty估算问题...

最近我的好友朱校博士在他的wechat公主号（名称：PKPD，里面内容非常丰富）写了一篇关于NONMEM $COV中R矩阵的文章，我顺手在文章中给他加了一些其他的估算方法，把其中其他方法的部分搬过来与大家讨论一下。

$COVARIANCE在NONMEM中意义很重要，当你完成建模后估算出参数值，这个参数值估的准不准肯定是modeler非常关心的部分，而这也是$COVARIANCE在NONMEM里的功能。其他关于$COVARIANCE的内容请移步朱博士的公主号查看，再次强烈推荐。

首先是NONMEM里$COV的稳定性，这点确实为业内人士所诟病，尤其是当OMEGA block的维度高于5的时候，高维诅咒简直惨不忍睹。关于这一点，NONMEM的开发者Robert Bauer博士也表示无可奈何，只能对于特定情况特定地分析。因此，我导师给我的一个课题就是关于如何避免$COV的不稳定性，我们给出了若干种不同的方法规避，希望对大家有所帮助。

1. 先看看NONMEM内部的救济措施。上文中朱校提到的IMP与SIR现在已经内置在NONMEM中了，这个方法避免了求高维函数的二阶导数矩阵，用蒙特卡洛法来求处参数关于OFV的uncertainty。我在此给出两个推荐的setting。

a. $ESTIMATION MAXEVALS=9000 METHOD=1 INTER NONINFETA=1

$ESTIMATION METHOD=IMPMAP MAPITER=0 ISAMPLE=5000 NITER=15 EONLY=1

$COVARIANCE UNCON。

第一行是主估算单元，负责给出参数的估计值；IMPMAP过程负责根据主估算单元的结果进一步得出“精确”的高维分布信息。其中ISAMPLE与NITER可以根据需要自行调整；而MAPITER表示Expectation Maximum算法中的最大化过程，在求covariance矩阵时可以省略节省时间。EONLY=1表示不改变主估算单元的结果，单纯计算其精确的后验分布。

b. $ESTIMATION MAXEVALS=9000 METHOD=1 INTER NONINFETA=1

$COVARIANCE UNCONDITIONAL SIRSAMPLE=1000 MATRIX=R SIRNITER=3 SIRDF=N

其中第一行依然是主估算单元，第二行不再是传统的$COV而是加持了SIR非参数法。其中SIRSAMPLE是sampling的取样数目，可以根据需要调整。SIRNITER是迭代次数，次数越多结果越稳定。（注意：NONMEM是不会进行resampling过程的，只会计算出各个样本的weight）SIRDF=N表示使用自由度为N的T分布进行抽样，对于分布的扭曲偏度较高的有效。

c. 这两种内部方法的稳定性很好，但是缺点是耗时较长。

2. BOOTSTRAP。这个方法过于有名所以没什么好说的，我个人认为是求解模型参数的uncertainty的“黄金方法”，稳定有效但略慢，PSN里同样内置的专门的函数，使用很方便。

3. Stochastic simulation estimation (SSE)法。顾名思义就是用原模型simulate一组新的数据，然后用原模型进行拟合计算一组新的参数估计值。具体使用时，可以将原模型的$COV去掉，只计算估计值。如果进行了1000次SSE过程，我们将得到1000组参数的估计值，然后只要进行最简单的统计就能得到$COV相同的结果，当然耗时同样很长。该方法与BOOTSTRAP不相上下，我个人的经验是对于复杂模型容易产生local minima的，BOOTSTRAP胜出，因为SSE容易得到一大串非最优的参数估计；而对于数据内subject少或者存在一些特殊受试者的数据，SSE胜出；总的来说还是BOOTSTRAP稍胜一筹。

4. OPTIMAL DESIGN法。这个方法暴风迅速（NONMEM测试中的$DESIGN通常在5秒内结束，快到怀疑人生），对于大部分稳定模型的参数uncertainty估计有着谜一般准确度（我测试了10个模型，8个成功），其基本原理来自Cramer-Rao不等式，见我的PPT。FIM就是费舍信息矩阵，它有个好处是可以不需要估算模型就可以得到，而且它的逆矩阵是VARIANCE-COVARAINCE三明治矩阵的下限，因此很显然，如果模型得到的是OFV最小的参数估计，那大概率它的FIM的逆矩阵就能用来估计$COV的结果。具体的原理阐述请参考文献[1]。如果碰到$COV失败的模型且计算的OMEGA block是非sigular的话，强烈推荐尝试这种方法，工具可以使用Andrew Hooker开发的多平台工具包POPED。（小声BB：Bauer已经将optimal design内置入NONMEM的7.5测试版，我已测试过，非常好用，全新的$DESIGN不久将会与我们见面。）

5. 快速公式法。你没有听错，真的有个简单公式可以用来快速检测参数的population variance的uncertainty，这也是我课题的一个小结果。公式的使用前提是参数的random effect（通常就是NONMEM code里的ETA）最好是正态分布，当然偏一点也没问题，附上PPT一张。对于variance，只要X和Y相同即可。我最常使用的形式是：

当你已经得到了某个参数的OMEGA估算值

，它的uncertainty就是如公式所计算的，可以作为快速诊断工具（Quick Analysis）。而那个effective number则是得到EBE之后，在phi文件中找到其对应的ETC，然后用它除以

后得到一个小于1的数（我的课题内命名为个体shrinkage），再将数据内的subject数目减去所有subject的个体shrinkage相加之和即可得到。。（更新一下）因为我有一个朋友刚刚使用了这个公式，所以想多解释一句，这个公式估算的属于uncertainty的下限，一般情况都是这个公式的估算值约等于$COV的估算结果，一旦出现公式估算小于$COV很多的情况，就是这个参数的分布也许有扭曲的问题或者这个参数的shrinkage较大，需要特别留意，因此这个公式也算个小诊断工具吧。这个公式非常有趣，我在课题内尝试过几个不同方向居然推导出了相同的公式，在此也感叹一下数学之美。

参考文献：

1. Nyberg J, Ueckert S, Strömberg EA, Hennig S, Karlsson MO, Hooker AC. PopED: an extended, parallelized, nonlinear mixed effects models optimal design tool. Computer methods and programs in biomedicine. 2012 Nov 1;108(2):789-805.