The R code for the core analyses is included in the file below.
r code for adaptation.doc
# 7 predictor model in individual patient data full <- lrm(STATUS~SEX+AGE10+MI+CHF+ISCHEMIA+LUNG+RENAL,data=AAA,x=T,y=T) # Estimate shrinkage factor full.validate <- validate(full, B=200) # Apply shrinkage factor on model coefficients full.shrunk.coefficients <- full.validate["Slope","index.corrected"] * full$coefficients # Adjust intercept; alternatively, use as offset variable full.shrunk.coefficients[1] <- full.shrunk.coefficients[1] + full.validate["Intercept","index.corrected"] # Full model, penalized estimation penalty <- pentrace(full,penalty=c(0.5,1,2,3,4,6,8,12,16,24),maxit=25) full.penalized <- update(full, penalty=penalty$penalty) # Adaptation method; start with univariate literature coefficients # analyzed with random effects meta-analysis lit.coefficients <- c(0.3606,0.788,1.034,1.590,1.514,0.8502,1.302) lit.se <- c(0.176,0.112,0.317,0.4109,0.378,0.2367,0.2495) # Univariate coefs in individual patient data, followed by simple adaptation method fit.ind.x <- rep(0,7) for (i in 1:7) {fit.ind.x <- lrm.fit(y=full$y,x=full$x[,i]) # Adaptation method 1, identical to proposal Greenland adapt.coef[i] <- full$coefficients[i+1] + 1 * (lit.coefficients[i] - fit.ind.x$coefficients[2]) adapt.var[i] <- full$var[i+1,i+1] - fit.ind.x$var[2,2] + lit.se[i]^2 } # We note that the variance has decreased considerably sqrt(adapt.var) / sqrt(diag(full$var)[2:8]) # Round the coefficients after shrinkage 0.9 cat("\n",round(9*adapt.coef,0)) # we estimate the intercept, first scale age around 70 years (7 decades) X <- as.matrix(full$x) X[,2] <- X[,2] – 7 # age rescaled # offset with scores / 10 offset.AAA <- X %*% c(.3,.3,.2,.8,.7,.6,1.1) fit.offset <- lrm.fit(y=full$y, offset=offset.AAA) # intercept given these scores for predictors fit.offset$coef[1] # result: -3.48 fit.offset$coef[1] + log((0.05/0.95)/(18/220)) # result: -3.92 # there is a function ‘bootcor.uni.mult’ which performs bootstrapping to obtain the correlation between univariate and multivariable logistic regression coefficients in the individual patient data; this correlation is required for Adaptation method 2 # Obtain correlations uni – mult, and shrinkage factor full.uni.mult.cor <- bootcor.uni.mult(full,B=200,maxit=10,trim=0.1,save.indices = T) # Adaptation method 2, with c optimal for (i in 1:7) { fit.ind.x <- lrm.fit(y=full$y,x=full$x[,i],maxit=25) adapt.factors[i] <- (full.uni.mult.cor$r[i] * sqrt(full$var[i+1,i+1]) * sqrt(fit.ind.x$var[2,2]) ) / (lit.se[i]^2 + fit.ind.x$var[2,2]) adapt.coefficients[i] <- full$coefficients[i+1] + adapt.factors[i] * (lit.coefficients[i] - fit.ind.x$coefficients[2]) adapt.var[i] <- full$var[i+1,i+1] * (1-(full.uni.mult.cor$r[i]^2 * fit.ind.x$var[2,2] / (lit.se[i]^2 + fit.ind.x$var[2,2]))) } # end loop over 7 predictors ## Recalibrate literature coefficients with one calibration factor ## # make offset based on univariate literature coefficients offset.AAA2 <- X %*% lit.coefficients fit.offset2 <- lrm.fit(y=full$y, x=offset.AAA2) fit.offset2$stats round(fit.offset2$coef[2],2) # result: 0.69 # scores with this overall recalibration cat("\n",round(10*fit.offset2$coef[2]*lit.coefficients,0),"\n") # use these scores / 10 as offset for intercept offset.AAA3 <- X %*% c(.3,.5,.7,1.1,1.0,.6,.9) fit.offset3 <- lrm.fit(y=full$y, offset=offset.AAA3) fit.offset3$coef[1] # -4.05 fit.offset3$coef[1] + log((0.05/0.95)/(18/220)) # -4.49