评估R中已发布和新的风险评分的校准然后绘制校准曲线

如何解决评估R中已发布和新的风险评分的校准然后绘制校准曲线

我目前正在预测疾病状况时对两个风险评分的校准（或错误校准程度）（二进制结果：1 =患有疾病，0 =无疾病）。一个是基于逻辑回归创建的，而另一个（称为弗雷明汉CVD风险评分）是使用cox比例风险模型创建的。

我在文献中看到了两个示例，一种方法涉及仅使用线性预测变量和intecept进行logistic回归估算预测概率：

例如glm(disease_status ~ -1 + offset(intercept_and_linear_predictor),data = df,family = 'binomial'))

方法2涉及使用绝对风险预测（以前的研究中未提供代码，因此下面的代码反映了我认为他们根据文章所进行的操作）。

我的问题是，当使用不同方法（一种逻辑回归；另一种Cox比例风险模型）计算风险得分时，如何使用方法2评估校准？

这是我到目前为止的内容： 加载软件包：

#load pROC
library(tidyverse)
library(rms)
library(Hmisc)
library(knitr)
library(broom)
library(pander)
library(ggbeeswarm)
library(gridExtra)
library(grid)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(viridis)
library(CalibrationCurves)

生成数据：

#generate df with random numbers
set.seed(123)
df <- data.frame(disease_status = rbinom(n=100,size=1,prob=0.20),sex = rbinom(n=100,prob=0.50),years_to_diagnosis = rnorm(100,mean=3.2,sd=1),new_risk_score_linear_predictor = rnorm(100,mean=2,sd=2),FRS_linear_predictor_males =rnorm(100,mean=6,FRS_linear_predictor_females =rnorm(100,mean=5,sd=1.2))

方法1

# based on https://darrendahly.github.io/post/homr/
# Add intercept to linear predictor
df$new_score_intercept <- 3.76 + df$new_risk_score

# fit model with disease as outcome
new_score_orig <-     glm(disease_status ~ -1 + offset(new_score_intercept),family = 'binomial')

方法2

# Calculating absolute risk 
# FRS CVD risk 
df$FRS_absolute_risk <-   exp(df$FRS_linear_predictor_males - 23.9802)
df$FRS_final_females   <-   1 -  0.88936^exp(df$FRS_absolute_risk)

df$FRS_absolute_risk <-   exp(df$FRS_linear_predictor_females - 26.1931)
df$FRS_final_females   <-   1 -  0.95012^exp(df$FRS_absolute_risk)

df$FRS_absolute_risk_combined <- ifelse(df$sex==1,df$FRS_final_males,df$FRS_final_females)
summary(df$FRS_absolute_risk_combined)

df$new_risk_score_absolute_risk       <- (exp(df$new_score_intercept)) / (1 + exp(df$new_score_intercept))

new_score_orig <-     glm(disease_status ~ -1 + offset(new_risk_score_absolute_risk),family = 'binomial')
FRS_orig <-  glm(disease_status ~ -1 + offset(df$FRS_absolute_risk),family = 'binomial')

用于校准的循环

models <- c('new_score_orig','FRS_orig')
for(x in models){
 model <- get(x)
 # get predicted probabilities from each model and create a new variable in test dataset
 df[paste(x,'pred',sep="_")] <- predict(model,type = "response")
 log_pred <- as.data.frame(predict(model,type='response'))
 assign(paste('log_pred',x,sep="_"),log_pred)
 cat('Brier score: log_pred',val.prob(log_pred[,1],df$disease_status)['Brier'],'\n') 
 x <- val.prob.ci.2(log_pred[,df$disease_status,smooth = "rcs",CL.smooth = "fill",logistic.cal=FALSE) 
 print('----------------------------------------------------------------')
}

真的很感谢您提出的任何建议-我完全感到困惑！