Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
128,074 characters · 30 sections · 129 citation commands
Conventional and Fuzzy Data Envelopment Analysis with deaR
Data envelopment analysis (DEA) Charnes1978 is a non-parametric technique used to measure the relative efficiency of a homogeneous set of decision making units (DMUs) that use multiple inputs to obtain multiple outputs. Using mathematical programming methods, DEA allows the identification of the best practice frontier (efficient frontier). DMUs that form the best practice frontier are qualified as efficient, while DMUs that move away from the frontier are inefficient.
DEA has been used to evaluate efficiency in many different fields such as education, agriculture and farm, banking, health, transportation, public administration, etc. Recently, Emrouznejad2018 have compiled a list of more than $10000$ articles related to DEA. Along with the methodological advances and development of new DEA models and practical applications, software (commercial and non-commercial) has also appeared to facilitate its use by both researchers and practitioners. Daraio2019 perform a review of DEA software available to assess efficiency and productivity and compare their different options. Among the commercial software, MaxDEA (\url{http://maxdea.com/MaxDEA.htm}), DEA-Solver-PRO Cooper2007 and DEAFrontier Zhu2014 stand out for the diversity of DEA models implemented. Also, some of the most popular software, such as SAS, GAMS or Stata, include modules to estimate efficiency using basic DEA models. Moreover, Alvarez2020 have written a package for MATLAB that includes several DEA models like radial DEA, directional distance function, Malmquist index and cross-efficiency among others.
As far as non-commercial software is concerned, one of the most widely used is undoubtedly DEAP, developed by Coelli1996. However, this software is very limited in terms of the variety of DEA models it can solve: basic models by Charnes1978 and Banker1984 (referred to as CCR and BCC models), Malmquist index (FGNZ decomposition due to Fare1994) and cost/revenue models.
Nowadays, the use of non-commercial software to apply DEA models involves the use of R packages. The first R packages for measuring efficiency and productivity using DEA were FEAR by Wilson2008 (it is distributed under license and has recently been updated after many years without doing so) and Benchmarking by Bogetoft2022, which complements the book by the same authors Bogetoft2011. These two packages contain methods that allow the use of different technology assumptions and different efficiency measures (radial measures, super-efficiency, additive models, cost efficiency, etc.), but it is worth noting that they include routines to apply the bootstrap methods described by Simar1998. Without being exhaustive, other relevant packages available in R for analysing efficiency and productivity using the DEA technique are:
Unlike the R packages cited above, deaR Coll2022 has a wide variety of models implemented that allow the user to apply both conventional and fuzzy models, which consider imprecision or uncertainty in the data. Moreover, the package includes super-efficiency, cross-efficiency, Malmquist index and bootstrapping models. The functions in deaR are designed to be easily used by a non-expert R user.
This paper is organized as follows. In Section (ref) we describe how to use the deaR package in 3 steps: introducing data, running a model and extracting the results. We also present the main S3 data classes and review some basic radial efficiency models. In Section (ref) we review other conventional efficiency models such as multiplier, free disposal hull, directional, non-radial, additive or SBM models. In Section (ref) we introduce some special features on data, namely non-controllable, non-discretionary and undesirable variables. Section (ref) is devoted to the main super-efficiency models: radial, SBM and additive models. In Section (ref) we present the cross-efficiency analysis, while in Section (ref) we review some popular fuzzy models: Kao-Liu, Guo-Tanaka and possibilistic models. We show the Malmquist productivity indices in Section (ref), and the bootstrapping methodology in Section (ref). Finally, we give an example on how to perform a non-parametric metafrontier analysis in Section (ref) and provide some conclusions in Section (ref).
Throughout this paper, we consider $\mathcal{D}=\left\{ \textrm{DMU}_1, \ldots ,\textrm{DMU}_n \right\} $ a set of $n$ DMUs with $m$ inputs and $s$ outputs. Matrices $X=(x_{ij})$ and $Y=(y_{rj})$ are the input and output data matrices, respectively, where $x_{ij}$ and $y_{rj}$ denote the $i$-th input and $r$-th output of the $j$-th DMU. We also assume that $x_{ij}$ and $y_{rj}$ are all positive, i.e., greater than $0$. Nevertheless, this restriction is not strictly considered in deaR package, and we can run models with negative or zero data, but results should be interpreted with care. In general, we denote vectors by bold-face letters and they are considered as column vectors unless otherwise stated. The elements of a vector are denoted by the same letter as the vector, but unbolded and with subscripts. The $0$-vector is denoted by $\bm{0}$ and the context determines its dimension.
First of all, deaR is a package available on CRAN and it can be installed with
In order to perform a DEA study over DMUs, we need to solve optimization problems. Specifically, all the models implemented in deaR involve linear programming problems. There are many options to solve linear programs using R, as shown in the Optimization Task View (linear programming section) but, in deaR, we use the lpSolve package Berkelaar2020. In fact, deaR can be considered as a wrapper of lpSolve.
The workflow in deaR has the following steps (see Figure (ref)):
Therefore, performing a DEA analysis with deaR can be divided into $3$ phases: introducing data, running a model and extracting the results. We will describe each phase in detail below.
Datasets in DEA are usually formatted in a spreadsheet table, in which one of the columns (typically the first one) corresponds to the DMUs identification labels, the next $m$ columns correspond to the different inputs, and the following $s$ columns correspond to the outputs. This type of data arrangement will be referred to as standard DEA dataset.
The package contains more than 20 datasets from different books and research papers, that can be used to reproduce the results obtained in those sources (use data(package = "deaR") to get the list of all available datasets). For example, the Fortune500 dataset is a standard DEA dataset containing 15 companies from the 1995's Fortune 500 list Zhu2014:
In this standard DEA dataset, the first column (“Company”) contains the names of the different firms included. The next three columns (“Assets”, “Equity” and “Employees”) define the inputs, and the last two columns (“Revenue” and “Profit”) correspond to the outputs.
The package deaR first needs to read the dataset and identify the DMUs, inputs, outputs and all other relevant information about the data. To do so, we will make use of the function make_deadata. For a standard DEA dataset, the syntax is straightforward. By default, the function assumes that the names of the DMUs are in the first column (although this can be changed with the parameter dmus) and therefore, we will only need to declare the number of inputs (ni) and the number of outputs (no):
Instead of defining the number of inputs/outputs with the parameters ni/no, we can explicitly define which columns contain the inputs and which the outputs, either by a number (position of the column) or by a column name. For instance, we can define a smaller DEA dataset with two inputs and one output without having to subset the dataframe:
In this example, the DEA dataset have “Equity” and “Employees” as inputs (columns 3 and 4) and “Profit” as output (defined by the column name). Alternatively, we can explicitly define the input and output data matrices, with DMUs in columns:
In this case, we only consider “Equity”, “Employees” and “Profit” of the first three DMUs. Moreover, if names are not provided by dimnames, then they are automatically generated as “DMU1”, “DMU2”, “Input1”, “Input2”, etc.
It is important to note that negative or zero data are not allowed by some models. To solve this problem, it is recommended to translate the base point of the inputs/outputs with negative or zero data in order to get only positive values. Nevertheless, depending on the nature of the data and the model, this may not be appropriate in some cases. Moreover, if there are data with very different orders of magnitude, then it is also recommended to redefine the units of measure in order to prevent ill-posed linear problems.
Finally, the resulting value of the make_deadata function is an object of class deadata which is a list with the following fields: \dirtree{ .1 deadata (for a model with $n$ DMUs, $m$ inputs and $s$ outputs.). .2 input: input data matrix of size $m\times n$.. .2 output: output data matrix of size $s\times n$.. .2 dmunames: character vector containing the names of the DMUs.. .2 nc_inputs: \textrm{integer vector identifying the non-controllable inputs (or \texttt{NULL})}.. .2 nc_outputs: \textrm{integer vector identifying the non-controllable outputs (or \texttt{NULL})}.. .2 nd_inputs: \textrm{integer vector identifying the non-discretionary inputs (or \texttt{NULL})}.. .2 nd_outputs: \textrm{integer vector identifying the non-discretionary outputs (or \texttt{NULL})}.. .2 ud_inputs: \textrm{integer vector identifying the undesirable inputs (or \texttt{NULL})}.. .2 ud_outputs: \textrm{integer vector identifying the undesirable outputs (or \texttt{NULL})}.. } Special features like non-controllable, non-discretionary and undesirable inputs/outputs are explained in Section (ref).
In general, a DMU is efficient if there is not any feasible activity in a given production possibility set “better than” the DMU, in the sense that consumes less inputs and produces more outputs. Hence, an efficiency model first establishes the production possibility set and then it checks if the DMU is efficient or not. If the DMU is inefficient, then the model usually gives a score and a target for improving the activity.
Once the data has been read and we have an object of class deadata, we can proceed to select and run a model. The package deaR has quite a wide range of different models available. Table (ref) lists all the model functions included in deaR. In this section we shall illustrate the use of those functions with classical examples of the basic radial CCR and BCC models contained in function model_basic.
The first model we are going to introduce is the so called CCR model Charnes1978,Charnes1979,Charnes1981. This model assumes that the production possibility set, i.e., the set of feasible activities defined by the set $\mathcal{D}$ of DMUs, is under constant returns to scale (CRS) and given by
where $X$ and $Y$ are the input and output data matrices, respectively, and $\bm{\lambda}=(\lambda_1,\ldots,\lambda_n)^\top$ is a column vector. Hence, $\textrm{DMU}_o\in\mathcal{D}$ is efficient if and only if there is no $(\mathbf{x},\mathbf{y})\in P$ such that $x_{io}\geq x_i$ and $y_{ro}\leq y_r$ with at least one strict inequality.
The CCR model can be either input or output-oriented. In the former case (see (ref) (a)) we look for determining the maximal proportionate reduction of inputs allowed by the production possibility set, while maintaining the current output level of $\textrm{DMU}_o$. On the other hand, in the output-oriented case (see (ref) (b)), we want to find the maximal proportionate increase of outputs while keeping the current input consumption of $\textrm{DMU}_o$:
where $\mathbf{x}_o=(x_{1o},\ldots,x_{mo})^\top$ and $\mathbf{y}_o=(y_{1o},\ldots,y_{so})^\top$ are column vectors. In a second stage, with our knowledge of the optimal objectives $\theta ^*$ or $\eta ^*$, we solve the following linear program, (ref) (a) or (ref) (b), in order to find the max-slack solution:
where the weights $\mathbf{w}^-$ and $\mathbf{w}^+$ are positive row vectors. In the input-oriented CCR model, $\textrm{DMU}_o$ is efficient if and only if $\theta ^*=1$ and $\omega ^*=0$. If $\textrm{DMU}_o$ is inefficient, then $0<\theta ^*\leq 1$ is the efficiency score and $(X\bm{\lambda}^*,Y\bm{\lambda}^*)$ is the target, that can be interpreted as the projection of $\textrm{DMU}_o$ onto the efficient frontier. Note that there can be inefficient DMUs with $\theta ^*=1$, called weakly efficient. On the other hand, in the output-oriented CCR model, we have $\eta ^*=1/\theta ^*$.
The BCC model Banker1984 considers the production possibility set under variable returns to scale (VRS),
where $\mathbf{e}=(1,\ldots ,1)$ is a row vector. Oriented BCC models are constructed by adding $\mathbf{e}\bm{\lambda} =1$ to the constraints of (ref) and (ref). Efficiency scores and targets are defined analogously to the CCR model.
In the production possibility set, the returns to scale condition can be changed to non-increasing (NIRS), non-decreasing (NDRS) or generalized returns to scale (GRS). In these cases, the condition $\mathbf{e}\bm{\lambda}=1$ is replaced by $0\leq \mathbf{e}\bm{\lambda}\leq 1$ (NIRS), $\mathbf{e}\bm{\lambda}\geq 1$ (NDRS) or $L\leq \mathbf{e}\bm{\lambda}\leq U$ (GRS), with $0\leq L\leq 1$ and $U\geq 1$. These conditions are added to the constraints of (ref) and (ref) in order to build models with different returns to scale.
The syntax of model_basic is very flexible, and contains a great deal of parameters allowing the user to run different models from within the same function. The main parameters are:
Other optional parameters of interest are:
Now, we can run a model for the Fortune500 dataset, which was already defined in the object dataFortune of class deadata. For instance, for the input-oriented CCR model, we can use
while for the BCC model of the same characteristics, we would write
although orientation = "io" and rts = "crs" are not necessary because they are the default values. Moreover, note that “datadea =” is not necessary in the first field because datadea is always the first parameter in model functions.
The results delivered by any model_xxx function is an object of class dea which is basically a list containing the information regarding the data, the call to function model_xxx and the results obtained for each DMU: \dirtree{ .1 dea. .2 modelname: name of the model.. .2 orientation: orientation of the model.. .2 rts: returns to scale of the model.. .2 DMU: \textrm{results of the model for each evaluated DMU}.. .2 data: \textrm{the object of class \texttt{deadata} to which the model has been applied}.. .2 dmu_eval: \textrm{evaluated DMUs}.. .2 dmu_ref: \textrm{evaluation reference set (with respect to which DMUs we have evaluated)}.. .2 maxslack: \textrm{logical parameter indicating if the max-slack solution has been computed}.. .2 weight_slack_i: \textrm{weight vector for input slacks in the max-slack solution}.. .2 weight_slack_o: \textrm{weight vector for output slacks in the max-slack solution}.. } Other specific parameters for some models (such as parameters \texttt{L} and \texttt{U} for generalized returns to scale, translation vectors \texttt{vtrans_i} and \texttt{vtrans_o} for undesirable variables, or \texttt{orientation_param}, with the input and output directions in directional models) are also stored because we want the class \texttt{dea} object to contain all the information about the model in question so that the results can be replicated.
The field DMU is itself a list containing, for each one of the evaluated DMUs, all the results obtained by the model. Namely,
In order to easily obtain those aforementioned results, there are several functions designed to extract the required information from the dea class object: efficiencies, lambdas, slacks, targets and multipliers. For example, in the case of the CCR model applied to the \texttt{Fortune500} dataset,
shows the efficiency scores stored in ccrFortune. It is important to remark that function efficiencies returns the scores (i.e., optimal objective values) of the model, that may not always be interpreted as efficiency scores.
Moreover, there are some other functions such as references or eff_dmus, whose parameter is also a dea object. Function references returns a list with the reference set for each inefficient DMU. Note that the reference set of a DMU is formed by all the efficient DMUs that appear in the linear combination that conforms its target. On the other hand, function \texttt{eff_dmus} returns an array with the efficient DMUs evaluated by the corresponding model.
Alternatively, instead of running all the aforementioned functions to extract the results, we can use the function summary.dea, which is a specific method for dea class objects and can be invoked with the generic function summary. For instance,
returns all the results stored in ccrFortune as a list of data frames. Otherwise, if returnList = FALSE (by default), then all these data frames are column-wise merged into a single data frame. Anyway, since exportExcel = TRUE (by default), the results in res are also exported to an Excel file named “ResultsDEAYYYYmmdd_HH.MM.SS.xls”, where \texttt{YYYYmmdd_HH.MM.SS} represents the current system date and time. This default name can be changed with the parameter \texttt{filename}.
Finally, we can make use of the function plot, which needs a dea class object in order to make some plots, depending on the model. For example, plot(ccrFortune) returns the plots shown in Figure (ref).
In this section, we are going to review the rest of the DEA models implemented in functions of the form model_xxx. Analogously to model_basic, most of these functions also use as input parameters datadea, orientation, rts, dmu_eval, \texttt{dmu_ref}, \texttt{maxslack}, \texttt{weight_slack_i}, \texttt{weight_slack_o} and \texttt{returnlp}.
All basic radial models can be applied in multiplier form Charnes1962, solving the dual of the linear problem in the first stage. For example, the input and output-oriented GRS models in multiplier form are given by
respectively, where $\mathbf{e}=(1,\ldots,1)$ is a row vector, $\mathbf{v}$ and $\mathbf{u}$ are row vectors interpreted as input and output “weights” respectively, and $\xi _L$, $\xi _U$ are the multipliers associated to the returns to scale constraints. Models with different returns to scale can be deduced from (ref) by taking $\xi _L=\xi _U=0$ (CRS), $L=U=1$ (VRS), $\xi _L=0$, $U=1$ (NIRS) or $L=1$, $\xi _U=0$ (NDRS). In any case, $\textrm{DMU}_o$ is efficient if and only if the optimal objective of (ref) (a) or (b) is equal to $1$ and there exists at least one optimal solution with positive optimal weights $\mathbf{v}^*>\mathbf{0}$, $\mathbf{u}^*>\mathbf{0}$. For this reason, the zeroes in the non-negativity conditions of $\mathbf{v}$ and $\mathbf{u}$ are usually replaced by a positive non-Archimedean infinitesimal $\epsilon $.
Multiplier models are applied using model_multiplier. For example,
solves the CRS input-oriented multiplier model applied to the Fortune500 dataset, stores the multipliers in a list named multFortune and shows the efficiency scores. Note that rts = "crs" and orientation = "io" are not necessary because they are the default values. Parameter epsilon is the non-Archimedean infinitesimal $\epsilon$ whose default value is $0$. It is important to remark that a too high positive value for epsilon can significantly alter the results and produce infeasibilities. Precisely, in our example, “General Motors” and “Ford Motor” have \texttt{NA} efficiency scores because \texttt{epsilon} is too high and the model becomes infeasible. On the other hand, we can also use \texttt{model_basic} with parameter \texttt{compute_multiplier = TRUE}, but in this case we can not set parameter \texttt{epsilon}, which is taken as $0$.
Free disposal hull (FDH) models consider the production possibility set
where $\mathbf{x}_j=(x_{1j},\ldots,x_{mj})^\top$ and $\mathbf{y}_j=(y_{1j},\ldots,y_{sj})^\top$ are column vectors. In fact, it is equivalent to consider VRS with $\lambda_j\in \left\{ 0,1\right\} $ binary variables Thrall1999, Cherchye2000, Deprins2006. In deaR, these models are applied using model_fdh, considering the production possibility set (ref) into the models implemented in model_basic, including directional models (see Section (ref)). For example, we can replicate some results in Mamizadeh2013 for the input-oriented case:
Note that the efficiency score of DMU sc16 is not $1$, contrary to the results shown in Mamizadeh2013.
The orientation can also be generalized, leading to directional models as described in Chambers1996,Chambers1998. The associated linear program for CRS and its second stage are given by
respectively, where $\mathbf{g}=(-\mathbf{g}^-,\mathbf{g}^+)\neq \mathbf{0}$ is a preassigned direction (with $\mathbf{g}^-\in \mathbb{R}^m$ and $\mathbf{g}^+\in\mathbb{R}^s$ column vectors), while the weights $\mathbf{w}^-$ and $\mathbf{w}^+$ are positive row vectors. Different returns to scale can be easily considered by adding the corresponding constraints: $\mathbf{e}\bm{\lambda}=1$ (VRS), $0\leq \mathbf{e}\bm{\lambda}\leq 1$ (NIRS), $\mathbf{e}\bm{\lambda}\geq 1$ (NDRS) or $L\leq \mathbf{e}\bm{\lambda}\leq U$ (GRS), with $0\leq L\leq 1$ and $U\geq 1$. Efficient DMUs are those with optimal objectives $\beta ^*=0$ and $\omega ^*=0$. For inefficient DMUs, targets are also given by $(X\bm{\lambda}^*,Y\bm{\lambda}^*)$. If $\mathbf{g}=(-\mathbf{x}_o,\mathbf{0})$ then the model is input-oriented and $\beta ^*=1-\theta ^*$. On the other hand, if $\mathbf{g}=(\mathbf{0},\mathbf{y}_o)$ then the model is output-oriented and $\beta ^*=\eta ^*-1$. Moreover, if $\mathbf{g}=(-\mathbf{x}_o,\mathbf{y}_o)$ then the model is non-oriented and $\beta ^*$ coincides with the generalized Farrell measure Briec1997.
In function model_basic, directional models are selected by setting parameter orienta\-tion = "dir". Then, we can set the input and output directions $\mathbf{g}^-,\mathbf{g}^+$ of each DMU by means of parameters dir_input and dir_output respectively which, in general, are matrices of size [number of inputs/outputs]$\times$[number of DMUs in dmu_eval]. If dir_input or \texttt{dir_output} are omitted, then they are assumed to be the input or output matrices (of DMUs in \texttt{dmu_eval}), respectively. Moreover, if \texttt{dir_input} or \texttt{dir_output} are vectors of length [number of inputs/outputs] then the same directions are applied to all evaluated DMUs. Finally, if \texttt{dir_input} or \texttt{dir_output} are scalars, then the same constant directions are applied to all inputs/outputs and all evaluated DMUs.
In the following example, considering the Fortune500 dataset, we compute the generalized Farrell measures and later, we apply an input-oriented directional model that fully contracts the first input, contracts the second input by half of the usual, and does not contract the third input, all for CRS:
\paragraph{Range directional models.} Other particular cases of directional models are the range directional models (RDM), in which $g^-_i=x_{io}-\min\left\{ x_{i1},\ldots ,x_{in}\right\} $ for $i=1,\ldots ,m$, and $g^+_r=\max\left\{ y_{r1},\ldots ,y_{rn}\right\} -y_{ro}$ for $r=1,\ldots ,s$, under VRS Portela2004. These models are designed to deal with negative data and they are applied using model_rdm. Input and output-oriented versions are also considered by taking $\mathbf{g}^+=\mathbf{0}$ and $\mathbf{g}^-=\mathbf{0}$, respectively. In these cases, parameter orientation must be equal to "io" or "oo", respectively, instead of the default value "no". For example,
returns the optimal $\beta ^*$ scores. Moreover, the inverse range directional models (IRDM) are constructed by substituting the non-zero components of the RDM directions by their inverses. They can be applied by setting parameter irdm = TRUE.
The non-radial models defined by Fare1978 allow non-proportional reductions/augmentations in inputs/outputs. The input and output-oriented CRS versions are given by
respectively, where $\mathbf{e}=(1,\ldots,1)$ is a row vector of the adequate dimension, $\bm{\theta}=(\theta_1,\ldots ,\theta _m)^\top $ and $\bm{\eta}=(\eta_1,\ldots ,\eta _s)^\top $ are column vectors, while $\textrm{diag}(\bm{\theta})$ and $\textrm{diag}(\bm{\eta})$ are diagonal matrices. Different returns to scale can be easily considered by adding the corresponding constraints. A second stage for the input and output-oriented CRS versions are needed in order to find the max-slack solution, respectively:
where the weights $\mathbf{w}^-$ and $\mathbf{w}^+$ are positive row vectors. $\textrm{DMU}_o$ is efficient if and only if the optimal objectives $\bar{\theta}^*=1$ (or $\bar{\eta}^*=1$) and $\omega ^{+*}=0$ (or $\omega ^{-*}=0$). These non-radial models are applied using model_nonradial. For the max-slack solution, the weights $\mathbf{w}^-$ or $\mathbf{w}^+$ are introduced by means of parameter weight_slack, that can be a value, a vector of length [number of inputs/outputs], or a matrix of size [number of inputs/outputs]$\times$[number of DMUs in dmu_eval]. By default, these weights are set to $1$. For example, we can replicate the results in Wu2011:
Note that after reading the data, a warning message appears because there are data with very different orders of magnitude. Nevertheless, the model is correctly applied and, in this case, there is no need to rescale the data.
\paragraph{DEA/preference structure models.} Non-radial models are generalized into DEA/pref\-er\-ence structure models Zhu1996, that replace the objectives in (ref) by a weighted sum and remove the constraints $\bm{\theta}\leq \mathbf{1}$ or $\bm{\eta}\geq \mathbf{1}$. These models are applied using model_deaps. The weights of the objective, called preference weights, are introduced by means of the parameter weight_eff. Moreover, if the logical parameter restricted_eff is TRUE (by default) then constraints $\bm{\theta}\leq \mathbf{1}$ or $\bm{\eta}\geq \mathbf{1}$ are not removed. For example, we can apply the input-oriented VRS version to the Fortune500 dataset, with objective function $\frac{1}{6}\left( \theta _1+2\theta _2+3\theta _3\right) $:
Additive models Charnes1985 do not distinguish between orientations and do not need a second stage. The CRS version is given by
where the weights $\mathbf{w}^-$ and $\mathbf{w}^+$ are positive row vectors. Different returns to scale can be easily considered by adding the corresponding constraints. Hence, $\textrm{DMU}_o$ is efficient if and only if the optimal objective $\omega ^*=0$. Although weights must be positive and orientations are not considered, if $\mathbf{w}^+=\mathbf{0}$ then the model can be interpreted as input-oriented and if $\mathbf{w}^-=\mathbf{0}$ then it can be interpreted as output-oriented. But you have to take into account that if the weight of a slack is zero, then this slack is not taken into account by the objective function and hence, inefficient (weakly efficient) DMUs can get $\omega ^*=0$.
Additive models are applied using model_additive. The weights are introduced by parameters weight_slack_i and weight_slack_o. Moreover, parameter orientation can be either NULL (by default), "io" ($\mathbf{w}^+=\mathbf{0}$) or \texttt{"oo"} ($\mathbf{w}^-=\mathbf{0}$).
It is important to note that (ref) is not unit-invariant in general. Nevertheless, there are particular cases that are unit-invariant, like the measure of inefficiency proportions (MIP) model Cooper1999. This model takes the weights $w_i^-=1/x_{io}$ and $w^+_r=1/y_{ro}$ under VRS. For example,
for the Fortune500 dataset. Another important particular case is the range adjusted measure (RAM) of inefficiencies model Cooper1999,Cooper2001, that can be solved for the Fortune500 dataset using the following script:
Another family of additive models, called additive-Min, are developed by Aparicio2007 in order to find the closest targets to the efficient frontier. The CRS version is given by:
Different returns to scale can be considered easily adding the corresponding constraints. However, independently of the returns to scale, these models can produce non-monotonic scores.
Additive-Min models are applied using model_addmin. Program (ref) can be solved by the “MILP” method proposed by Aparicio2007 or the ”maximal friends” method proposed by Tone2010. We can choose the method by means of the parameter method, that can be equal to "milp" or "mf". We have to take into account that the “MILP” method is faster but very problematic numerically. Moreover, for the “MILP” method under non-constant returns to scale, a modification proposed by Zhu2018 is implemented.
In order to apply the “MILP” method, we have to compute a set of extreme efficient DMUs Charnes1991, i.e. DMUs spanning the facets of the efficient frontier that cannot be expressed as a positive linear combination of the other DMUs. We can find a set of extreme efficient DMUs from a deadata object by means of the function extreme_efficient, and pass this result to model_addmin through the parameter extreff. On the other hand, if we do not previously compute a set of extreme efficient DMUs, it is computed internally by model_addmin. For example, we compute the CRS additive-Min scores of the Fortune 500 dataset:
On the other hand, for applying the “maximal friends” method, we have to previously compute the maximal friends subsets Tone2010, i.e. the facets of the efficient frontier, by means of function maximal_friends (see Section (ref) for more details) and pass the result to model_addmin through the parameter maxfr.
Slacks-based measure (SBM) of efficiency models Tone2001 provide an efficiency score and the CRS version is given by
where the weights $w_i^-$, $w_r^+$ are positive with $\sum _{i=1}^mw_i^-=m$ and $\sum _{r=1}^sw_r^+=s$. Note that (ref) is expressed in a unit-invariant form and, moreover, it can be linearized using the Charnes-Cooper transformation. $\textrm{DMU}_o$ is efficient if and only if the optimal objective $\rho ^*=1$, i.e., the optimal slacks are all zero. The input and output-oriented CRS versions are given by:
respectively. In general, $\rho ^*_I\geq \rho ^*$ and $\rho ^*_O\geq \rho ^*$. Note that oriented SBM models (ref) do not serve to find efficient DMUs by their own because there can be inefficient DMUs with $\rho ^*_I=1$ or $\rho ^*_O=1$. As usual, different returns to scale can be easily considered in (ref) and (ref) by adding the corresponding constraints. SBM models are applied using model_sbmeff. Parameter orientation can be "no" (non-oriented, by default), "io" (input-oriented) or "oo" (output-oriented). The weights are introduced by means of weight_input and \texttt{weight_output}, whose default values are $1$. Moreover, according to Tone2001, SBM models in function \texttt{model_sbmeff} are automatically adapted to deal with zeros in data. We have to note that, for the specific case of zeros in output data, the Case 2 of Tone2001 is applied, but taking $1/100$ instead of $1/10$.
For example, we can replicate the results in Tone2001:
The original SBM efficiency model given by (ref) evaluates the inefficiency of a DMU referring to the efficient activity of the form $ \left( \mathbf{x}_o-\mathbf{s}^-,\mathbf{y}_o+\mathbf{s}^+\right) $ that produces the lowest score $\rho$. Hence, efficient targets may be far away from $\textrm{DMU}_o$ and they could be inappropriate. To overcome this issue, Tone2010 among others proposed to search the efficient activity of the form $ \left( \mathbf{x}_o-\mathbf{s}^-,\mathbf{y}_o+\mathbf{s}^+\right) $ that produces the highest score $\rho$, leading to the SBM-Max efficiency model:
Program (ref) is solved by means of the “maximal friends” (facets of the efficient frontier) technique Tone2010. Nevertheless, you may be careful because the SBM-Max efficiency model can produce non-monotonic scores. This model is applied setting the parameter kaizen = TRUE in model_sbmeff. For example, we can compare the CRS versions of the SBM-Max and the SBM-Min (original) models:
Moreover, you can find the maximal friends subsets of a given set of DMUs by means of the function maximal_friends. The result is a list with all the facets of the efficient frontier and the DMUs that compose them. Moreover, you can pass this result to function model_sbmeff through the parameter maxfr. For example,
Parameter silent in functions model_sbmeff and maximal_friends allows to hide the progress messages from the computation of the maximal friends.
The CRS cost, revenue and profit efficiency models Coelli2005 are given, respectively, by
where $\mathbf{c}$ and $\mathbf{p}$ are row vectors with the unit prices of inputs and outputs, respectively. Restricted versions of the cost and revenue efficiency models are given by adding the constraints $\mathbf{x}\leq \mathbf{x}_o$ to (ref) (a) and $\mathbf{y}\geq \mathbf{y}_o$ to (ref) (b). The cost, revenue and profit efficiency scores are given, respectively, by
$\textrm{DMU}_o$ is considered to be cost, revenue or profit efficient if its respective efficiency score is equal to $1$. Moreover, all efficiency scores are between $0$ and $1$, except for the case $\mathbf{p}\mathbf{y}_o<\mathbf{c}\mathbf{x}_o$, in which the profit efficiency score can be negative or $\geq 1$ (in this case, the higher the score, the greater the inefficiency). Different returns to scale can be easily considered by adding the corresponding constraints.
These models are applied using model_profit. Unit prices $\mathbf{c}$ and $\mathbf{p}$ are introduced by parameters price_input and price_output, respectively. As usual, they can be a value, a vector or a matrix of size [number of inputs/outputs]$\times$[number of DMUs in dmu_eval]. Restricted versions are considered setting parameter restricted_optimal = TRUE (by default). For example, the CRS restricted models can be solved for the Coelli_1998 dataset with this script:
We can determine with function make_deadata whether any of the inputs or outputs are either non-controllable, non-discretionary or undesirable variables (see Cooper2007 and Zhu2014 for more details about these special features). It can be done with the parameters nc_inputs/nc_outputs, nd_inputs/nd_outputs and \texttt{ud_inputs/ud_outputs}, which are integer numbers denoting the position of non-controllable, non-discretionary and undesirable inputs/outputs, respectively. For example, let us assume that the “Employees” input of the \texttt{Fortune500} dataset cannot be controlled by the decision maker, then we can flag this input as non-controllable by
Note that “nc_inputs = 3” does not mean “the third column is non-controllable”, but rather “the third input (Employees) is non-controllable”.
Non-controllable variables can not change their values. The input-oriented CCR model (ref) (a) and its second stage (ref) (a) are adapted:
where the superscripts $C$ and $NC$ refers to “controllable” and “non-controllable” respectively. Analogously, the output-oriented and directional models can be adapted, as well as the other non-radial models. Other returns to scale can be considered by adding the corresponding constraints.
Non-discretionary variables are exogenously fixed and therefore, it is not possible to vary them at the discretion of management Cooper2007. The input-oriented CCR model (ref) (a) and its second stage (ref) (a) are adapted:
where the superscripts $D$ and $ND$ refers to “discretionary” and “non-discretionary” respectively. Analogously, the output-oriented and directional models can be adapted. Other returns to scale can be considered by adding the corresponding constraints. For example, we can replicate the results in Ruggiero2007, where the second input is non-discretionary:
Another model that can be adapted for non-discretionary variables is the non-radial model by Fare1978. For example, the input-oriented CRS non-radial model (ref) (a) with its second stage (ref) (a):
Other non-radial models such as additive or SBM are not affected by non-discretionary variables.
An output is undesirable if producing less quantity of this output leads to more efficiency. By extension, an input is said to be undesirable if it behaves contrary to the other usual inputs, i.e., consuming more quantity of this input leads to more efficiency. Note that “undesirable” inputs are in fact “desirable”. A more grammatically correct denomination would be “good inputs” and “bad outputs”, but the term “undesirable” prevails in the literature. The production possibility set under CRS is defined by
where the superscripts $g$ and $b$ refers to “good” and “bad” respectively. Other returns to scale can be considered by adding the corresponding constraints. A $\textrm{DMU}_o$ is efficient in the presence of undesirable inputs/outputs if there is no vector $\left( \mathbf{x}^g,\mathbf{x}^b,\mathbf{y}^g,\mathbf{y}^b \right) \in P$ such that $\mathbf{x}^g_o\leq \mathbf{x}^g$, $\mathbf{x}^b_o\geq \mathbf{x}^b$, $\mathbf{y}^g_o\leq \mathbf{y}^g$, $\mathbf{y}^b_o\geq \mathbf{y}^b$ with at least one strict inequality.
In general, an undesirable output can be treated as an input, and vice versa. Nevertheless, this does not reflect the true production process and there can be interpretation issues in most models. In this case, models must be adapted to undesirable inputs/outputs. For example, a modified version of (ref) (non-oriented weighted SBM efficiency model under CRS) is
where $\mathbf{s}^-=(\mathbf{s}^{g-},\mathbf{s}^{b-})$, $\mathbf{s}^+=(\mathbf{s}^{g+},\mathbf{s}^{b+})$, $\mathbf{x}_o=(\mathbf{x}_o^{g},\mathbf{x}_o^{b})$, $\mathbf{y}_o=(\mathbf{y}_o^{g},\mathbf{y}_o^{b})$ Tone2003,Cooper2007,Tone2021. Note that some interpretation issues appear in this model because good inputs could generate negative efficiency scores.
Modified versions of the directional CRS model (ref) (a) and its second stage (ref) (b) are given in Fare2004:
Other returns to scale can be considered by adding the corresponding constraints.
For non-directional radial models under VRS, undesirable inputs/outputs are treated as proposed by Seiford2002. This technique consists of transform each undesirable input/output in this way:
where $\mathbf{u},\mathbf{v}$ are translation vectors that allow $^gx_{ij}$ and $^by_{rj}$ be positive. Usually, the “max + 1” translation is applied, i.e., $u_i=\max_{j}\left\{ x_{ij}^g\right\}+1$ and $v_r=\max_{j}\left\{ y_{rj}^b\right\}+1$. The VRS condition is advisable in order to assure translation invariance of the model. In function model_basic, parameters vtrans_i and vtrans_o correspond to the translation vectors $\mathbf{u}$ and $\mathbf{v}$, respectively. If vtrans_i[i] is NA, then it applies the “max + 1” translation to the i-th undesirable input. If vtrans_i is a scalar, then it applies the same constant translation to all undesirable inputs. If \texttt{vtrans_i} is \texttt{NULL}, then it applies the “max + 1” translation to all undesirable inputs (analogously for outputs). For example, we can replicate some results in Hua2007, where the third output is undesirable and the output translation parameter is set to $1500$:
Finally, function undesirable_basic transforms a deadata class object with undesirable inputs/outputs according to Seiford2002, making use of parameters vtrans_i and vtrans_o. This function also works with deadata_fuzzy class objects (see Section (ref)).
Efficiency models evaluate inefficient DMUs usually providing a score, but they do not discriminate between efficient DMUs. On the other hand, super-efficiency models precisely evaluate efficient DMUs in order to rank them.
Radial super-efficiency models consist of extracting the evaluated $\mathrm{DMU}_o$ from the evaluation reference set Andersen1993. For example, the input-oriented CCR super-efficiency model with its second stage is given by
where $X_{-o},Y_{-o}$ are the input and output data matrices, respectively, defined by $\mathcal{D}-\left\{ \textrm{DMU}_o\right\}$, $\bm{\lambda}=(\lambda_1,\ldots ,\lambda_{o-1},\lambda_{o+1},\ldots ,\lambda_n)^{\top}$, while the weights $\mathbf{w}^-$ and $\mathbf{w}^+$ are positive row vectors. Other returns to scale can be considered, but infeasibility problems may appear. Radial super-efficiency models can be applied using model_supereff.
With respect to non-radial models, the SBM super-efficiency (SSBM) models consist of projecting the evaluated $\mathrm{DMU}_o$ onto the part of the production possibility set defined by $\mathcal{D}-\left\{ \textrm{DMU}_o\right\}$ that consumes more inputs and produces less outputs than $\mathrm{DMU}_o$ Tone2002. For example, the non-oriented CRS version is given by
where $\mathbf{w}^-=\left( w_1^-,\ldots ,w_m^-\right) $ and $\mathbf{w}^+=\left( w_1^+,\ldots ,w_s^+\right) $ are weights, and $\mathbf{t}^-$, $\mathbf{t}^+$ are called super-slacks. Other returns to scale and orientations can be considered Tone2002. SSBM models are applied using the function model_sbmsupereff, in which parameters weight_slack_i and weight_slack_o stands for $\mathbf{w}^-$ and $\mathbf{w}^+$, respectively. As usual, they can be a value ($1$ by default), a vector or a matrix of size [number of inputs/outputs]$\times$[number of DMUs in dmu_eval]. For example, we can replicate the results in Tone2002, where the CCR super-efficiency and the input-oriented SSBM models are compared:
Du2010 adapted SSBM models to additive models, changing the objective $\delta $ in (ref) by $\mathbf{w}^-\mathbf{t}^-+\mathbf{w}^+\mathbf{t}^+$. These additive super-efficiency models are considered to be input-oriented if $\mathbf{w}^+=\mathbf{0}$, and output-oriented if $\mathbf{w}^-=\mathbf{0}$; if both weights are non-zero, they are non-oriented. Additive super-efficiency models can be applied using model_addsupereff, in which, as in SSBM models, parameters weight_slack_i and weight_slack_o stands for $\mathbf{w}^-$ and $\mathbf{w}^+$, respectively. By default, $w_i^-=1/x_{io}$ and $w_r^+=1/y_{ro}$, making the model unit invariant. Moreover, for comparison purposes, model_addsupereff also returns the corresponding score $\delta $ of the optimal solution. In fact, function efficiencies returns these $\delta $ scores. For example, we can verify that, taking the default weights, the input-oriented SSBM and additive super-efficiency models are equivalent (note that it does not hold for other orientations in general):
Cross-efficiency models evaluate the efficiency score of a DMU using the optimal weights of the other DMUs Doyle1994. For example, the input and output-oriented GRS cross-efficiency of $\textrm{DMU}_k$ based on the weights of $\textrm{DMU}_o$ are given by
where $\mathbf{v}_o^*,\mathbf{u}_o^*,\xi _{Lo}^*,\xi _{Uo}^*$ are optimal according to programs (ref) (a) or (b), respectively. Expressions for different returns to scale can be deduced from (ref) by taking $\xi _{Lo}^*=\xi _{Uo}^*=0$ (CRS), $L=U=1$ (VRS), $\xi _{Lo}^*=0$, $U=1$ (NIRS) or $L=1$, $\xi _{Uo}^*=0$ (NDRS). Then, the cross-efficiency score of $\textrm{DMU}_k$ is given by the average of the $k$-th column of matrix $E$:
Alternatively, $e_k$ can be computed excluding $E_{kk}$, i.e., without self-appraisal. Other interesting scores are the averages of the rows of $E$:
that can be also computed with or without self-appraisal. Finally, the Maverick index of $\textrm{DMU}_k$ is given by
Since optimal weights may not be unique, values of cross-efficiencies may vary depending on which optimal weights are chosen. In order to avoid this arbitrariness, two additional linear methods (II and III) with two different formulations (aggressive and benevolent) are implemented Doyle1994. For example, the aggressive input-oriented GRS method II applied to $\textrm{DMU}_o$ is given by:
and the corresponding method III:
The benevolent versions of (ref) and (ref) are given by maximizing instead of minimizing. Moreover, it is important to remark that the aggressive formulations of methods II and III under non constant returns to scale can lead to unbounded programs. In this case, deaR adds bound constraints automatically.
Finally, the correction proposed by Lim2015b can be applied in the input-oriented VRS model in order to fix negative cross-efficiency scores. This correction has been implemented analogously in the input-oriented NIRS and GRS models which can also give negative cross-efficiencies. For example, the corrected input-oriented GRS cross-efficiency of $\textrm{DMU}_k$ based on the weights of $\textrm{DMU}_o$ is given by
Method II is not affected by this correction, but method III can be adapted. In this way, program (ref) becomes
Cross-efficiency models can be applied using cross_efficiency. Apart from the usual datadea, orientation, rts, dmu_eval and dmu_ref, there are some other interesting parameters:
The output of cross_efficiency is a list with fields orientation, rts, L, U, selfapp, \texttt{correction}, \texttt{Arbitrary}, \texttt{M2_agg}, \texttt{M2_ben}, \texttt{M3_agg}, \texttt{M3_ben}, \texttt{data}, \texttt{dmu_eval}, \texttt{dmu_ref}, \texttt{epsi\-lon} and \texttt{modelname}. The results of arbitrary, method II (aggressive and benevolent) and method III (aggressive and benevolent) models are stored in \texttt{Arbitrary}, \texttt{M2_agg}, \texttt{M2_ben}, \texttt{M3_agg} and \texttt{M3_ben}, respectively. These fields have subfields \texttt{multiplier_input}, \texttt{multiplier\-_output}, \texttt{multiplier_rts} (for non constant returns to scale), \texttt{cross_eff}, \texttt{e}, \texttt{A} and \texttt{maverick}. Moreover, \texttt{efficiency} is stored in the field \texttt{Arbitrary}.
We can replicate the results in Golany1989:
For example, we can show the cross-efficiency scores (ref) for the benevolent formulation of method II:
Moreover, we can use function plot in order to visualize cross-efficiency matrices $E$ from different methods and formulations, as shown in Figure (ref).
Finally, we can replicate the results in Lim2015b and compare cross-efficiency scores for different methods:
Note that bound constraints are automatically added for unbounded methods under non constant returns to scale.
In all the DEA models seen so far, it is assumed that the data about the production process (i.e., inputs and outputs) are perfectly known, and thus, they can be considered deterministic. However, in practice, it is quite common that some degree of uncertainty is present in the data, and therefore methods that can deal with such vagueness must be defined.
There are mainly two approaches for dealing with non-deterministic data in DEA: stochastic and fuzzy models. The former uses statistical distributions to obtain a statistical characterization of the efficient frontier (see Olesen2016 and references therein for a review on stochastic DEA), while the latter uses fuzzy theory and membership functions to deal with the ambiguity in the data (see Hatami2011,Emrouznejad2014 for reviews on fuzzy DEA).
Currently, in the deaR package, three popular fuzzy models are implemented, namely: Kao-Liu, Guo-Tanaka and possibilistic models. Next, we shall describe those models and their implementation in deaR, after giving a brief introduction to fuzzy numbers.
A fuzzy set $A$ is defined by a function $\mu _A:\mathbb{R}\rightarrow \left[ 0,1\right] $, called membership function. For any $x\in \mathbb{R}$, the value $\mu _A(x)$ can be interpreted as the grade of membership of $x$ to $A$. Alternatively, a fuzzy set is completely determined by the so-called $\alpha$-cuts (or $h$-levels) which are defined by $A^{\alpha}=\left\{ x\in \mathbb{R}\,\,|\,\,\mu_A(x) \geq \alpha \right\} $ for $\alpha \in \left] 0,1\right]$, and $A^0=\overline{\left\{ x\in \mathbb{R}\,\,|\,\,\mu_A(x)>0\right\}}$, where the overline denotes clausure.
Fuzzy numbers are a particular case of fuzzy sets verifying:
In deaR, we only consider a special type of fuzzy numbers called trapezoidal, whose membership functions have a trapezoidal shape as shown in Figure (ref) (a). Thus, a trapezoidal fuzzy number is defined by the parameters $(\texttt{mL}, \texttt{mR}, \texttt{dL}, \texttt{dR})$. If $\texttt{mL}=\texttt{mR}$, the fuzzy number is called triangular; moreover, if $\texttt{dL}=\texttt{dR}$, it is called symmetric. Note that a crisp number is a degenerated case of symmetric triangular fuzzy number for which $\texttt{dL}=0$.
Let us assume that our DEA dataset contains some inputs/outputs which are trapezoidal fuzzy numbers. Then, for each one of those inputs/outputs, four columns must be defined, and they can be read with the make_deadata_fuzzy function. Parameters $(\texttt{mL}, \texttt{mR}, \texttt{dL}, \texttt{dR})$ corresponding to input variables are introduced by inputs.mL, inputs.mR, inputs.dL and inputs.dR, respectively. These parameters are numeric vectors of length $m$ (number of inputs), specifying in which column of the dataset the corresponding quantity is located. If any parameter is not defined for a given input, the corresponding position of the vector must be filled with NA. Analogously for outputs.
For example, Kao_Liu_2003 dataset Kao2003 contains one crisp input and five outputs, two of them (3rd and 5th) being triangular fuzzy numbers. This dataset should be read as follows:
Note that, in this example, only inputs.mL is needed because the input is crisp. Then, inputs.mR is taken equal to inputs.mL, and inputs.dL $=$ inputs.dR $=0$ automatically. On the other hand, outputs.mR is also not necessary because fuzzy outputs are triangular. The 1st, 2nd and 4th outputs are crisp, hence the 1st, 2nd and 4th entries of vectors \texttt{outputs.dL} and \texttt{outputs.dR} are set to \texttt{NA}. The 3rd and 5th entries of these vectors contain the column positions of the corresponding parameters (e.g., the 8th column of \texttt{Kao_Liu_2003} dataset contains the \texttt{dL} parameter of the 3rd output).
For the sake of completeness we are going to read, in the next example, a mixed fuzzy dataset with different types of the aforementioned fuzzy numbers as variables. The dataset FuzzyExample is included in deaR and contains 5 DMUs with 3 inputs and 3 outputs, whose types are given in Table (ref). Moreover, in this example, inputs and outputs corresponding to DMUs A and D are crisp numbers, but we need to write them as fuzzy numbers, i.e., $\texttt{mR}=\texttt{mL}$, and $\texttt{dL}=\texttt{dR}=0$ (see Table (ref)). The data reading is as follows:
Finally, make_deadata_fuzzy returns a deadata_fuzzy class object, whose structure is the same as a deadata object but including subfields mL, mR, dL, \texttt{dR} in the \texttt{input} and \texttt{output} fields. The \texttt{deadata_fuzzy} objects are passed to \texttt{modelfuzzy_xxx} functions in order to perform a DEA fuzzy analysis, as we are going to see in the next sections.
Kao-Liu models Kao2000a,Kao2000b,Kao2003 consist on applying an existing model in the worst and best scenarios for the evaluated DMUs at each $\alpha$-cut. The worst scenario for a DMU occurs when it consumes the largest possible input amounts and produces the smallest possible output amounts, while, on the contrary, the rest of the DMUs consume the smallest possible input amounts and produce the largest possible output amounts. On the other hand, the best scenario for a DMU occurs when it consumes the smallest possible input amounts and produces the largest possible output amounts, while the rest of the DMUs consume the largest possible input amounts and produce the smallest possible output amounts.
Hence, Kao-Liu models are in fact “metamodels” because another usual (crisp) model is needed. In deaR, Kao-Liu models are applied using modelfuzzy_kaoliu and the underlying model is selected by means of the parameter kaoliu_modelname, whose possible values are “basic”, “additive”, “addsupereff”, “deaps”, “fdh”, “multiplier”, “nonradial”, “profit”, “rdm”, “sbmeff”, “sbmsupereff” and “supereff”. Specific parameters for these models, such as orientation or returns to scale, can be also introduced.
Finally, $\alpha$-cuts are selected by parameter alpha, that is a numeric vector with the $\alpha$-cuts in $\left[ 0,1\right] $. Alternatively, if alpha $>1$, it determines the number of $\alpha$-cuts, equispatially distributed in $\left[ 0,1\right] $. For example,
applies Kao-Liu using the input-oriented BCC model. Note that alpha = 11 would produce the same $\alpha$-cuts. The function modelfuzzy_kaoliu returns an object of class dea_fuzzy containing all the information and parameters. The specific results of the applied submodel are stored in the field alphacut. Inside alphacut there are fields for each $\alpha$-cut and, inside these fields, the corresponding input/output data along with the results of each DMU are stored in the field DMU, as shown in Figure (ref).
One of the difficulties about fuzzy efficiency models is the representation of the efficiency scores. This is particularly cumbersome in Kao-Liu models, since the scores themselves can be non-trapezoidal fuzzy numbers. In the package deaR, some plot methods for objects of class dea_fuzzy are also implemented in function plot. For example, we can represent the results obtained from the Kao-Liu BCC model applied to the Leon_2003 dataset:
The results are depicted in Figure (ref). There, different types of DMUs can be found. For instance, F and H are completely inefficient, because no $\alpha$-cut contains the value 1 for the efficiency. B, D, and E are efficient only for some $\alpha$, but not for all, because 1 is contained in the $\alpha$-cuts only for sufficiently small values of $\alpha$. Finally, A, C, and G are efficient for all possible values of $\alpha$, since 1 is contained in all $\alpha$-cuts. However it is noteworthy to point out that A and C are “crisp-efficient” for large enough values of $\alpha$, in the sense that the $\alpha$-cut intervals are reduced to $\{1\}$ from a certain value of $\alpha$, while G always presents uncertainty.
Fuzzy models for symmetric triangular data under constant returns to scale were proposed by Guo2001. These models are implemented in modelfuzzy_guotanaka, specifically the input and output-oriented versions of the model in Guo2001. The fuzzy efficiencies are calculated using Guo2001. According to their notation, the $\alpha$-cuts are called $h$-levels, and the (crisp) relative efficiencies and multipliers for the level $h = 1$ are obtained from the multiplier model (model_multiplier). It is important to remark that the optimal solutions of the Guo-Tanaka models are not unique in general.
We can replicate the results in Guo2001:
The resulting scores are (non-symmetric) triangular fuzzy numbers for each $h$-level and can be extracted as usual with function efficiencies and plotted with plot (see Figure (ref)).
On the other hand, the dea_fuzzy object returned by the model stores the results in the field hlevel, whose structure is similar to the structure of alphacut shown in Figure (ref). The difference is that there are not Worst and Best fields, and hence, the fields DMU$_1$, \texttt{DMU$_2$}, $\ldots$ (in which there are stored the results about efficiencies and multipliers) hang directly from the field \texttt{DMU}.
As a complement to Guo-Tanaka models, we have implemented cross-efficiency fuzzy models (arbitrary, aggressive and benevolent formulations) in function cross_efficiency_fuzzy for dea_fuzzy objects returned by modelfuzzy_guotanaka. For example,
Alternatively, we can execute a Guo-Tanaka model internally, producing the same result:
Possibilistic fuzzy DEA models proposed by Leon2003 represent a generalization of the basic radial models to the fuzzy framework. By means of modelfuzzy_possibilistic, we can replicate the results in Leon2003:
Note that, as in Guo-Tanaka models, the $\alpha$-cuts are called $h$-levels but, in this case, efficiency scores are crisp numbers for each $h$-level. The results are stored in the field hlevel of the dea_fuzzy object returned by the model, and can be extracted as usual with functions efficiencies and lambdas. Moreover, function plot also works (see Figure (ref)).
Classical DEA models provide a still picture of the performance of the DMUs. However, the activities of the DMUs often vary over time and therefore other types of models are needed in order to analyse said time variation. One of the most popular techniques is the Malmquist methodology, in which the DMUs are evaluated with respect to efficient frontiers corresponding to different time periods and different returns to scale.
Let $X^t,Y^t$ be the input and output data matrices, where $t=1,\ldots ,T$ denotes the time period. We define the CRS production possibility set at $t$ in the contemporary (usual) way by $P^t=P\left( X^t,Y^t\right) $, according to (ref). Analogously, we define the VRS production possibility set $P_B^t=P_B\left( X^t,Y^t\right) $ according to (ref).
Given an activity $\left( \mathbf{x},\mathbf{y} \right) \in \mathbb{R}_{>0}^{m+s}$, we define input and output-oriented “distance” functions at $t=1,\ldots,T$ by
respectively. Note that for a DMU at $t$, (ref) (a) and (b) are the CCR efficiency score $\theta ^*$ and the inverse of $\eta ^*$ of the DMU (see (ref)), respectively. Moreover, we can define the VRS versions of (ref), $D_{IB}^t$ and $D_{OB}^t$ considering $P_B^t$ instead of $P^t$. In the following equations, $D$ can be $D_I$ or $D_O$, depending on the orientation.
According to Fare1994, the Malmquist index of DMU$_o$ at $t<T$ is given by
The CRS decomposition of (ref) is given by
with
being the technical change and the efficiency change of DMU$_o$ at $t$, respectively. On the other hand, the VRS decomposition of (ref) is given by
with $TC_o^t$ given by (ref) and
being the pure efficiency change and the scale change of DMU$_o$ at $t$, respectively.
According to Ray1997,Grifell1999, the Malmquist index of DMU$_o$ at $t<T$ is given by
The VRS decomposition of (ref) is given by (ref), with $PECH_o^t$ given by (ref) and
According to Grifell1999, a generalized Malmquist index is defined by (ref) with $TC_o^t$ given by (ref), $PECH_o^t$ given by (ref) and
Moreover, according to Fare1997 a biased Malmquist index can be computed from a biased technical change, given by
with
being the magnitude of technical change of DMU$_o$ at $t$, and
being the output and input bias indices of DMU$_o$ at $t$, respectively. Different returns to scale can be considered in the computation of the biased technical change, using $D$ (CRS) or $D_B$ (VRS) in (ref), (ref) and (ref).
All these indices can be computed in a sequential way (instead of the usual contemporary way) by considering the production possibility sets by $P^{\leq t}=P(X^{\leq t},Y^{\leq t})$ (CRS) and $P_B^{\leq t}=P_B(X^{\leq t},Y^{\leq t})$ (VRS), with $X^{\leq t},Y^{\leq t}$ being the input and output data matrices, respectively, of all the DMUs at time periods $\leq t$ Shestalova2003. The new corresponding “distance” functions can be denoted by $D^{\leq t}$ and $D_B^{\leq t}$, and they must be used instead of $D^t$ and $D_B^t$.
Finally, the technical change can be defined in a global way Pastor2005. Its CRS version is given by
where $D^{\leq T}$ denotes the “distance” function computed using a global production possibility set $P^{\leq T}$ considering all the DMUs of all time periods. A VRS version of (ref) is constructed by taking $D_B$ instead of $D$.
In order to read Malmquist datasets, the function make_malmquist can deal with datasets in wide or long formats. Wide format datasets are characterized by having the data for different time periods in different columns. In this case, parameter nper must contain the number of time periods and arrangement must be \texttt{"horizontal"} (by default). On the other hand, \emph{long} format datasets have the data for different time periods in the same column and hence, an extra column specifying to which time period the data belongs is required. In this case, parameter \texttt{percol} must contain the position of the column with the time periods and \texttt{arrangement} must be \texttt{"vertical"}. The rest of parameters are analogous to those in \texttt{make_deadata} function. For example,
for a wide format dataset, and
for a long format dataset. In both cases, the result is a list with the different time periods. In turn, within each time period, the data is stored as a deadata object.
Once we have read a dataset, we compute Malmquist and other indices using function malmquist_index. Parameter datadealist must contain the resulting list from function make_malmquist. As usual, parameters dmu_eval and dmu_ref indicate which DMUs are evaluated and with respect to which DMUs they are evaluated, respectively. The orientation is given by orientation, that can take values \texttt{"io"} (input-oriented, by default) or \texttt{"oo"} (output-oriented). Parameter \texttt{rts} indicate which decomposition is applied to the Malmquist index: (ref) for \texttt{"crs"} (by default) or (ref) for \texttt{"vrs"}.
Parameter type1 determines the way in which we compute the production possibility sets: "cont" (contemporary, by default), "seq" (sequential) or "glob" (global). On the other hand, parameter type2 determines the definition of the indices: "fgnz" Fare1994 (by default), \texttt{"rd"} Ray1997,Grifell1999, \texttt{"gl"} (generalized) Grifell1999 or \texttt{"bias"} (biased) Fare1997. Finally, \texttt{tc_vrs} is a logical parameter (\texttt{FALSE} by default) indicating if the biased technical change given by (ref) is computed under VRS or not.
Table (ref) shows expressions used for different parameters combinations, apart from orien\-tation. Expressions for type1 = "seq" are computed considering sequential production possibility sets ($P^{\leq t},P_B^{\leq t},\ldots$) instead of contemporary ones ($P^t,P_B^t,\ldots$). Expressions (ref), (ref), (ref) and (ref) can be computed under CRS or VRS. In the case type1 = "cont" or "seq", type2 = "bias" and rts = "vrs", expressions (ref) (ref) and (ref) are computed under CRS if \texttt{tc_vrs = FALSE} (by default) or under VRS if \texttt{tc_vrs = TRUE}. Moreover, if \texttt{type1 = "glob"} then \texttt{type2} is irrelevant.
Apart from the values of parameters datadealist, dmu_eval, dmu_ref, orientation, rts, type1, \texttt{type2} and \texttt{tc_vrs}, the list returned by \texttt{malmquist_index} contains all the indices involved in the computations: \texttt{mi}, \texttt{ec}, \texttt{tc}, \texttt{pech}, \texttt{sech}, \texttt{obtech}, \texttt{ibtech} and/or \texttt{matech}. Moreover, all the efficiency scores involved are stored in the field \texttt{eff_all}, so users can build their own indices:
where * can be crs or vrs (note that $D$ is replaced by $D_B$ in the vrs case). Moreover, if type1 = "seq" then $D^t$ and $D^{t+1}$ are replaced by $D^{\leq t}$ and $D^{\leq t+1}$, respectively.
We can replicate the results in Wang2011:
Moreover, we can also replicate the results in Grifell1999:
The bootstrap sampling method allows us to analyse the sensitivity of efficiency scores relative to variations in the efficient frontier. We have implemented the bootstrapping algorithm proposed by Simar1998 for basic radial models in function bootstrap_basic. Parameters datadea, orientation, rts, L and U acts as usual. On the other hand, parameter \texttt{B} indicates the number of bootstrap iterations ($2000$ by default) and \texttt{alpha} is a value between $0$ and $1$ ($0.05$ by default) determining the confidence intervals. Moreover, parameter \texttt{h} represents the bandwidth of smoothing windows. By default, \texttt{h} $ = 0.014$ but the optimal bandwidth factor can also be calculated following the proposals of Silverman1998 and Daraio2007. So, \texttt{h = "h1"} is the optimal bandwidth referred to as “robust normal-reference rule” Daraio2007, \texttt{h = "h2"} is the value of \texttt{h1} but with the factor $0.9$ instead of $1.06$, \texttt{h = "h3"} is the value of \texttt{h1} adjusted for scale and sample size Daraio2007, and \texttt{h = "h4"} is the bandwidth provided by a Gaussian kernel density estimate.
The result is a list with the following fields of interest:
We can replicate the results in Simar1998, but with $100$ iterations:
One of the grounding hypotheses in DEA is that the DMUs must be comparable, not only in the sense that all DMUs should use the same inputs to produce the same outputs, but also that their production processes should be defined under the same technology framework. However, there are cases in which the objective is precisely to compare the efficient frontiers of groups of DMUs operating under different technologies, and to determine the so-called “technology gap” between those frontiers and the theoretical potential frontier of all the DMUs as a whole, which is called the metafrontier.
There are several approaches to the metafrontier analysis. On one hand, the stochastic metafrontier analysis Battese2002 assumes a parametric relationship between inputs and outputs. On the other hand, there is the non-parametric metafrontier analysis, that can be either concave Battese2004 or non-concave Tiedemann2011.
In this section, we are going to illustrate with an example, how to compute the non-parametric metafrontiers using deaR. It is noteworthy that there is no specific function to do this, because it can be readily computed using the dmu_eval and dmu_ref parameters in the model_xxx functions (see Section (ref)). This example considers three groups of DMUs with one input X and one output Y operating under variable returns to scale (BCC model) with input orientation. Figure (ref) depicts the DMUs of this example, together with the three efficient frontiers and the resulting non-parametric metafrontiers.
First we create a synthetic dataset and we define the grouping:
The following code makes use of the dmu_eval and dmu_ref parameters in order to evaluate the efficiency score of each DMU with respect to the three different efficient frontiers:
The efficiency of a DMU with respect to the non-parametric non-concave metafrontier is the minimum of the aforementioned efficiency scores evaluated for that DMU, after removing NAs for infeasible problems:
Finally, the efficiency with respect to the non-parametric concave metafrontier is the usual efficiency score considering all the DMUs as the evaluation reference set:
deaR package allows both researchers and practitioners to make use of a wide variety of DEA models. Among the conventional DEA models, the user can choose between radial and non-radial models (directional, additive, SBM, etc.), with different orientations and returns to scale, as well as consider variables with special features (non-controllable, non-discretionary or undesirable inputs/outputs). In addition, the package includes super-efficiency, cross-efficiency, Malmquist index and bootstrapping models. On the other hand, the versatility of deaR allows the user, for example, to calculate the generalized Farrell measure or to perform a metafrontier analysis.
A new feature unique to this package is the inclusion of several fuzzy DEA models. Currently, it contains the Kao-Liu, Guo-Tanaka and possibilistic models. Specifically, the Kao-Liu model can be applied to all conventional DEA models implemented in deaR, both efficiency and super-efficiency.
In an effort to improve the communication of the results, the package includes novel graphical representations such as the references graph (see Figure (ref)), which aims to show the relationships among the DMUs that determine the efficient frontier and the inefficient DMUs. Also, in the case of the cross-efficiency analysis, the heat maps of the methods used are shown (see Figure (ref)). Moreover, the different fuzzy DEA models incorporate specific visualizations.
Finally, the architecture of the package makes it easy to upgrade with new features and models, e.g., stochastic or network DEA. Currently, we are also working on the design of new plots to aid the visualization of results, as we think this is a shortcoming in the DEA efficiency analysis.