关于数据挖掘:R语言基于copula的贝叶斯分层混合模型的诊断准确性研究附代码数据

最近咱们被客户要求撰写对于贝叶斯分层混合模型的钻研报告，包含一些图形和统计输入。

在对诊断测试准确性的零碎评估中，统计分析局部旨在预计测试的均匀（跨钻研）敏感性和特异性及其变异性以及其余测量。灵敏度和特异性之间往往存在负相关，这表明须要相干数据模型。因为用户，剖析在统计上具备挑战性

解决两个摘要统计，
必须思考敏感性和特异性之间的相关性，
必须思考到钻研中的敏感性和特异性的异质性
应该容许纳入协变量。

本教程介绍并演示了用于诊断准确性钻研的荟萃剖析的分层混合模型。在层次结构的第一级中，给定每个钻研的灵敏度和特异性，两个二项分布用于别离形容患病和衰弱个体中真阳性和真阳性样本数的变动。在第二级，咱们应用二元散布模仿未察看到的敏感性和特异性。尽管应用了分层模型，但 meta 剖析的重点在于钻研中的汇总平均值，而在给定的钻研估算中很少。

应用来自两个先前公布的 meta 剖析的数据集来演示这些办法：

尿液中端粒酶的诊断准确性作为诊断原发性膀胱癌的肿瘤标志物，因为它是一个有问题的数据集，其相干参数估计为 - 1 并且没有协变量而引起收敛问题（Glas et al.2003）
比拟病毒检测（应用 HC2 检测）的敏感性和特异性与反复细胞学查看对具备宫颈病变的女性进行分类，以检测潜在的宫颈癌前病变（Arbyn 等，2013）。第二个数据集用于证实具备一个协变量的元回归，该协变量能够天然地扩大到包含几个协变量。

因为其灵活性和 MCMC 模仿的应用，简单建模通常能够在贝叶斯框架内更容易地实现。通过管制先验散布，贝叶斯推断能够躲避可识别性问题，而没有先验散布的频率推理中的数值迫近算法可能会因识别性问题而陷入困境。然而，贝叶斯办法通常须要统计专业知识和急躁，因为 MCMC 模仿是计算密集型的。相同，最频繁的办法已被蕴含在规范“程序”中，这些程序须要较少的统计常识和编程技能。此外，频率论办法通过最大似然预计（MLE）进行优化，与 MCMC 模仿相比，其具备更短的运行工夫。

JAGS（Plummer 等人 2003）是 Stan 的另一种可扩大的通用采样引擎。扩大 JAGS 须要晓得 C++ 动态链接库（DLL）模块。依据教训，配置和构建模块是一项艰巨而繁琐的工作，尤其是在 Windows 操作系统中。上述毛病加上这样的事实，Stan即便从较差的初始值开始，往往会以较少的迭代收敛。

为了评估链的模型收敛和平稳性，有必要查看潜在的比例缩减因子，无效样本大小（ESS），MCMC 误差和参数的跟踪图。当所有链达到目标后验散布时，预计后验方差靠近于链方差，使得两者的比率靠近 1，表明链是稳固，可能已达到目标散布。无效的样本大小示意实际上无关某个参数的信息量。当样本主动相干时，冀望参数后验散布的信息少于样本独立时的信息。因为模仿了后验散布，因而近似值有可能偏离一定量;MCMC 误差靠近 0 示意可能已达到目标散布。

Watanabe-Alkaike 信息准则（WAIC）（Watanabe 2010）是一种最近的模型比拟工具，用于测量拟合模型在贝叶斯框架中的预测精度，用于比拟模型。WAIC 能够被视为对 Deviance Information Criterion（DIC）的改良，只管风行，但它曾经存在一些问题（Plummer 2008）。WAIC 是一个齐全贝叶斯工具，十分靠近贝叶斯穿插验证，对从新参数化不变，可用于简略以及分层和混合模型。

（Glas 等，2003）系统地回顾了细胞学和其余标志物（包含端粒酶）对膀胱癌初步诊断的敏感性和特异性。他们报告说，端粒酶有敏感性和特异性别离为 0.75 ，和0.86。他们得出结论，端粒酶不够灵活，不宜用于日常应用。

将数据加载到 R 环境中并生成以下输入

##    ID TP  TN FN FP
## 1   1 25  25  8  1
## 2   2 17  11  4  3
## 3   3 88  31 16 16
## 4   4 16  80 10  3
## 5   5 40 137 17  1
## 6   6 38  24  9  6
## 7   7 23  12 19  0
## 8   8 27  18  6  2
## 9   9 14  29  3  3
## 10 10 37   7  7 22

ID是钻研的标识符，DIS是患病的数量，TP是真阳性 NonDis 的数量，是衰弱 TN 的数量，是真阴性的数量。

（Arbyn 等人，2013 年）对人乳头瘤病毒检测的准确性进行了 Cochrane 评估，并反复细胞学剖析，对宫颈涂片进行查看以诊断宫颈癌前病变。他们 SAS 应用 BRMA 模型进行 METADAS 了 10 项钻研，其中应用了两种测试。他们用于 HC2 和反复细胞学的相对灵敏度别离为 0.909 [0.857, 0.944] 和0.715 [0.629, 0.788]。这些数据用于演示如何在回归设置中扩大截距模型。将数据加载到 R 环境中并生成以下输入


##    Test         StudyID  TP   FP  TN FN
## 1  RepC  Andersson 2005   6   14  28  4
## 2  RepC   Bergeron 2000   8   28  71  4
## 3  RepC Del Mistro 2010  20  191 483  7
## 4  RepC Kulasingam 2002  20   74 170  6
## 5  RepC     Lytwyn 2000   4   20  26  2
## 6  RepC      Manos 1999  48  324 570 15
## 7  RepC  Monsonego 2008  10   18 168 15
## 8  RepC      Morin 2001  14  126 214  5
## 9  RepC  Silverloo 2009  24   43 105 10
## 10 RepC    Solomon 2001 227 1132 914 40
## 11  HC2  Andersson 2005   6   17  25  4
## 12  HC2   Bergeron 2000  10   38  61  2
## 13  HC2 Del Mistro 2010  27  154 566  2
## 14  HC2 Kulasingam 2002  23  115 129  3
## 15  HC2     Lytwyn 2000   4   19  33  1
## 16  HC2      Manos 1999  58  326 582  7
## 17  HC2  Monsonego 2008  22  110  72  2
## 18  HC2      Morin 2001  17   88 253  2
## 19  HC2  Silverloo 2009  34   65  81  2
## 20  HC2    Solomon 2001 256 1050 984 11

Test是一个解释变量，显示分类测试的类型，StudyID是钻研标识符，TP是真阳性数量，TN是真阴性的数量，FN是假阴性的数量。

默认状况下，chains = 3 `cores = 3。从下面的代码中，从每个 3 链中抽取 1000 样本，抛弃第一个样本，而后使得每个链具备 900 个“预烧期”后抽取。种子值 seed = 3 指定随机数生成器以容许后果的再现性，cores = 3容许通过应用3` 核来并行处理链，每个链一个核。

上面的迹线图显示了链和收敛。

接下来，取得如下的模型概要预计


## 95% 置信区间的后验边缘平均值和中位数敏感性和特异性
##            Parameter      Mean      Lower    Median    Upper  n_eff  Rhat
## MUse[1]  Sensitivity  0.756762  6.904e-01  0.756036  0.81658 1196.9 1.001
## MUsp[1]  Specificity  0.798289  6.171e-01  0.813517  0.90640  704.4 1.004
## ktau[1]  Correlation -0.820176 -9.861e-01 -0.876343 -0.33334  269.2 1.015
## Varse[1]   Var(Sens)  0.006198  8.321e-06  0.005047  0.01947  165.7 1.007
## Varsp[1]   Var(Spec)  0.048111  1.357e-02  0.041060  0.12204  169.5 1.007
##
##
## 模型特色
##
## Copula function: gauss, sampling algorithm: NUTS(diag_e)
##
## Formula(1):  MUse ~ 1
## Formula(2):  MUsp ~ 1
## Formula(3):  Omega ~ 1
## 3 chain(s)each with iter=28000; warm-up=1000; thin=30.
## post-warmup draws per chain=900;total post-warmup draws=2700.
#### 模型的预测准确性
##
## Log point-wise predictive density (LPPD): -38.0607
## Effective number of parameters: 7.5807
## Watanabe-Akaike information Criterion (WAIC): 91.2828

从下面的输入，所述元剖析灵敏度 MUse[1] 和特异性 MUsp[1] 是0.7568 [0.6904, 0.8166]和 0.7983 [0.6171, 0.9064]。灵敏度和特异性别离为0.0062 [0, 0.0195] 和0.0048 [0.0136, 0.1220]。Kendall 在敏感性和特异性之间的 tau 相关性预计为-0.8202 [-0.9861, -0.3333]。

以下命令生成一系列森林图。

## $G1


## 
## $G2


## 
## $G3
## Warning: Removed 2 rows containing missing values (geom_errorbar).

$G1是钻研特异性敏感性和特异性（品红色点）及其相应的 95%置信区间（黑线）的图。$G2是后验钻研敏感性和特异性及其相应的 95%置信区间（黑线）的图。

$G3是后验钻研敏感性和特异性及其相应的 95%置信区间（黑线）的图。还给出了钻研特异性的灵敏度和特异性（品红点）及其相应的 95%置信区间（粗灰线）。

如上图所示，总体均匀敏感性和特异性存在“膨胀”：后验钻研的预计取决于全局预计，因而也取决于所有其余钻研。

接下来，通过创立如下列表来筹备数据

在 data 块中，指定了数据集中变量的维度和名称，此处 Ns 指数据集中的钻研数量。该 parameters 块引入了待预计的未知参数。etarho; 示意 Fisher 氏变换的关联参数的模式的标量，mul示意的灵敏度和特异性在对数下的平均值为核心的察看值，其中随机效应是矢量零。

在 transformed parameters 块中进一步转换参数。在 model 块中定义所有参数和数据似然的先验散布。最初，在 generated quantities 块中，loglik是计算 WAIC 所需的对数似然向量。

接下来，stan调用函数将代码转换为C++，编译代码并从后验散布中提取样本，如下所示提取参数估计，并应用以下代码进一步查看链收敛和自相干


## Inference for Stan model: 61572683b29d52354783115614fab729.
## 3 chains, each with iter=5000; warmup=1000; thin=10; 
## post-warmup draws per chain=400, total post-warmup draws=1200.
## 
##               mean se_mean     sd    2.5%     50%   97.5% n_eff   Rhat
## MU[1]       0.7525  0.0018 0.0517  0.6323  0.7562  0.8415   796 0.9999
## MU[2]       0.7908  0.0034 0.1095  0.5273  0.8094  0.9539  1045 1.0008
## mu[1]       0.7668  0.0013 0.0388  0.6869  0.7688  0.8369   891 0.9990
## mu[2]       0.8937  0.0027 0.0753  0.6943  0.9115  0.9825   789 0.9992
## rho        -0.9311  0.0070 0.1353 -0.9996 -0.9813 -0.5626   372 1.0077
## Sigma[1,1]  0.3376  0.0091 0.2918  0.0579  0.2554  0.9851  1026 1.0023
## Sigma[1,2] -1.2291  0.0272 0.8765 -3.4195 -1.0031 -0.2724  1040 0.9991
## Sigma[2,1] -1.2291  0.0272 0.8765 -3.4195 -1.0031 -0.2724  1040 0.9991
## Sigma[2,2]  5.6827  0.1282 4.1931  1.4720  4.6330 16.9031  1070 1.0002
## 
## Samples were drawn using NUTS(diag_e) at Mon Oct 09 09:19:55 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

所述元剖析灵敏度（MU[1]）和特异性（MU[2]）和 95%置信区间是 0.7525[0.6323, 0.8415] 和0.7908[0.5273, 0.9539]。这与作者以两种形式发表的文章（0.75 [0.66,0.74]和 0.86 [0.71,0.94]）不同。作者将规范双变量正态分布拟合到 logit 转换的敏感性和特异性值，在钻研中容许钻研之间的异质性，并疏忽了更高层次的分层模型。因而，作者必须应用 0.5 的连续性校对，这是分层模型中没有遇到的问题。

下图显示除了 Clayton copula 模型之外，大多数拟合模型的链条混合成果令人满意，简直没有自相干。

所有拟合散布预计的均匀灵敏度和特异性如下表所示。


## Warning: 6 (30.0%) p_waic estimates greater than 0.4.
## We recommend trying loo() instead.
##       Model   Parameter      Mean      Lower    Median      Upper    n_eff
## 1  Gaussian Sensitivity  0.756762  6.904e-01  0.756036  8.166e-01 1196.857
## 2  Gaussian Specificity  0.798289  6.171e-01  0.813517  9.064e-01  704.379
## 3  Gaussian Correlation -0.820176 -9.861e-01 -0.876343 -3.333e-01  269.179
## 4  Gaussian   Var(Sens)  0.006198  8.321e-06  0.005047  1.947e-02  165.705
## 5  Gaussian   Var(Spec)  0.048111  1.357e-02  0.041060  1.220e-01  169.508
## 6       C90 Sensitivity  0.751379  6.913e-01  0.753546  8.098e-01   25.638
## 7       C90 Specificity  0.807051  6.549e-01  0.821119  9.069e-01  119.897
## 8       C90 Correlation -0.528340 -9.766e-01 -0.725178 -4.020e-18    4.111
## 9       C90   Var(Sens)  0.004885  3.400e-04  0.003297  1.955e-02   11.615
## 10      C90   Var(Spec)  0.045694  1.556e-02  0.038049  1.020e-01  137.149
## 11     C270 Sensitivity  0.757528  6.877e-01  0.761163  8.210e-01  273.236
## 12     C270 Specificity  0.803502  6.328e-01  0.811740  9.097e-01 1081.987
## 13     C270 Correlation -0.697493 -9.827e-01 -0.808717 -3.332e-06   40.012
## 14     C270   Var(Sens)  0.006662  2.667e-04  0.005293  2.027e-02  556.055
## 15     C270   Var(Spec)  0.044767  1.268e-02  0.037922  1.112e-01 1098.815
## 16      FGM Sensitivity  0.759407  6.891e-01  0.761931  8.174e-01 2475.208
## 17      FGM Specificity  0.802588  6.453e-01  0.812498  9.045e-01 2293.332
## 18      FGM Correlation -0.174538 -2.222e-01 -0.222221  2.222e-01  785.016
## 19      FGM   Var(Sens)  0.005390  7.425e-07  0.004181  1.813e-02 1019.633
## 20      FGM   Var(Spec)  0.041890  1.177e-02  0.036671  9.997e-02 2479.371
## 21    Frank Sensitivity  0.756683  6.855e-01  0.758340  8.152e-01 2686.631
## 22    Frank Specificity  0.808239  6.472e-01  0.818777  9.110e-01 1910.561
## 23    Frank Correlation -0.706819 -8.550e-01 -0.692019  1.000e+00 2700.000
## 24    Frank   Var(Sens)  0.006678  5.896e-04  0.005280  2.140e-02 2699.766
## 25    Frank   Var(Spec)  0.042067  1.201e-02  0.035908  1.039e-01 1937.653
## 26     BRMA Sensitivity  0.752531  6.323e-01  0.756181  8.415e-01  796.037
## 27     BRMA Specificity  0.790796  5.273e-01  0.809420  9.539e-01 1044.902
## 28     BRMA Correlation -0.822353 -9.824e-01 -0.876654 -3.804e-01  238.268
## 29     BRMA  Var(lSens)  0.337556  5.792e-02  0.255387  9.851e-01 1025.609
## 30     BRMA  Var(lSpec)  5.682692  1.472e+00  4.632967  1.690e+01 1070.481
##      Rhat  WAIC
## 1  1.0014 91.28
## 2  1.0044 91.28
## 3  1.0154 91.28
## 4  1.0072 91.28
## 5  1.0069 91.28
## 6  1.1047 91.40
## 7  1.0304 91.40
## 8  1.3707 91.40
## 9  1.1005 91.40
## 10 1.0311 91.40
## 11 1.0096 90.75
## 12 1.0001 90.75
## 13 1.0407 90.75
## 14 1.0024 90.75
## 15 0.9999 90.75
## 16 0.9998 97.37
## 17 0.9996 97.37
## 18 1.0070 97.37
## 19 1.0034 97.37
## 20 0.9999 97.37
## 21 0.9994 90.55
## 22 0.9997 90.55
## 23    NaN 90.55
## 24 0.9990 90.55
## 25 0.9992 90.55
## 26 0.9999 86.76
## 27 1.0008 86.76
## 28 1.0218 86.76
## 29 1.0023 86.76
## 30 1.0002 86.76

后果以图形形式出现如下

下面显示，BRMA 模型和高斯 copula 双变量 β 预计的相关性更加极其。另一个极其是模型 FGM copula 双变量 β 的预计，这是因为 FGM copula 中关联参数的束缚，其中值在 | 2/9 | 内。

在上图中g1，只管相干构造存在差别，但五个双变量 β 散布的边际均匀灵敏度和特异性与 95%置信区间的轻微差别相当。

在没有预计艰难的状况下，上表显示了 Pearson 预计的相关性-0.8224[-0.9824, -0.3804]。这是因为贝叶斯办法不受样本量的影响，因而可能解决具备较少问题的小样本量的状况。

基本上，所有六个模型在第一级层次结构中是等价的，并且在指定“钻研特异性”敏感性和特异性的先验散布不同。因而，模型应具备雷同数量的参数，在这种状况下，比拟预测密度是有意义的。在查看时，来自五个基于 copula 的模型的对数预测密度实际上是等效的（min=-38.77, max=-37.89）然而参数的无效数量有点不同（min=7.25, max=9.92）。

ascus数据集有 Test 作为协变量。应用协变量是有意义的，以钻研其对敏感性和特异性（包含相关性）的联结散布的影响。以下将基于 copula 的二元 beta 二项分布拟合到 ascus 数据。

下图显示了实用于 ascus 数据的所有六个模型的迹线图，其中所有参数（包含相干参数（BRMA 除外））都被建模为协变量的函数。除了基于 Clayton copula 的双变量 β 的状况外，存在适当的链混合和收敛。

从基于 copula 的双变量 β 散布来看，很显著 HC2 和反复细胞学中的敏感性和特异性之间的相关性是不同的。


## Warning: 19 (47.5%) p_waic estimates greater than 0.4.
## We recommend trying loo() instead.
##       Model Test   Parameter     Mean   Lower     Median      Upper
## 9  Gaussian  HC2 Correlation -0.43812 -0.9984 -6.959e-01  9.847e-01
## 10 Gaussian Repc Correlation -0.91991 -0.9997 -9.643e-01 -6.103e-01
## 23      C90  HC2 Correlation -0.06588 -0.7610 -1.039e-17 -7.624e-19
## 24      C90 Repc Correlation -0.85157 -0.9804 -9.120e-01 -4.906e-01
## 37     C270  HC2 Correlation -0.03038 -0.6452 -7.000e-18 -1.782e-18
## 38     C270 Repc Correlation -0.77847 -0.9757 -7.058e-01 -5.394e-01
## 51      FGM  HC2 Correlation -0.07618 -0.2222 -2.215e-01  2.222e-01
## 52      FGM Repc Correlation -0.19819 -0.2222 -2.222e-01  1.894e-01
## 65    Frank  HC2 Correlation -0.48806 -0.8140 -4.497e-01  1.000e+00
## 66    Frank Repc Correlation -0.73784 -0.8627 -7.275e-01  1.000e+00
## 81     BRMA Both Correlation -0.84808 -0.9839 -8.980e-01 -4.497e-01
##       n_eff   Rhat  WAIC
## 9   154.238 1.0066 236.4
## 10   24.342 1.0542 236.4
## 23   30.089 1.0690 235.7
## 24    8.532 1.1128 235.7
## 37   76.410 1.0326 227.5
## 38    2.945 1.4613 227.5
## 51 2422.935 1.0007 245.1
## 52 2550.145 0.9997 245.1
## 65 2700.000    NaN 238.3
## 66 2700.000    NaN 238.3
## 81  102.983 1.0254 233.7

Clayton90模型具备最低的 WAIC。

因而，这个例子表明，查看模型的拟合度和合理性是否充沛是至关重要的，而不是自觉地依赖信息规范来抉择最适宜数据的规范。

从上面绘制的后验绝对敏感性和特异性，所有收敛的模型通常认为反复细胞学比 HC2 敏感性低，而特异性没有显着损失。

基于 Copula 的模型提供了极大的灵活性和易用性，但它们须要审慎应用。尽管本文中应用的 copula 具备吸引力，因为它们在数学上易于解决，但（Mikosch 2006）和（Genest and Remillard 2006）指出，从数据中估算 copula 可能很艰难。此外，copula 模型背地的概念略微简单一些，因而须要统计专业知识来了解和编程，因为它们尚未作为统计软件中的规范程序。

在本文中，简要探讨了几种用于诊断准确性钻研的 metat 统计模型。

在评估 meta 剖析的敏感性和特异性以及相关性时，模型之间存在一些差别。因而，有必要进一步钻研某些参数的影响，例如钻研数量，样本量和联结散布的指定对 meta 剖析的预计。

提出的贝叶斯模型应用 copula 来构建二元 β 散布，该模型预计特定钻研的敏感性和特异性，具备特定的随机效应值。

在 ASCUS 分类数据中，基于拟合模型的论断与作者得出的论断统一：HC2 比反复细胞学查看更敏感但更轻微，并且没有显著低于特异性巴氏涂片诊断宫颈癌前病变的女性。

尽管 BRMA 对于两个数据集都具备最低的 WAIC，但咱们依然倡议应用双变量 β 散布对灵敏度和特异性进行建模，因为它们能够间接提供 meta 剖析预计。

最受欢迎的见解

1.Python 中的 Apriori 关联算法 - 市场购物篮剖析

2.R 语言绘制生存曲线预计 | 生存剖析 | 如何 R 作生存曲线图

3. 用关联规定数据挖掘摸索药物配伍中的法则

4. 通过 Python 中的 Apriori 算法进行关联规定开掘

5. 用关联规定数据挖掘摸索药物配伍中的法则

6. 采纳 SPSS Modeler 的 Web 简单网络对所有腧穴进行剖析

7.R 语言如何在生存剖析与 COX 回归中计算 IDI，NRI 指标

8.R 语言如何找到患者数据中具备差别的指标？（PLS—DA 剖析）

9.R 语言中的生存剖析 Survival analysis 早期肺癌患者 4 例

关于数据挖掘:R语言基于copula的贝叶斯分层混合模型的诊断准确性研究附代码数据

原文链接：http://tecdat.cn/?p=3060

原文出处：拓端数据部落公众号

介绍

荟萃 Meta 剖析的统计办法

推理框架和软件

模型诊断

模型比拟和抉择

数据集

端粒酶数据

ASCUS 分类数据

截距模型

模型比拟

Meta 回归

探讨

论断

十分感谢您浏览本文，有任何问题请在上面留言！