05-Multiplier.Rmd

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

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

# Multiplier Models

## Introduction

\index{Multiplier model|(}

Every linear program has a related alter ego referred to as the
\index{Dual values} dual. By duality, the two models have the same
optimal objective function values. In DEA, the multiplier models are
simply the dual of the envelopment models. In terms of the matrix
representations of DEA given in chapter 3, the A matrix is transposed,
rows give way to columns and columns become rows. The right hand side
values appear in the objective function while the previous objective
function coefficients become right hand side values. More could be said
about this but this will be deferred for people interested in duality or
more algorithmic aspects of linear programming.

Let's turn our attention instead to intuitively developing the DEA
multiplier model through a few steps.

## The Two Output, No Input Model

Assume that you, Alfred, are one of six restaurant managers of various
franchisees of the popular TexMex chain, Chiposter. Each restaurant
branch is operating under different conditions and are operated by your
colleagues: Barb, Chris, Don, Ed, and Fran. Each restaurant manager is
achieving outcomes (outputs.) At the end of the year, you want to make a
case for having done well.

```{r multiplier-rest-Y-data }
YRest1 <- matrix(c(25, 14, 42, 18, 19, 14, 3500, 2100, 1050,
                   4200, 2500, 1500), ncol=2,
                 dimnames=list(c("Alfred", "Barb","Chris", 
                                 "Don", "Ed","Fran"),
                               c("$y_1$", "$y_2$")))

ND <- nrow(YRest1); NY <- ncol(YRest1) # Define data size
```

```{r, include=FALSE}
library(kableExtra)
```

```{r, echo=TRUE}
kbl (YRest1, caption="Restaurant Data", booktabs=T, escape=F)  |>
  kable_styling(latex_options = "HOLD_position")
```

You can select any way of weighting the two outputs with just two
conditions:

-   No weight can be negative
-   No one can receive a weighted score above 100%

Given these two limitations, your goal is to make your "Alfred score" as
high as possible.

While it is possible to experiment with different weights, this can also
be framed as an optimization model. The goal or objective is to find the
best Alfred score, using weights ($u_1$, and $u_2$) subject to the
constraints that no one's scores is above 100%.

Let's take the restaurant data and formulate the model.

$$
\begin{split}
\begin{aligned}
    \text {max }   & 25 u_1 + 3500 u_2 \\
    \text{s.t.:  } & 25 u_1 + 3500 u_2 \leq 1.0\\
                   & 14 u_1 + 2100 u_2 \leq 1.0\\
                   & 42 u_1 + 1050 u_2 \leq 1.0\\
                   & 18 u_1 + 4200 u_2 \leq 1.0\\
                   & 19 u_1 + 2500 u_2 \leq 1.0\\
                   & 14 u_1 + 1500 u_2 \leq 1.0\\
                   & u_1, u_2\geq 0  \
  \end{aligned}
\end{split}
$$

Next let's implement this in \index{ompr} `ompr` which may serve as a
convenient refresher as well. As usual, we will load the needed
packages.

```{r, message=FALSE}
library(ROI, quietly=TRUE)              
library(ROI.plugin.glpk, quietly=TRUE)  
library(ompr, quietly=TRUE)             
library(ompr.roi, quietly=TRUE)     
```

```{r}

result <- MIPModel ()                         |>
  add_variable (u1, type="continuous", lb=0)  |>
  add_variable (u2, type="continuous", lb=0)  |>
  set_objective (25*u1+3500*u2, "max")        |>
  add_constraint (25*u1+3500*u2 <= 1.0)       |>
  add_constraint (14*u1+2100*u2 <= 1.0)       |>
  add_constraint (42*u1+1050*u2 <= 1.0)       |>
  add_constraint (28*u1+4200*u2 <= 1.0)       |>
  add_constraint (19*u1+2500*u2 <= 1.0)       |>
  add_constraint (14*u1+1500*u2 <= 1.0)       |>
  solve_model(with_ROI(solver = "glpk"))

u1_result <- get_solution(result, u1) 
u2_result <- get_solution(result, u2) 
alfred_score <- 25 * u1_result + 3500 * u2_result
names(alfred_score)<-""
```

```{r, echo=FALSE}
kbl (cbind(alfred_score, u1_result, u2_result), 
     escape=F, booktabs = T,
     col.names = c("Score", "$u_1$", "$u_2$"),
     caption="Alfred's Score and Results")       |>
  kable_styling(latex_options = "HOLD_position")
```

We could repeat the process for each person, by simply changing the
objective function.

```{r}
result <- MIPModel ()                         |>
  add_variable (u1, type="continuous", lb=0)  |>
  add_variable (u2, type="continuous", lb=0)  |>
  set_objective (14*u1+2100*u2, "max")        |>
  add_constraint (25*u1+3500*u2 <= 1.0)       |>
  add_constraint (14*u1+2100*u2 <= 1.0)       |>
  add_constraint (42*u1+1050*u2 <= 1.0)       |>
  add_constraint (28*u1+4200*u2 <= 1.0)       |>
  add_constraint (19*u1+2500*u2 <= 1.0)       |>
  add_constraint (14*u1+1500*u2 <= 1.0)       |>
  solve_model(with_ROI(solver = "glpk"))

u1_result <- get_solution(result, u1) 
u2_result <- get_solution(result, u2) 
barb_score <- 14 * u1_result + 2100 * u2_result
names(barb_score)<-""
```

```{r, echo=FALSE}
kbl (cbind(barb_score, u1_result, u2_result), 
     escape=F, booktabs = T, digits = 5,
     col.names = c("Score", "$u_1$", "$u_2$"),
     caption="Barb's Score and Results")           |>
  kable_styling(latex_options = "HOLD_position")
```

Let's do it for Chris.

```{r}
result <- MIPModel ()                         |>
  add_variable (u1, type="continuous", lb=0)  |>
  add_variable (u2, type="continuous", lb=0)  |>
  set_objective (42*u1+1050*u2, "max")        |>
  add_constraint (25*u1+3500*u2 <= 1.0)       |>
  add_constraint (14*u1+2100*u2 <= 1.0)       |>
  add_constraint (42*u1+1050*u2 <= 1.0)       |>
  add_constraint (28*u1+4200*u2 <= 1.0)       |>
  add_constraint (19*u1+2500*u2 <= 1.0)       |>
  add_constraint (14*u1+1500*u2 <= 1.0)       |>
  solve_model(with_ROI(solver = "glpk"))

u1_result <- get_solution(result, u1) 
u2_result <- get_solution(result, u2) 
chris_score <- 42 * u1_result + 1050 * u2_result
names(chris_score)<-""
```

```{r, echo=FALSE}
kbl (cbind(chris_score, u1_result, u2_result), 
     escape=F, booktabs = T, digits=5,
     col.names = c("Score", "$u_1$", "$u_2$"),
     caption="Chris's Score and Results")|>
  kable_styling(latex_options = "HOLD_position")
```

This process of building the model for each person individually and hard
coding the data into the model makes it difficult to maintain,
generalize, and apply to other cases.

Let's now generalize this by building the model algebraically. Let's
define the data, $y_{r,j}$ to be the value of output *r* for manager
*j*. Therefore, in this application, $y_{1,1}=25$ and $y_{2,1}=3500$
reflects Alfred (*j=1*) has outputs *1* and *2* of 25 and 3500
respectively.

We can then rewrite the formulation to the following. The summation is
over the two outputs.

$$
\begin{split}
\begin{aligned}
    \text {max } & \sum_{r=1}^{2} u_r y_{r,1}\\
    \text{s.t.:  } &  \sum_{r=1}^{2} u_r y_{r,1} \leq 1.0 \quad \text{[Alfred]}\\
     &  \sum_{r=1}^{2} u_r y_{r,2} \leq 1.0  \quad \text{[Barb]}\\
     &  \sum_{r=1}^{2} u_r y_{r,3} \leq 1.0  \quad \text{[Chris]}\\
     &  \sum_{r=1}^{2} u_r y_{r,4} \leq 1.0  \quad \text{[Don]}\\
     &  \sum_{r=1}^{2} u_r y_{r,5} \leq 1.0  \quad \text{[Ed]}\\
     &  \sum_{r=1}^{2} u_r y_{r,6} \leq 1.0  \quad \text{[Fran]}\\
     &                  u_1, u_2\geq 0  
  \end{aligned}
\end{split}
(\#eq:Ch5NLP6PeopleHArdCode)
$$

Now we will extend this further to reflect that the constraint of being
less than or equal to 1.0 is the same for every manager. A key notation
to learn is the $\forall$ symbol which means to repeat for all possible
values of the index.

$$
\begin{split}
\begin{aligned}
    \text {max }    & \sum_{r=1}^{2} u_r y_{r,1}\\
    \text{s.t.:   } &  \sum_{r=1}^{2} u_r y_{r,j} 
                          \leq 1.0  \; \forall \; j\\
                    &  u_r \geq 0  \; \forall \; r
  \end{aligned}
\end{split}
$$

The $\forall\;j$ in the constraints indicates that the constraint should
be repeated for all possible values of j. In our case it would have a
constraint for $j=1$, (Alfred), then another constraint with $j=2$
(Barb), etc. The $\forall\;r$ in variable lower bounds is used to repeat
the non-negativity of each weighting variable (in this case, $u_1$ and
$u_2$).

The above formulation is nearly generalized. Let's replace the last
hard-coded items remaining. The first is that the formulation only
calculates the score for first manager. To generalize this, let's
replace the *1* in the objective function with *k*. This would allow us
to calculate the optimal score for any manager, *k*. Secondly, the
summations each assume that there are only two outputs. Let's replace
this with $N^Y$ to serve as a count of the number of outputs.

$$
\begin{split}
\begin{aligned}
    \text {max }    & \sum_{r=1}^{N^Y} u_r y_{r,k}\\
    \text{s.t.:   } &  \sum_{r=1}^{N^Y} u_r y_{r,j} 
                          \leq 1.0 \; \forall \; j\\
                    & u_r \geq 0  \; \forall \; r
  \end{aligned}
\end{split}
$$

This formulation is now a linear programming model that will find an
optimal score for manager *k* regardless of the number of outputs and
the number of other managers being considered. This approach could be
used in a variety of peer evaluation or self-evaluation applications
where each unit is given the opportunity to evaluate themselves relative
to their peers. This model will always give a person their highest
possible score. Note that this model does not incorporate the resources
used to achieve these outcomes. For that, we can use a traditional ratio
idea of output over input.

## The Ratio Model

Now, let's further extend this model to incorporate inputs. Not only are
the restaurant managers producing different outcomes, they are also
using different resources, say capital and labor. Let's denote capital
investment as $x_1$ and labor in full-time-equivalents as $x_2$.

We can frame the question of doing a two-input, two-output study in the
same way as we did for the two output example earlier. Instead of a
simple score, the manager can have a ratio of weighted outputs over
weighted inputs.

For the two-input, two-output case, the formulation is given below.

$$
\begin{split}
\begin{aligned}
    \text {max } & \frac{\sum_{r=1}^{2} u_r y_{r,k}} {\sum_{i=1}^{2} v_i x_{i,k} } \\
    \text{s.t.:   } & \frac{\sum_{r=1}^{2} u_r y_{r,j}} {\sum_{i=1}^{2} v_i x_{i,j} }
                          \leq 1 \; \forall \; j\\
                    & u_r, v_i\geq 0  \; \forall \; r,\; i
  \end{aligned}
\end{split}
$$

We have two sets of variables now. The first is $u_r$ which is the
weight on the *r*'th output. The second is $v_i$ which is the weight on
the *i*'th input.

For clarification, we make a small variation from traditional DEA
notation. Normally the number of inputs is *m* and the number of outputs
is *s*. To make the code more readable, I will use a convention of `NX`
instead of *m* to refer to the number of inputs (x's) and `NY` to be the
number of outputs (y's) instead of *s*. Also, *n* is used to denote the
number of Decision Making Units (DMUs) and therefore I'll use `ND` to
indicate that in the R code. In the mathematical formulations,
superscripts are used to differentiate the different N's resulting in
$N^X$, $N^Y$, and $N^D$ for `NX`, `NY`, and `ND` or *m*, *s*, and *n*
respectively.

$$
\begin{split}
\begin{aligned}
    \text {max } & \frac{\sum_{r=1}^{N^Y} u_r y_{r,k}} {\sum_{i=1}^{N^X} v_i x_{i,k} } \\
    \text{s.t.: } & \frac{\sum_{r=1}^{N^Y} u_r y_{r,j}} {\sum_{i=1}^{N^X} v_i x_{i,j} }
                          \leq 1 \; \forall \; j\\
                  &   u_r, \;v_i\geq 0  \; \forall \; r,i
  \end{aligned}
\end{split}
$$

This isn't a linear program because we are dividing functions of
variables by functions of variables. We need to make a few
transformations. First, we clear the denominator of each of the
constraints and collect the terms on resulting in the following
formulation.

$$
\begin{split}
 \begin{aligned}
    \text {max }   & \frac{\sum_{r=1}^{N^Y} u_r y_{r,k}} {\sum_{i=1}^{N^X} v_i x_{i,k} } \\
    \text{s.t.:  } & \sum_{r=1}^{N^Y} u_r y_{r,j} - \sum_{i=1}^{N^X} v_i x_{i,k} 
                          \leq 0 \; \forall \; j\\
                   & u_r, \; v_i\geq 0  \; \forall \; r,i
  \end{aligned}
 \end{split}
$$

Alas, it is still a \index{Nonlinear function} nonlinear objective
function. Fortunately, we can take advantage of the structure of the
problem. There are an infinite number of possible combinations of
numerators and denominators that can give the same ratio. For example,
if a unit were to receive a 90% efficiency score, this could be with a
numerator of nine and a denominator of 10 or a numerator of 90 and a
denominator of 100. Each of these pairs satisfies the constraints and is
therefore equally valid. One clever trick solves the non-linearity by
setting the denominator to be a constant. The next step is to select a
normalizing value for the objective function. Let's set the denominator
equal to one. In this case, we simply add a constraint,
$\sum_{i=1}^{N^X} v_i x_{i,k}$, to the linear program and replace the
denominator in the objective function to one means that we can rewrite
the objective function using only the numerator.

$$
\begin{split}
 \begin{aligned}
    \text {max   } & \sum_{r=1}^{N^Y} u_r y_{r,k} \\
    \text{s.t.:  } & \sum_{i=1}^{N^X} v_i x_{i,k} = 1 \\
    & \sum_{r=1}^{N^Y} u_r y_{r,j} - \sum_{i=1}^{N^X} v_i x_{i,j} 
                          \leq 0 \; \forall \; j\\
                   &  u_r, \; v_i\geq 0  \; \forall \; r,i
  \end{aligned}
 \end{split}
$$

This linear program can then be implemented in an linear programming
system.

## Creating the LP - The Algebraic Approach

We will implement this using the `ompr` package again, as we did in
chapters 2 and 4.

We're going to use our data from earlier of grocery stores.

```{r, echo=FALSE}
x <- matrix(c(10, 20, 30, 50, 22, 29,
              75, 100, 220, 480, 290, 210), byrow=FALSE, ncol=2, 
            dimnames=list(LETTERS[1:6],c("x1","x2")))
y <- matrix(c(75,100,300,400, 280, 120),ncol=1,
            dimnames=list(LETTERS[1:6],"y"))

storenames<- c("Al\'s Pantry", "Bob\'s Mill", 
               "Trader Carrie\'s", "Dilbertson\'s",
               "Ed\'s Eggshop", "Flo\'s Farmacopia")
temp<-cbind(storenames,x,y)
colnames(temp)<-c("Store Name", '"Employees \n (x1)"', 
                  "Floor Space (x2)", "Sales (y)")
kbl (temp, booktabs=TRUE, escape=FALSE)

ND <- nrow(x); NX <- ncol(x); NY <- ncol(y) # Define data size
xdata       <- x [1:ND,]   # Call it xdata
dim (xdata) <- c (ND,NX)   # structure data correctly
ydata       <- y [1:ND,]
dim (ydata) <- c (ND,NY)
```

Remember the inputs are named as "x1" and"x2" to represent the staff
(full-time equivalent employees) and the store size ($m^2$)
respectively. The output is sales and is named "y1". As usual, I'm
copying the source data into objects, `xdata` and `ydata`, for actual
operation.

```{r Structure-Results}

mm1.efficiency <- matrix(rep(-1.0, ND), nrow=ND)
mm1.lambda     <- matrix(rep(-1.0, ND^2), nrow=ND)
mm1.vweight    <- matrix(rep(-1.0, ND*NX), nrow=ND) 
mm1.uweight    <- matrix(rep(-1.0, ND*NY), nrow=ND) 
mm1.xslack     <- matrix(rep(-1.0, ND*NX), nrow=ND) 
mm1.yslack     <- matrix(rep(-1.0, ND*NY), nrow=ND) 

mm1names <- TRA::DEAnames(NX, NY, ND)
```

We are now ready to do the analysis. In chapter 2 we used the
envelopment model to examine this case. Now we will use the multiplier
model.

```{r eval=TRUE}
for (k in 1:ND) {

  result <- MIPModel() |>
  add_variable(vweight[i], i = 1:NX, type = "continuous", lb = 0)   |>
  add_variable(uweight[r], r = 1:NY, type = "continuous", lb = 0)   |>
  set_objective(sum_expr(uweight[r] * ydata[k,r], r = 1:NY), "max") |>
  add_constraint(sum_expr(vweight[i] * xdata[k,i], i = 1:NX) == 1)  |>
  add_constraint((sum_expr(uweight[r] * ydata[j,r], r = 1:NY)-
                    sum_expr(vweight[i] * xdata[j,i], i = 1:NX)) 
                 <= 0, j = 1:ND)
  result

  result <- solve_model(result, with_ROI(solver = "glpk", 
                                         verbose = FALSE))
  mm1.efficiency[k] <- objective_value (result) 

  # Get the output weights
  tempvweight <- get_solution(result, vweight[i])
  mm1.vweight[k,] <- tempvweight[,3]

  # Get the input weights
  tempuweight <- get_solution(result, uweight[i])
  mm1.uweight[k,] <- tempuweight[,3]
  templambda <- as.matrix(get_row_duals(result))
  mm1.lambda[k,] <- templambda [-1]
  # Drops the first dual value since that 
  #    constraint since is for setting denom=1
 }

mm1.combined <- cbind (mm1.efficiency, 
                       mm1.vweight, mm1.uweight,
                       mm1.lambda)
```

```{r}
kbl (mm1.combined , booktabs=T, escape=F, digits=4,
     col.names=c("$\\theta^{CRS}$", mm1names$VnamesLX,
                 mm1names$UnamesLX,
                 mm1names$LambdanamesbyletterLX),
     caption =c("Multiplier Model Results from Six Store Example")) |>
  kable_styling (latex_options = c("HOLD_position"))
```

The weights from the multiplier model often have multiple alternative
solutions-in other words, different values for
\index{Decision variables} decision variables that give the same
objective function value. In optimization, this situation is called
\index{Multiple optima} multiple optima. In DEA this means that the
weights may be different depending upon the particular settings used for
solving the linear program. The result is that you should be careful in
over interpreting the weights as well as in how you interpret results
when trying to reproduce analyses.

\index{Multiplier model|)}

### Envelopment Multiplier Model Duality Relationship

\index{Dual values} Duality is a fundamental characteristic of linear
programming. Every linear program has a dual program that is a bit of
mirror image of the original. The multiplier model and the
\index{Envelopment model} envelopment models are fundamentally connected
by duality. The multiplier model of DEA is considered the dual of the
envelopment model and the reverse is also true.

Consider the basic structure of an LP as discussed in Chapter 3.

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

A linear program can be converted to its dual by:

-   Creating a variable for every constraint in the original LP
-   Creating a constraint for every variable in the original LP
-   Transposing the *A* matrix, sometimes denoted as $A^{\tau}$
-   Using the *B* vector as the objective function coefficients
-   Using the *C* vector as the right hand side values of the
    constraints.
-   Changing the objective function from a *min* to a *max*
-   Reverse the directions of the constraint inequalities

$$
\begin{split}
 \begin{aligned}
    \text{max  }   & BY \\
    \text{s.t.:  } & A^{\tau} Y \leq C\\ 
                       & Y \geq 0  \\
  \end{aligned}
   \end{split}
$$

We won't dwell further on this transformation and will leave it to the
reader to explore further but this is why we can extract dual values of
the envelopment model to get the multiplier weights from the multiplier
model. Also, it is why we can extract the envelopment's $\lambda$ values
from the multiplier model's results.

## Weight Restrictions in the Multiplier Model

\index{Weight restrictions|(}

$$
\begin{split}
 \begin{aligned}
    \text {max   } & \sum_{r=1}^{N^Y} u_r y_{r,k} \\
    \text{s.t.:  } & \sum_{i=1}^{N^X} v_i x_{i,k} = 1 \\
    & \sum_{r=1}^{N^Y} u_r y_{r,j} - \sum_{i=1}^{N^X} v_i x_{i,j} 
                          \leq 0 \; \forall \; j\\
                   &  u_r, \; v_i\geq 0  \; \forall \; r,i
  \end{aligned}
 \end{split}
$$

Weight restrictions can take a number of forms and are easy to visualize in the multiplier model.  

For example, we could simply indicate that the second output is more valuable or challenging to produce than the first.  We could then add a constraint that $u_1\leq u_2$.  If we had more specific information, we might even put a range on that relationship such as the weight of output must be at least 10% of the first output and no more than double that of the first output.  This would be represented as $0.1\leq \frac{u_2}{u_1} \leq 2.0$ which is readily linearized by muliplying out the $u_1$ term.

Another approach is to say that the component scores for each unit must follow a certain relationship. We might say that amount of the weighted output from the second output must be at least 10% of that of the first output but no more than double that of the first weighted output or $0.1\leq \frac{u_2y_{2,k}}{u_1y_{1,k}} \leq 2.0$.

Great care must be used in applying weight restrictions. As it removes a fundamental characteristic of weighting flexibility from DEA. Also certain weight restriction approaches may risk making the analysis of certain units infeasible. 

Let's build a model with and without the weight restrictions that an employee, $x_1$ is  considered at least twice as expensive as a unit of floor space, $x_2$ or $v_1\geq v_2$.    


```{r Structure-MM2-Wrs-Results}

mm2.efficiency <- matrix(rep(-1.0, ND), nrow=ND)
mm2.lambda     <- matrix(rep(-1.0, ND^2), nrow=ND)
mm2.vweight    <- matrix(rep(-1.0, ND*NX), nrow=ND) 
mm2.uweight    <- matrix(rep(-1.0, ND*NY), nrow=ND) 
mm2.xslack     <- matrix(rep(-1.0, ND*NX), nrow=ND) 
mm2.yslack     <- matrix(rep(-1.0, ND*NY), nrow=ND) 

mm2names <- TRA::DEAnames(NX, NY, ND)
```

We are now ready to do the analysis. In chapter 2 we used the
envelopment model to examine this case. Now we will use the multiplier
model.

```{r eval=TRUE}
for (k in 1:ND) {

  mm2result <- MIPModel ()                                           |>
  add_variable (vweight[i], i = 1:NX, type = "continuous", lb = 0)   |>
  add_variable (uweight[r], r = 1:NY, type = "continuous", lb = 0)   |>
  set_objective (sum_expr(uweight[r] * ydata[k,r], r = 1:NY), "max") |>
  add_constraint (sum_expr(vweight[i] * xdata[k,i], i = 1:NX) == 1)  |>
  add_constraint ((sum_expr(uweight[r] * ydata[j,r], r = 1:NY)-
                    sum_expr(vweight[i] * xdata[j,i], i = 1:NX)) 
                 <= 0, j = 1:ND)                                     |>         
  add_constraint (vweight[1]>=2*vweight[2])  # Weight Restriction
  
  mm2result

  mm2result <- solve_model(mm2result, with_ROI(solver = "glpk", 
                                         verbose = FALSE))
  mm2.efficiency[k] <- objective_value (mm2result) 

  # Get the output weights
  tempvweight <- get_solution(mm2result, vweight[i])
  mm2.vweight[k,] <- tempvweight[,3]

  # Get the input weights
  tempuweight <- get_solution(mm2result, uweight[i])
  mm2.uweight[k,] <- tempuweight[,3]
  templambda <- as.matrix(get_row_duals(mm2result))
  mm2.lambda[k,] <- templambda [2:7]
  # Drops the first and last row duals as those constraints
  #   since the first is for setting demon=1, last is for 
  #   weight restrictions.
 }

mm2.combined <- cbind (mm1.efficiency, mm2.vweight, mm1.uweight,
                       mm2.efficiency, mm2.vweight, mm2.uweight)
```

```{r, include=FALSE}
kbl (mm2.combined , booktabs=T, escape=F, digits=5,
     col.names=c("$\\theta^{CRS}$", mm1names$VnamesLX,
                 mm1names$UnamesLX,
                 "$\\theta^{CRS}$", mm1names$VnamesLX,
                 mm1names$UnamesLX),
     caption =c("Multiplier Model Results from Grocery Store Example")) |>
  add_header_above(c("Original Results" = 4, 
                     "With Weight Restrictions" = 4))                   |>  
  kable_styling (latex_options = c("HOLD_position","scale_down"))
```

In this case, the impact of weight restrictions is relatively mild as it decreased the efficiency of store B modestly and store F by only a small amount. A key insight is that weight restrictions will not improve any unit's efficiency score and can only decrease it or at best leave it unchanged. Also, weights are dependent on the unit of measure of the input or output. 

Other approaches can be used for restrictions including the so-called cone ratio model for use in the envelopment model.  More generally, ordinal weight restrictions can be employed as a data pre-processing step and thereby be employed in many DEA models. This will be demonstrated with an extended example in Chapter 8.

\index{Weight restrictions|)}

## Output-Oriented Multiplier Model

Recall that the input-oriented multiplier model had a ratio of weighted
outputs over weighted inputs.  \index{Output-oriented|(}

$$
\begin{split}
\begin{aligned}
    \text {max } & \frac{\sum_{r=1}^{N^Y} u_r y_{r,k}} {\sum_{i=1}^{N^X} v_i x_{i,k} } \\
    \text{s.t.: } & \frac{\sum_{r=1}^{N^Y} u_r y_{r,j}} {\sum_{i=1}^{N^X} v_i x_{i,j} }
                          \leq 1 \; \forall \; j\\
                  &   u_r, \;v_i\geq 0  \; \forall \; r,i
  \end{aligned}
\end{split}
$$

The output-oriented model is then simply the reciprocal of this along
with changing the *max* to a *min* and changing the direction of the
inequality constraints.

$$
\begin{split}
\begin{aligned}
    \text {min } &  \frac{\sum_{i=1}^{N^X} v_i x_{i,k}} {\sum_{r=1}^{N^Y} u_r y_{r,k} } \\
    \text{s.t.: } & \frac{\sum_{i=1}^{N^X} v_i x_{i,j}} {\sum_{r=1}^{N^Y} u_r y_{r,j} }
                          \geq 1 \; \forall \; j\\
                  &   u_r, \;v_i\geq 0  \; \forall \; r,i
  \end{aligned}
\end{split}
$$

The linearization follows in a similar process to the input-oriented model.

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

\index{Output-oriented|)}

## Variable Returns to Scale

\index{Variable returns to scale|(} 

Formulate and implement a model Implement model. By duality, adding a
constraint to the primal model, is equivalent to adding a variable for
in dual. If we consider the envelopment model to be the primal and the
multiplier to be its dual, then we are adding a variable. We'll term
this variable $\mu_0$ to reflect the returns to scale associated with
this DMU. If we are implementing a VRS model, then the equality
constraint indicates that $\mu_0$ is an unrestricted (or free) variable.

The input-oriented multiplier model is shown below.

$$
\begin{split}
 \begin{aligned}
    \text {max   } & \sum_{r=1}^{N^Y} u_r y_{r,k} +\mu_k\\
    \text{s.t.:  } & \sum_{i=1}^{N^X} v_i x_{i,k} = 1 \\
    & \sum_{r=1}^{N^Y} u_r y_{r,j} - \sum_{i=1}^{N^X} v_i x_{i,j} + \mu_k
                          \leq 0 \; \forall \; j\\
                   &  u_r, \; v_i\geq 0  \; \forall \; r,i\\
                   & \mu_k \;\text{ free}
  \end{aligned}
 \end{split}
$$

The output-oriented VRS multiplier model is the following.  The returns to scale term is now denoted as *nu* or $\nu_k$ for unit *k*.  Note that it is common in some DEA multiplier models to use the Greek symbols nu ($\nu$) and mu ($\mu$) as the input and output multipliers instead of *v* and *u* and then use the alternate as the returns to scale terms.  In this book I use *u* and *v* consistently  to make it more accessible for reader and limit confusion between what appear to be visually similar terms *v* and $\nu$ to this this section.

$$ 
\begin{split} 
 \begin{aligned}     
    \text {min }  &  \sum_{i=1}^{N^X} v_i x_{i,k}+\nu_k\\
    \text{s.t.: } &  \sum_{r=1}^{N^Y} u_r y_{r,k} =1 \\
                  &  \sum_{i=1}^{N^X} v_i x_{i,j} - \sum_{r=1}^{N^Y} u_r y_{r,j} + \nu_k \geq 0 \; \forall \; j\\
                  &  u_r, \;v_i\geq 0  \; \forall \; r,i\\
                  & \nu_k \;\text{ free}
  \end{aligned} 
 \end{split} 
$$

### Numerical Example Multiplier 

```{r Structure-VRS-MM-Results}

mm3.efficiency <- matrix(rep(-1.0, ND), nrow=ND)
mm3.lambda     <- matrix(rep(-1.0, ND^2), nrow=ND)
mm3.vweight    <- matrix(rep(-1.0, ND*NX), nrow=ND) 
mm3.uweight    <- matrix(rep(-1.0, ND*NY), nrow=ND) 
mm3.xslack     <- matrix(rep(-1.0, ND*NX), nrow=ND) 
mm3.yslack     <- matrix(rep(-1.0, ND*NY), nrow=ND) 
mm3.mu         <- matrix(rep(-1.0, ND), nrow=ND)
mm3names <- TRA::DEAnames(NX, NY, ND)
```

```{r eval=TRUE}
for (k in 1:ND) {

  mm3result <- MIPModel ()                                           |>
  add_variable (vweight[i], i = 1:NX, type = "continuous", lb = 0)   |>
  add_variable (uweight[r], r = 1:NY, type = "continuous", lb = 0)   |>
  add_variable (mu, type = "continuous")                             |>
  set_objective (sum_expr(uweight[r] * ydata[k,r], r = 1:NY)+mu, 
                 "max")                                              |>
  add_constraint (sum_expr(vweight[i] * xdata[k,i], i = 1:NX) == 1)  |>
  add_constraint ((sum_expr(uweight[r] * ydata[j,r], r = 1:NY)-
                    sum_expr(vweight[i] * xdata[j,i], i = 1:NX))+mu 
                 <= 0, j = 1:ND)                                             

  mm3result

  mm3result <- solve_model(mm3result, with_ROI(solver = "glpk", 
                                         verbose = FALSE))
  mm3.efficiency[k] <- objective_value (mm3result) 

  # Get the output weights
  tempvweight <- get_solution(mm3result, vweight[i])
  mm3.vweight[k,] <- tempvweight[,3]
  mm3.mu[k] <- get_solution(mm3result,mu)
  # Get the input weights
  tempuweight <- get_solution(mm3result, uweight[i])
  mm3.uweight[k,] <- tempuweight[,3]
  templambda <- as.matrix(get_row_duals(mm2result))
  mm3.lambda[k,] <- templambda [2:7]
  # Drops the first and last row duals as those contstraints
  #   since the first is for setting demon=1, last is for 
  #   weight restrictions.
 }

mm3.combined <- cbind (mm1.efficiency, mm2.vweight, mm1.uweight,
                       mm3.efficiency, mm3.mu, mm3.vweight, mm3.uweight)
```

```{r echo=FALSE}
kbl (mm3.combined , booktabs=T, escape=F, digits=5,
     col.names=c("$\\theta^{CRS}$", mm1names$VnamesLX,
                 mm1names$UnamesLX,
                 "$\\theta^{VRS}$", "$\\mu_k$", mm1names$VnamesLX,
                 mm1names$UnamesLX),
     caption =c("Multiplier Model Results from Grocery Store Example")) |>
  add_header_above(c("Original Results" = 4, 
                     "With Weight Restrictions" = 5))                   |>  
  kable_styling (latex_options = c("HOLD_position","scale_down"))
```

These results match those from the envelopment model's six grocery store case in chapter 2.  The scale factor, $mu_k$ indicates whether the unit is operating at most productive scale size ($\mu_k=0$), at too small a size where growing could result in greater returns ($\mu_k<0$) or at too large a size where it is suffering from diseconomies of scale and may be better off shrinking ($\mu_k>0$).Of course these results are all relative to the data available. 

### Other Returns to Scale Models

Restrictions on $\nu$ or $\mu$ can be made to multiplier model to allow for other returns scale assumptions. \index{Increasing returns to scale} \index{Decreasing returns to scale} Recall that from Chapter 2, Increasing returns to scale is more formally described as non-decreasing returns to scale and that Decreasing returns to scale is more formally described as non-increasing returns to scale.

```{r, echo=FALSE}
RTS_Table  <- matrix(c("CRS", "$\\text{None or }\\sum_{j=1}^{n} \\lambda_j  \\geq 0$", "$\\mu_k =0$", "$\\nu_k =0$", 
                       "VRS", "$\\sum_{j=1}^{n} \\lambda_j = 1$", "$\\mu_k \\text{ free}$", "$\\nu_k \\text{ free}$",
                       "DRS/NIRS", "$\\sum_{j=1}^{n} \\lambda_j \\leq 1$", "$\\mu_k \\leq 0$", "$\\nu_k \\geq 0$", 
                       "IRS/NDRS", "$\\sum_{j=1}^{n} \\lambda_j \\geq 1$","$\\mu_k \\geq 0$", "$\\nu_k \\leq 0$"),
                     byrow=TRUE, ncol=4)
kbl (RTS_Table , booktabs=T, escape=F, 
     col.names=c("Returns to Scale", "Envelopment", "Input-Oriented", "Output-Oriented"),
     caption =c("Returns to Scale Constraints for Envelopment and Multiplier Bounds"))  |>
  kable_styling (latex_options = c("HOLD_position","scale_down"))
```

\index{Variable returns to scale|)}

## Slacks in the Multiplier Model

Just as the envelopment model can be extended to allow for identifying weakly efficient units, *i.e.* units that are radially efficient but have a positive nonradial slack, we can also extend the multiplier model. Recall that we used $\epsilon$ as a non-Archimedean infinitesimal to indicate that the priority was primarily to minimize $\theta$. \index{non-Archimedean infinitesimal} \index{Slack|(}

$$
\begin{split}
 \begin{aligned}
    \text{minimize  }   & \theta - \epsilon ( \sum_{i} s^x_i + \sum_{r} s^y_r)\\
    \text{subject to } & \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 \; j,i,r
  \end{aligned}
 \end{split}
$$

The corresponding change in the multiplier model is to put a lower bound of $\epsilon$ on the input and output weights.  


$$
\begin{split}
 \begin{aligned}
    \text {max   } & \sum_{r=1}^{N^Y} u_r y_{r,k} \\
    \text{s.t.:  } & \sum_{i=1}^{N^X} v_i x_{i,k} = 1 \\
    & \sum_{r=1}^{N^Y} u_r y_{r,j} - \sum_{i=1}^{N^X} v_i x_{i,j} 
                          \leq 0 \; \forall \; j\\
                   &  u_r, \; v_i\geq \epsilon  \; \forall \; r,i
  \end{aligned}
 \end{split}
$$

Unfortunately finite approximations to $\epsilon$ will cause problems in the multiplier model just like in the envelopment model. Again, a second phase approach could be used to address this issue.  

$$
\text{Phase 1:          }
\begin{split}
 \begin{aligned}
    \text {max   } & \theta_k=\sum_{r=1}^{N^Y} u_r y_{r,k} \\
    \text{s.t.:  } & \sum_{i=1}^{N^X} v_i x_{i,k} = 1 \\
    & \sum_{r=1}^{N^Y} u_r y_{r,j} - \sum_{i=1}^{N^X} v_i x_{i,j} 
                          \leq 0 \; \forall \; j\\
                   &  u_r, \; v_i\geq 0  \; \forall \; r,i
  \end{aligned}
 \end{split}
$$

$$
\text{Phase 2:          }
\begin{split}
 \begin{aligned}
    \text {min   } & Q_k \\
    \text{s.t.:  } & \sum_{i=1}^{N^X} v_i x_{i,k} = 1 \\
    & \sum_{r=1}^{N^Y} u_r y_{r,j} - \sum_{i=1}^{N^X} v_i x_{i,j} 
                          \leq 0 \; \forall \; j\\
                   & \sum_{r=1}^{N^Y} u_r y_{r,k} =\theta_k^* &\text{Keep fixed from Phase 1}\\
                   &  u_r, \; v_i\geq Q_k  \; \forall \; r,i &\text{Increase smallest weight if possible}\\
                   & Q_k \geq 0
  \end{aligned}
 \end{split}
$$

If $\theta_k^*=1$ and $Q_k>0$ then unit is *k* is weakly effficient.  \index{Slack|)}

## Super-Efficiency in the Multiplier Model

The multiplier model can be modified to account for super-efficiency in a similar way to that of the envelopment model.

In the ratio model, we don't require that unit *k* limits its own score to 1.0.  This means that we simply eliminate the ratio constraint where $j=k$. \index{Super-efficiency}

$$
\begin{split}
 \begin{aligned}
    \text {max }   & \frac{\sum_{r=1}^{N^Y} u_r y_{r,k}} {\sum_{i=1}^{N^X} v_i x_{i,k} } \\
    \text{s.t.:  } & \frac{\sum_{r=1}^{N^Y} u_r y_{r,j}} { \sum_{i=1}^{N^X} v_i x_{i,k} }
                          \leq 1 \; \forall \; j, \; j\neq k\\
                   & u_r, \; v_i\geq 0  \; \forall \; r,i
  \end{aligned}
 \end{split}
$$

We can then make the same change to the corresponding linear program version of the ratio model.

$$ 
\begin{split} 
 \begin{aligned}     
    \text {min }  &  \sum_{i=1}^{N^X} v_i x_{i,k}+\nu_k\\
    \text{s.t.: } &  \sum_{r=1}^{N^Y} u_r y_{r,k} =1 \\
                  &  \sum_{i=1}^{N^X} v_i x_{i,j} - \sum_{r=1}^{N^Y} u_r y_{r,j} + \nu_k \geq 0 \; \forall \; j, \; j\neq k\\
                  &  u_r, \;v_i\geq 0  \; \forall \; r,i\\
                  & \nu_k \;\text{ free}
  \end{aligned} 
 \end{split} 
$$