dea.sbm: Slacks-Based Measure (SBM) of Efficiency
In additiveDEA: Additive Data Envelopment Analysis Models
Description
Usage
Arguments
Details
Value
Author(s)
References
See Also
Examples
Description
Calculate additive Data Envelopment Analysis (DEA) efficiency with the Slacks-Based Measure (SBM) of efficiency (Tone, 2001)
Usage
1
2
Arguments
base
A data frame with N rows and S+M columns, where N is the number of Decision-Making Units (DMUs), S is the number of
outputs and M is the number of inputs.
noutput
The number of outputs produced by the DMUs. All DMUs must produce the same number of outputs.
fixed
A numeric vector containing column indices for fixed (non-controllable) outputs and/or inputs (if any) in the data.
Defaults to NULL.
rts
Returns to scale specification. 1 for constant returns to scale and 2 (default) for variable returns to scale.
bound
A data frame with N rows and S+M columns containing user-defined bounds on the slacks of each DMU. If bounds are
supplied by the user in cases where some outputs and/or inputs are fixed, values should be 0 for these fixed
variables. Same for slacks that do not require bounds. Defaults to NULL.
whichDMUs
Numeric vector specifying the line numbers of the DMUs for which efficiency should be calculated. Defaults to NULL,
i.e. by default efficiency is calculated for all DMUs in the dataset.
print.status
Defaults to FALSE. If the solution of the linear program is NA for one or more DMUs, print.status can be set to TRUE to find out why.
Details
The Slacks-Based Measure (SBM) of efficiency (Tone, 2001) is an additive DEA model that maximizes the sum of input and output slacks for each DMU. Unlike other additive DEA models, SBM's objective function has a ratio form, with input slacks summed in the numerator and output slacks summed in the denominator. These sums in the ratio result in an aggregate measure of all inefficiencies that a DMU may exhibit in its inputs and outputs and can thus be seen as a holistic measure of (Pareto-Koopmans) efficiency. SBM weighs input and output slacks in the objective function with the respective inputs used and outputs produced by each DMU. SBM is units invariant, that is, the value of the aggregate inefficiency score of a DMU is independent of the units in which the inputs and outputs are measured, as long as these units are the same for every DMU. For each DMU, SBM returns an efficiency score that is bounded by 0 and 1.
Value
If print.status=FALSE, returns a data frame with N rows and 1+N+S+M columns. Column 1 reports the inefficiency score of each DMU. Columns 2 to N+1 contain the lambda values (intensity variables) indicating the DMU(s) that serve as reference for each DMU. Columns 1+N+1 to 1+N+S report the optimal output slacks for each DMU. Columns 1+N+S+1 to 1+N+S+M report the optimal input slacks for each DMU. When the linear program of a DMU is infeasible, the row of this DMU in the data frame will have NA values.
If print.status=TRUE, returns a list with two elements. The first element is the aforementioned data frame. The second element is a numeric vector indicating the solution status of the linear program (e.g. 0: solution found; 2: problem is infeasible; 5: problem needs scaling; etc. See https://CRAN.R-project.org/package=lpSolveAPI).
Author(s)
Andreas Diomedes Soteriades, andreassot10@yahoo.com
References
Tone K. (2001) A slacks-based measure of efficiency in data envelopment analysis. European Journal of Operational Research, 130, 498–509
See Also
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
# Twelve DMUs, 2 inputs, 2 outputs
# (see Table 1.5 in Cooper et al., 2007):
base <- data.frame(
y1= c(100,150,160,180,94,230,220,152,190,250,260,250),
y2= c(90,50,55,72,66,90,88,80,100,100,147,120),
x1= c(20,19,25,27,22,55,33,31,30,50,53,38),
x2= c(151,131,160,168,158,255,235,206,244,268,306,284))
# Example 1: Get inefficiency scores,
# lambdas and slacks for all DMUs:
dea.sbm(base, noutput= 2, rts= 1)
dea.sbm(base, noutput= 2, rts= 2)
# Example 2: Same as above, but consider output y2 and input x1 as fixed:
dea.sbm(base, noutput= 2, rts= 1, fixed= c(2,3))
dea.sbm(base, noutput= 2, rts= 2, fixed= c(2,3))
# Example 3: Impose an upper bound to all slacks:
# results with bounds
dea.sbm(base, noutput= 2, rts= 1, bound= base/12)$eff
# removing bounds allows for larger inefficiencies:
dea.sbm(base, noutput= 2, rts= 1, bound= NULL)$eff
# check solution status of linear programs when y2 and x1 are fixed
# and y1, x2 slacks are bounded:
bound <- base
fixed <- c(2,3)
bound[fixed] <- 0
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 1, print.status= TRUE)[[2]]
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 2, print.status= TRUE)[[2]]
# Example 4: Get inefficiency scores for DMUs 11 and 12:
bound <- base
fixed <- c(2,3)
bound[fixed] <- 0
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 1, whichDMUs= c(11, 12))$eff
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 2, whichDMUs= c(11, 12))$eff
## The function is currently defined as
function (base, noutput, fixed = NULL, rts = 2, bound = NULL,
whichDMUs = NULL, print.status = FALSE)
{
s <- noutput
m <- ncol(base) - s
n <- nrow(base)
ifelse(!is.null(whichDMUs), nn <- length(whichDMUs), nn <- n)
if (is.null(whichDMUs)) {
whichDMUs <- 1:n
}
re <- data.frame(matrix(0, nrow = nn, ncol = 1 + n + s +
m))
names(re) <- c("eff", paste("lambda", 1:n, sep = ""), paste("slack.y",
1:s, sep = ""), paste("slack.x", 1:m, sep = ""))
slacks <- diag(s + m)
slacks[1:s, ] <- -slacks[1:s, ]
type <- rep("=", s + m)
if (!is.null(fixed)) {
slacks[, fixed] <- 0
type[fixed[fixed <= s]] <- ">="
type[fixed[fixed > s]] <- "<="
}
k <- 0
A.aux <- cbind(t(base), slacks)
index.fixed.y <- fixed[which(fixed %in% 1:s)]
index.fixed.x <- fixed[which(fixed %in% (s + 1):(s + m))]
S <- s - length(index.fixed.y)
M <- m - length(index.fixed.x)
if (print.status == TRUE) {
ifelse(!is.null(whichDMUs), solution.status <- rep(0,
nn), solution.status <- rep(0, n))
}
if (rts == 2 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 2 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
A.bound <- matrix(0, nrow = nrows, ncol = n + s +
m)
kk <- 0
for (j in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1 + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (!is.null(fixed)) {
re <- re[, -(1 + n + fixed)]
}
if (print.status == TRUE) {
reList <- list()
reList[[1]] <- re
reList[[2]] <- solution.status
names(reList[[2]]) <- whichDMUs
re <- reList
}
return(re)
}
Example output
Loading required package: lpSolveAPI
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.0000000 0.0000000 0 0 0 0 0
2 1.0000000 0.0000000 1.0000000 0 0 0 0 0
3 0.8286846 0.3183246 0.8544503 0 0 0 0 0
4 1.0000000 0.0000000 0.0000000 0 1 0 0 0
5 0.7287893 0.2169014 0.9295775 0 0 0 0 0
6 0.6882808 0.8502618 0.9664921 0 0 0 0 0
7 0.8796459 0.6732984 1.0178010 0 0 0 0 0
8 0.7732857 1.1505759 0.2462827 0 0 0 0 0
9 0.9040795 0.8090909 0.7272727 0 0 0 0 0
10 0.7681086 0.7801047 1.1465969 0 0 0 0 0
11 0.8647961 0.9332547 1.2601415 0 0 0 0 0
12 0.9334602 0.8636364 1.0909091 0 0 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.y2 slack.x1
1 0 0 0 0 0 0.00000 0.000000 0.0000000
2 0 0 0 0 0 0.00000 0.000000 0.0000000
3 0 0 0 0 0 0.00000 16.371728 2.3989529
4 0 0 0 0 0 0.00000 0.000000 0.0000000
5 0 0 0 0 0 67.12676 0.000000 0.0000000
6 0 0 0 0 0 0.00000 34.848168 19.6314136
7 0 0 0 0 0 0.00000 23.486911 0.1958115
8 0 0 0 0 0 0.00000 35.865969 3.3091099
9 0 0 0 0 0 0.00000 9.181818 0.0000000
10 0 0 0 0 0 0.00000 27.539267 12.6125654
11 0 0 0 0 0 22.34670 0.000000 10.3922170
12 0 0 0 0 0 0.00000 12.272727 0.0000000
slack.x2
1 0.000000
2 0.000000
3 0.000000
4 0.000000
5 3.473239
6 0.000000
7 0.000000
8 0.000000
9 26.554545
10 0.000000
11 0.000000
12 10.681818
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.00000000 0.0000000 0 0.0000000 0 0 0
2 1.0000000 0.00000000 1.0000000 0 0.0000000 0 0 0
3 0.8500731 0.14196891 0.6870466 0 0.0000000 0 0 0
4 1.0000000 0.00000000 0.0000000 0 1.0000000 0 0 0
5 0.7392282 0.40000000 0.6000000 0 0.0000000 0 0 0
6 0.7409135 0.03584672 0.0000000 0 0.2088999 0 0 0
7 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 1
8 0.7918668 0.38133333 0.4080000 0 0.0000000 0 0 0
9 0.9465246 0.29333333 0.1600000 0 0.0000000 0 0 0
10 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 0
11 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 0
12 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.y2 slack.x1
1 0 0 0 0 0.0000000 0 0.00000 0.000000
2 0 0 0 0 0.0000000 0 0.00000 0.000000
3 0 0 0 0 0.1709845 0 12.64767 2.609326
4 0 0 0 0 0.0000000 0 0.00000 0.000000
5 0 0 0 0 0.0000000 36 0.00000 2.600000
6 0 0 0 0 0.7552534 0 18.89740 19.943140
7 0 0 0 0 0.0000000 0 0.00000 0.000000
8 0 0 0 0 0.2106667 0 0.00000 7.616000
9 0 0 0 0 0.5466667 0 0.00000 0.320000
10 0 0 1 0 0.0000000 0 0.00000 0.000000
11 0 0 0 1 0.0000000 0 0.00000 0.000000
12 0 0 0 0 1.0000000 0 0.00000 0.000000
slack.x2
1 0.00000
2 0.00000
3 0.00000
4 0.00000
5 19.00000
6 0.00000
7 0.00000
8 35.14133
9 23.49333
10 0.00000
11 0.00000
12 0.00000
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.00000000 0.0000000 0 0.0000000 0 0 0
2 1.0000000 0.00000000 1.0000000 0 0.0000000 0 0 0
3 0.8733333 0.00000000 1.2213740 0 0.0000000 0 0 0
4 1.0000000 0.00000000 0.0000000 0 1.0000000 0 0 0
5 0.5705672 0.21690141 0.9295775 0 0.0000000 0 0 0
6 0.7877124 0.00000000 1.9465649 0 0.0000000 0 0 0
7 0.8247267 0.03098592 1.7042254 0 0.0000000 0 0 0
8 0.6503483 0.00000000 1.3488372 0 0.1744186 0 0 0
9 0.8168174 0.56338028 0.9859155 0 0.0000000 0 0 0
10 0.8146766 0.00000000 2.0458015 0 0.0000000 0 0 0
11 0.8993112 0.62711864 0.0000000 0 1.2577684 0 0 0
12 0.8802817 0.53521127 1.4366197 0 0.0000000 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.x2
1 0 0 0 0 0 0.00000 0.000000
2 0 0 0 0 0 0.00000 0.000000
3 0 0 0 0 0 23.20611 0.000000
4 0 0 0 0 0 0.00000 0.000000
5 0 0 0 0 0 67.12676 3.473239
6 0 0 0 0 0 61.98473 0.000000
7 0 0 0 0 0 38.73239 7.067606
8 0 0 0 0 0 81.72093 0.000000
9 0 0 0 0 0 14.22535 29.774648
10 0 0 0 0 0 56.87023 0.000000
11 0 0 0 0 0 29.11017 0.000000
12 0 0 0 0 0 19.01408 14.985915
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.0000000 0.0000000 0 0.000000e+00 0 0 0
2 1.0000000 0.0000000 1.0000000 0 0.000000e+00 0 0 0
3 0.8958333 0.0000000 0.6666667 0 3.333333e-01 0 0 0
4 1.0000000 0.0000000 0.0000000 0 1.000000e+00 0 0 0
5 0.5870344 0.2080537 0.4429530 0 3.489933e-01 0 0 0
6 0.9389356 0.0000000 0.0000000 0 2.857143e-01 0 0 0
7 1.0000000 0.0000000 0.0000000 0 1.718923e-11 0 0 1
8 0.7151124 0.0000000 0.0000000 0 8.933333e-01 0 0 0
9 0.8962792 0.2608696 0.1739130 0 0.000000e+00 0 0 0
10 1.0000000 0.0000000 0.0000000 0 0.000000e+00 0 0 0
11 1.0000000 0.0000000 0.0000000 0 0.000000e+00 0 0 0
12 1.0000000 0.0000000 0.0000000 0 0.000000e+00 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.x2
1 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 0.000000
2 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 0.000000
3 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 16.666667
4 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 0.000000
5 0 0 0.0000000 0.0000000 0.000000e+00 56.067114 9.926174
6 0 0 0.7142857 0.0000000 0.000000e+00 0.000000 15.571429
7 0 0 0.0000000 0.0000000 2.062707e-11 0.000000 0.000000
8 0 0 0.0000000 0.1066667 0.000000e+00 36.533333 23.280000
9 0 0 0.0000000 0.0000000 5.652174e-01 3.478261 21.304348
10 0 0 1.0000000 0.0000000 0.000000e+00 0.000000 0.000000
11 0 0 0.0000000 1.0000000 0.000000e+00 0.000000 0.000000
12 0 0 0.0000000 0.0000000 1.000000e+00 0.000000 0.000000
[1] 1.0000000 1.0000000 0.8663182 1.0000000 0.8461538 0.8992118 0.8896748
[8] 0.8461538 0.9257886 0.8767537 0.9232046 0.9339153
[1] 1.0000000 1.0000000 0.8286846 1.0000000 0.7287893 0.6882808 0.8796459
[8] 0.7732857 0.9040795 0.7681086 0.8647961 0.9334602
1 2 3 4 5 6 7 8 9 10 11 12
0 0 0 0 0 0 0 0 0 0 0 0
1 2 3 4 5 6 7 8 9 10 11 12
0 0 0 0 0 0 0 0 0 0 0 0
[1] 0.8993112 0.8802817
[1] 1 1
function (base, noutput, fixed = NULL, rts = 2, bound = NULL,
whichDMUs = NULL, print.status = FALSE)
{
s <- noutput
m <- ncol(base) - s
n <- nrow(base)
ifelse(!is.null(whichDMUs), nn <- length(whichDMUs), nn <- n)
if (is.null(whichDMUs)) {
whichDMUs <- 1:n
}
re <- data.frame(matrix(0, nrow = nn, ncol = 1 + n + s +
m))
names(re) <- c("eff", paste("lambda", 1:n, sep = ""), paste("slack.y",
1:s, sep = ""), paste("slack.x", 1:m, sep = ""))
slacks <- diag(s + m)
slacks[1:s, ] <- -slacks[1:s, ]
type <- rep("=", s + m)
if (!is.null(fixed)) {
slacks[, fixed] <- 0
type[fixed[fixed <= s]] <- ">="
type[fixed[fixed > s]] <- "<="
}
k <- 0
A.aux <- cbind(t(base), slacks)
index.fixed.y <- fixed[which(fixed %in% 1:s)]
index.fixed.x <- fixed[which(fixed %in% (s + 1):(s + m))]
S <- s - length(index.fixed.y)
M <- m - length(index.fixed.x)
if (print.status == TRUE) {
ifelse(!is.null(whichDMUs), solution.status <- rep(0,
nn), solution.status <- rep(0, n))
}
if (rts == 2 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 2 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
A.bound <- matrix(0, nrow = nrows, ncol = n + s +
m)
kk <- 0
for (j in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1 + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (!is.null(fixed)) {
re <- re[, -(1 + n + fixed)]
}
if (print.status == TRUE) {
reList <- list()
reList[[1]] <- re
reList[[2]] <- solution.status
names(reList[[2]]) <- whichDMUs
re <- reList
}
return(re)
}
additiveDEA documentation built on May 2, 2019, 7:31 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Author(s) References See Also Examples
Description
Calculate additive Data Envelopment Analysis (DEA) efficiency with the Slacks-Based Measure (SBM) of efficiency (Tone, 2001)
Usage
1 2 |
Arguments
base |
A data frame with N rows and S+M columns, where N is the number of Decision-Making Units (DMUs), S is the number of outputs and M is the number of inputs. |
noutput |
The number of outputs produced by the DMUs. All DMUs must produce the same number of outputs. |
fixed |
A numeric vector containing column indices for fixed (non-controllable) outputs and/or inputs (if any) in the data. Defaults to NULL. |
rts |
Returns to scale specification. 1 for constant returns to scale and 2 (default) for variable returns to scale. |
bound |
A data frame with N rows and S+M columns containing user-defined bounds on the slacks of each DMU. If bounds are supplied by the user in cases where some outputs and/or inputs are fixed, values should be 0 for these fixed variables. Same for slacks that do not require bounds. Defaults to NULL. |
whichDMUs |
Numeric vector specifying the line numbers of the DMUs for which efficiency should be calculated. Defaults to NULL, i.e. by default efficiency is calculated for all DMUs in the dataset. |
print.status |
Defaults to FALSE. If the solution of the linear program is NA for one or more DMUs, print.status can be set to TRUE to find out why. |
Details
The Slacks-Based Measure (SBM) of efficiency (Tone, 2001) is an additive DEA model that maximizes the sum of input and output slacks for each DMU. Unlike other additive DEA models, SBM's objective function has a ratio form, with input slacks summed in the numerator and output slacks summed in the denominator. These sums in the ratio result in an aggregate measure of all inefficiencies that a DMU may exhibit in its inputs and outputs and can thus be seen as a holistic measure of (Pareto-Koopmans) efficiency. SBM weighs input and output slacks in the objective function with the respective inputs used and outputs produced by each DMU. SBM is units invariant, that is, the value of the aggregate inefficiency score of a DMU is independent of the units in which the inputs and outputs are measured, as long as these units are the same for every DMU. For each DMU, SBM returns an efficiency score that is bounded by 0 and 1.
Value
If print.status=FALSE, returns a data frame with N rows and 1+N+S+M columns. Column 1 reports the inefficiency score of each DMU. Columns 2 to N+1 contain the lambda values (intensity variables) indicating the DMU(s) that serve as reference for each DMU. Columns 1+N+1 to 1+N+S report the optimal output slacks for each DMU. Columns 1+N+S+1 to 1+N+S+M report the optimal input slacks for each DMU. When the linear program of a DMU is infeasible, the row of this DMU in the data frame will have NA values.
If print.status=TRUE, returns a list with two elements. The first element is the aforementioned data frame. The second element is a numeric vector indicating the solution status of the linear program (e.g. 0: solution found; 2: problem is infeasible; 5: problem needs scaling; etc. See https://CRAN.R-project.org/package=lpSolveAPI).
Author(s)
Andreas Diomedes Soteriades, andreassot10@yahoo.com
References
Tone K. (2001) A slacks-based measure of efficiency in data envelopment analysis. European Journal of Operational Research, 130, 498–509
See Also
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 | # Twelve DMUs, 2 inputs, 2 outputs
# (see Table 1.5 in Cooper et al., 2007):
base <- data.frame(
y1= c(100,150,160,180,94,230,220,152,190,250,260,250),
y2= c(90,50,55,72,66,90,88,80,100,100,147,120),
x1= c(20,19,25,27,22,55,33,31,30,50,53,38),
x2= c(151,131,160,168,158,255,235,206,244,268,306,284))
# Example 1: Get inefficiency scores,
# lambdas and slacks for all DMUs:
dea.sbm(base, noutput= 2, rts= 1)
dea.sbm(base, noutput= 2, rts= 2)
# Example 2: Same as above, but consider output y2 and input x1 as fixed:
dea.sbm(base, noutput= 2, rts= 1, fixed= c(2,3))
dea.sbm(base, noutput= 2, rts= 2, fixed= c(2,3))
# Example 3: Impose an upper bound to all slacks:
# results with bounds
dea.sbm(base, noutput= 2, rts= 1, bound= base/12)$eff
# removing bounds allows for larger inefficiencies:
dea.sbm(base, noutput= 2, rts= 1, bound= NULL)$eff
# check solution status of linear programs when y2 and x1 are fixed
# and y1, x2 slacks are bounded:
bound <- base
fixed <- c(2,3)
bound[fixed] <- 0
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 1, print.status= TRUE)[[2]]
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 2, print.status= TRUE)[[2]]
# Example 4: Get inefficiency scores for DMUs 11 and 12:
bound <- base
fixed <- c(2,3)
bound[fixed] <- 0
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 1, whichDMUs= c(11, 12))$eff
dea.sbm(base, noutput= 2, bound= bound,
fixed= c(2, 3), rts= 2, whichDMUs= c(11, 12))$eff
## The function is currently defined as
function (base, noutput, fixed = NULL, rts = 2, bound = NULL,
whichDMUs = NULL, print.status = FALSE)
{
s <- noutput
m <- ncol(base) - s
n <- nrow(base)
ifelse(!is.null(whichDMUs), nn <- length(whichDMUs), nn <- n)
if (is.null(whichDMUs)) {
whichDMUs <- 1:n
}
re <- data.frame(matrix(0, nrow = nn, ncol = 1 + n + s +
m))
names(re) <- c("eff", paste("lambda", 1:n, sep = ""), paste("slack.y",
1:s, sep = ""), paste("slack.x", 1:m, sep = ""))
slacks <- diag(s + m)
slacks[1:s, ] <- -slacks[1:s, ]
type <- rep("=", s + m)
if (!is.null(fixed)) {
slacks[, fixed] <- 0
type[fixed[fixed <= s]] <- ">="
type[fixed[fixed > s]] <- "<="
}
k <- 0
A.aux <- cbind(t(base), slacks)
index.fixed.y <- fixed[which(fixed %in% 1:s)]
index.fixed.x <- fixed[which(fixed %in% (s + 1):(s + m))]
S <- s - length(index.fixed.y)
M <- m - length(index.fixed.x)
if (print.status == TRUE) {
ifelse(!is.null(whichDMUs), solution.status <- rep(0,
nn), solution.status <- rep(0, n))
}
if (rts == 2 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 2 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
A.bound <- matrix(0, nrow = nrows, ncol = n + s +
m)
kk <- 0
for (j in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1 + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (!is.null(fixed)) {
re <- re[, -(1 + n + fixed)]
}
if (print.status == TRUE) {
reList <- list()
reList[[1]] <- re
reList[[2]] <- solution.status
names(reList[[2]]) <- whichDMUs
re <- reList
}
return(re)
}
|
Example output
Loading required package: lpSolveAPI
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.0000000 0.0000000 0 0 0 0 0
2 1.0000000 0.0000000 1.0000000 0 0 0 0 0
3 0.8286846 0.3183246 0.8544503 0 0 0 0 0
4 1.0000000 0.0000000 0.0000000 0 1 0 0 0
5 0.7287893 0.2169014 0.9295775 0 0 0 0 0
6 0.6882808 0.8502618 0.9664921 0 0 0 0 0
7 0.8796459 0.6732984 1.0178010 0 0 0 0 0
8 0.7732857 1.1505759 0.2462827 0 0 0 0 0
9 0.9040795 0.8090909 0.7272727 0 0 0 0 0
10 0.7681086 0.7801047 1.1465969 0 0 0 0 0
11 0.8647961 0.9332547 1.2601415 0 0 0 0 0
12 0.9334602 0.8636364 1.0909091 0 0 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.y2 slack.x1
1 0 0 0 0 0 0.00000 0.000000 0.0000000
2 0 0 0 0 0 0.00000 0.000000 0.0000000
3 0 0 0 0 0 0.00000 16.371728 2.3989529
4 0 0 0 0 0 0.00000 0.000000 0.0000000
5 0 0 0 0 0 67.12676 0.000000 0.0000000
6 0 0 0 0 0 0.00000 34.848168 19.6314136
7 0 0 0 0 0 0.00000 23.486911 0.1958115
8 0 0 0 0 0 0.00000 35.865969 3.3091099
9 0 0 0 0 0 0.00000 9.181818 0.0000000
10 0 0 0 0 0 0.00000 27.539267 12.6125654
11 0 0 0 0 0 22.34670 0.000000 10.3922170
12 0 0 0 0 0 0.00000 12.272727 0.0000000
slack.x2
1 0.000000
2 0.000000
3 0.000000
4 0.000000
5 3.473239
6 0.000000
7 0.000000
8 0.000000
9 26.554545
10 0.000000
11 0.000000
12 10.681818
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.00000000 0.0000000 0 0.0000000 0 0 0
2 1.0000000 0.00000000 1.0000000 0 0.0000000 0 0 0
3 0.8500731 0.14196891 0.6870466 0 0.0000000 0 0 0
4 1.0000000 0.00000000 0.0000000 0 1.0000000 0 0 0
5 0.7392282 0.40000000 0.6000000 0 0.0000000 0 0 0
6 0.7409135 0.03584672 0.0000000 0 0.2088999 0 0 0
7 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 1
8 0.7918668 0.38133333 0.4080000 0 0.0000000 0 0 0
9 0.9465246 0.29333333 0.1600000 0 0.0000000 0 0 0
10 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 0
11 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 0
12 1.0000000 0.00000000 0.0000000 0 0.0000000 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.y2 slack.x1
1 0 0 0 0 0.0000000 0 0.00000 0.000000
2 0 0 0 0 0.0000000 0 0.00000 0.000000
3 0 0 0 0 0.1709845 0 12.64767 2.609326
4 0 0 0 0 0.0000000 0 0.00000 0.000000
5 0 0 0 0 0.0000000 36 0.00000 2.600000
6 0 0 0 0 0.7552534 0 18.89740 19.943140
7 0 0 0 0 0.0000000 0 0.00000 0.000000
8 0 0 0 0 0.2106667 0 0.00000 7.616000
9 0 0 0 0 0.5466667 0 0.00000 0.320000
10 0 0 1 0 0.0000000 0 0.00000 0.000000
11 0 0 0 1 0.0000000 0 0.00000 0.000000
12 0 0 0 0 1.0000000 0 0.00000 0.000000
slack.x2
1 0.00000
2 0.00000
3 0.00000
4 0.00000
5 19.00000
6 0.00000
7 0.00000
8 35.14133
9 23.49333
10 0.00000
11 0.00000
12 0.00000
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.00000000 0.0000000 0 0.0000000 0 0 0
2 1.0000000 0.00000000 1.0000000 0 0.0000000 0 0 0
3 0.8733333 0.00000000 1.2213740 0 0.0000000 0 0 0
4 1.0000000 0.00000000 0.0000000 0 1.0000000 0 0 0
5 0.5705672 0.21690141 0.9295775 0 0.0000000 0 0 0
6 0.7877124 0.00000000 1.9465649 0 0.0000000 0 0 0
7 0.8247267 0.03098592 1.7042254 0 0.0000000 0 0 0
8 0.6503483 0.00000000 1.3488372 0 0.1744186 0 0 0
9 0.8168174 0.56338028 0.9859155 0 0.0000000 0 0 0
10 0.8146766 0.00000000 2.0458015 0 0.0000000 0 0 0
11 0.8993112 0.62711864 0.0000000 0 1.2577684 0 0 0
12 0.8802817 0.53521127 1.4366197 0 0.0000000 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.x2
1 0 0 0 0 0 0.00000 0.000000
2 0 0 0 0 0 0.00000 0.000000
3 0 0 0 0 0 23.20611 0.000000
4 0 0 0 0 0 0.00000 0.000000
5 0 0 0 0 0 67.12676 3.473239
6 0 0 0 0 0 61.98473 0.000000
7 0 0 0 0 0 38.73239 7.067606
8 0 0 0 0 0 81.72093 0.000000
9 0 0 0 0 0 14.22535 29.774648
10 0 0 0 0 0 56.87023 0.000000
11 0 0 0 0 0 29.11017 0.000000
12 0 0 0 0 0 19.01408 14.985915
eff lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7
1 1.0000000 1.0000000 0.0000000 0 0.000000e+00 0 0 0
2 1.0000000 0.0000000 1.0000000 0 0.000000e+00 0 0 0
3 0.8958333 0.0000000 0.6666667 0 3.333333e-01 0 0 0
4 1.0000000 0.0000000 0.0000000 0 1.000000e+00 0 0 0
5 0.5870344 0.2080537 0.4429530 0 3.489933e-01 0 0 0
6 0.9389356 0.0000000 0.0000000 0 2.857143e-01 0 0 0
7 1.0000000 0.0000000 0.0000000 0 1.718923e-11 0 0 1
8 0.7151124 0.0000000 0.0000000 0 8.933333e-01 0 0 0
9 0.8962792 0.2608696 0.1739130 0 0.000000e+00 0 0 0
10 1.0000000 0.0000000 0.0000000 0 0.000000e+00 0 0 0
11 1.0000000 0.0000000 0.0000000 0 0.000000e+00 0 0 0
12 1.0000000 0.0000000 0.0000000 0 0.000000e+00 0 0 0
lambda8 lambda9 lambda10 lambda11 lambda12 slack.y1 slack.x2
1 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 0.000000
2 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 0.000000
3 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 16.666667
4 0 0 0.0000000 0.0000000 0.000000e+00 0.000000 0.000000
5 0 0 0.0000000 0.0000000 0.000000e+00 56.067114 9.926174
6 0 0 0.7142857 0.0000000 0.000000e+00 0.000000 15.571429
7 0 0 0.0000000 0.0000000 2.062707e-11 0.000000 0.000000
8 0 0 0.0000000 0.1066667 0.000000e+00 36.533333 23.280000
9 0 0 0.0000000 0.0000000 5.652174e-01 3.478261 21.304348
10 0 0 1.0000000 0.0000000 0.000000e+00 0.000000 0.000000
11 0 0 0.0000000 1.0000000 0.000000e+00 0.000000 0.000000
12 0 0 0.0000000 0.0000000 1.000000e+00 0.000000 0.000000
[1] 1.0000000 1.0000000 0.8663182 1.0000000 0.8461538 0.8992118 0.8896748
[8] 0.8461538 0.9257886 0.8767537 0.9232046 0.9339153
[1] 1.0000000 1.0000000 0.8286846 1.0000000 0.7287893 0.6882808 0.8796459
[8] 0.7732857 0.9040795 0.7681086 0.8647961 0.9334602
1 2 3 4 5 6 7 8 9 10 11 12
0 0 0 0 0 0 0 0 0 0 0 0
1 2 3 4 5 6 7 8 9 10 11 12
0 0 0 0 0 0 0 0 0 0 0 0
[1] 0.8993112 0.8802817
[1] 1 1
function (base, noutput, fixed = NULL, rts = 2, bound = NULL,
whichDMUs = NULL, print.status = FALSE)
{
s <- noutput
m <- ncol(base) - s
n <- nrow(base)
ifelse(!is.null(whichDMUs), nn <- length(whichDMUs), nn <- n)
if (is.null(whichDMUs)) {
whichDMUs <- 1:n
}
re <- data.frame(matrix(0, nrow = nn, ncol = 1 + n + s +
m))
names(re) <- c("eff", paste("lambda", 1:n, sep = ""), paste("slack.y",
1:s, sep = ""), paste("slack.x", 1:m, sep = ""))
slacks <- diag(s + m)
slacks[1:s, ] <- -slacks[1:s, ]
type <- rep("=", s + m)
if (!is.null(fixed)) {
slacks[, fixed] <- 0
type[fixed[fixed <= s]] <- ">="
type[fixed[fixed > s]] <- "<="
}
k <- 0
A.aux <- cbind(t(base), slacks)
index.fixed.y <- fixed[which(fixed %in% 1:s)]
index.fixed.x <- fixed[which(fixed %in% (s + 1):(s + m))]
S <- s - length(index.fixed.y)
M <- m - length(index.fixed.x)
if (print.status == TRUE) {
ifelse(!is.null(whichDMUs), solution.status <- rep(0,
nn), solution.status <- rep(0, n))
}
if (rts == 2 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 2 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
A.bound <- matrix(0, nrow = nrows, ncol = n + s +
m)
kk <- 0
for (j in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
add.constraint(lpmodel, xt = c(-1, rep(1, n), rep(0,
s + m)), type = "=", rhs = 0)
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + 1 + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1 + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & is.null(bound)) {
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (rts == 1 & !is.null(bound)) {
index <- which(colSums(bound) != 0)
nrows <- length(index)
A.bound <- matrix(0, nrow = nrows, ncol = n + s + m)
kk <- 0
for (i in index) {
kk <- kk + 1
A.bound[kk, n + index[kk]] <- 1
}
A.bound <- cbind(0, A.bound)
for (i in whichDMUs) {
k <- k + 1
lpmodel <- make.lp(nrow = 0, ncol = 1 + n + s + m)
A <- cbind(-t(base)[, i], A.aux)
A.bound[, 1] <- as.numeric(-bound[i, index])
A.bound[A.bound[, 1] == 0, ] <- 0
A <- rbind(A, A.bound)
xt <- as.numeric((1/S) * (1/base[i, 1:s]))
if (!is.null(fixed)) {
xt[index.fixed.y] <- 0
}
xt <- c(1, rep(0, n), xt, rep(0, m))
add.constraint(lpmodel, xt = xt, type = "=", rhs = 0)
for (j in 1:(s + m)) {
add.constraint(lpmodel, xt = A[j, ], type = type[j],
rhs = 0)
}
for (l in 1:kk) {
add.constraint(lpmodel, xt = A[s + m + l, ],
type = "<=", rhs = 0)
}
set.rhs(lpmodel, b = c(1, rep(0, s + m + l)))
obj <- as.numeric((-1/M) * (1/base[i, (s + 1):(s +
m)]))
if (!is.null(fixed)) {
obj[index.fixed.x - s] <- 0
}
obj <- c(1, rep(0, n + s), obj)
set.objfn(lpmodel, obj = obj)
x <- solve(lpmodel)
if (print.status == TRUE) {
solution.status[k] <- x
}
re[k, ] <- c(get.objective(lpmodel), tail(get.primal.solution(lpmodel),
n + s + m)/get.primal.solution(lpmodel)[1 + 1 +
s + m + sum(colSums(bound) != 0) + 1])
if (x != 0) {
re[k, ] <- rep(NA, ncol(re))
}
}
}
if (!is.null(fixed)) {
re <- re[, -(1 + n + fixed)]
}
if (print.status == TRUE) {
reList <- list()
reList[[1]] <- re
reList[[2]] <- solution.status
names(reList[[2]]) <- whichDMUs
re <- reList
}
return(re)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.