PRAPaperAnalysis.Rmd

---
title: "Panel Reactive Antibodies (PRAs) and Their Association with Blood Product Exposures in Pediatric Patients Undergoing Cardiac Transplantation"
author:
  - 'Dr. Nina Guzzetta, Dr. Susana Cruz Beltran'
  - 'Statisticians: Dr. Renee Moore, Andrea Lane'
date: "3/23/2021"
geometry: margin=1cm
output:
  pdf_document:
    toc: false
    number_sections: false
header-includes:
    - \usepackage{caption}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE)
knitr::opts_chunk$set(fig.width=6, fig.height=4) 
library(readxl)
library(epitrix)
library(dplyr)
library(rockchalk)
library(kableExtra)
library(summarytools)
library(ggplot2)
library(ggmosaic)
library(tidyverse)
library(berryFunctions)
library(caret)
library(glmnet)

options(knitr.table.format = "latex")

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00")

```

\captionsetup[table]{labelformat=empty}

```{r upload}
df.sheet1 <- read_excel("Transplant data.final.withUNOS.xlsx", sheet = 1,
                        col_types = c("text","text","date","text","numeric",
                                      "text","text","text", "date","date",
                                      "text","text","text","numeric","text",
                                    "text",rep("numeric",5),"text",rep("numeric",16),"text"))

df.sheet2 <- read_excel("Transplant data.final.withUNOS.xlsx", sheet = 2,
                        col_types = c("text","text","date","text","numeric",
                                      "text","text","text", "date","date",
                                "text","text","text","numeric",rep("numeric",6),"text",
                                    "text",rep("numeric",5),"text",rep("numeric",16),"text"))

df.neg <- df.sheet1[-1,]
df.pos <- df.sheet2[-1,]

names(df.neg) <- clean_labels(names(df.neg))
names(df.pos) <- clean_labels(names(df.pos))

names(df.neg)[1] <- "Subj_num"
names(df.pos)[1] <- "Subj_num"

names(df.neg)[25] <- "prior_rbc"
names(df.neg)[26] <- "prior_ffp"
names(df.neg)[27] <- "prior_plt"
names(df.neg)[28] <- "prior_cryo"

names(df.neg)[30] <- "dos_rbc"
names(df.neg)[31] <- "dos_ffp"
names(df.neg)[32] <- "dos_plt"
names(df.neg)[33] <- "dos_cryo"

names(df.neg)[35] <- "post_rbc"
names(df.neg)[36] <- "post_ffp"
names(df.neg)[37] <- "post_plt"
names(df.neg)[38] <- "post_cryo"

names(df.neg)[39] <- "unos"


names(df.pos)[31] <- "prior_rbc"
names(df.pos)[32] <- "prior_ffp"
names(df.pos)[33] <- "prior_plt"
names(df.pos)[34] <- "prior_cryo"

names(df.pos)[36] <- "dos_rbc"
names(df.pos)[37] <- "dos_ffp"
names(df.pos)[38] <- "dos_plt"
names(df.pos)[39] <- "dos_cryo"

names(df.pos)[41] <- "post_rbc"
names(df.pos)[42] <- "post_ffp"
names(df.pos)[43] <- "post_plt"
names(df.pos)[44] <- "post_cryo"

names(df.pos)[45] <- "unos"

df.all <- bind_rows("1" = df.neg, "2" = df.pos, .id="pra_status")

```

```{r clean}
df.all$age_unit <- substr(df.all$age_years,nchar(df.all$age_years),nchar(df.all$age_years))
df.all$age_num <- substr(df.all$age_years,1,nchar(df.all$age_years)-1)

df.all$age <- ifelse(df.all$age_unit == "m", as.numeric(df.all$age_num)/12, ifelse(df.all$age_unit == "w", as.numeric(df.all$age_num)/52, as.numeric(df.all$age_num)))

df.all$abo_incompatible_nomiss <- ifelse(is.na(df.all$abo_incompatible),0,1)
df.all$prior_mech_nomiss <- ifelse(is.na(df.all$prior_mech_support),0,1)

df.all$plasma_c <- ifelse(df.all$plasmapheresis %in% c("Yes","yes"),1,0)
df.all$rej_c <- ifelse(df.all$rejection=="N",0,1)

df.all$plasma_c <- factor(df.all$plasma_c, levels = c(0,1),
                          labels = c("No","Yes"))
df.all$rej_c <- factor(df.all$rej_c, levels = c(0,1), labels = c("No","Yes"))
df.all$ecmo_c <- factor(df.all$ecmo_support_postop, levels = c("no","yes"),
                        labels = c("no","yes"))
df.all$diag_c <- df.all$diagnosis
df.all$diag_c[df.all$diag_c=="LV failure"] <- "Dilated CM"
df.all$diag_c[df.all$diag_c=="RV Failure" | df.all$diag_c=="Shone's complex- Coarct - Multiple VSD"] <- "Other"

df.all$diag_c <- factor(df.all$diag_c, 
      levels = c("ARVD","Dilated CM","Failed SV","Restrictive CM","Other"),
      labels = c("ARVD","Dilated CM","Failed SV","Restrictive CM","Other"))
df.all$donor_c <- factor(df.all$donor_specific_ab, levels = c("N","Y"),
                         labels = c("No","Yes"))
df.all$gender_c <- factor(df.all$gender, levels = c("F","M"), labels = c("Female","Male"))
df.all$race_c <- factor(df.all$race, 
        levels = c("White","African American","Hispanic","Asian","Declined"),             labels = c("White","African American","Hispanic","Asian","No response"))
df.all$abo_incompatible_nomiss <- factor(df.all$abo_incompatible_nomiss, 
                                         levels = c(0,1), labels = c("No","Yes"))
df.all$prior_mech_nomiss <- factor(df.all$prior_mech_nomiss, levels = c(0,1),
                                         labels = c("No","Yes"))

df.all$unos_new <- ifelse(df.all$unos=="1A",0,1)
df.all$unos_f <- factor(df.all$unos_new,levels = c(0,1),labels=c("1A","1B"))

df.all <- df.all[!is.na(df.all$prior_donor_exposures),]
```

```{r functions}
my2x2 <- function(X,Y,col=TRUE){
  perc <- ifelse(col,2,1)
  mytab <- table(X,Y)
  proptab <- format(round(prop.table(mytab,perc)*100,1),nsmall = 1,trim=TRUE)
  mytab2 <- matrix(paste0(mytab," (",proptab,")"),nlevels(X),2)
  mytab2 <- rbind(mytab2, colSums(mytab))
  mytab <- cbind(mytab2, c(rowSums(mytab),sum(rowSums(mytab))))
  return(mytab)
}

mypval <- function(X,Y,symbol=TRUE){
  if (class(X)=="numeric"){
    test <- "Wilcoxon rank-sum"
    pval <- wilcox.test(X[df.all$pra_status==1],X[df.all$pra_status==2])$p.value
    pval <- ifelse(pval == 1,0.99,pval)
    pval <- format(round(pval,2), nsmall=2)
    pval <- ifelse(pval==0,"<0.0001",pval)
  } else if (class(X)=="factor"){
  
  tab <- table(X,Y)
  exp.tab <- chisq.test(tab)$expected
  test <- ifelse(all(exp.tab > 5),"Chi Square","Fisher's Exact")
  pval <- ifelse(all(exp.tab > 5),chisq.test(tab,correct = FALSE)$p.value,fisher.test(tab)$p.value)
  pval <- ifelse(pval == 1,0.99,pval)
  pval <- format(round(pval,2), nsmall=2)
  pval <- ifelse(pval=="0.00","<0.01",pval)
  if(symbol){pval <- ifelse(test=="Chi Square",
                 paste0(pval,""),
                 paste0(pval,"^"))}
   }
  return(list(pval = pval, test = test))
}

varlevels <- function(X){
  if(class(X)=="numeric"){
    nlev <- 1
  } else if (class(X)=="factor"){
    nlev <- nlevels(X)
  }
  return(nlev)
}

mystats <- function(X,Y){
  if (class(X)=="numeric"){
    tab <- matrix(NA,1,3)
    all_1 <- format(round(median(X,na.rm = TRUE),1),nsmall = 1)
    all_2 <- paste0("(",format(round(quantile(X,na.rm = TRUE)[2],1),nsmall=1),", ",
                    format(round(quantile(X,na.rm = TRUE)[4],1),nsmall=1),")")
    
    strat1_1 <- format(round(median(X[Y==1],na.rm = TRUE),1),nsmall=1)
    strat2_1 <- format(round(median(X[Y==2],na.rm = TRUE),1),nsmall=1)
    
    strat1_2 <- paste0("(",format(round(quantile(X[Y==1],na.rm = TRUE)[2],1),nsmall=1),", ",
                       format(round(quantile(X[Y==1],na.rm = TRUE)[4],1),nsmall=1),")")
    strat2_2 <- paste0("(",format(round(quantile(X[Y==2],na.rm = TRUE)[2],1),nsmall = 1),", ",
                       format(round(quantile(X[Y==2],na.rm = TRUE)[4],1),nsmall = 1),")")
    
    cell1 <- paste(all_1,all_2)
    cell2 <- paste(strat1_1,strat1_2)
    cell3 <- paste(strat2_1,strat2_2)
    
    tab <- cbind(cell1,cell2,cell3)
  } else if (class(X)=="factor"){
    nlev <- nlevels(X)
    all_1 <- c()
    all_2 <- c()
    strat1_1 <- c()
    strat2_1 <- c()
    strat1_2 <- c()
    strat2_2 <- c()
    tab <- matrix(NA,nlev,3)
    for(i in 1:nlev){
    tab1 <- table(X)
    all_1[i] <- paste(tab1[i])
    all_2[i] <- paste0("(",format(round(prop.table(tab1)[i]*100,1),nsmall = 1),")")
    
    tab2 <- table(X,Y)
    strat1_1[i] <- tab2[i,1]
    strat1_2[i] <- paste0("(",format(round(prop.table(tab2,2)[i,1]*100,1),nsmall = 1),")")
    
    strat2_1[i] <- tab2[i,2]
    strat2_2[i] <- paste0("(",format(round(prop.table(tab2,2)[i,2]*100,1),nsmall = 1),")")
      
    }
    cell1 <- paste(all_1,all_2)
    cell2 <- paste(strat1_1,strat1_2)
    cell3 <- paste(strat2_1,strat2_2)

    tab <- cbind(cell1,cell2,cell3)

  }
  return(tab)
}

mystats2 <- function(X,Y){
  if(class(X)=="numeric"&&class(Y)=="factor"){
    tab <- matrix(NA,1,3)
    all_1 <- format(round(median(X,na.rm = TRUE),1),nsmall = 1)
    all_2 <- paste0("(",format(round(quantile(X,na.rm = TRUE)[2],1),nsmall=1),", ",
                    format(round(quantile(X,na.rm = TRUE)[4],1),nsmall=1),")")
    
    strat1_1 <- format(round(median(X[Y==levels(Y)[1]],na.rm = TRUE),1),nsmall=1)
    strat2_1 <- format(round(median(X[Y==levels(Y)[2]],na.rm = TRUE),1),nsmall=1)
    
    strat1_2 <- paste0("(",format(round(quantile(X[Y==levels(Y)[1]],na.rm = TRUE)[2],1),nsmall=1),", ",
                       format(round(quantile(X[Y==levels(Y)[1]],na.rm = TRUE)[4],1),nsmall=1),")")
    strat2_2 <- paste0("(",format(round(quantile(X[Y==levels(Y)[2]],na.rm = TRUE)[2],1),nsmall = 1),", ",
                       format(round(quantile(X[Y==levels(Y)[2]],na.rm = TRUE)[4],1),nsmall = 1),")")
    
    cell1 <- paste(all_1,all_2)
    cell2 <- paste(strat1_1,strat1_2)
    cell3 <- paste(strat2_1,strat2_2)

    tab <- cbind(cell1,cell2,cell3)
  }
  
  if(class(X)=="factor"&&class(Y)=="factor"){
    nlev <- nlevels(X)
    all_1 <- c()
    all_2 <- c()
    strat1_1 <- c()
    strat2_1 <- c()
    strat1_2 <- c()
    strat2_2 <- c()
    tab <- matrix(NA,nlev,3)
    for(i in 1:nlev){
    tab1 <- table(X)
    all_1[i] <- paste(tab1[i])
    all_2[i] <- paste0("(",format(round(prop.table(tab1)[i]*100,1),nsmall = 1),")")
    
    tab2 <- table(X,Y)
    strat1_1[i] <- tab2[i,1]
    strat1_2[i] <- paste0("(",format(round(prop.table(tab2,2)[i,1]*100,1),nsmall = 1),")")
    
    strat2_1[i] <- tab2[i,2]
    strat2_2[i] <- paste0("(",format(round(prop.table(tab2,2)[i,2]*100,1),nsmall = 1),")")
      
    }
    cell1 <- paste(all_1,all_2)
    cell2 <- paste(strat1_1,strat1_2)
    cell3 <- paste(strat2_1,strat2_2)

    tab <- cbind(cell1,cell2,cell3)

  }
  if(class(X)=="numeric"&&class(Y)=="numeric"){
    
    tab <- cor.test(df.blood$prior_donor_exposures,df.blood$cpra)$estimate

  }

  return(tab)
}

otherstats <- function(X,Y){
  if(class(X)=="factor"&&class(Y)=="numeric"){
  stattab <- t(mystats2(Y,X))
  pval <- format(round(wilcox.test(Y[X==levels(X)[1]],Y[X==levels(X)[2]])$p.value,2),nsmall=2)
  pval <- ifelse(pval == 1,0.99,pval)
  pval <- ifelse(pval=="0.00","<0.01",pval)
  pval <- rep(pval,nlevels(X))
  tab <- cbind(stattab,pval)
  }
  if(class(X)=="numeric"&&class(Y)=="numeric"){
      stattab <- format(round(cor.test(X,Y,method = "spearman")$estimate,2),nsmall=2)
      pval <- format(round(cor.test(X,Y,method = "spearman")$p.value,2),nsmall = 2)
      tab <- cbind(stattab,pval)
  }
  
  return(tab)
}

```

\clearpage
\section{Tables}
```{r table1}
varlist <- c("age","weight","no_of_prior_surgeries","gender_c","race_c","prior_mech_nomiss","diag_c","unos_f","donor_time","cpb_time","abo_incompatible_nomiss","plasma_c","icu_length_of_stay_day", "vis", "time_to_extubation_hr","rej_c", "donor_c","ecmo_c")
# var_levels <- lapply(df.all[varlist], FUN = varlevels)
# df_levels <- data.frame(matrix(unlist(var_levels), nrow=length(varlist),byrow = T))

myfulltab <- lapply(df.all[varlist], FUN = mystats, Y = df.all$pra_status)

mat.tab <- as.matrix(unlist(myfulltab[[1]]))
for (i in 2:length(myfulltab)){
  mat.tab <- rbind(mat.tab,as.matrix(unlist(myfulltab[[i]])))
}

mypvals <- lapply(df.all[varlist], function(x)
                  rep(mypval(x,df.all$pra_status)$pval,varlevels(x)))
# mytest <- lapply(df.all[varlist], function(x)
#                   rep(mypval(x,df.all$pra_status)$test,varlevels(x)))
# mytestmat <- as.matrix(unlist(mytest))
mypvalmat <- as.matrix(unlist(mypvals))
mat.tab.p <- cbind(mat.tab,mypvalmat)
all.tab <- data.frame(mat.tab.p)

names(all.tab)[2] <- paste0(names(all.tab)[2],
  footnote_marker_symbol(1, "latex"))
rownames(all.tab) <- c("Age (Years)", "Weight at transplantation (kg)", 
                      "Number of previous surgeries","Female", "Male", "White",
                      "African American", "Hispanic","Asian", "No response",  
                       "No", "Yes","ARVD","Dilated CM","Failed SV",
                      "Restrictive CM","Other","1A","1B","Donor cross clamp time (min)",
                      "CPB time (min)","No ","Yes ","No  ","Yes  ",
                      "Length of ICU stay (days)","Vasoactive inotrope score",
                      "Time to extubation (hrs)","No   ","Yes   ","No    ",
                      "Yes    ","No     ","Yes     ")

demonly <- matrix(all.tab[1:25,1],25,1)
rownames(demonly) <- rownames(all.tab)[1:25]
kable(demonly, "latex", booktabs = T, align = "l", col.names=c("All patients"),caption = "Table 1: Demographic, Preoperative, and Intraoperative Variables for all patients (N=64)") %>%
  kable_styling(latex_options = "hold_position") %>% 
  group_rows("Demographic and preoperative",1,19,bold=TRUE) %>% 
  group_rows("Sex", 4, 5, bold = FALSE) %>% 
  group_rows("Race", 6, 10, bold = FALSE) %>% 
  group_rows("Mechanical support prior to transplantation", 11, 12, bold = FALSE) %>% 
  group_rows("Intraoperative",20,25,bold=TRUE) %>% 
  group_rows("Diagnosis",13,17,bold=FALSE) %>% 
  group_rows("UNOS status",18,19,bold=FALSE) %>% 
  group_rows("ABO incompatible",22,23,bold=FALSE) %>% 
  group_rows("Plasmapheresis",24,25,bold=FALSE) %>% 
  footnote(general = "Continuous variables expressed as median (Q1,Q3); Categorical variables expressed as n(%)",threeparttable=T)



kable(all.tab[1:25,], "latex", booktabs = T, align = "l", col.names = c("All patients (n=64)", "HLA ab- (n=35)", "HLA ab+ (n=29)", "p value"), caption = "Table 2a: Demographic, Preoperative, and Intraoperative Variables by binary HLA ab status") %>%
  kable_styling(latex_options = "hold_position") %>% 
  group_rows("Demographic and preoperative",1,19,bold=TRUE) %>% 
  group_rows("Sex", 4, 5, bold = FALSE) %>% 
  group_rows("Race", 6, 10, bold = FALSE) %>% 
  group_rows("Mechanical support prior to transplantation", 11, 12, bold = FALSE) %>% 
  group_rows("Intraoperative",20,25,bold=TRUE) %>% 
  group_rows("Diagnosis",13,17,bold=FALSE) %>% 
  group_rows("UNOS status",18,19,bold=FALSE) %>% 
  group_rows("ABO incompatible",22,23,bold=FALSE) %>% 
  group_rows("Plasmapheresis",24,25,bold=FALSE) %>% 
  collapse_rows(columns = 5,valign = "top", latex_hline = "none") %>% 
  footnote(general = "Continuous variables expressed as median (Q1,Q3), Wilcoxon rank sum test; Categorical variables expressed as n(%), Chi square test or Fisher's exact test^", threeparttable=T)
  


```


\clearpage

```{r bloodprodclean}
##create binary variables for did/did not have prior exposures
df.blood <- df.all[!is.na(df.all$prior_donor_exposures),]
df.blood$prior_bin <- factor(ifelse(df.blood$prior_donor_exposures>0,1,0),levels = c(0,1))
df.blood$rbc_bin <- factor(ifelse(df.blood$prior_rbc>0,1,0),levels = c(0,1))
df.blood$ffp_bin <- factor(ifelse(df.blood$prior_ffp>0,1,0),levels = c(0,1))
df.blood$plt_bin <- factor(ifelse(df.blood$prior_plt>0,1,0),levels = c(0,1))
df.blood$cryo_bin <- factor(ifelse(df.blood$prior_cryo>0,1,0),levels = c(0,1))
df.blood$cell <- df.blood$prior_rbc + df.blood$prior_plt
df.blood$acell <- df.blood$prior_ffp + df.blood$prior_cryo
df.blood$cell_bin <- factor(ifelse(df.blood$cell>0,1,0),levels = c(0,1))
df.blood$acell_bin <- factor(ifelse(df.blood$acell>0,1,0),levels = c(0,1))
df.blood$pra_bin <- factor(ifelse(df.blood$pra_status==1,0,1),levels=c(0,1))

```

\clearpage
```{r allcpratab}
varlist <- c("age","weight","no_of_prior_surgeries","gender_c","prior_mech_nomiss","donor_time","cpb_time","abo_incompatible_nomiss","plasma_c")

myfulltab2 <- lapply(df.blood[varlist], FUN = otherstats, Y = df.blood$cpra)

mat.tab2 <- as.matrix(unlist(myfulltab2[[1]]))
for (i in 2:length(myfulltab2)){
  mat.tab2 <- rbind(mat.tab2,as.matrix(unlist(myfulltab2[[i]])))
}

mat.tab2 <- mat.tab2[-c(4,7,12,15),]
rownames(mat.tab2) <- c("Age (Years)", "Weight at transplantation (kg)",
                      "Number of previous surgeries","Female", "Male",
                       "No", "Yes","Donor cross clamp time (min)",
                      "CPB time (min)","No ","Yes ","No  ","Yes  ")


kable(mat.tab2, "latex", booktabs = T, align = "l", col.names = c("HLA ab","p value"), caption = "Table 2b: Pre, Intra, and Postoperative variables by continuous HLA ab") %>%
  kable_styling(latex_options = "hold_position") %>% 
  group_rows("Demographic and preoperative",1,7,bold=TRUE) %>% 
  group_rows("Sex",4,5,bold=FALSE) %>% 
  group_rows("Mechanical support prior to transplantation",6,7,bold=FALSE) %>% 
  group_rows("Intraoperative",8,13,bold=TRUE) %>% 
  group_rows("ABO incompatible",10,11,bold=FALSE) %>% 
  group_rows("Plasmapheresis",12,13,bold=FALSE) %>% 
  collapse_rows(columns = 3,valign = "top", latex_hline = "none") %>% 
  footnote(general = "Continuous variables expressed as Spearman rank correlation $r$; Categorical variables expressed as median (Q1,Q3), Wilcoxon rank sum test",threeparttable=T,escape=F)

```


\clearpage

```{r studyobjectives}

varlist <- c("prior_donor_exposures","prior_rbc","prior_plt","prior_ffp","prior_cryo","prior_bin","rbc_bin","plt_bin","ffp_bin","cryo_bin")

myfulltab <- lapply(df.blood[varlist], FUN = mystats2, Y = df.blood$pra_bin)

mat.tab <- as.matrix(unlist(myfulltab[[1]]))
for (i in 2:length(myfulltab)){
  mat.tab <- rbind(mat.tab,as.matrix(unlist(myfulltab[[i]])))
}

mypvals <- lapply(df.blood[varlist], function(x)
                  rep(mypval(x,df.blood$pra_bin)$pval,varlevels(x)))
mypvalmat <- as.matrix(unlist(mypvals))
mat.tab.p <- cbind(mat.tab,mypvalmat)
df.tab <- data.frame(mat.tab.p)



varlist2 <- c("prior_donor_exposures","prior_rbc","prior_plt","prior_ffp","prior_cryo","prior_bin","rbc_bin","plt_bin","ffp_bin","cryo_bin")

myfulltab2 <- lapply(df.blood[varlist], FUN = otherstats, Y = df.blood$cpra)

mat.tab2 <- as.matrix(unlist(myfulltab2[[1]]))
for (i in 2:length(myfulltab2)){
  mat.tab2 <- rbind(mat.tab2,as.matrix(unlist(myfulltab2[[i]])))
}

mat.tab2 <- mat.tab2[-c(6,9,12,15,18),]

rownames(df.tab) <- rownames(mat.tab2) <- c("Total prior donor exposures ","Red blood cells ","Platelets ","Fresh frozen plasma ","Cryoprecipitate ","No","Yes","No ","Yes ","No  ","Yes  ","No   ","Yes   ","No    ","Yes    ")

kable(df.tab, "latex", booktabs = T, align = "l", col.names = c("All patients (n=64)", "HLA ab- (n=35)", "HLA ab+ (n=29)", "p value"), caption = "Table 3a: Blood products and binary HLA ab status") %>%
  kable_styling(latex_options = "hold_position") %>% 
   group_rows("Blood products - binary n (%)",6,15,bold=TRUE) %>% 
   group_rows("Blood products - continuous median(Q1,Q3)", 1, 5, bold = TRUE) %>% 
   collapse_rows(columns = 5,valign = "top", latex_hline = "none") %>% 
   row_spec(5, hline_after = TRUE) %>%
   group_rows("Total prior donor exposures",6,7,bold=FALSE) %>% 
  group_rows("Red blood cells",8,9,bold=FALSE) %>%
  group_rows("Platelets",10,11,bold=FALSE) %>%
  group_rows("Fresh frozen plasma",12,13,bold=FALSE) %>%
  group_rows("Cryoprecipitate",14,15,bold=FALSE) %>%
  footnote(general = "Continuous blood products: Wilcoxon rank sum test; Binary blood products: Chi square test or Fisher's exact test^",threeparttable=T)

```
\clearpage

```{r bloodprod2}

kable(mat.tab2, "latex", booktabs = T, align = "l", col.names = c("HLA ab", "p value"), caption = "Table 3b: Blood products and continuous HLA ab") %>%
  kable_styling(latex_options = "hold_position") %>% 
    group_rows("Blood products - binary: median (Q1,Q3)",6,15,bold=TRUE) %>% 
    group_rows("Blood products - continuous: correlation r", 1, 5, bold = TRUE) %>% 
    collapse_rows(columns = 3,valign = "top", latex_hline = "none") %>% 
    row_spec(5, hline_after = TRUE) %>%
    group_rows("Total prior donor exposures",6,7,bold=FALSE) %>% 
    #group_rows("RBCs",8,9,bold=FALSE) %>% 
   group_rows("Platelets",10,11,bold=FALSE) %>%
   group_rows("Fresh frozen plasma",12,13,bold=FALSE) %>%
   group_rows("Cryoprecipitate",14,15,bold=FALSE) %>%
  footnote(general = "Continuous blood product variables expressed as Spearman rank correlation $r$; Binary blood product variables expressed as median (Q1,Q3), Wilcoxon rank sum test",threeparttable=T,escape=F)

```

\clearpage
```{r postop}
kable(all.tab[26:34,], "latex", booktabs = T, align = "l", col.names = c("All patients (n=64)", "HLA ab- (n=35)", "HLA ab+ (n=29)", "p value"), caption = "Table 4: Postoperative Variables by binary HLA ab status") %>%
  kable_styling(latex_options = "hold_position") %>% 
  group_rows("Rejection",4,5,bold=FALSE) %>% 
  group_rows("Donor specific antibody",6,7,bold=FALSE) %>% 
  group_rows("ECMO",8,9,bold = FALSE) %>% 
  collapse_rows(columns = 5,valign = "top", latex_hline = "none") %>% 
  footnote(general = "Continuous variables expressed as median (Q1,Q3), Wilcoxon rank sum test; Categorical variables expressed as n(%), Chi square test or Fisher's exact test^", threeparttable=T)

```




```{r bloodprodtabs}
make2x2table <- function(X,Y,mycols,myrows,caption){
  tab <- my2x2(X,Y)
  pval <- mypval(X,Y,FALSE)
  foot <- paste0("p = ",pval$pval," (",pval$test,")")
  rownames(tab) <- myrows
  rowspec <- nlevels(X)
  mytable <- kable(tab, col.names = mycols, row.names = TRUE, "latex", booktabs = T,
                   caption=caption) %>%
  kable_styling(latex_options = "hold_position", full_width = FALSE) %>%
  footnote(general = foot) %>%
  row_spec(rowspec, hline_after = TRUE) %>%
  column_spec(c(2,4), border_left  = TRUE)

  return(mytable)
}
colnames.all <- c("HLA ab-", "HLA ab+", "Total")

```

```{r alleles}
df.alleles <- df.blood[df.blood$pra_status==2,]

df.alleles$allele_c <- ifelse(df.alleles$class_i==1 & df.alleles$class_ii==1,3,
                              ifelse(df.alleles$class_i==1 & df.alleles$class_ii==0,1,2))

df.alleles$allele_c <- factor(df.alleles$allele_c, levels=c(1,2,3),
                              labels = c("1","2","both"))

```

```{r tab17}
make2x2table(df.alleles$allele_c,df.alleles$donor_c,
             mycols = c("No donor \nspecific antibody","Donor specific \nantibody","Total"),
             myrows=c("Class I","Class II","Class I and II","Total"),
             "Table 5: Allele class and donor specific antibody")
```

\newpage 

```{r scatter}
plot(df.blood$age, df.blood$cpra,main="Relationship between age and continuous HLA ab \n (r=-0.04)",xlab="Age",ylab="HLA ab")
abline(lm(df.blood$cpra~df.blood$age), col="red")

```

```{r scatter2}
plot(df.blood$weight, df.blood$cpra,main="Relationship between weight and continuous HLA ab \n(r=-.002)",xlab="Weight",ylab="HLA ab")
abline(lm(df.blood$cpra~df.blood$weight), col="red")
```

```{r scatter 3}
plot(df.blood$prior_donor_exposures,df.blood$cpra,main="Relationship between total blood product exposures \nand continuous HLA ab (r=0.27)",xlab="Total blood product exposures",ylab="HLA ab")
abline(lm(df.blood$cpra~df.blood$prior_donor_exposures), col="red")


```


```{r univar,eval=FALSE}
unimod <- function(var){
    mod <- glm(pra_bin ~ var,data=df.blood,family="binomial")
    coefs <- round(summary(mod)$coefficients[,1],2)
    coefs.p <- matrix(paste0(coefs," (",round(summary(mod)$coefficients[,2],2),")"),nrow=2,ncol=1)
    oddsrat <- matrix(paste0(round(exp(coefs),2)," (",round(exp(confint(mod))[,1],2),", ",round(exp(confint(mod))[,2],2),")"),nrow=2,ncol=1)
    pval <- matrix(rep(round(summary(mod)$coefficients[2,4],3),2),nrow=2,ncol=1)
    
    tab <- cbind(coefs.p,oddsrat,pval)

return(tab)
}


varlist <- c("prior_bin","rbc_bin","plt_bin","ffp_bin","cryo_bin","prior_donor_exposures","prior_rbc","prior_plt","prior_ffp","prior_cryo")

modlist <- lapply(df.blood[varlist], FUN=unimod)

mat.tab <- as.matrix(unlist(modlist[[1]]))
for (i in 2:length(modlist)){
  mat.tab <- rbind(mat.tab,as.matrix(unlist(modlist[[i]])))
}

rownames(mat.tab) <- c("Intercept","Total","Intercept ","Red blood cells","Intercept  ","Platelets","Intercept   ","Fresh frozen plasma","Intercept    ","Cryoprecipitate","Intercept     ","Total ","Intercept      ","Red blood cells ","Intercept       ","Platelets ","Intercept        ","Fresh frozen plasma ","Intercept         ","Cryoprecipitate ")

kable(mat.tab, "latex", booktabs = T, align = "l", col.names = c("Estimate (SE)","OR (95% CI)" ,"p value"), caption = "Table 4: Univariate results for blood product exposures") %>%
  kable_styling(latex_options = "hold_position") %>% 
  group_rows("Blood products - binary",1,10,bold=TRUE) %>% 
  group_rows("Blood products - count", 11, 20, bold = TRUE) %>% 
  row_spec(c(2,4,6,8,10,12,14,16,18), hline_after = TRUE) %>%
  collapse_rows(columns = 4,valign = "bottom", latex_hline = "none") %>% 
  footnote(general = "Absence of blood product is the refernce group for all binary variables")



cormat <- matrix(cor(df.blood[,c("prior_donor_exposures","prior_rbc","prior_plt","prior_ffp","prior_cryo")]),5,5)
rownames(cormat) <- c("Total","RBCs","Plt","FFP","Cryo")
kable(cormat, "latex", booktabs = T, align = "l", col.names = c("Total","RBCs","Plt","FFP","Cryo"), caption = "Table 5: Correlation matrix for blood products") %>%
  kable_styling(latex_options = "hold_position")


demvarlist <- c("age","weight","no_of_prior_surgeries","gender_c","prior_mech_nomiss")

modlist2 <- lapply(df.blood[demvarlist], FUN=unimod)

mat.tab2 <- as.matrix(unlist(modlist2[[1]]))
for (i in 2:length(modlist2)){
  mat.tab2 <- rbind(mat.tab2,as.matrix(unlist(modlist2[[i]])))
}

mat.tab2r <- insertRows(as.data.frame(mat.tab2),8,c("--","--",mat.tab2[7,3]))
mat.tab2r <- insertRows(as.data.frame(mat.tab2r),11,c("--","--",mat.tab2[10,3]))

rownames(mat.tab2r) <- c("Intercept","Age","Intercept ","Weight","Intercept  ","Number of previous surgeries","Intercept   ","Male (ref)","Female","Intercept    ","No (ref)","Yes")

kable(mat.tab2r, "latex", booktabs = T, align = "l", col.names = c("Estimate (SE)","OR (95% CI)" ,"p value"), caption = "Table 7: Univariate results for other demographic variables") %>%
  kable_styling(latex_options = "hold_position") %>%
  group_rows("Sex",7,9,bold=FALSE) %>%
  group_rows("Prior mechanical support", 10, 12, bold = FALSE) %>%
  row_spec(c(2,4,6,9), hline_after = TRUE) %>%
  collapse_rows(columns = 4,valign = "bottom", latex_hline = "none")


```


```{r model, eval=FALSE}
df.blood$race2 <- factor(ifelse(df.blood$race_c=="Hispanic"|df.blood$race_c=="Asian"|df.blood$race_c=="No response","H/A/NR",df.blood$race_c))

mod <- glm(pra_bin ~ prior_bin + age + weight + prior_mech_nomiss +  no_of_prior_surgeries + gender_c + race2 + diag_c, data=df.blood, family="binomial",na.action = na.pass)

#step.aic <- step(mod)

finalmod <- glm(pra_bin ~ prior_bin + prior_mech_nomiss + race2 + diag_c,data=df.blood,family="binomial")


coefs <- summary(finalmod)$coefficients[,1]
coefs.se <- summary(finalmod)$coefficients[,2]
oddsrat <- exp(coefs)
or.se <- exp(confint(finalmod))
pvals <- summary(finalmod)$coefficients[,4]

mymat <- round(cbind(coefs,coefs.se,oddsrat,or.se,pvals),3)

tabmat <- matrix(0,nrow(mymat),3)
tabmat[,1] <- paste0(mymat[,1]," (",mymat[,2],")")
tabmat[,2] <- paste0(mymat[,3]," (",mymat[,4],",",mymat[,5],")")
tabmat[,3] <- mymat[,6]

rownames(tabmat) <- c("Intercept","Total blood products","Mech support","Race (AA)","Race (H/A/NR)","Diagnosis (Dilated CM)","Diagnosis (Failed SV)","Diagnosis (Restrictive CM)","Diagnosis (Other)")
kable(tabmat, "latex", booktabs = T, align = "l", col.names = c("Estimate (SE)","Odds Ratio (95% CI)","p value"), caption = "Final Model") %>%
  kable_styling(latex_options = "hold_position")

library(pROC)
test_prob = predict(finalmod, type = "response")
test_roc = roc(df.blood$pra_bin ~ test_prob, plot = TRUE, print.auc = TRUE)

```

```{r modelparts, eval=FALSE}
mod2 <- glm(pra_bin ~ rbc_bin + plt_bin + ffp_bin + cryo_bin + prior_mech_nomiss + race2 + diag_c, data=df.blood, family="binomial",na.action = na.pass)

coefs <- summary(mod2)$coefficients[,1]
coefs.se <- summary(mod2)$coefficients[,2]
oddsrat <- exp(coefs)
or.se <- exp(confint(mod2))
pvals <- summary(mod2)$coefficients[,4]

mymat <- round(cbind(coefs,coefs.se,oddsrat,or.se,pvals),3)

tabmat <- matrix(0,nrow(mymat),3)
tabmat[,1] <- paste0(mymat[,1]," (",mymat[,2],")")
tabmat[,2] <- paste0(mymat[,3]," (",mymat[,4],",",mymat[,5],")")
tabmat[,3] <- mymat[,6]

rownames(tabmat) <- c("Intercept","Red blood cells","Platelets","Fresh frozen plasma","Cryoprecipitate","Mech support","Race (AA)","Race (H/A/NR)","Diagnosis (Dilated CM)","Diagnosis (Failed SV)","Diagnosis (Restrictive CM)","Diagnosis (Other)")
kable(tabmat, "latex", booktabs = T, align = "l", col.names = c("Estimate (SE)","Odds Ratio (95% CI)","p value"), caption = "Model with Individual Blood Components") %>%
  kable_styling(latex_options = "hold_position")

test_prob2 = predict(mod2, type = "response")
test_roc2 = roc(df.blood$pra_bin ~ test_prob2, plot = TRUE, print.auc = TRUE)



```