Reliability of LDF in determining the the pulp vitality

Ibrahim Hassan

2024-01-23

Load the necessary libraries

library(gtsummary)
library(tidyverse)
library(readxl)
library(corrr)
library(skimr)
library(pROC)
library(ggpubr)
library(flextable)
library(gt)
library(scales)
library(xtable)
library(pillar)
library(Hmisc)
library(kableExtra)
library(paletteer)
library(ggthemes)
library(cutpointr)
library(ggsci)
library(hrbrthemes)
library(easystats)
library(reshape2)
library(psych)
library(randomForest)

Look at the data

glimpse(df)
## Rows: 82
## Columns: 12
## $ Patient         <chr> "Trauma", "Trauma", "Trauma", "Trauma", "Trauma", "Tra…
## $ N               <chr> "C1", "C1.1", "E1", "E1.1", "J1", "J1.1", "T1", "T1.1"…
## $ Letter          <chr> "C", "C", "E", "E", "J", "J", "T", "T", "F", "F", "J",…
## $ number          <chr> "1", "1.1", "1", "1.1", "1", "1.1", "1", "1.1", "1", "…
## $ age             <dbl> 12, 12, 12, 12, 14, 14, 13, 13, 30, 30, 14, 14, 56, 56…
## $ gender          <chr> "Female", "Female", "Female", "Female", "Female", "Fem…
## $ Tooth           <chr> "UL1", "UL1", "UR1", "UR1", "UL2", "UL2", "UL1", "UL1"…
## $ LDF             <dbl> 8.2, 7.7, 5.0, 4.8, 3.0, 2.6, 0.5, 0.4, 0.9, 9.8, 2.2,…
## $ S.Dev           <dbl> 1.0, 0.9, 0.8, 0.9, 1.4, 0.2, 0.5, 0.3, 0.4, 1.5, 0.8,…
## $ S.Err           <dbl> 0.03, 0.03, 0.02, 0.03, 0.04, 0.01, 0.02, 0.01, 0.01, …
## $ `LDF Max range` <dbl> 11.2, 10.1, 6.9, 8.1, 10.0, 3.2, 4.2, 2.1, 2.7, 15.2, …
## $ `Min range`     <dbl> 5.8, 5.0, 2.9, 2.8, 0.4, 2.0, 0.0, 0.0, 0.1, 6.0, 12.1…
glimpse(val)
## Rows: 41
## Columns: 2
## $ LDF1 <dbl> 8.2, 5.0, 3.0, 0.5, 0.9, 2.2, 8.5, 6.6, 7.9, 8.9, 8.1, 7.2, 8.7, …
## $ LDF2 <dbl> 7.7, 4.8, 2.6, 0.4, 9.8, 1.1, 7.9, 6.5, 10.3, 10.2, 7.8, 7.2, 6.9…

Characteristics table

df %>%
  select(Patient,age,gender, LDF,Tooth) %>%
  tbl_summary(
    by = Tooth,
    type = all_continuous() ~ "continuous2",
    digits = list(all_continuous2() ~ 2, 
                  all_categorical() ~ c(0,2)),
    statistic = all_continuous2() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})"),
    label = list(age ~ "Age", gender ~ "Gender", LDF ~ "Mean LDF")
    ) %>%
    add_p(list(
        all_continuous2() ~ "aov",
        all_categorical() ~ "chisq.test"
        )
    ) %>%
    bold_p( t = 0.05, q = FALSE)
Characteristic UL1, N = 381 UL2, N = 21 UR1, N = 401 UR2, N = 21 p-value2
Patient



0.054
    Control 26 (68.42%) 0 (0.00%) 32 (80.00%) 2 (100.00%)
    Negative 4 (10.53%) 0 (0.00%) 4 (10.00%) 0 (0.00%)
    Trauma 8 (21.05%) 2 (100.00%) 4 (10.00%) 0 (0.00%)
Age



0.5
    Range 12.00 - 59.00 14.00 - 14.00 12.00 - 59.00 14.00 - 14.00
    Mean (±SD) 25.05 (±13.86) 14.00 (±0.00) 24.50 (±13.72) 14.00 (±0.00)
    Median (IQR) 26.00 (13.00, 29.50) 14.00 (14.00, 14.00) 22.50 (13.00, 28.50) 14.00 (14.00, 14.00)
Gender



0.5
    Female 24 (63.16%) 2 (100.00%) 26 (65.00%) 2 (100.00%)
    Male 14 (36.84%) 0 (0.00%) 14 (35.00%) 0 (0.00%)
Mean LDF



0.012
    Range 0.40 - 27.60 2.60 - 3.00 0.40 - 10.30 9.60 - 17.80
    Mean (±SD) 7.87 (±6.79) 2.80 (±0.28) 4.76 (±3.34) 13.70 (±5.80)
    Median (IQR) 7.45 (2.58, 9.10) 2.80 (2.70, 2.90) 5.05 (1.18, 7.75) 13.70 (11.65, 15.75)
1 n (%)
2 Pearson’s Chi-squared test; One-way ANOVA

Boxplot for Age between patient groups

df %>% 
   ggboxplot(
    y = "age", 
    rug = TRUE, 
    x = "Patient", 
    fill = "Patient", 
    ylab = "Age",
    xlab = "",
    title = "", 
    palette = "npg", 
    ggtheme = theme_pubr(),
    bxp.errorbar = TRUE
  ) +
    geom_pwc(method = "t.test", label = "p.adj.format",  bracket.nudge.y = -0.08)+
    stat_compare_means(method = "aov",label.y.npc = .7, label.x.npc = .5) +
    scale_fill_jama()+
    labs_pubr()

Boxplot for mean LDF between patient groups

df %>% 
   ggboxplot(
    y = "LDF", 
    rug = TRUE, 
    x = "Patient", 
    fill = "Patient", 
    ylab = "Mean LDF",
    xlab = "",
    title = "", 
    palette = "npg", 
    ggtheme = theme_pubr(),
    bxp.errorbar = TRUE
  ) +
    geom_pwc(method = "t.test", label = "p.adj.format",  bracket.nudge.y = -0.08)+
    stat_compare_means(method = "aov",label.y.npc = .7, label.x.npc = .5 ) +
    scale_fill_jama()+
    labs_pubr()

Reliability test: Cronbach’s alpha

alpha(val)[["total"]][["raw_alpha"]]
## [1] 0.8830882

ROC curve between Trauma and control

New dataframe

TC <-
    df %>%
    filter(Patient != "Negative") %>%
    mutate(P = case_when(Patient == "Trauma" ~ 1, .default = 0))

ROC plot

pROC::roc(
    TC$Patient,
    pos_class = "Trauma",
    direction= ">",
    TC$LDF,
    plot = TRUE,

    legacy.axes=TRUE,
    print.auc = TRUE,
    lwd = 4,
    col = "steelblue",
    ) 

## 
## Call:
## roc.default(response = TC$Patient, predictor = TC$LDF, direction = ">",     plot = TRUE, pos_class = "Trauma", legacy.axes = TRUE, print.auc = TRUE,     lwd = 4, col = "steelblue")
## 
## Data: TC$LDF in 60 controls (TC$Patient Control) > 14 cases (TC$Patient Trauma).
## Area under the curve: 0.5923

ROC plot (False & True Positive Percentages)

pROC::roc(
    TC$Patient,
    TC$LDF,
    pos_class = "Trauma",
    direction= ">",
    plot = TRUE,
    legacy.axes=TRUE,
    print.auc = TRUE,
    percent = TRUE,
    lwd = 4,
    col = "steelblue",
    xlab="False Positive Percentage",
    ylab="True Postive Percentage"
    )

## 
## Call:
## roc.default(response = TC$Patient, predictor = TC$LDF, percent = TRUE,     direction = ">", plot = TRUE, pos_class = "Trauma", legacy.axes = TRUE,     print.auc = TRUE, lwd = 4, col = "steelblue", xlab = "False Positive Percentage",     ylab = "True Postive Percentage")
## 
## Data: TC$LDF in 60 controls (TC$Patient Control) > 14 cases (TC$Patient Trauma).
## Area under the curve: 59.23%

Optimal cutoff point

cp_TC <- 
    cutpointr(
        TC,
        LDF,
        Patient,
        pos_class = "Trauma",
        direction= "<=",
        method = maximize_metric,
        metric = sum_sens_spec
        )
summary(cp_TC)
## Method: maximize_metric 
## Predictor: LDF 
## Outcome: Patient 
## Direction: <= 
## 
##     AUC  n n_pos n_neg
##  0.5923 74    14    60
## 
##  optimal_cutpoint sum_sens_spec    acc sensitivity specificity tp fn fp tn
##                 5        1.3262 0.6757      0.6429      0.6833  9  5 19 41
## 
## Predictor summary: 
##     Data Min.  5% 1st Qu. Median     Mean 3rd Qu.   95% Max.       SD NAs
##  Overall  0.4 0.4   2.525   6.55 6.736486   8.775 17.80 27.6 5.652823   0
##  Control  0.4 0.4   2.650   7.00 6.640000   8.825 17.80 18.0 4.818791   0
##   Trauma  0.4 0.4   2.225   4.85 7.150000   7.700 26.04 27.6 8.592013   0

ROC curve between Negative and control

New dataframe

NC <-
    df %>%
    filter(Patient != "Trauma") %>%
    mutate(P = case_when(Patient == "Negative" ~ 1, .default = 0))

ROC plot (Sensitivity & Specificity)

pROC::roc(
    NC$Patient,
    NC$LDF,
    pos_class = "Negative",
    direction= ">",
    plot = TRUE,
    legacy.axes=TRUE,
    print.auc = TRUE,
    lwd = 4,
    col = "steelblue"
    ) 

## 
## Call:
## roc.default(response = NC$Patient, predictor = NC$LDF, direction = ">",     plot = TRUE, pos_class = "Negative", legacy.axes = TRUE,     print.auc = TRUE, lwd = 4, col = "steelblue")
## 
## Data: NC$LDF in 60 controls (NC$Patient Control) > 8 cases (NC$Patient Negative).
## Area under the curve: 0.6865

ROC plot (False & True Positive Percentages)

pROC::roc(
    NC$Patient,
    NC$LDF,
    pos_class = "Negative",
    direction= ">",
    plot = TRUE,
    legacy.axes=TRUE,
    print.auc = TRUE,
    percent = TRUE,
    lwd = 4,
    col = "steelblue",
    xlab="False Positive Percentage",
    ylab="True Postive Percentage"
    )

## 
## Call:
## roc.default(response = NC$Patient, predictor = NC$LDF, percent = TRUE,     direction = ">", plot = TRUE, pos_class = "Negative", legacy.axes = TRUE,     print.auc = TRUE, lwd = 4, col = "steelblue", xlab = "False Positive Percentage",     ylab = "True Postive Percentage")
## 
## Data: NC$LDF in 60 controls (NC$Patient Control) > 8 cases (NC$Patient Negative).
## Area under the curve: 68.65%

Optimal cutoff point

cp_NC <- 
    cutpointr(
        NC,
        LDF,
        Patient,
        direction= "<=",
        pos_class = "Negative",
        method = maximize_metric,
        metric = sum_sens_spec
        )
summary(cp_NC)
## Method: maximize_metric 
## Predictor: LDF 
## Outcome: Patient 
## Direction: <= 
## 
##     AUC  n n_pos n_neg
##  0.6865 68     8    60
## 
##  optimal_cutpoint sum_sens_spec    acc sensitivity specificity tp fn fp tn
##               4.2        1.5583 0.7059       0.875      0.6833  7  1 19 41
## 
## Predictor summary: 
##      Data Min.    5% 1st Qu. Median     Mean 3rd Qu.    95% Max.       SD NAs
##   Overall  0.4 0.435   2.025   6.25 6.211765   8.800 17.625 18.0 4.771480   0
##   Control  0.4 0.400   2.650   7.00 6.640000   8.825 17.800 18.0 4.818791   0
##  Negative  0.9 0.970   1.175   1.85 3.000000   3.375  7.840  9.8 2.972253   0