在科学网—发表期刊学术文章-我作为作者的故事 - 谢钢的博文 (sciencenet.cn) 里我写道，Haas教授的微生物（带来的健康）风险的定量评估的统计分析框架QMRA的核心内容是所谓的贝塔-泊阿松混合分布剂量-反应模型（the beta-Poisson dose–response model）。这个模型是基于这样几个假设条件：（1）风险人员摄入了致病源微生物的实际数量/剂量是一个随机变量并服从一个已知平均剂量d的泊阿松分布；（2）任何一个致病源微生物只要到达了风险人员体内某个部位就能100%导致感染；（3）任何一个致病源微生物能够到达了风险人员体内某个部位的概率服从参数值为a, b的贝塔分布。贝塔-泊阿松混合分布剂量-反应模型可以在数学上表示为：P_I(d) = 1 – 1F1(a,a+b,-d), 这里P_I(d)代表在平均摄入剂量d时的被感染的概率，1F1(., ., .)是库玛合流超越几何函数（the Kummer confluent hypergeometric function）。因此，贝塔-泊阿松混合分布剂量-反应模型没有精确的解析解只有数值计算的近似解。为了方便没有该超越几何函数数值计算能力的研究人员也能最大限度地应用该模型，一个普遍使用的贝塔-泊阿松混合分布剂量-反应模型近似公式是P_I(d) = 1 – (1 + d/b)^(-a)，已有研究发现，要保证该近似公式的结果在合理范围之内必须满足两个条件：（1）a << b; (2) b >> 1。然而，究竟a应该比b大多少才算足够大，b值要大到多大满足要求，没有人给出具体的量化标准。我因为对这个在进行水微生物质量健康风险定量评估所普遍使用的贝塔-泊阿松混合分布剂量-反应模型近似公式的应用条件提出了具体的数值判断标准，因而在Haas教授做学科编辑(area editor in Microbial Risk Assessment)的Risk Analysis杂志上发表了相关的文章。

在2018年9月Haas教授邀请我做一篇关于针对轮状病毒(Rotavirus)的健康风险评估的新提出的一个剂量-反应模型的文章的同行评审人，虽然知道这完全是无偿自愿的工作，但我认为这是保证专业期刊学术水平的有效的做法之一，我接受了邀请。Risk Analysis不是一个‘开放获取’(Open Access)类的专业期刊，因此发表文章的作者是不必给所谓的‘出版费’的。而Risk Analysis的同行评审过程是典型的‘双盲评审’ -即仅仅根据杂志编辑所提供的信息投稿人与评审人都不可能知道对方是谁。让我没有料到的是针对这篇投稿文章我这个同行评审人从18年9月到20年3月连续做了三轮四次评审，而更令我感到奇怪的是最终我却没能在Risk Analysis上找到这篇文章发表的信息。

为了使我能比较具体地把这个评审过程报告给大家，我必须先解释一下所涉及的健康风险评估的剂量-反应模型的一些技术细节。Haas教授提出的的微生物污染物健康风险的定量评估的统计分析框架QMRA的最核心的内容就是贝塔-泊阿松混合分布剂量-反应模型：P_I(d) = 1 – 1F1(a,a+b,-d)。在实际应用时因为求解其中的库玛合流超越几何函数的困难，研究人员往往用近似公式模型P_I(d) = 1 – (1 + d/b)^(-a) 替代精确模型来分析样本数据；当我们把贝塔-泊阿松混合分布剂量-反应模型的第三个假设条件简化为‘任何一个致病源微生物能够到达了风险人员体内某个部位的概率是一个固定值r’，该剂量-反应模型则变成了更简单的‘指数剂量-反应模型’ （Exponential dose-response model）: P_I(d) = 1 – exp(-rd)。请注意到Haas教授在他的QMRA的书中也描述了一些纯粹从数学上按最佳拟合样本数据点的角度提出的其它几种健康风险剂量-反应模型，但这个贝塔-泊阿松混合分布剂量-反应模型的导出是以微生物导致的健康风险的物理机制过程为基础得到的-即前述的三项假设条件。在同样的能够拟合样本数据的前提下，一个以相关的物理机制为基础的数量分析模型当然是最好的。这是为什么Haas教授首先提出的贝塔-泊阿松混合分布剂量-反应模型能够经得起时间考验的根本原因。

我作为评审人的文章的作者（们）提出了一个新的剂量-反应模型P_I(d) = （1-ø）[1 – exp(-rd)]，其中ø代表风险人员自身具备的免疫能力（以百分比表示，0表示完全没有免疫力，1表示完全不会被感染）。文章的作者（们）用他们提出的新的剂量-反应模型对一组最常被引用的轮状病毒(Rotavirus)人体实验的样本数据进行了拟合分析并与标准的贝塔-泊阿松混合分布剂量-反应模型族的模型做了比较。结果表明他们的新模型的拟合度最好，所以他们认为这个研究结果值得在Risk Analysis上发表。

之前大家都用贝塔-泊阿松混合分布剂量-反应模型P_I(d) = 1 – 1F1(a,a+b,-d)来分析拟合反应模型对这一组最常被引用的轮状病毒(Rotavirus)人体实验的样本数据。这篇投稿文章提出了一个新的剂量-反应模型P_I(d) = （1-ø）[1 – exp(-rd)]其实是在原来的指数剂量-反应模型P_I(d) = 1 – exp(-rd)加了一个考虑风险人员的免疫能力的因子项（1-ø）。因为贝塔-泊阿松混合分布剂量-反应模型与新的剂量-反应模型同为二参数的模型（前者为a, b；后者为r, ø），新的模型的拟合度稍好一点，至少对这一组样本数据而言新模型应该被认为是一个改进（尽管程度有限）。但是，投稿的文章在介绍已经存在的贝塔-泊阿松混合分布剂量-反应模型的时候居然没有把Haas教授的QMRA的经典著作作为参考文献之一。不仅如此，文章作者还声称现有的剂量-反应模型没有把健康风险的形成机制考虑进来，而他们的模型则考虑了免疫因素是重大改进。显然，文章的作者在文献综述上做的完全不到位。我在评审意见中指出，他们的模型与Haas教授的模型都考虑了健康风险的形成机制，区别只是角度不同而已；从模型的角度就是参数化的方式不同而已。显然，文章作者根本没有能理解贝塔-泊阿松混合分布剂量-反应模型是怎么推导出来的。另外很明显，文章作者缺乏基本的统计学的基本训练。比如，把两个不同的概率分布函数（分别为二参数的分布），仅仅因为用了相同的两个希腊字母来表示其分布函数参数就认为这两对参数彼此对等(equivalent) -其英文原文为：“In this form, the approximate beta-Poisson shares the same parameters as the exact beta-Poisson.” (这是两组不同的参数，当然不是‘the same parameters同样的参数’！)；还有就是投稿文章的一些统计分析结果的图是明显的牛头不对马嘴的。基于上述理由，我的总的评审意见是在适当回应了我的要求修改的意见后可以发表。

2019年3月，编辑部把作者们的修改稿给我做第二次的评审，对修改稿我只提出了仍然存在的一个次要的数学方面的作者们显然没有弄明白的问题要他们作出明确及正确的表述（我在评审意见中把我认为的正确答案也提供给了他们）。编辑部的决定是：作者完成了小修改(minor revision)后即可发表。

奇怪的是，到了2019年9月，我居然又收到Haas教授邀请让我对这篇稿子评审，不过不是做为修改稿评审，而是做为作者们的同一篇文章的重新投稿来评审。在这一稿里，作者们用了统计假设检验p-值大于0.05来证明他们的模型的正确。读过我的科学网—统计显著性问题的历史由来及最新进展 - 谢钢的博文 (sciencenet.cn)的读者可以想象的到我的评审意见，我强烈要求作者们好好学习有关的文献，修改了统计假设检验的错误用法之后再来投稿。另外一点我指出要他们改正的是，在几个最重要的关于他们的剂量-反应模型是如何推导出来的公式的出处的问题上，必须要给出具体的参考文献（我甚至提醒他们，公式的出处是否就是Haas教授的QMRA书的某一页的哪个公式）；否则，作者们必须提供原始的证明。仅仅把重要公式罗列出来是不可以接受的做法。

时间来到了2020年的三月，我收到了Haas教授邀请让我对这篇稿子进行第三轮的评审，像第二轮一样是做为作者们的同一篇文章的重新投稿来评审。这一次，文章的作者们加进了用贝叶斯统计的方法(MCMC)来应用他们的‘新’的剂量-反应模型来分析样本数据。令我费解的是，作者们居然用R直接写MCMC的模型拟合及参数估计的程序代码，而不是用现成贝叶斯统计分析软件语言，比如在R里调用BUGS程序做贝叶斯统计分析。因为除非你是专门搞贝叶斯统计分析的又有很强的数学与计算机背景，自己写的MCMC的模型拟合及参数估计的程序代码基本上是质量不高的（比如模拟计算结果无法收敛或程序的缺陷未能排除，等等。）而且，他们的结论是，在模型表现的好坏比较时贝叶斯分析结果在某些模型上做不出来。无奈，我只好花了几个小时把他们所需要用到的有关的用于分析（参数估计与判断模型拟合程度）贝塔-泊阿松混合分布剂量-反应模型族的传统的频率法(frequentist approach)特殊的R函数代码全写了下来。在提供评审意见时附上这些R代码请编辑部转交供投稿文章作者们参考。我同时提醒他们，贝叶斯统计分析在他们的研究中完全是多此一举，因为他们用的是non-informative prior，这样算出来的结果(posterior distribution)理论上应该是等同于frequentist approach 的结果。如果他们真的明白贝叶斯分析的好处何在，则应该考虑用informative prior来导出posterior distribution。我特别提醒他们我在用贝叶斯统计方法构建贝塔-泊阿松混合分布剂量-反应模型的方面也做过一点研究，请他们参考我的研究成果（如本博文最后的截图所示）。

在2021年，我有空的时候耐心查阅了有关的文献资料，发现其实这篇投稿文章提出的剂量-反应模型P_I(d) = （1-ø）[1 – exp(-rd)]不是他们第一个提出的，作者们的‘创新’大概只是在于把它用在了轮状病毒的数据上并且指出了其拟合程度比贝塔-泊阿松混合分布剂量-反应模型的结果稍好。剂量-反应模型P_I(d) = （1-ø）[1 – exp(-rd)]之前就已经被提出了：Schmidt2014: Norovirus Dose–Response: Are Currently Available Data Informative Enough to Determine How Susceptible Humans Are to Infection from a Single Virus? Risk Analysis， Volume 35, Issue 7 p. 1364-1383, https://doi.org/10.1111/risa.12323。当然，这里还有一种可能就是Schmidt是这篇投稿文章的共同作者之一，因为Risk Analysis实行的是双盲评审，对此我无从知晓。有鉴于此，我干脆就把我写的R代码整理之后放在了我的researchgate的个人网页上供有需要的研究人员自由取用。

这三轮四次的投稿文章评审工作其实是耗费了我不少的时间精力，但我认为是值得的，不论是在提高自己的对相关问题的认识水平方面，还是对保证专业期刊发表的文章的质量上，我的这一点点付出都是值得的。

########################################################################

# Special R functions were developed (frequentist approach) for estimation and assessment of 

#     Dose-Response models used in risk analysis  (in particular in QMRA = Quantitative Microbial Risk Assessment)

########################################################################

# Everyone is welcome to use the following r code for any non-profit purposes.

# Please acknowledge the ownership/originality to the author: Gang Xie (John), March 2022, Office

of Research Services and Graduate Studies,

# Charles Sturt University, NSW, Australia wherever is necessary.

# All intellectual property rights are reserved (e.g., License CC BY 4.0).

# License CC BY 4.0 allows you, for free, to copy and redistribute the material in any medium or

format;

# or to remix, transform, and build upon the material for any purpose, even commercially as long as

the author's

# (or the institution that the author affiliated to) original intellectual property ownership is

acknowledged/recognised.

#

# Rotavirus data were used to demonstrate the usage of these functions.

#

# References:

# Haas2014: Haas, C.N., Rose, J.B., & Gerba, C.P. (2014). Quantitative microbial risk assessment, 2nd

edition, Hoboken, NJ: John Wiley & Sons.

# Schmidt2014: Norovirus Dose–Response: Are Currently Available Data Informative Enough to

Determine How Susceptible Humans Are to Infection

# from a Single Virus? Risk Analysis， Volume 35, Issue 7 p. 1364-1383,

https://doi.org/10.1111/risa.12323

######################################################################

用于分析（参数估计与判断模型拟合程度）贝塔-泊阿松混合分布剂量-反应模型族的特殊的R函数代码：根据不同的微生物污染物的样本数据可派生变化出六种专门的剂量-反应模型。有兴趣的读者可从我的Researchgate个人网页上下载详细完整的R代码：(PDF) R functions for QMRA Dose-Response models (researchgate.net)

下面这张截图是我在2015年底在黄金海岸（澳洲，昆士兰州）举办的一个有关贝叶斯统计的国际会议上发表的一张科学海报(Poster)的一部分。