Replicating published results

library(pgt)

The estimators with a numerical example printed in their source paper are validated against it: the weak-G-disposability model of Hampf and Rødseth (2015) and its directional representation against Rødseth (2025) Tables 2 and 3, and the by-production model against the analytic Example 1 of Murty, Russell and Levkoff (2012). These replications run as unit tests on every check; this vignette reproduces them in the open. The remaining estimators carry analytic hand-computed checks in the test suite rather than published-table replications; the final section shows the internal-consistency check for the materials-balance cost model.

Rødseth (2025), Table 2: minimal controlled emissions

The pigfarms data is Table 1 of Rødseth (2025): five pig-finishing farms with feed and piglets as material inputs, meat as the good output, and controlled, uncontrolled and abatement emission columns for nitrogen. The paper prints no nitrogen coefficients for the inputs or the meat. The good-output coefficient \(v = 7/6\) follows from farms A and B, which share inputs: their meat differs by 3 and their uncontrolled emissions by 3.5. Each farm’s feed coefficient then closes its account exactly to the uncontrolled column.

data(pigfarms)
pigfarms
#>   farm feed piglet labor capital meat controlled uncontrolled abatement
#> 1    A   21      2    50     100   10         16         16.7       0.7
#> 2    B   21      2    50     100    7         16         20.2       4.2
#> 3    C   21      2    50     100    7         20         20.2       0.2
#> 4    D   24      3    50     100   11         20         20.4       0.4
#> 5    E   24      3    50     100   10         21         21.6       0.6
v <- 7 / 6
u_feed <- (pigfarms$uncontrolled + v * pigfarms$meat) / pigfarms$feed
tech <- pgt_tech(
  x = pigfarms[, c("feed", "piglet", "labor", "capital")],
  y = pigfarms$meat,
  b = pigfarms$controlled,
  u = cbind(u_feed, 0, 0, 0),
  v = v,
  a = pigfarms$abatement,
  id = pigfarms$farm
)
fit <- pgt(tech, model = "wgd")
data.frame(farm = fit$results$id, b_star = fit$results$b_star,
           z_star = fit$results$z_star, a_star = fit$results$a_star)
#>   farm b_star z_star a_star
#> 1    A     16   16.7    0.7
#> 2    B     16   20.2    4.2
#> 3    C     16   20.2    4.2
#> 4    D     20   20.4    0.4
#> 5    E     16   16.7    0.7

The minimal controlled emissions are \((16, 16, 16, 20, 16)\), which match the paper’s Table 2. Because the abatement output is observed, the fit also reports the projection’s implied uncontrolled emission z_star and abatement a_star, and for every farm these match the paper’s Table 2 uncontrolled and abatement entries. The retained content of the meat decides the peer of farms B and C: at \(v = 0\) farms A and B would tie as peers at \(b = 16\), but the three extra units of meat that farm A produces retain 3.5 units of nitrogen, so farm B is the unique peer. Setting \(v = 0\) (for example by using the uncontrolled column as the only carrier input, \(u = 1\)) leaves the minima unchanged but makes the targets of farms B and C solver-dependent.

Rødseth (2025), Table 3: the directional representation

The directional model (Eq. 13, the indirect form of the directional programme of Eq. 12, with abatement fixed at each farm’s own level, which Proposition 4 equates with the factorially determined multi-output form of Eq. 14) jointly expands the good output and contracts the bad output. It uses the same technology.

fdmo <- pgt(tech, model = "fdmo")
r3 <- fdmo$results
comparison <- data.frame(
  farm = r3$id,
  good_eff = round(r3$good_eff, 3),
  paper_good = c(0, 0, 3, 0, 0.8),
  bad_eff = round(r3$bad_eff, 3),
  paper_bad = c(0, 0, 3.5, 0, 1.0),
  maximal_y = round(r3$maximal_y, 3)
)
comparison
#>   farm good_eff paper_good bad_eff paper_bad maximal_y
#> 1    A    0.000        0.0   0.000       0.0    10.000
#> 2    B    0.000        0.0   0.000       0.0     7.000
#> 3    C    3.000        3.0   3.500       3.5    10.000
#> 4    D    0.000        0.0   0.000       0.0    11.000
#> 5    E    0.789        0.8   0.921       1.0    10.789

Farms A, B and D are on the frontier (gross inefficiency 0), and farm C expands its good output by 3 and cuts its bad output by 3.5, a gross inefficiency of 6.5: farms A to D reproduce Table 3. Farm E is the one divergence, shown side by side above. The paper prints \((\theta_y, \theta_b) = (0.8, 1.0)\); on the printed data the package returns \((0.789, 0.921)\), the exact optimum of the programme. The paper’s data are printed to one decimal, so the divergence is conditional on treating them as exact. Common coefficients \(v = 7/6\), \(u_{\text{feed}} = 31/25\) and \(u_{\text{piglet}} = 349/300\), with the printed inputs, meat and controlled emissions, give uncontrolled emissions and abatement that round to the printed columns, and a fitted Table 3 whose every row, farm E included, rounds to the printed one:

u_c <- c(feed = 31 / 25, piglet = 349 / 300, labor = 0, capital = 0)
x_pig <- as.matrix(pigfarms[, names(u_c)])
z_c <- as.vector(x_pig %*% u_c) - v * pigfarms$meat
a_c <- z_c - pigfarms$controlled
c(uncontrolled = all(round(z_c, 1) == pigfarms$uncontrolled),
  abatement = all(round(a_c, 1) == pigfarms$abatement))
#> uncontrolled    abatement 
#>         TRUE         TRUE
rc <- pgt(pgt_tech(x = x_pig, y = pigfarms$meat, b = pigfarms$controlled,
                   u = u_c, v = v, a = a_c, id = pigfarms$farm),
          model = "fdmo")$results
data.frame(farm = rc$id, gross = round(rc$gross, 3),
           good_eff = round(rc$good_eff, 3),
           bad_eff = round(rc$bad_eff, 3),
           maximal_y = round(rc$maximal_y, 3))
#>   farm gross good_eff bad_eff maximal_y
#> 1    A 0.000    0.000   0.000    10.000
#> 2    B 0.000    0.000   0.000     7.000
#> 3    C 6.500    3.000   3.500    10.000
#> 4    D 0.000    0.000   0.000    11.000
#> 5    E 1.785    0.824   0.961    10.824

These coefficients are one set consistent with the printed data, not necessarily the paper’s own; the farm E row can be confirmed only with the unrounded data.

Murty, Russell and Levkoff (2012): the by-production example

Their Example 1 is a constant-returns technology with one input, one good output \(y\) and one bad output, written \(z\) in the paper (our \(b\)), solved analytically in the paper.

mrl <- pgt_tech(
  x = matrix(c(1, 1, 1, 2, 2), ncol = 1, dimnames = list(NULL, "x")),
  y = c(2, 3 / 2, 2 / 3, 3, 2),
  b = c(4, 1, 2, 5, 3),
  u = 1, polluting = "x",
  id = paste0("DMU", 1:5)
)
bp <- pgt(mrl, model = "byprod", returns = "crs")
rb <- bp$results
data.frame(id = rb$id, output_eff = round(rb$output_eff, 4),
           emission_eff = round(rb$efficiency, 4), mean_eff = round(rb$mean_eff, 4))
#>     id output_eff emission_eff mean_eff
#> 1 DMU1     1.0000       0.2500   0.6250
#> 2 DMU2     0.7500       1.0000   0.8750
#> 3 DMU3     0.3333       0.5000   0.4167
#> 4 DMU4     0.7500       0.4000   0.5750
#> 5 DMU5     0.5000       0.6667   0.5833

For DMU3 the output efficiency is \(1/3\), the emission efficiency \(1/2\) and their mean, the by-production efficiency index, \(5/12\); for DMU2 they are \(3/4\), \(1\) and \(7/8\). These match the analytic scores in the paper.

Coelli, Lauwers and Van Huylenbroeck (2007): coefficients and the EE = TE x EAE decomposition

This section is an internal-consistency check, not a replication of a printed table: the paper’s farm-level data are not published, so the four farms below are invented. What comes from the paper are the phosphorus contents (their footnote 16): feed 0.0124, piglet 0.0117 and meat 0.0117 kg P per kg. The check confirms that the implemented cost model decomposes environmental efficiency as \(EE = TE \times EAE\) with these coefficients; the model’s LPs are verified against hand-computed solutions in the test suite.

farms <- data.frame(
  feed = c(200, 210, 190, 205),
  piglet = c(20, 22, 18, 21),
  meat = c(100, 100, 95, 102)
)
tech_c <- pgt_tech(
  x = as.matrix(farms[, c("feed", "piglet")]),
  y = farms$meat,
  b = 0.0124 * farms$feed + 0.0117 * farms$piglet - 0.0117 * farms$meat,
  u = c(feed = 0.0124, piglet = 0.0117),
  v = 0.0117
)
mbc <- pgt(tech_c, model = "mb_cost", returns = "crs")
rc <- mbc$results
data.frame(id = rc$id, EE = round(rc$efficiency, 4),
           TE = round(rc$te, 4), EAE = round(rc$eae, 4))
#>   id     EE     TE    EAE
#> 1  1 0.9955 1.0000 0.9955
#> 2  2 0.9442 0.9524 0.9914
#> 3  3 1.0000 1.0000 1.0000
#> 4  4 0.9885 0.9951 0.9934
all.equal(rc$efficiency, rc$te * rc$eae)
#> [1] TRUE

The identity \(EE = TE \times EAE\) holds to machine precision. (The identity holds by construction of the reported components, so this check guards the wiring, not the mathematics; the LP solutions themselves are checked by hand in the test suite.)

References