#INTRODUCTION# Install packages

# Loading required libraries
if (!requireNamespace("readr", quietly = TRUE)) {
  install.packages("readr")
}
if (!requireNamespace("dplyr", quietly = TRUE)) {
  install.packages("dplyr")
}
if (!requireNamespace("kableExtra", quietly = TRUE)) {
  install.packages("kableExtra")
}
if (!requireNamespace("gt", quietly = TRUE)) {
  install.packages("gt")
}
if (!requireNamespace("flextable", quietly = TRUE)) {
  install.packages("flextable")
}
if (!requireNamespace("tidyverse", quietly = TRUE)) {
  install.packages("tidyverse")
}
if (!requireNamespace("broom", quietly = TRUE)) {
  install.packages("broom")
}
if (!requireNamespace("pscl", quietly = TRUE)) {
  install.packages("pscl")
}

library(readr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(gt)
library(flextable)
## 
## Attaching package: 'flextable'
## The following objects are masked from 'package:kableExtra':
## 
##     as_image, footnote
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::compose()         masks flextable::compose()
## ✖ dplyr::filter()          masks stats::filter()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ dplyr::lag()             masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom)
library(pscl)
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.

Import data

# Importing the data
data <- read_csv("MER_Diagnostic.csv")
## Rows: 640 Columns: 33
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (29): WorkingHours, LocationKNH, LocaWardAEOther, LocationReporting, Loc...
## dbl  (4): IDNumber, form_number, PatientAge, Agegroup2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Create a group vector based on the diagnostic variable

# Create a group vector based on the diagnostic variable
data <- data %>%
  mutate(Group = ifelse(Diagnostic == "YES", "Diagnostics", "Nondiagnostic"))

#TABLE 1 PATIENT CHARACTERISTICS#

data$Group <- as.factor(data$Group)
data$PatientGender <- as.factor(data$PatientGender)
data$Agegroup2 <- as.factor(data$Agegroup2)
data$PatientDiagnosis <- as.factor(data$PatientDiagnosis)
data$Whoisreporting <- as.factor(data$Whoisreporting)
data$Outcome <- as.factor(data$Outcome)

# Contingency table for Agegroup2
contingency_Agetable <- table(data$Agegroup2, data$Group)
print(contingency_Agetable)
##    
##     Diagnostics Nondiagnostic
##   1          21            44
##   2          72            83
##   3          98           143
##   4          22            31
##   5           2            14
# Chi-square test of independence for Agegroup2
chi_square_test_age <- chisq.test(contingency_Agetable)
print(chi_square_test_age)
## 
##  Pearson's Chi-squared test
## 
## data:  contingency_Agetable
## X-squared = 9.3135, df = 4, p-value = 0.05372
# Contingency table for PatientGender
contingency_PatientGendertable <- table(data$PatientGender, data$Group)
print(contingency_PatientGendertable)
##         
##          Diagnostics Nondiagnostic
##   Female         105           169
##   Male           135           189
# Chi-square test of independence for PatientGender
chi_square_test_gender <- chisq.test(contingency_PatientGendertable)
print(chi_square_test_gender)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_PatientGendertable
## X-squared = 0.55932, df = 1, p-value = 0.4545
# Contingency table for PatientDiagnosis
contingency_PatientDiagnosistable <- table(data$PatientDiagnosis, data$Group)
print(contingency_PatientDiagnosistable)
##           
##            Diagnostics Nondiagnostic
##   Medical           82           108
##   Oncology          19            30
##   surgical         127           197
# Chi-square test of independence for PatientDiagnosis
chi_square_test_diagnosis <- chisq.test(contingency_PatientDiagnosistable)
print(chi_square_test_diagnosis)
## 
##  Pearson's Chi-squared test
## 
## data:  contingency_PatientDiagnosistable
## X-squared = 0.84558, df = 2, p-value = 0.6552
# Contingency table for Whoisreporting
contingency_Whoisreportingtable <- table(data$Whoisreporting, data$Group)
print(contingency_Whoisreportingtable)
##             
##              Diagnostics Nondiagnostic
##   Doctor               0             3
##   Nurse              248           366
##   Pharmasist           0             2
# Chi-square test of independence for Whoisreporting
chi_square_test_who <- chisq.test(contingency_Whoisreportingtable)
## Warning in chisq.test(contingency_Whoisreportingtable): Chi-squared
## approximation may be incorrect
print(chi_square_test_who)
## 
##  Pearson's Chi-squared test
## 
## data:  contingency_Whoisreportingtable
## X-squared = 3.3695, df = 2, p-value = 0.1855
# Contingency table for Outcome
contingency_Outcometable <- table(data$Outcome, data$Group)
print(contingency_Outcometable)
##       
##        Diagnostics Nondiagnostic
##   Good         210           329
##   Poor          46            55

#Table 2: Distribution by of error reporting in relation to location and time error happened#

add_totals_and_proportions <- function(contingency_table) {
  # Add totals for rows and columns
  contingency_with_totals <- addmargins(contingency_table, FUN = list(Total = sum))
  
  # Calculate row proportions
  row_proportions <- prop.table(contingency_table, margin = 1)
  
  list(
    table_with_totals = contingency_with_totals,
    row_proportions = row_proportions
  )
}

# Create and print the contingency table for LocaWardAEOther
contingency_LocaWardAEOthertable <- table(data$LocaWardAEOther, data$Group)
results_LocaWardAEOther <- add_totals_and_proportions(contingency_LocaWardAEOthertable)
## Margins computed over dimensions
## in the following order:
## 1: 
## 2:
print("Contingency Table for LocaWardAEOther with Totals:")
## [1] "Contingency Table for LocaWardAEOther with Totals:"
print(results_LocaWardAEOther$table_with_totals)
##        
##         Diagnostics Nondiagnostic Total
##   AE            102           122   224
##   Other          37            79   116
##   Ward          105           170   275
##   Total         244           371   615
print("Row Proportions for LocaWardAEOther:")
## [1] "Row Proportions for LocaWardAEOther:"
print(round(results_LocaWardAEOther$row_proportions, 2))
##        
##         Diagnostics Nondiagnostic
##   AE           0.46          0.54
##   Other        0.32          0.68
##   Ward         0.38          0.62
# Create and print the contingency table for WorkingHours
contingency_WorkingHourstable <- table(data$WorkingHours, data$Group)
results_WorkingHours <- add_totals_and_proportions(contingency_WorkingHourstable)
## Margins computed over dimensions
## in the following order:
## 1: 
## 2:
print("Contingency Table for WorkingHours with Totals:")
## [1] "Contingency Table for WorkingHours with Totals:"
print(results_WorkingHours$table_with_totals)
##        
##         Diagnostics Nondiagnostic Total
##   NO            131           153   284
##   YES           125           231   356
##   Total         256           384   640
print("Row Proportions for WorkingHours:")
## [1] "Row Proportions for WorkingHours:"
print(round(results_WorkingHours$row_proportions, 2))
##      
##       Diagnostics Nondiagnostic
##   NO         0.46          0.54
##   YES        0.35          0.65

#Table 3 Types of Medical Errors# ## Diagnostic, Medication, Treatment, Preventive and others##

# Calculate proportions for a binary variable
calculate_proportions <- function(df, variable) {
  prop_table <- prop.table(table(df[[variable]]))
  return(as.data.frame(prop_table))
}

# Calculate proportions for each error type
diagnostic_proportions <- calculate_proportions(data, "Diagnostic")
treatment_proportions <- calculate_proportions(data, "Treatment")
medication_proportions <- calculate_proportions(data, "Medication")
preventive_proportions <- calculate_proportions(data, "Preventive")
others_proportions <- calculate_proportions(data, "Others")

# Combine all proportions into one data frame
proportions_df <- data.frame(
  Type = c("Diagnostic", "Diagnostic", "Treatment", "Treatment", "Medication", "Medication", "Preventive", "Preventive", "Others", "Others"),
  Response = rep(c("NO", "YES"), 5),
  Proportion = c(
    diagnostic_proportions$Freq,
    treatment_proportions$Freq,
    medication_proportions$Freq,
    preventive_proportions$Freq,
    others_proportions$Freq
  )
)

# Create a flextable to display the proportions
flex_table <- flextable(proportions_df) %>%
  set_header_labels(
    Type = "Type of Error",
    Response = "Response",
    Proportion = "Proportion"
  ) %>%
  add_header_row(values = c("Univariate Analysis of Types of Errors"), colwidths = 3) %>%
  theme_vanilla() %>%
  autofit()

# Display the flextable
flex_table

Univariate Analysis of Types of Errors

Type of Error

Response

Proportion

Diagnostic

NO

0.6000000

Diagnostic

YES

0.4000000

Treatment

NO

0.4578125

Treatment

YES

0.5421875

Medication

NO

0.8218750

Medication

YES

0.1781250

Preventive

NO

0.9296875

Preventive

YES

0.0703125

Others

NO

0.7968750

Others

YES

0.2031250

#Figure 1 Sub-Classify Diagnostic errors#

# Filter out missing values and calculate the proportion of each diagnostic-specific error type
error_proportions <- data %>%
  filter(!is.na(DiagsticSpecific)) %>%
  count(DiagsticSpecific) %>%
  mutate(Proportion = n / sum(n) * 100)

# Create the bar graph using ggplot2
ggplot(error_proportions, aes(x = DiagsticSpecific, y = Proportion, fill = DiagsticSpecific)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = paste0(round(Proportion, 1), "%")), vjust = -0.5) +
  labs(
    title = "Proportion of Diagnostic-Specific Errors",
    x = "Type of Diagnostic Error",
    y = "Proportion (%)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

#Table 4 Contributing factors to diagnostic errors#

# Create a group vector based on the diagnostic variable
data <- data %>%
  mutate(Group = ifelse(Diagnostic == "YES", "Diagnostics", "Nondiagnostic"))

# Make the vectors as factors
data$Group <- as.factor(data$Group)
data$ContrCommPlus <- as.factor(data$ContrCommPlus)
data$contrPatient <- as.factor(data$contrPatient)
data$contStaff <- as.factor(data$contStaff)
data$CombiEquipment <- as.factor(data$CombiEquipment)
data$LackofPolicies <- as.factor(data$LackofPolicies)

# Contingency table for Communication
contingency_Communicatable <- table(data$ContrCommPlus, data$Group)
print(contingency_Communicatable)
##      
##       Diagnostics Nondiagnostic
##   NO          155           254
##   YES         101           130
# Chi-square test of independence for Communication
chi_square_test_Comm <- chisq.test(contingency_Communicatable)
print(chi_square_test_Comm)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_Communicatable
## X-squared = 1.8518, df = 1, p-value = 0.1736
# Contingency table for Patientfactor
contingency_Patientfactorstable <- table(data$contrPatient, data$Group)
print(contingency_Patientfactorstable)
##      
##       Diagnostics Nondiagnostic
##   NO          222           344
##   YES          34            40
# Chi-square test of independence for Patientfactor
chi_square_test_Patientfactor <- chisq.test(contingency_Patientfactorstable)
print(chi_square_test_Patientfactor)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_Patientfactorstable
## X-squared = 0.96839, df = 1, p-value = 0.3251
# Contingency table for Stafffactors
contingency_Stafffactorstable <- table(data$contStaff, data$Group)
print(contingency_Stafffactorstable)
##      
##       Diagnostics Nondiagnostic
##   NO          130           197
##   YES         126           187
# Chi-square test of independence for Stafffactors
chi_square_test_Stafffactors <- chisq.test(contingency_Stafffactorstable)
print(chi_square_test_Stafffactors)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_Stafffactorstable
## X-squared = 0.0023449, df = 1, p-value = 0.9614
# Contingency table for Equipment factors
contingency_Equipmenttable <- table(data$CombiEquipment, data$Group)
print(contingency_Equipmenttable)
##      
##       Diagnostics Nondiagnostic
##   NO          221           319
##   YES          35            65
# Chi-square test of independence for Equipment factors
chi_square_test_Equipment <- chisq.test(contingency_Equipmenttable)
print(chi_square_test_Equipment)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_Equipmenttable
## X-squared = 1, df = 1, p-value = 0.3173
# Contingency table for Lack of Policies
contingency_LackofPoliciestable <- table(data$LackofPolicies, data$Group)
print(contingency_LackofPoliciestable)
##      
##       Diagnostics Nondiagnostic
##   NO          192           311
##   YES          64            73
# Chi-square test of independence for Lack of Policies
chi_square_test_LackofPolicies <- chisq.test(contingency_LackofPoliciestable)
print(chi_square_test_LackofPolicies)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_LackofPoliciestable
## X-squared = 2.929, df = 1, p-value = 0.087

#Table 5 Assocoation of factors with likelyhood of reporting diagnostic errors#

# Convert the diagnostic variable to a binary factor
data <- data %>%
  mutate(DiagnosticError = ifelse(Diagnostic == "YES", 1, 0))

# Set the reference levels
data$LocaWardAEOther <- relevel(factor(data$LocaWardAEOther), ref = "Ward")
data$WorkingHours <- relevel(factor(data$WorkingHours), ref = "YES")
data$LocationKNH <- relevel(factor(data$LocationKNH), ref = "NO")

# Perform logistic regression
logistic_model <- glm(DiagnosticError ~ LocaWardAEOther + WorkingHours + LocationKNH, 
                      data = data, 
                      family = binomial)

# Summarize the logistic regression results
summary(logistic_model)
## 
## Call:
## glm(formula = DiagnosticError ~ LocaWardAEOther + WorkingHours + 
##     LocationKNH, family = binomial, data = data)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -1.4904     0.3315  -4.496 6.93e-06 ***
## LocaWardAEOtherAE      0.4630     0.2022   2.289  0.02205 *  
## LocaWardAEOtherOther  -0.2838     0.2355  -1.205  0.22831    
## WorkingHoursNO         0.3408     0.1718   1.984  0.04731 *  
## LocationKNHYES         0.8935     0.3014   2.964  0.00303 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 826.16  on 614  degrees of freedom
## Residual deviance: 806.76  on 610  degrees of freedom
##   (25 observations deleted due to missingness)
## AIC: 816.76
## 
## Number of Fisher Scoring iterations: 4
# Calculate odds ratios and confidence intervals
logistic_model_results <- tidy(logistic_model, exponentiate = TRUE, conf.int = TRUE)
print(logistic_model_results)
## # A tibble: 5 × 7
##   term                 estimate std.error statistic   p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)             0.225     0.332     -4.50   6.93e-6    0.116     0.426
## 2 LocaWardAEOtherAE       1.59      0.202      2.29   2.21e-2    1.07      2.37 
## 3 LocaWardAEOtherOther    0.753     0.236     -1.20   2.28e-1    0.471     1.19 
## 4 WorkingHoursNO          1.41      0.172      1.98   4.73e-2    1.00      1.97 
## 5 LocationKNHYES          2.44      0.301      2.96   3.03e-3    1.37      4.48
# Calculate McFadden's R-squared
r2_mcfadden <- pR2(logistic_model)["McFadden"]
## fitting null model for pseudo-r2
cat("McFadden's R-squared: ", r2_mcfadden, "\n")
## McFadden's R-squared:  0.02347439

OUTPUT “McFadden’s R-squared ###cat(“McFadden’s R-squared:”, r2_mcfadden, “”)### McFadden’s R-squared: 0.02347439

#END# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX