Multiple pollutants and heterogeneous coefficients

library(pgt)

A technology can carry several pollutants, each with its own materials-balance account, and the material flow coefficients can vary across units (heterogeneous input quality; Eder 2022, Rødseth 2025). This vignette demonstrates both features and the staged decomposition, on a small synthetic example.

A two-pollutant technology

Forty plants burn coal and gas. Each fuel carries CO2 potential and SO2 potential, with different coefficients; a share of the CO2 potential is retained in the product (\(v_{CO2} = 0.1\)), none of the SO2. b is a matrix with one named column per pollutant, u a named list with one element per pollutant, and v a named vector.

set.seed(42)
L <- 40
coal <- runif(L, 20, 60)
gas <- runif(L, 10, 40)
y <- runif(L, 15, 45)
u_co2 <- c(coal = 2.5, gas = 1.4)
u_so2 <- c(coal = 0.020, gas = 0.001)
cap_co2 <- 2.5 * coal + 1.4 * gas - 0.1 * y
cap_so2 <- 0.020 * coal + 0.001 * gas
b_co2 <- cap_co2 * runif(L, 0.55, 0.95)
b_so2 <- cap_so2 * runif(L, 0.55, 0.95)
# three plants under-report scrubbing and breach their SO2 account
b_so2[1:3] <- cap_so2[1:3] * runif(3, 1.05, 1.25)

tech2 <- pgt_tech(
  x = cbind(coal = coal, gas = gas),
  y = y,
  b = cbind(co2 = b_co2, so2 = b_so2),
  u = list(co2 = u_co2, so2 = u_so2),
  v = c(co2 = 0.1, so2 = 0),
  group = factor(rep(c("A", "B"), each = L / 2))
)
tech2
#> Pollution-generating technology
#>   DMUs: 40   inputs: 2   good outputs: 1   bad outputs: 2
#>   inputs: coal, gas
#>   material flow coefficients [co2]: u = [2.5, 1.4], v = 0.1
#>   material flow coefficients [so2]: u = [0.020, 0.001], v = 0
#>   groups: A (20), B (20)
#>   materials balance: 3 of 80 DMU-pollutant accounts violate the identity (see mb_check())

mb_check() audits one account per plant and pollutant; the three seeded SO2 violations surface immediately.

mb <- mb_check(tech2)
table(subset(as.data.frame(mb), violated)$pollutant)
#> 
#> so2 
#>   3

Every pollutant’s cap constrains the projection

pgt(model = "wgd") selects one pollutant for the objective through pollutant, but the caps of every pollutant constrain the peer mix (see vignette("models", "pgt") for the collapsed cap rows). A peer that looks attractive on CO2 can be unusable because leaning on it would breach the evaluated plant’s SO2 account.

fit_co2 <- pgt(tech2, model = "wgd", pollutant = "co2")
summary(fit_co2)
#> pgt fit: model = wgd, returns = vrs, peers = all
#>   pollutant: co2
#>   DMUs: 40   failed LPs: 0
#> 
#> Efficiency (b*/b):
#>     0%    25%    50%    75%   100% 
#> 0.2716 0.3897 0.5354 0.8206 1.0000 
#> 
#> By group:
#>  group  n n_na median   mean median_dual_output
#>      A 20    0 0.4644 0.5236              2.645
#>      B 20    0 0.6411 0.6751              2.716

Dropping the SO2 account (a single-pollutant technology on the same data) weakly lowers every b_star, since the projection loses one set of constraints:

tech1 <- pgt_tech(
  x = cbind(coal = coal, gas = gas), y = y, b = b_co2, u = u_co2,
  v = 0.1, group = factor(rep(c("A", "B"), each = L / 2))
)
fit1 <- pgt(tech1, model = "wgd")
summary(fit1$results$b_star - fit_co2$results$b_star)
#>       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
#> -1.421e-14  2.842e-14  4.263e-14  8.148e-13  6.573e-14  1.027e-11

The three SO2 violators show the documented boundary behaviour: a plant violating another pollutant’s identity loses its self-reference, so its CO2 score can exceed 1 (its MB-consistent projection emits more CO2 than the plant reports) or its programme can be infeasible. Neither happens for plants whose accounts all hold.

head(cbind(fit_co2$results[c("id", "efficiency", "status")],
           so2_violated = mb$violated[mb$pollutant == "so2"]), 5)
#>   id efficiency status so2_violated
#> 1  1  0.4613612      0         TRUE
#> 2  2  0.3388767      0         TRUE
#> 3  3  0.8210681      0         TRUE
#> 4  4  0.4116247      0        FALSE
#> 5  5  0.6013043      0        FALSE

DMU-specific coefficients

Coefficients can vary by plant, for example with coal quality (Eder 2022): pass an \(L \times N\) matrix per pollutant instead of a vector. Here the first twenty plants burn higher-carbon coal.

U_co2 <- matrix(u_co2, L, 2, byrow = TRUE,
                dimnames = list(NULL, c("coal", "gas")))
U_co2[1:20, "coal"] <- 2.8
tech_het <- pgt_tech(
  x = cbind(coal = coal, gas = gas), y = y,
  b = pmin(b_co2, (2.8 * coal + 1.4 * gas - 0.1 * y) * 0.95),
  u = U_co2, v = 0.1
)
fit_het <- pgt(tech_het, model = "wgd")
median(fit_het$results$efficiency, na.rm = TRUE)
#> [1] 0.5354456

Each plant’s cap now uses its own coefficients, so two plants with the same fuel bill face different materials-balance ceilings.

The staged decomposition

With a technology group, pgt_decompose(type = "rodseth") decomposes the full weak-G-disposability score in three stages: own-group peers (technical efficiency), all peers (technology gap), and all peers with the input constraints dropped (input-mix component), an exact multiplicative identity for rows with all stages feasible.

dec <- pgt_decompose(tech2, type = "rodseth", pollutant = "co2")
#> Warning: 1 of 40 DMUs have infeasible or failed stage LPs; the affected
#> components are NA. This flags DMUs violating a materials-balance identity; run
#> mb_check().
summary(dec)
#> pgt decomposition: type = rodseth, returns = vrs, DMUs = 40
#>   components: te x technology x ae (product = total)
#> 
#> Group medians (components, total, attribution shares):
#>  group  n n_na    te technology ae total share_te share_technology  share_ae
#>      A 20    1 0.492      0.929  1 0.464    0.932           0.0578 -2.69e-16
#>      B 20    0 0.691      0.933  1 0.591    0.860           0.0859 -3.03e-16
#> 
#> (share_* are medians of DMU-level shares; they need not sum to 1.
#>  n_na counts DMUs with any NA component.)

The caveats of the estimator carry over: a plant violating another pollutant’s account can have te and total above 1, and infeasible stages return NA (see ?pgt_decompose).

References