optimal_k: Compute approximate optimal k

In vlad: Variable Life Adjusted Display and Other Risk-Adjusted Quality Control Charts

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.

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

Value

Returns a single value which is the approximate optimal k.

Author(s)

Philipp Wittenberg

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))

## End(Not run)

vlad documentation built on Feb. 15, 2021, 5:12 p.m.