03-IO_Env_Matrix.Rmd

---
output:
  pdf_document: default
  html_document: default
---

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

# Matrix Implementations of DEA

```{r Install-support-libraries, include=FALSE}
library (kableExtra)
```


In the previous chapter we created the DEA linear program algebraically using the `ompr` package. Before `ompr` was available, we needed to use another approach for most DEA packages such as \index{TFDEA} `TFDEA` and \index{MultiplierDEA} `MultiplierDEA`, which consisted of using matrices. This approach requires more "bookkeeping". Both approaches work but I think the algebraic approach is likely to be easier to write, debug, and extend in the future. For solving, we will use the widely adopted \index{lpSolveAPI} `lpSolveAPI` package, but this same matrix-based approach holds other LP solver engines such as C-Plex, Gurobi, and CBC.

People interested in creating their own package or working on contributing to a package that uses a matrix approach instead of an algebraic approach may find this chapter helpful. For people that want to stick to the algebraic implementation and using \index{ompr} `ompr`, you can skip this chapter.

To get it ready for matrix formulating, we need to review the standard form of linear programs which means that only numbers can be on the right hand side of the inequalities.
The standard form of linear programming can be described as the following.

$$
\begin{split}
 \begin{aligned}
    \text{min  }   & CX \\
    \text{s.t.:  } & AX \geq B\\ 
                       & X \geq 0  
  \end{aligned}
   \end{split}
$$

In this formulation, *X* is a vector of variables. This is a little confusing for people doing DEA since in DEA *x* represents \index{Inputs} inputs. In an envelopment DEA model, the vector *X* has length equal to the number of units examined plus one more for theta and the number of inputs and the number of outputs to cover slack variables. *C* is a vector of data with the same number of elements as the *X* vector. *A* is a matrix of data and *B* is a vector describing the constraint right hand sides.

Recall that in the previous chapter we had presented the standard LP model of the input-oriented envelopment model.

$$
\begin{split}
\begin{aligned}
    \text{min    } & \theta \\
    \text{s.t.:  } & \sum_{j=1}^{N^D} x_{i,j}\lambda_j \leq \theta x_{i,k} \; \forall \; i\\
                   & \sum_{j=1}^{N^D} y_{r,j}\lambda_j \geq  y_{r,k} \;\forall \; r\\
                   & \lambda_j \geq 0  \; \forall \; j
  \end{aligned}
   \end{split}
$$

The first step is to simply move the expression with a variable from the right hand side of the input inequality to the right-hand side.

$$
\begin{split}
 \begin{aligned}
    \text{min    } & \theta \\
    \text{s.t.:  } & \sum_{j=1}^{N^D} x_{i,j}\lambda_j - \theta x_{i,k} \leq 0 \; \forall \; i\\
                   & \sum_{j=1}^{N^D} y_{r,j}\lambda_j \geq  y_{r,k} \; \forall \; r\\
                   & \lambda_j \geq 0  \; \forall \; j
  \end{aligned}
   \end{split}
$$

Now we will convert the model's \index{Constraints!Input} input and \index{Constraints!Output} output constraints from inequalities into equalities by explicitly defining slack variables. This isn't necessary just yet but will be used later in the chapter. Note that when we do so, we will be adding more variable to the *X* vector.

$$
\begin{split}
 \begin{aligned}
    \text{min   } & \theta \\
    \text{s.t.: } & \sum_{j=1}^{N^D} x_{i,j}\lambda_j - \theta x_{i,k} + s^x_i = 0 \;\forall \; i\\
                  & \sum_{j=1}^{N^D} y_{r,j}\lambda_j - s^y_r =  y_{r,k} \; \forall \; r\\
                  & \lambda_j , \;s^x_i, \;s^y_r \geq 0  \; \forall \; i,\;r,\;j
  \end{aligned}
   \end{split}
  (\#eq:Ch3LPCCRIOE-Slack)
$$

Now, let's populate the information for our example from the beginning of chapter 2.

```{r  Defining-Problem-Data, echo=FALSE }
  x <- matrix(c(10,20,30,50),ncol=1,dimnames=list(LETTERS[1:4],"x"))

  y <- matrix(c(75,100,300,400),ncol=1,dimnames=list(LETTERS[1:4],"y"))
  
  ND <- nrow(x); NX <- ncol(x); NY <- ncol(y); 
  
  kbl (cbind(x,y),booktabs=T,
      caption="Input-Output Data")
```

The next step is to incorporate the data into the model. Since there is only one input and one output, there are only two constraints and two slack variables. To keep a close connection to the application, we will refer to the specific envelopment variables by DMU name rather than number. For example, this lists $\lambda_A$ rather than $\lambda_1$.

$$
\begin{split}
 \begin{aligned}
    \text{min    } & 0\lambda_A + 0\lambda_B + 0\lambda_C + 0\lambda_D + 
                          1\theta + 0s^x_1 + 0s^y_1\\
    \text{s.t.:  } & 10\lambda_A + 20\lambda_B + 30\lambda_C + 50\lambda_D -
                          20\theta + 1s^x_1 + 0s^y_1 = 0\\
                   & 75\lambda_A + 100\lambda_B + 300\lambda_C + 400\lambda_D + 
                           0\theta + 0s^x_1 - 1s^y_1 = 100\\
                   & \lambda_A, \; \lambda_B, \; \lambda_C, \; \lambda_D, \; \theta, \; s^x_1, \; s^y_1 \geq 0  \;
  \end{aligned}
   \end{split}
$$

Notice that this formulation now includes seven \index{Decision variables} decision variables. The linear programming convention is to define the vector of all variables to be *X*. In this case, $X=[\lambda_A, \lambda_B, \lambda_C, \lambda_D, \theta, s^x_1, s^y_1 ]^\tau$. The $\tau$ is a transpose operator that changes a vector from a column vector to a row vector or row to column. In this case, the *X* is a column vector but the elements are shown as a horizontal row vector to save space and then flipped using the transpose operator. Using the *X* vector for decision variables allows us to separate the data from the variables into the following.

$$
\begin{split}
 \begin{aligned}
    \text{min }   & \begin{bmatrix}
         0 & 0 & 0 & 0 & 1 & 0 & 0
     \end{bmatrix} X \\
    \text{s.t.:  } & \begin{bmatrix}
         10 & 20 & 30 & 50 & -20 & 1 & 0 & \\ 
          75 & 100 & 300 & 400  & 0 & 0 & -1 &
     \end{bmatrix}  X \begin{matrix}
=
\end{matrix} \begin{bmatrix}
0\\ 
100
\end{bmatrix}\\ 
   & X \geq 0  
  \end{aligned}
\end{split}
$$

Now, let's show the data elements. The *C* vector contains the coefficients of the objective function.

$$
\begin{split}
    C = \begin{bmatrix}
         0 & 0 & 0 & 0 & 1 & 0 & 0
     \end{bmatrix}  
\end{split}
$$

The *A* matrix contains the variable's coefficients from the constraints.

$$
\begin{split}
    A = \begin{bmatrix}
         10 & 20 & 30 & 50 & -20 & 1 & 0 & \\ 
          75 & 100 & 300 & 400  & 0 & 0 & -1 &
     \end{bmatrix}  \\
\end{split}
$$

Lastly, the *B* vector contains the right hand sides of the constraints.

$$
\begin{split}
     B = \begin{bmatrix}
0\\ 
100
\end{bmatrix} 
\end{split}
$$

The bulk of the LP is then defined by passing along the numerical A matrix and the vectors of B and C. We then just need to indicate the sense of the objective function (min or max) and the directions of the constraints.

As models get larger, the full numerical representation becomes difficult for humans to read even if it is a very good format for computers. If you stare at that too long, you may go "bug-eyed." Personally, I really like the algebraic representations. The \index{lpSolve} `lpSolve`package in R is capable of taking an LP model in this form and solving it. Now we just need to structure the data in R to put it in the same format.

```{r warning=FALSE, message=FALSE}
library (lpSolve)
Amatrix <- rbind(c(10, 20, 30, 50, -20, 1, 0),
                 c(75, 100, 300, 400, 0, 0, -1))
Bvector <- c(0, 100)
Cvector <- c(0, 0, 0, 0, 1, 0, 0)
constraintdir <- c("=", "=")

res3_1 <- lp ("min", Cvector, Amatrix, constraintdir, Bvector)

res3_1
```

Now let's look at the optimal values of the \index{Decision variables} decision variables.

```{r}
format(res3_1$solution, digits=4)
```

The solver in this case knew nothing of our model's structure so the results in terms of decision variable values would need be organized based on how we structured the model to pass into \index{lpSolve} `lpSolve`.

```{r}
res3_1.solution <- as.matrix(res3_1$solution)
rownames (res3_1.solution) <- c("L_A", "L_B", "L_C", "L_D", "theta", 
                                 "s_x", "s_y")
kbl (t(res3_1.solution), booktabs=T, digits=4,
        caption="Optimal solution values for DMU A (CCR-IO") |>
  kable_styling(latex_options = "HOLD_position")
```

Let's demonstrate the use of \index{LaTeX} LaTeX in table captions. There are several things to bear in mind when using LaTeX rendered in column labels.

-   The LaTeX must be wrapped in \$ symbols when used in-line for either column labels or in-line text.
-   A double backslash is needed to show symbols
-   When the plain text version is displayed such as by running a code chunk, the plain text version will display the LaTeX commands making for a harder to read table.
-   Using the proper mathematical notation in the table gives a nice and rigorous display.
-   Superscripts and subscripts are triggered with a `^` or `_` respectively but this is done for just immediately following character. To have multiple characters in a superscript or subscript, they need to be wrapped in curly brackets.  

```{r, results = 'asis'}
res3_1.solution <- as.matrix(res3_1$solution)
rownames (res3_1.solution) <- c("$\\lambda_A$", "$\\lambda_B$", 
                                "$\\lambda_C$", "$\\lambda_D$", 
                                "$\\theta^{CRS}$", 
                                "$s^x_1$", "$s^y_1$")
kbl (t(res3_1.solution), booktabs=T, escape=F, digits=4,
        caption=
  "Optimal solution values for DMU A with LaTeX column labels")|>
  kable_styling(latex_options = "HOLD_position")
```

## Creating the LP Incrementally

Another approach is to incrementally build up the model. Let's use the \index{lpSolveAPI} `lpSolveAPI` package to demonstrate this way of doing the analysis.

```{r}
library(lpSolveAPI)
```

Now it gets more involved. We are defining a function for calculating DEA with options for \index{Returns to scale} returns to scale\index{Orientation}, and two different options for \index{Super-efficiency} super-efficiency.

For this walk through, we will focus on running the LP for calculating DEA for only one unit. We typically would do a loop for every unit but the entire `for` loop would need to be inside a single R chunk of this document. Instead, we will walk through how the data structure works. If you understand this, extending it to other units with a `for` loop is straightforward.

First, let's declare each of the needed basic data structures. First, determine the number of DMUs, number of inputs, and number of outputs.

```{r Defining-Problem-Size}
x <- matrix(c(10,20,30,50), ncol=1,
            dimnames=list(LETTERS[1:4],"x"))

y <- matrix(c(75,100,300,400), ncol=1,
            dimnames=list(LETTERS[1:4],"y"))
  
ND <- nrow(x); NX <- ncol(x); NY <- ncol(y); 

kbl(cbind(x,y), booktabs=T,
    caption="Input-Output Data")
```

Now we will declare the data structures for inputs and outputs that will contain the data that we use. It is a little redundant to create new objects for the x and y data but there are advanced cases when it is helpful.

```{r Creating-Alternative-Data}
xdata <- x[1:ND,] 
dim(xdata) <- c(ND,NX) 
ydata <- y[1:ND,]
dim(ydata) <- c(ND,NY)
```

Now we can define the structure of the results that we will be calculating.

This is overkill since it has space for results of every unit instead of just the one that we are analyzing now.

```{r Declaring-Structures-of-Results}
results.efficiency <- matrix(rep(-1.0, ND), nrow=ND, ncol=1)
results.lambda     <- matrix(rep(-1.0, ND^2), nrow=ND,ncol=ND)
results.vweight    <- matrix(rep(-1.0, ND*NX), nrow=ND,ncol=NX) 
results.uweight    <- matrix(rep(-1.0, ND*NY), nrow=ND,ncol=NY) 
results.xslack     <- matrix(rep(-1.0, ND*NX), nrow=ND,ncol=NX) 
results.yslack     <- matrix(rep(-1.0, ND*NY), nrow=ND,ncol=NY) 
```

Now it is time for more model options.

This is where we start creating the populating the linear programming model. This model declaration creates a problem with the needed number of \index{Variables} columns (variables) for standard DEA envelopment models.

```{r Create-lp-opject}
# This creates an empty LP object that we will later populate
    lpe <- make.lp(0,(ND+1+NX+NY))   
```

The name `lpe` is meant to indicate a linear programming object for the envelopment model. There are two basic DEA models: the \index{Envelopment model} envelopment model and the \index{Multiplier model} multiplier model. They generate the same results but are structured very differently. We'll return to the multiplier model at a later time.

We easily display what our empty `lpe` object looks like. We are initializing it without any constraints but columns for each variable that we expect to use:

-   Columns for each lambda, *i.e.* $\lambda_1$, $\lambda_2$, ..., $\lambda_{N^D}$\
-   A column for the efficiency, $\theta$
-   Columns for each input slack and each output slack

```{r Display-lp-Empty}
lpe 
```

The empty LP isn't very interesting but building it up step by step shows what we are doing and why.

This is an important step where we say that we are picking a particular unit to examine. Since our earlier graphical analysis pointed to unit *B*, we will also look at unit 2. Changing to a different unit is as simple as changing the number 2. Of course, even easier is to use a `for` loop to run through doing it for every unit but that simple change is something that we will also leave for later. The current format lets me continue to explain how DEA works one chunk of R code at a time.

```{r Which-Unit}
k <- 2   # Declare a particular unit instead of 
         # using a for loop such    for (k in 1:ND)     
```

The next thing that we will do is to create the objective function.

In this case, the "1" corresponds to the column for the $\theta$ variable of the input-oriented model. The input-oriented model is a min. The output-oriented model is a max. Let's make sure to set the model up accordingly.

```{r}
lpcontrol <- lp.control(lpe, sense="min") 
              # Set to a minimize function

set.objfn(lpe, c(rep(0, ND), 1, rep(0,NX+NY)))  
              # Minimize the theta variable.
```

The above commands have two important functions. The first merely sets it to be a *minimize* objective function. The second puts the relevant 0's and 1's in the objective function.

The `lp.control` function returns a long list of options being used for the linear programming solver. The defaults are fine for now and it can take a full page to list all the settings so I just store it into an object, `lpcontrol` so it doesn't display it for now.

```{r Display-lp-with-Obj}
lpe 
```

Now we will prepare for adding names to elements of our data matrix. This is a little tricky so don't worry if these commands are initially confusing. Also, if you are fluent in R, you can probably come up with more elegant ways of doing this.

```{r Add_names_to_LP}

lambda_names <- paste("L", 1:ND, sep = "")  
# We'll use L to represent Lambda
  
eff_name <- "theta"   # Pick the appropriate variable name

xslack_names <- paste ("sx", 1:NX, sep = "")  
     # These are slack variables instead of inequalities
yslack_names <- paste ("sy", 1:NY, sep = "")

ColNames <- c(lambda_names, eff_name, xslack_names, yslack_names)

x_constraint_names <- paste("Con_X_", 1:NX, sep = "")
y_constraint_names <- paste("Con_Y_", 1:NY, sep = "")
    
RowNames <- NULL
 
```

The following step is more complicated but the core idea is to populate a data matrix to describe the constraints for a DEA model. The first `for` loop create the row(s) corresponding to input constraints.

```{r Build-lp-x-constraints}
for (i in 1:NX) {
  add.constraint(lpe, c(xdata[,i],-1*xdata[k,i], 
                        if(i>1) rep(0,i-1),+1,
                        if(i<NX) rep(0,NX-i),
                        rep(0,NY)),"=", 0)
          }
RowNames <- c(RowNames, x_constraint_names )
dimnames(lpe) <- list(RowNames, ColNames)  
   # Add X Constraint names to LP

lpe 
```

Now the rows \index{Constraints!Rows} (constraints) are created for the outputs.

```{r Build-lp-y-constraints}

for (r in 1:NY) {
  add.constraint (lpe, c(ydata[,r], 0, rep(0,NX), 
                         if(r>1) rep(0,r-1),-1,
                         if(r<NY) rep(0,NY-r)),"=", 
                  ydata[k,r]) 
          }

RowNames <- c(RowNames, y_constraint_names )
dimnames(lpe) <- list(RowNames, ColNames)  
# Add X Constraint names to LP

lpe 
```

Set the lower limits for all variables to be 0.0

```{r}
set.bounds (lpe, lower = rep(0.0, ND+1+NX+NY))  
```

All the variables should be non-negative. (Technically the radial efficiency scores, $\theta$ is usually declared to be un-restricted in sign but this won't affect the solution.)

```{r Display-lp-Non-Negativity-constraints }
lpe 
```

## Solving and Processing Results

Finally, solve the LP to calculate the efficiency score.

```{r}
solve.lpExtPtr (lpe)
```

Simply solving doesn't tell us anything very useful. It just gives an error code of 0 (or at least we hope to get that-anything else is worth digging into!) An error code of 0 means that it solved the problem. When it does so, it adds results to the `lpe` object. Now we need to extract the results from the solution in a form that we can use.

We have four separate steps. In the input-oriented model, we read out the objective function value and place it in the *k*'th place of `results.efficiency`.

The following step sets *Phase1* equal to the same efficiency score. We will use this value in a second LP that we creatively refer to as *Phase2* in order to make sure that we have an optimal target. This is consistent with the two-phase approach discussed in Chapter 2.

Next we get \index{Dual values!Weights} input dual weights and place them in `results.vweight`. These are useful because the dual weights give us the solution to the DEA multiplier model. We can also think of them roughly as the "prices" or "values" placed on inputs and outputs that make the unit appear as good as possible. (It is tricky that all rows are offset from how you might count them because row 1 is considered to be the objective function.)

The fourth item gets the dual weight from the returns to scale constraint. Notice that this is extracted from row 2 of the dual solution.

The same steps are repeated for the output-oriented model but sign corrections may need to be made by multiplying these items by *-1*, you can check for yourself if this needs to be done.

The last step of extracting output weights is common for both orientations since it does not require a sign correction.

```{r Ch3-Collect-Phase1-Results}
results.efficiency[k] <-  get.objective(lpe)      
Phase1 <- get.objective(lpe)
results.lambda[k,] <- get.variables(lpe)[1:(ND)] 
    # Put lambda values in matrix            

results.vweight[k,] <- -1 * get.dual.solution(lpe)[2:(1+NX)] 
    # Get input weights, v
      
results.uweight[k,] <- get.dual.solution(lpe)[(2+NX):(1+NX+NY)]
    # Get output weights, u
DMUnames <- list(c(LETTERS[1:ND]))
```

Let's look at our results!

```{r}
Xnames <- lapply(list(rep("X",NX)),paste0,1:NX)
Ynames <- lapply(list(rep("Y",NY)),paste0,1:NY)
Vnames <- lapply(list(rep("v",NX)),paste0,1:NX)
Unames <- lapply(list(rep("u",NY)),paste0,1:NY)
SXnames <- lapply(list(rep("sx",NX)),paste0,1:NX)
SYnames <- lapply(list(rep("sy",NY)),paste0,1:NY)
dimnames (results.lambda) <- c(DMUnames,DMUnames)
dimnames (results.xslack) <- c(DMUnames,SXnames)
dimnames (results.yslack) <- c(DMUnames,SYnames)
dimnames (results.vweight) <- c(DMUnames,Vnames)
dimnames (results.uweight) <- c(DMUnames,Unames)

results.efficiency [k]  # Display efficiency score for unit k
results.lambda  [k,]    # Display lambda variables
results.vweight [k,]    # Display input weights
results.uweight [k,]    # Display output weights
```

```{r, echo=FALSE}
kbl (cbind(results.efficiency [k],t(results.lambda  [k,]),
           results.vweight [k,],results.uweight [k,]),
     booktabs=T, escape=F, digits=4, 
     caption="DEA Efficiency Scores and Lambdas", 
     col.names =  c("$\\theta^{CRS}$","$\\lambda_A$", "$\\lambda_B$", 
                                "$\\lambda_C$", "$\\lambda_D$", 
                                "$v_1$", "$u_1$") )|>
  add_header_above(c("Efficiency" = 1, "Envelopment Variables" = 4, 
                     "Weights" = 2))|>
  kable_styling(latex_options = "HOLD_position")
```

## Slack Maximization

\index{Slack|(}

This is a little more advanced. In analyses where our only purpose is to find the efficiency score, we could stop here. In cases where we want to find the *best* target for an inefficient unit, we need to go one step further and say that we want to also maximize the slack given that same level of efficiency. The simple one input, one output case with constant returns to scale won't demonstrate a problem so it is only shown here at this time for completeness and the particularly adventurous reader.

Let's explain a little bit of slack maximization. Phase 1 gives simple radial efficiency scores but does not necessarily give the best target for efficiency. The purpose here is to ensure that the $\lambda$ values actually point to the best possible target of performance for each unit.

Now we need to do a step that is called slack maximization. We will need to add a constraint that the efficiency variable ($\theta$) is held fixed to the value found in Phase 1.

The linear program is modified in the following manner. We now want to maximize the sum of both input slack(s) and output slack(s).

$$
\begin{split}
 \begin{aligned}
    \text{max  }  & \sum_{i}s^x_i + \sum_{r}s^y_r\\
    \text{s.t.  } & \sum_{j=1}^{N^D} \lambda_j  = 1\\
                       & \sum_{j=1}^{N^D} x_{i,j}\lambda_j - \theta x_{i,k} + s^x_i = 0 \; \forall \; i\\
                       & \sum_{j=1}^{N^D} y_{r,j}\lambda_j - s^y_r =  y_{r,k} \; \forall \; r\\
                       & \theta = \theta^*\\
                       & \lambda_j , s^x_i, s^y_r \geq 0  \; \forall \; j
  \end{aligned}
   \end{split}
  (\#eq:Ch3LPCCRIOE-MaxSlacks)
$$

The following are the changes to the LP object.

```{r}
   # Maximize sum of slacks
add.constraint (lpe,c(rep(0,ND),1,rep(0,NX+NY)),"=",Phase1)
   # Hold Efficiency constant
lpcontrol <- lp.control(lpe, sense="max")  
   # Change from min to max
set.objfn (lpe,c(rep(0,ND+1),rep(1,NX+NY))) 
   # Maximize sum of slacks
```

We are now ready to solve phase 2. We can then extract the values of the $\lambda$ variables and the slack variables.

```{r}
solve.lpExtPtr (lpe)                                        
results.lambda[k,] <- get.variables(lpe)[1:(ND)]            
results.efficiency[k] <- get.variables(lpe)[ND+1]
results.xslack[k,] <- get.variables(lpe)[(ND+2):(ND+1+NX)]  
results.yslack[k,] <- get.variables(lpe)[(ND+2+NX):(ND+1+NX+NY)] 
```

Again, let's look over the results. We will arrange them in the order used for building the vector of \index{Decision variables} decision variables.

```{r}
kbl(cbind(t(results.lambda[k,]), results.efficiency[k], 
    results.xslack[k,], results.yslack[k,]), 
    booktabs = T, digits = 2, escape=F,
     caption="Analysis results for DMU B",
    col.names =  c("$\\lambda_A$", "$\\lambda_B$", 
                   "$\\lambda_C$", "$\\lambda_D$", 
                   "$\\theta^{CRS}$","$s_x$", "$s_y$") )|>
  add_header_above(c("Envelopment Variables" = 4, 
                     "Efficiency" = 1, "Slacks" = 2))|>
  kable_styling(latex_options = "HOLD_position")
```

Now, let's put it all together to summarize the results.

What does this mean? We are saying that unit `r k` (B) has an efficiency of `r results.efficiency[k]`. This is obtained by using a vector $\lambda$=`r results.lambda[k,]`. Also, the weights on input(s) that correspond to this score was `r results.vweight[k,]` and for outputs `r results.uweight[k,]`.

In this case there was no slack among the inputs or outputs to be maximized but the same approach could be used for other applications.

\index{Slack|)}

## Exercises

If you have followed along this far - you deserve a medal! If you work your way through by cutting and pasting code fragments, you should be able to reproduce my exact results. At this point you are ready for a few challenges:

1.  Graphically determine the efficiency score and lambdas for the other units. *Difficulty: Decaf cup of coffee*
2.  Interpret the solution in terms of who is doing well, who is doing poorly, and who should be learning from whom. *Difficulty: Cup of coffee*
3.  Add a fifth unit, E, that produces 400 units output using 30 units of input. Graphically evaluate all five units for their efficiency scores and lambda values.
4.  Interpret the solution in terms of who is doing well, who is doing poorly, and who should be learning from whom. *Difficulty: Cup of coffee*
5.  Examine another unit using the R code (hint: change where k is set.) *Difficulty: Decaf cup of coffee*
6.  Interpret the solution in terms of who is doing well, who is doing poorly, and who should be learning from whom. *Difficulty: Cup of coffee*
7.  Wrap a `for` loop around the model to examine every unit. *Difficulty: Cup of coffee*
8.  Use a bigger data set (more inputs, outputs, and units.) *Difficulty: Cup of coffee*
9.  Validate results against other DEA packages (ex. `Benchmarking`, `nonparaeff`) *Difficulty: Pot of coffee*
10. Construct an example where Phase 2 increases positive slacks from Phase 1. *Difficulty: Pot of coffee*
11. Create "cool" graphs or plots of results. *Difficulty: It depends...*