vlad (version 0.2.0)

optimal_k: Compute optimal k

Description

Compute optimal k.

Usage

optimal_k(QA, df, coeff, yemp = TRUE)

Arguments

QA

Double. Defines the performance of a surgeon with the odds ratio ratio of death Q.

df

Data Frame. First column Parsonnet Score and second column outcome of each operation.

coeff

Numeric Vector. Estimated coefficients \(\alpha\) and \(\beta\) from the binary logistic regression model. For more information see details.

yemp

Logical. If TRUE use observed outcome value, if FALSE use estimated binary logistc regression model.

Value

Returns a single value which is the approximate optimal k for a set of given Parsonnet scores.

Details

Formula deterioration: $$ k{det} = \frac{Q{A} - 1 - log(Q{A})}{log(Q{A})}\bar{p} , Q{A} > 1 $$ Formula improvement: $$ k{imp} = \frac{1 - Q{A} + log(Q{A})}{log(Q{A})}\bar{p} , Q{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("dplyr")
library(vlad)
data("cardiacsurgery", package = "spcadjust")

## preprocess data to 30 day mortality and subset phase I (In-control) of surgeons 2
S2I <- cardiacsurgery %>% rename(s = Parsonnet) %>%
  mutate(y = ifelse(status == 1 & time <= 30, 1, 0),
        phase = factor(ifelse(date < 2*365, "I", "II"))) %>%
  filter(phase == "I", surgeon == 2) %>% select(s, y)

coeff <- coef(glm(y ~ s, data = S2I, family = "binomial"))

## (Deterioration)
optimal_k(QA = 2, df = S2I, coeff = coeff, yemp = FALSE)

## manually find optimal k for detecting deterioration
QA <- 2
pbar <- mean(sapply(S2I[, 1], gettherisk, coef = coeff))
kopt <- pbar * ( QA - 1 - log(QA) ) / log(QA)

all.equal(kopt, optimal_k(QA = 2, df = S2I, coeff = coeff, yemp = FALSE))

## (Improvement)
optimal_k(QA = 1/2, df = S2I, coeff = coeff, yemp = FALSE)

## manually find optimal k for detecting improvement
QA <- 1/2
pbar <- mean(sapply(S2I[, 1], gettherisk, coef = coeff))
kopt <- pbar * ( 1 - QA + log(QA) ) / log(QA)

all.equal(kopt, optimal_k(QA = 1/2, df = S2I, coeff = coeff, yemp = FALSE))
# }