12-Software-Comp.Rmd

---
output:
  pdf_document: default
  html_document: default
editor_options: 
  markdown: 
    wrap: sentence
---

```{r, echo=FALSE, eval=FALSE}
library(bookdown); library(rmarkdown); rmarkdown::render("12-Software-Comp.Rmd", "pdf_book")
```

# Benchmarking the Benchmarkers

```{r setupch12, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(kableExtra)
library(DiagrammeRsvg, quietly=TRUE)
library(DiagrammeR, quietly=TRUE)
library(rsvg, quietly=TRUE)
library(htmltools, quietly=TRUE)
library(TRA)
```

```{r, echo=FALSE, eval=FALSE}
library(bookdown); library(rmarkdown); rmarkdown::render("12-Software-Comp.Rmd", "pdf_book")
```

## Introduction

Over the years, a variety of software packages have been developed for doing Data Envelopment Analysis.
The \index{Benchmarking} Benchmarking package was one of the earliest and most full-featured packages available.
It continues to be the most downloaded R package for doing DEA to this day.

Let's review data on the current frequency of downloads for DEA packages.

```{r package_download_stats, echo=FALSE}
library("ggplot2")
library("dlstats")

dea_packages <- c("Benchmarking", "nonparaeff", "DJL", "rDEA", "deaR", "MultiplierDEA" )

pkg_dl_data <- cran_stats(packages = dea_packages)

```

```{r plot_dea_package_use, echo=FALSE}

ggplot(pkg_dl_data, aes(end, downloads, group=package, color=package)) +
    geom_line() + geom_point(aes(shape=package)) + 
    labs(title = "Monthly Downloads of DEA Packages",
       subtitle = "Wide choice of DEA packages",
       caption = "Source: CRAN", 
       x = "Month", y = "Downloads") 

```

The chapter "Computational Aspects of DEA" by Iqbal Ali in Data Envelopment Analysis, edited by Charnes, Cooper, Lewin, and Seiford does an interesting job with exploring computational issues in early DEA implementations.
Among the key insights is that a finite approximation for $\epsilon$ should not be used to maximize slacks.
Iqbal Ali tests $\epsilon =10^{-5}$ to $\epsilon=10^{-8}$ demonstrating how DEA implementations can return incorrect values or even unbounded solutions.

Any reliable DEA package should be tested to ensure that they do not rely on finite values for $\epsilon$ and instead implement a two-phase approach as discussed earlier.

Other issues considered by Iqbal Ali include:

-   Which DEA models are implemented?
-   Is an anticycling technique employed?
-   Testing against ill-conditioned data

Most R packages for doing DEA rely on a standard linear programming solver.

```{r loadhelperfiles, echo=FALSE }
source('Helperfiles.R')
#knitr::read_chunk('Helperfiles.R')
#<<poscolfunct>>   
   # This reads in a chunk that defines the poscol function
   # This function will filter out columns that are zero
   # More precisely, it factors out column with
   # column sums that are zero.  This is helpful
   # for tables of lambda values in DEA.
source('Helperfiles.R')
#knitr::read_chunk('Helperfiles.R')
#<<DrawIOdiagramfunction>>   
```

## Demonstrate Each Package

Create a common dataset, perhaps the ill-conditioned dataset from Iqbal Ali's paper to test $\epsilon$ and computational values.

## Testing with Ill-Conditioned Data

Data scaled with wide ranges may highlight problems in the DEA implementation of a package or a weakness in the underlying solver's numerical methods.
In the late 1990s I received a phone call from an Economics professor that had done her PhD doing DEA but was at a complete loss as to why her simple DEA input-oriented efficiency scores were as large as 47 when they should have ranged from just zero to 1.0.
She was was using SAS-OR.
I asked her the magnitudes of the inputs and outputs and found that they were in units of dollars but ranged to the millions or billions.
I suggested simply rescaling and that immediately solved the problem.
Another random call from a student team of Michigan State students had the same problem while using the Excel Solver.
I expect that the current linear programming solvers used by R DEA packages are robust enough to avoid scaling problems but is another item that should be checked.

As discussed earlier, a common problem in some earlier DEA implementations was the use of finite values of $\epsilon$ causing computational issues and solvers that were susceptible to poorly scaled data.

Testing of large scale DEA problems dates back to the work of Dr. Dick Barr and Dr. Matt Durchholz at Southern Methodist University in the late 1980s and early 1990s.
Since then Jose Dula and others have continued this work.
Since DEA requires a series of simple linear programs, in general they don't represent a major computational burden.
Any modern DEA package should be able to handle 1000 DMU problems with ease.
Problems with tens of thousands to millions or more could be considered but careful consideration should be given to the impact of using an outlier based approach such as DEA in such a large dataset.
A six-sigma type of extreme value for positive performance of a particular DMU could cause a tremendous distortion of the efficiency frontier rendering the scores for many of the other DMUs inappropriate.

For demonstrating DEA, let's use a data set from Charnes, Cooper, and Li (1990) that was later also used Iqbal Ali for his computational performance.
This data set consists of three inputs and three outputs to be examined with an input-oriented CCR (CRS) envelopment model.

For the sake of illustration purposes, we will only use a select portion of the dataset to to start with.

```{r create-bad-data, echo=FALSE}
x <- matrix(c(483.01, 1397736, 616961,
              371.95,  355509, 385453,
              268.23,  685584, 341941,
              202.02,  452713, 117424,
              197.93,  471650, 112634,
              178.96,  423124, 189743,
              148.04,  367012,  97004,
              184.93,  408311, 111904,
              123.33,  245542,  91861,
              116.91,  305316,  91710,
              129.62,  295812,  92409,
              106.26,  198703,  53499), ncol=3, byrow=TRUE,
           dimnames=list(LETTERS[1:12],c("Labor", "WF", "INV")))

y <- matrix(c(6785798, 1594957, 1088699, 
              2505984,  545140,  835745,
              2292025,  406947,  473600,
              1158016,  135939,  336165,
              1244124,  204909,  317709,
              1187130,  190178,  605037,
               658910,   86514,  239760,
               993238, 1411954,  353896,
               854188,  135327,  239360,
               606743,   78357,  208188,
               736545,  114365,  298112,
               454584,   67154,  233733), ncol=3, byrow=TRUE,
            dimnames=list(LETTERS[1:12],c("GIOF", "PT", "RS")))

kbl (cbind(y,x), caption="Complex Data from Charnes, Cooper, and Li (1990)", booktabs = T, digits = 4) |>
  kable_styling(latex_options = c("hold_position"))
```

```{r CCL_1990_IO_Model}
XFigNames <- c("Labor", "WF", "INV")
YFigNames <- c("GIOF", "PT", "RS")
Figure <- DrawIOdiagram (XFigNames,YFigNames, '"\nCCR\nIO\n "')  

tmp <- capture.output(rsvg_png(charToRaw (export_svg(Figure)),
                               'CCL_1990_IO_diagram.png'))
```

![Input-Output Model from Charnes, Cooper, Li, (1990)](CCL_1990_IO_diagram.png){#fig:CCL1990}

For each package, we will examine an input-oriented, constant returns to scale model.
We will select slack maximization whenever it is available.
Results are given for each analysis demonstrate the structure of the the lambda (envelopment variable) values returned as well as the efficiency scores.

Results are summarized at the end after all analyses to illustrate the efficiency score.

A few things to take note of:

-   `suppressPackageStartupMessages` is useful to cut down on warnings about duplicated functions among other items. After all, it isn't surprising that most of the packages have a function named `dea`.
-   Functions are called from specific packages using the `::` such as `Benchmarking::dea` to call the `dea` from the `Benchmarking` package.
-   Parameters are passed by name into the DEA function to avoid ambiguity and to clarify the different naming \index{Conventions!Naming} conventions used in each package.

## The Benchmarking Package

```{r run-benchmarking}
suppressPackageStartupMessages(library (Benchmarking))

res_benchmarking <- Benchmarking::dea (X=x, Y=y, RTS="crs", 
                                       ORIENTATION="in", 
                                       SLACK=TRUE, DUAL=FALSE)

kbl (cbind(res_benchmarking$eff,poscol(res_benchmarking$lambda)),
        caption="Results from the Benchmarking package", 
     booktabs = T, digits = 4) |>
  kable_styling(latex_options = 
                  c("hold_position"))
```

The results from Benchmarking are simple and straightforward.

-   DMU names are used for column names in the efficiency scores but in the lambda (envelopment variables), they are numbered as L1, L2, etc. and names are not used in rows.

```{r run-nonparaeff}

suppressPackageStartupMessages(library ("nonparaeff"))

nonparaeff_data <- data.frame(y,x)  # Format for data entry
res_nonparaeff <- nonparaeff::dea(base = nonparaeff_data,
                                  noutput = 3,        # count from left
                                  orientation = 1,    # 1 = IO
                                  rts = 1,            # 1 = CRS
                                  onlytheta = FALSE)

kbl (res_nonparaeff$eff, booktabs = T, digits = 4,
     caption="Efficiency Scores from nonparaeff")

res_lambda_nonparaeff<- res_nonparaeff[,1:13]

kbl (cbind(res_nonparaeff$eff,poscol(res_lambda_nonparaeff)),
     booktabs = T, digits = 4,
     caption="Results from the Benchmarking package") |>
  kable_styling(latex_options = 
                  c("hold_position", "scale_down"))
```

Some items to note about `nonparaeff`.

-   The input and outputs are passed together with outputs first, then inputs and the number of outputs specified.
-   The $\lambda$ variables are returned separately for each DMU.
-   Returned objects lose the row and column names.

## The rDEA Package

```{r run-rDEA}
suppressPackageStartupMessages(library (rDEA))

res_rDEA <- rDEA::dea (XREF=x, YREF=y, X=x, Y=y,
                               model="input", RTS="constant")

kbl (cbind(res_rDEA$thetaOpt, res_rDEA$lambda), booktabs = T, digits = 4, 
     caption="Results from Benchmarking") |>
kable_styling (latex_options = c("hold_position", "scale_down"))
```

Some notable characteristics of using the rDEA package are the following.

-   The rDEA allows for passing input prices in order to do a cost-minimization analysis.\
-   The model parameter combines both orientation and model decisions.
-   Explicitly allows for separating the production technology (XREF and YREF) from the units to be estimated.

## The DJL Package

The DJL package comes from the name "Distance Measure Based Judgement and Learning" but is coincidentally also the initials of the primary author, Dr. Dong-Joon Lim.
The emphasis in this package is on a variety efficiency models from the perspective of distance metrics.
With this it provides a variety of models for accommodating super-efficiency related issues - particularly infeasibility under certain situations.

```{r run-DJL}
suppressPackageStartupMessages(library (DJL))

res_DJL <- DJL::dm.dea (xdata=x, ydata=y, 
                        rts="crs", orientation="i")

kbl(cbind(res_DJL$eff, res_DJL$lambda), booktabs = T, digits = 4,
       caption="Efficiency Results from DJL")|>
kable_styling (latex_options = c("hold_position", "scale_down"))

```

## The `multiplierDEA` Package

The `multiplierDEA` package was developed by Aurobindh Kalathil Puthanpura.
It was developed based on an interest in the multiplier DEA model, applying weight restrictions in the multiplier model, and revisiting cross-efficiency.

```{r run-MultiplierDEA}

suppressPackageStartupMessages(library (MultiplierDEA))

res_MultiplierDEA <- MultiplierDEA::DeaMultiplierModel (x=x, y=y, 
                               rts="crs", orientation="input")

kbl(cbind (res_MultiplierDEA$Efficiency, res_MultiplierDEA$Lambda),
    booktabs = T, digits = 4,
       caption="Efficiency Resfults from multiplierDEA") |>
kable_styling (latex_options = c("hold_position", "scale_down"))

```

## The `deaR` Package

The \index{deaR package} `deaR` package has a rich collection of data sets included in the package.
It also includes a broader collection of DEA variants than most of the other packages including:

-   Malmquist Productivity Index,
-   Undesirable input and output models,
-   Fuzzy DEA,
-   Cross-Efficiency,
-   Bootstrapping DEA.

Alas, the data structure is quite different from the other DEA packages making direct comparisons somewhat more difficult.
A comparison of deaR will be added in the future.

```{r run-deaR}

suppressPackageStartupMessages(library (deaR))

#res_deaR <- deaR::model_basic(orientation="io", rts="crs", 
#                              maxslack=TRUE, 
# Note:  Data format is complicated...

#kable(res_deaR$eff, caption="Results from deaR")
#kable(res_deaR$lambda, caption="Results from deaR")

```

## Summarizing the Results from Different DEA Packages

Let's now pull all the results together.
Ideally, all of the efficiency scores should match.

First, let's review the varying commands to illustrate the different syntax for the calls.

```{r DEA-commands, eval=FALSE}

Benchmarking::dea (X=x, Y=y, RTS="crs", ORIENTATION="in", 
                                       SLACK=TRUE, DUAL=FALSE)

DJL::dm.dea (xdata=x, ydata=y, rts="crs", orientation="i")

MultiplierDEA::DeaMultiplierModel (x=x, y=y, 
                               rts="crs", orientation="input")

nonparaeff::dea(base = nonparaeff_data,noutput = 3, # count from left
                orientation = 1,    # 1 = IO
                rts = 1,            # 1 = CRS
                onlytheta = FALSE)

rDEA::dea (XREF=x, YREF=y, X=x, Y=y,model="input", RTS="constant")

```

```{r combined-results, echo=FALSE}

combined_eff <- cbind (res_benchmarking$eff, res_nonparaeff$eff,
                       res_rDEA$thetaOpt, res_DJL$eff, 
                       res_MultiplierDEA$Efficiency)

colnames (combined_eff) <- c("Benchmarking", "nonparaeff", "rDEA", 
                             "DJL", "MultiplierDEA")

kbl (combined_eff, booktabs = T, digits = 6, 
     caption = "Comparison of Efficiency Scores")      |>
  kable_styling(latex_options=c("hold_position","scale_down"))
```

```{r combined-results-using-kable}
combined_eff2 <- cbind (res_nonparaeff$eff-res_benchmarking$eff,
                       res_rDEA$thetaOpt-res_benchmarking$eff, 
                       res_DJL$eff-res_benchmarking$eff, 
                       res_MultiplierDEA$Efficiency-res_benchmarking$eff)

colnames (combined_eff2) <- c("nonparaeff", "rDEA", 
                             "DJL", "MultiplierDEA")

kbl (combined_eff2, digits=16, booktabs=T,
     caption="Difference from Benchmarking Package Results") |>
  kable_styling(latex_options=c("hold_position"))
```

The result is that the largest magnitude differences between Benchmarking's results and that of the other packages tested was a positive `r max(combined_eff2)` and the most negative (smaller than Benchmarking's efficiency scores) was `r min(combined_eff2)`.
These values are quite well behaved and more precise than most data used for analyses.

```{r, echo=FALSE}
comparative_deviations <- rbind (colMeans(combined_eff2),
                      colMeans(abs(combined_eff2)))

rownames(comparative_deviations) <- 
  c("Mean Difference", "Mean Absolute Deviation")

kbl (summary(combined_eff2),  booktabs=T,
     caption="Summary of Differences between Benchmarking and other Packages") |>
  kable_styling(latex_options="hold_position")
```

Computational performance could be further compared in many ways such as the following.

-   Additional comparisons could be made for the slacks to ensure that non-radial slacks are maximized.
-   Larger data sets could be tested.
-   Other forms of ill-conditioned data could be tested.
-   Computational speed could be tested with data sets of over 1000.
-   Returns to scale other than constant returns to scale or CRS could be used.
-   Detailed comparisons of lambda variables could be examined.
-   Testing could be done to to see if small numerical anomalies arise and affect results or if rounding processes to address this cause other problems. For example, LP solvers may return lambda values such as $-10^{-16}$ rather than zero or $0.99999998$ rather than $1.0$.
-   Examining DEA bootstrapping.

## DEA Package Dependencies

Let's explore another view of DEA packages by examining the packages that they depend upon.
A key aspect is what optimization engine each package uses.

```{r DEA-package-dependencies, eval=FALSE}

library("miniCRAN")

dea_packages_cran <- c("Benchmarking", "nonparaeff", "DJL", "rDEA") 
#solver_packages <- c("lpsolve", "lpsolveAPI", "GLPK", "Rglpk",
#                      "glpkAPI", "ompr", "ROI","rsymphony")
# tags <- c("ggplot2", "data.table", "plyr", "knitr", "shiny", 
#           "xts", "lattice")
pkgDep(dea_packages_cran , suggests = FALSE, 
       enhances = FALSE)
        #, availPkgs = cranJuly2014

dg <- makeDepGraph(dea_packages_cran, enhances = FALSE)
set.seed(1)
plot(dg, legendPosition = c(-1, -1), vertex.size = 10, cex = 0.7)
```

This graph seems to indicate that Benchmarking is relatively independent of other packages despite being a powerful and well accepted tool.
Other packages appear to rely on a complex web of packages unrelated to DEA.
Interestingly, it appears that \index{lpSolveAPI} `lpSolveAPI` is only used by `DJL` and `Benchmarking` while the other packages depend on other solvers.

It would be more interesting to filter this to only include solver engines such as 
\index{lpSolveAPI} `lpSolveAPI`, \index{lpSolve} `lpSolve`, \index{glpk} `glpk`, 
and \index{symphony} `symphony` along with the DEA packages.

## Summary of Package Features

Comparative features Columns are used for each package Rows are:

-   Package Name
-   Package Author(s)
-   Version and date
-   Main LP package used? (ex. LPSolve API, GLPK, etc.)
-   All Traditional Models (IO/OO, CRS/VRS/DRS/IRS)?
-   Both Envelopment & Multiplier?
-   Slack Maximization?
-   Malmquist?
-   Weight Restrictions?
-   Bad outputs/Good Inputs?
-   Non-discretionary Inputs/Outputs?
-   Directional Distance Functions?
-   Window Analysis?
-   Cross-Efficiency?
-   Special Features?

## Future Work