Rで予測モデルの内的妥当性検証を行うためのパッケージはいくつか有るが、自分がほしい解析を一括で扱ってくれるいい感じのパッケージがなくて困っていた。
特に、c-statistics, calibration slope, calibration interceptのbias-corrected estimationを最尤推定法だけでなくshrinkage estimation (lasso, ridge)で評価できるパッケージは存在しない。
そこで、既存のパッケージと組み合わせていい感じのBias-corredted estimatonをしてくれる関数を自作してみた。
使用したパッケージは
pacman::p_load( tidyverse, #データ処理で必須 rms, #Harrellの回帰分析とにかく便利全部盛り glmnet, #lasso, ridgeには必要不可欠 CalibrationCurves # 予測性能指標を一括して算出してくれる関数が入っている )
まずは各種推定法における予測性能指標のデータフレームを算出する関数を自作する
MLE regression
MLE.pred<- function(var,data,newdata){ #var=variable, data= modeling data, newdata=predicted data form<-as.formula(paste("death30~",paste(var,collapse="+"))) #death30の部分は二値アウトカムなので使うときに自由に設定してください model<-lrm(form,data=data) res<-valProbggplot( logit=predict(model, newdata=newdata), y=as.numeric(newdata$death30)-1, pl=T) res<- res$stats %>% t() %>% as.data.frame() return(res) }
Lasso regression
lasso.pred<- function(var,data,newdata){ #var=variable, data= modeling data, newdata=predicted data form<-as.formula(paste("death30~",paste(var,collapse="+"))) mat<-model.matrix(form, data=data)[,-1] newmat<-model.matrix(form, data=newdata)[,-1] #set.seed(1) lasso.cv<-cv.glmnet(mat, data$death30, alpha = 1,family = "binomial") res<-valProbggplot( p=predict(lasso.cv, newx=newmat, type="response", s="lambda.min"), y=as.numeric(newdata$death30)-1, pl=T) res<- res$stats %>% t() %>% as.data.frame() return(res) }
ridge regression
ridge.pred<- function(var,data,newdata){ #var=variable, data= modeling data, newdata=predicted data form<-as.formula(paste("death30~",paste(var,collapse="+"))) mat<-model.matrix(form, data=data)[,-1] newmat<-model.matrix(form, data=newdata)[,-1] #set.seed(1) ridge.cv<-cv.glmnet(mat, data$death30, alpha = 0,family = "binomial") res<-valProbggplot( p=predict(ridge.cv, newx=newmat, type="response", s="lambda.min"), y=as.numeric(newdata$death30)-1, pl=T) res<- res$stats %>% t() %>% as.data.frame() return(res) }
ちなみにCalibrationCurvesパッケージのval.prob.ci.2を使っても予測性能指標は出るのだが、関数のエラーなのか、calibraiton plotの出力をoffにする(pl=F)にするとエラーが出てしまい、pl=TにするとBootstrap時にBootstrap回数枚のplot画像が出力されてしまうので、かわりにvalProbggplotを使用。
なぜかこちらだと結果を変数に格納すればplotが出ないので便利。各種指標の計算方法はval.prob.ci.2と一緒。
説明変数のベクトル、dataを入れてあげれば以下のようなデータフレームで出力される。 元々Bootstrapを想定してnewdataを変数に組み込んだ形の関数にしているが、apparent estimationをしたければ、data=newdataにしてあげればOK。
3つの推定方法(MLE, lasso, ridge)について、Bootstrap resampling法でBias-corrected estimateする自作関数
boot.bias.corr.estimation<- function(orig.data,var,B=5,method){ ## Adjustment of the optimism by bootstrap N <- dim(orig.data)[1] #patient number res.i1 <- res.i2 <- list() # MLE estimation if (method=="MLE"){ for(i in 1:B){ #bootstrap sample set.seed(i) bs.i <- sample(1:N, N, replace=TRUE) d.i <- orig.data[bs.i,] #bootstrap modelでのapparent prediction ability res.i1[[i]]<-MLE.pred(var,data=d.i,newdata=d.i) #bootstrap modelでoriginal dataのpredictionを外挿した結果 res.i2[[i]]<-MLE.pred(var,data=d.i,newdata = orig.data) #print(paste(i,"th bootstrap iteration is completed.", sep="")) } } # lasso estimation else if (method=="lasso"){ for(i in 1:B){ #bootstrap sample set.seed(i) bs.i <- sample(1:N, N, replace=TRUE) d.i <- orig.data[bs.i,] #bootstrap modelでのapparent prediction ability res.i1[[i]]<-lasso.pred(var,data=d.i,newdata=d.i) #bootstrap modelでoriginal dataのpredictionを外挿した結果 res.i2[[i]]<-lasso.pred(var,data=d.i,newdata = orig.data) #print(paste(i,"th bootstrap iteration is completed.", sep="")) } } # ridge estimation else if (method=="ridge"){ for(i in 1:B){ #bootstrap sample set.seed(i) bs.i <- sample(1:N, N, replace=TRUE) d.i <- orig.data[bs.i,] #bootstrap modelでのapparent prediction ability res.i1[[i]]<-ridge.pred(var,data=d.i,newdata=d.i) #bootstrap modelでoriginal dataのpredictionを外挿した結果 res.i2[[i]]<-ridge.pred(var,data=d.i,newdata = orig.data) #print(paste(i,"th bootstrap iteration is completed.", sep="")) } } res.i1.df <- do.call(rbind, res.i1) res.i2.df <- do.call(rbind, res.i2) # Extract the relevant metrics and compute the difference opt1 <- list() for (i in 1:B) { # Assuming the result data frames have columns named "AUC", "slope", and "intercept" auc_diff <- res.i1[[i]]$`C (ROC)` - res.i2[[i]]$`C (ROC)` slope_diff <- res.i1[[i]]$Slope - res.i2[[i]]$Slope intercept_diff <- res.i1[[i]]$Intercept - res.i2[[i]]$Intercept brier_diff <- res.i1[[i]]$Brier - res.i2[[i]]$Brier opt1[[i]] <- data.frame(AUC = auc_diff, slope = slope_diff, intercept = intercept_diff, brier = brier_diff) } opt1.df <- do.call(rbind, opt1) save(opt1.df, file="opt1.df.RData") # bias-corrected estimation original.res<- if(method=="MLE"){ MLE.pred(var,data=orig.data,newdata=orig.data) %>% select("AUC"="C (ROC)","slope"="Slope","intercept"="Intercept","brier"="Brier") } else if(method=="lasso"){ lasso.pred(var,data=orig.data,newdata=orig.data) %>% select("AUC"="C (ROC)","slope"="Slope","intercept"="Intercept","brier"="Brier") } else if(method=="ridge"){ ridge.pred(var,data=orig.data,newdata=orig.data) %>% select("AUC"="C (ROC)","slope"="Slope","intercept"="Intercept","brier"="Brier") } opt1df.mean<- opt1.df %>% summarise_all(mean) bias.corrected.estimation<- original.res - opt1df.mean result<-rbind(original.res, bias.corrected.estimation) %>% mutate(rowname=c("original","bias_corrected"),.before=AUC) return(list(optimism.data=opt1.df, optimisim.hist=hist(opt1.df), result=result)) }
以下のように、original data, 説明変数のベクトル(var)、Bootstrap回数(B)、推定方法(method)を定義すると、
apparent estimateとbias-corrected estimateが比較されたデータフレーム形式で算出される。今はc-statistics, calibration slope, intercept, brier scoreを算出しているが、内部コードをいじれば簡単にECIなども算出することはできる。
boot.bias.corr.estimation(orig.data=d,var=model1,B=500,method="MLE")
また、同時にoptimismの全データとヒストグラムも表示されるようにしたので、あとからbias-corrected estimateの95%信頼区間を算出することもできる。
色々Githubとか調べてみたけどニーズにハマるパッケージが見つけられず数日浪費していたが、最終的に自作するのが一番だと気づいた。これで予測モデル研究がだいぶはかどる。 もう少し頑張ればBootstrapしたCalibration plotの作成もできそうだけど、ひとまずはここまで。
