【关键词】 危险;队列研究 　　Monte-Carlo simulation estimation of rate ratio based on person year

　　ZHENG Ying-jie，JIANG Qing-wu，ZHANG Zhi-jie.

　　Department of Epidemiology，School of Public Health，Fudan University，Shanghai 200032，China 　　【Abstract】 Objective To estimate the distribution of person-year based rate ratio by using Monte-Carlo simulation.Methods According to the example given，rate ratio（RR）model was constructed and simulated with Latin Hypercube Sampling.The sampled frequency distribution was fitted and95%confidence interval of RR was calculated with several methods.Results The simulated RR frequency distribution was fitted good among Pearson5，Lognorm，Gamma and InvGauss distribution with Pear-son5best.95%confidence interval of RR was similar among the different methods，but both interval of simulated RR were more narrower.Conclusions The method of combining Monte-Carlo simulation with Latin Hypercube Sampling could be used to estimate the distribution of the rate ratio based on per-son-year.It was expected to extend to much more complicated parameter estimation. 　　【Key words】 Risk;Cohort studies 　　Poisson过程用于描述单位时间（面积、容积、里程等）内可连续、稳定地发生某事件的情形，由3个参数组成 ［1］ :暴露”时间”t、观察结局数α和单位时间平均发生结局数λ。已知其中的2个参数，即可以通过Gamma或Poisson模型来估计第3个参数（图1）。在队列研究中，确定暴露与结局（发病或死亡）之间的关系常常依赖基于人时的相对危险度（rate ratio，RR），RR是暴露组和非暴露组的人时结局密度之比。由于暴露组和非暴露组均为相应总体的一个随机样本，假定其结局的出现是相互独立的，则可通过Poisson过程来描述。因此，RR实质上是两个Gamma分布的比值。 　　图1 描述Poisson过程的3个参数之间的关系及其适用的统计模型（略） 　　Figure1 The relationship and applicable statistical models of three parameters in poisson process　　　　由于Gamma分布函数本身的复杂性，对基于2个Gamma分布比值的RR的精确估计非常困难，因此在实际工作中，常常用计算相对简单的近似法，如Wald法或Score法。本文应用Monte-Carlo模拟，结合拉丁超立方抽样技术（latin hypercube sampling，LHS） ［2］ 实现对RR的分布估计。

　　1 材料与方法 　　1.1 @ Risk软件简介 @ Risk Version4.5是Mi-crosoft Excel的一个插件，是基于Monte-Carlo模拟技术的分析工具，通过建模结合使用概率分布函数、统计量函数和Excel相应功能，实现定量风险分析和参数估计，支持最后的决策。 　　1.2 RR在 @ Risk中的建模和模拟 在队列研究中，最常见的数据整理见表1，观察人时和结局例数是最基本的两个数据。在 @ Risk软件中，用Gamma分布构建结局密度模型以及RR分布模型，经过Monte-Carlo模拟抽样，获得RR的频率分布或其相关统计量，并可进一步拟合RR的分布。应用拟合分布、统计量函数和常用方法（Wald法和Score法）进行RR的95%可信区间的估计和比较。 　　表1 某危险因素与某结局的队列研究中@Risk的模型构建 （略）

　　Table1 The model construction between some risk factor and the events during a cohort study in@Risk 　　抽样策略是Monte-Carlo模拟的核心。 @ Risk软件提供2种抽样技术:Monte Carlo抽样和LHS抽样。Monte-Carlo抽样是完全随机的，当低概率事件对结果产生主要影响时，聚集现象（较大概率的数值更多地被抽到）变得非常明显，需要更多的抽样次数才能解决。而LHS抽样方法是近年发展的一种更有效的抽样技术，采用分层无返回随机抽样，与Monte-Carlo抽样相比，能够以较少的抽样次数、更短的时间实现估计参数的收敛。同时，拉丁超立方抽样方法有助于解决低概率分布问题。 　　因此本研究中所有结果均采用在 @ Risk中采用LHS抽样技术、抽样次数为10000次时的结果，应用χ 2 检验进行分布的拟合优度检验。 　　2 结果 　　下面的例子引自Stampfer等 ［3］ 的绝经后妇女使用激素与发生冠心病危险的队列研究，随访32317名绝经后的妇女，其中激素使用组和不使用组分别观察54308.7和51477.5人年，分别发生30和60例冠心病。 　　2.1 结局密度概率分布 由于结局的出现事件可通过Poisson过程来描述，并且假定暴露组和非暴露组均为相应总体的一个随机样本，因此绝经后妇女 激素使用组和不使用组每10万人年冠心病平均发生数λ的真值应当分别落在其相应的概率分布范围内（图2）。 　　图2 每10万人年冠心病平均发生数λ的概率分布 （略）

　　Figure2 The probability distribution of mean number of coronary heart disease per100000person yearλ 　　2.2 RR频率分布、拟合分布及其拟合优度检验 以结局概率分布作为输入模型，通过LHS抽样和Monte-Carlo模拟，获得RR的频率分布，并采用12种概率分布进行拟合（表2、图3）。由表2可见，Pearson5、Lognorm、Gamma和InvGauss分布满足拟合要求，实际上这4种分布几乎完全重叠（图形未提供，从表3的95%可信区间界值也可看出），但以Pearson5分布拟合最佳，其χ 2 值最小，为54.04。 　　图3 RR的频数分布和拟合分布（略） 　　Figure3 The frequency and simulated distribution of RR 　　2.3 RR的区间估计及与常用方法的比较 由于获得了RR的拟合概率分布，RR的各种统计量如均值、标准差、方差、任意百分位值、峰度系数和偏度系数等等，均可获得。实际上，在无法获得理想的拟合概率分布时，通过由 @ Risk提供的统计量函数，在Monte-Carlo模拟的同时，也可获得上述统计量。模拟的RR值95%可信区间与Wald法或Score法结果大致相当，但其上限值和下限值均略小。 　　表2 RR频率分布及其拟合分布的拟合优度检验（略） 　　Table2 The goodness-of-fit test of the frequency and simulated distribution of RR 　　表3 RR的区间估计的Monte-Carlo模拟法和常用方法的比较 （略）

　　Table3 The comparison of confidential interval of RR by Monte-Carlo simulation method and other in-common-use methods 　　3 讨论 　　近年来，流行病学意识到显著性检验或New-man-Pearson假设检验的局限性，因此常常进行参数的区间估计。由于参数估计常常涉及多个统计过程，需经多步建模而获得的参数模型往往非常复杂或实际上难以通过数学推导获得，因此通常在假定满足一定条件的情形下，应用近似法等来计算，在实际工作中应用广泛。但单是进行统计过程的每步点值估计，往往会造成严重的后果，任何一个统计过程中点值估计的偏差，都会导致最后的参数估计出现较大的错误。因此，要避免上述问题，就应当对各个统计过程中所有可能的结果有一个整体的掌握，以概率分布来替代点估计来描述建模过程中所有情形的不确定性，通过Monte-Carlo模拟技术达到整个统计过程中所有变量的所有信息，即包括变量的所有可能取值和这些值的发生可能性的测定结果，实现参数的精确估计 ［1］ 。基于Excel软件的插件-@ Risk软件是实现对多个统计过程建立的参数模型进行分布估计的良好工具。 　　以往基于人时的RR区间估计可分为两类:大样本资料常常采用近似的Score法或Wald法，其前提是假定在一定条件下RR符合对数正态分布，但这两种方法的标准差计算有所差别，Wald法被认为是Score法的一种近似;而小样本方法则常常首先应用二项分布来估计全部病例中的暴露比例，通过暴露比例与RR的关系，来推算RR的区间，并使用 中位P值以获得更高精度的参数估计。由于RR系两个Gamma分布的比值，因此要估计其分布异常困难。本文应用 @ Risk软件对例子中的资料进行简单的建模和Monte-Carlo模拟，获得对RR的分布估计。结果表明有4个分布（Pearson5、Lognorm、Gamma和InvGauss）满足拟合要求，但以Pearson5分布拟合最佳。因此，应用该拟合分布进行后续的统计量估计变得异常简单。另外，由于数据本身的复杂多样性，在不能获得较为满意的拟合分布时，可通过 @ Risk软件中提供的统计量函数同样可实现参数的估计。从RR区间估计的结果可以看出，模拟估计的结果与常用方法计算结果大致相当，这与本次例子的大样本有关。 　　应用 @ Risk软件模拟了统计过程的真实抽样，这有许多的优点:①构建模型的过程简单，容易理解，规避了许多复杂的数学公式;②可获得丰富的参数信息，如频率分布、拟合分布、各种统计量及其图表等等，为进一步的模型建立奠定基础;③建模后的所有工作均由计算机完成;④与近似法不同，本文方法只是对建模过程中使用的每个统计模型有应用条件的要求，而没有其他方面的条件要求，如近似法的大样本要求等;⑤采用的LHS抽样方法是近年发展的一种更有效的抽样技术，采用分层无返回随机抽样，与完全随机的Monte-Carlo抽样相比，能够以较少的抽样次数、更短的时间实现参数估计的收敛;同时，LHS抽样方法有助于解决低概率分布问题 ［2］ 。总之，本文以基于人年的RR分布估计为例，应用Monte-Carlo模拟结合拉丁超立方抽样技术，以概率分布的形式真实地模拟了统计过程中多个模型的抽样过程，从而实现最终的参数分布估计。 @ Risk4.5软件提供了37种的参数和非参数概率模型，为实际工作中参数的分布估计奠定基础，实际上本文的思路可应用于许多更为复杂的参数估计。

　　【参考文献】 　　［1］ Vose D.Risk analysis:a quantitative guide［M］.2nd.West Sus-sex:John Wiley&Sons Ltd，2001. 　　［2］ Iman RL，Conover WJ.Small sample sensitivity analysis tech-niques for computer models，with an application to risk assessment ［J］.Communications in Statistics，1980，A9，1749-1842. 　　［3］ Stampfer MJ，Willett WC，Colditz GA，et al.A prospective study of postmenopausal hormones and coronary heart disease ［J］.N Engl J Med，1985，313（17）:1044-1049. 　　（葛阜阳校）