The Malmquist-Luenberger index (Chung, Färe and Grosskopf 1997)
measures productivity change when part of the output is undesirable, by
scoring each unit on a directional output distance function, the output
counterpart of the directional input distance function of Chambers,
Chung and Färe (1996), that expands the good output and contracts the
bad output together. The adjacent-period version can be infeasible
across periods; the global index of Oh (2010) pools every period into
one reference technology, which removes the infeasibility and makes the
index circular. It does not address the separate inconsistency shown by
Aparicio, Pastor and Zofío (2013), in which the technical-change
component signals regress although the frontier moves toward more good
and less bad output. pgt_ml() computes the global index
under the pollution-generating technology and splits it into efficiency
change and best-practice gap change, \(GML =
EC \times BPC\), with \(GML >
1\) signalling growth.
This estimator is experimental: the global Malmquist-Luenberger index under a materials-balance technology is not yet settled in the literature, and the default technologies do not re-impose the per-DMU materials-balance cap on the projected point.
The technology argument selects the reference
technology. The default "wd" imposes weak disposability of
the bad output in the Kuosmanen (2005) form. With
returns = "crs" its abatement weights are redundant and,
for one bad output, it is the constant-returns technology of Chung, Färe
and Grosskopf (1997, Eq. 3.9), which Oh (2010) also uses; both papers
index several bad outputs jointly, whereas pgt_ml() indexes
one pollutant at a time; the variable-returns default extends it. The
"input_envelope" and "envelope" options treat
the bad output as reducible to the peer emission envelope and are
exploratory companions; the input-keeping variant matches the shape of
the input-fixed benchmark of pgt() rather than its Eq. 6
programme.
The panel below has a frontier that improves between the two periods: the same output is produced with 15 percent less emission.
set.seed(1)
n <- 15
d1 <- data.frame(id = 1:n, y = runif(n, 5, 10), x = runif(n, 8, 12),
b = runif(n, 2, 6))
d2 <- data.frame(id = 1:n, y = d1$y * 1.03, x = d1$x, b = d1$b * 0.85)
d <- rbind(cbind(d1, period = 1), cbind(d2, period = 2))
panel <- pgt_tech(x = d[, "x", drop = FALSE], y = d$y, b = d$b,
period = d$period, id = d$id)
ml <- pgt_ml(panel)
summary(ml)
#> pgt global Malmquist-Luenberger index (technology = wd, returns = vrs)
#>
#> Distribution (GML > 1 is productivity growth):
#> 0% 25% 50% 75% 100%
#> GML 1 1.033 1.0343 1.0459 1.0948
#> EC 1 1.000 1.0000 1.0000 1.0000
#> BPC 1 1.033 1.0343 1.0459 1.0948
#>
#> Geometric mean by period transition:
#> to n gml ec bpc
#> 2 15 1.036 1 1.036The median index exceeds 1, and best-practice gap change carries it: the frontier moved out while each unit kept its position relative to the frontier, so efficiency change stays near 1.
The Philippine rice panel in the frontier package spans
eight years. Treating fertiliser nitrogen as the pollutant potential
gives a nitrogen-based productivity index.
data("riceProdPhil", package = "frontier")
r <- riceProdPhil
uN <- 0.46
rice_panel <- pgt_tech(
x = as.matrix(r[, c("AREA", "LABOR", "NPK", "OTHER")]),
y = r$PROD,
b = uN * r$NPK,
u = c(AREA = 0, LABOR = 0, NPK = uN, OTHER = 0),
v = 0,
period = r$YEARDUM,
id = as.character(r$FMERCODE)
)
rice_ml <- pgt_ml(rice_panel)
summary(rice_ml)
#> pgt global Malmquist-Luenberger index (technology = wd, returns = vrs)
#>
#> Distribution (GML > 1 is productivity growth):
#> 0% 25% 50% 75% 100%
#> GML 0.5786 0.9047 1 1.1360 1.7404
#> EC 0.6317 0.9410 1 1.0743 1.6022
#> BPC 0.5981 0.9205 1 1.1041 1.6719
#>
#> Geometric mean by period transition:
#> to n gml ec bpc
#> 2 43 0.9833 1.0296 0.9551
#> 3 43 1.0483 1.0607 0.9883
#> 4 43 1.0044 0.9807 1.0242
#> 5 43 0.9339 0.9721 0.9607
#> 6 43 1.0129 0.9739 1.0401
#> 7 43 0.8930 1.0120 0.8823
#> 8 43 1.2184 0.9881 1.2331Efficiency scores are point estimates from a finite sample of DMUs.
boot_pgt() gauges their sampling variability by re-solving
each unit against random subsamples of the reference set (m out of n,
without replacement). Subsampling gives valid intervals under minimal
assumptions (Politis and Romano 1994; Politis, Romano and Wolf 1999);
for DEA, Kneip, Simar and Wilson (2008) prove the m-out-of-n bootstrap,
which draws with replacement, consistent under variable returns for the
score of a fixed point, and Simar and Wilson (2011) note that the
resampling can also be done without replacement, although in their
simulations resampling with replacement gave coverage closer to nominal.
boot_pgt() covers cross-sectional efficiency scores only;
no inference is currently available for the pgt_ml() index.
The example takes the sixty plants of one year of the synthetic
steeldemo panel, an independent cross-section.
data(steeldemo)
steel60 <- steeldemo[steeldemo$year == 2021, ]
tech <- pgt_tech(
x = steel60[, c("coal_coke", "other_fuel", "raw_material", "flux",
"capture_energy")],
y = steel60$production, b = steel60$emissions, a = steel60$captured,
v = 0.01467, x_abate = "capture_energy",
group = steel60$route, id = steel60$plant
)
bt <- boot_pgt(tech, model = "wgd", B = 40, seed = 1)
#> Warning: 1 of 60 DMUs have infeasible LPs in more than half of their subsample
#> replicates; their bounds are returned as NA (see per_dmu$n_ok).
head(bt$per_dmu)
#> id group estimate lower upper se n_ok
#> 1 P001 Integrated 0.1258220 0.00000000 0.1258220 0.07672506 40
#> 2 P002 Integrated 0.1536810 0.00000000 0.1534836 0.06553732 40
#> 3 P003 Integrated 0.2801537 0.02650546 0.2801537 0.09771562 40
#> 4 P004 Integrated 0.7870091 0.60591234 0.7870091 0.05612814 40
#> 5 P005 Integrated 0.3612169 0.12117032 0.3612169 0.08785406 40
#> 6 P006 Integrated 0.1719616 0.00000000 0.1715463 0.05198277 40Each unit is evaluated against the subsample frontier alone; the unit
itself is not appended to the subsample. Replicates in which its
programme is infeasible are dropped and counted in
per_dmu$n_ok, and the point estimate can lie above its own
upper endpoint, because the interval targets the true score, which the
estimate overestimates.
The intervals follow the subsampling recipe: the subsample scores are recentred at the point estimate and rescaled by the relative convergence rate \((m/L)^\kappa\), so the interval extends downward from the estimate, matching the direction of the frontier bias (a smaller reference sample can only move the estimated frontier inward). The intervals remain heuristic: resampling theory for data envelopment analysis covers the radial and directional technologies, not the materials-balance programs here, so \(\kappa\) borrows the DEA rate as an approximation. Before relying on the intervals, check that their width is stable in the subsample size (Simar and Wilson 2011).
boot_pgt_sensitivity(tech, model = "wgd", B = 25,
m_grid = c(15, 25, 40), seed = 1)
#> Warning: 2 of 60 DMUs have infeasible LPs in more than half of their subsample
#> replicates; their bounds are returned as NA (see per_dmu$n_ok).
#> Warning: 1 of 60 DMUs have infeasible LPs in more than half of their subsample
#> replicates; their bounds are returned as NA (see per_dmu$n_ok).
#> m median_width median_width_pos mean_se
#> 1 15 0.18617425 0.18617425 0.09873924
#> 2 25 0.15961667 0.15961667 0.08235490
#> 3 40 0.08075829 0.08075829 0.03232624A width that changes little across the grid supports the chosen subsample size; a strong trend means the size, not the data, is driving the interval.
When the technology carries groups, boot_pgt() also
reports intervals for the group mean scores. These are a descriptive
stability band, not supported inference: the replicate statistic
averages all group members against a common subsample frontier, so
replicate variation reflects movement of the frontier only, and the
per-DMU convergence rate is reused for the mean (Kneip, Simar and Wilson
2015 give the theory for aggregates of efficiency scores).
The heuristic label invites the obvious question: how close to
nominal do these intervals get when the truth is known? The package
ships a Monte Carlo study (inst/simulations/coverage.R) for
the envelope model with exponential inefficiency and two known
frontiers, the affine \(b_{\min}(y) = 2 +
0.5y\) and the curved \(b_{\min}(y) = 1
+ 0.05y^2\), where the true score \(\theta = b_{\min}(y)/b\) is the estimand.
Two hundred replicates per cell, B = 100, nominal level 95
per cent. A unit gets no interval when fewer than B/2 of
its subsample programmes are feasible, so the study reports
availability, coverage among units with an interval (“given interval”)
and the share of all units whose interval exists and covers
(“delivered”). The table shows the cells at the defaults \(m = \mathrm{round}(L^{0.7})\) and \(\kappa = 2/3\) from the precomputed results
shipped with the package.
| frontier | L | available | given interval | delivered | frontier (theta > 0.9) | interior | group means |
|---|---|---|---|---|---|---|---|
| affine | 50 | 0.973 | 0.965 | 0.939 | 0.978 | 0.960 | 0.950 |
| affine | 100 | 0.981 | 0.984 | 0.965 | 0.991 | 0.981 | 0.978 |
| curved | 50 | 0.973 | 0.928 | 0.903 | 0.953 | 0.919 | 0.948 |
| curved | 100 | 0.981 | 0.923 | 0.905 | 0.946 | 0.913 | 0.917 |
On the affine frontier the intervals over-cover, and on the curved
one they under-cover by two to three points at both sample sizes. The
missing intervals fall on the units with the largest output. The shipped
results also vary m and kappa at
L = 100: a larger subsample or a faster rate narrows the
intervals and lowers coverage, to 0.77 under kappa = 1 on
the curved frontier. One simulation design is evidence, not proof: the
DEA rate remains borrowed for the materials-balance programmes, which is
why the intervals keep their heuristic label and why
boot_pgt_sensitivity() belongs in the workflow.
Aparicio, J., Pastor, J. T., & Zofío, J. L. (2013). On the inconsistency of the Malmquist-Luenberger index. European Journal of Operational Research, 229(3), 738-742. doi:10.1016/j.ejor.2013.03.031
Chambers, R. G., Chung, Y., & Färe, R. (1996). Benefit and distance functions. Journal of Economic Theory, 70(2), 407-419. doi:10.1006/jeth.1996.0096
Chung, Y. H., Färe, R., & Grosskopf, S. (1997). Productivity and undesirable outputs: A directional distance function approach. Journal of Environmental Management, 51(3), 229-240. doi:10.1006/jema.1997.0146
Kneip, A., Simar, L., & Wilson, P. W. (2008). Asymptotics and consistent bootstraps for DEA estimators in nonparametric frontier models. Econometric Theory, 24(6), 1663-1697. doi:10.1017/S0266466608080651
Kneip, A., Simar, L., & Wilson, P. W. (2015). When bias kills the variance: Central limit theorems for DEA and FDH efficiency scores. Econometric Theory, 31(2), 394-422. doi:10.1017/S0266466614000413
Kuosmanen, T. (2005). Weak disposability in nonparametric production analysis with undesirable outputs. American Journal of Agricultural Economics, 87(4), 1077-1082. doi:10.1111/j.1467-8276.2005.00788.x
Oh, D.-h. (2010). A global Malmquist-Luenberger productivity index. Journal of Productivity Analysis, 34(3), 183-197. doi:10.1007/s11123-010-0178-y
Politis, D. N., & Romano, J. P. (1994). Large sample confidence regions based on subsamples under minimal assumptions. The Annals of Statistics, 22(4), 2031-2050. doi:10.1214/aos/1176325770
Politis, D. N., Romano, J. P., & Wolf, M. (1999). Subsampling. Springer, New York. doi:10.1007/978-1-4612-1554-7
Simar, L., & Wilson, P. W. (2011). Inference by the m out of n bootstrap in nonparametric frontier models. Journal of Productivity Analysis, 36(1), 33-53. doi:10.1007/s11123-010-0200-4