このマークダウンではNRIとIDIについてまとめます。

はじめに

予測モデルの性能は，いくつかの尺度で評価することができる．

  • 総合指標
  • 識別能
  • 較正能

総合指標には、較正能と識別能が含まれており、例えば\(R^2\)はそのような指標である。
識別能とは、事象のある人とない人を分けるモデルの能力。ROC解析は識別力を評価する。
また、NRI上やIDIも識別力の指標となる。

較正能は、観測された結果とモデルの予測との間の一致度である。
キャリブレーションプロットとホスマー・レメショ適合度検定が較正能を評価する。

識別能はこれまではACUやC統計量によって定量化し評価されてきた。
AUCは、予測されるリスクのすべての可能なカットポイントにおけるモデルの感度と特異性を考慮したものである。

しかし、予測モデルは、しばしば、患者を診断や治療法の決定に対応するリスクカテゴリーに分類するために使用される。このため、絶対リスク推定値に基づく臨床的なリスクカテゴリーを適切に割り当てる能力に応じてモデルを比較するという考えが生まれた.

リスクモデルを更新することにより、個人は最初のモデルとは異なるカテゴリーに移動する可能性がある。リスクカテゴリーを変更となった（例えば低いから高いへ）個人の割合は、再分類割合と呼ばれる。再分類は、リスク予測モデルの改良が医療上の意思決定に与える影響の程度を評価するものである。 すべての再分類が正しいとは限らないため、正しい再分類の割合や再分類の改善割合などの代替指標が開発された。

NRI（Net reclassification improvement）は、従来の診断・予測に新しい因子を加える、もしくは別の診断・予測指標を用いたときに、再分類がどれくらいになるか、表す指標である。

NRIのような指標は非常に有効であり、マーカーや予測モデルの評価に関するガイダンス文書では、本格的な費用対効果分析に先立つステップとして採用されている。

NRI（Net reclassification improvement）

NRIには二種類ある。

Category-based NRI
- アウトカムの予測確率をあるカットオフ値でカテゴリーに区切る
- 意味のあるカテゴリー（例：ガイドラインで治療方針が変わるカットオフ値）を設定し，calibrationが良好なモデルなら臨床に直接適用できる情報が得られる
- 臨床に意味のないカットオフ値を用いた場合の解釈は困難 (Ann Intern Med.2014;160:122-131.)
- 意味のある閾値がない場合，イベント発生割合で区切る方法が提案されている(Stat Med. 2017;36:4455‒67)

category-free (continuous) NRI
- アウトカムの予測確率をカテゴリーを作らず連続的に扱う
- 使用を推奨される状況はほとんどない（イベントありで予測確率が0.1％でも上昇すれば正しく再分類されたことになる）
- category-free NRIの大きさを拡大解釈してはいけない（通常category-based NRIよりもか なり大きくなる）

IDI Integrated Discrimination Improvement

IDIはイベントあり、なしの人の平均予測確率の差である。

条件付き期待値は、イベントありとイベントなしの標本平均で推定される、つまり以下の式で表せる

\(IDI=(IS_{new} − IS_{old}) − (IP_{new} − IP_{old})\)
IS: Integrated Sn, IP: Integrated 1-Sp.

\(\hat{IDI}=(\tilde{\hat{p_{new,events}}} − \tilde{\hat{p_{old,events}}}) − (\tilde{\hat{p_{new,nonevents}}} − \tilde{\hat{p_{old,nonevents}}})\)

\(\hat{IDI}=(\tilde{\hat{p_{new,events}}} − \tilde{\hat{p_{new,nonevents}}}) − (\tilde{\hat{p_{old,events}}} − \tilde{\hat{p_{old,nonevents}}})\)

$={i in event} $
$={j in nonevent} $

C統計量(AUC)。平均化された感度 - カットポイント間で均等に重み付けしない
- リスク再分類を評価しない

IDI: 平均化された感度
- カットポイントを均等に重み付け
- 特異度の調整はC統計量と異なる
- 潜在的に重要な予測因子をより明確にすることができる

カテゴリフリーNRI: 予測確率において正確な移動をする被験者の割合
- イベントと非イベントNRIの性能はかなり異なる可能性がある

Loading data

データは、Hosmer and Lemeshow Applied Logistic Regression:第3版）より引用。
aplore3パッケージの中のlowbirthweightのデータを使用する.

library(aplore3)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.9
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## Fix cases 10 and 39 to make them identical to Dr. Orav's dataset
lbw <- as.tibble(lowbwt)
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
lbw[c(10,39),"BWT"] <- c(2655,3035)


## Recoding
lbw %>% mutate(racecat=race,
               low=case_when(low==">= 2500 g"~0,
                             TRUE~1),
               ftvcat=case_when(ftv=="None"~"None",
                                    ftv=="One"~"Normal",
                                    TRUE~"Many"),
               preterm=as.factor(case_when(ptl=="None"~0,
                                 TRUE~1)))->lbw

summary(lbw)
##        id             low              age             lwt           race   
##  Min.   :  4.0   Min.   :0.0000   Min.   :14.00   Min.   : 80.0   White:96  
##  1st Qu.: 68.0   1st Qu.:0.0000   1st Qu.:19.00   1st Qu.:110.0   Black:26  
##  Median :123.0   Median :0.0000   Median :23.00   Median :121.0   Other:67  
##  Mean   :121.1   Mean   :0.3122   Mean   :23.24   Mean   :129.8             
##  3rd Qu.:176.0   3rd Qu.:1.0000   3rd Qu.:26.00   3rd Qu.:140.0             
##  Max.   :226.0   Max.   :1.0000   Max.   :45.00   Max.   :250.0             
##                                                                             
##  smoke            ptl        ht        ui             ftv           bwt      
##  No :115   None     :159   No :177   No :161   None     :100   Min.   : 709  
##  Yes: 74   One      : 24   Yes: 12   Yes: 28   One      : 47   1st Qu.:2414  
##            Two, etc.:  6                       Two, etc.: 42   Median :2977  
##                                                                Mean   :2945  
##                                                                3rd Qu.:3475  
##                                                                Max.   :4990  
##                                                                              
##       BWT        racecat      ftvcat          preterm
##  Min.   :2655   White:96   Length:189         0:159  
##  1st Qu.:2750   Black:26   Class :character   1: 30  
##  Median :2845   Other:67   Mode  :character          
##  Mean   :2845                                        
##  3rd Qu.:2940                                        
##  Max.   :3035                                        
##  NA's   :187

ロジスティック回帰モデルの適合

アウトカムは新生児の低出生体重児の2値化指標。3つのモデルを作成。

null: 切片のみ lwt: 母体体重を予測因子とする lwt.racecat: 母体重量と母体人種を予測変数とするモデル

logistic.model.list <-
    list(null        = glm(low ~ 1,             data = lbw, family = binomial),
         lwt         = glm(low ~ lwt,           data = lbw, family = binomial),
         lwt.racecat = glm(low ~ lwt + racecat, data = lbw, family = binomial)
         )

予測値をデータセットに取り込む

lbw <- within(lbw, {
  ## Obtain fitted values from two models
  fitted.lwt             <- fitted(logistic.model.list[["lwt"]])
  fitted.lwt.racecat     <- fitted(logistic.model.list[["lwt.racecat"]])
  
  ## Changes in predicted probability
  change                 <- factor(sign(fitted.lwt.racecat - fitted.lwt),
                                   levels=c(-1,0,1),labels=c("Down","Unchanged","Up"))
  
  ## Mark classification
  fitted.lwt.pos         <- as.numeric(fitted.lwt         >= 0.4)
  fitted.lwt.racecat.pos <- as.numeric(fitted.lwt.racecat >= 0.4)
  
  ## Changes in classification
  reclass               <- factor(sign(fitted.lwt.racecat.pos - fitted.lwt.pos),
                                  levels = c(-1,0,1), labels = c("Down","Unchanged","Up"))
})

ROC analysis using pROC package

library(pROC)
## Type 'citation("pROC")' for a citation.
## 
##  次のパッケージを付け加えます: 'pROC'
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      cov, smooth, var
## Create two ROC curves
roc.lwt         <- roc(low ~ fitted.lwt, data = lbw)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc.lwt.racecat <- roc(low ~ fitted.lwt.racecat, data = lbw)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Plot ROCs
junk <- plot(roc.lwt,         lty = 1, col = 1, add = F, legacy.axes = T)
junk <- plot(roc.lwt.racecat, lty = 2, col = 2, add = T)
legend(0.4, 0.2, lty = 1:2, col = 1:2, legend = c("lwt","lwt.racecat"), bty = "n")

## Compare ROC DeLong’s test for two correlated ROC curves

対応しているのでDeLongで検定

## Test the increment
roc.test(roc.lwt, roc.lwt.racecat)
## 
##  DeLong's test for two correlated ROC curves
## 
## data:  roc.lwt and roc.lwt.racecat
## Z = -0.96694, p-value = 0.3336
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
##  -0.10359626  0.03514776
## sample estimates:
## AUC of roc1 AUC of roc2 
##   0.6131030   0.6473272

Best cutoff by closest to topleft method

## Best cutoff by closest to topleft method
coords(roc = roc.lwt , x= "best" , input = "threshold", best.method = "closest.topleft")

Discrimination boxplot and slopes using PredictABEL package

平行な2つの箱ひげ図の分離は、右の2予測因子モデルでより良くなっている。
判別の傾きは、非ケースとケースの平均を結ぶ直線の傾きである。

library(PredictABEL)
## Warning: パッケージ 'PredictABEL' はバージョン 4.1.3 の R の下で造られました
lbw<-as.data.frame(lbw) #Only data frame is readable

layout(matrix(1:2, ncol = 2))
slope.lwt <-
    plotDiscriminationBox(data = lbw, cOutcome = 2, plottitle = "lwt",
                          predrisk = fitted(logistic.model.list[["lwt"]]),
                          labels   = c("Without disease", "With disease"))

slope.lwt.racecat <-
    plotDiscriminationBox(data = lbw, cOutcome = 2, plottitle = "lwt.racecat",
                          predrisk = fitted(logistic.model.list[["lwt.racecat"]]),
                          labels   = c("Without disease", "With disease"))

c(slope.lwt = slope.lwt, slope.lwt.racecat = slope.lwt.racecat)
## $slope.lwt.Discrim_Slope
## [1] 0.031
## 
## $slope.lwt.racecat.Discrim_Slope
## [1] 0.058

判別の傾きは、有病者と無病者の平均予測確率の間で形成される傾きである。

2つの異なるモデルからの傾きの差がIDIである。

NRI and IDI using PredictABEL package

  • NRI（Net Reclassification Improvement）
  • IDI（Integrated Discrimination Improvement）

NRIは以下のように定義される。

\(NRI = [P(up|D = 1) - P(down|D = 1)] - [P(up|D = 0) - P(down|D = 0)]\)

カテゴリー別NRI: 与えられたカットオフを越えた再分類のみを考慮し、上方または下方への移動は連続NRIで考慮される。

これは事実上、良い再分類割合の合計から悪い再分類の合計を差し引いたものである。

カテゴリー別NRIでは、実証のためにカットオフ値を0.4に任意に設定した。

カットオフの設定は既存研究や、臨床的に決める。

reclassification(data = lbw, cOutcome = 2,
                 predrisk1 = fitted(logistic.model.list[["lwt"]]),
                 predrisk2 = fitted(logistic.model.list[["lwt.racecat"]]),
                 cutoff = c(0, 0.4, 1)
                 )
##  _________________________________________
##  
##      Reclassification table    
##  _________________________________________
## 
##  Outcome: absent 
##   
##              Updated Model
## Initial Model [0,0.4) [0.4,1]  % reclassified
##       [0,0.4)     104      16              13
##       [0.4,1]       5       5              50
## 
##  
##  Outcome: present 
##   
##              Updated Model
## Initial Model [0,0.4) [0.4,1]  % reclassified
##       [0,0.4)      34      16              32
##       [0.4,1]       3       6              33
## 
##  
##  Combined Data 
##   
##              Updated Model
## Initial Model [0,0.4) [0.4,1]  % reclassified
##       [0,0.4)     138      32              19
##       [0.4,1]       8      11              42
##  _________________________________________
## 
##  NRI(Categorical) [95% CI]: 0.1357 [ -0.0138 - 0.2853 ] ; p-value: 0.0753 
##  NRI(Continuous) [95% CI]: 0.3588 [ 0.0572 - 0.6604 ] ; p-value: 0.0197 
##  IDI [95% CI]: 0.0277 [ 0.0028 - 0.0526 ] ; p-value: 0.02938

手動で計算してみよう。
NRIは
\(NRI = [P(up|D = 1) - P(down|D = 1)] - [P(up|D = 0) - P(down|D = 0)]\)
なので、
Outcome presenceのUpdate modelでは16/(34+16+3+6),
Outcome presenceのInitial modelでは3/(34+16+3+6),
Outcome absentのUpdate modelでは16/(105+16+5+5),
Outcome absentのInitial modelでは5/(105+16+5+5),

よってカテゴリーRNIは以下で計算される。

## Categorical NRI manual calculation
(16 - 3) / (34 + 16 + 3 + 6) - (16 - 5) / (104 + 16 + 5 + 5)
## [1] 0.1357236

Functionにして計算するなら次のようになる。

## Implementation as a function
catNriElements <- with(lbw, by(fitted.lwt.racecat.pos - fitted.lwt.pos, low,
                                function(x) {
                                    ## Probability of up/down/unchanged
                                    propChanges <- prop.table(table(sign(x)))
                                    ## up - down
                                    propChanges["1"] - propChanges["-1"]
                                }))
catNriElements["1"] - catNriElements["0"]
##         1 
## 0.1357236

Continuous NRIでは以下。 Continuous NRI: \[[(up - down)/(up + down + unchanged) in low = 1] - [(up - down)/(up + down + unchanged) in low = 0]\]

## Continuous NRI: [(up - down)/(up + down + unchanged) in low = 1] - [(up - down)/(up + down + unchanged) in low = 0] in terms of predicted probability changes

xtabs(~ change +low, lbw)
##            low
## change       0  1
##   Down      74 23
##   Unchanged  0  0
##   Up        56 36

なので手動計算は

(36 - 23) / (23 + 0 + 36) - (56 - 74) / (74 + 0 + 56)
## [1] 0.3588005
## Implementation as a function
contNriElements <- with(lbw, by(fitted.lwt.racecat - fitted.lwt, low,
                                function(x) {
                                    ## Probability of up/down/unchanged
                                    propChanges <- prop.table(table(sign(x)))
                                    ## up - down
                                    propChanges["1"] - propChanges["-1"]
                                }))
contNriElements["1"] - contNriElements["0"]
##         1 
## 0.3588005

NRIの結果の書き方

詳細は論文を参照(Ann Intern Med. 2014;160:122-131)
- event, nonevent NRIを分けて報告する
- Reclassification tableはNRIよりも情報が多い
- 単位: event, nonevent NRIは％表記してもよいが，overall NRI(-2から2までとりうる）は単位がないので，％表記は不適切

  • Calibration: 比較するモデルのcalibrationに関する情報を提示する

結果の解釈

異なる対象者や異なるアウトカム，異なるカットオフ値を用いた場合のNRIを直接比較することに基づいて強い結論を導いてはならない NRIは予測ではなく，あくまで再分類を評価している Ann Intern Med. 2014;160:122-131.

IDI

IDIは予測確率の平均値の比較なので以下で求まる。

## IDI: Compare mean predicted probability among low = 1 and low = 0, and subtract
idiElements <- with(lbw, by(fitted.lwt.racecat - fitted.lwt, low, mean))
idiElements["1"] - idiElements["0"]
##         1 
## 0.0276851

これらはreclassificationで全て求まるので楽

Graphical representation of NRI

非ケースの場合（結果＝0）、予測確率が減少するのは良いことである。
ケース(結果=1)の場合、予測される確率を上がると良いと言える。 これを図示してみる

library(reshape2)
## 
##  次のパッケージを付け加えます: 'reshape2'
##  以下のオブジェクトは 'package:tidyr' からマスクされています: 
## 
##      smiths
lbw.melt <- melt(lbw[,c("id","low","change","reclass","fitted.lwt","fitted.lwt.racecat")],
                 id.vars = c("id","low","change","reclass"))

lbw.melt

Categorical NRI

ggplot(lbw.melt, aes(x = variable, y = value, group = id, color = change)) +
    geom_point() +
    geom_line() +
    facet_grid(. ~ low) +
    scale_color_manual(name = "Reclassification",
                       values = c("Up"="red", "Down"="blue"),
                       limits = c("Up","Down")) +
    scale_x_discrete(name = "Models") +
    scale_y_continuous(name = "Predicted probability", limit = c(0,1)) +
    labs(title = "Faceted by outcome")

非ケースの場合（アウトカム=0）、ダウンリクラシフィケーション（青）は良好。 ケース(結果=1)の場合、アップ再分類(赤)は良いことです。

その逆は誤った再分類である。灰色の線で示された患者は再分類されなかったもの。

Continuous NRI

ggplot(lbw.melt, aes(x = variable, y = value, group = id, color = reclass)) +
    geom_point() +
    geom_line() +
    facet_grid(. ~ low) +
    scale_color_manual(name = "Reclassification",
                       values = c("Up"="red", "Down"="blue", "Unchanged"="gray"),
                       limits = c("Up","Unchanged","Down")) +
    scale_x_discrete(name = "Models") +
    scale_y_continuous(name = "Predicted probability", breaks = c(0, 0.4, 1), limit = c(0,1)) +
    labs(title = "Faceted by outcome")

Calibration plot and Hosmer-Lemeshow goodness of fit test

plotCalibration(data = lbw, 
                cOutcome = 2, 
                predRisk = fitted(logistic.model.list[[3]]), 
                groups = 10)

## $Table_HLtest
##                total meanpred meanobs predicted observed
## [0.0589,0.168)    19    0.125   0.105      2.37        2
## [0.1681,0.225)    21    0.202   0.190      4.25        4
## [0.2255,0.255)    17    0.237   0.294      4.03        5
## [0.2548,0.274)    19    0.264   0.211      5.01        4
## [0.2738,0.299)    19    0.285   0.421      5.42        8
## [0.2987,0.338)    19    0.321   0.316      6.11        6
## [0.3383,0.372)    23    0.358   0.261      8.23        6
## [0.3716,0.415)    15    0.390   0.200      5.85        3
## [0.4152,0.483)    20    0.443   0.600      8.85       12
## [0.4829,0.597]    17    0.523   0.529      8.89        9
## 
## $Chi_square
## [1] 7.619
## 
## $df
## [1] 8
## 
## $p_value
## [1] 0.4715

NRIの解釈

モデルのAUCROCがほぼ同じ場合 予測因子の追加によって予後予測が実質的に改善されることはない
再分類は予測因子の追加が医学的判断を変えることを示唆する

この場合は、再分類がどの程度正しいか（個人が正しい方向に移動したかどうか）を調査し、２つのモデルを用いて再分類された個人がどのくらい正しい再分類であるかを明らかにする。

よくあるのは、再分類は悪い再分類とほぼ同数の良い再分類をもたらしている時に、AUROCは同等だがNRIが改善するというようなことがおきる。リスクファクターが予測を向上させない場合，再分類は単に更新されたモデルが異なる誤りをすることを示すだけである．

新しいリスクファクターの増分値を評価することは、逐次的なプロセスであり、AUCは、モデルの予測性能の改善に関する第一印象に非常に適した尺度である。更新されたモデルのAUCが元のモデルと顕著な差がない場合、再分類、健康への影響、費用対効果をさらに評価することは保証されない。再分類は、予測性能の改善によって裏付けられた場合にのみ、臨床的有用性を示す。再分類は、AUCの分析に取って代わるものではなく、むしろそれに続くものであるべきである。

使用における注意

最近では、複数の研究者がNRIの現在の用途を見直し、不適切な使用に関する懸念を表明している。 Improvement of risk prediction by genomic profiling: reclassification measures versus the area under the receiver operating characteristic curve. American Journal of Epidemiology 2010; 172(3):353–361. DOI: 10.1093/aje/kwq122.
Assessment of improved prediction beyond traditional risk factors: when does a difference make a difference? Circulation: Cardiovascular Genetics 2010; 3(1):3–5. DOI: 10.1161/CIRCGENETICS.110.938092.
Clinically relevant measures of fit? A note of caution. American Journal of Epidemiology 2012; 176(6):488–491. DOI: 10.1093/aje/kws208.
Does the net reclassification improvement help us evaluate models and markers? Annals of Internal Medicine 2014; 160(2):136–137. DOI: 10.7326/M13-2841.Does the net reclassification improvement help us evaluate models and markers? Annals of Internal Medicine 2014; 160(2):136–137. DOI: 10.7326/M13-2841.
Net reclassification indices for evaluating risk prediction instruments: a critical review. Epidemiology 2014; 25(1):114–121. DOI: 10.1097/EDE.0000000000000018.

予測モデルを治療の意思決定に用いる場合、予測されたリスクはカテゴリーであれば、高いリスクと低いリスクに分類される。

これらの代替指標は、正しい再分類の定義において大きく異なっている。再分類が正しいとみなされるのは、再分類されたグループで観察されたリスクが新しいリスクカテゴリー内に収まるか、それに近づくときだけでなく、病気を発症する人がより高いリスクカテゴリーに変わり、発症しない人がより低いリスクカテゴリーに変わるときでもある。再分類の割合は、選択した閾値の値と数に依存するため、臨床的に適切なリスク閾値が利用できる場合にのみ、再分類を評価することが有意義である。 。

Pencinaらは、NRIと統合的識別性改善（IDI）の紹介論文で、「IDI（ひいてはNRIも）はモデルのキャリブレーションに依存する」ことを指摘し、「識別スロープはモデルのキャリブレーションに依存するという欠点があり、IDIにも同様の影響があるかもしれない」「新しいマーカー追加後のモデルの性能を評価する場合、キャリブレーションで改善（あるいは少なくとも他の指標が改善すれば悪影響がない）されているかどうかを確認することが不可欠」と指摘している。

HildenとGerdsは、予測モデルへのマーカーの付加価値を評価する際に、「統合識別改善（IDI）とNRIに頼らないこと」と「これらの指標は信頼できない指針を提供すること」を警告している。 A note on the evaluation of novel biomarkers: do not rely on integrated discrimination improvement and net reclassification index. Statistics in Medicine 2014, DOI: 10.1002/sim.5804.

彼らは理論的な例と洗練された数学的展開を用いて、「IDIとNRIが予測性能の向上を測定するために使われる場合、キャリブレーションが不十分なモデルが有利に見えることがあり、シミュレーション研究では、実際にデータを生成したモデル（したがって、可能な限り最高のモデル）でも、測定情報を追加せずに改善できる」と結論付けている。このことから、IDIと連続NRIは、AUCやBrierスコアとは対照的に、「非正規」である測定基準の範疇に入る。

シミュレーションや数学的な発展がなくても、HildenやGerdsと同じ結論に達することは容易である。例えば、研究対象者全員の予測リスクを単純に2倍にすると、通常はかなりの割合が他のリスクカテゴリーに再分類され、予測リスクの差が膨らむことになる。これらの知見の妥当性に関する疑問は、HildenとGerdsの論文では納得のいく回答が得られていない

Continuous NRIへの批判

Continuous NRIは、臨床的な有用性を評価するためには最も適切でない可能性がある。特に、連続的な純再分類改善は、モデルの改善というよりも効果の指標であり、時に不規則な挙動を示すことがある。これを改善の指標として用いるには、注意が必要である。また、Continuous NRIの検定と統合的な識別改善の検定は、入れ子モデルにおける尤度比検定に類似しており、過大に解釈される可能性がある。リスク層での再分類は、閾値が必要であるが、治療方針の変更の可能性を検討することができ、臨床的により適切であると思われる。 American Journal of Epidemiology, Volume 176, Issue 6, 15 September 2012, Pages 488-491,https://doi.org/10.1093/aje/kws208

References