In ampel-leipzig/ameld: Data and Model of End-Stage Liver Disease used in the AMPEL Project

Sys.setenv(LANGUAGE = "en")
library("ameld")
library("survival")

data(eldd)
data(eldr)

Authors: r packageDescription("ameld")[["Author"]]
Last modified: r file.info("ameld.Rmd")$mtime
Compiled: r date()

The ameld Package

The ameld R package provides a dataset, eldd, of patients evaluated for liver transplantation at the University Hospital of Leipzig from r paste(strftime(attr(eldd, "period"), "%B %Y"), collapse = " to "). eldr contains the reference limits for the laboratory measurements used in eldd.

Access the Data

The datasets could be loaded as follows:

library("ameld")
data(eldd)
data(eldr)

colnames(eldd)
knitr::kable(head(eldd[, 1:10]))
knitr::kable(head(eldr))

Survival Estimates

The ameld package provides a plot_surv function that is similar to survival::plot.survfit() but uses different defaults:

library("survival")

srv <- Surv(eldd$DaysAtRisk, eldd$Deceased)

srvfit <- survfit(srv ~ 1)

plot_surv(
    srvfit,
    main = "Kaplan-Meier Survival Estimate"
)

plot_surv(
    srvfit,
    cumhaz = TRUE,
    main = "Cumulative Hazard"
)

A more complex example of the survival plot with risk tables could be generated by combining plot_surv and plot_table:

## timepoints of interest
times <- c(0, 30, 90)

## calculate risk tables
sm <- summary(srvfit, times = times)
nrisk <- as.matrix(sm$n.risk)
ncumevents <- as.matrix(cumsum(sm$n.event))
rownames(nrisk) <- rownames(ncumevents) <- times

## keep old graphic parameters and restore them afterwards
old.par <- par(no.readonly = TRUE)

layout(matrix(1:3, nrow = 3), height = c(5, 1, 1))
par(cex.main = 2)
plot_surv(
    srvfit,
    main = "Kaplan-Meier Survival Estimate to day 90",
    times = times,
    xmax = 90
)
par(mar = c(5.1, 4.1, 1.1, 2.1))
plot_table(
    nrisk, at = times, main = "Number at risk",
    xaxis = FALSE, cex.text = 1.5, ylabels = FALSE
)
par(mar = c(5.1, 4.1, 1.1, 2.1))
plot_table(
    ncumevents, at = times, main = "Cumulative number of events",
    xaxis = FALSE, cex.text = 1.5, ylabels = FALSE
)

par(old.par)

Z(log) Transformation

The $z(log)$-transformation was suggested in @hoffmann2017. It is similar to common $z$-transformation but standardizes laboratory measurements by their respective reference or normal values. An R implementation is provided by the zlog package [@zlog].

library("zlog")

## transform reference data.frame for zlog
r <- eldr[c("Code", "AgeDays", "Sex", "LowerLimit", "UpperLimit")]
names(r) <- c("param", "age", "sex", "lower", "upper")
r$age <- r$age / 365.25
r <- set_missing_limits(r)

## we just want to standardize laboratory values
cn <- colnames(eldd)
cnlabs <- cn[grepl("_[SCEFQ1]$", cn)]
zeldd <- eldd
zeldd[c("Age", "Sex", cnlabs)] <- zlog_df(eldd[, c("Age", "Sex", cnlabs)], r)

We could use the (mean) $z(log)$ values to get a first impression of the influence of these values on survival/non-survival:

## divide data.frame by dead/alive
s <- split(zeldd[cnlabs], zeldd$Deceased)
names(s) <- c("survived", "dead")

## calculate mean standardized lab values
s <- lapply(s, colMeans, na.rm = TRUE)
o <- order(s$dead)

## comparison plot
col <- palette.colors(2)

## keep old graphic parameters and restore them afterwards
old.par <- par(no.readonly = TRUE)

par(mar = c(7.1, 4.1, 4.1, 2.1))

plot(
    s$dead[o], type = "b", pch = 20, lwd = 2, col = col[1],
    axes = FALSE, ann = FALSE
)
lines(s$survived[o], type = "b", pch = 20, lwd = 2, col = col[2])
legend(
    "bottomright",
    legend = c("dead", "survived"),
    col = col, lwd = 2, pch = 20, bty = "n"
)
title(
    main = "Mortality Status vs Mean Standardized Laboratory Values", adj = 0
)
title(xlab = "Laboratory Measurements", adj = 1L, line = 5)
title(ylab = "Mean Standardized Values", adj = 1L)
r <- range(unlist(s))
axis(
    2, at = seq(from = floor(r[1]), to = ceiling(r[2])),
    lwd.ticks = 0, col = "#808080"
)
axis(
    1, at = seq_along(o), labels = names(s$dead[o]), las = 2,
    lwd.ticks = 0L, col = "#808080"
)
par(old.par)

Observed vs Expected Mortality

We divide our cohort into 5 MELD-score risk categories as described in @wiesner2003. Our observed 90 day mortality is higher than the predicted one, except in the lowest category.

eldd$MELD <- meld(
    creatinine = as_metric(eldd$CRE_S, "creatinine"),
    bilirubin = as_metric(eldd$BILI_S, "bilirubin"),
    inr = eldd$INR_C,
    dialysis = eldd$Dialysis,
    cause = "other"
)

# Mortality rates and categories as reported in Wiesner et al. 2003
mr <- c(1.9, 6.0, 19.6, 52.6, 71.3) / 100

mcat <- cut(
    eldd$MELD,
    breaks = c(-Inf, seq(10, 40, by=10), Inf),
    labels = c(
        paste0("[", floor(min(eldd$MELD, na.rm = TRUE)), ",9]"),
        "[10,20)", "[20,30)", "[30,40)",
        paste0( "[40,", ceiling(max(eldd$MELD, na.rm = TRUE)), ")")
    ),
    right = FALSE
)

tbl <- observed_vs_expected_mortality(srv, time = 90, f = mcat, expected = mr)
tbl[c("ObservedMortality", "ExpectedMortality")] <-
    tbl[c("ObservedMortality", "ExpectedMortality")] * 100

knitr::kable(
    tbl,
    col.names =
        c(
            "Observed deaths (n)",
            "Expected deaths (n)",
            "Standardized mortality ratio (SMR)",
            "Observed mortality (%)",
            "Expected mortality (%)"
        ),
    caption = paste0(
        "Observed vs MELD-expected 90 day mortality. ",
        "MELD mortality values are taken from @wiesner2003. ",
        "All patients censored before day 90 are ignored for ",
        " the calculation of the MELD-expected deaths. ",
        "SMR, Standardized mortality ratio = observed deaths/expected deaths."
    ),
    digits = 1
)

Acknowledgment

This work is part of the AMPEL (Analysis and Reporting System for the Improvement of Patient Safety through Real-Time Integration of Laboratory Findings) project.

This measure is co-funded with tax revenues based on the budget adopted by the members of the Saxon State Parliament.

Session Information

sessionInfo()

References

ampel-leipzig/ameld documentation built on Aug. 23, 2024, 7:31 p.m.