CPM Practical4 performance

Model performance

予測は、相対的な効果（OR、RRなど）指標で表現されるだけではなりません。
私たちが主に興味を持っているのは、絶対的な効果です。
例えば、

“この特性とこの特性を持つ人が（結果を）得る確率はどのくらいか？”

などです。

予測モデルを使用するには、少なくともモデル開発に使用した対象については、これらの予測が有効であることが望まれます（内部検証）。

また、予後予測モデルの開発に使用したサンプル以外の他の集団におけるモデルの性能や予測の妥当性を推論しなければなりません（外部検証）
　　 多くの場合、新しいマーカー（遺伝子マーカーや画像マーカーなど）が、対象となる疾患を予測する上で重要かどうかを判断する必要があります。マーカーとアウトカムの「統計的に特異な」関連性だけでは十分ではなく、新しいマーカーは、モデルの予測能力を向上させなければなりません。また、マーカーは多くの候補マーカー（潜在的な遺伝子マーカーなど）の中から選ばれることが多いため、この新しいマーカーの性能を過大評価してしまう恐れがあります。

これらの理由から代替モデル（旧モデルと新マーカーを使用したモデル）を選択するための予測性能（パフォーマンス）の指標が必要です。　　

従来の統計的アプローチでは、予測が実際の結果にどれだけ近いかを測定していた。 例えば、R2、Brier scoreがあります。
現在、予測性能を推定するのに重要な指標としてはDiscliminationとCalibrationがあります。 ちなみに、R2やBrierはDiscriminationの指標になります。

Discrimination

Discrimination（識別能）はアウトカムを持つ人が持たない人よりもリスク予測が高いかどうかを確認することです. 二値の結果の場合、予測値p（0と1の間の確率）は、閾値tに従って二値化される必要があります。これにより、p<tの場合は0、p≧tの場合は1が割り当てられ、これによりSensitivity、Specificityなどを計算することができます。

• Sensitivity、Specificity
  

• ROCは、1-Specificityに対するSensitivityをプロット
  

• AUCとC score
  

• Discrimination slope(Yates の判別係数): 結果が出ている人と出ていない人の平均予測値の差 　　 # Calibration Calibration（較正能）は予測確率がα％の人たちが、平均してα％の人たちが実際にアウトカムがおきているかをチェックするものです.つまり予測値が真の値に近いかを評価します。二分法の場合、真の事象はY=0または1であるため、グループを作成するし予測された確率pは真値Yに近いかを確認します。

• キャリブレーションプロット: 予測平均値と真の平均値のプロットは対角線上にあるのが望ましいです

Y=0,1 真の結果
p= 予測される結果の確率
　　 とした場合に、logit(p)は線形予測式であらわされます。
キャリブレーションプロットを回帰線で表すと

\(logit(Y) = α_{cal} + β_{cal}×logit(p)\)

\(β_{cal}\)は校正の傾き よくキャリブレーションされたモデルは、pはYに近い値になると予想されるため \(α_{cal}\)、切片が0
\(β_{cal}\)、傾きは1
になるはずです

• キャリブレーションスロープ
  予測値は最尤法を用いて最適なフィット感を得ているため、モデルの開発（見かけ上のキャリブレーション）に使用したデータセットについては、 \(β_{cal}=１\)
  となるはずです。

内部検証では、

optimismのため\(β_{cal}<1\)

となりますこれは、モデルが誤キャリブレーションされていることを示し、観察された集団と同様の新しい集団において予測変数の効果を過大評価することになります。

また、モデルの係数に必要なshrinkを示しています。
外部検証では、\(β_{cal}\)は潜在的なオーバーフィッティングを反映するだけでなく、予測変数の効果の違いも反映します

他の性能評価として
- モデルの臨床的有用性には、偽陽性と真陽性の相対的価値の判断が含まれること
- 決定曲線分析
があります。

R package

この章ではROCR, pROCパッケージを使用してROC曲線をフィットさせ、AUCを推定します。
rmsパッケージには、内部検証のためのコマンドが含まれています。
foreignパッケージでは、他のソフトウェア（ここではSPSS）からデータファイルをインポートすることができます。

rm(list=ls())
library(ROCR)
library(pROC)
library(rms)
library(foreign)
library(tidyverse)

Dataset

今回は、精巣がんのデータを使用します。
転移性の非セミノーマ性精巣癌に対する化学療法後には、残存腫瘍（成熟した奇形腫や生存している癌細胞）がまだ存在している可能性があるため、最初の転移の残骸を取り除くために外科的切除が一般的に受け入れられている治療法です。この実習では、残存腫瘍の有無を予測するロジスティック回帰モデルを開発し、検証します。

データセットt821は、変数NEC（1は壊死; 良性の組織、0は腫瘍）に関する4つの研究からのデータを含んでいます。
データは以下から引用されています：
Steyerberg EW et al. European Journal of Cancer, 30 (9), 1994, 1231-1239;
Steyerberg EW et al. 1995 13(5):1177-87;
Vergouwe Y et al. J Urol. 2001 Jan;165(1):84-8、
再解析は
Steyerberg EW et al. Statist. Med. 2001; 20:3847-3859.

まず、データを読み込みます。

t821  <- read.spss('t821.sav',use.value.labels=F, to.data.frame=T)

このデータセットには、NECを予測できる5つの変数が含まれています：
- 原発腫瘍のテラトーマ要素
- 化学療法前のAFPとHCGのレベル
- 化学療法後の腫瘤の大きさ
- 腫瘤の大きさの減少
です。
NECを調査した研究にしたがって、3つのデータセットD1からD3に分割します。

D1  <- t821[t821$STUDY==1,]
D2  <- t821[t821$STUDY==3,] 
D3  <- t821[t821$STUDY==2 | t821$STUDY==4,] 

Developing the model and estimation of the in-sample discrimination and calibration

まず，データセットD1を見てみましょう

head(D1)

ここでは，D1の5つの予測変数(TER,PREAFP,PREHCG,SQPOST,REDUC10)を用いてlrmでモデルをフィットさせます。glmだとrmsとの組み合わせが悪いからです。

x=Tと入力するとexpanded design matrixがxに返ってきます
y=Tと入力すると、応答変数（欠測を除く）がyという名前で返されます。

model1 <- lrm(NEC ~ TER+PREAFP+PREHCG+SQPOST+REDUC10,
            data=D1,x=T,y=T,linear.predictors=T)
model1

## Logistic Regression Model
##  
##  lrm(formula = NEC ~ TER + PREAFP + PREHCG + SQPOST + REDUC10, 
##      data = D1, x = T, y = T, linear.predictors = T)
##  
##                         Model Likelihood    Discrimination    Rank Discrim.    
##                               Ratio Test           Indexes          Indexes    
##  Obs           544    LR chi2     186.90    R2       0.389    C       0.818    
##   0            299    d.f.             5    g        1.760    Dxy     0.636    
##   1            245    Pr(> chi2) <0.0001    gr       5.810    gamma   0.637    
##  max |deriv| 8e-09                          gp       0.317    tau-a   0.315    
##                                             Brier    0.174                     
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept -2.1072 0.5797 -3.63  0.0003  
##  TER        0.9955 0.2082  4.78  <0.0001 
##  PREAFP     0.8592 0.2281  3.77  0.0002  
##  PREHCG     0.5541 0.2190  2.53  0.0114  
##  SQPOST    -0.0737 0.0658 -1.12  0.2626  
##  REDUC10    0.2644 0.0495  5.35  <0.0001 
## 

参考にglmでモデルを作っておきます

gm1 <- glm(NEC ~ TER+PREAFP+PREHCG+SQPOST+REDUC10,
            data=D1,family=binomial)

summary(gm1)

## 
## Call:
## glm(formula = NEC ~ TER + PREAFP + PREHCG + SQPOST + REDUC10, 
##     family = binomial, data = D1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4076  -0.8362  -0.3334   0.8835   2.4087  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.10715    0.57975  -3.635 0.000278 ***
## TER          0.99546    0.20816   4.782 1.73e-06 ***
## PREAFP       0.85916    0.22812   3.766 0.000166 ***
## PREHCG       0.55414    0.21897   2.531 0.011386 *  
## SQPOST      -0.07371    0.06579  -1.120 0.262551    
## REDUC10      0.26442    0.04946   5.346 9.00e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 748.78  on 543  degrees of freedom
## Residual deviance: 561.87  on 538  degrees of freedom
## AIC: 573.87
## 
## Number of Fisher Scoring iterations: 5

AUC

さて、AUCはどのくらいでしょうか？
それでは、2つの異なるアプローチでROC曲線をプロットします。
まず、モデルの線形予測変数にlogitの逆数をとることによって、予測される確率を計算する必要があります。
model1$linear.predictorsに線形の予測変数が格納されているのでそれをplogis関数で予測確率にもどします。

prob=plogis(model1$linear.predictors)
#prob=plogis(predict(model1,type=c("lp"))) #これは上と同じ

GLMでロジスティック回帰の場合は

D1$prob1 <- predict(gm1,type=c("response"))

で予測確率がもとまります

pROC package

次に、pROCパッケージを使って曲線をプロットします。

rocは、pROCパッケージの主な機能です。　　 この関数はROC曲線を構築し、クラス “roc”のリストである “roc”オブジェクトを返します。このオブジェクトは、print、plot、または関数 auc, ci, smooth.roc, coordsに渡すことができます。
さらに、2つのrocオブジェクトをroc.testで比較することができます。

roc(response, predictor)でROC曲線を生成できます。
D1$NECは１，０変数でアウトカム(response)です。probがpredictorです。

rocmodel1 <- roc(D1$NEC, prob)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

plot(rocmodel1)

rocmodel1$auc

## Area under the curve: 0.818

glmモデルを指定して書く場合は

roc1 <- roc(NEC~prob1, data = D1)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

(AUC1 <- roc1$auc)      # AUC estimate

## Area under the curve: 0.818

ci(roc1)        # 95% CI of AUC

## 95% CI: 0.7833-0.8527 (DeLong)

plot(roc1)  # plot the ROC curve

ROCR package

ROCR パッケージを使ってプロットしてみます。
まずはprediction {ROCR}を使って予測オブジェクトを生成します。
ROCRを使った分類器の評価は、すべて予測オブジェクトの作成から始まります。この関数は，入力データ（ベクトル，行列，データフレーム，リストの形式がある）を標準的な形式に変換するために使用されます．

Ex prediction(predictions, labels, label.ordering = NULL)

model1$yにはlrmで作成した応答変数の行列が格納されています。

pred<-prediction(prob,model1$y)
pred

## A prediction instance
##   with 544 data points

これでpredには予測インスタンス（確率とアウトカムがリストになっている）が生成されています。

performance {ROCR}はpredictorの性能をみるための関数です。
performance(prediction.obj, measure, x.measure = “cutoff”, …)
prediction.objにpredictionを入れmesureで評価指標をいれます
tpr:True positive rate. P(Yhat = + | Y = +). Estimated as: TP/P.
fpr:False positive rate. P(Yhat = + | Y = -). Estimated as: FP/N.
です

ROC curves:measure=“tpr”, x.measure=“fpr”. の設定でRCOが生成されます。

perf <- performance(pred, measure = "tpr", x.measure = "fpr")

プロットします。

plot(perf, col=rainbow(7), main="ROC curve", xlab="1-Specificity (FPR)", 
     ylab="Sensitivity (TPR)")    
abline(0, 1) 

discrimination slopeは、壊死のあるものとないものの予測確率の平均値の差です。

mean(prob[D1$NEC==1])-mean(prob[D1$NEC==0])

## [1] 0.300757

ちなみにAUCの比較は以下の用に行う roc.test(roc1, roc2, method=“delong”) # test of the difference of AUC between roc1 and roc2

Calibration

サンプル内のキャリブレーションの傾きはどのくらいですか？
validateのSlopeにCalibration slopeの傾きが記載されます。
optimismを補正したmのがindex,correctedなslopeにあたる。

予測値は最尤法を用いて最適なフィットを得ているので、開発モデルでは1になります。 しかし、optimismを含んだ見かけ上のキャリブレーションです。

model1Bvalidation<-validate(model1, B=200) 
model1Bvalidation

##           index.orig training   test optimism index.corrected   n
## Dxy           0.6361   0.6385 0.6299   0.0086          0.6275 200
## R2            0.3890   0.3936 0.3826   0.0110          0.3779 200
## Intercept     0.0000   0.0000 0.0053  -0.0053          0.0053 200
## Slope         1.0000   1.0000 0.9738   0.0262          0.9738 200
## Emax          0.0000   0.0000 0.0067   0.0067          0.0067 200
## D             0.3417   0.3473 0.3350   0.0122          0.3295 200
## U            -0.0037  -0.0037 0.0004  -0.0040          0.0004 200
## Q             0.3454   0.3510 0.3347   0.0163          0.3291 200
## B             0.1737   0.1722 0.1757  -0.0035          0.1772 200
## g             1.7596   1.7873 1.7324   0.0549          1.7047 200
## gp            0.3167   0.3175 0.3138   0.0037          0.3130 200

なおブートストラップしていますが開発データでみるのはorigの値です。

Calibration plotを手書きで書いてみます。 数理研の野間先生のCodeをこのデータにフィットさせています。

## Calibration plot and Hosmer-Lemeshow test
# yは真の値
# yhatは期待値
# gはグラフの分割（10分割がデフォルト）

calib <- function(y, yhat, g=10) { # create calib function
  cutyhat <- cut(yhat, breaks = quantile(yhat, probs=seq(0, 1, 1/g)), include.lowest=TRUE)  #予測値を10分割してcutyhatに格納
  obs <- xtabs(cbind(1 - y, y) ~ cutyhat)
  expect <- xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
  tot <- xtabs(rep(1,times=length(y)) ~ cutyhat)
  chisq <- sum((obs - expect)^2/expect)
  P <- 1 - pchisq(chisq, g - 2)

  O <- obs[,2]/tot
  E <- expect[,2]/tot
  R <- data.frame(O,as.vector(E))
  colnames(R) <- c("cutyhat","Obs","Exp")
  
  plot(as.vector(O)~as.vector(E), xlim=0:1, ylim=0:1, ylab="Observed Frequency", xlab="Predicted Probability")   # calibration plot
  abline(0,1, lty=2)
  
  lm1 <- summary(lm(as.vector(O)~as.vector(E)))
  
  return(list(datatable=R, calib.line=lm1, Chisq=chisq, p.value=P))
  
}

calib(y=D1$NEC, yhat=fitted(gm1))

## $datatable
##             cutyhat        Obs        Exp
## 1  [0.00407,0.0709] 0.01785714 0.04350387
## 2    (0.0709,0.164] 0.07547170 0.11883063
## 3     (0.164,0.263] 0.24074074 0.21000073
## 4     (0.263,0.345] 0.27272727 0.30291883
## 5     (0.345,0.436] 0.46296296 0.38975436
## 6     (0.436,0.536] 0.56896552 0.49191920
## 7     (0.536,0.632] 0.59615385 0.59639749
## 8     (0.632,0.736] 0.66037736 0.68244295
## 9     (0.736,0.847] 0.74074074 0.78972480
## 10    (0.847,0.945] 0.87272727 0.88836093
## 
## $calib.line
## 
## Call:
## lm(formula = as.vector(O) ~ as.vector(E))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.04642 -0.03084 -0.01474  0.02195  0.07799 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.004341   0.029608   0.147    0.887    
## as.vector(E) 0.989246   0.056155  17.616  1.1e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04839 on 8 degrees of freedom
## Multiple R-squared:  0.9749, Adjusted R-squared:  0.9717 
## F-statistic: 310.3 on 1 and 8 DF,  p-value: 1.102e-07
## 
## 
## $Chisq
## [1] 6.011022
## 
## $p.value
## [1] 0.6459972

val.prob(prob, D1$NEC, m=nrow(D1)/10, cex=.5)

##           Dxy       C (ROC)            R2             D      D:Chi-sq 
##  6.360385e-01  8.180192e-01  3.889772e-01  3.417348e-01  1.869037e+02 
##           D:p             U      U:Chi-sq           U:p             Q 
##            NA -3.676471e-03  2.273737e-13  1.000000e+00  3.454113e-01 
##         Brier     Intercept         Slope          Emax           E90 
##  1.736997e-01 -5.896206e-12  1.000000e+00  5.115574e-02  3.864631e-02 
##          Eavg           S:z           S:p 
##  2.355920e-02  1.541363e-01  8.775023e-01

さらに単純にrms::calibrateで書く場合
glmでなくlmrでないと対応していない

cali<-calibrate(model1, method = "boot", B=100,
                predy = seq(.01,.99,length=100))
plot(cali)

## 
## n=544   Mean absolute error=0.026   Mean squared error=0.00088
## 0.9 Quantile of absolute error=0.041

ちょっとださい気もします。

Internally validated discrimination and calibration

上記のcalibrationの値もoptimismの対象となります。
つまり、実際の性能はサンプルから得られたものよりも低いということです。

optimismを修正するために、ブートストラップを使用して、モデルの性能を新たに推定します。

rmsパッケージのvalidationコマンドがこれを行います。

model1Bvalidation<-validate(model1, B=200) 
model1Bvalidation

##           index.orig training    test optimism index.corrected   n
## Dxy           0.6361   0.6422  0.6298   0.0124          0.6236 200
## R2            0.3890   0.3975  0.3821   0.0154          0.3736 200
## Intercept     0.0000   0.0000 -0.0090   0.0090         -0.0090 200
## Slope         1.0000   1.0000  0.9658   0.0342          0.9658 200
## Emax          0.0000   0.0000  0.0092   0.0092          0.0092 200
## D             0.3417   0.3515  0.3346   0.0169          0.3248 200
## U            -0.0037  -0.0037  0.0003  -0.0040          0.0003 200
## Q             0.3454   0.3551  0.3343   0.0209          0.3245 200
## B             0.1737   0.1715  0.1757  -0.0042          0.1779 200
## g             1.7596   1.8021  1.7293   0.0728          1.6867 200
## gp            0.3167   0.3195  0.3136   0.0059          0.3108 200

内部で検証されたAUCはどのくらいでしょうか？
この関数は、AUCの代わりにSomer’s DをDxyとして返します（理由は不明ですが、これは実に厄介です）。
AUCを得るためには、Somer’s D / 2 + 0.5という式を使います。

#get internally validated AUC -  Somer's D / 2 + .5
model1Bvalidation[1,]/2+0.5

##      index.orig        training            test        optimism index.corrected 
##       0.8180329       0.8211159       0.8148937       0.5062222       0.8118107 
##               n 
##     100.5000000

較正用の傾きはどれくらいでしょうか？

#Calibration internally validated
model1Bvalidation[4,] #4列目がslopeなので

##      index.orig        training            test        optimism index.corrected 
##       1.0000000       1.0000000       0.9658207       0.0341793       0.9658207 
##               n 
##     200.0000000

index.correctedの値をみればOKです。

ちなみに、手書きでブートストラップしたければ以下のようにやります。

## Adjustment of the optimism by bootstrap

B <- 1000   # number of resampling

N <- dim(D1)[1]

AUC.i1 <- AUC.i2 <- numeric(B)

for(i in 1:B){
  bs.i <- sample(1:N, N, replace=TRUE)
  D1.i <- D1[bs.i,]
  gm1.i <- gm1 <- glm(NEC ~ TER+PREAFP+PREHCG+SQPOST+REDUC10,
            data=D1.i,
            family=binomial)
  D1.i$prob1 <- predict(gm1.i,type=c("response"))
  AUC.i1[i] <- roc(NEC ~ prob1, data = D1.i)$auc
  D1$prob.i <- predict(gm1.i, newdata=D1,type = c("response"))
  AUC.i2[i] <- roc(NEC ~ prob.i, data = D1)$auc
  print(paste(i,"th bootstrap iteration is completed.", sep=""))
  }

opt1 <- AUC.i1 - AUC.i2
summary(opt1)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.061566 -0.006982  0.005884  0.004976  0.016556  0.050675

hist(opt1)

lam1 <- mean(opt1)      # estimate of the optimism
(AUC.c1 <- AUC1 - lam1)     # bias corrected AUC estimate

## [1] 0.8130565

External validation using the dataset D2

D2としてコード化したインディアナ州のクリニックのデータを使って、model1の性能を評価します。

head(D2)

dim(D2)

## [1] 273  13

先ほどと同じモデルを使いますが、今度はそれをD2での予測に適用します。　　 　　 まず線形予測値を計算し、続いて逆ロジット関数を使ってD2患者の壊死の確率を計算します。

#use the same model as before, but now apply it to make predictions in D2
linpredD2<-predict(model1,D2)
predD2 <- plogis(linpredD2)

ここで、新しいAUCとcalibration slopeを計算します。どれくらいでしょうか？ 開発データセットD1の時と比べて、判別やキャリブレーションが悪くなったのか、良くなったのか？

val.prob {rms}は二値事象に対する予測確率を検証するのに役立ちます． 予測確率pまたは予測対数オッズlogitのセットと，予測pまたはlogitの開発に使用されなかった二値アウトカムyのベクトルが与えられると，val.probは以下の指標と統計量を計算します．

• pとyの間のSomersのD_{xy}
  
• 順位相関[2(C-.5)，C=ROC領域]
  
• Nagelkerke-Cox-Snell-Maddala-MageeのR二乗指数
  
• 判別指数D [(ロジスティックモデルL.R.カイ二乗-1)/n]
  
• L.R.カイ二乗，そのP値
  
• 信頼性のない指数U
  
• 信頼性のないことを検定するための2d.f.のカイ二乗
  
• 非信頼性を検定するための2 d.f.のカイ二乗(H0: intercept=0, slope=1)、そのP値
  
• 品質指数Q
  
• Brierスコア(pとyの平均二乗差)
  
• Intercept
  
• Slope
  
• E_{max}=予測確率と loessで校正された確率の最大絶対差
  
• Eavg
  
• 同平均
  
• E90 , 0.9分位点
• 校正精度のSpiegelhalter Z-test、その両側P値

#AUC and calibration slope
val.prob(predD2, D2$NEC, m=20, cex=.5)

##          Dxy      C (ROC)           R2            D     D:Chi-sq          D:p 
##  0.570331285  0.785165643  0.249804398  0.186605912 51.943414017           NA 
##            U     U:Chi-sq          U:p            Q        Brier    Intercept 
##  0.007487938  4.044207019  0.132376716  0.179117974  0.160631870 -0.158494642 
##        Slope         Emax          E90         Eavg          S:z          S:p 
##  0.740150183  0.162098288  0.046780716  0.027916813  1.466675457  0.142464363

m=20はグループ化された比率が必要な場合、グループごとの平均オブザベーション数
cexはレジェンドの文字サイズ

DxyはSomers D C(ROC)がAUCになります。

モデルのキャリブレーションについてはどうでしょうか？ Slopeが0.74で1に比べるとだいぶ悪いです。

Re-calibration of a model for a new setting

model1のキャリブレーションを改善し、より良い予測をするために、モデルの切片と傾きの両方を再キャリブレーションします。

model1を再校正するために適合させる新しい回帰モデルを書きます

#re-calibrate intercept and slope
recabmodel1<-lrm(NEC~linpredD2,data=D2)

AUCはどのくらいになりましたか？
モデルはよくキャリブレーションされたか？

#check that now slope is 1
val.prob(plogis(recabmodel1$linear.predictors), D2$NEC, m=20, cex=.5)

##           Dxy       C (ROC)            R2             D      D:Chi-sq 
##  5.703313e-01  7.851656e-01  2.673319e-01  2.014199e-01  5.598762e+01 
##           D:p             U      U:Chi-sq           U:p             Q 
##            NA -7.326007e-03  1.136868e-13  1.000000e+00  2.087459e-01 
##         Brier     Intercept         Slope          Emax           E90 
##  1.599225e-01  6.977180e-17  1.000000e+00  5.027273e-02  4.266982e-02 
##          Eavg           S:z           S:p 
##  1.926467e-02 -2.181190e-01  8.273364e-01

Additional work: Extend the model with LDH as predictor

追加作業です。
LDHを予測因子としてモデルを拡張する 潜在的に重要なマーカーの1つである精巣乳酸脱水素酵素(LDH)は、モデルの改善に使用できます。まず、model1をLDHを含むmodel2にアップデートします。

# Extend the model with LDH #
model2 <- update(model1, ~ . + LNLDHST)

次に、内部で検証された性能を推定します。AUC（楽観性の補正）はどのくらいですか？そのモデルはよくキャリブレーションされているでしょうか？

#Internally validated characteristics
model2Bvalidation<-validate(model2, B=200)

#get internally validated AUC -  Somer's D / 2 + .5
model2Bvalidation[1,]/2+0.5

##      index.orig        training            test        optimism index.corrected 
##       0.8387141       0.8407255       0.8352555       0.5054700       0.8332441 
##               n 
##     100.5000000

#Calibration internally validated
model1Bvalidation[4,]

##      index.orig        training            test        optimism index.corrected 
##       1.0000000       1.0000000       0.9658207       0.0341793       0.9658207 
##               n 
##     200.0000000

あなたはどちらのモデルを選びますか？model1とmodel2のどちらがいいですか？ データセットD3は、model2を外部から評価するために使うことができます。

#### Externally validated characteristics
linpredD3<-predict(model2,D3)
predD3 <- plogis(linpredD3)

#AUC and calibration slope
valD3<-val.prob(predD3, D3$NEC, m=20, cex=.5)

print(valD3)

##          Dxy      C (ROC)           R2            D     D:Chi-sq          D:p 
##  0.647120943  0.823560471  0.376233495  0.319475439 89.494696661           NA 
##            U     U:Chi-sq          U:p            Q        Brier    Intercept 
##  0.002963517  2.820894089  0.244034165  0.316511923  0.165273693 -0.156465211 
##        Slope         Emax          E90         Eavg          S:z          S:p 
##  0.824818732  0.091738900  0.064673109  0.033788668  1.188036701  0.234818951

明らかにモデルの再キャリブレーションが必要です。

####recalibrate intercept and slope
recabmodel2<-lrm(NEC~linpredD3,data=D3)

#check that now slope is 1
val.prob(plogis(recabmodel2$linear.predictors), D3$NEC, m=20, cex=.5)

##           Dxy       C (ROC)            R2             D      D:Chi-sq 
##  6.471765e-01  8.235883e-01  3.862287e-01  3.296592e-01  9.231559e+01 
##           D:p             U      U:Chi-sq           U:p             Q 
##            NA -7.220217e-03  0.000000e+00  1.000000e+00  3.368794e-01 
##         Brier     Intercept         Slope          Emax           E90 
##  1.640637e-01  6.374558e-17  1.000000e+00  4.154895e-02  3.163878e-02 
##          Eavg           S:z           S:p 
##  1.666310e-02 -1.128883e-01  9.101191e-01

Scoring system

係数を点数化して予測モデル化していく方法を試します

gm2 <- glm(NEC ~ TER+PREAFP+PREHCG,data=D1,family=binomial)
library(broom)
tidy(gm2)

beta係数(estimate)が

TER1.17
PREAFP0.97
PREHCG0.48

なので便宜上それぞれ2,2,1点として係数を作成。
結果がどの程度あてはまりを持つかはこの後のステップで検証します

scoreという変数をTER, PREAFP, PREHCGの足し算で作成します

D1%>%mutate(score=2*TER+2*PREAFP+1*PREHCG)->D1 #mutate関数でscoreを作成

table(D1$score,D1$NEC)

##    
##       0   1
##   0 116  34
##   1  35  10
##   2  80  73
##   3  44  63
##   4  15  21
##   5   9  44

では、実際にDiscrimination、Calibrationを見ていきましょう。 念の為必要なパッケージを明示的に示しておきます。 polsplineは先にinstallしてください。

library(rms) #calibration
library(polspline)
library(ResourceSelection) #Hosmer-lemeshow

## ResourceSelection 0.3-5   2019-07-22

library(pROC) #AUC

glmでモデルを作ってみます。

gm_score<-glm(NEC ~ score,
            data=D1,family=binomial)
summary(gm_score)

## 
## Call:
## glm(formula = NEC ~ score, family = binomial, data = D1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7888  -1.0879  -0.7026   1.0468   1.7435  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.27301    0.16059  -7.927 2.24e-15 ***
## score        0.52946    0.06385   8.293  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 748.78  on 543  degrees of freedom
## Residual deviance: 667.57  on 542  degrees of freedom
## AIC: 671.57
## 
## Number of Fisher Scoring iterations: 4

AUC

これまでの説明と同様にAUCを求めてみます。

予測値を格納します。

D1$prob_score <- predict(gm_score,type=c("response"))

pROCを使って、AUCをプロットします。

roc_score <- roc(NEC~prob_score, data = D1)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

(AUC1 <- roc_score$auc)     # AUC estimate

## Area under the curve: 0.7109

ci(roc_score)       # 95% CI of AUC

## 95% CI: 0.6687-0.7531 (DeLong)

plot(roc_score) # plot the ROC curve

Calibration

rmsではlrmとolsモデルしか対応しておらず、glmだとerrorになってします。
このためlrmで再度モデルを書きます。
ちなみにglmでなくGlmにするとrmsでも対応します。

f1<-lrm(NEC~ score,
            data=D1,x=T,y=T)
cal<-calibrate(f1, group=D1$NEC)

plot(cal)

## 
## n=544   Mean absolute error=0.015   Mean squared error=0.00034
## 0.9 Quantile of absolute error=0.018

やっぱりちょっとださいので、 lrmでの予測確率を変数に格納し、Calibration plotを描きます。

prob_score=plogis(f1$linear.predictors)
val.prob(prob_score, D1$NEC, m=20, cex=.5, pl=T)

## Warning in min(xx[xx > upper]): min の引数に有限な値がありません: Inf を返します

##           Dxy       C (ROC)            R2             D      D:Chi-sq 
##  4.218688e-01  7.109344e-01  1.855107e-01  1.474428e-01  8.120891e+01 
##           D:p             U      U:Chi-sq           U:p             Q 
##            NA -3.676471e-03  5.684342e-13  1.000000e+00  1.511193e-01 
##         Brier     Intercept         Slope          Emax           E90 
##  2.123582e-01  2.418950e-16  1.000000e+00  1.463054e-02  9.581675e-03 
##          Eavg           S:z           S:p 
##  7.626502e-03  1.084086e-02  9.913504e-01

Hosmer-Lemshow test

適合度の検定としてhosmer-Lemshow testがあります。
観測された事象率がモデル母集団のサブグループでの期待される事象率に適合するかどうかを評価します。 十分位数で分割した適合リスク値をサブグループとする。サブグループでの期待された事象率と観測された事象率が類似するモデルは “well calibrated” と呼ばれる。

hoslem.test(gm_score$y,gm_score$fitted.values)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  gm_score$y, gm_score$fitted.values
## X-squared = 3.5786, df = 8, p-value = 0.893

一般的にHosmer-Lemeshow検定ではp値が0.05より大きければ適合しているとみなします。

あとはValidationの手続きは上と同様に行っていけば良い。