Data envelopment analysis estimates a production frontier as the smallest set containing the data that satisfies whichever axioms you are willing to assume, and measures each unit’s distance to it. This vignette fits the four estimators in the package to one simulated data set and shows what each of them answers.
The package bundles the Program Follow Through data of Charnes, Cooper and Rhodes (1981): 70 US primary school sites, five inputs describing family background and school resources, three outputs measuring reading, mathematics and self-esteem. This is the application the CCR model was published with.
data(charnes1981)
x <- charnes1981[, paste0("x", 1:5)]
y <- charnes1981[, paste0("y", 1:3)]
dea(x, y, rts = "crs", orientation = "in")
#> --- Data envelopment analysis ---
#> model: radial (Debreu-Farrell)
#> technology: CRS, input orientation
#> DMUs: 70 inputs: 5 outputs: 3 (0.021 sec)
#>
#> theta (input-oriented, 1 = on the frontier)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.7883 0.9002 0.9404 0.9378 1.0000 1.0000
#>
#> efficient DMUs: 19 of 70Nineteen of seventy sites score 1. Before reading anything into that, note the dimensions: five inputs and three outputs is a problem in eight dimensions, and a large efficient set is what eight dimensions produce whatever the data says. The last section of this vignette is about exactly that.
Real data cannot tell you whether an estimator is right, because the
right answer is not known. dea_sim() draws from a
technology whose distance functions have a closed form, so every
estimate below can be scored against a truth that was written down
before the estimator ran.
sim <- dea_sim(n = 150, p = 2, q = 1, returns = 0.9, seed = 1)
sim
#> --- Simulated DEA design ---
#> n = 150 inputs = 2 outputs = 1 elasticity of scale = 0.9
#> input support: [1, 2]
#> inefficiency: exp, E[u] = 0.3
#> true input efficiency theta: 0.1212 / 0.7842 / 0.9994 (min / median / max)
#> attainable MSE slope, vrs: -1 crs: -1.333sim$theta and sim$phi are the true input-
and output-oriented Farrell efficiencies.
fit <- dea(sim$x, sim$y, rts = "vrs", orientation = "in")
fit
#> --- Data envelopment analysis ---
#> model: radial (Debreu-Farrell)
#> technology: VRS, input orientation
#> DMUs: 150 inputs: 2 outputs: 1 (0.009 sec)
#>
#> theta (input-oriented, 1 = on the frontier)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.5664 0.8033 0.8876 0.8689 0.9531 1.0000
#>
#> efficient DMUs: 22 of 150The score is the factor by which every input could be cut at
unchanged output. peers() says which units span the piece
of frontier each DMU is measured against, and slacks() says
what is left over once the radial cut has been made.
head(peers(fit), 4)
#> dmu peer lambda
#> 1 1 86 0.6086542
#> 2 1 116 0.2484980
#> 3 1 106 0.1428478
#> 4 2 116 0.7687440
head(slacks(fit), 4)
#> sx_x1 sx_x2 sy_y1
#> 1 0 0 0.0000000
#> 2 0 0 0.2822239
#> 3 0 0 0.0000000
#> 4 0 0 0.0000000Those two are related. A DMU can be radially efficient — no
common factor cuts every input — and still be dominated,
because one input alone could come down. print() reports
both counts when they differ, and fit$efficient is the
stricter, Pareto–Koopmans judgement.
dea_sbm() folds that distinction into a single number:
the average fraction of each input being wasted, against the average
fraction of each output forgone. It equals 1 exactly on the
Pareto–Koopmans efficient set.
dea_ddf() moves inputs down and outputs up at the same
time, along a direction you choose. It measures inefficiency additively,
so 0 is on the frontier and larger is worse — and, unlike the two above,
it does not divide by the DMU’s own data, so zeros and negative values
are fine.
ddf <- dea_ddf(sim$x, sim$y, direction = "both", rts = "vrs")
summary(ddf$beta)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.00000 0.02913 0.07734 0.09634 0.14886 0.40068The radial models are its two one-sided special cases:
dea_rts(x, y, orientation = "out")
#> --- Returns to scale (output orientation) ---
#> DMUs: 70 inputs: 5 outputs: 3
#>
#> classification
#>
#> irs crs drs
#> 13 19 38
#>
#> scale efficiency
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.8733 0.9819 0.9939 0.9843 1.0000 1.0000
#>
#> at most productive scale size (SE = 1): 19 of 70Scale efficiency is the part of the distance to the constant-returns frontier that comes from operating at the wrong size rather than from being badly run. Note that the classification depends on the orientation, and is meant to: it describes the frontier point a DMU is benchmarked against, and the two orientations project to different points.
Every score above is biased. The estimated frontier is spanned by the observed DMUs, so it lies inside the true one, and every efficiency estimate is pulled toward 1. The bias is not a rounding error:
c(mean_estimate = mean(fit$eff), mean_truth = mean(sim$theta))
#> mean_estimate mean_truth
#> 0.8688975 0.7455112dea_boot() estimates that bias by resampling and returns
a corrected estimate with a confidence interval.
set.seed(1)
b <- dea_boot(fit, B = 100, seed = 1, progress = FALSE)
b
#> --- Simar-Wilson bootstrap ---
#> model: radial, vrs, input orientation
#> B = 100 bandwidth = 0.04764 (silverman) 0.48 sec
#>
#> mean bias: 0.02517 mean se: 0.01386
#> bias correction worthwhile (|bias|/se > 1/sqrt(3)) for 150 of 150 DMUs
#>
#> first rows
#> dmu eff bias_corr lower upper
#> 1 1 0.9143 0.8928 0.8746 0.9111
#> 2 2 0.7401 0.7253 0.7027 0.7376
#> 3 3 0.7950 0.7788 0.7638 0.7919
#> 4 4 0.7155 0.7057 0.6949 0.7123
#> 5 5 0.8798 0.8616 0.8394 0.8762
#> 6 6 0.8704 0.8577 0.8439 0.8693
#> 7 7 0.6847 0.6682 0.6540 0.6798
#> 8 8 0.9558 0.9404 0.9262 0.9527
#> 9 9 0.8138 0.8033 0.7901 0.8117
#> 10 10 0.9553 0.9305 0.8998 0.9522B = 100 keeps this vignette quick; 2000 is the usual
recommendation.
The correct_worthwhile column applies Simar and Wilson’s
own warning: the correction removes a bias and adds the variance of the
estimate of that bias, so it is a net loss where the bias is small
relative to the noise. It is reported rather than applied silently.
Slowly, and more slowly with every variable you add.
dea_rate() gives the slope that log mean squared error can
attain against log sample size:
c(`1 in, 1 out` = dea_rate(1, 1, "vrs"),
`2 in, 2 out` = dea_rate(2, 2, "vrs"),
`4 in, 4 out` = dea_rate(4, 4, "vrs"))
#> 1 in, 1 out 2 in, 2 out 4 in, 4 out
#> -1.3333333 -0.8000000 -0.4444444At four inputs and four outputs the slope is −4/9: a hundredfold
increase in sample size buys about a factor of 8 in mean squared error,
where a parametric estimator would buy 100. This is the strongest
practical argument for keeping the number of variables small, and it is
why the package warns when n < 3(p+q).
Which brings the Program Follow Through analysis back into view.
Seventy units clears the n >= 3(p+q) rule comfortably —
70 against 24 — but the rule is about whether the scores mean anything
at all, not about whether they are precise:
A slope of −0.5 on eight dimensions. The 1981 analysis is still the right analysis of that data; it is simply an analysis whose efficient set is large because the space is large, and the nineteen sites scoring 1 should be read as “not dominated by any observed combination” rather than as “best practice”.
Nineteen sites score 1 and the model cannot separate them. That is not a failure of the data — it is what “not dominated by any observed combination” means in eight dimensions. But it leaves the question a Program Follow Through evaluator actually had, which of these sites to learn from, unanswered.
Cross-efficiency answers it by replacing self-appraisal with peer appraisal. Ordinary DEA lets every site choose the weights that flatter it most; here each site is also scored under every other site’s chosen weights, and the average is its cross-efficiency.
ben <- dea_cross(x, y, secondary = "benevolent")
summary(ben)
#> --- Cross-efficiency ---
#> technology: CRS, input oriented
#> secondary goal: benevolent
#> DMUs: 70 raters: 70 (0.011 sec)
#>
#> cross-efficiency (mean appraisal by every rater)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.7330 0.8121 0.8623 0.8574 0.8880 0.9831
#>
#> own DEA score, for comparison
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.7883 0.9002 0.9404 0.9378 1.0000 1.0000
#>
#> DMUs scoring 1 on their own weights: 19 of 70; ties in the cross ranking: 0
#>
#> ranked by cross-efficiency (first 10 of 70)
#> rank dmu cross own maverick spread
#> 1 1 52 0.9831 1 0.0172 0.2133
#> 2 2 44 0.9712 1 0.0296 0.2663
#> 3 3 62 0.9641 1 0.0372 0.3670
#> 4 4 22 0.9597 1 0.0420 0.3451
#> 5 5 47 0.9594 1 0.0423 0.2453
#> 6 6 17 0.9576 1 0.0443 0.2313
#> 7 7 58 0.9503 1 0.0523 0.1561
#> 8 8 20 0.9358 1 0.0686 0.1435
#> 9 9 27 0.9205 1 0.0863 0.2810
#> 10 10 21 0.9189 1 0.0883 0.2212All seventy are ranked, with no ties. The maverick
column is the gap between what a site awards itself and what its peers
award it: a large value marks a site that looks efficient only under
weights nobody else would choose.
The weights are not unique, and this matters enough to report twice. An efficient site has an entire face of optimal weight vectors, all giving it a score of 1 and all giving other sites different scores. A cross-efficiency computed from whichever vertex the solver stopped at is therefore an artefact. Doyle and Green’s secondary goals bracket it — hold each site’s own score fixed and then either maximize or minimize the average score it awards everyone else:
agg <- dea_cross(x, y, secondary = "aggressive")
range(ben$eff - agg$eff)
#> [1] 0.03250341 0.11695412A ranking that survives from one end of that bracket to the other is in the data. One that does not was in the solver. Reporting a single cross-efficiency without saying which secondary goal produced it is the common practice and it is not a defensible one.
Everything above measures distance to the frontier without asking what anything costs. That is the right default, because prices are usually unavailable. When they are available they answer a question technical efficiency cannot: a site can sit exactly on the frontier and still be spending its budget on the wrong things.
w <- c(1.4, 0.9, 2.1, 1.2, 1.0) # one price list, faced by every site
ce <- dea_cost(x, y, w, rts = "crs")
ce
#> --- Cost efficiency ---
#> technology: CRS
#> DMUs: 70 inputs: 5 outputs: 3 (0.004 sec)
#>
#> cost efficiency (1 = best)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.6229 0.7533 0.7931 0.8087 0.8465 1.0000
#>
#> = technical x allocative
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> technical 0.7883 0.9002 0.9404 0.9378 1.0000 1
#> allocative 0.7226 0.8051 0.8516 0.8616 0.9186 1
#>
#> fully cost efficient: 2 of 70 (19 are technically efficient, of which 17 use the wrong mix for their prices)Cost efficiency factors exactly into the two:
\[CE_o \;=\; \underbrace{\theta_o}_{\text{technical}} \times \underbrace{AE_o}_{\text{allocative}}\]
and the print output above makes the point of the exercise: of the nineteen sites that are technically efficient, only two are also buying the cheapest mix that would keep them there. The other seventeen are on the frontier and still overspending.
sum(ce$technical > 1 - 1e-9) # technically efficient
#> [1] 19
sum(ce$eff == 1) # and allocatively efficient too
#> [1] 2dea_revenue() does the same on the output side.
dea_profit() uses both price vectors, and reports the
Nerlovian profit gap rather than a ratio — observed profit is routinely
zero or negative, so a ratio is undefined exactly where the question is
most interesting. That measure adds into technical plus
allocative rather than multiplying, and it needs variable returns:
maximum profit over a cone is unbounded as soon as one unit is
profitable.
Every radial score can be read off either of two programs. The envelopment form asks which combination of other units dominates this one; the multiplier form asks what prices would make this unit look as good as possible. Strong duality makes them agree on the score, but only the second returns the weights:
fit <- dea(x, y, rts = "crs", orientation = "in", multipliers = TRUE)
round(head(multipliers(fit), 4), 4)
#> v_x1 v_x2 v_x3 v_x4 v_x5 u_y1 u_y2 u_y3
#> 1 0.0000 0.0000 0.0207 0.0000 0.0000 0 0.0000 0.0241
#> 2 0.0008 0.0877 0.0000 0.0000 0.0154 0 0.0035 0.0301
#> 3 0.0016 0.0000 0.0000 0.0255 0.0039 0 0.0000 0.0325
#> 4 0.0000 0.1363 0.0066 0.0000 0.0000 0 0.0000 0.0552Each row satisfies \(v'x_o = 1\) in the caller’s own units, and values no site above break-even. Zeros are common and are the substantive finding: a zero weight is an input the site has chosen to declare irrelevant to its own assessment. Restricting that freedom is what assurance regions are for, and they are not in this version.