interceptAdj.R
In predtools: Prediction Model Tools

## ---- include = FALSE---------------------------------------------------------

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)


## ----setup--------------------------------------------------------------------

library(predtools)
library(magrittr)
library(dplyr)
library(ggplot2)


## -----------------------------------------------------------------------------

data(dev_data)
data(val_data)


## ----echo=FALSE---------------------------------------------------------------
set.seed(1)
val_data$y_alt <- ifelse(val_data$y == 0, 0,
                         ifelse(runif(n = nrow(val_data)) <= 0.5, 0, 1))

val_data %>%
  select(y, y_alt) %>%
  summary() %>%
  knitr::kable()


## ----echo=FALSE---------------------------------------------------------------

knitr::kable(dev_data[1:7,])


## -----------------------------------------------------------------------------

reg<-glm(y~sex+age+severity+comorbidity,data=dev_data,family=binomial(link="logit"))
summary(reg)


## ----echo=FALSE---------------------------------------------------------------

cfs <- coefficients(reg)
str<-paste0(round(cfs[1],4),"+",paste0(round(cfs[-1],4),"*",names(cfs[-1]),collapse="+"))
str_risk_model <- gsub("+-", "-", str, fixed = T)


## -----------------------------------------------------------------------------

dev_data$pred <- predict.glm(reg, type = 'response')
val_data$pred <- predict.glm(reg, newdata = val_data, type = 'response')

calibration_plot(data = dev_data, obs = "y", pred = "pred", title = "Calibration plot for development data")
calibration_plot(data = val_data, obs = "y_alt", pred = "pred", y_lim = c(0, 0.6),
                 title = "Calibration plot for validation data")


## -----------------------------------------------------------------------------

odds_correction_factor <- odds_adjust(p0 = mean(dev_data$y), p1 = mean(val_data$y_alt), v = var(dev_data$pred))
odds_correction_factor


## -----------------------------------------------------------------------------

dev_data$pred <- predict.glm(reg, type = 'response')
val_data$pred <- predict.glm(reg, newdata = val_data, type = 'response')

val_data$odds_adj <- (val_data$pred / (1 - val_data$pred)) * odds_correction_factor
val_data$pred_adj <- val_data$odds_adj / (1 + val_data$odds_adj)

val_data$id <- c(1 : nrow(val_data))
val_data_long <- reshape(data = val_data, direction = "long", varying = c("pred", "pred_adj"), v.name = "preds",
                         idvar = "id", timevar = "Method", times = c("Primitive", "Adjusted"))

calibration_plot(data = val_data, obs = "y_alt", pred = "pred_adj",
                 title = "Calibration plot for development data - after recalibration")
calibration_plot(data = val_data_long, obs = "y_alt", pred = "preds", group = "Method",
                 title = "Calibration plot for development data - before and after recalibration")