vlad (version 0.2.2)

optimal_k: Compute approximate optimal k

Description

Compute approximate optimal k.

Usage

optimal_k(pmix, RA, yemp = FALSE)

Arguments

pmix

Data Frame. A three column data frame. First column is the operation outcome. Second column are the predicted probabilities from the risk model. Third column can be either the predicted probabilities from the risk model or average outcome.

RA

Double. Odds ratio of death under the alternative hypotheses. Detecting deterioration in performance with increased mortality risk by doubling the odds Ratio RA = 2. Detecting improvement in performance with decreased mortality risk by halving the odds ratio of death RA = 1/2. Odds ratio of death under the null hypotheses is 1.

yemp

Logical. If TRUE, use emirical outcome values, else use model.

Value

Returns a single value which is the approximate optimal k.

Details

Formula deterioration: $$ k{det} = \frac{R{A} - 1 - log(R{A})}{log(R{A})}\bar{p} , R{A} > 1 $$ Formula improvement: $$ k{imp} = \frac{1 - R{A} + log(R{A})}{log(R{A})}\bar{p} , R{A} < 1 $$

References

Wittenberg P, Gan FF, Knoth S (2018). A simple signaling rule for variable life-adjusted display derived from an equivalent risk-adjusted CUSUM chart. Statistics in Medicine, 37(16), pp 2455--2473.

Examples

# NOT RUN {
library(vlad)
library(dplyr)
data("cardiacsurgery", package = "spcadjust")

## preprocess data to 30 day mortality
SALL <- cardiacsurgery %>% rename(s = Parsonnet) %>%
  mutate(y = ifelse(status == 1 & time <= 30, 1, 0),
         phase = factor(ifelse(date < 2*365, "I", "II")))
SI <- subset(SALL, phase == "I")
GLM <- glm(y ~ s, data = SI, family = "binomial")
pi1 <- predict(GLM, type = "response", newdata = data.frame(s = SI$s))
pmix <- data.frame(SI$y, pi1, pi1)

## (Deterioration)
optimal_k(pmix = pmix, RA = 2)

## manually find optimal k for detecting deterioration
RA <- 2
pbar <- mean(pmix$pi1)
kopt <- pbar * ( RA - 1 - log(RA) ) / log(RA)

all.equal(kopt, optimal_k(pmix = pmix, RA = 2))

## (Improvement)
optimal_k(pmix = pmix, RA = 1/2)

## manually find optimal k for detecting improvement
RA <- 1/2
pbar <- mean(pmix$pi1)
kopt <- pbar * ( 1 - RA + log(RA) ) / log(RA)

all.equal(kopt, optimal_k(pmix = pmix, RA = 1/2))
# }