Read Data from RedCap

  • Read data in R and an R dataframe named newdata

  • Data Variables used to score GRAIL

    • age_admission (numeric)
    • sex (categorical; Options: Male, Female)
    • alb_admit (numeric)
    • creatinine_admission (numeric)
    • bun_admit (numeric)

  • Limits on Labs (Check the data)

    • 0 < Cr < 0.1 then Cr=0.1; Cr > 10 then Cr=10
    • 0 < BUN < 1 then BUN=1; BUN > 219 then BUN=219
    • 0 < Alb < 1 then Alb=1; Alb > 6 then Alb=6

Create work Directory

Save GRAIL Model (R Object)

  • Save the file “Equationgrail.rda” in the work directory

GRAIL Model Fitted

  • Read saved fitted model as R object from a file
  • Print ANOVA Table

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = mGFR_cap ~ Age_mGFR + Sex + logAlb + rcs(logCr, 
    3) + rcs(logBUN, 3), data = df, x = TRUE, y = TRUE, scale = TRUE)
Model Likelihood
Ratio Test		 Discrimination
Indexes		 Rank Discrim.
Indexes
Obs 1235 LR χ2 1275.56 R2 0.644 ρ 0.770
Distinct Y 151 d.f. 7 R21235 0.644
Y0.5 96 Pr(>χ2) <0.0001 R27,1235 0.642
max |∂log L/∂β| 0.0004 Score χ2 1557.72 R21234.3 0.644
Pr(>χ2) <0.0001 R27,1234.3 0.642
g 2.691
gr 14.746
|Pr(Y ≥ median)-½| 0.287
β S.E. Wald Z Pr(>|Z|)
Age_mGFR  -0.0455  0.0057 -8.00 <0.0001
Sex=Male   1.3339  0.1121 11.90 <0.0001
logAlb   1.2588  0.2556 4.93 <0.0001
logCr  -4.6981  0.4111 -11.43 <0.0001
logCr'  -0.3406  0.6373 -0.53 0.5931
logBUN  -0.5138  0.2663 -1.93 0.0537
logBUN'  -1.2204  0.3697 -3.30 0.0010

Calculate the Score from the Model

  • Applied to saved model to the newdata to calculate the score
png 
  2 
png 
  2 
png 
  2 
fit6 <- glm(outcome1_kidney~Q10_GRAIL_30,data = Grail_Patients_Data_E0,family = "binomial")
prob_fit6 <- predict(fit6, type = "response")


dt<-Grail_Patients_Data_E0
dt$prob_fit6 <- prob_fit6

roc_mod6 <- roc(response=dt$outcome1_kidney, predictor=prob_fit6)
roc_mod6_AUC<-round(roc_mod6$auc, digits = 3)

g6 <- ggplot(dt, aes(d = outcome1_kidney, m = prob_fit6)) +   ##### prob_fit6 was prob_fit5 CHECK #########
  geom_roc(n.cuts = 0, color = "#0071BF") +
  style_roc() +
  theme(text=element_text(size=20, face="bold"), #change font size of all text
        axis.text=element_text(size=20), #change font size of axis text
        axis.title=element_text(size=20, face="bold"), #change font size of axis titles
        plot.title=element_text(size=14, face="bold"), #change font size of plot title
        legend.text=element_text(size=20, face="bold"), #change font size of legend text
        legend.title=element_text(size=20)) + #change font size of legend title
  geom_abline(intercept = 0, slope = 1, linetype = 'dashed') +
  labs(x = "1 - Specificity", y = "Sensitivity") 

g6<- g6 + annotate("text", x=0.40, y=0.95, size=6, fontface="bold",
                   label=paste("AUC Q10-GRAIL (Cutoff 30) = ", roc_mod6_AUC))



png(filename="C:\\Users\\to909\\Desktop\\Projects\\GRAIL code\\AUC Q10-GRAIL (Cutoff 30) new.png")
g6
dev.off()
png 
  2

renal event

fit_renal <- glm(renal_event~Q10_GRAIL_30,data = Grail_Patients_Data_E0,family = "binomial")
prob_fit_renal <- predict(fit_renal, type = "response")


dt<-Grail_Patients_Data_E0
dt$prob_fit_renal <- prob_fit_renal

roc_mod_renal <- roc(response=dt$renal_event, predictor=prob_fit_renal)
roc_mod_renal_AUC<-round(roc_mod_renal$auc, digits = 3)

g_renal <- ggplot(dt, aes(d = renal_event, m = prob_fit_renal)) +   ##### prob_fit_renal was prob_fit5 CHECK #########
  geom_roc(n.cuts = 0, color = "#0071BF") +
  style_roc() +
  theme(text=element_text(size=20, face="bold"), #change font size of all text
        axis.text=element_text(size=20), #change font size of axis text
        axis.title=element_text(size=20, face="bold"), #change font size of axis titles
        plot.title=element_text(size=14, face="bold"), #change font size of plot title
        legend.text=element_text(size=20, face="bold"), #change font size of legend text
        legend.title=element_text(size=20)) + #change font size of legend title
  geom_abline(intercept = 0, slope = 1, linetype = 'dashed') +
  labs(x = "1 - Specificity", y = "Sensitivity") 

g_renal<- g_renal + annotate("text", x=0.40, y=0.95, size=6, fontface="bold",
                   label=paste("AUC Q10-GRAIL (Cutoff 30) = ", roc_mod_renal_AUC))



png(filename="C:\\Users\\to909\\Desktop\\Projects\\GRAIL code\\AUC Q10-GRAIL (Cutoff 30) renal.png")
g_renal
dev.off()
png 
  2 
library(epiR)
tab1 = as.table(table(Grail_Patients_Data_E0$Q10_GRAIL_30, Grail_Patients_Data_E0$outcome1_kidney))
tab1
     
         0    1
  >30  298   27
  ≤30 1308  191
Grail_Patients_Data_E0$Q10_GRAIL_30<-factor(Grail_Patients_Data_E0$Q10_GRAIL_30,
  levels = c("≤30", ">30"), ordered = TRUE) 
Grail_Patients_Data_E0$outcome1_kidney<-factor(Grail_Patients_Data_E0$outcome1_kidney,
  levels = c("1", "0"), ordered = TRUE) 

tab1 = as.table(table(Grail_Patients_Data_E0$Q10_GRAIL_30, Grail_Patients_Data_E0$outcome1_kidney))
tab1
     
         1    0
  ≤30  191 1308
  >30   27  298
rval <- epi.tests(tab1, method="exact", digits = 3, conf.level = 0.95)
print(rval)
          Outcome +    Outcome -      Total
Test +          191         1308       1499
Test -           27          298        325
Total           218         1606       1824

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.822 (0.803, 0.839)
True prevalence *                      0.120 (0.105, 0.135)
Sensitivity *                          0.876 (0.825, 0.917)
Specificity *                          0.186 (0.167, 0.205)
Positive predictive value *            0.127 (0.111, 0.145)
Negative predictive value *            0.917 (0.881, 0.945)
Positive likelihood ratio              1.076 (1.018, 1.137)
Negative likelihood ratio              0.667 (0.462, 0.964)
False T+ proportion for true D- *      0.814 (0.795, 0.833)
False T- proportion for true D+ *      0.124 (0.083, 0.175)
False T+ proportion for T+ *           0.873 (0.855, 0.889)
False T- proportion for T- *           0.083 (0.055, 0.119)
Correctly classified proportion *      0.268 (0.248, 0.289)
--------------------------------------------------------------
* Exact CIs
#summary(rval)
library(epiR)

tab1 = as.table(table(Grail_Patients_Data_E0$Median_GRAIL_30, Grail_Patients_Data_E0$outcome1_kidney))
tab1
     
        1   0
  >30  52 685
  ≤30 166 921
Grail_Patients_Data_E0$Median_GRAIL_30<-factor(Grail_Patients_Data_E0$Median_GRAIL_30,
  levels = c("≤30", ">30"), ordered = TRUE) 
Grail_Patients_Data_E0$outcome1_kidney<-factor(Grail_Patients_Data_E0$outcome1_kidney,
  levels = c("1", "0"), ordered = TRUE) 

tab1 = as.table(table(Grail_Patients_Data_E0$Median_GRAIL_30, Grail_Patients_Data_E0$outcome1_kidney))
tab1
     
        1   0
  ≤30 166 921
  >30  52 685
rval <- epi.tests(tab1, method="exact", digits = 3, conf.level = 0.95)
print(rval)
          Outcome +    Outcome -      Total
Test +          166          921       1087
Test -           52          685        737
Total           218         1606       1824

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.596 (0.573, 0.619)
True prevalence *                      0.120 (0.105, 0.135)
Sensitivity *                          0.761 (0.699, 0.816)
Specificity *                          0.427 (0.402, 0.451)
Positive predictive value *            0.153 (0.132, 0.175)
Negative predictive value *            0.929 (0.909, 0.947)
Positive likelihood ratio              1.328 (1.219, 1.446)
Negative likelihood ratio              0.559 (0.438, 0.714)
False T+ proportion for true D- *      0.573 (0.549, 0.598)
False T- proportion for true D+ *      0.239 (0.184, 0.301)
False T+ proportion for T+ *           0.847 (0.825, 0.868)
False T- proportion for T- *           0.071 (0.053, 0.091)
Correctly classified proportion *      0.467 (0.443, 0.490)
--------------------------------------------------------------
* Exact CIs
tab1 = as.table(table(Grail_Patients_Data_E0$CKDEpi.rf_30, Grail_Patients_Data_E0$outcome1_kidney))
tab1
     
        1   0
  >30  80 980
  ≤30 138 626
Grail_Patients_Data_E0$CKDEpi.rf_30<-factor(Grail_Patients_Data_E0$CKDEpi.rf_30,
  levels = c("≤30", ">30"), ordered = TRUE) 
Grail_Patients_Data_E0$outcome1_kidney<-factor(Grail_Patients_Data_E0$outcome1_kidney,
  levels = c("1", "0"), ordered = TRUE) 

tab1 = as.table(table(Grail_Patients_Data_E0$CKDEpi.rf_30, Grail_Patients_Data_E0$outcome1_kidney))
tab1
     
        1   0
  ≤30 138 626
  >30  80 980
rval <- epi.tests(tab1, method="exact", digits = 3, conf.level = 0.95)
print(rval)
          Outcome +    Outcome -      Total
Test +          138          626        764
Test -           80          980       1060
Total           218         1606       1824

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.419 (0.396, 0.442)
True prevalence *                      0.120 (0.105, 0.135)
Sensitivity *                          0.633 (0.565, 0.697)
Specificity *                          0.610 (0.586, 0.634)
Positive predictive value *            0.181 (0.154, 0.210)
Negative predictive value *            0.925 (0.907, 0.940)
Positive likelihood ratio              1.624 (1.443, 1.828)
Negative likelihood ratio              0.601 (0.503, 0.719)
False T+ proportion for true D- *      0.390 (0.366, 0.414)
False T- proportion for true D+ *      0.367 (0.303, 0.435)
False T+ proportion for T+ *           0.819 (0.790, 0.846)
False T- proportion for T- *           0.075 (0.060, 0.093)
Correctly classified proportion *      0.613 (0.590, 0.635)
--------------------------------------------------------------
* Exact CIs

#renal

tab1 = as.table(table(Grail_Patients_Data_E0$CKDEpi.rf_30, Grail_Patients_Data_E0$renal_event))
tab1
     
        0   1
  ≤30 494 270
  >30 862 198
Grail_Patients_Data_E0$CKDEpi.rf_30<-factor(Grail_Patients_Data_E0$CKDEpi.rf_30,
  levels = c("≤30", ">30"), ordered = TRUE) 
Grail_Patients_Data_E0$renal_event<-factor(Grail_Patients_Data_E0$renal_event,
  levels = c("1", "0"), ordered = TRUE) 

tab1 = as.table(table(Grail_Patients_Data_E0$CKDEpi.rf_30, Grail_Patients_Data_E0$renal_event))
tab1
     
        1   0
  ≤30 270 494
  >30 198 862
rval <- epi.tests(tab1, method="exact", digits = 3, conf.level = 0.95)
print(rval)
          Outcome +    Outcome -      Total
Test +          270          494        764
Test -          198          862       1060
Total           468         1356       1824

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.419 (0.396, 0.442)
True prevalence *                      0.257 (0.237, 0.277)
Sensitivity *                          0.577 (0.531, 0.622)
Specificity *                          0.636 (0.609, 0.661)
Positive predictive value *            0.353 (0.319, 0.388)
Negative predictive value *            0.813 (0.788, 0.836)
Positive likelihood ratio              1.584 (1.426, 1.758)
Negative likelihood ratio              0.666 (0.594, 0.745)
False T+ proportion for true D- *      0.364 (0.339, 0.391)
False T- proportion for true D+ *      0.423 (0.378, 0.469)
False T+ proportion for T+ *           0.647 (0.612, 0.681)
False T- proportion for T- *           0.187 (0.164, 0.212)
Correctly classified proportion *      0.621 (0.598, 0.643)
--------------------------------------------------------------
* Exact CIs
tab1 = as.table(table(Grail_Patients_Data_E0$Q10_GRAIL_30, Grail_Patients_Data_E0$renal_event))
tab1
     
         1    0
  ≤30  398 1101
  >30   70  255
Grail_Patients_Data_E0$Q10_GRAIL_30<-factor(Grail_Patients_Data_E0$Q10_GRAIL_30,
                                            levels = c("b   $30", ">30"), ordered = TRUE) 
Grail_Patients_Data_E0$renal_event<-factor(Grail_Patients_Data_E0$renal_event,
                                               levels = c("1", "0"), ordered = TRUE) 

tab1 = as.table(table(Grail_Patients_Data_E0$Q10_GRAIL_30, Grail_Patients_Data_E0$renal_event))
tab1
        
           1   0
  b\t$30   0   0
  >30     70 255
rval <- epi.tests(tab1, method="exact", digits = 3, conf.level = 0.95)
print(rval)
          Outcome +    Outcome -      Total
Test +            0            0          0
Test -           70          255        325
Total            70          255        325

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.000 (0.000, 0.011)
True prevalence *                      0.215 (0.172, 0.264)
Sensitivity *                          0.000 (0.000, 0.051)
Specificity *                          1.000 (0.986, 1.000)
Positive predictive value *            NaN (0.000, 1.000)
Negative predictive value *            0.785 (0.736, 0.828)
Positive likelihood ratio              NaN (NaN, NaN)
Negative likelihood ratio              1.000 (1.000, 1.000)
False T+ proportion for true D- *      0.000 (0.000, 0.014)
False T- proportion for true D+ *      1.000 (0.949, 1.000)
False T+ proportion for T+ *           NaN (0.000, 1.000)
False T- proportion for T- *           0.215 (0.172, 0.264)
Correctly classified proportion *      0.785 (0.736, 0.828)
--------------------------------------------------------------
* Exact CIs
library(epiR)

tab1 = as.table(table(Grail_Patients_Data_E0$Median_GRAIL_30, Grail_Patients_Data_E0$renal_event))
tab1
     
        1   0
  ≤30 327 760
  >30 141 596
Grail_Patients_Data_E0$Median_GRAIL_30<-factor(Grail_Patients_Data_E0$Median_GRAIL_30,
  levels = c("≤30", ">30"), ordered = TRUE) 
Grail_Patients_Data_E0$renal_event<-factor(Grail_Patients_Data_E0$renal_event,
  levels = c("1", "0"), ordered = TRUE) 

tab1 = as.table(table(Grail_Patients_Data_E0$Median_GRAIL_30, Grail_Patients_Data_E0$renal_event))
tab1
     
        1   0
  ≤30 327 760
  >30 141 596
rval <- epi.tests(tab1, method="exact", digits = 3, conf.level = 0.95)
print(rval)
          Outcome +    Outcome -      Total
Test +          327          760       1087
Test -          141          596        737
Total           468         1356       1824

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.596 (0.573, 0.619)
True prevalence *                      0.257 (0.237, 0.277)
Sensitivity *                          0.699 (0.655, 0.740)
Specificity *                          0.440 (0.413, 0.466)
Positive predictive value *            0.301 (0.274, 0.329)
Negative predictive value *            0.809 (0.778, 0.836)
Positive likelihood ratio              1.247 (1.156, 1.345)
Negative likelihood ratio              0.685 (0.590, 0.797)
False T+ proportion for true D- *      0.560 (0.534, 0.587)
False T- proportion for true D+ *      0.301 (0.260, 0.345)
False T+ proportion for T+ *           0.699 (0.671, 0.726)
False T- proportion for T- *           0.191 (0.164, 0.222)
Correctly classified proportion *      0.506 (0.483, 0.529)
--------------------------------------------------------------
* Exact CIs
tab1 = as.table(table(Grail_Patients_Data_E0$CKDEpi.rf_30, Grail_Patients_Data_E0$renal_event))
tab1
     
        1   0
  ≤30 270 494
  >30 198 862
Grail_Patients_Data_E0$CKDEpi.rf_30<-factor(Grail_Patients_Data_E0$CKDEpi.rf_30,
  levels = c("≤30", ">30"), ordered = TRUE) 
Grail_Patients_Data_E0$renal_event<-factor(Grail_Patients_Data_E0$renal_event,
  levels = c("1", "0"), ordered = TRUE) 

tab1 = as.table(table(Grail_Patients_Data_E0$CKDEpi.rf_30, Grail_Patients_Data_E0$renal_event))
tab1
     
        1   0
  ≤30 270 494
  >30 198 862
rval <- epi.tests(tab1, method="exact", digits = 3, conf.level = 0.95)
print(rval)
          Outcome +    Outcome -      Total
Test +          270          494        764
Test -          198          862       1060
Total           468         1356       1824

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.419 (0.396, 0.442)
True prevalence *                      0.257 (0.237, 0.277)
Sensitivity *                          0.577 (0.531, 0.622)
Specificity *                          0.636 (0.609, 0.661)
Positive predictive value *            0.353 (0.319, 0.388)
Negative predictive value *            0.813 (0.788, 0.836)
Positive likelihood ratio              1.584 (1.426, 1.758)
Negative likelihood ratio              0.666 (0.594, 0.745)
False T+ proportion for true D- *      0.364 (0.339, 0.391)
False T- proportion for true D+ *      0.423 (0.378, 0.469)
False T+ proportion for T+ *           0.647 (0.612, 0.681)
False T- proportion for T- *           0.187 (0.164, 0.212)
Correctly classified proportion *      0.621 (0.598, 0.643)
--------------------------------------------------------------
* Exact CIs
table(master$outcome1_kidney,useNA = "always")

   0    1 <NA> 
1667  229  145 
table(master$outcome1_kidney1,useNA = "always")

   0    1 <NA> 
2013   27    1 
table(master$outcome1_kidney2,useNA = "always")

   0    1 <NA> 
1685  205  151 
table(master$renal_event,useNA = "always")

   0    1 <NA> 
1499  542    0 
library(gtsummary)

names(newdata)
 [1] "Age_mGFR"             "sex"                  "alb_admit"           
 [4] "creatinine_admission" "bun_admit"            "outcome1_kidney"     
 [7] "MELD_Na_baseline"     "outcome1_kidney1"     "outcome1_kidney2"    
[10] "renal_event"          "logAlb"               "logCr"               
[13] "logBUN"               "Sex"                  "SEXn"                
[16] "Median_GRAIL_02"      "Q10_GRAIL_02"         "CKDEpi.rf"           
[19] "CKDEpi.rf_15"         "CKDEpi.rf_20"         "CKDEpi.rf_30"        
[22] "Median_GRAIL_15"      "Median_GRAIL_20"      "Median_GRAIL_30"     
[25] "Q10_GRAIL_10"         "Q10_GRAIL_15"         "Q10_GRAIL_20"        
[28] "Q10_GRAIL_30"        
library(table1)
table1(~Age_mGFR+Sex+Median_GRAIL_02+MELD_Na_baseline+Q10_GRAIL_02+CKDEpi.rf+as.factor(renal_event),data=newdata)
Overall
(N=2041)
Age_mGFR
Mean (SD) 61.0 (12.0)
Median [Min, Max] 62.0 [23.0, 96.0]
Sex
Female 779 (38.2%)
Male 1262 (61.8%)
Median_GRAIL_02
Mean (SD) 32.0 (26.7)
Median [Min, Max] 26.3 [2.00, 150]
MELD_Na_baseline
Mean (SD) 25.4 (8.52)
Median [Min, Max] 26.0 [-4.00, 52.0]
Missing 62 (3.0%)
Q10_GRAIL_02
Mean (SD) 18.5 (19.1)
Median [Min, Max] 12.2 [2.00, 124]
CKDEpi.rf
Mean (SD) 39.2 (24.1)
Median [Min, Max] 33.7 [2.67, 131]
as.factor(renal_event)
0 1499 (73.4%)
1 542 (26.6%)