问题描述
我目前正在预测疾病状况时对两个风险评分的校准（或错误校准程度）（二进制结果：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('----------------------------------------------------------------')
}
真的很感谢您提出的任何建议-我完全感到困惑！
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）