1. Install the relevant packages

Hide code cell source
   suppressWarnings({
     invisible({
       rm(list = ls())

       if (!require(icdpicr, quietly = TRUE, warn.conflicts = FALSE)) install.packages('icdpicr')
       if (!require(dplyr, quietly = TRUE, warn.conflicts = FALSE)) install.packages('dplyr')
       if (!require(readr, quietly = TRUE, warn.conflicts = FALSE)) install.packages('readr')
       if (!require(tidyr, quietly = TRUE, warn.conflicts = FALSE)) install.packages('tidyr')

       library(icdpicr, quietly = TRUE, warn.conflicts = FALSE)
       library(dplyr, quietly = TRUE, warn.conflicts = FALSE)
       library(readr, quietly = TRUE, warn.conflicts = FALSE)
       library(tidyr, quietly = TRUE, warn.conflicts = FALSE)
     })
   })
  1. Get some documentation

   ls("package:icdpicr")
  1. 'cat_trauma'
  2. 'injury'
  1. Study the syntax

   print(cat_trauma) 
Hide code cell output
function (df, dx_pre, icd10, i10_iss_method, calc_method = 1, 
    verbose = FALSE) 
{
    if (!is.data.frame(df)) 
        stop("First argument must be a dataframe")
    if (NROW(df) == 0) 
        stop("Data contains no observations. It must contain at least one row.")
    if (!is.character(dx_pre)) 
        stop("Second argument must be a character string")
    if (make.names(dx_pre) != dx_pre) 
        stop("Second argument must be a valid variable name in R")
    if (!(calc_method %in% c(1, 2))) 
        stop("calc_method must be either 1 or 2")
    if (!(icd10 %in% c(TRUE, FALSE, "cm", "base"))) 
        stop("icd10 must be TRUE, FALSE, 'cm', or 'base'")
    if (icd10 == FALSE) 
        i10_iss_method <- ""
    if (i10_iss_method == "roc_max") 
        stop("The roc_max option has been depricated. Please use roc_max_NIS, roc_max_TQIP, roc_max_NIS_only, or roc_max_TQIP_only instead.")
    if ((icd10 != FALSE) && !(i10_iss_method %in% c("roc_max_NIS", 
        "roc_max_TQIP", "roc_max_NIS_only", "roc_max_TQIP_only", 
        "gem_max", "gem_min"))) 
        stop("i10_iss_menthod must be roc_max_NIS, roc_max_TQIP, roc_max_NIS_only, roc_max_TQIP_only, gem_max, or gem_min.")
    regex_dx <- paste0("^", dx_pre, "([0-9]+)$")
    dx_colnames <- grep(regex_dx, names(df), value = TRUE)
    dx_nums <- as.numeric(sub(regex_dx, "\\1", dx_colnames))
    num_dx <- length(dx_nums)
    if (num_dx == 0) 
        stop("No variables with prefix found in data")
    df <- data.frame(df)
    if (isTRUE(icd10)) 
        icd10 <- "cm"
    if (icd10 %in% c("base", "cm")) {
        etab <- rbind(etab_s1, i10_ecode)
        ntab <- switch(i10_iss_method, roc_max_NIS = rbind(ntab_s1, 
            .select_i10_data("NIS", icd10)), roc_max_TQIP = rbind(ntab_s1, 
            .select_i10_data("TQIP", icd10)), roc_max_NIS_only = rbind(ntab_s1, 
            .select_i10_data("NIS_only", icd10)), roc_max_TQIP_only = rbind(ntab_s1, 
            .select_i10_data("TQIP_only", icd10)), gem_max = rbind(ntab_s1, 
            i10_map_max), gem_min = rbind(ntab_s1, i10_map_min))
    }
    else {
        ntab <- ntab_s1
        etab <- etab_s1
    }
    for (i in dx_nums) {
        dx_name <- paste0(dx_pre, i)
        df_ss <- df[, dx_name, drop = FALSE]
        df_ss$n <- 1:NROW(df_ss)
        df_ss[, dx_name] <- sub("\\.", "", df_ss[, dx_name])
        if (icd10 == TRUE & i10_iss_method == "roc_max") {
            i9_valid <- c("8", "9", "E")
            i10_valid <- c("S", "T", "U", "V", "W", "X", "Y")
            df_ss[, dx_name] <- ifelse(substr(df_ss[, dx_name], 
                1, 1) %in% c(i9_valid, i10_valid), df_ss[, dx_name], 
                NA)
            process_i10 <- function(s) {
                stopifnot(is.character(s) | is.na(s))
                ret_val <- NA
                s <- sub("\\.", "", s)
                if (!substr(s, 1, 1) %in% c("S", "T", "U", "V", 
                  "W", "X", "Y")) {
                  ret_val <- s
                }
                else if (nchar(s) < 7 & !grepl("X", substr(s, 
                  2, nchar(s)))) {
                  ret_val <- s
                }
                else if (nchar(s) != 7) {
                  ret_val <- ""
                }
                else if (substr(s, 7, 7) != "A") {
                  ret_val <- ""
                }
                else if (substr(s, 5, 5) == "X") {
                  ret_val <- substr(s, 1, 4)
                }
                else if (substr(s, 6, 6) == "X") {
                  ret_val <- substr(s, 1, 5)
                }
                else {
                  ret_val <- substr(s, 1, 6)
                }
                return(ret_val)
            }
            df_ss[, dx_name] <- sapply(df_ss[, dx_name], process_i10)
        }
        temp <- merge(df_ss, ntab, by.x = dx_name, by.y = "dx", 
            all.x = TRUE, all.y = FALSE, sort = FALSE)
        temp <- temp[order(temp$n), ]
        temp <- temp[, c("severity", "issbr")]
        if (calc_method == 2) {
            temp[which(temp$severity == 6), "severity"] <- 5
        }
        names(temp) <- paste0(c("sev_", "issbr_"), i)
        df <- .insert_columns(df, dx_name, temp)
    }
    body_regions <- unique(i10_map_max$issbr)
    issbr_names <- gsub("/", "_", body_regions)
    for (i in body_regions) {
        temp <- df[, grepl("sev_", names(df)), drop = FALSE] * 
            (1 * (df[, grepl("issbr_", names(df))] == i))
        df[, paste0("mxaisbr_", gsub("/", "", i))] <- apply(temp, 
            1, function(row) {
                row <- ifelse(row == 0, NA, row)
                if (all(is.na(row))) {
                  maxaisbr <- 0
                }
                else if (all(row == 9, na.rm = TRUE)) {
                  maxaisbr <- 9
                }
                else {
                  maxaisbr <- max(c(0, row[row != 9]), na.rm = TRUE)
                }
                return(maxaisbr)
            })
    }
    c9to0 <- function(x) ifelse(x == 9, 0, x)
    df$maxais <- apply(df, 1, function(row) {
        row <- row[grepl("mxaisbr", names(row))]
        if (all(is.na(row))) {
            maxais <- as.numeric(NA)
        }
        else if (max(c9to0(row), na.rm = TRUE) == 0) {
            maxais <- max(row, na.rm = TRUE)
        }
        else {
            maxais <- max(c9to0(row), na.rm = TRUE)
        }
        return(maxais)
    })
    df$maxais <- as.numeric(df$maxais)
    df$riss <- apply(df, 1, function(row) {
        temp <- row[grepl("^mxaisbr", names(row))]
        temp <- as.numeric(c9to0(temp))
        sum(temp[order(-temp)[1:3]]^2)
    })
    df[df$maxais == 6, "riss"] <- 75
    df[df$maxais == 9, "riss"] <- NA
    df$niss <- apply(df, 1, function(row) {
        temp <- row[grepl("^sev_", names(row))]
        temp <- as.numeric(temp)
        temp <- ifelse(is.na(temp) | temp == 9, 0, temp)
        sum(temp[order(-temp)[1:3]]^2)
    })
    df[df$maxais == 6, "niss"] <- 75
    df[df$maxais == 9, "niss"] <- NA
    ecode_colnames <- paste0("ecode_", 1:4)
    df[, ecode_colnames] <- NA
    ecode_regex <- paste0("^", etab$dx, collapse = "|")
    df[, ecode_colnames] <- t(apply(df, 1, function(row) {
        row <- sub("\\.", "", row)
        row_ecodes <- stringr::str_extract(as.character(unlist(row)), 
            ecode_regex)
        row_ecodes <- na.omit(row_ecodes)
        row_ecodes[1:4]
    }))
    for (i in 1:4) {
        col_name <- paste("ecode_", i, sep = "")
        df_ss <- df[, col_name, drop = FALSE]
        df_ss$n <- 1:NROW(df_ss)
        df_ss[, col_name] <- sub("\\.", "", df_ss[, col_name])
        temp <- merge(df_ss, etab, by.x = col_name, by.y = "dx", 
            all.x = TRUE, all.y = FALSE, sort = FALSE)
        temp <- temp[order(temp$n), ]
        temp <- temp[, c("mechmaj", "mechmin", "intent")]
        names(temp) <- paste(c("mechmaj", "mechmin", "intent"), 
            i, sep = "")
        df <- .insert_columns(df, col_name, temp)
    }
    if (stringr::str_detect(i10_iss_method, "NIS|TQIP") && icd10 %in% 
        c("cm", "base")) {
        if (verbose) 
            print("Calculating mortality prediction")
        coef_df <- .select_i10_coef(prefix = stringr::str_extract(i10_iss_method, 
            "NIS|TQIP"), icd10)
        stopifnot(max(coef_df$intercept, na.rm = TRUE) == min(coef_df$intercept, 
            na.rm = TRUE))
        intercept <- max(coef_df$intercept, na.rm = TRUE)
        coef_df <- coef_df[!is.na(coef_df$effect), ]
        effect_hash <- coef_df$effect
        names(effect_hash) <- coef_df$dx
        calc_mortality_prediction <- function(dx) {
            x <- sum(effect_hash[sub("\\.", "", dx)], na.rm = TRUE) + 
                intercept
            1/(1 + exp(-x))
        }
        mat <- as.matrix(df[, grepl(paste0("^", dx_pre), names(df))])
        df$Pmort <- apply(mat, 1, calc_mortality_prediction)
    }
    rownames(df) <- 1:nrow(df)
    df
}
<bytecode: 0x7fc620b3d770>
<environment: namespace:icdpicr>
  1. What is the value of dx_pre?

  print(injury)
Hide code cell output
# A tibble: 100,477 × 11
   dx1      dx2      dx3      dx4      dx5   dx6   dx7   dx8   dx9   dx10   died
   <chr>    <chr>    <chr>    <chr>    <chr> <chr> <chr> <chr> <chr> <chr> <int>
 1 S72.342A NA       NA       NA       NA    NA    NA    NA    NA    NA        0
 2 S05.22XA NA       NA       NA       NA    NA    NA    NA    NA    NA        0
 3 S00.01XA S00.03XA S00.11XA S00.12XA S00.… S00.… S00.… S01.… S02.… S80.…     0
 4 S21.119A NA       NA       NA       NA    NA    NA    NA    NA    NA        0
 5 S82.191A NA       NA       NA       NA    NA    NA    NA    NA    NA        0
 6 S22.42XA NA       NA       NA       NA    NA    NA    NA    NA    NA        0
 7 S92.052A S92.065A S92.325A S92.335A S92.… S93.… NA    NA    NA    NA        0
 8 S02.112A S06.5X0A S12.090A S12.100A S12.… S20.… S20.… S22.… S22.… S22.…     0
 9 S00.03XA S22.058A S22.068A S22.078A S22.… S30.… S42.… S62.… S81.… S82.…     0
10 S61.411A S62.624B S62.626B S66.300A S66.… S66.… S66.… NA    NA    NA        0
# ℹ 100,467 more rows
  1. As such, the script should first use a sample:

   # use subset of injury given its size
   inj = injury[1:100,1:3]
   df_score = cat_trauma(df=inj, dx_pre="dx", icd10=TRUE, i10_iss_method="roc_max_NIS", calc_method=1, verbose=FALSE)

   # visualize the output
   df_score[1:9,1:9]; df_score[1:9,10:14]; df_score[1:9,15:18]; df_score[1:9,34:35]
Hide code cell output
A data.frame: 9 × 9
dx1sev_1issbr_1dx2sev_2issbr_2dx3sev_3issbr_3
<chr><int><chr><chr><int><chr><chr><int><chr>
1S72.342A1ExtremitiesNA NANA NA NANA
2S05.22XA1Face NA NANA NA NANA
3S00.01XA1Head/Neck S00.03XA 2Head/Neck S00.11XA 2Face
4S21.119A3Chest NA NANA NA NANA
5S82.191A3ExtremitiesNA NANA NA NANA
6S22.42XA3Chest NA NANA NA NANA
7S92.052A1ExtremitiesS92.065A 1ExtremitiesS92.325A 1Extremities
8S02.112A1Face S06.5X0A 4Head/Neck S12.090A 4Head/Neck
9S00.03XA2Head/Neck S22.058A 1Chest S22.068A 3Chest
A data.frame: 9 × 5
mxaisbr_Generalmxaisbr_HeadNeckmxaisbr_Facemxaisbr_Extremitiesmxaisbr_Chest
<dbl><dbl><dbl><dbl><dbl>
100010
200100
302200
400003
500030
600003
700010
804100
902003
A data.frame: 9 × 4
mxaisbr_Abdomenmaxaisrissniss
<dbl><dbl><dbl><dbl>
101 1 1
201 1 1
302 8 9
403 9 9
503 9 9
603 9 9
701 1 3
8041733
9031314
A data.frame: 9 × 2
intent4Pmort
<chr><dbl>
1NA0.013857916
2NA0.014412795
3NA0.017673701
4NA0.026078454
5NA0.027839150
6NA0.024364215
7NA0.005516156
8NA0.037222518
9NA0.021037152
  1. What is the value of the ISS?

  print(df_score$iss)
NULL
  1. What is the value of the NISS?

   print(df_score$niss)
Hide code cell output
  [1]  1  1  9  9  9  9  3 33 14  6 22  1  1  2  1  4  1  6  1  6  1  9  1  3 22
 [26]  6 10  6  1 30  4  9 11 33 17  1 17  1 17  3  2  6  8  9  5  4  1 21  1 11
 [51] 20 11  6  1 26  9  3  2 27 11  6 19  1 11  6  3  3 41  2 10  9 11  1  1  1
 [76]  1 18 22 10  3 19  6  1  9  9 14  9  5  4 33 21 22  1  4  2  1 13  1  9 35
  1. What is the value of the RTS?

   print(df_score$rts)
NULL

\(\vdots\)

  1. Edit code to run on the full dataset

    # use subset of injury given its size (replace "injury" with your datafile of interest)
    df_score = cat_trauma(df=injury, dx_pre="dx", icd10=TRUE, i10_iss_method="roc_max_NIS", calc_method=1, verbose=FALSE)
    
    # visualize the output
    df_score[1:9,1:9]; df_score[1:9,10:14]; df_score[1:9,15:18]; df_score[1:9,34:35]