---
title: Baseline Creatinine Determination Method Impacts Association Between Acute Kidney Injury and Clinical Outcomes
author: 


output: 
  officedown::rdocx_document:
    reference_docx: word-style-reference-manuscript.docx
bibliography: AKI_CKD.bib
csl: american-medical-association-brackets.csl

header-includes:
- \usepackage{placeins}

---
```r
library(officedown)
```


W. Alton Russell, MS^1^, David Scheinker, PhD\*^1,2,3^, Scott M. Sutherland, MD\*^2,3^


<br>

^1^Department of Management Science and Engineering, Stanford University, Stanford, CA, USA

^2^Stanford School of Medicine, Stanford, CA, USA

^3^Lucile Packard Children’s Hospital Stanford, Palo Alto, CA, USA

\* Co-senior authors

<br>

# Corresponding author

Scott M. Sutherland, MD 

300 Pasteur Drive, Room G-306

Stanford, CA 94304

650-387-5286

suthersm@stanford.edu

#####

```{r setup, include = FALSE}
# library(officedown)
library(ggplot2)
library(GGally)
library(data.table)
library(pROC)
library(gridExtra)
library(flextable)
library(officer)
theme_set(theme_bw())


run_ROC = FALSE

prop.with.CI <- function(num, denom, alpha = 0.05){
  prop <- num / denom
  z = qnorm(1-alpha/2)
  CI <- prop + c(-1,1)*z*sqrt(prop*(1-prop)/num)
  return(list(prop, CI[1], CI[2]))
}



incidence <- function(baseline, data.table = dt){
  n.pos <- sum(data.table[, ifelse(AKI.48 == 1, 1, ifelse( max.SCr.first.7.days / get(baseline) >= 1.5, 1, 0))])
  results <- prop.with.CI(num = n.pos, denom = nrow(data.table))
  #print(paste(n.pos, nrow(dt)))
  return(list(baseline, unlist(results[1]), unlist(results[2]), unlist(results[3])))
}

incidence.on.first <- function(baseline){
  n.pos <- sum(dt[, ifelse(AKI.48 == 1, 1, ifelse( first.IP.SCr / get(baseline) >= 1.5, 1, 0))])
  results <- prop.with.CI(num = n.pos, denom = nrow(dt))
  #print(paste(n.pos, nrow(dt)))
  return(list(baseline, unlist(results[1]), unlist(results[2]), unlist(results[3])))
}

# Equations for imputation
schwartz_bed_inv <- function(height, gfr = 107){
  return(.413*height/gfr)
}

schwartz_scr_inv <- function(height, gfr = 107){
  return ( ( (42.3*(height/100)^0.78) / gfr ) ^(1/0.78) )
}

fas_age <- function(age, is.female){
  return(
    ifelse(is.female == 1,
           0.23 + 0.034*age - 0.0018*age^2 + 0.00017*age^3 - 0.0000051*age^4,
           0.21 + 0.057*age - 0.0075*age^2 + 0.00064*age^3 -  0.000016*age^4)
         )
}

fas_height <- function(height){
  height = height/100
  return(
    3.94 - 13.4*height + 17.6*height^2 - 9.84*height^3 + 2.04*height^4
  )
}

predict.outcome <- function(baseline, outcome, threshold = 1.5, full = TRUE){
  #print(dt)
  if(outcome == "CKD3plus"){
    AKI.pos <- dt[CKD3plus != -1 , ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / get(baseline) >= threshold, 1, 0))]
    outcome.pos <- dt[CKD3plus != -1 , get(outcome)]

  } else {
    AKI.pos <- dt[, ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / get(baseline) >= threshold, 1, 0))]
    outcome.pos <- dt[ , get(outcome)]
  }
  
  TP <- sum( AKI.pos*outcome.pos )
  TN <- sum((1-AKI.pos)*(1-outcome.pos))
  FP <- sum(AKI.pos*(1-outcome.pos))
  FN <- sum((1-AKI.pos)*outcome.pos)
  Sn.with.CI <- prop.with.CI(TP, TP + FN)
  Sp.with.CI <- prop.with.CI(TN, TN + FP)
  
  if (full == TRUE){
    return(list(baseline, outcome, threshold, TP, TN, FP, FN, 
              unlist(Sn.with.CI[1]), unlist(Sn.with.CI[2]), unlist(Sn.with.CI[3]), 
              unlist(Sp.with.CI[1]), unlist(Sp.with.CI[2]), unlist(Sp.with.CI[3])))
  } else {
    return(c(sum(AKI.pos)/length(AKI.pos), unlist(Sn.with.CI[1]), unlist(Sp.with.CI[1])))
  }
}

predict.outcome.PreSCr <- function(baseline, outcome, threshold = 1.5, full = TRUE){

  if(outcome == "CKD3plus"){
    AKI.pos <- dt.with_preSCr[CKD3plus != -1 , ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / get(baseline) >= threshold, 1, 0))]
    outcome.pos <- dt.with_preSCr[CKD3plus != -1 , get(outcome)]

  } else {
    AKI.pos <- dt.with_preSCr[, ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / get(baseline) >= threshold, 1, 0))]
    outcome.pos <- dt.with_preSCr[ , get(outcome)]
  }
  
  TP <- sum( AKI.pos*outcome.pos )
  TN <- sum((1-AKI.pos)*(1-outcome.pos))
  FP <- sum(AKI.pos*(1-outcome.pos))
  FN <- sum((1-AKI.pos)*outcome.pos)
  Sn.with.CI <- prop.with.CI(TP, TP + FN)
  Sp.with.CI <- prop.with.CI(TN, TN + FP)
  
  if (full == TRUE){
    return(list(baseline, outcome, threshold, TP, TN, FP, FN, 
              unlist(Sn.with.CI[1]), unlist(Sn.with.CI[2]), unlist(Sn.with.CI[3]), 
              unlist(Sp.with.CI[1]), unlist(Sp.with.CI[2]), unlist(Sp.with.CI[3])))
  } else {
    return(c(sum(AKI.pos)/length(AKI.pos), unlist(Sn.with.CI[1]), unlist(Sp.with.CI[1])))
  }
}

predict.outcome.over12m <- function(baseline, outcome, threshold = 1.5, full = TRUE){

  if(outcome == "CKD3plus"){
    AKI.pos <- dt.with_preSCr_over12m[CKD3plus != -1 , ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / get(baseline) >= threshold, 1, 0))]
    outcome.pos <- dt.with_preSCr_over12m[CKD3plus != -1 , get(outcome)]

  } else {
    AKI.pos <- dt.with_preSCr_over12m[, ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / get(baseline) >= threshold, 1, 0))]
    outcome.pos <- dt.with_preSCr_over12m[ , get(outcome)]
  }
  
  TP <- sum( AKI.pos*outcome.pos )
  TN <- sum((1-AKI.pos)*(1-outcome.pos))
  FP <- sum(AKI.pos*(1-outcome.pos))
  FN <- sum((1-AKI.pos)*outcome.pos)
  Sn.with.CI <- prop.with.CI(TP, TP + FN)
  Sp.with.CI <- prop.with.CI(TN, TN + FP)
  
  if (full == TRUE){
    return(list(baseline, outcome, threshold, TP, TN, FP, FN, 
              unlist(Sn.with.CI[1]), unlist(Sn.with.CI[2]), unlist(Sn.with.CI[3]), 
              unlist(Sp.with.CI[1]), unlist(Sp.with.CI[2]), unlist(Sp.with.CI[3])))
  } else {
    return(c(sum(AKI.pos)/length(AKI.pos), unlist(Sn.with.CI[1]), unlist(Sp.with.CI[1])))
  }
}

#For correlation matrix plots
lowerfun <- function(data,mapping){
  ggplot(data = data, mapping = mapping)+
    geom_point(alpha = 0.1)+
    scale_x_continuous(limits = c(0, 1.2), breaks = c(0, 0.3, 0.6, 0.9))+
    scale_y_continuous(limits = c(0, 1.2), breaks = c(0, 0.3, 0.6, 0.9, 1.2))+
    geom_abline(slope = 1, intercept = 0, color = "darksalmon", linetype = "dashed", size = 1)+
    geom_smooth()
}

#For correlation matrix plots
diagfun <- function(data, mapping){
  ggplot(data = data, mapping = mapping) +
    geom_density(fill = "turquoise") +
    scale_x_continuous(limits = c(0, 1.2), breaks = c(0, 0.3, 0.6, 0.9))
    #geom_label(label = "test", aes(y = 0.6, x = 0.6))
}
dt.bldisp <- data.table(
  AKI_def = factor(c("AKI.Schwartz.bedside.120", "AKI.Schwartz.bedside.107", "AKI.FAS.height", "AKI.FAS.age", "AKI.FAS.age.tbl", "AKI.bl.min", "AKI.lowest.FAS", "AKI.bl.second.min", "AKI.second.lowest.FAS", "AKI.bl.prox", "AKI.prox.FAS", "AKI.bl.mean", "AKI.mean.FAS")),
  baseline = c("Schwartz.bedside.120", "Schwartz.bedside.107", "FAS.height", "FAS.age", "FAS.age.tbl", "bl.min", "lowest.FAS", "bl.second.min", "second.lowest.FAS", "bl.prox", "prox.FAS", "bl.mean", "mean.FAS"),
  bl.disp = factor(c("Schwartz 120", "Schwartz 107", 
                     "FAS height", "FAS", "FAS table", 
                     "Lowest", "Lowest-FAS", "Second Lowest", "Second-Lowest-FAS", "Proximal", 
                     "Proximal-FAS", "Mean", "Mean-FAS"), 
                   levels = c("Schwartz 120", "Schwartz 107", 
                     "FAS height", "FAS", "FAS table", 
                     "Lowest", "Lowest-FAS", "Second Lowest", "Second-Lowest-FAS", "Proximal", 
                     "Proximal-FAS", "Mean", "Mean-FAS")),
  bl.disp.first = factor(c("Schwartz 120", "Schwartz 107", 
                           "FAS height", "FAS", "FAS table", 
                           "Lowest-Schwartz", "Lowest-FAS", 
                           "Second Lowest-Schwartz", "Second-Lowest-FAS", 
                           "Proximal-Schwartz", "Proximal-FAS", 
                           "Mean-Schwartz", "Mean-FAS"), 
                   levels = c("Schwartz 120", "Schwartz 107", 
                           "FAS height", "FAS", "FAS table", 
                           "Lowest-Schwartz", "Lowest-FAS", 
                           "Second Lowest-Schwartz", "Second-Lowest-FAS", 
                           "Proximal-Schwartz", "Proximal-FAS", 
                           "Mean-Schwartz", "Mean-FAS")),
  impute.type = factor(c("Schwartz imputation", "Schwartz imputation", 
                         "FAS imputation", "FAS imputation", 
                         "FAS imputation", 
                         "Schwartz imputation", "FAS imputation", 
                         "Schwartz imputation", "FAS imputation", 
                         "Schwartz imputation", "FAS imputation", 
                         "Schwartz imputation", "FAS imputation"),
            levels = c("Schwartz imputation", "FAS imputation")),
  empiric.type = factor(c("Imputed", "Imputed", "Imputed", "Imputed", "Imputed", 
                          "Lowest-1st", "Lowest-1st", 
                          "Second Lowest-1st", "Second Lowest-1st", 
                          "Proximal-1st", "Proximal-1st", 
                          "Mean-1st", "Mean-1st"
),
                   levels = c("Imputed", "Lowest-1st",
                              "Second Lowest-1st", 
                              "Proximal-1st", "Mean-1st"))
)

#Generates diagonal lines for ROC plots
diag_fun <- function(x) -1*x+1

#Read data
dt <- fread("AKI_processed_dataframe.csv")
#dt[ , bl.imp := NULL]
n_pt <- nrow(dt)
n_pt_CKD <- nrow(dt[CKD3plus != -1])

percent.with.prelabs <- nrow(dt[n.prelabs > 0])/n_pt
percent.mult.prelabs <- nrow(dt[n.prelabs > 1])/n_pt

dt[ , age.yr := adm.age/(60*60*24*365)]
dt[ , is.female := ifelse(Sex.Title == "FEMALE", 1, ifelse(Sex.Title == "MALE", 0, -1))]

# Add imputed baselines to the data.table
dt[ , Schwartz.bedside.107 := schwartz_bed_inv(Huse, gfr = 107.3)]
dt[ , Schwartz.SCr.107 := schwartz_scr_inv(Huse, gfr = 107.3)]
dt[ , Schwartz.bedside.120 := schwartz_bed_inv(Huse, gfr = 120)]
dt[ , Schwartz.SCr.120 := schwartz_scr_inv(Huse, gfr = 120)]
dt[ , FAS.age := fas_age(age.yr, is.female)]
dt[ , FAS.height := fas_height(Huse)]
dt[ , lowest.FAS := ifelse(n.prelabs == 0, FAS.age, bl.min)]
dt[ , second.lowest.FAS := ifelse(n.prelabs == 0, FAS.age, bl.second.min)]
dt[ , mean.FAS := ifelse(n.prelabs == 0, FAS.age, bl.mean)]
dt[ , prox.FAS := ifelse(n.prelabs == 0, FAS.age, bl.prox)]

dt[ , GFR.min := 0.413*(Huse/max.SCr)]
dt[ , AKI.bl.min := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / bl.min >= 1.5, 1, 0))]
dt[ , AKI.bl.second.min := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / bl.second.min >= 1.5, 1, 0))]
dt[ , AKI.bl.mean := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / bl.mean >= 1.5, 1, 0))]
dt[ , AKI.bl.prox := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / bl.prox >= 1.5, 1, 0))]

dt[ , AKI.Schwartz.bedside.107 := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / Schwartz.bedside.107 >= 1.5, 1, 0))]
dt[ , AKI.Schwartz.SCr.107 := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / Schwartz.SCr.107 >= 1.5, 1, 0))]
dt[ , AKI.Schwartz.bedside.120 := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / Schwartz.bedside.120 >= 1.5, 1, 0))]
dt[ , AKI.Schwartz.SCr.120 := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / Schwartz.SCr.120 >= 1.5, 1, 0))]
dt[ , AKI.FAS.age := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / FAS.age >= 1.5, 1, 0))]
dt[ , AKI.FAS.height := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / FAS.height >= 1.5, 1, 0))]
dt[ , AKI.lowest.FAS := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / lowest.FAS >= 1.5, 1, 0))]
dt[ , AKI.second.lowest.FAS := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / second.lowest.FAS >= 1.5, 1, 0))]
dt[ , AKI.mean.FAS := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / mean.FAS >= 1.5, 1, 0))]
dt[ , AKI.prox.FAS := ifelse(AKI.48 == 1, 1, ifelse(max.SCr.first.7.days / prox.FAS >= 1.5, 1, 0))]

# Add "raw AKI" score for ROC analysis
dt[ , raw.AKI.Schwartz.bedside.120 := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/Schwartz.bedside.120)]
dt[ , raw.AKI.Schwartz.bedside.107 := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/Schwartz.bedside.107)]
dt[ , raw.AKI.FAS.age := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/FAS.age)]
dt[ , raw.AKI.FAS.height := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/FAS.height)]
dt[ , raw.AKI.bl.min := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/bl.min)]
dt[ , raw.AKI.bl.second.min := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/bl.second.min)]
dt[ , raw.AKI.bl.mean := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/bl.mean)]
dt[ , raw.AKI.bl.prox := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/bl.prox)]
dt[ , raw.AKI.lowest.FAS := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/lowest.FAS)]
dt[ , raw.AKI.second.lowest.FAS := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/second.lowest.FAS)]
dt[ , raw.AKI.mean.FAS := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/mean.FAS)]
dt[ , raw.AKI.prox.FAS := ifelse(AKI.48 == 1, 100, max.SCr.first.7.days/prox.FAS)]


#Exclusion counts
n_tot_stay <- nrow(dt)

bl.imp <- c("Schwartz.bedside.120", "Schwartz.bedside.107", "FAS.height", "FAS.age")
#bl.imp.107 <- c("Schwartz.SCr.107", "Schwartz.bedside.107", "FAS.height", "FAS.age")


pct.diagnosed <- function(dt, stage, bl){
  raw.AKI <- paste0("raw.AKI.", bl)
  if(stage == 1){
    AKI.vect <- dt[ , (max.SCr.first.7.days/get(bl) >= 1.5) | (AKI.48 == TRUE) ]
  } else if (stage == 2){
    AKI.vect <- dt[ , (max.SCr.first.7.days/get(bl) >= 2.0) ]
  } else if (stage == 3){
    AKI.vect <- dt[ , (max.SCr.first.7.days/get(bl) >= 3.0) | 
                      (GFR.min <= 35) |
                      (max.SCr >= 4.0)]
  }
  
  
  return(
    paste0(
      sum(AKI.vect & dt[, AKI.Diagnosis == 1]), 
      "/", 
      sum(AKI.vect), 
      "\n(", 
      sprintf("%.1f%%", 
              100*sum(AKI.vect & dt[, AKI.Diagnosis == 1])/sum(AKI.vect)
              ), 
      ")"
      )
  )
}


```

```{r preprocess_data, include=FALSE, eval=FALSE}

library(data.table)
library(jsonlite)


####
# READ IN DATA
####

# h = heights (for calculating post-stay EGFR)
# labs = SCr (for calculating post-stay EGFR )
# inputs contains all data for the stay

h <- readRDS("G:/My Drive/AKI and CKD/July 2019 data/heights.rds")
h <- data.table(h)
labs <- readRDS("G:/My Drive/AKI and CKD/July 2019 data/labs.rds")
labs <- data.table(labs)
input <- fromJSON("G:/My Drive/AKI and CKD/July 2019 data/input_hashed.json")
input <- data.table(input)
print("Initial number of labs")

colnames(labs) <- c("lab.date", "lab.time", "lab.name", "lab.value", "value", "DI.MRN")

time.to.seconds <- function(y){
  
  sapply(strsplit(y,":"),
  function(x) {
    x <- as.numeric(x)
    (x[1]*60+x[2])*60
    }
)
}

####
# Preprocess labs
####

#Create date/time field. Adding 4070880422 to convert to EPIC format, which is seconds since 1840-12-31
labs[ , dt := as.numeric(as.POSIXct(lab.date)) + time.to.seconds(lab.time) - as.numeric(as.POSIXct("1840-12-31"))]

labs[,c("lab.time", "lab.value", "lab.date") := NULL]

#Delete lab values that are too large or equal to 0
labs <- labs[value<=20]
labs <- labs[value>0]
print("labs after eliminating those with unreasonable value")
print(nrow(labs))

# # ###check for close together lab results
# labsClose <- rbind(labs)
# setkey(labsClose, DI.MRN, dt)
# labsClose[order(dt) , lastval := shift(value), by=DI.MRN]
# labsClose[order(dt) , lastval.dt := shift(dt), by=DI.MRN]
# labsClose[ , dtDif := dt - lastval.dt]
# labsClose[ , hrs.Dif := dtDif/(60*60)]
# labsSame <- labsClose[ hrs.Dif == 0 ]

# Eliminate any labs that were entered as at the exact same second as others
labs[ , dups := .N, by=c("DI.MRN", "dt")]
labs <- labs[dups == 1]
print("labs after eliminating those with duplicate timestamp")

####
# Preprocess heights
####

colnames(h) <- c("encounter", "enc.date", "enc.date2", "h", "DI.MRN")


h[ , dt := as.numeric(as.POSIXct(enc.date)) - as.numeric(as.POSIXct("1840-12-31"))]



###convert height to h.inch and h.cm
h$h <- as.character(h$h)
h$ft <- sapply(strsplit(h$h,"'"), `[`, 1)
h$ft <- as.numeric(h$ft)
h$inch <- sapply(strsplit(h$h,"'"), `[`, 2)
h$inch <- gsub("\"", "",h$inch)
h$inch <- as.numeric(h$inch)
h$h.inch <- h$ft*12+h$inch
h$h.cm <- h$h.inch*2.54
h[, c("ft","inch", "h","h.inch", "encounter", "enc.date", "enc.date2"):=NULL]

# h[ , dups := .N, by=c("DI.MRN", "dt")]
# table(h$dups)
# h_dups <- h[dups > 1]
# h <- h[dups == 1]
# h2 <- cbind(h)
#If multiple heights recorded at same second, take the average of them (usually they're identical anyway)
h[ , h.cm := mean(h.cm), by=c("DI.MRN", "dt")]
h <- unique(h, by = c("DI.MRN", "dt"))

print("labs after eliminating those with duplicate timestamp")
print(nrow(labs))

####
# Pull sex, age from input to labs and h
####

## Format input to match dt, MRN
setnames(input, "SURF_MRN", "DI.MRN")


### Add age in days and sex to h (height table) and labs

#First read in birthdays
bdays <- unique(input[, c("DI.MRN", "Birth.Date")])
setkey(input, DI.MRN)
setkey(h, DI.MRN)
setkey(labs, DI.MRN)

# #heights
# h <- h[DI.MRN %in% input$DI.MRN]
# h <- h[bdays][is.na(dt)==FALSE]
#labs
labs <- labs[DI.MRN %in% input$DI.MRN]
#After eliminating labs with MRN that doesn't match our input table
print(nrow(labs))
labs <- labs[bdays][is.na(dt)==FALSE]

h <- h[DI.MRN %in% input$DI.MRN]
h <- h[bdays][is.na(dt)==FALSE]

#Calculate age
labs[ , age := dt - Birth.Date]
labs[ , Birth.Date := NULL]
h[ , age := dt - Birth.Date]
h[ , Birth.Date := NULL]
#Pull in sex
sex <- unique(input[, c("DI.MRN", "Sex.Title")])
setkey(labs, DI.MRN)
labs <- labs[sex][is.na(dt)==FALSE]

setkey(h, DI.MRN)
h <- h[sex][is.na(dt)==FALSE]

rm(bdays); rm(sex)


# Rm height and lab values for patients too young or too old
h <- h[age >= 90*24*60*60]
h <- h[age <= 7000*24*60*60]
labs <- labs[age >= 90*24*60*60]
labs <- labs[age <= 7000*24*60*60]
#After applying age filter
print(nrow(labs))

####
# Create regression model for imputing missing heights
####

#get height model
h.male = h[Sex.Title == "MALE"]
Hmod.male <- lm(h.cm~poly(age,3), data=h.male)
summary(Hmod.male)

h.female = h[Sex.Title == "FEMALE"]
Hmod.female <- lm(h.cm~poly(age,3), data=h.female)
summary(Hmod.female)
rm(h.female); rm(h.male)

# # plot
# library(ggplot2)
# ggplot(h, aes(x=age, y=h.cm, color=Sex.Title)) + geom_point(size=.05)+ ylim(0, 200) +
#  geom_line(aes(y = predict(Hmod.male, data.frame(age = h$age))), size = 1.5, linetype = 'solid', color = "royalblue4")+
#  geom_line(aes(y = predict(Hmod.female, data.frame(age = h$age))), size = 1.5, linetype = 'solid', color = "red4")

###add predicted height to labs (does not consider any actual height reading)
labs[, hPred.male := predict(Hmod.male, data.frame(age = labs$age)) ]
labs[ , hPred.female := predict(Hmod.female, data.frame(age = labs$age)) ]
labs[ , hPred := ifelse(Sex.Title == "MALE", hPred.male, hPred.female) ]
labs[, c("hPred.male","hPred.female") := NULL]


####
# Pull correct height into labs
####

# If no heigh in data, impute. If heigh available but more than 90 days before/after lab, adjust with regression. If nearest height is within 90 days, use.

### Add adjusted height reading to labs
# take nearest height to labs append in labsWH table
#labs <- labs[DI.MRN %in% input$DI.MRN]
labs[,join_dt := dt]
h[,join_dt := dt]

labs[,lab_dt := dt]
h[,h_dt := dt]

labs[,lab_age := age]
h[,h_age := age]

setkey(labs, DI.MRN, join_dt)
setkey(h, DI.MRN, join_dt)
labsWH <- h[labs, roll = "nearest"]
labsWH <- labsWH[ is.na(h.cm) == FALSE]
labsWH <- labsWH[ , c("DI.MRN", "h.cm", "Sex.Title", "h_dt", "h_age", "lab_dt", "lab_age")]
# hWlabs <- h[labs, roll = "nearest"]
# hWlabs <- hWlabs[ DI.MRN %in% h$DI.MRN]
#In labsWH calculate height adjusted: height value + [Hmod(labAge) - Hmod(heightAge)]
labsWH[ , sec_fr_h := lab_dt - h_dt]
labsWH[ , Hadjust.male := h.cm + predict(Hmod.male, data.frame(age = lab_age)) - predict(Hmod.male, data.frame(age = h_age)) ]
labsWH[ , Hadjust.female := h.cm + predict(Hmod.female, data.frame(age = lab_age)) - predict(Hmod.female, data.frame(age = h_age)) ]
labsWH[ , Hadjusted := ifelse(Sex.Title == "MALE", Hadjust.male, Hadjust.female) ]
labsWH[ , Huse := ifelse(abs(sec_fr_h) < 90*24*60*60, h.cm, Hadjusted) ]


#pull Huse back into labs
Huse_tbl <- labsWH[ , c("DI.MRN", "lab_dt", "Huse")]

Huse_tbl[ , join_dt := lab_dt]
setkey(Huse_tbl, DI.MRN, join_dt)
setkey(labs, DI.MRN, join_dt)


labs <- merge(labs, Huse_tbl, by=c("DI.MRN", "join_dt"), all.x = TRUE)

labs[ , H := ifelse(!is.na(Huse), Huse, hPred)]
labs[,c("join_dt","lab.name","hPred","lab_dt.x","lab_dt.y", "lab_age") := NULL]

rm(labsWH); rm(Huse_tbl)

####
# Pull correct height into input
####

#Get height at admission
input[ , c("h.inch", "h.dt") := tstrsplit(Adm.Height, '\\|')]
input[ , h.dt := as.numeric(h.dt)]
input[ , h.inch := as.numeric(h.inch)]
input[ , h.cm := h.inch*2.54]
input[ , adm.age := Admit.Date.Time - Birth.Date]
input[ , h.age := h.dt - Birth.Date]

#MOVING LATER TO GET EXCLUSION CRITERIA COUNTS
#input <- input[adm.age >= 90*24*60*60]
#input <- input[adm.age <= 7000*24*60*60]

input[ , sec_fr_h := h.dt - Admit.Date.Time]

input[ , Hadjust.male := h.cm + predict(Hmod.male, data.frame(age = adm.age)) - predict(Hmod.male, data.frame(age = h.age)) ]
input[ , Hadjust.female := h.cm + predict(Hmod.female, data.frame(age = adm.age)) - predict(Hmod.female, data.frame(age = h.age)) ]

input[ , Hadjusted := ifelse(Sex.Title == "MALE", Hadjust.male, Hadjust.female) ]

input[, H.regress := ifelse(Sex.Title == "MALE", predict(Hmod.male, data.frame(age = adm.age)),
                            predict(Hmod.female, data.frame(age = adm.age)))]


input[ , Huse := ifelse(abs(sec_fr_h) < 90*24*60*60, h.cm, Hadjusted) ]
input[ , Huse := ifelse(is.na(Huse), H.regress, Huse)]

input[ , Hsource := ifelse(is.na(h.cm), "Imputed", 
                           ifelse(abs(sec_fr_h) < 90*24*60*60, "Measured",
                                  "Adjusted"))]

input[ , c("h.dt", "sec_fr_h", "h.inch", "h.age", "Hadjust.male", "Hadjust.female", "Hadjusted", "H.regress") := NULL]



####
# Assign CKD3+ labels
####

# Calculate eGFR in labs
# Define post-stay window for each stay
# Assign post-stay CKD labels
# If 2+ eGFR > 60, CKD3 == 1, 
# CKD3plus = 1 (positive), 0 (negative), or -1 (unlabeled)

labs[ , eGFR := 0.413*H/value ]
labs[ , CKD3plus := ifelse(eGFR < 60, 1, 0)]

input[ , post.start := Discharge.Date.Time + 90*24*60*60]
input[ , post.end := Discharge.Date.Time + 450*24*60*60]


#count labs in post-window for given CSN matching the MRN
setkey(input, DI.MRN, post.start, post.end)
labs[, dt2:=dt]
postlabs <- foverlaps(labs, input, 
                      by.x = c("DI.MRN","dt", "dt2"), 
                      by.y = c("DI.MRN","post.start","post.end"))

postlabs <- postlabs[!is.na(SURF_CSN)]

postlabs <- postlabs[, c("SURF_CSN", "dt", "CKD3plus")]


###add countPostLabs to stays ###
postlabs[, count := .N , by=SURF_CSN]
input[ , countPostLabs := postlabs$count[match(input$SURF_CSN, postlabs$SURF_CSN)]]
input[, countPostLabs := ifelse(is.na(countPostLabs), 0, countPostLabs)]

postlabs$count <- NULL


### Calculate CKD presence and add to stays. 
postlabs[ , countPositives := sum(CKD3plus),by = SURF_CSN]
input$CKD3plus <-0
input[ ,CKD3plus := ifelse(postlabs$countPositives[match(input$SURF_CSN,postlabs$SURF_CSN)] > 1,
                           1,
                           0)]
input[,CKD3plus := ifelse(countPostLabs < 2, -1, CKD3plus)]
input[ , c("post.start","post.end","countPostLabs") := NULL]

table(input$CKD3plus)

######
#
# Analysis of how close together post-stay labs were for CKD+
#
######
input.CKDpos <- input[CKD3plus == 1]
postlabs <- postlabs[CKD3plus > 0]

first_pos <- postlabs[ , min(dt), by = SURF_CSN]
last_pos <- postlabs[ , max(dt), by = SURF_CSN]
colnames(first_pos) <- c("SURF_CSN", "dt_first_CKD_lab")
colnames(last_pos) <- c("SURF_CSN", "dt_last_CKD_lab")

input.CKDpos <- input.CKDpos[first_pos, on = "SURF_CSN", nomatch = 0]

input.CKDpos <- input.CKDpos[last_pos, on = "SURF_CSN", nomatch = 0]
input.CKDpos[ , days_between_CKDlabs := (dt_last_CKD_lab - dt_first_CKD_lab)/(60*60*24)]

ggsave("cdf_time_between_poststay_GFR.png",
  ggplot(input.CKDpos)+
    stat_ecdf(aes(x = days_between_CKDlabs))+
    scale_x_continuous(limits = c(0, 360))+
    scale_y_continuous(labels = scales::percent)+
    xlab("Interval between qualifying post-hospitalization SCr (days)")+
    ylab("Cumulative % of CKD-positive patients"),
  width = 5, height = 5, units = "in")

1 - ecdf(input.CKDpos$days_between_CKDlabs)(c(1, 7, 30, 60, 90))



fwrite(input, "input_labled_Aug2019v3.csv")



# rm(list = ls())


```

```{r generate_dataframe, include=FALSE, eval=FALSE}

library(data.table)



#Read data
dt <- fread("input_labled_Aug2019v3.csv")

#Exclusion counts
n_tot_stay <- nrow(dt)
n_old <- dt[adm.age >= 19*365*24*60*60, .N]
n_young <- dt[adm.age< 90*24*60*60, .N]
n_missing_SCr <- dt[ nchar(Creatinine) == 0, .N]
n_renal <- dt[ Renal == 1, .N]


#Apply exclusion criteria
dt <- dt[adm.age <= 19*365*24*60*60]
dt <- dt[adm.age >= 90*24*60*60]

dt <- dt[ nchar(Creatinine) > 0]

dt <- dt[ Renal == 0]

# 26 728 stays

#Remove superfluous columns

dt <- dt[ , IP.Mortality := ifelse(Disch.Disp == "EXPIRED", 1, 0)]
dt <- dt[ , c("Adm.Height", "Admit.Date.Time", "Discharge.Date.Time", "Creatinine", "AKI.Diagnosis", "Deceased", "CKD3plus", "IP.Mortality", "adm.age", "Huse", "Sex.Title", "Race.Title", "Hsource")]

# Add imputed baseline
dt[ , bl.imp := (.413/120)*Huse]

dt[ , n.iplabs := 0]
dt[ , bl.min := -5.0]
dt[ , bl.mean := -5.0]
dt[ , bl.prox := -5.0]
dt[ , bl.second.min := -5.0]
# Loop through all stays and 
#  - generate table of labs with timestamps
#  - use those in 92-day prestay window to generate min/second min/mean/prox baselines
#  - use those from the stay to diagnose AKI under 48-hr absolute increase criteria
#  - use each baseline to apply relative increase criteria to SCr within 7 days of admit
#  - use each bl to diagnose AKI on first
#  - identify max inpatient SCr
for(row in 1: nrow(dt)){
  
  #Get dates for stay
  adm.temp = dt[row, Admit.Date.Time]
  dsc.temp = dt[row, Discharge.Date.Time]
  
  #Create temporary labs table for this stay 
  temp = tstrsplit(tstrsplit(dt[row, Creatinine], split = '\\;'), split = '\\|')
  
  labs = data.table(
    lab = unlist(temp[1]),
    date = unlist(temp[2])
  )
  labs[, lab := as.numeric(lab)]
  labs[, date := as.numeric(date)]
  labs <- labs[ !is.na(lab)]
  labs <- unique(labs)
  labs <- labs[date >= adm.temp - 92*24*60*60]
  # Add period (pre/ip/post)
  labs[ , period := ifelse(date < adm.temp, 'pre', ifelse(date < dsc.temp, 'ip', 'post'))]
  
  
  ip.labs <- labs[period == 'ip']
  dt[row , n.prelabs := nrow(labs[period == 'pre']) ]
  dt[row , n.iplabs := nrow(labs[period == 'ip']) ]

  if (dt[row , n.iplabs ] > 0){
    #add baselines
    if(dt[row , n.prelabs] > 0 ){
      dt[row , bl.min := labs[ period == 'pre', min(lab)]]
      dt[row , bl.mean := labs[ period == 'pre', mean(lab)]]
      dt[row , bl.prox := labs[ period == 'pre'][which.max(date), lab]]
      dt[row , bl.second.min := sort(labs[period == 'pre', lab])[2]]
    }
    
    dt[row, bl.min := ifelse(is.na(bl.min)| bl.min == -5.0, bl.imp, bl.min)]
    dt[row, bl.mean := ifelse(is.na(bl.mean) | bl.mean == -5.0, bl.imp, bl.mean)]
    dt[row, bl.prox := ifelse(is.na(bl.prox) | bl.prox == -5.0, bl.imp, bl.prox)]
    dt[row, bl.second.min := ifelse(is.na(bl.second.min) | bl.second.min == -5.0, bl.min, bl.second.min)] 
    
    #Diagnose AKI based on 48hr absolute increase
    j = 1
    has.AKI.48 = 0
    while (has.AKI.48 == 0 & j <= nrow(ip.labs)){
      this.date = ip.labs[j, date]
      this.val = ip.labs[j, lab]
      inc.48 = 0
      for (k in nrow(ip.labs)) {
        if(this.date - ip.labs[k, date] < 48*60*60 ){
          inc.48 = max(inc.48, this.val - ip.labs[k, lab])
          #print(inc.48)
        }
      }
      if(inc.48 >= 0.3) { has.AKI.48 <- 1 }
      j = j + 1
    }
    dt[row, AKI.48 := has.AKI.48]
    
    #USe baselines to diagnose AKI based on relative increase
    dt[row, max.SCr.first.7.days := max(ip.labs[date <= adm.temp + 7*24*60*60]$lab)]
    dt[row, max.SCr := max(ip.labs$lab)]
    dt[row, first.IP.SCr := ip.labs[which.min(date), lab]]
    
    # dt[row, AKI.rel.imp := has.AKI.rel(baseline = bl.imp, labs = ip.labs[date <= adm.temp + 7*24*60*60]$lab, threshold = 0.5) ]
    # dt[row, AKI.rel.min := has.AKI.rel(baseline = bl.min, labs = ip.labs[date <= adm.temp + 7*24*60*60]$lab, threshold = 0.5) ]
    # dt[row, AKI.rel.second.min := has.AKI.rel(baseline = bl.second.min, labs = ip.labs[date <= adm.temp + 7*24*60*60]$lab, threshold = 0.5) ]
    # dt[row, AKI.rel.mean := has.AKI.rel(baseline = bl.mean, labs = ip.labs[date <= adm.temp + 7*24*60*60]$lab, threshold = 0.5) ]
    # dt[row, AKI.rel.prox := has.AKI.rel(baseline = bl.prox, labs = ip.labs[date <= adm.temp + 7*24*60*60]$lab, threshold = 0.5) ]
    # 
    # #Combine AKI.48 and rel increase
    # dt[row, AKI.imp := max(AKI.48, AKI.rel.imp)]
    # dt[row, AKI.min := max(AKI.48, AKI.rel.min)]
    # dt[row, AKI.second.min := max(AKI.48, AKI.rel.second.min)]
    # dt[row, AKI.mean := max(AKI.48, AKI.rel.mean)]
    # dt[row, AKI.prox := max(AKI.48, AKI.rel.prox)]
  }
}

# total 26 728

#Remove patients with no IP labs
dt <- dt[n.iplabs > 0]
# 20 450

#remive patients with inpatient SCr after the first 7 days
dt <- dt[max.SCr.first.7.days != -Inf]
# 20 373

dt[ , c("Adm.Height", "Creatinine", "Admit.Date.Time", "Discharge.Date.Time", "Deceased") := NULL]

fwrite(dt, "AKI_processed_dataframe.csv")

```

```{r fig_inclusion_diag, include = FALSE, eval=FALSE}
#Read data
dt_all <- fread("input_labled_Aug2019v3.csv")
n_old <- dt_all[adm.age >= 19*365*24*60*60, .N]
n_young <- dt_all[adm.age< 90*24*60*60, .N]
n_missing_SCr <- dt_all[ nchar(Creatinine) == 0, .N]
n_renal <- dt_all[ Renal == 1, .N]


library(DiagrammeR)
library(rsvg)
library(DiagrammeRsvg)

inclusion_flow <- grViz("digraph flowchart {
      # node definitions with substituted label text
      node [fontname = Helvetica, shape = rectangle]
      tab1 [label = '@@1']
      tab2 [label = '@@2']
      tab3 [label = '@@3']
      tab4 [label = '@@4']
      tab5 [label = '@@5']

      # edge definitions with the node IDs
      tab1 -> tab2 -> tab3; tab1 -> tab4; tab2 -> tab5;
      }

      [1]: paste0(format(n_tot_stay, big.mark=','),  ' patient stays')
      [2]: paste0(format(n_pt, big.mark=','), ' stays included\\n(incidence, inpatient mortality,\\n and diagnosis code analyses)')
      [3]: paste0(format(n_pt_CKD, big.mark=',')  , ' stays included\\n(CKD analysis)')
      [4]: paste0(format(n_tot_stay - n_pt, big.mark=','), ' stays excluded', '\\n', format(n_young, big.mark=','), ' for age < 3 months', '\\n', format(n_old, big.mark=','), ' for age > 18 years', '\\n', format(n_missing_SCr, big.mark=','), ' for < 1 SCr in first 7 days', '\\n', format(n_renal, big.mark=','), ' for renal diagnosis')
      [5]: paste0(format(n_pt - n_pt_CKD, big.mark=',')  , ' stays excluded \\n from CKD analysis \\n for < 2 post-stay SCr')
      ")



png::writePNG(rsvg(charToRaw(export_svg(inclusion_flow))), "inclusion_flow.png")


```

```{r fas_age_table, include = FALSE, eval = FALSE}
FAS_age_table <- fread("FAS_age_table.csv")

dt[ , age.yr.int := floor(age.yr)]
dt[ , age.yr.int := ifelse(age.yr.int == 0, 1, age.yr.int)]

dt <- dt[FAS_age_table, on = c("Sex.Title", "age.yr.int"), nomatch = 0]

ggplot(melt(dt[, c("age.yr", "FAS.age", "FAS.age.tbl", "Sex.Title")], id.vars = c("age.yr", "Sex.Title")))+
  geom_line(aes(x = age.yr, y = value, color = variable, linetype = Sex.Title))



```

```{r impute_analysis, include = FALSE}

# dt.small <- dt[sample(.N, 100)]




# PAIRWISE CORRELATION MATRIX

dt.imp <- dt[ , ..bl.imp]


setnames(dt.imp, old = dt.bldisp$baseline, new = as.character(dt.bldisp$bl.disp), skip_absent = TRUE)



# plt.pairs.imp <- ggpairs(data = dt.imp,
#                     diag = list(continuous = diagfun),
#                     lower = list(continuous = lowerfun)
#                     ) + theme_bw()
# 
# ggsave("imputed_bl_pairwise_fig.png", plt.pairs.imp, width = 6, height = 5, units = "in")

# MEAN AND STANDARD DEVIATION

dt.imputed.bls <- melt(dt[ , ..bl.imp], measure.vars = bl.imp, variable.name = "baseline")[dt.bldisp, on = "baseline", nomatch = 0]

plt.mean.sd.imputed <- 
  ggplot(dt.imputed.bls, aes(col = bl.disp))+
         stat_summary(aes(x = bl.disp, y = value), fun.data = function(y){
           data.frame(y = mean(y),
               ymin = mean(y) - sd(y), 
               ymax = mean(y) + sd(y))})+
         ylab("Mean and standard deviation baseline SCr (mg/dL)")+
         xlab("Baseline imputation method")+
  coord_flip()+
  theme(legend.position = "none")
  
ggsave("baseline_mean_sd_imp_fig.png", plt.mean.sd.imputed, width = 6, height = 3, units = "in")


# AGREEMENT BY HEIGHT AND AGE

dt.plt <- rbind(
  cbind(facet = "By height (cm)", melt(dt[, c("Huse", ..bl.imp)], id.vars = "Huse")),
  cbind(facet = "By age (years)", melt(dt[, c("age.yr", ..bl.imp)], id.vars = "age.yr")),
  use.names = FALSE
)
setnames(dt.plt, "variable", "baseline")
dt.plt <- dt.plt[dt.bldisp, on = "baseline", nomatch = 0]

plt.imputed.by.age.height <- ggplot(dt.plt, aes(x = Huse, y = value, color = bl.disp)) +
  facet_grid( cols = vars(facet), scales = "free_x") +
  geom_point(alpha = 0.05, size = 0.15)+
  geom_smooth(se = FALSE)+
  theme(legend.position = "bottom")+
  ylab("Baseline SCr (mg/dL)")+
  xlab("")+scale_color_discrete(name = "Imputation method") +
  coord_cartesian(ylim = c(0, 1.6))

ggsave("imputed_by_age_height_fig.png", height = 5, width = 6, units = "in")


#
# Incidence
#

dt.inc.imputed <- data.table(
  baseline <- character(),
  inc <- numeric(),
  inc.LB <- numeric(),
  inc.UB <- numeric()
)

for(bl in bl.imp){
  dt.inc.imputed <- 
    rbind(dt.inc.imputed,
        incidence(bl), use.names = FALSE)
}

colnames(dt.inc.imputed) <- c("baseline", "inc", "inc.LB", "inc.UB")

dt.inc.imputed <- dt.inc.imputed[dt.bldisp, on = "baseline", nomatch = 0]


p.inc.imputed <- 
  ggplot(dt.inc.imputed, aes(x=bl.disp, y=inc, ymin = inc.LB, ymax = inc.UB, col = bl.disp))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  theme(legend.position="none")+
  scale_y_continuous(labels = scales::percent)+
  ylab("Incidence")+
  xlab("Baseline Determination Method")+
  coord_flip()

ggsave("p_inc_imputed_fig.png", plot = p.inc.imputed, device = "png", width = 6, height = 3, units = "in")


#dt.inc.imputed[ , baseline := factor(baseline, levels = c("bl.imp","bl.min","bl.second.min", "bl.prox", "bl.mean")) ]

#
# Clinical outcomes
#

dt.imp.outcomes <- data.table(
  baseline <- character(),
  outcome <- character(),
  threshold <- numeric(),
  TP <- numeric(),
  TN <- numeric(),
  FP <- numeric(),
  FN <- numeric(),
  Sn <- numeric(),
  Sn.LB <- numeric(),
  Sn.UB <- numeric(),
  Sp <- numeric(),
  Sp.LB <- numeric(),
  Sp.UB <- numeric()
)

colnames(dt.imp.outcomes) <- c("baseline", "Outcome", "Threshold", "TP", "TN", "FP", "FN", "Sn", "Sn.LB", "Sn.UB", "Sp", "Sp.LB", "Sp.UB")
  
outcomes <- c("AKI.Diagnosis", "CKD3plus", "IP.Mortality")
for(bl in bl.imp){
  for(outcome in outcomes){
     dt.imp.outcomes <- 
       rbind(dt.imp.outcomes,
             predict.outcome(bl, outcome), use.names = FALSE)
  }
}


Outcome.labs <- c("AKI ICD Code", "Post-hospitalization CKD", "Inpatient Mortality")
names(Outcome.labs) <- c("AKI.Diagnosis", "CKD3plus", "IP.Mortality")
point.type.labs <- c("Sensitivity", "Specificity")
names(point.type.labs) <- c("Sn", "Sp")

dt.outcomes.melted <- melt(dt.imp.outcomes, id.vars = c("baseline", "Outcome"),
                           measure.vars = c("Sn", "Sp"),
                           variable.name = "point.type",
                           value.name = "point.est")


dt.outcomes.melted.LB <- melt(dt.imp.outcomes, id.vars = c("baseline", "Outcome"),
                           measure.vars = c("Sn.LB", "Sp.LB"),
                           variable.name = "point.type",
                           value.name = "LB")

dt.outcomes.melted.UB <- melt(dt.imp.outcomes, id.vars = c("baseline", "Outcome"),
                           measure.vars = c("Sn.UB", "Sp.UB"),
                           variable.name = "point.type",
                           value.name = "UB")

dt.outcomes.melted <- cbind(dt.outcomes.melted, dt.outcomes.melted.LB$LB, dt.outcomes.melted.UB$UB)
colnames(dt.outcomes.melted) <- c("baseline", "Outcome", "point.type", "point.est", "LB", "UB")


dt.imp.outcomes.melted <- dt.outcomes.melted[dt.bldisp, on = "baseline", nomatch = 0]

p.facet_sn_sp <- ggplot(dt.imp.outcomes.melted, aes(x = bl.disp, y = point.est, ymin = LB, ymax = UB, col = bl.disp))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  facet_grid(Outcome~point.type, scales="free", labeller = labeller(Outcome = Outcome.labs, point.type = point.type.labs))+
  theme(legend.position="none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
  scale_y_continuous(labels = scales::percent)+ylab("")+xlab("Baseline determination method")

ggsave("sens_spec_imp_fig.png", plot = p.facet_sn_sp, device = "png", width = 6, height = 7, units = "in")



#
# AUROC with CI with pROC package
# 

dt.AUROC.imp <- data.table(Outcome = character(),
                       baseline = character(),
                       AUROC = numeric(),
                       AUROC.LB = numeric(),
                       AUROC.UB = numeric())

#baselines_short <- c("Schwartz", "FAS.age", "FAS.height")
for (outcome in outcomes){
  for (baseline in bl.imp){
    if (outcome == "CKD3plus"){
      CIs <- eval(
        parse(text = paste0("ci.auc(roc(data = dt[CKD3plus != -1], response = ", 
                            outcome , ", predictor = raw.AKI.", baseline , "))")))
    } else {
      CIs <- eval(
        parse(text = paste0("ci.auc(roc(data = dt, response = ", 
                            outcome , ", predictor = raw.AKI.", baseline , "))")))
    }
    dt.AUROC.imp <- rbind(dt.AUROC.imp,
          list(outcome, baseline, CIs[2], CIs[1], CIs[3]))
  }
}


dt.AUROC.imp <- dt.AUROC.imp[dt.bldisp, on = "baseline", nomatch=0]

p.facet_AUROC_imp <- ggplot(dt.AUROC.imp, 
                            aes(x = bl.disp, 
                                y = AUROC, ymin = AUROC.LB, ymax = AUROC.UB, 
                                col = bl.disp))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  facet_grid(cols = vars(Outcome))+
  theme(legend.position="none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
  scale_y_continuous(labels = scales::percent)+ylab("")+xlab("")

ggsave("AUROC_imp_fig.png", plot = p.facet_AUROC_imp, device = "png", width = 6, height = 3, units = "in")
  
if(run_ROC == TRUE){
  
  #
  # ROC curves
  #
  
  n=300
  T.min = 0
  T.max = 15

  #
  # Coded AKI diagnosis
  #
  
  
  dt.ROC.code <- data.table(
    baseline = rep(bl.imp, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.imp)))
  
  
  dt.ROC.code[ , c("incidence", "Sensitivity", "Specificity") 
               := data.table(t(mapply(predict.outcome, 
                                       baseline = baseline,
                                       threshold = threshold, 
                                       outcome = "AKI.Diagnosis",
                                       full = FALSE)))]
  dt.ROC.code <- dt.ROC.code[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.code <-  
    ggplot(dt.ROC.code, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
      geom_line(size = 1)+
      stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
      scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      coord_fixed()+
      theme_light()+
      ylab("Sensitivity")+
      xlab("Specificity")+
      geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
      theme(legend.position =  "none")
  
  #
  # Post-hospitalization CKD
  #
  
  dt.ROC.CKD <- data.table(
    baseline = rep(bl.imp, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.imp)))
  
  
  dt.ROC.CKD[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "CKD3plus",
                                       full = FALSE)))]
  dt.ROC.CKD <- dt.ROC.CKD[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.CKD <- ggplot(dt.ROC.CKD, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")
  
  #
  # IP mortality
  #
  
  dt.ROC.Mort <- data.table(
    baseline = rep(bl.imp, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.imp)))
  
  
  dt.ROC.Mort[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "IP.Mortality",
                                       full = FALSE)))]
  dt.ROC.Mort <- dt.ROC.Mort[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.Mort <- ggplot(dt.ROC.Mort, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")
  
  
  dt.ROC.imp <- rbind(cbind(dt.ROC.code, Outcome = "AKI ICD code"),
                      cbind(dt.ROC.CKD, Outcome ="Post-hospitalization CKD"),
                      cbind(dt.ROC.Mort, Outcome = "Inpatient Mortality"))
  
  
  ROC.panel.imp <- ggplot(dt.ROC.imp, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    facet_grid(cols = vars(Outcome))+
    stat_function(fun = diag_fun, linetype = "dotted", size = 1, color = "darkred") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")+
    labs(color = "Baseline")
  

  ggsave("ROC_panel_imputed_fig.png", width = 6, height = 4, units = "in")
  
  
      

  ggsave("ROC_and_AUC_imp_fig.png", 
         grid.arrange(ROC.panel.imp, p.facet_AUROC_imp, nrow = 2), 
         width = 6, height = 5, units = "in")


}




```

```{r with_preadmit_analysis, include = FALSE}



dt.with_preSCr <- dt[n.prelabs > 0 ]
bl.preadmit <- c("bl.min","bl.second.min", "bl.prox", "bl.mean",
            "Schwartz.bedside.120", "FAS.age")


bl.empirical <- c("bl.min","bl.second.min", "bl.prox", "bl.mean")

# dt.small <- dt[sample(.N, 100)]


# PAIRWISE CORRELATION MATRIX




dt.bl.preadmit<- dt.with_preSCr[, ..bl.preadmit]

setnames(dt.bl.preadmit, old = dt.bldisp$baseline, new = as.character(dt.bldisp$bl.disp), skip_absent = TRUE)



# plt.pairs.imp <- ggpairs(data = dt.bl.preadmit,
#                     diag = list(continuous = diagfun),
#                     lower = list(continuous = lowerfun)
#                     ) + theme_bw()
# 
# ggsave("bl_pairwise_with_preSCr_fig.png", plt.pairs.imp, width = 8, height = 9, units = "in")

# MEAN AND STANDARD DEVIATION

dt.preSCr.bls <- melt(dt[ , ..bl.preadmit], measure.vars = bl.preadmit, variable.name = "baseline")[dt.bldisp, on = "baseline", nomatch = 0]

plt.preSCr.mean.sd.imputed <- 
  ggplot(dt.preSCr.bls, aes(col = bl.disp, x=bl.disp, y=value))+
  stat_summary(aes(col = bl.disp), fun.data=mean_sdl, geom="pointrange")+
  geom_violin(alpha = 0.3, col = NA, aes(fill = bl.disp)) + 
         ylab("Mean and standard deviation baseline SCr (mg/dL)")+
         xlab("Baseline determination method")+
  coord_flip(ylim = c(0, 1.5))+
  theme(legend.position = "none")

ggsave("baseline_mean_sd_preSCr_fig.png", plt.preSCr.mean.sd.imputed, width = 6, height = 4, units = "in")



# Mean and standard deviation GFR corresponding to baseline for empirical methods

dt.with_preSCr[ , paste0("gfr.", bl.empirical) := lapply(.SD, function(x) 0.413*Huse/x), .SDcols = bl.empirical]

gfr.bl <- paste0("gfr.", bl.empirical)

dt.bl.gfr.plt <- melt(dt.with_preSCr[ , ..gfr.bl])
dt.bl.gfr.plt <- dt.bl.gfr.plt[ , "baseline" := substr(variable, 5, length(variable))][dt.bldisp, on = "baseline", nomatch = 0]


p.bl.GFR <- ggplot(dt.bl.gfr.plt, aes(x=bl.disp, y=value)) + 
  geom_violin(alpha = 0.3, col = NA, aes(fill = bl.disp)) + 
  stat_summary(aes(col = bl.disp), fun.data=mean_sdl, geom="pointrange")+
  geom_hline(yintercept = c(107.3, 120), color = c("dodgerblue4", "darkred"), linetype = "dashed")+
  annotate("text", y = c(92, 130), x = c(4.5, 4.5), label =  c("107.3", "120"), color = c("dodgerblue4", "darkred"))+
  ylab(bquote('Baseline eGFR, mL/min/1.73m'^2))+xlab("Empirical baseline method")+
  scale_y_continuous() + coord_flip(ylim = c(0,300)) + theme(legend.position = "none")

ggsave("bl_gfr_empirical.png",p.bl.GFR, width = 6, height = 4, units = "in" )


#
# Incidence
#

dt.inc.preSCr <- data.table(
  baseline <- character(),
  inc <- numeric(),
  inc.LB <- numeric(),
  inc.UB <- numeric()
)

for(bl in bl.preadmit){
  dt.inc.preSCr <- 
    rbind(dt.inc.preSCr,
        incidence(bl, dt.with_preSCr), use.names = FALSE)
}

colnames(dt.inc.preSCr) <- c("baseline", "inc", "inc.LB", "inc.UB")

dt.inc.preSCr <- dt.inc.preSCr[dt.bldisp, on = "baseline", nomatch = 0]


p.inc.preSCr <- 
  ggplot(dt.inc.preSCr, aes(x=bl.disp, y=inc, ymin = inc.LB, ymax = inc.UB, col = bl.disp))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  theme(legend.position="none")+
  scale_y_continuous(labels = scales::percent)+
  ylab("Incidence")+
  xlab("Baseline Determination Method")+
  coord_flip()

ggsave("p_inc_preScr_fig.png", plot = p.inc.preSCr, device = "png", width = 5, height = 3, units = "in")


#dt.inc.imputed[ , baseline := factor(baseline, levels = c("bl.imp","bl.min","bl.second.min", "bl.prox", "bl.mean")) ]

#
# Clinical outcomes
#

dt.preSCr.outcomes <- data.table(
  baseline <- character(),
  outcome <- character(),
  threshold <- numeric(),
  TP <- numeric(),
  TN <- numeric(),
  FP <- numeric(),
  FN <- numeric(),
  Sn <- numeric(),
  Sn.LB <- numeric(),
  Sn.UB <- numeric(),
  Sp <- numeric(),
  Sp.LB <- numeric(),
  Sp.UB <- numeric()
)

colnames(dt.preSCr.outcomes) <- c("baseline", "Outcome", "Threshold", "TP", "TN", "FP", "FN", "Sn", "Sn.LB", "Sn.UB", "Sp", "Sp.LB", "Sp.UB")
  
outcomes <- c("AKI.Diagnosis", "CKD3plus", "IP.Mortality")
for(bl in bl.preadmit){
  for(outcome in outcomes){
     dt.preSCr.outcomes <- 
       rbind(dt.preSCr.outcomes,
             predict.outcome.PreSCr(bl, outcome), use.names = FALSE)
  }
}


# Outcome.labs <- c("AKI ICD Code", "Post-hospitalization CKD", "Inpatient Mortality")
# names(Outcome.labs) <- c("AKI.Diagnosis", "CKD3plus", "IP.Mortality")
# point.type.labs <- c("Sensitivity", "Specificity")
# names(point.type.labs) <- c("Sn", "Sp")
# 
# dt.outcomes.melted <- melt(dt.preSCr.outcomes, id.vars = c("baseline", "Outcome"),
#                            measure.vars = c("Sn", "Sp"),
#                            variable.name = "point.type",
#                            value.name = "point.est")
# 
# 
# dt.outcomes.melted.LB <- melt(dt.preSCr.outcomes, id.vars = c("baseline", "Outcome"),
#                            measure.vars = c("Sn.LB", "Sp.LB"),
#                            variable.name = "point.type",
#                            value.name = "LB")
# 
# dt.outcomes.melted.UB <- melt(dt.preSCr.outcomes, id.vars = c("baseline", "Outcome"),
#                            measure.vars = c("Sn.UB", "Sp.UB"),
#                            variable.name = "point.type",
#                            value.name = "UB")
# 
# dt.outcomes.melted <- cbind(dt.outcomes.melted, dt.outcomes.melted.LB$LB, dt.outcomes.melted.UB$UB)
# colnames(dt.outcomes.melted) <- c("baseline", "Outcome", "point.type", "point.est", "LB", "UB")
# 
# 
# dt.preSCr.outcomes.melted <- dt.outcomes.melted[dt.bldisp, on = "baseline", nomatch = 0]
# 
# p.facet_sn_sp <- ggplot(dt.preSCr.outcomes.melted, aes(x = bl.disp, y = point.est, ymin = LB, ymax = UB, col = bl.disp))+
#   geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
#   facet_grid(Outcome~point.type, scales="free", labeller = labeller(Outcome = Outcome.labs, point.type = point.type.labs))+
#   theme(legend.position="none",
#         axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
#   scale_y_continuous(labels = scales::percent)+ylab("")+xlab("Baseline determination method")
# 
# ggsave("sens_spec_imp_fig.png", plot = p.facet_sn_sp, device = "png", width = 6, height = 7, units = "in")


#
# AUROC with CI with pROC package
# 


dt.AUROC.preSCr <- data.table(Outcome = character(),
                       baseline = character(),
                       AUROC = numeric(),
                       AUROC.LB = numeric(),
                       AUROC.UB = numeric())




for (outcome in outcomes){
  for (baseline in bl.preadmit){
    if (outcome == "CKD3plus"){
      CIs <- eval(parse(text = paste0("ci.auc(roc(data = dt.with_preSCr[CKD3plus != -1], response = ", outcome , ", predictor = raw.AKI.", baseline , "))")))
    } else {
      CIs <- eval(parse(text = paste0("ci.auc(roc(data = dt.with_preSCr, response = ", outcome , ", predictor = raw.AKI.", baseline , "))")))
    }
    dt.AUROC.preSCr <- rbind(dt.AUROC.preSCr,
          list(outcome, baseline, CIs[2], CIs[1], CIs[3]))
  }
}


dt.AUROC.preSCr <- dt.AUROC.preSCr[dt.bldisp, on = "baseline", nomatch=0]

dt.AUROC.preSCr[ , Outcome := factor(
  ifelse(Outcome == "AKI.Diagnosis", "AKI ICD code", ifelse(Outcome == "IP.Mortality", "Inpatient mortality", "Post-hospitalization CKD"))
  , levels = c("AKI ICD code", "Post-hospitalization CKD", "Inpatient mortality"))]



p.facet_AUROC_preSCr <- ggplot(dt.AUROC.preSCr, aes(x = bl.disp, y = AUROC, ymin = AUROC.LB, ymax = AUROC.UB, col = bl.disp))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  facet_grid(cols = vars(Outcome))+
  theme(legend.position="none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
  scale_y_continuous(labels = scales::percent)+ylab("")+xlab("")+
  geom_hline(yintercept = 0.5, color = "darkred", linetype = "dotted", size = 1)+
  ylab("AUROC")

ggsave("AUROC_preSCr_fig.png", plot = p.facet_AUROC_preSCr, device = "png", width = 6, height = 3, units = "in")

#
# ROC curves
#

if(run_ROC == TRUE){
  n=401
  T.min = 0
  T.max = 8
  
  # Coded AKI diagnosis
  
  dt.ROC.code <- data.table(
    baseline = rep(bl.preadmit, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.preadmit)))
  
  
  dt.ROC.code[ , c("incidence", "Sensitivity", "Specificity") :=
                 data.table(t(mapply(predict.outcome.PreSCr, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "AKI.Diagnosis",
                                       full = FALSE)))]
  dt.ROC.code <- dt.ROC.code[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.code <-  
    ggplot(dt.ROC.code, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
      geom_line(size = 1)+
      stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
      scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      coord_fixed()+
      theme_light()+
      ylab("Sensitivity")+
      xlab("Specificity")+
      geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")
  
  # Post-hospitalization CKD
  
  dt.ROC.CKD <- data.table(
    baseline = rep(bl.preadmit, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.preadmit)))
  
  
  dt.ROC.CKD[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome.PreSCr, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "CKD3plus",
                                       full = FALSE)))]
  dt.ROC.CKD <- dt.ROC.CKD[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.CKD <- ggplot(dt.ROC.CKD, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")
  
  
  # IP mortality
  
  dt.ROC.Mort <- data.table(
    baseline = rep(bl.preadmit, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.preadmit)))
  
  
  dt.ROC.Mort[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome.PreSCr, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "IP.Mortality",
                                       full = FALSE)))]
  dt.ROC.Mort <- dt.ROC.Mort[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.Mort <- ggplot(dt.ROC.Mort, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")
  
  
  dt.ROC.preSCr <- rbind(cbind(dt.ROC.code, Outcome = "AKI ICD code"),
                      cbind(dt.ROC.CKD, Outcome ="Post-hospitalization CKD"),
                      cbind(dt.ROC.Mort, Outcome = "Inpatient mortality"))
  dt.ROC.preSCr[ , Outcome := factor(Outcome, levels = c("AKI ICD code", "Post-hospitalization CKD", "Inpatient mortality"))]
  
  
  ROC.panel.preSCr <- ggplot(dt.ROC.preSCr, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    facet_grid(cols = vars(Outcome))+
    stat_function(fun = diag_fun, linetype = "dotted", size = 1, color = "darkred") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")+
    labs(color = "Baseline")
  

  
  ggsave("ROC_panel_preSCr_fig.png", ROC.panel.preSCr, width = 6, height = 4, units = "in")
  
  
    dt.ROC.stages <- dt.ROC.preSCr[threshold %in% c(1.5, 2, 3), ][ , AKI_stage := ifelse(threshold == 1.5, "1+", ifelse(threshold == 2, "2+", "3"))]
    
    ROC_with_stages <- ggplot(dt.ROC.preSCr, 
                              aes(x = Specificity, y = Sensitivity, color = bl.disp))+
      geom_line(size = 1, alpha = 0.3)+
      facet_grid(rows = vars(Outcome))+
      stat_function(fun = diag_fun, linetype = "dotted", size = 1, color = "darkred") +
      scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      coord_fixed()+
      ylab("Sensitivity")+
      xlab("Specificity")+
      geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
      labs(color = "Baseline")+
      geom_point(data = dt.ROC.stages, 
                 aes(x = Specificity, y = Sensitivity, color = bl.disp, shape = AKI_stage), size = 4, alpha = 0.6)+
      theme(legend.position = "right")
    
    ggsave("ROC_with_stages.png", ROC_with_stages, 
           width = 6, height = 8.1, units = "in")
      

      ggsave("ROC_and_AUC_preSCr_fig.png", 
             grid.arrange(ROC.panel.preSCr, p.facet_AUROC_preSCr, nrow = 2), 
             width = 6, height = 5, units = "in")


    }










# Nadir inpatient GFR


dt.with_preSCr.GFR.melted <-melt(dt.with_preSCr,  id.vars = c("GFR.min"),
                           measure.vars = paste0("AKI.", bl.preadmit),
                           variable.name = "AKI_def",
                           value.name = "AKI_pos")


dt.with_preSCr.GFR.melted[ , AKI_pos := ifelse(AKI_pos == 1, "AKI+", "AKI-")]

dt.with_preSCr.GFR.melted <- dt.with_preSCr.GFR.melted[ AKI_pos == "AKI+"]

dt.with_preSCr.GFR.melted <- dt.with_preSCr.GFR.melted[dt.bldisp, on = "AKI_def", nomatch = 0]


p.nadir.GFR <- ggplot(dt.with_preSCr.GFR.melted, aes(x=bl.disp, y=GFR.min, color = bl.disp)) + 
  geom_violin(alpha = 0.3, col = NA, aes(fill = bl.disp))  + 
  stat_summary(fun.data=mean_sdl, geom="pointrange")+
  geom_hline(yintercept = 90, color = "red", linetype = "dashed")+
  annotate("text", y = 145, x = 6.3, label =  bquote('90 mL/min/1.73m'^2), color = "red")+
  ylab(bquote('Nadir eGFR, mL/min/1.73m'^2))+xlab("Baseline determination method")+
  coord_flip(ylim = c(0,200))+ theme(legend.position = "none")


ggsave("Nadir_gfr_fig.png", p.nadir.GFR, width = 4, height = 4, units = "in")

dt.frac.over.gfr <- dt.with_preSCr.GFR.melted[, sum(GFR.min > 90)/.N, by = bl.disp ]

```

```{r over12m, include = FALSE}


dt.with_preSCr_over12m <- dt[n.prelabs > 0 & age.yr >= 1 ]



#dt.inc.imputed[ , baseline := factor(baseline, levels = c("bl.imp","bl.min","bl.second.min", "bl.prox", "bl.mean")) ]

#
# Clinical outcomes
#
outcomes <- c("AKI.Diagnosis", "CKD3plus", "IP.Mortality")

dt.preSCr.outcomes <- data.table(
  baseline <- character(),
  outcome <- character(),
  threshold <- numeric(),
  TP <- numeric(),
  TN <- numeric(),
  FP <- numeric(),
  FN <- numeric(),
  Sn <- numeric(),
  Sn.LB <- numeric(),
  Sn.UB <- numeric(),
  Sp <- numeric(),
  Sp.LB <- numeric(),
  Sp.UB <- numeric()
)


if(run_ROC == TRUE){
  n=100
  T.min = 0
  T.max = 8
  
  # Coded AKI diagnosis
  
  dt.ROC.code <- data.table(
    baseline = rep(bl.preadmit, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.preadmit)))
  
  
  dt.ROC.code[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome.over12m, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "AKI.Diagnosis",
                                       full = FALSE)))]
  dt.ROC.code <- dt.ROC.code[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.code <-  
    ggplot(dt.ROC.code, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
      geom_line(size = 1)+
      stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
      scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      coord_fixed()+
      theme_light()+
      ylab("Sensitivity")+
      xlab("Specificity")+
      geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")
  
  # Post-hospitalization CKD
  
  dt.ROC.CKD <- data.table(
    baseline = rep(bl.preadmit, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.preadmit)))
  
  
  dt.ROC.CKD[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome.over12m, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "CKD3plus",
                                       full = FALSE)))]
  dt.ROC.CKD <- dt.ROC.CKD[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.CKD <- ggplot(dt.ROC.CKD, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")
  
  
  # IP mortality
  
  dt.ROC.Mort <- data.table(
    baseline = rep(bl.preadmit, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.preadmit)))
  
  
  dt.ROC.Mort[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome.over12m, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "IP.Mortality",
                                       full = FALSE)))]
  dt.ROC.Mort <- dt.ROC.Mort[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.Mort <- ggplot(dt.ROC.Mort, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")
  
  
  dt.ROC.preSCr <- rbind(cbind(dt.ROC.code, Outcome = "AKI ICD code"),
                      cbind(dt.ROC.CKD, Outcome ="Post-hospitalization CKD"),
                      cbind(dt.ROC.Mort, Outcome = "Inpatient mortality"))
  dt.ROC.preSCr[ , Outcome := factor(Outcome, levels = c("AKI ICD code", "Post-hospitalization CKD", "Inpatient mortality"))]
  
  
  ROC.panel.over12mo <- ggplot(dt.ROC.preSCr, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    facet_grid(cols = vars(Outcome))+
    stat_function(fun = diag_fun, linetype = "dotted", size = 1, color = "black", alpha = 0.4) +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")+
    labs(color = "Baseline")
  
}

#
# AUROC with CI with pROC package
# 


dt.AUROC.preSCr_over12m <- data.table(Outcome = character(),
                       baseline = character(),
                       AUROC = numeric(),
                       AUROC.LB = numeric(),
                       AUROC.UB = numeric())





for (outcome in outcomes){
  for (baseline in bl.preadmit){
    if (outcome == "CKD3plus"){
      CIs <- eval(parse(text = paste0("ci.auc(roc(data = dt.with_preSCr_over12m[CKD3plus != -1], response = ", outcome , ", predictor = raw.AKI.", baseline , "))")))
    } else {
      CIs <- eval(parse(text = paste0("ci.auc(roc(data = dt.with_preSCr_over12m, response = ", outcome , ", predictor = raw.AKI.", baseline , "))")))
    }
    dt.AUROC.preSCr_over12m <- rbind(dt.AUROC.preSCr_over12m,
          list(outcome, baseline, CIs[2], CIs[1], CIs[3]))
  }
}

dt.AUROC.preSCr_over12m <- dt.AUROC.preSCr_over12m[dt.bldisp, on = "baseline", nomatch=0]

dt.AUROC.preSCr_over12m[ , Outcome := factor(
  ifelse(Outcome == "AKI.Diagnosis", "AKI ICD code", ifelse(Outcome == "IP.Mortality", "Inpatient mortality", "Post-hospitalization CKD"))
  , levels = c("AKI ICD code", "Post-hospitalization CKD", "Inpatient mortality"))]



p.facet_AUROC_over12mo <- ggplot(dt.AUROC.preSCr_over12m, aes(x = bl.disp, y = AUROC, ymin = AUROC.LB, ymax = AUROC.UB, col = bl.disp))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  facet_grid(cols = vars(Outcome))+
  theme(legend.position="none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
  scale_y_continuous(labels = scales::percent)+ylab("")+xlab("")+
  geom_hline(yintercept = 0.5, color = "darkred", linetype = "dotted", size = 1)+
  ylab("AUROC")

# ggsave("AUROC_preSCr_fig.png", plot = p.facet_AUROC_preSCr, device = "png", width = 6, height = 3, units = "in")

if(run_ROC == TRUE){
  ggsave("ROC_and_AUC_over12mo_fig.png", 
         grid.arrange(ROC.panel.over12mo, p.facet_AUROC_over12mo, nrow = 2), 
         width = 6, height = 5, units = "in")
}



```

```{r empirical_fist_analysis, include = FALSE}
bl.emp_first <- c("bl.min", "bl.mean", "bl.prox", "bl.second.min", 
                  "Schwartz.bedside.120", "FAS.age",
                  "lowest.FAS", "second.lowest.FAS", "mean.FAS", "prox.FAS")


#
# Incidence
#

dt.inc.empfirst <- data.table(
  baseline <- character(),
  inc <- numeric(),
  inc.LB <- numeric(),
  inc.UB <- numeric()
)

for(bl in bl.emp_first){
  dt.inc.empfirst <- 
    rbind(dt.inc.empfirst,
        incidence(bl, dt), use.names = FALSE)
}

colnames(dt.inc.empfirst) <- c("baseline", "inc", "inc.LB", "inc.UB")

dt.inc.empfirst <- dt.inc.empfirst[dt.bldisp, on = "baseline", nomatch = 0]


p.inc.empfirst <- 
  ggplot(dt.inc.empfirst, aes(x=bl.disp.first, y=inc, ymin = inc.LB, ymax = inc.UB, col = empiric.type))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  theme(legend.position="none")+
  scale_y_continuous(labels = scales::percent)+
  ylab("Incidence")+
  xlab("Empirical-first method")+
  coord_flip()

ggsave("p_inc_empfirst_fig.png", plot = p.inc.empfirst, device = "png", width = 6, height = 3, units = "in")


#dt.inc.imputed[ , baseline := factor(baseline, levels = c("bl.imp","bl.min","bl.second.min", "bl.prox", "bl.mean")) ]

#
# Clinical outcomes
#

dt.empfirst.outcomes <- data.table(
  baseline <- character(),
  outcome <- character(),
  threshold <- numeric(),
  TP <- numeric(),
  TN <- numeric(),
  FP <- numeric(),
  FN <- numeric(),
  Sn <- numeric(),
  Sn.LB <- numeric(),
  Sn.UB <- numeric(),
  Sp <- numeric(),
  Sp.LB <- numeric(),
  Sp.UB <- numeric()
)

colnames(dt.empfirst.outcomes) <- c("baseline", "Outcome", "Threshold", "TP", "TN", "FP", "FN", "Sn", "Sn.LB", "Sn.UB", "Sp", "Sp.LB", "Sp.UB")
  
outcomes <- c("AKI.Diagnosis", "CKD3plus", "IP.Mortality")
for(bl in bl.emp_first){
  for(outcome in outcomes){
     dt.empfirst.outcomes <- 
       rbind(dt.empfirst.outcomes,
             predict.outcome(bl, outcome), use.names = FALSE)
  }
}


Outcome.labs <- c("AKI ICD Code", "Post-hospitalization CKD", "Inpatient Mortality")
names(Outcome.labs) <- c("AKI.Diagnosis", "CKD3plus", "IP.Mortality")
point.type.labs <- c("Sensitivity", "Specificity")
names(point.type.labs) <- c("Sn", "Sp")

dt.outcomes.melted <- melt(dt.empfirst.outcomes, id.vars = c("baseline", "Outcome"),
                           measure.vars = c("Sn", "Sp"),
                           variable.name = "point.type",
                           value.name = "point.est")


dt.outcomes.melted.LB <- melt(dt.empfirst.outcomes, id.vars = c("baseline", "Outcome"),
                           measure.vars = c("Sn.LB", "Sp.LB"),
                           variable.name = "point.type",
                           value.name = "LB")

dt.outcomes.melted.UB <- melt(dt.empfirst.outcomes, id.vars = c("baseline", "Outcome"),
                           measure.vars = c("Sn.UB", "Sp.UB"),
                           variable.name = "point.type",
                           value.name = "UB")

dt.empfirst.outcomes.melted <- cbind(dt.outcomes.melted, dt.outcomes.melted.LB$LB, dt.outcomes.melted.UB$UB)
colnames(dt.empfirst.outcomes.melted) <- c("baseline", "Outcome", "point.type", "point.est", "LB", "UB")


dt.empfirst.outcomes.melted <- dt.empfirst.outcomes.melted[dt.bldisp, on = "baseline", nomatch = 0]

p.facet_sn_sp <- ggplot(dt.empfirst.outcomes.melted, aes(x = bl.disp.first, y = point.est, ymin = LB, ymax = UB, col = empiric.type))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  facet_grid(Outcome~point.type, scales="free", labeller = labeller(Outcome = Outcome.labs, point.type = point.type.labs))+
  theme(legend.position="none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
  scale_y_continuous(labels = scales::percent)+ylab("")+xlab("Empirical-first method")

ggsave("sens_spec_empfirst_fig.png", plot = p.facet_sn_sp, device = "png", width = 6, height = 7, units = "in")





#
# AUROC with CI with pROC package
# 

dt.AUROC.emp_first <- data.table(Outcome = character(),
                       baseline = character(),
                       AUROC = numeric(),
                       AUROC.LB = numeric(),
                       AUROC.UB = numeric())



#baselines_short <- c("Schwartz", "FAS.age", "FAS.height")
for (outcome in outcomes){
  for (baseline in bl.emp_first){
    if (outcome == "CKD3plus"){
      CIs <- eval(
        parse(text = paste0("ci.auc(roc(data = dt[CKD3plus != -1], response = ", 
                            outcome , ", predictor = raw.AKI.", baseline , "))")))
    } else {
      CIs <- eval(
        parse(text = paste0("ci.auc(roc(data = dt, response = ", 
                            outcome , ", predictor = raw.AKI.", baseline , "))")))
    }
    dt.AUROC.emp_first <- rbind(dt.AUROC.emp_first,
          list(outcome, baseline, CIs[2], CIs[1], CIs[3]))
  }
}


dt.AUROC.emp_first <- dt.AUROC.emp_first[dt.bldisp, on = "baseline", nomatch=0]

dt.AUROC.emp_first[ , Outcome := factor(
  ifelse(Outcome == "AKI.Diagnosis", "AKI ICD code", ifelse(Outcome == "IP.Mortality", "Inpatient mortality", "Post-hospitalization CKD"))
  , levels = c("AKI ICD code", "Post-hospitalization CKD", "Inpatient mortality"))]

p.facet_AUROC_emp_first <- ggplot(dt.AUROC.emp_first, 
                            aes(x = bl.disp, 
                                y = AUROC, ymin = AUROC.LB, ymax = AUROC.UB, 
                                col = bl.disp))+
  geom_pointrange(position=position_dodge(width=0.30), shape = 4)+
  facet_grid(cols = vars(Outcome))+
  theme(legend.position="none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
  scale_y_continuous(labels = scales::percent)+ylab("")+xlab("")

ggsave("AUROC_emp_first_fig.png", plot = p.facet_AUROC_emp_first, device = "png", width = 6, height = 3, units = "in")
  
if(run_ROC == TRUE){
  
  #
  # ROC curves
  #
  
  n=300
  T.min = 0
  T.max = 15

  #
  # Coded AKI diagnosis
  #
  
  
  dt.ROC.code <- data.table(
    baseline = rep(bl.emp_first, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.emp_first)))
  
  
  dt.ROC.code[ , c("incidence", "Sensitivity", "Specificity") 
               := data.table(t(mapply(predict.outcome, 
                                       baseline = baseline,
                                       threshold = threshold, 
                                       outcome = "AKI.Diagnosis",
                                       full = FALSE)))]
  dt.ROC.code <- dt.ROC.code[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.code <-  
    ggplot(dt.ROC.code, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
      geom_line(size = 1)+
      stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
      scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
      coord_fixed()+
      theme_light()+
      ylab("Sensitivity")+
      xlab("Specificity")+
      geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
      theme(legend.position =  "none")
  
  #
  # Post-hospitalization CKD
  #
  
  dt.ROC.CKD <- data.table(
    baseline = rep(bl.emp_first, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.emp_first)))
  
  
  dt.ROC.CKD[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "CKD3plus",
                                       full = FALSE)))]
  dt.ROC.CKD <- dt.ROC.CKD[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.CKD <- ggplot(dt.ROC.CKD, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")
  
  #
  # IP mortality
  #
  
  dt.ROC.Mort <- data.table(
    baseline = rep(bl.emp_first, each = n),
    threshold = rep(seq(from=T.min, to =T.max, by = (T.max-T.min)/(n-1) ), length(bl.emp_first)))
  
  
  dt.ROC.Mort[ , c("incidence", "Sensitivity", "Specificity") := data.table(t(mapply(predict.outcome, 
                                       baseline = baseline, 
                                       threshold = threshold, 
                                       outcome = "IP.Mortality",
                                       full = FALSE)))]
  dt.ROC.Mort <- dt.ROC.Mort[dt.bldisp, on = "baseline", nomatch=0]
  
  ROC.Mort <- ggplot(dt.ROC.Mort, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    stat_function(fun = diag_fun, linetype = "dashed", size = 1, color = "black") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    theme_light()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")
  
  
  dt.ROC.emp_first <- rbind(cbind(dt.ROC.code, Outcome = "AKI ICD code"),
                      cbind(dt.ROC.CKD, Outcome ="Post-hospitalization CKD"),
                      cbind(dt.ROC.Mort, Outcome = "Inpatient Mortality"))
  
  dt.ROC.emp_first[ , Outcome := factor(Outcome,
                                           levels = c("AKI ICD code", 
                                                      "Post-hospitalization CKD",
                                                      "Inpatient Mortality"))]
  ROC.panel.emp_first <- ggplot(dt.ROC.emp_first, aes(x = Specificity, y = Sensitivity, color = bl.disp))+
    geom_line(size = 1)+
    facet_grid(cols = vars(Outcome))+
    stat_function(fun = diag_fun, linetype = "dotted", size = 1, color = "darkred") +
    scale_x_reverse(limits = c(1,0), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    scale_y_continuous(limits = c(0,1), breaks = c(0, .2, .4, .6, .8, 1), minor_breaks = NULL)+
    coord_fixed()+
    ylab("Sensitivity")+
    xlab("Specificity")+
    geom_rect(xmin = -1, xmax = 0, ymin = 0, ymax = 1, fill = NA, color = "black")+
    theme(legend.position =  "none")+
    labs(color = "Baseline")
  

  ggsave("ROC_panel_emp_first_fig.png", width = 6, height = 4, units = "in")
  
  
      

  ggsave("ROC_and_AUC_emp_first_fig.png", 
         grid.arrange(ROC.panel.emp_first, p.facet_AUROC_emp_first, nrow = 2), 
         width = 6, height = 5, units = "in")


}



```

``` {r incidence_by_demo, include = FALSE}
dt[, age.bin := ifelse(age.yr < 1, "3 months to 1 year",
                       ifelse(age.yr < 6, "1 to 5.9 years",
                              ifelse(age.yr < 13, "6 to 12.9 years",
                                     "13 to 18 years")))]

AKI.all <- paste0("AKI.",bl.emp_first)



dt.inc_by_demo <- data.table(
  Type =character(),
  Group = character(),
  N = numeric(),
  AKI.Schwartz.bedside.120 = numeric(),
  AKI.bl.min = numeric(),
  AKI.bl.second.min = numeric(),
  AKI.bl.mean = numeric(),
  AKI.bl.prox = numeric(),
  AKI.FAS.age = numeric(),
  AKI.lowest.FAS = numeric(),
  AKI.second.lowest.FAS = numeric(),
  AKI.mean.FAS = numeric(),
  AKI.prox.FAS = numeric()
)

AKI.all <- colnames(dt.inc_by_demo)[4:13]

dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "",  Group = "All", 
                          N = dt[ , .N]),
                          dt[ , lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )
#SEX
dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Sex",  Group = "Female", 
                          N = dt[Sex.Title =="FEMALE", .N]),
                          dt[Sex.Title =="FEMALE", lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )

dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Sex",  Group = "Male", 
                          N = dt[Sex.Title =="MALE", .N]),
                          dt[Sex.Title =="MALE", lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )

# RACE
dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Primary race",  Group = "Asian", 
                          N = dt[Race.Title =="ASIAN", .N]),
                          dt[Race.Title =="ASIAN", lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )

dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Primary race",  Group = "Black", 
                          N = dt[Race.Title =="BLACK OR AFRICAN AMERICAN", .N]),
                          dt[Race.Title =="BLACK OR AFRICAN AMERICAN", lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )

dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Primary race",  Group = "Native American", 
                          N = dt[Race.Title =="AMERICAN INDIAN OR ALASKA NATIVE", .N]),
                          dt[Race.Title =="AMERICAN INDIAN OR ALASKA NATIVE", lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )

dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Primary race",  Group = "Pacific islander", 
                          N = dt[Race.Title =="NATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER", .N]),
                          dt[Race.Title =="NATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER", lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )


dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Primary race",  Group = "White", 
                          N = dt[Race.Title =="WHITE OR CAUCASIAN", .N]),
                          dt[Race.Title =="WHITE OR CAUCASIAN", lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )





other_races <- c("OTHER", "DECLINES TO STATE","UNKNOWN"  )
dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Primary race",  Group = "Other", 
                          N = dt[Race.Title %in% other_races, .N]),
                          dt[Race.Title %in% other_races, lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )
# Age

for (bin in c("3 months to 1 year", "1 to 5.9 years", "6 to 12.9 years", "13 to 18 years")){
  dt.inc_by_demo <- rbind(dt.inc_by_demo,
                        append(list(
                          Type = "Age",  Group = bin, 
                          N = dt[age.bin == bin, .N]),
                          dt[age.bin == bin, lapply(.SD, function(x) 100*sum(x)/.N), .SDcols =AKI.all]
                          )
                        )
}


dt.inc_by_demo[ , group_N := paste0(Group, " (", trimws(format(N, big.mark = ',')), ")")]

dt.inc_by_demo[, c("Group", "N") := NULL]
setcolorder(dt.inc_by_demo, neworder = c("Type", "group_N"))



col_labels <- c("Group (N)", "Imputed", "Lowest", "Second lowest", "Mean", "Proximal",
                "Imputed", "Lowest", "Second lowest", "Mean", "Proximal")
names(col_labels) <- c("group_N", AKI.all)

## FORMAT FLEXTABLE

#Incidence.by.demo <- fread("incidence_by_demo.csv")
dt.inc_by_demo <- as_grouped_data(dt.inc_by_demo, groups = c("Type"))


t1 <- as_flextable(dt.inc_by_demo[-1,])
t1 <- compose(t1, i = ~ !is.na(Type), j = "group_N",
              value = as_paragraph(as_chunk(Type)))



t1 <- set_header_labels(t1, values = col_labels)
t1 <- add_header_row(t1, values = c("Group (N)", 
                                    rep("Schwartz 120 imputation", 5),
                                    rep("FAS imputation", 5)))
t1 <- merge_h(t1, part = "header")
t1 <- merge_v(t1, part = "header")

t1 <- theme_box(t1)
t1 <- fontsize(t1, size = 10, part = "all")
t1 <- font(t1, fontname = "Times", part = "all")
t1 <- align(t1, align = "left", part = "body")
t1 <- align(t1, align = "center", part = "header")

t1 <- colformat_num(t1, j = 2:11, digits = 1, suffix = "%")

t1 <- bold(t1, j=1, i = ~ !is.na(Type), bold = TRUE, part = "body")
t1 <- width(t1, j = 1, width = 1.5)
t1 <- height(t1, height = 0.4, part = "body")
t1 <- height(t1, height = 0.4, part = "header")
t1 <- valign(t1, valign = "top", part = "all")
t1 <- align(t1, i = c(2, 5, 12), align = "left", part = "body")
t1 <- bg(t1, i = c(2, 5, 12), bg = "#DDDDDD", part = "body")
t1 <- width(t1, j = 2:11, width = 0.75)
t1 <- width(t1, j = 1, width = 1.5)

```

``` {r AKI_diag_by_stage, include = FALSE}

dt.AKI_diag_by_stage <- data.table(
  baseline = bl.emp_first
)
dt.AKI_diag_by_stage[ , `AKI stage 1+` := lapply(baseline, function(X) pct.diagnosed(dt, 1, X))]
dt.AKI_diag_by_stage[ , `AKI stage 2+` := lapply(baseline, function(X) pct.diagnosed(dt, 2, X))]
dt.AKI_diag_by_stage[ , `AKI stage 3` := lapply(baseline, function(X) pct.diagnosed(dt, 3, X))]

dt.AKI_diag_by_stage <- dt.AKI_diag_by_stage[dt.bldisp[ , c("baseline", "bl.disp.first")], on = "baseline", nomatch = 0]
dt.AKI_diag_by_stage[, baseline := NULL]

setcolorder(dt.AKI_diag_by_stage, "bl.disp.first")
colnames(dt.AKI_diag_by_stage)[1] <- "Baseline method"

fwrite(dt.AKI_diag_by_stage, file = "AKI_diag_by_stage.csv")


```

# Abstract

**Background:** The current consensus definition for acute kidney injury (AKI) does not specify how baseline serum creatinine should be determined. Studies have used heterogeneous approaches. We assessed how baseline determination impacted AKI incidence and the association between AKI and clinical outcomes.

**Methods:** We retrospectively applied empirical (based on measured serum creatinine) and imputed (based on age or height) baseline estimation methods to pediatric patients discharged between 2014 and 2019 from an academic hospital. Using each method, we estimated AKI incidence and assessed the area under the ROC curve (AUROC) for AKI as a predictor of three clinical outcomes: application of an AKI billing code (proxy for more clinically overt disease), inpatient mortality, and post-hospitalization chronic kidney disease.

**Results:** Incidence was highly variable across baseline methods (12.2% – 26.7%). Incidence was highest when the lowest pre-admission creatinine was used if available and the Schwartz bedside equation was used to impute one otherwise. AKI was more predictive of application of an AKI billing code when baseline was imputed universally, regardless of pre-admission values (AUROC 80.7% – 84.9%) than it was with any empirical approach (AUROC 64.5% – 76.6%). AKI was predictive of post-hospitalization CKD when using universal imputation baseline methods (AUROC 67.0% – 74.6%); AKI was not strongly predictive of post-hospitalization CKD when using empirical baseline methods (AUROC 46.4% – 58.5%).  Baseline determination method did not affect the association between AKI and inpatient mortality.

**Conclusions:** Method of baseline determination influences AKI incidence and the association between AKI and clinical outcomes, illustrating the need for standard criteria. Imputing the baseline for all patients, even when preadmission creatinine is available, may identify a more clinically relevant subset of the disease.

#####

# Key words

Acute kidney injury, AKI, children, serum creatinine, epidemiology

<br> 

# Declarations

**Funding:** This research and manuscript submission were unfunded.

**Conflicts:** DS is an advisor to Carta Healthcare, a healthcare analytics company. WAR and SMS report no relevant conflicts of interest.

**Ethics/Consent:** Research using retrospective patient data obtained through routine care performed with waiver of consent.

**Data and materials:** Not applicable.

**Code availability:** Upon acceptance we will upload all code to an open, DOI-indexed online repository.

**Authors’ contributions:** WAR collected and analyzed the data and drafted and revised the manuscript. SMS conceived of the analysis, advised on methods, and critically revised the manuscript. DS advised on methods and critically revised the manuscript.


#####

# Introduction

Acute kidney injury (AKI) is a common complication among hospitalized patients that is associated with increased morbidity and mortality [@Kaddourah2017; @Coca2012; @Sutherland2013; @Sutherland2015; @Hoste2015]. AKI is currently defined according to the Kidney Disease: Improving Global Outcomes (KDIGO) criteria, which diagnose and stage AKI based upon increases in serum creatinine (creatinine) and decrements in urine output (UOP) [@KDIGO2012]. KDIGO defines AKI as an increase in creatinine of at least 1.5 times the baseline value, but it does not unambiguously specify a baseline determination method. As a result, application of the baseline concept varies between studies and across organizations. In practice, one of the most commonly employed methods is an "empirical-first" approach where preadmission creatinine is used as the baseline if measured and imputation is performed if not. Although inconsistency in baseline determination method has been shown to influence estimates of AKI incidence, it is not currently known whether the baseline determination method affects the association between AKI and clinical outcomes [@Siew2015; @Siew2010]. We analyzed the impact of baseline creatinine determination on the relationship between AKI and inpatient mortality, application of an AKI billing code (as a proxy for more severe, clinically overt AKI episodes), and the subsequent development of chronic kidney disease (CKD).  We hypothesized that the method used to determine a patient’s baseline creatinine would affect the relationship between the development of AKI and subsequent AKI-related outcomes and also generate heterogeneous epidemiologic findings.

#####

# Methods

This was a retrospective study of patients admitted to Lucile Packard Children’s Hospital at Stanford University between May 2014 and May 2019. Included patients were between 3 months and 18 years of age at hospital admission.  Patients were excluded due to the absence of a creatinine measurement in the first seven hospital days and if they had an ICD diagnosis code for CKD, end stage renal disease (ESRD), or kidney transplantation during or before the index hospitalization. Following the KDIGO definition, patients were considered AKI positive during the inpatient stay if creatinine increased at least 150% over baseline during the first seven days of hospitalization (relative increase criterion) or if creatinine increased at least 0.3 mg/mL within 48 hours during hospitalization (absolute increase criterion) [@KDIGO2012].

<br>

Four baseline imputation methods and four empirical baseline methods were evaluated (see Supplementary Table S1 for equations). Two of the imputed methods back-calculated creatinine from the original Schwartz bedside equation [@Schwartz2009] by assuming a reference GFR. The **Schwartz 120** method used  a reference GFR of 120 mL/min/1.73 m^2^, as is common practice [@Kaddourah2017; @Zappitelli2008], and the **Schwartz 107** method used a reference GFR of 107.3 mL/min/1.73 m^2^, the median GFR among healthy children age 4-18 [@Pottel2012]. The third and fourth imputation methods relied upon the pediatric full age spectrum (FAS) eGFR equations which assign a baseline based on age and sex (**FAS** method) or based on height alone (**FAS height**); both FAS eGFR equations assume a normal baseline GFR of 107.3 mL/min/1.73 m^2^ [@Hoste2014]. The four empirical baseline methods used pre-admission creatinines measured in the three months prior to admission: using the lowest pre-admission value available (**lowest**), using the second lowest pre-admission value available (**second lowest**), using the value most proximal to admission date (**proximal**), and using the mean value of all creatinine values available in the three months prior to admission (**mean**). If no height was recorded within 90 days of admission, a sex-specific polynomial regression was used to impute the patient’s height as described in the online supplement. Because the empirical methods require a preadmission creatinine, we also evaluated "empirical-first" versions of each of the empirical methods, where imputation was performed for any patients who did not have a preadmission creatinine value available, as is commonly done in practice [@Kaddourah2017].

<br>

## Analysis of universal baseline imputation methods

The first analysis was an evaluation of universal imputation using each of the four baseline imputation methods. The distribution of baseline creatinine resulting from each imputation method and the correlation between methods were assessed. Incidence of AKI was estimated using each imputation method. The predictiveness of AKI for three clinical outcomes -- inpatient mortality, receipt of an AKI International Classification of Diseases (ICD) 9/10 billing code, and post-hospitalization CKD -- was assessed for each imputation method. Post-hospitalization CKD was defined as an eGFR below 60 mL/min/1.73m^2^ on two separate occasions between 90 and 450 days after discharge. Patients' eGFR was calculated using the bedside Schwartz equation [@Schwartz2009], and only patients with at least two post-discharge creatinine measures were included in the CKD analysis. To assess predictive performance, sensitivity and specificity were calculated using the standard KDIGO definition, which uses a relative increase threshold of at least 50% over baseline within the first 7 days of hospitalization. To gain a more comprehensive understanding of the impact of baseline determination method, receiver operator characteristic (ROC) analysis was conducted by varying the relative increase threshold to investigate the predictiveness of AKI for clinical outcomes across the different baseline methods (rationale for ROC analysis is further discussed in the Online Supplement).


<br>

## Comparison of empirical and imputed methods for patients with preadmission creatinine

In the second analysis, the four empirical baseline methods were compared to two of the most commonly applied imputation methods-- the Schwartz 120 and FAS methods. Since empirical methods cannot be applied to patients who are missing a preadmission creatinine, this analysis was restricted to the subset of patients for whom creatinine was measured in the 90-day pre-admission period. Baseline distribution, correlation across baseline methods, and ability to predict clinical outcomes were analyzed. Because patients with preadmission creatinine do not reflect a general patient population, AKI incidence and sensitivity and specificity for clinical outcomes were not reported for this cohort. Nadir inpatient GFR for patients classified as AKI-positive was compared across baseline determination methods.

<br>

## Comparison of "empirical-first" methods and universal imputation

Two universal imputation methods (**Schwartz 120** and **FAS**) were compared to eight empirical-first methods that combine an empirical methods (**lowest**, **second lowest**, **mean**, or **proximal**) for patients with measured preadmission creatinine and an imputed method (Schwartz 120 or FAS) for other patients. These empirical-first methods were"" **lowest-Schwartz**, **second lowest-Schwartz**, **mean-Schwartz**, **proximal-Schwartz**, and then **lowest-FAS**, **second lowest-FAS**, **mean-FAS**, and **proximal-FAS**. AKI incidence as well as sensitivity and specificity for clinical outcomes were computed for each  baseline method. Incidence was also analyzed by age, sex, and race.


#####

# Results


## Patient population

`r format(n_pt, big.mark = ',')` patients were included in the analysis, of which `r format(n_pt_CKD, big.mark = ',')` (`r sprintf("%.1f%%", 100*n_pt_CKD/n_pt)`) patients had two post-discharge creatinine measures and were included in the post-hospitalization CKD analysis (Figure 1). `r sprintf("%.1f%%", 100* percent.with.prelabs)` of patients had at least one creatinine measured in the three months prior to admission and `r sprintf("%.1f%%", 100*percent.mult.prelabs)` of patients had more than one. `r format(dt[ Hsource == "Measured", .N], big.mark = ',')` (`r sprintf("%.1f%%", 100* dt[ Hsource == "Measured", .N]/dt[,.N])`) patients had a height measurement within 3 months of admission. Height was imputed for the `r format(dt[ Hsource == "Imputed", .N], big.mark = ',')` (`r sprintf("%.1f%%", 100* dt[ Hsource == "Imputed", .N]/dt[,.N])`) patients with no height measurement on record and for the `r format(dt[ Hsource == "Adjusted", .N], big.mark = ',')` (`r sprintf("%.1f%%", 100* dt[ Hsource == "Adjusted", .N]/dt[,.N])`) patients whose nearest height measurement was more than 3 months from their admission date. `r format(dt[ , sum(AKI.48)], big.mark = ',')` (`r sprintf("%.1f%%", 100*dt[ , sum(AKI.48)/.N])`) patients had a documented increase in creatinine of more than 0.3 mg/dL within 48 hours during their stay, meeting the absolute increase criterion of the KDIGO AKI definition. In the ROC analyses presented below, these patients were classified as AKI-positive regardless of relative increase threshold or baseline determination method. Across the entire cohort, the inpatient mortality rate was `r sprintf("%.1f%%", 100*sum(dt$IP.Mortality)/nrow(dt))`, and an AKI billing code was applied to `r sprintf("%.1f%%", 100*sum(dt$AKI.Diagnosis)/nrow(dt))` of the cohort.  Among the `r format(n_pt_CKD, big.mark = ',')` patients for whom two post-hospitalization creatinine measurements were available, the rate of post-hospitalization CKD was `r sprintf("%.1f%%", 100*sum(dt$CKD3plus == 1)/nrow(dt[CKD3plus != -1]))`.

<br>

## Analysis of universal imputation baseline determination methods

Of the universal imputation methods, the FAS method led to the largest mean  and standard deviation baseline creatinine (mean = `r format(mean(dt$FAS.age), big.mark = ',', digits = 2)` mg/dL; standard deviation `r format(dt[ , sd(FAS.age)], big.mark = ',', digits = 2)` mg/dL) and the Schwartz method with a reference GFR of 120 mL/min/1.73 m^2^ led to the smallest (mean `r format(mean(dt$Schwartz.bedside.120), big.mark = ',', digits = 2)` mg/dL; standard deviation `r format(dt[ , sd(Schwartz.bedside.120)], big.mark = ',', digits = 2)` mg/dL; Figure S1). The FAS height method was the least correlated with the other methods (correlation coefficient of 0.885 with the FAS method; 0.887 with the Schwartz methods; Figure S2). The FAS height method produced higher baseline creatinine values than the other methods for patients shorter than 75 cm (`r sprintf("%.1f%%", 100* dt[ Huse < 75, .N]/dt[,.N])` of patients) and patients taller than 175 cm (`r sprintf("%.1f%%", 100* dt[ Huse > 175, .N]/dt[,.N])` of patients; Figure S3). Incidence was similar for the FAS method (`r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "FAS.age", inc] )`, 95% CI `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "FAS.age", inc.LB] )` -- `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "FAS.age", inc.UB] )`), FAS height method (`r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "FAS.height", inc] )`, 95% CI `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "FAS.height", inc.LB] )` -- `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "FAS.height", inc.UB] )`), and Schwartz 107 method (`r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "Schwartz.bedside.107", inc] )`, 95% CI `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "Schwartz.bedside.107", inc.LB] )` -- `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "Schwartz.bedside.107", inc.UB] )`), all of which use the same reference GFR of 107 mL/min/1.73 m^2^. Incidence was significantly higher using the Schwartz method with a reference GFR of 120 mL/min/1.73 m^2^ (`r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "Schwartz.bedside.120", inc] )`, 95% CI `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "Schwartz.bedside.120", inc.LB] )` -- `r sprintf("%.1f%%", 100*dt.inc.imputed[baseline == "Schwartz.bedside.120", inc.UB] )`, Figure S4). In assessing how baseline imputation affected AKI as a predictor for the three clinical outcomes (application of an AKI ICD billing code, post-hospitalization CKD, and inpatient mortality), sensitivity was highest and specificity was lowest using the Schwartz 120 method, although this difference was not statistically significant in the case of sensitivity for inpatient mortality (Figure S5). Differences in AUROC were not statistically significant for all three outcomes (Figure S6).

<br>

## Comparison of empirical and universal imputation baseline determination methods among patients with measured preadmission creatinine

For the `r format(dt[n.prelabs>0, .N], big.mark = ',')` patients with one or more pre-admission creatinine available, correlation was high between all empirical methods (at least 0.82) but lower between empirical and imputed methods (coefficients between 0.50 and 0.57; Figure S7). Empirical methods led to larger standard deviations in baseline creatinine estimates compared to the imputed methods (Figure 2). While the imputed baseline methods implied a baseline eGFR of 107.3 or 120 mL/min/1.73 m^2^ by design, empirical baselines were associated with a wide range of baseline GFRs, particularly when the lowest pre-admission creatinine was used (Figure S8). Of the empirical methods, using the lowest pre-admission creatinine yielded the largest baseline eGFR on average (mean `r format(mean(dt.with_preSCr$gfr.bl.min), big.mark = ',', digits = 3)` mL/min/1.73 m^2^; standard deviation `r format(dt.with_preSCr[ , sd(gfr.bl.min)], big.mark = ',', digits = 2)`), and the mean method yielded the smallest (mean `r format(mean(dt.with_preSCr$gfr.bl.mean), big.mark = ',', digits = 3)` mL/min/1.73 m^2^; standard deviation `r format(dt.with_preSCr[ , sd(gfr.bl.mean)], big.mark = ',', digits = 2)`; Figure S8). Incidence was negatively correlated with average baseline creatinine; baseline methods that yielded lower baseline creatinine values yielded higher incidences. As an exception, the FAS method had the third-lowest incidence (Figure S9) despite yielding the highest baseline creatinine (Figure 2). 

<br>

In ROC analysis, AKI was more predictive of CKD and receipt of an ICD-coded diagnoses when the two universal imputation methods (FAS and Schwartz 120) were employed compared to when any of the empirical methods were used (Figure 3). Notably, AKI was weakly predictive, or not predictive at all, for CKD using the empirical baseline methods (AUROC < 60% for all methods; AUROC = `r sprintf("%.1f%%", 100* dt.AUROC.preSCr[Outcome == "Post-hospitalization CKD" & baseline == "bl.prox"]$AUROC)` [95% CI `r sprintf("%.1f%%", 100* dt.AUROC.preSCr[Outcome == "Post-hospitalization CKD" & baseline == "bl.prox"]$AUROC.LB)` -- `r sprintf("%.1f%%", 100* dt.AUROC.preSCr[Outcome == "Post-hospitalization CKD" & baseline == "bl.prox"]$AUROC.UB)`] using the proximal method) There was no significant difference in association between AKI and inpatient mortality across any of the imputed or empirical approaches.

<br> 

The distribution of nadir inpatient eGFR for AKI-positive patients varied across baseline determination methods (Figure 4). Many patients classified as AKI-positive never had an eGFR below 90 mL/min/1.73 m^2^ during their inpatient stay using the lowest (`r sprintf("%.1f%%", 100*dt.frac.over.gfr[bl.disp == "Lowest", V1])`) and second lowest (`r sprintf("%.1f%%", 100*dt.frac.over.gfr[bl.disp == "Second Lowest", V1])`) baseline method. This proportion was considerably lower using the mean (`r sprintf("%.1f%%", 100*dt.frac.over.gfr[bl.disp == "Mean", V1])`) or proximal (`r sprintf("%.1f%%", 100*dt.frac.over.gfr[bl.disp == "Proximal", V1])`) methods. Almost all AKI-positive patients had a nadir inpatient eGFR below 90 mL/min/1.73 m^2^ using the Schwartz 120 (`r sprintf("%.1f%%", 100*(1-dt.frac.over.gfr[bl.disp == "Schwartz 120", V1]))`) and FAS (`r sprintf("%.1f%%", 100*(1-dt.frac.over.gfr[bl.disp == "FAS", V1]))`) methods.

<br> 

## Comparison of "empirical-first" and universal imputation baseline determination methods

Among the "empirical-first" methods, incidence was lowest using the mean-FAS (`r sprintf("%.1f%%", 100*dt.inc.empfirst[baseline == "mean.FAS", inc])`) and proximal-FAS (`r sprintf("%.1f%%", 100*dt.inc.empfirst[baseline == "prox.FAS", inc])`) baseline methods and highest using the lowest-Schwartz method (`r sprintf("%.1f%%", 100*dt.inc.empfirst[baseline == "bl.min", inc])`; Table 1). Regardless of baseline determination method, AKI incidence was similar for male and female patients and highest for patients under 1 year of age (Table 1). 

<br> 

Sensitivity of AKI as a predictor was positively correlated with incidence (baseline methods that produced a higher AKI incidence led to higher sensitivity), and specificity was negatively correlated with incidence (baseline methods that produced a higher AKI incidence led to lower specificity) across all three outcomes (Figure S10; Figure S11). However, the Schwartz 120 method yielded higher sensitivity for ICD-coded diagnosis and post-hospitalization CKD than the lowest-Schwartz and second lowest-Schwartz methods despite yielding a lower incidence, and the FAS method yielded a higher sensitivity for ICD-coded diagnosis and post-hospitalization CKD than the second lowest-FAS method despite yielding a lower incidence. 


#####
               
# Discussion

Although the KDIGO consensus definition has improved our epidemiologic understanding of AKI, it does not currently specify an approach to baseline creatinine determination. Our analysis confirmed the findings of prior studies that the approach used to determine baseline greatly impacts AKI incidence in both adults and children [@Siew2015; @Siew2010; @Zappitelli2008]. Using the lowest available pre-admission creatinine, which is the current standard practice, resulted in the highest incidence, likely identifying patients with milder disease.  When the second lowest creatinine, the average creatinine, and the most proximal creatinine were used as the reference baseline, AKI incidence was lower in all three cases.  This suggests that pre-admission creatinine values are dynamic, even over a three-month period.  Furthermore, it raises concern that a single, spuriously low creatinine value may falsely identify a patient as having AKI.  Using an average value or the second lowest value might mitigate this issue, however, such an approach may be challenging to operationalize.  The two imputation methods used most ubiquitously in practice, Schwartz 120 and FAS, both resulted in a lower AKI incidence than using the lowest available creatinine; of the two, FAS led to the  lower incidence, likely due to use of a lower reference GFR (i.e. 107.3 vs. 120mL/min/1.73 m^2^).

<br>

More importantly, our findings demonstrate that the approach used to determine a patient’s baseline creatinine dramatically affects the association between AKI and clinical outcomes. Among patients with available pre-admission data, when compared with using the lowest available creatinine as the baseline, universally imputing a baseline creatinine led to a stronger association between AKI and the application of an ICD 9/10 AKI billing code, which we used as a proxy for more severe, clinically overt disease.  Additionally, both universal imputation methods (Schwartz 120 and FAS) resulted in a stronger association with the development of post-hospitalization CKD.  In practice, however, many hospitalized children do not have a pre-admission creatinine; indeed, approximately 40% of our population had no such value available. Currently, the standard approach is to use the baseline value if it is available and to extrapolate a baseline if one is not.  Notably, the aforementioned findings were true even when we analyzed our entire cohort.  For both receipt of an ICD 9/10 AKI code and the development of post-hospitalization CKD, the association was stronger with universally imputed baselines than when the lowest pre-admission creatinine was used if available (lowest-Schwartz and lowest-FAS methods; Figure 5).  With regard to inpatient mortality, there was no statistically significant difference across the empirical/imputation methods; mortality is a rare event and it is possible that our analysis was not sufficiently large to detect an effect. In adult populations where mortality is more common, baseline determination method had an impact on post-AKI mortality rates [@Siew2010].

<br>

These findings, when interpreted in combination with the higher AKI incidence identified, suggest that using the lowest pre-admission creatinine as the baseline identifies a subset of milder AKI that may be of lesser clinical relevance. Indeed, when nadir inpatient eGFR was analyzed, half of the patients considered AKI-positive using the lowest pre-admission SCr never experienced an abnormally reduced inpatient eGFR (i.e. < 90mL/min/1.73 m^2^). By comparison, imputed baseline methods classified more than 99% of patients with normal nadir eGFRs as AKI-negative. Of note, the sensitivity and specificity of AKI as a predictor of clinical outcomes using alternative baseline methods was highly correlated with incidence, necessitating the use of a measure that is incidence independent to achieve an unbiased comparison of predictive performance (i.e. AUROC); this is an important methodologic consideration for future studies. Stage 1, 2, and 3 AKI correspond to different cutoff points on the ROC curve, with a higher AKI stage yielding a higher specificity but lower sensitivity for clinical outcomes (Figure S14).

<br>

Imputed baseline methods differ fundamentally from empirical methods. Imputed baselines are calculated from few variables (height, age, sex) and therefore may not accurately reflect creatinine variability, which can be due to differences in body composition, comorbidities, and disease severity. Imputed baseline methods, however, can be used independent of pre-admission data.  In fact, when imputation is used, the relative increase criterion becomes mathematically equivalent to identifying all patients who experience a sufficiently low eGFR during the first week of hospitalization. For the Schwartz 120 method this threshold eGFR is the same for all patients, 80 mL/min/1.74 m^2^ (reference baseline eGFR of 120 mL/min/1.74 m^2^); for the height-independent FAS method, the threshold inpatient GFR below which patients meet the relative increase criteria for AKI depends on the patient’s age and sex. Notably, patients who are AKI-positive using imputed but not empirical methods, typically have higher pre-admission creatinines and may have a higher baseline risk of adverse outcomes (i.e. lower baseline GFR).  Alternatively, patients identified as AKI-positive using empirical but not imputed baseline methods typically exhibit abnormally low pre-admission creatinines, which may reflect malnutrition, volume overload, or measurement error.


<br>
         
While these findings are novel and intriguing, it is important to interpret our findings in the context of the study’s limitations.  First, this analysis is retrospective in nature.  As such, it is subject to the biases inherent to this methodology. Second, we used only the creatinine criteria from the KDIGO definition [@KDIGO2012].  Several studies demonstrate that use of the creatinine criteria in isolation underestimates AKI incidence [@Kellum2015; @Kaddourah2019]. Of note, those studies suggest that patients identified by the urine output criteria do have outcomes similar to those identified by the creatinine criteria; therefore, while inclusion of urine output criteria may have impacted AKI incidence, it is unlikely to have modified the association identified between AKI and outcomes. That said, it is certainly possible that our findings are not reflective of the entire AKI population. Third, we utilized ICD codes as a proxy for severe, clinically relevant AKI events.  Historically, ICD codes have been maligned for their lack of sensitivity, however, the best available data suggest that the codes have relatively high specificity and we believe they serve as an accurate proxy for disease severity [@Waikar2006; @Vlasschaert2011; @Siew2016]. Fourth, we assessed these baseline methods in patients as young as 3 months of age even though the imputation methods were developed in patients one year or older [@Hoste2014; @Schwartz2009]. This was done to reflect current clinical practice, and a sensitivity analysis where patients between 3 months and one year of age were excluded found similar results (Figure S9). Baseline GFR has been shown to be lower in patients younger than two years of age [@Hessey2017], and further research should explore baseline imputation for this age group. Fifth, in diagnosing post-hospitalization CKD we did not include any restrictions on the interval between qualifying post-stay labs, where eGFR was below 60 mL/min/1.73 m^2^, to avoid further restricting the number of patients with the labs necessary to assess post-hospitalization CKD. The formal diagnosis criteria for CKD require that qualifying labs be separated by 90 days or more [@KDIGO2012], a condition met by 60% of the patients we classified as developing post-hospitalization CKD (Figure S10). Finally, it is important to note that this analysis was done in a pediatric population where the incidence of pre-hospitalization CKD is low.  In other populations where the prevalence of CKD is higher (i.e., adult populations), it is likely that imputation of baseline serum creatinine will be less accurate at identifying true AKI events and may correlate less strongly with outcomes [@Sutherland2018; @Sutherland2017].


<br>

In summary, our analysis confirms prior data which demonstrate that the choice of baseline determination method has a significant impact on AKI incidence.  It expands on these findings to demonstrate that baseline determination method affects the relationship between AKI and subsequent outcomes as well.  The fact that preadmission creatinine values are commonly absent is a well-known challenge across AKI epidemiologic studies, particularly in children [@Siew2010; @Kaddourah2017].   Interestingly, our findings suggest that consideration of preadmission creatinine values does not strengthen the association between AKI and clinical outcomes. Indeed, outright disregard for preadmission creatinine when defining AKI led to a stronger association with clinical outcomes in our cohort. Notably, the strong association with post-hospitalization CKD suggests that using an imputed baseline could improve post-AKI CKD risk estimation, an important clinical challenge [@Nelson2019; @Hessey2019; @Hsu2016; @Sigurjonsdottir2018]. At a minimum, our findings underscore the importance of a uniform approach to baseline determination in AKI studies.  We believe that every effort should be made to standardize the concept and, until then, that AKI epidemiologic studies should be interpreted in the context of their baseline determination method.



<!---BLOCK_LANDSCAPE_START--->
# Tables


```{r table1, echo = FALSE, eval = TRUE}
t1
```

**Table 1** Incidence of AKI according to baseline determination method by patient demographic group. One patient whose sex was recorded as ‘unknown’ was excluded from the sex-based analysis. Patients whose primary race was recorded as ‘declines to state’ (52 patients) or ‘unknown’ (16 patients) were grouped with ‘other’. Abbreviations: AKI, acute kidney injury; FAS, full age spectrum.

<!---BLOCK_LANDSCAPE_STOP--->
#####

# Figures


![](inclusion_flow_Scott.png)

**Figure 1** Inclusion and exclusion criteria.  93,105 pediatric hospitalizations were initially evaluated. 72,839 were excluded leaving `r format(n_pt, big.mark=',')` admissions for the primary analysis. For the post-hospitalization CKD analysis, `r format(n_pt - n_pt_CKD, big.mark=',')` were excluded as they did not have at least two post-discharge creatinines leaving `r format(n_pt_CKD, big.mark=',')` hospitalizations for the final CKD analysis. Abbreviations: CKD, chronic kidney disease; SCr, serum creatinine



#####



![](baseline_mean_sd_preSCr_fig.png)
**Figure 2** Mean and standard deviation of baseline serum creatinine by determination method among patients with measured pre-admission creatinine.



#####

![](ROC_and_AUC_preSCr_fig.png)

**Figure 3** ROC curves (top) and AUROC (bottom) for AKI as a predictor of clinical outcomes based on baseline method for patients with measured pre-admission creatinine.

#####


![](Nadir_gfr_fig.png)

**Figure 4** Distribution of nadir estimated GFR for AKI+ patients under alternative baseline determination methods for patients with measured preadmission SCr.


#####


![](ROC_and_AUC_emp_first_fig.png)

**Figure 5** ROC curves (top) and AUROC (bottom) for AKI as a predictor of clinical outcomes based on baseline method for empirical-first and universal baseline determination methods.



#####

# References

<div id="refs"></div>
