---
title: "Multiple pollutants and heterogeneous coefficients"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Multiple pollutants and heterogeneous coefficients}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6,
                      fig.height = 4)
```

```{r setup}
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 multi-pollutant accounts, DMU-specific
coefficients, several intended outputs and the five-component
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.

```{r}
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
```

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

```{r}
mb <- mb_check(tech2)
table(subset(as.data.frame(mb), violated)$pollutant)
```

## Every pollutant's cap constrains the input-fixed benchmark

The faithful `"wgd"` programme of Eq. 6 uses only the selected
pollutant's data. Cross-account discipline lives in the input-fixed
benchmark: `pgt(model = "wgd_input_fixed")` 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.

```{r}
fit_co2 <- pgt(tech2, model = "wgd_input_fixed", pollutant = "co2")
summary(fit_co2)
```

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:

```{r}
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_input_fixed")
summary(fit1$results$b_star - fit_co2$results$b_star)
```

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.

```{r}
head(cbind(fit_co2$results[c("id", "efficiency", "status")],
           so2_violated = mb$violated[mb$pollutant == "so2"]), 5)
```

## 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.

```{r}
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_input_fixed")
median(fit_het$results$efficiency, na.rm = TRUE)
```

Each plant's cap now uses its own coefficients, so two plants with the
same fuel bill face different materials-balance ceilings. The input
coefficients do not enter the reduced form of the faithful `"wgd"`
programme (only the output coefficients $v_i$ do); producer-specific
$u$ instead surfaces in the `quality` component of the decomposition
below.

## The five-component decomposition

`pgt_decompose(type = "rodseth")` splits the weak-G-disposability
score into the five components of Rødseth (2025, Eq. 11): production
technical efficiency, coefficient quality, production allocative
efficiency, abatement technical efficiency and abatement allocative
efficiency (the stage programmes are stated in
`vignette("models", "pgt")`). A component collapses to exactly 1 when
the data cannot separate it. Here no abatement output is observed and
no pollution-control input is marked, so the two abatement components
are 1, while the heterogeneous coefficients of `tech_het` activate the
quality component.

```{r}
dec <- pgt_decompose(tech_het, type = "rodseth")
summary(dec)
```

The product of the five components equals the `"wgd"` efficiency
exactly, and technology groups are not required for this type.

## Several intended outputs

Good outputs can be a matrix with one named column per output, and the
retained content becomes $v_i'y_i$ with one coefficient per output.
Here the plants sell power and process heat, retaining different
carbon shares per unit, and the programme holds both outputs at the
evaluated plant's levels.

```{r}
y2 <- cbind(power = y, heat = runif(L, 5, 20))
tech_m <- pgt_tech(
  x = cbind(coal = coal, gas = gas),
  y = y2,
  b = b_co2,
  u = u_co2,
  v = c(0.10, 0.05)
)
fit_m <- pgt(tech_m, model = "wgd")
head(fit_m$results, 3)
```

The fit carries one column per output and per-output duals
(`dual_power`, `dual_heat`). Every model except `"fdmo"` handles
several outputs; `mac_curve()` needs a single-output fit, while
`shadow_prices()` reports every output dual.

## References

- Eder, A. (2022). Environmental efficiency measurement when producers
  control pollutants under heterogeneous conditions: a generalization
  of the materials balance approach. *Journal of Productivity
  Analysis*, 57(2), 157-176. doi:10.1007/s11123-021-00623-y
- Rødseth, K. L. (2025). On the development of a unified, nonparametric
  materials balance-based efficiency analysis model and its
  applications. *Journal of Productivity Analysis*, 64(3), 305-319.
  doi:10.1007/s11123-025-00768-0