introduction-to-qDEA.R
In qDEA: Quantile Data Envelopment Analysis

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----eval=FALSE---------------------------------------------------------------
# # Install from CRAN
# install.packages("qDEA")
# 
# # Load the package
# library(qDEA)

## ----echo=FALSE---------------------------------------------------------------
# For vignette building
library(qDEA)

## ----basic-dea----------------------------------------------------------------
# Load example data
data(CST11)
head(CST11)

# Prepare input and output matrices
X <- as.matrix(CST11$EMPLOYEES)
Y <- as.matrix(CST11$SALES_EJOR)

# Run output-oriented DEA with constant returns to scale
result_dea <- qDEA(X = X, Y = Y, 
                   orient = "out", 
                   RTS = "CRS",
                   qout = 0)  # Standard DEA (no outliers allowed)

# View efficiency scores
result_dea$effvals

## ----basic-qdea---------------------------------------------------------------
# Run qDEA allowing one outlier
qout <- 1/nrow(X)  # Proportion = 1/8 = 0.125

result_qdea <- qDEA(X = X, Y = Y, 
                    orient = "out", 
                    RTS = "CRS",
                    qout = qout)

# Compare DEA and qDEA efficiency scores
comparison <- data.frame(
  Store = CST11$STORE,
  DEA = round(result_dea$effvals, 3),
  qDEA = round(result_qdea$effvalsq, 3),
  Difference = round(result_qdea$effvalsq - result_dea$effvals, 3)
)
print(comparison)

## ----two-inputs---------------------------------------------------------------
# Load two-input, one-output data
data(CST21)
head(CST21)

# Prepare matrices
X <- as.matrix(CST21[, c("EMPLOYEES", "FLOOR_AREA")])
Y <- as.matrix(CST21$SALES)

# Input-oriented DEA with VRS
result_vrs <- qDEA(X = X, Y = Y,
                   orient = "in",
                   RTS = "VRS",
                   qout = 1/nrow(X))

# Display results
data.frame(
  Store = CST21$STORE,
  Employees = CST21$EMPLOYEES,
  Floor_Area = CST21$FLOOR_AREA,
  Sales = CST21$SALES,
  Efficiency = round(result_vrs$effvalsq, 3)
)

## ----two-outputs--------------------------------------------------------------
# Load one-input, two-output data
data(CST12)
head(CST12)

# Prepare matrices
X <- as.matrix(CST12$EMPLOYEES)
Y <- as.matrix(CST12[, c("CUSTOMERS", "SALES")])

# Output-oriented qDEA
result_outputs <- qDEA(X = X, Y = Y,
                       orient = "out",
                       RTS = "CRS",
                       qout = 0.15)  # Allow 15% outliers

# Results
data.frame(
  Store = CST12$STORE,
  Employees = CST12$EMPLOYEES,
  Customers = CST12$CUSTOMERS,
  Sales = CST12$SALES,
  DEA_Eff = round(result_outputs$effvals, 3),
  qDEA_Eff = round(result_outputs$effvalsq, 3)
)

## ----hospitals----------------------------------------------------------------
# Load hospital data
data(CST22)
head(CST22)

# Prepare matrices
X <- as.matrix(CST22[, c("DOCTORS", "NURSES")])
Y <- as.matrix(CST22[, c("OUT_PATIENTS", "IN_PATIENTS")])

# Run qDEA with 10% outliers allowed
result_hospital <- qDEA(X = X, Y = Y,
                        orient = "in",
                        RTS = "VRS",
                        qout = 0.10)

# Create results table
hospital_results <- data.frame(
  Hospital = CST22$HOSPITAL,
  Doctors = CST22$DOCTORS,
  Nurses = CST22$NURSES,
  Out_Patients = CST22$OUT_PATIENTS,
  In_Patients = CST22$IN_PATIENTS,
  Efficiency = round(result_hospital$effvalsq, 3)
)

print(hospital_results)

# Identify efficient hospitals
efficient <- hospital_results$Hospital[hospital_results$Efficiency >= 0.99]
cat("\nEfficient hospitals:", paste(efficient, collapse = ", "))

## ----input-orientation--------------------------------------------------------
result_in <- qDEA(X = X, Y = Y, orient = "in", RTS = "CRS", qout = 0.10)
cat("Input-oriented efficiencies:\n")
print(round(result_in$effvalsq, 3))

## ----output-orientation-------------------------------------------------------
result_out <- qDEA(X = X, Y = Y, orient = "out", RTS = "CRS", qout = 0.10)
cat("Output-oriented efficiencies:\n")
print(round(result_out$effvalsq, 3))

## ----graph-orientation--------------------------------------------------------
result_graph <- qDEA(X = X, Y = Y, orient = "inout", RTS = "VRS", qout = 0.10)
cat("Graph-oriented distances:\n")
print(round(result_graph$distvalsq, 3))

## ----rts-comparison-----------------------------------------------------------
# Constant Returns to Scale (CRS)
result_crs <- qDEA(X = X, Y = Y, orient = "in", RTS = "CRS", qout = 0.10)

# Variable Returns to Scale (VRS)
result_vrs_comp <- qDEA(X = X, Y = Y, orient = "in", RTS = "VRS", qout = 0.10)

# Decreasing Returns to Scale (DRS)
result_drs <- qDEA(X = X, Y = Y, orient = "in", RTS = "DRS", qout = 0.10)

# Increasing Returns to Scale (IRS)
result_irs <- qDEA(X = X, Y = Y, orient = "in", RTS = "IRS", qout = 0.10)

# Compare results
rts_comparison <- data.frame(
  Hospital = CST22$HOSPITAL,
  CRS = round(result_crs$effvalsq, 3),
  VRS = round(result_vrs_comp$effvalsq, 3),
  DRS = round(result_drs$effvalsq, 3),
  IRS = round(result_irs$effvalsq, 3)
)

print(rts_comparison)

## ----bootstrap, eval=FALSE----------------------------------------------------
# # Run qDEA with bootstrap (100 replications)
# # Note: This takes longer to run
# result_boot <- qDEA(X = X, Y = Y,
#                     orient = "in",
#                     RTS = "VRS",
#                     qout = 0.10,
#                     nboot = 100,
#                     seedval = 12345)
# 
# # Access bias-corrected estimates
# bias_corrected <- result_boot$BOOT_DATA$effvalsq.bc
# 
# # Compare original and bias-corrected
# comparison_boot <- data.frame(
#   Hospital = CST22$HOSPITAL,
#   Original = round(result_boot$effvalsq, 3),
#   Bias_Corrected = round(bias_corrected, 3),
#   Bias = round(result_boot$effvalsq - bias_corrected, 3)
# )
# 
# print(comparison_boot)

## ----peers--------------------------------------------------------------------
# Run qDEA to get peer information
data(CST11)
X <- as.matrix(CST11$EMPLOYEES)
Y <- as.matrix(CST11$SALES_EJOR)

result_peers <- qDEA(X = X, Y = Y,
                     orient = "out",
                     RTS = "VRS",
                     qout = 0.10)

# Access peer information
peers <- result_peers$PEER_DATA$PEERSq
print(peers)

## ----projections--------------------------------------------------------------
# Run with projection calculation
result_proj <- qDEA(X = X, Y = Y,
                    orient = "out",
                    RTS = "VRS",
                    qout = 0.10,
                    getproject = TRUE)

# Access projected values
X_target <- result_proj$PROJ_DATA$X0HATq
Y_target <- result_proj$PROJ_DATA$Y0HATq

# Compare actual vs. target
projection_comparison <- data.frame(
  Store = CST11$STORE,
  Actual_Input = X[,1],
  Target_Input = round(X_target[,1], 1),
  Actual_Output = Y[,1],
  Target_Output = round(Y_target[,1], 1),
  Efficiency = round(result_proj$effvalsq, 3)
)

print(projection_comparison)

## ----iterative----------------------------------------------------------------
# Run with multiple iterations
result_iter <- qDEA(X = X, Y = Y,
                    orient = "out",
                    RTS = "CRS",
                    qout = 0.15,
                    nqiter = 5)  # Allow up to 5 iterations

# Check number of iterations used
cat("Iterations completed:", result_iter$LP_DATA$qiter, "\n")

# Check actual proportion of outliers
cat("Actual outlier proportion:", round(result_iter$LP_DATA$qhat, 4), "\n")

## ----output-structure, eval=FALSE---------------------------------------------
# result <- qDEA(X = X, Y = Y, orient = "out", RTS = "CRS", qout = 0.10)
# 
# # Main efficiency scores
# result$effvals    # DEA efficiency scores
# result$effvalsq   # qDEA efficiency scores
# result$distvals   # DEA distance measures
# result$distvalsq  # qDEA distance measures
# 
# # Input data (for reference)
# result$INPUT_DATA
# 
# # Bootstrap data (if nboot > 0)
# result$BOOT_DATA
# 
# # Peer information
# result$PEER_DATA
# 
# # Projected values (if getproject = TRUE)
# result$PROJ_DATA
# 
# # LP solver information
# result$LP_DATA

## ----complete-analysis--------------------------------------------------------
# Load data
data(CST22)

# Prepare data
X <- as.matrix(CST22[, c("DOCTORS", "NURSES")])
Y <- as.matrix(CST22[, c("OUT_PATIENTS", "IN_PATIENTS")])

# Run comprehensive analysis
result_complete <- qDEA(
  X = X, 
  Y = Y,
  orient = "in",
  RTS = "VRS",
  qout = 0.10,
  nqiter = 3,
  getproject = TRUE
)

# Create comprehensive results table
complete_results <- data.frame(
  Hospital = CST22$HOSPITAL,
  Doctors = CST22$DOCTORS,
  Nurses = CST22$NURSES,
  Out_Patients = CST22$OUT_PATIENTS,
  In_Patients = CST22$IN_PATIENTS,
  DEA_Eff = round(result_complete$effvals, 3),
  qDEA_Eff = round(result_complete$effvalsq, 3),
  Target_Doctors = round(result_complete$PROJ_DATA$X0HATq[,1], 1),
  Target_Nurses = round(result_complete$PROJ_DATA$X0HATq[,2], 1)
)

print(complete_results)

# Summary statistics
cat("\n=== Summary Statistics ===\n")
cat("Mean DEA efficiency:", round(mean(result_complete$effvals), 3), "\n")
cat("Mean qDEA efficiency:", round(mean(result_complete$effvalsq), 3), "\n")
cat("Efficient units (qDEA ≥ 0.99):", 
    sum(result_complete$effvalsq >= 0.99), "out of", nrow(X), "\n")

# Potential input reductions
input_savings <- cbind(
  X - result_complete$PROJ_DATA$X0HATq
)
colnames(input_savings) <- c("Doctor_Savings", "Nurse_Savings")

cat("\nTotal potential input reductions:\n")
cat("Doctors:", round(sum(input_savings[,1]), 1), "\n")
cat("Nurses:", round(sum(input_savings[,2]), 1), "\n")

## ----visualization, fig.width=7, fig.height=5---------------------------------
# Create efficiency comparison plot
data(CST22)
X <- as.matrix(CST22[, c("DOCTORS", "NURSES")])
Y <- as.matrix(CST22[, c("OUT_PATIENTS", "IN_PATIENTS")])

result_viz <- qDEA(X = X, Y = Y, orient = "in", RTS = "VRS", qout = 0.10)

# Prepare data for plotting
plot_data <- data.frame(
  Hospital = CST22$HOSPITAL,
  DEA = result_viz$effvals,
  qDEA = result_viz$effvalsq
)

# Simple bar plot
barplot(
  t(as.matrix(plot_data[, c("DEA", "qDEA")])),
  beside = TRUE,
  names.arg = plot_data$Hospital,
  col = c("skyblue", "coral"),
  main = "Hospital Efficiency: DEA vs qDEA",
  ylab = "Efficiency Score",
  xlab = "Hospital",
  legend.text = c("DEA", "qDEA"),
  args.legend = list(x = "topright")
)
abline(h = 1, lty = 2, col = "red")

## ----session-info-------------------------------------------------------------
sessionInfo()