06-Cross.Rmd

---
output:
  html_document: default
  pdf_document: default
---
# Cross-Efficiency-DRAFT 

Introduction
-----------------------------

##The Ratio Model

Put yourself as a competitor trying to argue that you are the best in converting inputs into outputs among a set of other units.  You have data on what the competitors' inputs and outputs.  You can 

The Linear Programs for DEA
-----------------------------

On the other hand, what if it allowed for blending of units.  There are a few assumptions that we could make.  Let's start by saying that we can compare any particular products by rescaling (up or down) any other product as well as any combination of units.  

We'll start by creating a mathematical framework. Can you find a combination of units that produces at least as much output using less input?  Let's define the proportion of input needed as $\theta$.  A value of $\theta=1$ then means no input reduction can be found in order to produce that unit's level of output.  The blend of other units is described by a vector $\lambda$.  Another way to denote this is $\lambda_j$ is the specific amount of a unit _j_ used in setting the target for for performance for unit _k_.  Similarly, $x_j$ is the amount of input used by unit _j_ and $y_j$ is the amount of output produced by unit _j_.  

This can be easily expanded to the multiple input and multiple output case by defining $x_i,j$ to  be the amount of the _i_'th input used by unit _j_ and $y_r,j$ to be the amount of the _r_'th output produced by unit _j_.  For simplicity, this example will focus on the one input and one output case rather than the _m_ input and _s_ output case but the R code explicitly allows for $m,s>1$.  To make the code more readable, I will use which corresponds to _NX_ instead of _m_ to refer to the number of inputs (x's) and _NY_ to be the number of ouputs (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.  

We have two important 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.  

The multiplier model can be thought of as finding a weighting scheme for outputs over inputs that give you the best possible score while giving no one better than _1.0._



$$
\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{subject to } & \frac{\sum_{r=1}^{N^Y} u_r y_{r,k}} {\sum_{i=1}^{N^X} v_i x_{i,k} }
                          \leq 1 \forall \; j\\
                       u_r, v_i\geq 0  \; \forall \; r,i
  \end{aligned}
\end{split}
(\#eq:LPCCRIOM-Ratio-1a)
$$

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 consraints 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{subject to } & \sum_{r=1}^{N^Y} u_r y_{r,k} - \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}
(\#eq:LPCCRIOM-Ratio-2a)
$$

Now we will convert the problem input and output constraints from inequalities into equalities by explicitly defining slack variables.  

There are an infinite number of possible combinations of numerators and denominators that can give the same ratio.  The next step is to select 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.  


$$
\begin{split}
 \begin{aligned}
    \text {max }  \sum_{r=1}^{N^Y} u_r y_{r,k} \\
    \text{subject to } & \sum_{i=1}^{N^X} v_i x_{i,k} \\
    & \sum_{r=1}^{N^Y} u_r y_{r,k} - \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}
(\#eq:LPCCRIOM-NoSlack)  
$$

## Creating the LP - The Algebraic Approach

We will implement this using the ompr package again.

We're going to use our data from earlier.  For this example, we will use the dataset from Kenneth Baker's third edition of _Optimization Modeling with Spreadsheets_, pages 175-178, Example 5.3 titled "Hope Valley Health Care Association." In this case, a health care organization wants to benchmark six nursing homes against each other.

```{r multiplier_model }

library(pander, quietly=TRUE)           # For better tables
library(dplyr, quietly=TRUE)            # For data structure manipulation
library(ROI, quietly=TRUE)              # R Optimization Interface 
library(ROI.plugin.glpk, quietly=TRUE)  # Connection to glpk as solver
library(ompr, quietly=TRUE)             # Optimization Modeling using R
library(ompr.roi, quietly=TRUE)         # Connective tissue

  XBaker1 <- matrix(c(150, 400, 320, 520, 350, 320, .2, 0.7, 1.2, 2.0, 1.2, 0.7),
                  ncol=2,dimnames=list(LETTERS[1:6],c("x1", "x2")))

  YBaker1 <- matrix(c(14000, 14000, 42000, 28000, 19000, 14000, 3500, 21000, 10500, 
                    42000, 25000, 15000),
                  ncol=2,dimnames=list(LETTERS[1:6],c("y1", "y2")))

ND <- nrow(XBaker1); NX <- ncol(XBaker1); NY <- ncol(YBaker1); # Define data size

xdata      <-XBaker1 [1:ND,]  # Call it xdata
dim(xdata) <-c(ND,NX)  # structure data correctly
ydata      <-YBaker1[1:ND,]
dim(ydata) <-c(ND,NY)
  
```

Remember the inputs are hard coded as "x1" and"x2" to represent the _staff hours per day_ and the _supplies per day_ respectively.  The two outputs of _reimbursed patient-days_ and _privately paid patient-days_ are named "y1" and "y2".  

```{r}
YBaker1 <- matrix(c(14000, 14000, 42000, 28000, 19000, 14000, 3500, 21000, 10500, 
                    42000, 25000, 15000),
                  ncol=2,dimnames=list(LETTERS[1:6],c("y1", "y2")))

ND <- nrow(XBaker1); NX <- ncol(XBaker1); NY <- ncol(YBaker1); # Define data size

```

Note that I'm naming the data sets based on their origin and then loading them into xdata and ydata for actual operation.

```{r Structure_Results}
# Need to remember to restructure the results matrices.

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) 
```

We are now ready to do the analysis.  In Baker, the analysis was done using the multiplier model.  In chapter 2 we used the envelopment model to examine this case.  Now we will use the multiplier model using the DEAMultiplerModel package by Aurobindh Kalathil Puthanpura.  

```{r}
library (MultiplierDEA)

# Example from Kenneth R. Baker: Optimization Modeling with Spreadsheets, 
#    Third Edition, p. 176, John Wiley & Sons, Inc.
dmu <- c("A", "B", "C", "D", "E", "F")
x <- data.frame(c(150,400,320,520,350,320),c(0.2,0.7,1.2,2.0,1.2,0.7))
rownames(x) <- dmu
colnames(x)[1] <- c("StartHours")
colnames(x)[2] <- c("Supplies")
y <- data.frame(c(14,14,42,28,19,14),c(3.5,21,10.5,42,25,15))
rownames(y) <- dmu
colnames(y)[1] <- c("Reimbursed")
colnames(y)[2] <- c("Private")
# Calculate the efficiency score
result <- DeaMultiplierModel(x,y,"crs", "input")
# Examine the efficiency score for DMUs
pander(result$Efficiency, caption="CCR Efficiency Scores")
```

The efficiency scores match the earlier results and that of Baker.  

The package also includes a cross-efficiency calculation function so as to find the cross-efficiencies given the multiplier model.  Cross-efficiency calculated in this way tends to suffer from frequent cases of multiple optima.  Later in this chapter, we will explore how to handle these issues. For now, let's apply cross-efficiency to the same data set.  

```{r}
# Example from Kenneth R. Baker: Optimization Modeling with Spreadsheets, 
#    Third Edition, p. 176, John Wiley & Sons, Inc.
dmu <- c("A", "B", "C", "D", "E", "F")
x <- data.frame(c(150,400,320,520,350,320),c(0.2,0.7,1.2,2.0,1.2,0.7))
rownames(x) <- dmu
colnames(x)[1] <- c("StartHours")
colnames(x)[2] <- c("Supplies")
y <- data.frame(c(14,14,42,28,19,14),c(3.5,21,10.5,42,25,15))
rownames(y) <- dmu
colnames(y)[1] <- c("Reimbursed")
colnames(y)[2] <- c("Private")
# Calculate the efficiency score
result <- CrossEfficiency(x,y,"crs", "input")
# Examine the cross efficiency score for DMUs
pander(result$ce_ave, caption="Cross-Efficiency Scores")

```

As might be expected, the cross-efficiency values are all significantly lower than the original efficiency scores.  In cross-efficiency, no nursing home is able to pick an extreme and unique weighting scheme to their advantage.  It also tends to break the ties at 1.0 that frequently occur in DEA so it usually gives a unique ranking.  These features come at a significant cost though which will be discussed in more detail in the future.

## Implementing Cross-Efficiency

Let's revisit the Baker model again.  

```{r Baker_cross, eval=TRUE}

multiplierIO <- function (x,y) 
  { 

  ND <- nrow(x); NX <- ncol(y); NY <- ncol(y); # Define data size
  
  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) 

    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] * y[k,r], r = 1:NY), "max") %>%
    add_constraint(sum_expr(vweight[i] * x[k,i], i = 1:NX) == 1) %>%
    add_constraint((sum_expr(uweight[r] * y[j,r], r = 1:NY)-
                    sum_expr(vweight[i] * x[j,i], i = 1:NX)) 
                  <= 0, j = 1:ND)
    result

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

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

   # Get the weights- Input weights
    tempuweight <- get_solution(result, uweight[i])
    results.uweight[k,] <- tempuweight[,3]

  } # End of for k loop
  
  resultlist <- list(efficiency=results.efficiency, 
                     vweight=results.vweight, uweight=results.uweight)
 
  return(resultlist)  

} # End of function  

resfunc <- multiplierIO(x,y)
DMUnames <- list(c(LETTERS[1:ND]))
Vnames<- lapply(list(rep("v",NX)),paste0,1:NX)   # Input weight names: v1, ...
Unames<- lapply(list(rep("u",NY)),paste0,1:NY)   # Output weight names: u1, ...

dimnames(resfunc$efficiency) <- c(DMUnames,"CCR-IO")
dimnames(resfunc$vweight)    <- c(DMUnames,Vnames)
dimnames(resfunc$uweight)    <- c(DMUnames,Unames)

pander (cbind(resfunc$efficiency,resfunc$vweight,resfunc$uweight), 
        caption="Efficiency Scores and Weights")
```

You can think of the cross-efficiency by going back to revisit the original ratio model of the multiplier model.  Let's extend the definition of the input and output weights to reflect for which unit the weights were determined.  For example, $u_{r,k}$ is then the weight on output r when analyzed from the perspetive of DMU k and a similar interpretation applies to input weight.  

$$
\begin{split}
\begin{aligned}
    \text {max } \frac{\sum_{r=1}^{N^Y} u_{r,k} y_{r,k}} {\sum_{i=1}^{N^X} v_i,k x_{i,k} } \\
    \text{subject to } & \theta_{j,k}=\frac{\sum_{r=1}^{N^Y} u_{r,k} y_{r,k}} {\sum_{i=1}^{N^X} v_{i,k} x_{i,k} }
                          \leq 1 \forall \; j\\
                       u_r, v_i\geq 0  \; \forall \; r,i
  \end{aligned}
\end{split}
(\#eq:LPCCRIO-Cross1)
$$

Now the cross-efficiency for a unit j can be calculated as the following.  I refer to the value $\theta_{k,j}$ as the cross-evaluation score for unit _j_ from the perspective of the evaluation of unit _k_.

$$
\begin{split}
 \begin{aligned}
    \frac{\sum_{r=1}^{N^Y} u_{r,k} y_{r,k}} {\sum_{i=1}^{N^X} v_i,k x_{i,k} } \\
    \text{subject to } & \theta_{k,j} = \frac{\sum_{r=1}^
                 {N^Y} u_{r,k} y_{r,k}} {\sum_{i=1}^{N^X} v_{i,k} x_{i,k} }
                          \leq 1 \forall \; j\\
                       u_r, v_i\geq 0  \; \forall \; r,i
  \end{aligned}
   \end{split}
  (\#eq:Cross-Calc-Peer)
$$

Now, the cross-efficiency score for unit _j_ is simply the average of all of the cross-evaluation scores given to _j_ by the other units when they do their analysis.

$$
\begin{split}
 \begin{aligned}
    CE_j= \frac {\sum_{k=1}^{N^D}\theta_{k,j}} {N^D}
        = \frac {\sum_{k=1}^{N^D} {
                        \frac {\sum_{r=1}^{N^Y} u_{r,k} y_{r,j} } 
                       {\sum_{i=1}^{N^X} v_{i,k} x_{i,j} } }    }
                {N^D}
  \end{aligned}
\end{split}
(\#eq:Calc-Cross)
$$



```{r}
results.crosseval  <- matrix(rep(-1.0, ND*ND), nrow=ND,ncol=ND) 
results.crosseff   <- matrix(rep(-1.0, ND)) 

for (k in 1:ND) {
   for (j in 1:ND) {results.crosseval[k,j]<-sum(resfunc$uweight[k,]*y[j,])/sum(resfunc$vweight[k,]*x[j,])
                   }
}
for (k in 1:ND) {
  results.crosseff[k]<-mean(results.crosseval[,k])
}

dimnames(results.crosseval) <- c(DMUnames,DMUnames)
dimnames(results.crosseff) <- c(DMUnames)
pander(results.crosseval, caption = "Cross-evaluation Matrix")

#for (j in 1:ND) {results.crosseff[j]<-sum(results.crosseval[,j])/ND}
pander(results.crosseff, caption = "Cross-Efficiency Results")

```

Notice that the values on the diagonal, $\theta_{k,k}$ are the same as the original efficiency scores.  

We are going to be doing the cross-efficiency calculations often so let's convert this to a function.  

```{r}
calccrosseff <- function (x,y,vmatrix,umatrix) {

  ND <- nrow(x); NX <- ncol(y); NY <- ncol(y); # Define data size
  
  DMUnames <- rownames(x)
  
  if (is.null(DMUnames)) {DMUnames<-list(c(LETTERS[1:ND]))}  
  
  crosseval  <- matrix(rep(-1.0, ND*ND), nrow=ND, ncol=ND,
                       dimnames=list(DMUnames,DMUnames))
  crosseff   <- matrix(rep(-1.0, ND), nrow=ND, 
                       dimnames=list(DMUnames, "Cross-Eff") )

  for (k in 1:ND) {
     for (j in 1:ND) {crosseval[k,j]<-sum(umatrix[k,]*y[j,])/sum(vmatrix[k,]*x[j,])
                   }
   }
  
   for (j in 1:ND) {crosseff[j]<-sum(crosseval[,j])/ND}

  resultlist <- list(crosseval=crosseval, crosseff=crosseff)
 
  return(resultlist)
}

```

```{r}
rescross <- calccrosseff (x, y, resfunc$vweight, resfunc$uweight)

pander(rbind(rescross$crosseval,t(rescross$crosseff)), 
             caption="Cross-Evaluation Matrix and Cross-Evaluation Scores")

```

## Dealing with Cross-Efficiency's Multiple Optima

As discussed in the chapter on the Multiplier model, DEA often has multiple optima, particularly for efficient DMUs.  While this doesn't affect the efficiency scores, it can have a major impact on the cross-evaluation scores and therefore the cross-efficiency scores.

To examine this issue, we need to invoke a secondary objective function again in the same manner as we did for the envelopment model.

A variety of mechanisms exist for dealing with the issue of multiple optima.

```{r Baker_Cross_Eval_Sec, eval=TRUE}
# Implements secondary objective function to resolve multiple optima

multiplierIOSec <- function (x, y) 
  { 

  ND <- nrow(x); NX <- ncol(y); NY <- ncol(y); # Define data size
  
  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) 

    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] * y[k,r], r = 1:NY), "max") %>%
      add_constraint(sum_expr(vweight[i] * x[k,i], i = 1:NX) == 1) %>%
      add_constraint((sum_expr(uweight[r] * y[j,r], r = 1:NY)-
                    sum_expr(vweight[i] * x[j,i], i = 1:NX)) 
                  <= 0, j = 1:ND)
      result

      result2 <- 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] * sum(y[,r]), r = 1:NY), "max") %>%
          # Modified objective function for 
      add_constraint(sum_expr(vweight[i] * sum(x[,i]), i = 1:NX) == 1) %>%
      add_constraint(result, sum_expr(uweight[r] * sum(y[,r]), r = 1:NY)-
                    sum_expr(vweight[i] * sum(x[,i]), i = 1:NX)  
                  <= 0)
      add_constraint((sum_expr(uweight[r] * y[j,r], r = 1:NY)-
                    sum_expr(vweight[i] * x[j,i], i = 1:NX)) 
                  <= 0, j = 1:ND)

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

    # Get the weights - Output weights
    tempvweight <- get_solution(result2, vweight[i])
    results.vweight[k,] <- tempvweight[,3]

   # Get the weights- Input weights
    tempuweight <- get_solution(result2, uweight[i])
    results.uweight[k,] <- tempuweight[,3]

  } # End of for k loop
  
  resultlist <- list(efficiency=results.efficiency, 
                     vweight=results.vweight, uweight=results.uweight)
 
  return(resultlist)  

} # End of function  

```

Let's explore the impact of the secondary objective function.  The MultiplierDEA package has a cross-efficiency function but the Mal_Ben is a specialized version of the cross-efficiency function that implements the secondary goal of either malevolently trying to minimize other people's scores or benevelolently trying to maximize other people's scores.  Let's compare these to this to the arbitrary cross-efficiency function which does neither.  

```{r}

pander(cbind(x,y),caption="Data Passed to function")

#resfunc <- multiplierIOSec (x=xdata, y=ydata)

resfunc    <- CrossEfficiency (x, y, rts="crs", orientation="input")
resfuncmal <- Mal_Ben (x, y, rts="crs", orientation="input", phase="mal")
resfuncben <- Mal_Ben (x, y, rts="crs", orientation="input", phase="ben")

rownames(resfunc$ce_ave)<-"Single-Phase Cross-Efficiency"
rownames(resfuncmal$ce_ave)<-"Malevolent Cross-Efficiency"
rownames(resfuncben$ce_ave)<-"Benevolent Cross-Efficiency"

pander(rbind(resfunc$ce_ave, resfuncmal$ce_ave, resfuncben$ce_ave) )

```

Next, let's review the implementation of multiple optima in Auro's multiplier model.  To test it, we'll load data from Takamura and Tone's 2003 paper on sites for relocating the capital of Japan.

```{r}

XTT03 <- matrix(c(1,1,1,1,1,1,1),ncol=1,dimnames=list(LETTERS[1:7],"x"))

YTT03 <- matrix(c(5,7,8,4,9,10,4,10,10,7,8,4,2,7,3,3,5,3,4,10,7),
                ncol=3, dimnames=list(LETTERS[1:7], c("C1","C2","C3")))

pander (cbind(XTT03, YTT03), caption = "Input-Output Data for Takamura and Tone, 2003")
```

Comments about cross-efficiency functions from DEAmultiplier package.

* Mal_and_Ben is not defined as a function.
* CrossEfficiency function should allow secondary objective function choice of "Malevelent" or "Benevelent"
* Results from CrossEfficiency do not Takamura and Tone.  
* Naming convention from results of analysis.   Would be good if column names follow convention like my "helper file"-For example lambdas could use L_ as a prefix.
* data validation.  Ex. allow "Ben", "ben" "benevelent", "Benevelent".  Also, RTS could do "CRS" or "crs"

Let's try the Spring case

```{r}
XSpringL <- matrix(c(1,1,1,1,1,1,1,1),ncol=1,dimnames=list(LETTERS[1:8],"x"))

YSpringL <- matrix(c(77,10,88,13,6,35,75,3,75,4,34,33,10,55,61,2,80,60,2,44,27,33,50,90,39,44,81,69,88,77,33,8),ncol=4,dimnames=list(LETTERS[1:8],c("HW","MQ","Proj", "Exam")))

SpringLres <- CrossEfficiency(XSpringL, YSpringL)
```


## Future Issues for Cross-Efficiency Chapter

Issues to consider in the future for this chapter or other chapters:

* Multiple Optima
* Secondary Objective Functions
* Fixed Weighting nature of certain model dimensions