From 8bd70c3cb114a1f53e84f9def9772f66ae860946 Mon Sep 17 00:00:00 2001
From: iik1 In efficiency analysis, we often study firms that operate under
-fundamentally different technologies. Steel producers using electric arc
-furnaces (EAF) face a different production possibility set than those
-using the blast furnace-basic oxygen furnace (BF-BOF) route. Hospitals
-in rural areas face different constraints than urban ones. Banks in
-developing economies operate under different regulatory and
-technological environments than those in advanced economies. Following Battese, Rao, and O’Donnell (2004) and O’Donnell, Rao, and
+Battese (2008), we conceive of a single industry
+metatechnology \(T^*\): the set of all input-output
+combinations that are technically feasible in the industry. Each group
+of firms operates within a restricted subset \(T_j \subseteq T^*\) of this metatechnology,
+where the restrictions arise from regulation, the physical environment,
+resource endowments, or the cost of switching production systems. Groups
+do not possess fundamentally different technologies; they face different
+restrictions of a common metatechnology. Standard stochastic frontier analysis (SFA) or data envelopment
-analysis (DEA) applied to the pooled sample implicitly assumes all firms
-share the same technology – an assumption that may be unrealistic.
-Estimating separate frontiers for each group solves this problem but
-makes efficiency scores incomparable across groups: a firm that is 90%
-efficient relative to a less advanced group frontier may actually be
-less productive than a firm that is 70% efficient relative to a more
-advanced frontier. The metafrontier framework, introduced by Battese,
Rao, and O’Donnell (2004) and extended by Huang, Huang, and Liu (2014)
and O’Donnell, Rao, and Battese (2008), resolves this by: This is the default method: The stochastic metafrontier provides a variance-covariance
matrix: The DEA metafrontier computes: For SFA fits, technical efficiencies are computed with the Battese
+and Coelli (1988) conditional expectation estimator by default
+( The For a formatted summary table: The package provides four built-in plot types: Points below the 45-degree line indicate a technology gap (TE* <
TE). The vertical distance from the line reflects the TGR. The poolability test evaluates whether group-specific frontiers are
statistically different from a single pooled frontier: A significant result (small p-value) indicates that the technology
-groups have genuinely different production technologies, justifying the
-metafrontier approach. A significant result (small p-value) indicates that the group
+frontiers genuinely differ, that is, the groups face different
+restrictions of the common metatechnology, justifying the metafrontier
+approach. Every estimation stage of a metafrontier fit can be inspected with
+ The The package supports three distributional assumptions for the
one-sided inefficiency term \(u_i\) in
SFA: Since we used simulated data, we can compare estimated values against
the truth: The package supports panel data via the Battese-Coelli (1992) and
(1995) models. Use the The When group membership is unobserved, use
For additive efficiency decomposition, use DDF-based
metafrontier: Battese, G.E. and Coelli, T.J. (1988). Prediction of firm-level
+technical efficiencies with a generalized frontier production function
+and panel data. Journal of Econometrics, 38(3),
+387–399. Battese, G.E., Rao, D.S.P. and O’Donnell, C.J. (2004). A
metafrontier production function for estimation of technical
efficiencies and technology gaps for firms operating under different
@@ -937,6 +1009,10 @@ Jondrow, J., Lovell, C.A.K., Materov, I.S. and Schmidt, P.
+(1982). On the estimation of technical inefficiency in the stochastic
+frontier production function model. Journal of Econometrics,
+19(2–3), 233–238. O’Donnell, C.J., Rao, D.S.P. and Battese, G.E. (2008).
Metafrontier frameworks for the study of firm-level efficiencies and
technology ratios. Empirical Economics, 34(2),
diff --git a/inst/doc/malmquist.R b/inst/doc/malmquist.R
index 8480aae..e0e40a1 100644
--- a/inst/doc/malmquist.R
+++ b/inst/doc/malmquist.R
@@ -6,9 +6,11 @@ knitr::opts_chunk$set(
fig.height = 5
)
+
## ----setup--------------------------------------------------------------------
library(metafrontier)
+
## ----simulate-panel-----------------------------------------------------------
set.seed(42)
panels <- lapply(1:4, function(t) {
@@ -29,29 +31,35 @@ panel_data <- do.call(rbind, panels)
table(panel_data$group, panel_data$time)
+
## ----malmquist----------------------------------------------------------------
malm <- malmquist_meta(
log_y ~ log_x1 + log_x2,
data = panel_data,
group = "group",
time = "time",
+ id = "id",
orientation = "output",
rts = "crs"
)
malm
+
## ----summary------------------------------------------------------------------
summary(malm)
+
## ----results-table------------------------------------------------------------
head(malm$malmquist, 10)
+
## ----verify-identity----------------------------------------------------------
m <- malm$malmquist
complete <- complete.cases(m[, c("MPI", "TEC", "TGC", "TC")])
all.equal(m$MPI[complete], m$TEC[complete] * m$TGC[complete] * m$TC[complete])
+
## ----group-vs-meta------------------------------------------------------------
# Within-group: MPI_group = EC_group x TC_group
head(malm$group_malmquist)
@@ -59,21 +67,25 @@ head(malm$group_malmquist)
# Metafrontier: MPI_meta = EC_meta x TC_meta
head(malm$meta_malmquist)
+
## ----tgr-dynamics-------------------------------------------------------------
tgr_df <- malm$tgr
# Mean TGR by group and period
aggregate(cbind(TGR_from, TGR_to) ~ group, data = tgr_df, FUN = mean)
+
## ----tgc-by-group-------------------------------------------------------------
aggregate(TGC ~ group, data = tgr_df, FUN = mean)
+
## ----vrs-comparison-----------------------------------------------------------
malm_vrs <- malmquist_meta(
log_y ~ log_x1 + log_x2,
data = panel_data,
group = "group",
time = "time",
+ id = "id",
rts = "vrs"
)
@@ -85,6 +97,7 @@ data.frame(
na.rm = TRUE)
)
+
## ----produc-example, eval = FALSE---------------------------------------------
# library(plm)
# data("Produc", package = "plm")
@@ -94,10 +107,12 @@ data.frame(
# data = Produc,
# group = "region",
# time = "year",
+# id = "state",
# rts = "crs"
# )
# summary(malm_us)
+
## ----utility-example, eval = FALSE--------------------------------------------
# library(sfaR)
# data("utility", package = "sfaR")
@@ -107,6 +122,7 @@ data.frame(
# data = utility,
# group = "regu",
# time = "year",
+# id = "firm",
# rts = "vrs"
# )
# summary(malm_util)
diff --git a/inst/doc/malmquist.Rmd b/inst/doc/malmquist.Rmd
index 4bd6306..f25678a 100644
--- a/inst/doc/malmquist.Rmd
+++ b/inst/doc/malmquist.Rmd
@@ -24,9 +24,10 @@ library(metafrontier)
Standard Malmquist productivity indices measure productivity change over
time by decomposing it into efficiency change and technical change. When
-firms operate under different technologies, however, this decomposition
-misses an important dimension: changes in the *technology gap* between a
-group's frontier and the global best practice.
+groups of firms face different restrictions of a common metatechnology,
+however, this decomposition misses an important dimension: changes in
+the *technology gap* between a group's frontier and the global best
+practice.
The **metafrontier Malmquist TFP index** of O'Donnell, Rao, and Battese
(2008) extends the standard index with a three-way decomposition that
@@ -97,6 +98,7 @@ malm <- malmquist_meta(
data = panel_data,
group = "group",
time = "time",
+ id = "id",
orientation = "output",
rts = "crs"
)
@@ -104,6 +106,10 @@ malm <- malmquist_meta(
malm
```
+The `id` argument identifies firms across periods. Firms not observed
+in both periods of a consecutive pair are dropped from that comparison,
+and `malmquist_meta()` warns about how many observations were dropped.
+
## Detailed results
@@ -190,6 +196,7 @@ malm_vrs <- malmquist_meta(
data = panel_data,
group = "group",
time = "time",
+ id = "id",
rts = "vrs"
)
@@ -218,13 +225,16 @@ malm_us <- malmquist_meta(
data = Produc,
group = "region",
time = "year",
+ id = "state",
rts = "crs"
)
summary(malm_us)
```
Similarly, `sfaR::utility` provides electric utility data with a binary
-`regu` variable (regulated vs. deregulated) over 1986--1996:
+`regu` variable (retail-access vs. traditionally regulated states) over
+1986--1996, where `regu = 1` indicates that the state had enacted
+retail-access (restructuring) legislation:
```{r utility-example, eval = FALSE}
library(sfaR)
@@ -235,14 +245,59 @@ malm_util <- malmquist_meta(
data = utility,
group = "regu",
time = "year",
+ id = "firm",
rts = "vrs"
)
summary(malm_util)
```
+## Caveats
+
+A few limitations should be kept in mind when interpreting the results.
+
+### Interpretation of the index
+
+The Malmquist index, and by extension its metafrontier version, is not
+a proper index in the sense of O'Donnell (2012): it is not
+multiplicatively complete, and its binary comparisons are not
+transitive, so chained comparisons across several periods depend on the
+path taken. Interpretation should therefore be confined to binary
+period-to-period comparisons; avoid cumulating the indices into a
+multi-period productivity level series.
+
+### Cross-period infeasibility
+
+Under `rts = "vrs"`, `"drs"`, `"irs"`, or `"fdh"`, the cross-period
+distance functions (evaluating a period-$t$ observation against the
+period-$s$ frontier) may have no feasible solution. Infeasible cases
+are returned as `NA` (never `Inf`), counted, and reported, so summary
+statistics over the remaining observations should be read with the
+reduced sample in mind. Under `rts = "crs"` the programmes are always
+feasible.
+
+### Firm matching
+
+Observations are matched across consecutive periods via the `id`
+variable. Firms observed in only one period of a pair are dropped from
+that comparison, and `malmquist_meta()` warns about the number of
+dropped observations. With substantial attrition or entry, the index is
+computed on the surviving matched sample only.
+
+### SFA-based index
+
+With `method = "sfa"`, the metafrontier at each period is formed as the
+pointwise maximum of the fitted group frontiers, an approximation that
+the function announces with a message. For exact decompositions, prefer
+the default `method = "dea"`.
+
+
## References
+- O'Donnell, C.J. (2012). An aggregate quantity framework for measuring
+ and decomposing productivity change. *Journal of Productivity
+ Analysis*, 38(3), 255--272.
+
- O'Donnell, C.J., Rao, D.S.P. and Battese, G.E. (2008). Metafrontier
frameworks for the study of firm-level efficiencies and technology
ratios. *Empirical Economics*, 34(2), 231--255.
diff --git a/inst/doc/malmquist.html b/inst/doc/malmquist.html
index 6e12373..ed01111 100644
--- a/inst/doc/malmquist.html
+++ b/inst/doc/malmquist.html
@@ -347,10 +347,10 @@ Standard Malmquist productivity indices measure productivity change
over time by decomposing it into efficiency change and technical change.
-When firms operate under different technologies, however, this
-decomposition misses an important dimension: changes in the
-technology gap between a group’s frontier and the global best
-practice. The metafrontier Malmquist TFP index of O’Donnell,
Rao, and Battese (2008) extends the standard index with a three-way
decomposition that separates within-group dynamics from cross-group
@@ -416,25 +416,31 @@ The Similarly, Introduction to metafrontier
What is a metafrontier?
Estimate the metafrontier
#> ------------------
#> Method: sfa
#> Metafrontier: deterministic
-#> Groups: G1, G2, G3
-#> Total obs: 600
-#> G1: 200 obs
-#> G2: 200 obs
-#> G3: 200 obs
-#>
-#> Group log-likelihoods:
-#> G1: 35.49
-#> G2: 0.84652
-#> G3: -25.399
-#>
-#> Mean TGR by group:
-#> G1: 1
-#> G2: 0.7593
-#> G3: 0.6047Deterministic SFA metafrontier (Battese, Rao, and O’Donnell,
objective = "lp", the default). The
+alternative minimum sum of squared deviations criterion is available via
+objective = "qp"; both criteria are proposed by Battese,
+Rao, and O’Donnell (2004), and the methods vignette
+discusses them in detail.
fit_det <- metafrontier(
@@ -482,49 +501,58 @@ Deterministic SFA metafrontier (Battese, Rao, and O’Donnell,
#> Metafrontier: deterministic
#>
#> --- Group: G1 (n = 200) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 1.025588 0.060881 16.85 <2e-16 ***
-#> log_x1 0.493508 0.009658 51.10 <2e-16 ***
-#> log_x2 0.294569 0.009805 30.04 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: 35.49
-#>
-#> --- Group: G2 (n = 200) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 0.76703 0.06478 11.84 <2e-16 ***
-#> log_x1 0.48428 0.01239 39.10 <2e-16 ***
-#> log_x2 0.29703 0.01210 24.55 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: 0.84652
-#>
-#> --- Group: G3 (n = 200) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 0.53665 0.05513 9.734 <2e-16 ***
-#> log_x1 0.49950 0.01238 40.352 <2e-16 ***
-#> log_x2 0.28205 0.01186 23.782 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: -25.399
-#>
-#> --- Metafrontier ---
-#> Estimate
-#> (Intercept) 1.0256
-#> log_x1 0.4935
-#> log_x2 0.2946
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 1.025588 0.060881 16.846 < 2e-16 ***
+#> log_x1 0.493508 0.009658 51.097 < 2e-16 ***
+#> log_x2 0.294569 0.009805 30.042 < 2e-16 ***
+#> log_sigma_v -1.800325 0.155456 -11.581 < 2e-16 ***
+#> log_sigma_u -1.627386 0.313233 -5.195 2.04e-07 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: 35.49
+#>
+#> --- Group: G2 (n = 200) ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 0.76703 0.06478 11.84 < 2e-16 ***
+#> log_x1 0.48428 0.01239 39.10 < 2e-16 ***
+#> log_x2 0.29703 0.01210 24.55 < 2e-16 ***
+#> log_sigma_v -1.65583 0.16136 -10.26 < 2e-16 ***
+#> log_sigma_u -1.39987 0.28108 -4.98 6.35e-07 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: 0.84652
+#>
+#> --- Group: G3 (n = 200) ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 0.53665 0.05513 9.734 <2e-16 ***
+#> log_x1 0.49950 0.01238 40.352 <2e-16 ***
+#> log_x2 0.28205 0.01186 23.782 <2e-16 ***
+#> log_sigma_v -1.94321 0.14404 -13.490 <2e-16 ***
+#> log_sigma_u -0.89271 0.08650 -10.320 <2e-16 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: -25.399
#>
-#> --- Efficiency Decomposition ---
-#> Group Mean_TE Mean_TGR Mean_TE_star
-#> G1 0.8571 1.0000 0.8571
-#> G2 0.8249 0.7593 0.6263
-#> G3 0.7401 0.6047 0.4475
+#> --- Metafrontier ---
+#> Estimate
+#> (Intercept) 1.0256
+#> log_x1 0.4935
+#> log_x2 0.2946
#>
-#> --- Technology Gap Ratio Summary ---
-#> Group N Mean SD Min Q1 Median Q3 Max
-#> G1 200 1.0000 0.0000 1.0000 1.0000 1.0000 1.0000 1.0000
-#> G2 200 0.7593 0.0099 0.7386 0.7518 0.7592 0.7679 0.7799
-#> G3 200 0.6047 0.0125 0.5769 0.5947 0.6056 0.6143 0.6289
Stochastic metafrontier (Huang, Huang, and Liu, 2014)
@@ -556,52 +584,61 @@ Stochastic metafrontier (Huang, Huang, and Liu, 2014)
#> Metafrontier: stochastic
#>
#> --- Group: G1 (n = 200) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 1.025588 0.060881 16.85 <2e-16 ***
-#> log_x1 0.493508 0.009658 51.10 <2e-16 ***
-#> log_x2 0.294569 0.009805 30.04 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: 35.49
-#>
-#> --- Group: G2 (n = 200) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 0.76703 0.06478 11.84 <2e-16 ***
-#> log_x1 0.48428 0.01239 39.10 <2e-16 ***
-#> log_x2 0.29703 0.01210 24.55 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: 0.84652
-#>
-#> --- Group: G3 (n = 200) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 0.53665 0.05513 9.734 <2e-16 ***
-#> log_x1 0.49950 0.01238 40.352 <2e-16 ***
-#> log_x2 0.28205 0.01186 23.782 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: -25.399
-#>
-#> --- Metafrontier ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 0.784379 0.186836 4.198 2.69e-05 ***
-#> log_x1 0.493270 0.005918 83.349 < 2e-16 ***
-#> log_x2 0.289361 0.005832 49.613 < 2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: 96.292
-#>
-#> --- Efficiency Decomposition ---
-#> Group Mean_TE Mean_TGR Mean_TE_star
-#> G1 0.8571 1.2894 1.1051
-#> G2 0.8249 0.9794 0.8078
-#> G3 0.7401 0.7797 0.5771
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 1.025588 0.060881 16.846 < 2e-16 ***
+#> log_x1 0.493508 0.009658 51.097 < 2e-16 ***
+#> log_x2 0.294569 0.009805 30.042 < 2e-16 ***
+#> log_sigma_v -1.800325 0.155456 -11.581 < 2e-16 ***
+#> log_sigma_u -1.627386 0.313233 -5.195 2.04e-07 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: 35.49
+#>
+#> --- Group: G2 (n = 200) ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 0.76703 0.06478 11.84 < 2e-16 ***
+#> log_x1 0.48428 0.01239 39.10 < 2e-16 ***
+#> log_x2 0.29703 0.01210 24.55 < 2e-16 ***
+#> log_sigma_v -1.65583 0.16136 -10.26 < 2e-16 ***
+#> log_sigma_u -1.39987 0.28108 -4.98 6.35e-07 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: 0.84652
+#>
+#> --- Group: G3 (n = 200) ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 0.53665 0.05513 9.734 <2e-16 ***
+#> log_x1 0.49950 0.01238 40.352 <2e-16 ***
+#> log_x2 0.28205 0.01186 23.782 <2e-16 ***
+#> log_sigma_v -1.94321 0.14404 -13.490 <2e-16 ***
+#> log_sigma_u -0.89271 0.08650 -10.320 <2e-16 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: -25.399
+#>
+#> --- Metafrontier ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 0.784379 0.186836 4.198 2.69e-05 ***
+#> log_x1 0.493270 0.005918 83.349 < 2e-16 ***
+#> log_x2 0.289361 0.005832 49.613 < 2e-16 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: 96.292
#>
-#> --- Technology Gap Ratio Summary ---
-#> Group N Mean SD Min Q1 Median Q3 Max
-#> G1 200 1.2894 0.0099 1.2735 1.2808 1.2887 1.2973 1.3073
-#> G2 200 0.9794 0.0157 0.9428 0.9687 0.9796 0.9906 1.0185
-#> G3 200 0.7797 0.0111 0.7537 0.7717 0.7798 0.7876 0.8027vcov(fit_sto)
@@ -640,7 +677,9 @@ DEA-based metafrontier
#> Mean TGR by group:
#> G1: 0.9959
#> G2: 0.9221
-#> G3: NA
Efficiency scores
# Verify the fundamental identity: TE* = TE x TGR
all.equal(te_star, te * tgr)
#> [1] TRUE
+estimator = "bc88"). The Jondrow et al. (1982) estimator
+is computed and stored alongside it, so you can switch without
+refitting:te_jlms <- efficiencies(fit_det, type = "group", estimator = "jlms")
+cor(te, te_jlms)
+#> [1] 0.9999314Technology gap ratio
technology_gap_ratio() function returns TGR values
grouped by technology:tgr_by_group <- technology_gap_ratio(fit_det)
-lapply(tgr_by_group, summary)
-#> $G1
-#> Min. 1st Qu. Median Mean 3rd Qu. Max.
-#> 1 1 1 1 1 1
-#>
-#> $G2
-#> Min. 1st Qu. Median Mean 3rd Qu. Max.
-#> 0.7386 0.7518 0.7592 0.7593 0.7679 0.7799
-#>
-#> $G3
-#> Min. 1st Qu. Median Mean 3rd Qu. Max.
-#> 0.5769 0.5947 0.6056 0.6047 0.6143 0.6289tgr_by_group <- technology_gap_ratio(fit_det)
+lapply(tgr_by_group, summary)
+#> $G1
+#> Min. 1st Qu. Median Mean 3rd Qu. Max.
+#> 1 1 1 1 1 1
+#>
+#> $G2
+#> Min. 1st Qu. Median Mean 3rd Qu. Max.
+#> 0.7386 0.7518 0.7592 0.7593 0.7679 0.7799
+#>
+#> $G3
+#> Min. 1st Qu. Median Mean 3rd Qu. Max.
+#> 0.5769 0.5947 0.6056 0.6047 0.6143 0.6289tgr_summary(fit_det)
-#> Group N Mean SD Min Q1 Median Q3
-#> 1 G1 200 1.0000000 0.000000000 1.0000000 1.0000000 1.0000000 1.0000000
-#> 2 G2 200 0.7592884 0.009947071 0.7385586 0.7517742 0.7592422 0.7679385
-#> 3 G3 200 0.6047017 0.012456784 0.5769063 0.5947118 0.6056158 0.6142562
-#> Max
-#> 1 1.0000000
-#> 2 0.7798918
-#> 3 0.6289392tgr_summary(fit_det)
+#> Group N Mean SD Min Q1 Median Q3
+#> 1 G1 200 1.0000000 0.000000000 1.0000000 1.0000000 1.0000000 1.0000000
+#> 2 G2 200 0.7592884 0.009947071 0.7385586 0.7517742 0.7592422 0.7679385
+#> 3 G3 200 0.6047017 0.012456784 0.5769063 0.5947118 0.6056158 0.6142562
+#> Max
+#> 1 1.0000000
+#> 2 0.7798918
+#> 3 0.6289392Coefficients
-# Metafrontier coefficients
-coef(fit_det, which = "meta")
-#> (Intercept) log_x1 log_x2
-#> 1.0255883 0.4935078 0.2945688
-
-# Group-specific coefficients
-coef(fit_det, which = "group")
-#> $G1
-#> (Intercept) log_x1 log_x2
-#> 1.0255883 0.4935078 0.2945688
-#>
-#> $G2
-#> (Intercept) log_x1 log_x2
-#> 0.7670252 0.4842818 0.2970328
-#>
-#> $G3
-#> (Intercept) log_x1 log_x2
-#> 0.5366478 0.4994989 0.2820525# Metafrontier coefficients
+coef(fit_det, which = "meta")
+#> (Intercept) log_x1 log_x2
+#> 1.0255883 0.4935078 0.2945688
+
+# Group-specific coefficients
+coef(fit_det, which = "group")
+#> $G1
+#> (Intercept) log_x1 log_x2
+#> 1.0255883 0.4935078 0.2945688
+#>
+#> $G2
+#> (Intercept) log_x1 log_x2
+#> 0.7670252 0.4842818 0.2970328
+#>
+#> $G3
+#> (Intercept) log_x1 log_x2
+#> 0.5366478 0.4994989 0.2820525Model information
-
+
Visualisation
Efficiency scatter
-
-Hypothesis testing
Poolability test
poolability_test(fit_det)
-#>
-#> Likelihood Ratio Test for Poolability of Group Frontiers
-#>
-#> data: metafrontier(formula = log_y ~ log_x1 + log_x2, data = sim$data, group = "group", meta_type = "deterministic")
-#> LR = 504.71, df = 10, p-value < 2.2e-16poolability_test(fit_det)
+#>
+#> Likelihood Ratio Test for Poolability of Group Frontiers
+#>
+#> data: fit_det
+#> LR = 504.71, df = 10, p-value < 2.2e-16Convergence diagnostics
+check_convergence(), which returns one row per stage (each
+group frontier and the metafrontier itself) with the estimation method,
+the optimiser’s convergence code, and a logical convergence
+indicator:check_convergence(fit_det)
+#>
+#> Convergence of estimation stages
+#> --------------------------------
+#> stage method code converged note
+#> group: G1 MLE 0 TRUE
+#> group: G2 MLE 0 TRUE
+#> group: G3 MLE 0 TRUE
+#> metafrontier LP 0 TRUEsummary() method also prints a convergence block, so
+estimation problems are flagged even if you never call
+check_convergence() directly.Inefficiency distributions
# Half-normal (default): u ~ |N(0, sigma_u^2)|
-fit_hn <- metafrontier(log_y ~ log_x1 + log_x2,
- data = sim$data, group = "group",
- dist = "hnormal")
-
-# Truncated normal: u ~ N+(mu, sigma_u^2)
-fit_tn <- metafrontier(log_y ~ log_x1 + log_x2,
- data = sim$data, group = "group",
- dist = "tnormal")
-
-# Exponential: u ~ Exp(1/sigma_u)
-fit_exp <- metafrontier(log_y ~ log_x1 + log_x2,
- data = sim$data, group = "group",
- dist = "exponential")# Half-normal (default): u ~ |N(0, sigma_u^2)|
+fit_hn <- metafrontier(log_y ~ log_x1 + log_x2,
+ data = sim$data, group = "group",
+ dist = "hnormal")
+
+# Truncated normal: u ~ N+(mu, sigma_u^2)
+fit_tn <- metafrontier(log_y ~ log_x1 + log_x2,
+ data = sim$data, group = "group",
+ dist = "tnormal")
+
+# Exponential: u ~ Exp(1/sigma_u)
+fit_exp <- metafrontier(log_y ~ log_x1 + log_x2,
+ data = sim$data, group = "group",
+ dist = "exponential")Comparing true and estimated values
# True vs estimated metafrontier coefficients
-cbind(
- True = sim$params$beta_meta,
- Estimated = coef(fit_det, which = "meta")
-)
-#> True Estimated
-#> (Intercept) 1.0 1.0255883
-#> log_x1 0.5 0.4935078
-#> log_x2 0.3 0.2945688
-
-# True vs estimated mean TGR by group
-true_tgr <- tapply(sim$data$true_tgr, sim$data$group, mean)
-est_tgr <- tapply(fit_det$tgr, fit_det$group_vec, mean)
-cbind(True = true_tgr, Estimated = est_tgr)
-#> True Estimated
-#> G1 1.0000000 1.0000000
-#> G2 0.7788008 0.7592884
-#> G3 0.6065307 0.6047017
-
-# Correlation between true and estimated efficiency
-cor(sim$data$true_te, fit_det$te_group)
-#> [1] 0.80659
-cor(sim$data$true_te_star, fit_det$te_meta)
-#> [1] 0.9400828# True vs estimated metafrontier coefficients
+cbind(
+ True = sim$params$beta_meta,
+ Estimated = coef(fit_det, which = "meta")
+)
+#> True Estimated
+#> (Intercept) 1.0 1.0255883
+#> log_x1 0.5 0.4935078
+#> log_x2 0.3 0.2945688
+
+# True vs estimated mean TGR by group
+true_tgr <- tapply(sim$data$true_tgr, sim$data$group, mean)
+est_tgr <- tapply(fit_det$tgr, fit_det$group_vec, mean)
+cbind(True = true_tgr, Estimated = est_tgr)
+#> True Estimated
+#> G1 1.0000000 1.0000000
+#> G2 0.7788008 0.7592884
+#> G3 0.6065307 0.6047017
+
+# Correlation between true and estimated efficiency
+cor(sim$data$true_te, fit_det$te_group)
+#> [1] 0.8049637
+cor(sim$data$true_te_star, fit_det$te_meta)
+#> [1] 0.9395109Panel SFA Metafrontier
panel argument:# Simulate panel data
-panel_sim <- simulate_panel_metafrontier(
- n_groups = 2, n_firms_per_group = 20, n_periods = 5, seed = 42
-)
-
-# BC92: time-varying inefficiency u_it = u_i * exp(-eta*(t-T))
-fit_panel <- metafrontier(
- log_y ~ log_x1 + log_x2,
- data = panel_sim$data,
- group = "group",
- panel = list(id = "firm", time = "year"),
- panel_dist = "bc92"
-)
-summary(fit_panel)
-
-# The eta parameter captures time-varying inefficiency
-# eta > 0: inefficiency decreasing over time
-# eta < 0: inefficiency increasing over time# Simulate panel data
+panel_sim <- simulate_panel_metafrontier(
+ n_groups = 2, n_firms_per_group = 20, n_periods = 5, seed = 42
+)
+
+# BC92: time-varying inefficiency u_it = u_i * exp(-eta*(t-T))
+fit_panel <- metafrontier(
+ log_y ~ log_x1 + log_x2,
+ data = panel_sim$data,
+ group = "group",
+ panel = list(id = "firm", time = "year"),
+ panel_dist = "bc92"
+)
+summary(fit_panel)
+
+# The eta parameter captures time-varying inefficiency
+# eta > 0: inefficiency decreasing over time
+# eta < 0: inefficiency increasing over timeBootstrap Confidence Intervals for TGR
boot_tgr() function provides parametric and
nonparametric bootstrap confidence intervals for the technology gap
ratio:sim <- simulate_metafrontier(n_groups = 2, n_per_group = 100, seed = 42)
-fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data,
- group = "group", meta_type = "stochastic")
-
-# Nonparametric bootstrap (case resampling within groups)
-boot <- boot_tgr(fit, R = 499, type = "nonparametric", seed = 1)
-print(boot)
-
-# Observation-level CIs
-ci <- confint(boot)
-head(ci)
-
-# Group-level mean TGR CIs
-boot$ci_group
-
-# Parametric bootstrap (resample from estimated error distributions)
-boot_par <- boot_tgr(fit, R = 499, type = "parametric", seed = 1)sim <- simulate_metafrontier(n_groups = 2, n_per_group = 100, seed = 42)
+fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data,
+ group = "group", meta_type = "stochastic")
+
+# Nonparametric bootstrap (case resampling within groups)
+boot <- boot_tgr(fit, R = 499, type = "nonparametric", seed = 1)
+print(boot)
+
+# Observation-level CIs
+ci <- confint(boot)
+head(ci)
+
+# Group-level mean TGR CIs
+boot$ci_group
+
+# Parametric bootstrap (resample from estimated error distributions)
+boot_par <- boot_tgr(fit, R = 499, type = "parametric", seed = 1)Murphy-Topel Variance Correction
@@ -868,66 +936,70 @@ Murphy-Topel Variance Correction
uses fitted values from Stage 1 as regressors. This “generated
regressor” problem means naive standard errors understate uncertainty.
The Murphy-Topel (1985) correction adjusts for this:
-fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data,
- group = "group", meta_type = "stochastic")
-
-# Naive (uncorrected) standard errors
-vcov(fit)
-
-# Murphy-Topel corrected standard errors
-vcov(fit, correction = "murphy-topel")
-
-# Corrected confidence intervals
-confint(fit, correction = "murphy-topel")fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data,
+ group = "group", meta_type = "stochastic")
+
+# Naive (uncorrected) standard errors
+vcov(fit)
+
+# Murphy-Topel corrected standard errors
+vcov(fit, correction = "murphy-topel")
+
+# Corrected confidence intervals
+confint(fit, correction = "murphy-topel")Latent Class Metafrontier
latent_class_metafrontier():sim <- simulate_metafrontier(n_groups = 2, n_per_group = 100, seed = 42)
-
-# Fit with 2 latent classes
-lc <- latent_class_metafrontier(
- log_y ~ log_x1 + log_x2,
- data = sim$data, n_classes = 2,
- n_starts = 5, seed = 123
-)
-print(lc)
-summary(lc)
-
-# Select optimal number of classes via BIC
-bic_table <- select_n_classes(
- log_y ~ log_x1 + log_x2, data = sim$data,
- n_classes_range = 2:4, n_starts = 3, seed = 42
-)
-print(bic_table) # choose n_classes with lowest BICsim <- simulate_metafrontier(n_groups = 2, n_per_group = 100, seed = 42)
+
+# Fit with 2 latent classes
+lc <- latent_class_metafrontier(
+ log_y ~ log_x1 + log_x2,
+ data = sim$data, n_classes = 2,
+ n_starts = 5, seed = 123
+)
+print(lc)
+summary(lc)
+
+# Select optimal number of classes via BIC
+bic_table <- select_n_classes(
+ log_y ~ log_x1 + log_x2, data = sim$data,
+ n_classes_range = 2:4, n_starts = 3, seed = 42
+)
+print(bic_table) # choose n_classes with lowest BICDirectional Distance Functions (DDF)
sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, seed = 42)
-# Use raw (non-log) data for DEA
-sim$data$y <- exp(sim$data$log_y)
-sim$data$x1 <- exp(sim$data$log_x1)
-sim$data$x2 <- exp(sim$data$log_x2)
-
-fit_ddf <- metafrontier(
- y ~ x1 + x2, data = sim$data, group = "group",
- method = "dea", type = "directional", direction = "output"
-)
-summary(fit_ddf)
-
-# Additive decomposition: beta_meta = beta_group + ddf_tgr
-head(data.frame(
- beta_meta = fit_ddf$beta_meta,
- beta_group = fit_ddf$beta_group,
- ddf_tgr = fit_ddf$ddf_tgr
-))sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, seed = 42)
+# Use raw (non-log) data for DEA
+sim$data$y <- exp(sim$data$log_y)
+sim$data$x1 <- exp(sim$data$log_x1)
+sim$data$x2 <- exp(sim$data$log_x2)
+
+fit_ddf <- metafrontier(
+ y ~ x1 + x2, data = sim$data, group = "group",
+ method = "dea", type = "directional", direction = "output"
+)
+summary(fit_ddf)
+
+# Additive decomposition: beta_meta = beta_group + ddf_tgr
+head(data.frame(
+ beta_meta = fit_ddf$beta_meta,
+ beta_group = fit_ddf$beta_group,
+ ddf_tgr = fit_ddf$ddf_tgr
+))References
+
References
to estimating the metafrontier production function based on a stochastic
frontier framework. Journal of Productivity Analysis, 42(3),
241–254.Metafrontier Malmquist Productivity
Motivation
Computing the index
data = panel_data,
group = "group",
time = "time",
- orientation = "output",
- rts = "crs"
-)
-
-malm
-#>
-#> Metafrontier Malmquist TFP Index
-#> ================================
-#> Orientation: output
-#> RTS: crs
-#> Groups: G1, G2
-#> Periods: 1 -> 2 -> 3 -> 4
-#> Observations: 300
-#>
-#> Mean decomposition (M* = TEC x TGC x TC*):
-#> MPI = 1.085
-#> TEC = 1.123
-#> TGC = 1.127
-#> TC* = 0.9818id argument identifies firms across periods. Firms
+not observed in both periods of a consecutive pair are dropped from that
+comparison, and malmquist_meta() warns about how many
+observations were dropped.Detailed results
@@ -447,56 +453,58 @@ Detailed results
#>
#> Call:
#> malmquist_meta(formula = log_y ~ log_x1 + log_x2, data = panel_data,
-#> group = "group", time = "time", orientation = "output", rts = "crs")
-#>
-#> Orientation: output
-#> RTS: crs
-#> Groups: G1, G2
-#> Periods: 1 -> 2 -> 3 -> 4
-#> Observations: 300
-#>
-#> Overall means:
-#> MPI = 1.085
-#> TEC = 1.123
-#> TGC = 1.127
-#> TC* = 0.9818
-#>
-#> --- Three-Way Decomposition by Group ---
-#> M* = TEC x TGC x TC*
-#>
-#> Group: G1 (n = 150 )
-#> MPI TEC TGC TC
-#> 1.0726 1.2388 0.9833 0.9874
-#>
-#> Group: G2 (n = 150 )
-#> MPI TEC TGC TC
-#> 1.0962 1.0082 1.2703 0.9766
-#>
-#> --- By Period ---
+#> group = "group", time = "time", id = "id", orientation = "output",
+#> rts = "crs")
+#>
+#> Method: dea
+#> Orientation: output
+#> RTS: crs
+#> Groups: G1, G2
+#> Periods: 1 -> 2 -> 3 -> 4
+#> Observations: 300
+#>
+#> Overall means:
+#> MPI = 1.145
+#> TEC = 1.123
+#> TGC = 1.127
+#> TC* = 1.002
+#>
+#> --- Three-Way Decomposition by Group ---
+#> M* = TEC x TGC x TC*
+#>
+#> Group: G1 (n = 150 )
+#> MPI TEC TGC TC
+#> 1.1806 1.2388 0.9833 1.0273
+#>
+#> Group: G2 (n = 150 )
+#> MPI TEC TGC TC
+#> 1.1096 1.0082 1.2703 0.9767
#>
-#> Period 1 -> 2
-#> MPI TEC TGC TC
-#> 1.1697 1.0360 1.3314 0.9044
-#>
-#> Period 2 -> 3
-#> MPI TEC TGC TC
-#> 1.0470 1.5288 1.0757 0.6282
-#>
-#> Period 3 -> 4
-#> MPI TEC TGC TC
-#> 1.0364 0.8056 0.9732 1.4224
-#>
-#> --- Technology Gap Ratios ---
+#> --- By Period ---
+#>
+#> Period 1 -> 2
+#> MPI TEC TGC TC
+#> 1.3082 1.0360 1.3314 0.9567
+#>
+#> Period 2 -> 3
+#> MPI TEC TGC TC
+#> 1.0250 1.5288 1.0757 0.6300
+#>
+#> Period 3 -> 4
+#> MPI TEC TGC TC
+#> 1.1021 0.8056 0.9732 1.4194
#>
-#> Group: G1
-#> Mean TGR (from): 0.987
-#> Mean TGR (to): 0.9687
-#> Mean TGC: 0.9833
-#>
-#> Group: G2
-#> Mean TGR (from): 0.7011
-#> Mean TGR (to): 0.8234
-#> Mean TGC: 1.27Interpreting the decomposition
@@ -559,7 +567,7 @@ Within-group vs metafrontier Malmquist
# Within-group: MPI_group = EC_group x TC_group
head(malm$group_malmquist)
#> id group period_from period_to MPI_group EC_group TC_group
-#> 1 1 G1 1 2 NA 2.7812965 NA
+#> 1 1 G1 1 2 6.164525 2.7812965 2.2164215
#> 2 2 G1 1 2 1.757668 2.2260475 0.7895912
#> 3 3 G1 1 2 0.556865 0.6482012 0.8590929
#> 4 4 G1 1 2 1.018773 1.1385028 0.8948361
@@ -607,21 +615,24 @@ Returns to scale assumptions
data = panel_data,
group = "group",
time = "time",
- rts = "vrs"
-)
-
-# Compare mean MPI under CRS vs VRS
-data.frame(
- CRS = colMeans(malm$malmquist[, c("MPI", "TEC", "TGC", "TC")],
- na.rm = TRUE),
- VRS = colMeans(malm_vrs$malmquist[, c("MPI", "TEC", "TGC", "TC")],
- na.rm = TRUE)
-)
-#> CRS VRS
-#> MPI 1.0847509 1.018219
-#> TEC 1.1234751 1.004931
-#> TGC 1.1267757 1.014285
-#> TC 0.9818377 1.005744Using real-world panel data
@@ -637,12 +648,15 @@ Using real-world panel data
data = Produc,
group = "region",
time = "year",
- rts = "crs"
-)
-summary(malm_us)sfaR::utility provides electric utility data
-with a binary regu variable (regulated vs. deregulated)
-over 1986–1996:regu variable (retail-access
+vs. traditionally regulated states) over 1986–1996, where
+regu = 1 indicates that the state had enacted retail-access
+(restructuring) legislation:
library(sfaR)
data("utility", package = "sfaR")
@@ -651,13 +665,60 @@ Using real-world panel data
data = utility,
group = "regu",
time = "year",
- rts = "vrs"
-)
-summary(malm_util)A few limitations should be kept in mind when interpreting the +results.
+The Malmquist index, and by extension its metafrontier version, is +not a proper index in the sense of O’Donnell (2012): it is not +multiplicatively complete, and its binary comparisons are not +transitive, so chained comparisons across several periods depend on the +path taken. Interpretation should therefore be confined to binary +period-to-period comparisons; avoid cumulating the indices into a +multi-period productivity level series.
+Under rts = "vrs", "drs",
+"irs", or "fdh", the cross-period distance
+functions (evaluating a period-\(t\)
+observation against the period-\(s\)
+frontier) may have no feasible solution. Infeasible cases are returned
+as NA (never Inf), counted, and reported, so
+summary statistics over the remaining observations should be read with
+the reduced sample in mind. Under rts = "crs" the
+programmes are always feasible.
Observations are matched across consecutive periods via the
+id variable. Firms observed in only one period of a pair
+are dropped from that comparison, and malmquist_meta()
+warns about the number of dropped observations. With substantial
+attrition or entry, the index is computed on the surviving matched
+sample only.
With method = "sfa", the metafrontier at each period is
+formed as the pointwise maximum of the fitted group frontiers, an
+approximation that the function announces with a message. For exact
+decompositions, prefer the default method = "dea".
O’Donnell, C.J. (2012). An aggregate quantity framework for +measuring and decomposing productivity change. Journal of +Productivity Analysis, 38(3), 255–272.
O’Donnell, C.J., Rao, D.S.P. and Battese, G.E. (2008). Metafrontier frameworks for the study of firm-level efficiencies and technology ratios. Empirical Economics, 34(2), diff --git a/inst/doc/methods.R b/inst/doc/methods.R index c610eca..8aa3abd 100644 --- a/inst/doc/methods.R +++ b/inst/doc/methods.R @@ -6,9 +6,11 @@ knitr::opts_chunk$set( fig.height = 5 ) + ## ----setup-------------------------------------------------------------------- library(metafrontier) + ## ----det-example-------------------------------------------------------------- sim <- simulate_metafrontier( n_groups = 2, n_per_group = 300, @@ -30,11 +32,29 @@ coef(fit_det, which = "meta") # Group coefficients for comparison coef(fit_det, which = "group") + ## ----verify-envelop----------------------------------------------------------- meta_b0 <- coef(fit_det, which = "meta")[1] group_b0 <- sapply(coef(fit_det, which = "group"), `[`, 1) meta_b0 >= group_b0 + +## ----det-qp------------------------------------------------------------------- +fit_qp <- metafrontier( + log_y ~ log_x1 + log_x2, + data = sim$data, + group = "group", + meta_type = "deterministic", + objective = "qp" +) + +# LP and QP criteria typically give very similar coefficients +cbind( + LP = coef(fit_det, which = "meta"), + QP = coef(fit_qp, which = "meta") +) + + ## ----sto-example-------------------------------------------------------------- fit_sto <- metafrontier( log_y ~ log_x1 + log_x2, @@ -45,6 +65,7 @@ fit_sto <- metafrontier( summary(fit_sto) + ## ----sto-inference------------------------------------------------------------ # Variance-covariance matrix vcov(fit_sto) @@ -52,10 +73,12 @@ vcov(fit_sto) # Log-likelihood of the metafrontier model logLik(fit_sto) + ## ----tgr-range---------------------------------------------------------------- tgr_vals <- efficiencies(fit_sto, type = "tgr") summary(tgr_vals) + ## ----dea-example-------------------------------------------------------------- # CRS metafrontier fit_crs <- metafrontier( @@ -81,6 +104,7 @@ cbind( VRS = tapply(fit_vrs$tgr, fit_vrs$group_vec, mean) ) + ## ----compare-methods---------------------------------------------------------- # Compare TGR estimates across methods tgr_det <- tapply(fit_det$tgr, fit_det$group_vec, mean) @@ -96,9 +120,11 @@ comparison <- data.frame( ) round(comparison, 4) + ## ----poolability-------------------------------------------------------------- poolability_test(fit_det) + ## ----monte-carlo, eval=FALSE-------------------------------------------------- # # Monte Carlo: check parameter recovery over 100 replications # set.seed(1) diff --git a/inst/doc/methods.Rmd b/inst/doc/methods.Rmd index 6d716ba..2f45be4 100644 --- a/inst/doc/methods.Rmd +++ b/inst/doc/methods.Rmd @@ -43,7 +43,12 @@ The group-specific technical efficiency is: $$TE_{ij} = \exp(-u_{ij}) \in (0, 1]$$ -estimated via the Jondrow et al. (1982) conditional mean estimator. +estimated by default via the Battese and Coelli (1988) conditional +expectation estimator $E[\exp(-u_{ij}) \mid \varepsilon_{ij}]$ +(`estimator = "bc88"`). The Jondrow et al. (1982) estimator +$\exp(-E[u_{ij} \mid \varepsilon_{ij}])$ is also computed and stored, and +`efficiencies(fit, estimator = "jlms")` switches between the two without +refitting. ### 1.2 The metafrontier @@ -53,9 +58,14 @@ that: $$x'\beta^* \ge x'\beta_j \quad \text{for all } x \text{ and all } j$$ -That is, the metafrontier weakly dominates all group frontiers. It -represents the production technology available to firms with unrestricted -access to all technologies. +That is, the metafrontier weakly dominates all group frontiers. It is +the boundary of the industry **metatechnology** $T^*$: the set of all +input-output combinations that are technically feasible in the industry +(Battese, Rao and O'Donnell, 2004; O'Donnell, Rao and Battese, 2008). +Each group frontier bounds a restricted subset $T_j \subseteq T^*$ of +this metatechnology, where the restrictions reflect regulation, the +physical environment, resource endowments, or the cost of switching +production systems. ### 1.3 The efficiency decomposition @@ -76,15 +86,36 @@ that input mix; values below 1 indicate a technology gap. ### 2.1 Estimation -After obtaining group estimates $\hat\beta_j$ in Stage 1, the metafrontier -parameters $\hat\beta^*$ are estimated by solving: +After obtaining group estimates $\hat\beta_j$ in Stage 1, the +metafrontier parameters $\hat\beta^*$ are chosen so that the metafrontier +lies as close as possible to the group frontiers while enveloping all of +them. Battese, Rao and O'Donnell (2004) propose two optimisation +criteria, both implemented in the package via the `objective` argument. -$$\min_{\beta^*} \sum_{j=1}^{J} \sum_{i=1}^{n_j} \left(x_{ij}'\beta^* - x_{ij}'\hat\beta_j\right)^2$$ +**Minimum sum of absolute deviations** (`objective = "lp"`, the default): + +$$\min_{\beta^*} \sum_{j=1}^{J} \sum_{i=1}^{n_j} \left|x_{ij}'\beta^* - x_{ij}'\hat\beta_j\right|$$ $$\text{subject to: } x_{ij}'\beta^* \ge x_{ij}'\hat\beta_j \quad \forall\, i, j$$ -This is a convex quadratic program. The `metafrontier` package solves it -using `constrOptim()` from base R, which implements an adaptive barrier -algorithm for linearly constrained optimisation. +Because the envelopment constraints force every deviation to be +non-negative, the absolute values can be dropped and the problem reduces +to a linear programme (O'Donnell, Rao and Battese, 2008, Eqs. 23-25). +This LP is the primary presentation in O'Donnell, Rao and Battese (2008) +and is the package default; it is solved with `lpSolveAPI`. + +**Minimum sum of squared deviations** (`objective = "qp"`): + +$$\min_{\beta^*} \sum_{j=1}^{J} \sum_{i=1}^{n_j} \left(x_{ij}'\beta^* - x_{ij}'\hat\beta_j\right)^2$$ + +subject to the same envelopment constraints. This is a convex quadratic +programme, solved exactly with `quadprog` when that package is +available, with an adaptive-barrier fallback via `constrOptim()` from +base R otherwise. + +The two criteria weight deviations differently (the squared criterion +penalises large gaps between the metafrontier and a group frontier more +heavily), but in practice they typically produce very similar +coefficients. ### 2.2 Properties @@ -129,6 +160,25 @@ group_b0 <- sapply(coef(fit_det, which = "group"), `[`, 1) meta_b0 >= group_b0 ``` +To use the minimum sum of squared deviations criterion instead, set +`objective = "qp"`: + +```{r det-qp} +fit_qp <- metafrontier( + log_y ~ log_x1 + log_x2, + data = sim$data, + group = "group", + meta_type = "deterministic", + objective = "qp" +) + +# LP and QP criteria typically give very similar coefficients +cbind( + LP = coef(fit_det, which = "meta"), + QP = coef(fit_qp, which = "meta") +) +``` + ## 3. Stochastic metafrontier (Huang, Huang, and Liu, 2014) @@ -349,8 +399,9 @@ on the pooled sample and $LL_j$ are the group-specific log-likelihoods. poolability_test(fit_det) ``` -A significant test (p < 0.05) confirms that the groups operate under -different technologies and the metafrontier decomposition is warranted. +A significant test (p < 0.05) confirms that the group frontiers differ, +that is, the groups face different restrictions of the common +metatechnology, and the metafrontier decomposition is warranted. ## 8. Simulation for Monte Carlo studies @@ -397,6 +448,10 @@ The `simulate_metafrontier()` function supports: ## References +- Battese, G.E. and Coelli, T.J. (1988). Prediction of firm-level + technical efficiencies with a generalized frontier production function + and panel data. *Journal of Econometrics*, 38(3), 387--399. + - Battese, G.E., Rao, D.S.P. and O'Donnell, C.J. (2004). A metafrontier production function for estimation of technical efficiencies and technology gaps for firms operating under different technologies. diff --git a/inst/doc/methods.html b/inst/doc/methods.html index 35e4b8d..abb6654 100644 --- a/inst/doc/methods.html +++ b/inst/doc/methods.html @@ -363,8 +363,13 @@
The group-specific technical efficiency is:
\[TE_{ij} = \exp(-u_{ij}) \in (0, 1]\]
-estimated via the Jondrow et al. (1982) conditional mean -estimator.
+estimated by default via the Battese and Coelli (1988) conditional
+expectation estimator \(E[\exp(-u_{ij}) \mid
+\varepsilon_{ij}]\) (estimator = "bc88"). The
+Jondrow et al. (1982) estimator \(\exp(-E[u_{ij} \mid \varepsilon_{ij}])\) is
+also computed and stored, and
+efficiencies(fit, estimator = "jlms") switches between the
+two without refitting.
\[x'\beta^* \ge x'\beta_j \quad \text{for all } x \text{ and all } j\]
-That is, the metafrontier weakly dominates all group frontiers. It -represents the production technology available to firms with -unrestricted access to all technologies.
+That is, the metafrontier weakly dominates all group frontiers. It is +the boundary of the industry metatechnology \(T^*\): the set of all input-output +combinations that are technically feasible in the industry (Battese, Rao +and O’Donnell, 2004; O’Donnell, Rao and Battese, 2008). Each group +frontier bounds a restricted subset \(T_j +\subseteq T^*\) of this metatechnology, where the restrictions +reflect regulation, the physical environment, resource endowments, or +the cost of switching production systems.
After obtaining group estimates \(\hat\beta_j\) in Stage 1, the metafrontier -parameters \(\hat\beta^*\) are -estimated by solving:
+parameters \(\hat\beta^*\) are chosen +so that the metafrontier lies as close as possible to the group +frontiers while enveloping all of them. Battese, Rao and O’Donnell +(2004) propose two optimisation criteria, both implemented in the +package via theobjective argument.
+Minimum sum of absolute deviations
+(objective = "lp", the default):
\[\min_{\beta^*} \sum_{j=1}^{J} -\sum_{i=1}^{n_j} \left(x_{ij}'\beta^* - -x_{ij}'\hat\beta_j\right)^2\] \[\text{subject to: } x_{ij}'\beta^* \ge +\sum_{i=1}^{n_j} \left|x_{ij}'\beta^* - +x_{ij}'\hat\beta_j\right|\] \[\text{subject to: } x_{ij}'\beta^* \ge x_{ij}'\hat\beta_j \quad \forall\, i, j\]
-This is a convex quadratic program. The metafrontier
-package solves it using constrOptim() from base R, which
-implements an adaptive barrier algorithm for linearly constrained
-optimisation.
Because the envelopment constraints force every deviation to be
+non-negative, the absolute values can be dropped and the problem reduces
+to a linear programme (O’Donnell, Rao and Battese, 2008, Eqs. 23-25).
+This LP is the primary presentation in O’Donnell, Rao and Battese (2008)
+and is the package default; it is solved with
+lpSolveAPI.
Minimum sum of squared deviations
+(objective = "qp"):
\[\min_{\beta^*} \sum_{j=1}^{J} +\sum_{i=1}^{n_j} \left(x_{ij}'\beta^* - +x_{ij}'\hat\beta_j\right)^2\]
+subject to the same envelopment constraints. This is a convex
+quadratic programme, solved exactly with quadprog when that
+package is available, with an adaptive-barrier fallback via
+constrOptim() from base R otherwise.
The two criteria weight deviations differently (the squared criterion +penalises large gaps between the metafrontier and a group frontier more +heavily), but in practice they typically produce very similar +coefficients.
To use the minimum sum of squared deviations criterion instead, set
+objective = "qp":
fit_qp <- metafrontier(
+ log_y ~ log_x1 + log_x2,
+ data = sim$data,
+ group = "group",
+ meta_type = "deterministic",
+ objective = "qp"
+)
+
+# LP and QP criteria typically give very similar coefficients
+cbind(
+ LP = coef(fit_det, which = "meta"),
+ QP = coef(fit_qp, which = "meta")
+)
+#> LP QP
+#> (Intercept) 1.0416500 1.0416500
+#> log_x1 0.4950084 0.4950084
+#> log_x2 0.1926936 0.1926936fit_sto <- metafrontier(
- log_y ~ log_x1 + log_x2,
- data = sim$data,
- group = "group",
- meta_type = "stochastic"
-)
-
-summary(fit_sto)
-#>
-#> Metafrontier Model Summary
-#> ==========================
-#>
-#> Call:
-#> metafrontier(formula = log_y ~ log_x1 + log_x2, data = sim$data,
-#> group = "group", meta_type = "stochastic")
-#>
-#> Method: sfa
-#> Metafrontier: stochastic
-#>
-#> --- Group: G1 (n = 300) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 1.041650 0.054009 19.29 <2e-16 ***
-#> log_x1 0.495008 0.009173 53.97 <2e-16 ***
-#> log_x2 0.192694 0.008907 21.64 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: 22.685
-#>
-#> --- Group: G2 (n = 300) ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 0.58144 0.05155 11.28 <2e-16 ***
-#> log_x1 0.52876 0.01019 51.88 <2e-16 ***
-#> log_x2 0.20155 0.01037 19.44 <2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: -36.696
-#>
-#> --- Metafrontier ---
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 0.813734 0.208469 3.903 9.49e-05 ***
-#> log_x1 0.513997 0.005155 99.716 < 2e-16 ***
-#> log_x2 0.197514 0.004989 39.588 < 2e-16 ***
-#> ---
-#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-#> Log-likelihood: 181.35
-#>
-#> --- Efficiency Decomposition ---
-#> Group Mean_TE Mean_TGR Mean_TE_star
-#> G1 0.8433 1.1839 0.9984
-#> G2 0.7328 0.8306 0.6087
-#>
-#> --- Technology Gap Ratio Summary ---
-#> Group N Mean SD Min Q1 Median Q3 Max
-#> G1 300 1.1839 0.0320 1.1222 1.1576 1.1861 1.2087 1.2464
-#> G2 300 0.8306 0.0181 0.7949 0.8148 0.8300 0.8445 0.8693The stochastic metafrontier provides standard errors:
-# Variance-covariance matrix
-vcov(fit_sto)
-#> (Intercept) log_x1 log_x2
-#> (Intercept) 4.345943e-02 -8.269764e-05 -6.986322e-05
-#> log_x1 -8.269764e-05 2.657002e-05 1.867196e-06
-#> log_x2 -6.986322e-05 1.867196e-06 2.489219e-05
+fit_sto <- metafrontier(
+ log_y ~ log_x1 + log_x2,
+ data = sim$data,
+ group = "group",
+ meta_type = "stochastic"
+)
-# Log-likelihood of the metafrontier model
-logLik(fit_sto)
-#> 'log Lik.' 181.3504 (df=3)
+summary(fit_sto)
+#>
+#> Metafrontier Model Summary
+#> ==========================
+#>
+#> Call:
+#> metafrontier(formula = log_y ~ log_x1 + log_x2, data = sim$data,
+#> group = "group", meta_type = "stochastic")
+#>
+#> Method: sfa
+#> Metafrontier: stochastic
+#>
+#> --- Group: G1 (n = 300) ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 1.041650 0.054009 19.287 < 2e-16 ***
+#> log_x1 0.495008 0.009173 53.966 < 2e-16 ***
+#> log_x2 0.192694 0.008907 21.635 < 2e-16 ***
+#> log_sigma_v -1.698990 0.114126 -14.887 < 2e-16 ***
+#> log_sigma_u -1.524665 0.228427 -6.675 2.48e-11 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: 22.685
+#>
+#> --- Group: G2 (n = 300) ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 0.58144 0.05155 11.28 <2e-16 ***
+#> log_x1 0.52876 0.01019 51.88 <2e-16 ***
+#> log_x2 0.20155 0.01037 19.44 <2e-16 ***
+#> log_sigma_v -2.00652 0.16307 -12.30 <2e-16 ***
+#> log_sigma_u -0.87260 0.08019 -10.88 <2e-16 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: -36.696
+#>
+#> --- Metafrontier ---
+#> Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) 0.813734 0.208469 3.903 9.49e-05 ***
+#> log_x1 0.513997 0.005155 99.716 < 2e-16 ***
+#> log_x2 0.197514 0.004989 39.588 < 2e-16 ***
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> Log-likelihood: 181.35
+#>
+#> --- Efficiency Decomposition ---
+#> Group Mean_TE Mean_TGR Mean_TE_star
+#> G1 0.8477 1.1839 1.0036
+#> G2 0.7372 0.8306 0.6124
+#>
+#> --- Technology Gap Ratio Summary ---
+#> Group N Mean SD Min Q1 Median Q3 Max
+#> G1 300 1.1839 0.0320 1.1222 1.1576 1.1861 1.2087 1.2464
+#> G2 300 0.8306 0.0181 0.7949 0.8148 0.8300 0.8445 0.8693
+#>
+#> --- Convergence ---
+#> All estimation stages converged.The stochastic metafrontier provides standard errors:
+# Variance-covariance matrix
+vcov(fit_sto)
+#> (Intercept) log_x1 log_x2
+#> (Intercept) 4.345943e-02 -8.269764e-05 -6.986322e-05
+#> log_x1 -8.269764e-05 2.657002e-05 1.867196e-06
+#> log_x2 -6.986322e-05 1.867196e-06 2.489219e-05
+
+# Log-likelihood of the metafrontier model
+logLik(fit_sto)
+#> 'log Lik.' 181.3504 (df=3)"drs" / "irs" (decreasing / increasing
returns)# CRS metafrontier
-fit_crs <- metafrontier(
- log_y ~ log_x1 + log_x2,
- data = sim$data,
- group = "group",
- method = "dea",
- rts = "crs"
-)
-
-# VRS metafrontier
-fit_vrs <- metafrontier(
- log_y ~ log_x1 + log_x2,
- data = sim$data,
- group = "group",
- method = "dea",
- rts = "vrs"
-)
-
-# Compare mean TGR
-cbind(
- CRS = tapply(fit_crs$tgr, fit_crs$group_vec, mean),
- VRS = tapply(fit_vrs$tgr, fit_vrs$group_vec, mean)
-)
-#> CRS VRS
-#> G1 1.0000000 1.0000000
-#> G2 0.6014273 0.8180669# CRS metafrontier
+fit_crs <- metafrontier(
+ log_y ~ log_x1 + log_x2,
+ data = sim$data,
+ group = "group",
+ method = "dea",
+ rts = "crs"
+)
+
+# VRS metafrontier
+fit_vrs <- metafrontier(
+ log_y ~ log_x1 + log_x2,
+ data = sim$data,
+ group = "group",
+ method = "dea",
+ rts = "vrs"
+)
+
+# Compare mean TGR
+cbind(
+ CRS = tapply(fit_crs$tgr, fit_crs$group_vec, mean),
+ VRS = tapply(fit_vrs$tgr, fit_vrs$group_vec, mean)
+)
+#> CRS VRS
+#> G1 1.0000000 1.0000000
+#> G2 0.6014273 0.8180669# Compare TGR estimates across methods
-tgr_det <- tapply(fit_det$tgr, fit_det$group_vec, mean)
-tgr_sto <- tapply(fit_sto$tgr, fit_sto$group_vec, mean)
-tgr_dea <- tapply(fit_crs$tgr, fit_crs$group_vec, mean)
-true_tgr <- tapply(sim$data$true_tgr, sim$data$group, mean)
-
-comparison <- data.frame(
- True = true_tgr,
- Deterministic = tgr_det,
- Stochastic = tgr_sto,
- DEA_CRS = tgr_dea
-)
-round(comparison, 4)
-#> True Deterministic Stochastic DEA_CRS
-#> G1 1.0000 1.000 1.1839 1.0000
-#> G2 0.6703 0.702 0.8306 0.6014# Compare TGR estimates across methods
+tgr_det <- tapply(fit_det$tgr, fit_det$group_vec, mean)
+tgr_sto <- tapply(fit_sto$tgr, fit_sto$group_vec, mean)
+tgr_dea <- tapply(fit_crs$tgr, fit_crs$group_vec, mean)
+true_tgr <- tapply(sim$data$true_tgr, sim$data$group, mean)
+
+comparison <- data.frame(
+ True = true_tgr,
+ Deterministic = tgr_det,
+ Stochastic = tgr_sto,
+ DEA_CRS = tgr_dea
+)
+round(comparison, 4)
+#> True Deterministic Stochastic DEA_CRS
+#> G1 1.0000 1.000 1.1839 1.0000
+#> G2 0.6703 0.702 0.8306 0.6014poolability_test(fit_det)
-#>
-#> Likelihood Ratio Test for Poolability of Group Frontiers
-#>
-#> data: metafrontier(formula = log_y ~ log_x1 + log_x2, data = sim$data, group = "group", meta_type = "deterministic")
-#> LR = 442.35, df = 5, p-value < 2.2e-16A significant test (p < 0.05) confirms that the groups operate -under different technologies and the metafrontier decomposition is -warranted.
+poolability_test(fit_det)
+#>
+#> Likelihood Ratio Test for Poolability of Group Frontiers
+#>
+#> data: fit_det
+#> LR = 442.35, df = 5, p-value < 2.2e-16A significant test (p < 0.05) confirms that the group frontiers +differ, that is, the groups face different restrictions of the common +metatechnology, and the metafrontier decomposition is warranted.
The simulate_metafrontier() function generates data from
a known DGP, enabling parameter recovery studies:
# Monte Carlo: check parameter recovery over 100 replications
-set.seed(1)
-n_rep <- 100
-beta_hat <- matrix(NA, n_rep, 3)
-
-for (r in seq_len(n_rep)) {
- sim_r <- simulate_metafrontier(
- n_groups = 2, n_per_group = 200,
- tech_gap = c(0, 0.3),
- sigma_u = c(0.2, 0.3),
- sigma_v = 0.15
- )
- fit_r <- metafrontier(
- log_y ~ log_x1 + log_x2,
- data = sim_r$data,
- group = "group",
- meta_type = "deterministic"
- )
- beta_hat[r, ] <- coef(fit_r, which = "meta")
-}
-
-# Bias
-true_beta <- c(1.0, 0.5, 0.3)
-colMeans(beta_hat) - true_beta# Monte Carlo: check parameter recovery over 100 replications
+set.seed(1)
+n_rep <- 100
+beta_hat <- matrix(NA, n_rep, 3)
+
+for (r in seq_len(n_rep)) {
+ sim_r <- simulate_metafrontier(
+ n_groups = 2, n_per_group = 200,
+ tech_gap = c(0, 0.3),
+ sigma_u = c(0.2, 0.3),
+ sigma_v = 0.15
+ )
+ fit_r <- metafrontier(
+ log_y ~ log_x1 + log_x2,
+ data = sim_r$data,
+ group = "group",
+ meta_type = "deterministic"
+ )
+ beta_hat[r, ] <- coef(fit_r, which = "meta")
+}
+
+# Bias
+true_beta <- c(1.0, 0.5, 0.3)
+colMeans(beta_hat) - true_betaThe simulate_metafrontier() function supports:
n_groups)Battese, G.E. and Coelli, T.J. (1988). Prediction of firm-level +technical efficiencies with a generalized frontier production function +and panel data. Journal of Econometrics, 38(3), +387–399.
Battese, G.E., Rao, D.S.P. and O’Donnell, C.J. (2004). A metafrontier production function for estimation of technical efficiencies and technology gaps for firms operating under different diff --git a/man/as_metafrontier_model.Rd b/man/as_metafrontier_model.Rd index 1610f6c..57ce9c1 100644 --- a/man/as_metafrontier_model.Rd +++ b/man/as_metafrontier_model.Rd @@ -12,17 +12,34 @@ as_metafrontier_model(x, ...) \item{...}{additional arguments passed to methods.} } \value{ -A list with components: \code{coefficients}, \code{efficiency}, - \code{X}, \code{y}, \code{sigma_v}, \code{sigma_u}, \code{logLik}, - \code{hessian}, \code{n}, \code{dist}. +A list of class \code{"metafrontier_model"} with components: + \code{beta}, \code{te}, \code{X}, \code{y}, \code{sigma_v}, + \code{sigma_u}, \code{logLik}, \code{hessian}, \code{n}, + \code{dist}. } \description{ Generic function that extracts the components needed by \code{\link{metafrontier}} from a pre-fitted frontier model. -Methods are provided for \pkg{sfaR}, \pkg{frontier}, and -\pkg{Benchmarking} objects, as well as plain lists with the +Methods are provided for \pkg{sfaR} (\code{"sfacross"}), +\pkg{frontier} (\code{"frontier"}), and \pkg{Benchmarking} +(\code{"Farrell"}) objects, as well as plain lists with the required fields. } +\details{ +\code{metafrontier(models = ...)} calls this function internally on +each supplied model, so fitted \pkg{sfaR} or \pkg{frontier} objects +can be passed to \code{metafrontier()} directly. Manual conversion +is only needed for hand-built list models. Converting an object that +has already been converted is a no-op, so it is safe to pass +converted objects to \code{metafrontier()} as well. + +Note that \code{Benchmarking::dea()} (\code{"Farrell"}) objects do +not store the inputs, outputs, or frontier coefficients, so the +converted model carries only efficiency scores and cannot be used +with \code{metafrontier(models = ...)}; use the formula interface +with \code{method = "dea"} instead. Converting a Farrell object +therefore raises a warning. +} \examples{ # Using a named list: mod <- as_metafrontier_model(list( diff --git a/man/check_convergence.Rd b/man/check_convergence.Rd new file mode 100644 index 0000000..82e0658 --- /dev/null +++ b/man/check_convergence.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/diagnostics.R +\name{check_convergence} +\alias{check_convergence} +\alias{check_convergence.metafrontier} +\alias{check_convergence.lc_metafrontier} +\alias{check_convergence.malmquist_meta} +\title{Check Convergence of All Estimation Stages} +\usage{ +check_convergence(object, ...) + +\method{check_convergence}{metafrontier}(object, ...) + +\method{check_convergence}{lc_metafrontier}(object, ...) + +\method{check_convergence}{malmquist_meta}(object, ...) +} +\arguments{ +\item{object}{a fitted model object.} + +\item{...}{additional arguments passed to methods.} +} +\value{ +A data frame of class \code{"metafrontier_convergence"} + with one row per estimation stage and columns: + \describe{ + \item{stage}{stage label, e.g. \code{"group: G1"} or + \code{"metafrontier"}} + \item{method}{how the stage was estimated: \code{"MLE"}, + \code{"LP"}, \code{"QP"}, \code{"QP (barrier)"}, + \code{"DEA"}, or \code{"external"}} + \item{code}{the integer convergence code returned by the + optimiser (0 indicates success); \code{NA} for DEA stages + and externally fitted groups} + \item{converged}{logical convergence indicator; for DEA + stages \code{TRUE} unless any efficiency score is + \code{NA} (infeasible programme); \code{NA} for externally + fitted groups} + \item{note}{additional detail, e.g. the number of infeasible + DEA programmes} + } +} +\description{ +Reports the convergence status of every estimation stage of a +fitted metafrontier model: each group-level frontier and the +metafrontier itself. This makes it easy to verify that all +optimisers (MLE) and mathematical programmes (LP/QP) finished +successfully before interpreting technology gap ratios, +confidence intervals, or efficiency decompositions. +} +\details{ +For groups supplied via the \code{models} argument of +\code{\link{metafrontier}} the convergence status of the external +fitter is not tracked, and the corresponding rows carry +\code{NA} with an explanatory note. +} +\examples{ +sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, seed = 42) +fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data, + group = "group") +check_convergence(fit) + +} diff --git a/man/coef.metafrontier.Rd b/man/coef.metafrontier.Rd new file mode 100644 index 0000000..ecb4218 --- /dev/null +++ b/man/coef.metafrontier.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/methods.R +\name{coef.metafrontier} +\alias{coef.metafrontier} +\title{Extract Coefficients from a Metafrontier Model} +\usage{ +\method{coef}{metafrontier}(object, which = c("meta", "group"), extraPar = FALSE, ...) +} +\arguments{ +\item{object}{a \code{"metafrontier"} object.} + +\item{which}{character. \code{"meta"} (default) returns the +metafrontier coefficients; \code{"group"} returns a named list +of group-specific coefficient vectors.} + +\item{extraPar}{logical. If \code{TRUE}, auxiliary parameters are +included alongside the frontier coefficients. For +\code{which = "group"} the variance parameters are +back-transformed to their natural scale (\code{sigmaV}, +\code{sigmaU}), \code{mu} and \code{eta} are kept as estimated, +and heteroscedastic Z coefficients are labelled with their +column names. For \code{which = "meta"} the Stage 2 variance +parameters are appended for stochastic metafrontiers.} + +\item{...}{additional arguments (currently unused).} +} +\value{ +A named numeric vector (\code{which = "meta"}) or a named + list of numeric vectors (\code{which = "group"}). +} +\description{ +Extract Coefficients from a Metafrontier Model +} +\examples{ +sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, seed = 42) +fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data, + group = "group", meta_type = "stochastic") +coef(fit) +coef(fit, extraPar = TRUE) +coef(fit, which = "group", extraPar = TRUE) + +} diff --git a/man/efficiencies.Rd b/man/efficiencies.Rd index ccc74c8..e3ee0a2 100644 --- a/man/efficiencies.Rd +++ b/man/efficiencies.Rd @@ -7,7 +7,7 @@ \usage{ efficiencies(object, ...) -\method{efficiencies}{metafrontier}(object, type = c("meta", "group", "tgr"), ...) +\method{efficiencies}{metafrontier}(object, type = c("meta", "group", "tgr"), estimator = NULL, ...) } \arguments{ \item{object}{a fitted \code{"metafrontier"} object.} @@ -18,6 +18,17 @@ efficiencies(object, ...) \code{"group"} for efficiency relative to the group frontier, \code{"meta"} (default) for efficiency relative to the metafrontier, or \code{"tgr"} for the technology gap ratio.} + +\item{estimator}{optional character. Override the efficiency +estimator used at fit time: \code{"bc88"} for the Battese-Coelli +(1988) conditional expectation \eqn{E[\exp(-u)|\varepsilon]} or +\code{"jlms"} for \eqn{\exp(-E[u|\varepsilon])} (Jondrow et al., +1982). Both are stored on SFA fits, so no refitting is needed; +\code{type = "meta"} is recomputed as \eqn{TE \times TGR}. The +TGR itself does not depend on the estimator. Ignored (with a +warning) for DEA fits and externally fitted group models that do +not carry both estimators. Default \code{NULL} returns the +scores selected at fit time.} } \value{ A numeric vector of efficiency scores of length diff --git a/man/latent_class_metafrontier.Rd b/man/latent_class_metafrontier.Rd index ae08401..dfe0ce9 100644 --- a/man/latent_class_metafrontier.Rd +++ b/man/latent_class_metafrontier.Rd @@ -64,6 +64,12 @@ probabilities, class-specific frontier parameters, and the metafrontier. } \details{ +Latent class estimation is available for SFA-based metafrontiers +only: the EM posterior class probabilities require a parametric +observation-level likelihood, which DEA does not provide. For DEA +fits with observed groups, see \code{\link{poolability_test}}. + + The EM algorithm iterates between: \itemize{ \item \strong{E-step}: compute posterior class membership probabilities diff --git a/man/malmquist_meta.Rd b/man/malmquist_meta.Rd index 6edb126..00c1c61 100644 --- a/man/malmquist_meta.Rd +++ b/man/malmquist_meta.Rd @@ -9,10 +9,12 @@ malmquist_meta( data = NULL, group = NULL, time = NULL, + id = NULL, method = c("dea", "sfa"), dist = c("hnormal", "tnormal", "exponential"), + estimator = c("bc88", "jlms"), orientation = c("output", "input"), - rts = c("crs", "vrs", "drs", "irs"), + rts = c("crs", "vrs", "drs", "irs", "fdh"), control = list(), ... ) @@ -33,20 +35,35 @@ indicators.} that identifies time periods, or a vector of time indicators. Periods must be consecutive integers or sortable.} +\item{id}{optional. A character string naming the column in +\code{data} that identifies firms across periods, or a vector +of firm identifiers. When supplied, firms are matched across +consecutive periods by identifier within each group. When +\code{NULL} (default), firms are matched by row position within +each group, which is valid only for balanced panels sorted +identically in every period (see Details).} + \item{method}{character. \code{"dea"} (default) for DEA-based distance functions or \code{"sfa"} for SFA-based parametric -distance functions.} +distance functions (an approximation; see Details).} \item{dist}{character. Distribution of the inefficiency term when \code{method = "sfa"}: \code{"hnormal"} (default), \code{"tnormal"}, or \code{"exponential"}.} +\item{estimator}{character. Technical efficiency estimator used +when \code{method = "sfa"}: \code{"bc88"} (default) for the +Battese and Coelli (1988) estimator +\eqn{E[\exp(-u)|\varepsilon]}, or \code{"jlms"} for the Jondrow +et al. (1982) estimator \eqn{\exp(-E[u|\varepsilon])}. Passed to +the group SFA fitter.} + \item{orientation}{character. \code{"output"} (default) or \code{"input"}.} \item{rts}{character. Returns to scale assumption: -\code{"crs"} (default), \code{"vrs"}, \code{"drs"}, or -\code{"irs"}.} +\code{"crs"} (default), \code{"vrs"}, \code{"drs"}, +\code{"irs"}, or \code{"fdh"}.} \item{control}{a list of control parameters for the SFA optimiser.} @@ -61,7 +78,9 @@ An object of class \code{"malmquist_meta"}, a list \code{MPI} (metafrontier Malmquist TFP index), \code{TEC} (technical efficiency change), \code{TGC} (technology gap change), - \code{TC} (metafrontier technical change)} + \code{TC} (metafrontier technical change). The \code{id} + column holds the supplied firm identifiers when \code{id} + is given, and the within-group match position otherwise.} \item{group_malmquist}{data frame with the within-group Malmquist index decomposition: \code{MPI_group}, \code{EC_group}, \code{TC_group}} @@ -75,10 +94,17 @@ An object of class \code{"malmquist_meta"}, a list end period), and \code{TGC} (technology gap change, \code{TGR_to / TGR_from})} \item{call}{the matched function call} + \item{method}{the estimation method used (\code{"dea"} or + \code{"sfa"})} \item{orientation}{the orientation used} \item{rts}{the returns to scale assumption} \item{groups}{group labels} \item{periods}{time periods} + \item{n_infeasible}{total number of infeasible cross-period + DEA programs (always \code{0} for \code{method = "sfa"})} + \item{infeasible_by_period}{data frame with the number of + infeasible cross-period DEA programs per period pair + (\code{method = "dea"} only)} } } \description{ @@ -105,18 +131,54 @@ where: the shift of the global production possibility frontier } -Computation uses DEA-based distance functions. For each +\strong{Firm matching:} when \code{id} is supplied, firms are +matched across consecutive periods by identifier within each +technology group. Duplicated (id, period) combinations within a +group are an error. Observations without a within-group match in +the adjacent period, either because the panel is unbalanced or +because a firm switches group between periods, are dropped, and a +single consolidated warning reports the number dropped per period +pair. When \code{id} is \code{NULL}, firms are matched by row +position within each group; this is valid only for balanced +panels sorted identically in every period, so a message is +emitted as a reminder, and a warning is issued when group sizes +differ across a period pair (the unmatched observations are +dropped). Supplying \code{id} is recommended. + +\strong{DEA-based computation (\code{method = "dea"}):} for each consecutive pair of periods \eqn{(s, t)}, eight sets of LP -problems are solved: within-group and pooled efficiencies at -each period, plus cross-period evaluations for the geometric -mean formulation of technical change. - -\strong{Balanced panel assumption:} Firms are matched across -periods by position within each group. The data should contain -a balanced panel (the same firms observed in every period) with -consistent ordering. If group sizes differ across periods, -only the first \code{min(n_s, n_t)} firms per group are paired -and unmatched observations are silently dropped. +problems are solved: within-group and pooled efficiencies at each +period, plus cross-period evaluations for the geometric mean +formulation of technical change. Distances to the metafrontier +are exact distances to the pooled-data frontier, as in O'Donnell, +Rao and Battese (2008). + +\strong{SFA-based computation is an approximation +(\code{method = "sfa"}):} period-specific group SFA frontiers are +estimated, and each observation's metafrontier distance is +approximated by the pointwise maximum of the estimated group +frontier functions evaluated at its inputs; no enveloping +metafrontier is re-estimated. This coincides with the O'Donnell +et al. (2008) metafrontier wherever a single group frontier +dominates, but can understate the metafrontier where group +frontiers cross, which affects TGC and TC*. Prefer +\code{method = "dea"} when an exact decomposition is required. + +\strong{Infeasible cross-period programs:} under +\code{rts = "vrs"}, \code{"drs"}, \code{"irs"}, or \code{"fdh"}, +cross-period LPs can be genuinely infeasible because the +reference technology cannot reach the evaluated observation. Such +cases yield \code{NA} (never \code{Inf}), are excluded from the +reported means, and are counted in a single consolidated warning; +the counts are stored in the \code{n_infeasible} and +\code{infeasible_by_period} components. \code{rts = "crs"} avoids +the issue, as does the hyperbolic orientation available in +\code{\link{metafrontier}}. + +Note that the standard Malmquist index is not a \sQuote{proper} +(multiplicatively complete and transitive) TFP index in the sense +of O'Donnell (2012), so chained comparisons of index levels +across more than two periods should be avoided. } \examples{ # Simulate panel data for 2 groups, 3 time periods @@ -134,12 +196,13 @@ panels <- lapply(1:3, function(t) { }) panel_data <- do.call(rbind, panels) -# Compute metafrontier Malmquist index +# Compute metafrontier Malmquist index, matching firms by id malm <- malmquist_meta( log_y ~ log_x1 + log_x2, data = panel_data, group = "group", - time = "time" + time = "time", + id = "id" ) summary(malm) @@ -149,4 +212,9 @@ O'Donnell, C.J., Rao, D.S.P. and Battese, G.E. (2008). Metafrontier frameworks for the study of firm-level efficiencies and technology ratios. \emph{Empirical Economics}, 34(2), 231--255. \doi{10.1007/s00181-007-0119-4} + +O'Donnell, C.J. (2012). An aggregate quantity framework for +measuring and decomposing productivity change. +\emph{Journal of Productivity Analysis}, 38(3), 255--272. +\doi{10.1007/s11123-012-0275-1} } diff --git a/man/metafrontier-package.Rd b/man/metafrontier-package.Rd index 5511f8b..5dced0c 100644 --- a/man/metafrontier-package.Rd +++ b/man/metafrontier-package.Rd @@ -5,7 +5,7 @@ \alias{metafrontier-package} \title{metafrontier: Analysis of Metafrontier Models for Efficiency and Productivity} \description{ -Implements metafrontier production function models for estimating technical efficiencies and technology gaps for firms operating under different technologies. Supports both stochastic frontier analysis (SFA) and data envelopment analysis (DEA) based metafrontiers. Includes the deterministic metafrontier of Battese, Rao, and O'Donnell (2004) \doi{10.1023/B:PROD.0000012454.06094.29}, the stochastic metafrontier of Huang, Huang, and Liu (2014) \doi{10.1007/s11123-014-0402-2}, and the metafrontier Malmquist productivity index of O'Donnell, Rao, and Battese (2008) \doi{10.1007/s00181-007-0119-4}. Additional features include panel SFA with time-varying inefficiency, bootstrap confidence intervals for technology gap ratios, latent class metafrontier estimation via the EM algorithm, Murphy-Topel corrected standard errors, and 'ggplot2' visualisation methods. +Implements metafrontier production function models for estimating technical efficiencies and technology gaps for groups of firms that face different restrictions of a common underlying metatechnology (group-specific technologies in the sense of Battese, Rao, and O'Donnell, 2004). Supports both stochastic frontier analysis (SFA) and data envelopment analysis (DEA) based metafrontiers. Includes the deterministic metafrontier of Battese, Rao, and O'Donnell (2004) \doi{10.1023/B:PROD.0000012454.06094.29}, the stochastic metafrontier of Huang, Huang, and Liu (2014) \doi{10.1007/s11123-014-0402-2}, and the metafrontier Malmquist productivity index of O'Donnell, Rao, and Battese (2008) \doi{10.1007/s00181-007-0119-4}. The deterministic metafrontier can be identified by either the minimum sum of absolute deviations (LP) or the minimum sum of squared deviations (QP) criterion. Additional features include panel SFA with time-varying inefficiency, bootstrap confidence intervals for technology gap ratios, a DEA poolability permutation test, latent class metafrontier estimation via the EM algorithm, Murphy-Topel corrected standard errors, convergence diagnostics, import of pre-fitted models from external estimation engines ('sfaR', 'frontier', 'Benchmarking'), and 'ggplot2' visualisation methods. } \seealso{ Useful links: diff --git a/man/metafrontier.Rd b/man/metafrontier.Rd index 94ab89b..cc4e58b 100644 --- a/man/metafrontier.Rd +++ b/man/metafrontier.Rd @@ -12,13 +12,17 @@ metafrontier( meta_type = c("deterministic", "stochastic"), dist = c("hnormal", "tnormal", "exponential"), orientation = c("output", "input"), - rts = c("crs", "vrs", "drs", "irs"), + rts = c("crs", "vrs", "drs", "irs", "fdh"), models = NULL, panel = NULL, panel_dist = c("bc92", "bc95"), - type = c("radial", "directional"), + type = c("radial", "directional", "hyperbolic"), direction = c("proportional", "output", "input"), control = list(), + estimator = c("bc88", "jlms"), + objective = c("lp", "qp"), + engine = c("internal", "sfaR", "frontier", "Benchmarking"), + slack = FALSE, ... ) } @@ -57,13 +61,20 @@ Ignored when \code{method = "dea"}.} \item{rts}{character. Returns to scale for DEA: \code{"crs"} (constant, default), \code{"vrs"} (variable), \code{"drs"} -(decreasing), or \code{"irs"} (increasing). Ignored when +(decreasing), \code{"irs"} (increasing), or \code{"fdh"} (free +disposable hull, i.e. no convexity). Ignored when \code{method = "sfa"}.} \item{models}{an optional named list of pre-fitted group-specific -frontier models (objects from \pkg{sfaR}, \pkg{frontier}, or -\pkg{Benchmarking}). If provided, \code{formula}, \code{data}, -and \code{group} are ignored.} +frontier models (objects from \pkg{sfaR} or \pkg{frontier}, or +hand-built lists). Fitted model objects are converted +automatically via \code{\link{as_metafrontier_model}}, so no +manual conversion is required (pre-converting is harmless, the +conversion is idempotent). Farrell objects from +\pkg{Benchmarking} store neither coefficients nor data and +cannot be used here; use the formula interface with +\code{method = "dea"} instead. If \code{models} is provided, +\code{formula}, \code{data}, and \code{group} are ignored.} \item{panel}{an optional list with components \code{id} and \code{time} naming the panel identifier and time columns in @@ -76,12 +87,19 @@ or \code{"bc95"} (Battese and Coelli 1995, observation-specific mean). Only used when \code{panel} is non-NULL.} \item{type}{character. For DEA: \code{"radial"} (default) for -standard radial DEA or \code{"directional"} for directional -distance functions.} +standard radial DEA, \code{"directional"} for directional +distance functions, or \code{"hyperbolic"} for hyperbolic +(graph) efficiency, which contracts inputs and expands outputs +simultaneously.} -\item{direction}{character. Direction vector for DDF: -\code{"proportional"} (default), \code{"output"}, or -\code{"input"}. Only used when \code{type = "directional"}.} +\item{direction}{direction vector for DDF. Either a character +preset (\code{"proportional"} (default), \code{"output"}, or +\code{"input"}), a numeric vector of length m + s giving a +common direction (first m elements for inputs, last s for +outputs), or a numeric n x (m + s) matrix of firm-specific +directions. With numeric directions the ratio-based TGR is not +defined; the additive gap (\code{ddf_gap}) is reported instead. +Only used when \code{type = "directional"}.} \item{control}{a named list of control parameters passed to \code{\link[stats]{optim}}. Common options include @@ -89,6 +107,42 @@ distance functions.} \code{reltol} (relative convergence tolerance, default 1e-10), and \code{fnscale} (set to -1 internally for maximisation).} +\item{estimator}{character. Technical efficiency estimator for +SFA models: \code{"bc88"} (default) for the conditional +expectation \eqn{E[\exp(-u)|\varepsilon]} of Battese and Coelli +(1988), which is the consistent estimator of technical +efficiency, or \code{"jlms"} for +\eqn{\exp(-E[u|\varepsilon])} following Jondrow et al. (1982). +Both are stored on the fitted object; see +\code{\link{efficiencies}}. Ignored when \code{method = "dea"}.} + +\item{objective}{character. Identification criterion for the +deterministic metafrontier: \code{"lp"} (default) minimises the +sum of absolute deviations (a linear programme), \code{"qp"} +minimises the sum of squared deviations (a quadratic programme, +solved exactly via \pkg{quadprog} when available). Both criteria +are proposed in Battese, Rao, and O'Donnell (2004). Only used +when \code{method = "sfa"} and \code{meta_type = +"deterministic"}.} + +\item{engine}{character. Estimation backend for the group +frontiers: \code{"internal"} (default) uses the package's own +estimators; \code{"sfaR"} or \code{"frontier"} delegate the SFA +group frontiers to \code{\link[sfaR]{sfacross}} or +\code{\link[frontier]{sfa}} (cross-sectional, single-part +formulas only); \code{"Benchmarking"} delegates the DEA group +frontiers and the pooled metafrontier to +\code{\link[Benchmarking]{dea}} (radial only), using its +\code{XREF}/\code{YREF} external-reference facility for the +metafrontier stage. The metafrontier stage for SFA methods is +always estimated internally (the Murphy-Topel correction +requires the internal likelihood).} + +\item{slack}{logical. For radial DEA, compute second-stage input +and output slacks (with the radial score held fixed) against +both the group and the pooled reference sets. Default +\code{FALSE}.} + \item{...}{additional arguments passed to the group-level estimation functions.} } @@ -109,7 +163,15 @@ An object of class \code{"metafrontier"} (with subclass \item{groups}{group labels} \item{method}{estimation method used} \item{meta_type}{metafrontier type used} - \item{convergence}{convergence status} + \item{meta_convergence}{integer convergence code for the + metafrontier stage (0 = success; \code{\link[stats]{optim}} + codes for the stochastic metafrontier and the QP barrier + fallback; 0 for a successful LP or DEA solution). Each SFA + group model in \code{group_models} additionally carries its + own \code{convergence} code. Use + \code{\link{check_convergence}} to inspect all stages.} + \item{estimator, objective, engine, meta_solver}{the + estimation choices used for the fit} } } \description{ @@ -129,17 +191,33 @@ where \eqn{TE_i} is efficiency relative to the group frontier and \eqn{TGR_i} is the technology gap ratio measuring how close the group frontier is to the metafrontier. -The deterministic metafrontier (Battese, Rao, and O'Donnell, 2004) -is estimated by solving a linear program that minimises the total -envelope overshoot subject to the constraint that the metafrontier -envelops all group frontiers. BRO (2004) originally proposed a -constrained least-squares (QP) formulation; the LP yields the -tightest envelope and is solved via \pkg{lpSolveAPI}, with a QP -fallback via \code{constrOptim()} when the LP is infeasible. +The deterministic metafrontier is identified by one of the two +criteria proposed by Battese, Rao, and O'Donnell (2004), subject in +both cases to the constraint that the metafrontier envelops all +group frontiers: minimising the sum of absolute deviations, which +reduces to a linear programme because the envelope constraints +force every deviation to be non-negative (O'Donnell, Rao, and +Battese, 2008, Eqs. 23-25), or minimising the sum of squared +deviations, a convex quadratic programme. The LP +(\code{objective = "lp"}, the default) is solved via +\pkg{lpSolveAPI}; the QP (\code{objective = "qp"}) is solved +exactly via \pkg{quadprog} when available, with an adaptive-barrier +fallback via \code{constrOptim()}. The stochastic metafrontier (Huang, Huang, and Liu, 2014) replaces this with a second-stage SFA, providing a distributional framework for inference on the TGR. +\strong{Convergence and failure handling:} estimation stops with an +error only when no usable estimate exists (for example, when both +the BFGS and Nelder-Mead optimisers fail for a group frontier). +When an optimiser stops at a non-zero convergence code, the fitted +object is returned with a warning and the code is recorded; use +\code{\link{check_convergence}} or \code{summary()} to verify all +estimation stages before interpreting technology gap ratios, +confidence intervals, or productivity decompositions. Infeasible +DEA programmes yield \code{NA} efficiency scores, accompanied by a +warning and counted by \code{\link{check_convergence}}. + \strong{Note on standard errors (stochastic metafrontier):} The stochastic metafrontier is a two-stage estimator. Stage 2 treats the fitted group frontier values as data, so the reported standard diff --git a/man/poolability_test.Rd b/man/poolability_test.Rd index bcdb8d1..e22659a 100644 --- a/man/poolability_test.Rd +++ b/man/poolability_test.Rd @@ -4,19 +4,28 @@ \alias{poolability_test} \title{Test Poolability of Group Frontiers} \usage{ -poolability_test(object, ...) +poolability_test(object, B = 199, seed = NULL, ...) } \arguments{ \item{object}{a fitted \code{"metafrontier"} object with -\code{method = "sfa"}.} +\code{method = "sfa"} or \code{method = "dea"}.} + +\item{B}{integer. Number of permutation replicates for the DEA +permutation test (default 199). Ignored for SFA objects.} + +\item{seed}{integer or \code{NULL}. Random seed for the DEA +permutation test, for reproducibility. Ignored for SFA objects.} \item{...}{additional arguments (currently unused).} } \value{ A list of class \code{"htest"} with components: \describe{ - \item{statistic}{the LR test statistic} - \item{parameter}{degrees of freedom} + \item{statistic}{the test statistic (LR statistic for SFA; the + mean technology gap, \eqn{\bar{S} = \mathrm{mean}(1 - TGR)}, + for DEA)} + \item{parameter}{degrees of freedom (SFA) or the effective + number of permutation replicates (DEA)} \item{p.value}{p-value of the test} \item{method}{description of the test} } @@ -24,11 +33,12 @@ A list of class \code{"htest"} with components: \description{ Tests the null hypothesis that all groups share a common frontier (i.e., the metafrontier coincides with all group frontiers) against -the alternative that group-specific frontiers differ. Uses a -likelihood ratio test. +the alternative that group-specific frontiers differ. For SFA-based +metafrontiers a likelihood ratio test is used; for DEA-based +metafrontiers a permutation test is used. } \details{ -The LR statistic is: +\strong{Likelihood ratio test (SFA).} The LR statistic is: \deqn{LR = -2 [LL_{pooled} - \sum_j LL_j]} where \eqn{LL_{pooled}} is the log-likelihood of the pooled (single frontier) model and \eqn{LL_j} are the group-specific @@ -41,7 +51,27 @@ For \eqn{J} groups each with \eqn{p} frontier parameters plus distributional parameters, this equals \eqn{(J - 1) \times p_{total}} where \eqn{p_{total}} includes frontier coefficients, \eqn{\sigma_v}, and \eqn{\sigma_u} -(and \eqn{\mu} for truncated-normal). +(and \eqn{\mu} for truncated-normal). This test requires a +likelihood and is therefore only available for SFA-based +metafrontiers. + +\strong{Permutation test (DEA).} DEA has no likelihood, so the +poolability hypothesis is assessed by a permutation test. Under the +null of a single pooled technology, group labels are exchangeable: +reassigning observations to groups at random should not +systematically change the distance between the group frontiers and +the metafrontier. The observed statistic is the mean technology gap, +\eqn{S_{obs} = \mathrm{mean}(1 - TGR_i)}, and its null distribution +is approximated by refitting the metafrontier on \code{B} random +permutations of the group labels. The p-value is +\eqn{(1 + \#\{S_b \ge S_{obs}\}) / (B + 1)}, following the +aggregate-efficiency inference logic of Simar and Zelenyuk (2007). +The smoothed subsampling approach of Kneip, Simar, and Wilson (2016) +is the asymptotically rigorous alternative for testing hypotheses in +nonparametric production models; the permutation test offered here +is a computationally simple approximation. The default \code{B = 199} +is a pragmatic choice; p-values have resolution \eqn{1/(B + 1)}, so +increase \code{B} for finer resolution. } \examples{ set.seed(42) @@ -51,4 +81,22 @@ fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data, group = "group") poolability_test(fit) +\donttest{ +# DEA permutation test +fit_dea <- metafrontier(log_y ~ log_x1 + log_x2, + data = sim$data, group = "group", + method = "dea") +poolability_test(fit_dea, B = 99, seed = 1) +} + +} +\references{ +Simar, L. and Zelenyuk, V. (2007). Statistical inference for +aggregates of Farrell-type efficiencies. \emph{Journal of Applied +Econometrics}, 22(7), 1367--1394. \doi{10.1002/jae.991} + +Kneip, A., Simar, L. and Wilson, P.W. (2016). Testing hypotheses +in nonparametric models of production. \emph{Journal of Business & +Economic Statistics}, 34(3), 435--447. +\doi{10.1080/07350015.2015.1049747} } diff --git a/man/print.metafrontier.Rd b/man/print.metafrontier.Rd index eb63c88..d19ad86 100644 --- a/man/print.metafrontier.Rd +++ b/man/print.metafrontier.Rd @@ -15,7 +15,11 @@ Invisibly returns \code{x}. } \description{ -Print a Metafrontier Object +Prints a compact overview of a fitted metafrontier model: the +estimation method and metafrontier type, the efficiency estimator +and identification objective (where applicable), the groups and +their sample sizes, group log-likelihoods, mean technology gap +ratio by group, and a one-line convergence status. } \examples{ sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, seed = 42) diff --git a/man/select_n_classes.Rd b/man/select_n_classes.Rd index bac4d2a..a7a93ae 100644 --- a/man/select_n_classes.Rd +++ b/man/select_n_classes.Rd @@ -14,19 +14,22 @@ select_n_classes(formula, data, n_classes_range = 2:5, ...) \item{n_classes_range}{integer vector of class counts to try.} \item{...}{additional arguments passed to -\code{\link{latent_class_metafrontier}}.} +\code{\link{latent_class_metafrontier}}, notably \code{n_starts}.} } \value{ A data frame with columns \code{n_classes}, \code{BIC}, and \code{marginal_ll}. } \description{ -Select Number of Latent Classes via BIC -} -\details{ -Fits latent class metafrontier models for each value in -\code{n_classes_range} and returns BIC values. The optimal -number of classes minimises BIC. +Fits a latent class metafrontier model for each value in +\code{n_classes_range} and tabulates the Bayesian information +criterion (BIC) and marginal log-likelihood of each fit. The +optimal number of classes minimises BIC. Fits that fail are +silently dropped from the table. Because the EM algorithm can +converge to local optima, the ranking is sensitive to the number +of random starts: pass \code{n_starts} (forwarded to +\code{\link{latent_class_metafrontier}}) and increase it for a +more reliable comparison across class counts. } \examples{ \donttest{ diff --git a/man/simulate_metafrontier.Rd b/man/simulate_metafrontier.Rd index a0318df..6fa274b 100644 --- a/man/simulate_metafrontier.Rd +++ b/man/simulate_metafrontier.Rd @@ -12,7 +12,10 @@ simulate_metafrontier( tech_gap = NULL, sigma_u = NULL, sigma_v = 0.2, - seed = NULL + seed = NULL, + beta_groups = NULL, + input_means = NULL, + input_corr = NULL ) } \arguments{ @@ -26,12 +29,15 @@ groups. If a vector, must be of length \code{n_groups}.} \item{beta_meta}{numeric vector. Metafrontier coefficients (including intercept). Length must be \code{n_inputs + 1}. -Default: \code{c(1.0, 0.5, 0.3)}.} +Default: \code{c(1.0, seq(0.5, 0.2, length.out = n_inputs))}, +i.e. \code{c(1.0, 0.5, 0.2)} for the default two inputs. Ignored +when \code{beta_groups} is supplied.} \item{tech_gap}{numeric vector of length \code{n_groups}. The technology gap for each group, defined as the reduction in the intercept relative to the metafrontier. Default: evenly spaced -from 0 to 0.5.} +from 0 to 0.5. Ignored (with a warning) when \code{beta_groups} +is supplied.} \item{sigma_u}{numeric vector of length \code{n_groups}. Standard deviation of the half-normal inefficiency term for each group. @@ -42,6 +48,26 @@ term. Default: 0.2.} \item{seed}{integer or \code{NULL}. Random seed for reproducibility.} + +\item{beta_groups}{optional group-specific frontier coefficients, +including slopes: either an \code{n_groups} x \code{(n_inputs + 1)} +numeric matrix (one row per group) or a list of \code{n_groups} +numeric vectors of length \code{n_inputs + 1}. When supplied, it +replaces the intercept-shift construction based on +\code{tech_gap}; see Details. Default \code{NULL} (intercept-shift +design).} + +\item{input_means}{optional \code{n_groups} x \code{n_inputs} numeric +matrix of per-group mean log-input levels. When supplied, the +log-inputs for group \code{g} are drawn from a normal distribution +centred at \code{input_means[g, ]}; see Details. Default +\code{NULL} (identical uniform inputs across groups).} + +\item{input_corr}{optional \code{n_inputs} x \code{n_inputs} +correlation matrix for the log-inputs. When supplied, the +log-inputs are drawn from a multivariate normal distribution with +this correlation structure; see Details. Default \code{NULL} +(independent inputs).} } \value{ A list with components: @@ -53,8 +79,38 @@ A list with components: } \description{ Generates synthetic data from a known metafrontier data-generating -process. Useful for Monte Carlo simulations, package testing, and -teaching. +process. Each group frontier lies weakly below the metafrontier, +consistent with groups facing different restrictions of a common +metatechnology (Battese, Rao and O'Donnell, 2004). Useful for Monte +Carlo simulations, package testing, and teaching. +} +\details{ +By default the group frontiers share the metafrontier slopes and +differ only in their intercepts, so the true technology gap ratio +(TGR) is constant within each group and equals +\code{exp(-tech_gap[g])}. When \code{beta_groups} is supplied the +group frontiers may differ in their slopes, in which case no single +log-linear metafrontier envelops all groups: the tightest log-linear +envelope is then a pseudo-true quantity. The returned +\code{true_tgr} is instead computed observation by observation +against the pointwise maximum over the group frontiers, +\eqn{TGR_i = \exp(x_i^ op eta_g - \max_j x_i^ op eta_j)}, which is +guaranteed to lie in (0, 1]. The true group frontier for each firm +is \eqn{x_i^ op eta_g}, \code{true_te} is generated exactly as in the +default design, and \code{true_te_star = true_te * true_tgr}. In +this case \code{params$beta_meta} is \code{NULL} and +\code{params$beta_groups} holds the supplied coefficients. + +By default the log-inputs are drawn i.i.d. from a uniform +distribution on \code{[0, 5]}, identically across groups. Supplying +\code{input_means} and/or \code{input_corr} switches to normal +log-inputs with standard deviation \code{5 / sqrt(12)} (matching the +spread of the uniform draws), centred at \code{input_means[g, ]} +(2.5 for every group and input when \code{input_means} is +\code{NULL}). When \code{input_corr} is supplied the draws are +multivariate normal with that correlation matrix; when it is +\code{NULL} but \code{input_means} is given, the inputs are drawn +independently. } \examples{ sim <- simulate_metafrontier(n_groups = 3, n_per_group = 200, @@ -65,4 +121,10 @@ table(sim$data$group) # The true metafrontier coefficients sim$params$beta_meta +# Group-specific slopes: per-observation true TGR +sim2 <- simulate_metafrontier( + beta_groups = rbind(c(1.0, 0.5, 0.2), c(0.9, 0.6, 0.1)) +) +range(sim2$data$true_tgr) + } diff --git a/man/simulate_panel_metafrontier.Rd b/man/simulate_panel_metafrontier.Rd index beb6a22..d98a4c0 100644 --- a/man/simulate_panel_metafrontier.Rd +++ b/man/simulate_panel_metafrontier.Rd @@ -13,7 +13,8 @@ simulate_panel_metafrontier( sigma_u = 0.3, sigma_v = 0.2, eta = 0.05, - seed = NULL + seed = NULL, + attrition = 0 ) } \arguments{ @@ -34,6 +35,15 @@ simulate_panel_metafrontier( \item{eta}{numeric. Time-decay parameter for BC92.} \item{seed}{integer or NULL. Random seed.} + +\item{attrition}{numeric in [0, 0.5]. Probability that each +firm-period observation after a firm's first period is dropped +independently, producing an unbalanced panel. Every firm's first +period is always kept, so all firms remain in the data. The +attrition draws are made after all other random numbers, so +\code{attrition = 0} (the default) reproduces legacy balanced +datasets exactly for the same seed. The realised share of at-risk +observations dropped is stored in \code{params$attrition_share}.} } \value{ A list with components: @@ -41,7 +51,9 @@ A list with components: \item{data}{data frame with columns: \code{firm}, \code{year}, \code{group}, \code{log_y}, \code{log_x1}, \code{log_x2}, \code{true_te}, \code{true_u}, \code{true_v}} - \item{params}{list of true parameter values used in generation} + \item{params}{list of true parameter values used in generation, + including \code{attrition} and the realised + \code{attrition_share}} } } \description{ @@ -57,4 +69,8 @@ sim <- simulate_panel_metafrontier( head(sim$data) str(sim$params) +# An unbalanced panel with roughly 20\% attrition +sim_unbal <- simulate_panel_metafrontier(seed = 42, attrition = 0.2) +table(table(sim_unbal$data$firm)) + } diff --git a/man/summary.metafrontier.Rd b/man/summary.metafrontier.Rd index eab223c..9611635 100644 --- a/man/summary.metafrontier.Rd +++ b/man/summary.metafrontier.Rd @@ -12,10 +12,44 @@ \item{...}{additional arguments (currently unused).} } \value{ -An object of class \code{"summary.metafrontier"}. +An object of class \code{"summary.metafrontier"}: a list + with components + \describe{ + \item{call}{the matched call of the original fit} + \item{method}{estimation method (\code{"sfa"} or \code{"dea"})} + \item{meta_type}{metafrontier type (\code{"deterministic"} or + \code{"stochastic"})} + \item{groups}{character vector of group labels} + \item{nobs}{named vector of observation counts (total and per + group)} + \item{group_tables}{named list of coefficient matrices, one per + group, with columns \code{Estimate}, \code{Std. Error}, + \code{z value}, and \code{Pr(>|z|)} where standard errors are + available (empty list for DEA fits)} + \item{meta_table}{metafrontier coefficient matrix in the same + format, or \code{NULL} for DEA fits} + \item{tgr_summary}{data frame of TGR statistics by group, as + returned by \code{\link{tgr_summary}}} + \item{efficiency_summary}{data frame with mean TE, mean TGR, + and mean TE* by group} + \item{logLik_groups}{named vector of group log-likelihoods, or + \code{NULL}} + \item{meta_logLik}{Stage 2 log-likelihood of the stochastic + metafrontier, or \code{NULL}} + \item{convergence}{data frame with columns \code{stage}, + \code{code}, and \code{converged} recording the optimiser + status of each estimation stage, or \code{NULL} if + unavailable; see \code{\link{check_convergence}}} + } } \description{ -Summary of a Metafrontier Model +Computes group-level summaries of technical efficiency (TE), +technology gap ratio (TGR), and metafrontier efficiency (TE*), +full coefficient tables for each group frontier (including +variance parameters and, for BC92 panels, \code{eta}, all with +standard errors where a Hessian is available), the metafrontier +coefficient table (with Murphy-Topel corrected standard errors +where applicable), and a per-stage convergence table. } \examples{ sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, seed = 42) diff --git a/man/technology_gap_ratio.Rd b/man/technology_gap_ratio.Rd index a45858b..a2d253e 100644 --- a/man/technology_gap_ratio.Rd +++ b/man/technology_gap_ratio.Rd @@ -39,7 +39,11 @@ metafrontier need not envelop the group frontiers at every point. A TGR of 1 means the group frontier coincides with the metafrontier at that input mix. Values less than 1 indicate a -technology gap. +technology gap. Since each group technology is a restricted subset +of the common metatechnology (Battese, Rao and O'Donnell, 2004), +the gap reflects the restrictions a group faces (regulation, +environment, endowments) rather than a fundamentally different +technology. } \examples{ set.seed(42) diff --git a/man/tgr_summary.Rd b/man/tgr_summary.Rd index bd8f9bf..30ff5ab 100644 --- a/man/tgr_summary.Rd +++ b/man/tgr_summary.Rd @@ -16,5 +16,20 @@ A data frame with columns: Group, N, Mean, SD, Min, Q1, Median, Q3, Max. } \description{ -Prints a summary table of TGR statistics by group. +Prints a summary table of TGR statistics by group. The underlying +observation-level TGR values are the same as those returned by +\code{efficiencies(object, type = "tgr")}. +} +\examples{ +sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, + seed = 42) +fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data, + group = "group", method = "sfa", + meta_type = "deterministic") +tgr_summary(fit) + +} +\seealso{ +\code{\link{efficiencies.metafrontier}}, + \code{\link{technology_gap_ratio}}, \code{\link{boot_tgr}} } diff --git a/man/vcov.metafrontier.Rd b/man/vcov.metafrontier.Rd index db729e5..3aa7445 100644 --- a/man/vcov.metafrontier.Rd +++ b/man/vcov.metafrontier.Rd @@ -4,7 +4,13 @@ \alias{vcov.metafrontier} \title{Variance-Covariance Matrix for Metafrontier Coefficients} \usage{ -\method{vcov}{metafrontier}(object, correction = c("none", "murphy-topel"), ...) +\method{vcov}{metafrontier}( + object, + correction = c("none", "murphy-topel"), + which = c("meta", "group"), + extraPar = FALSE, + ... +) } \arguments{ \item{object}{a \code{"metafrontier"} object.} @@ -15,13 +21,38 @@ the Murphy and Topel (1985) correction for first-stage estimation uncertainty (the generated-regressor problem). Only available for stochastic metafrontiers.} +\item{which}{character. \code{"meta"} (default) returns the +metafrontier (Stage 2) variance-covariance matrix; +\code{"group"} returns a named list with one full +variance-covariance matrix per group (from the inverse negative +Hessian of the group MLE), with \code{NULL} entries for groups +without a stored Hessian.} + +\item{extraPar}{logical. If \code{TRUE} and \code{which = "meta"}, +the full Stage 2 matrix is returned, including the rows and +columns for the auxiliary parameters (raw MLE parameterisation, +e.g. \code{log_sigma_v}); the default returns only the block for +the frontier coefficients.} + \item{...}{additional arguments (currently unused).} } \value{ -A variance-covariance matrix, or \code{NULL} if unavailable. +A variance-covariance matrix (\code{which = "meta"}), a + named list of matrices (\code{which = "group"}), or \code{NULL} + if unavailable. } \description{ -Variance-Covariance Matrix for Metafrontier Coefficients +Returns the variance-covariance matrix of the Stage 2 +(metafrontier) coefficients, or, with \code{which = "group"}, the +per-group matrices from the Stage 1 maximum likelihood fits. +\code{NULL} is returned when no Stage 2 Hessian exists: the +deterministic metafrontier is fitted by LP/QP optimisation and has +no sampling variance in this framework, so \code{which = "meta"} +returns \code{NULL} with a warning; with \code{which = "group"}, +list entries are \code{NULL} for groups without a stored Hessian +(e.g. externally fitted models). DEA-based metafrontiers are +nonparametric and \code{vcov()} signals an error; use +\code{\link{boot_tgr}} for inference instead. } \examples{ sim <- simulate_metafrontier(n_groups = 2, n_per_group = 50, seed = 42) diff --git a/tests/testthat/test-bc88.R b/tests/testthat/test-bc88.R new file mode 100644 index 0000000..42703d0 --- /dev/null +++ b/tests/testthat/test-bc88.R @@ -0,0 +1,88 @@ +# Tests for the BC88 (Battese and Coelli, 1988) efficiency estimator + +.bc88_test_data <- function(n = 200, seed = 1) { + set.seed(seed) + x1 <- rnorm(n) + x2 <- rnorm(n) + u <- abs(rnorm(n, 0, 0.3)) + v <- rnorm(n, 0, 0.2) + data.frame(y = 1 + 0.5 * x1 + 0.3 * x2 + v - u, x1 = x1, x2 = x2) +} + +test_that("bc88 and jlms are both stored, bounded, and highly correlated", { + d <- .bc88_test_data() + f <- Formula::Formula(y ~ x1 + x2) + + gm <- metafrontier:::.fit_sfa_group(f, d, "hnormal", list()) + + expect_true(all(gm$efficiency_bc88 > 0 & gm$efficiency_bc88 < 1)) + expect_true(all(gm$efficiency_jlms > 0 & gm$efficiency_jlms < 1)) + expect_gt(cor(gm$efficiency_bc88, gm$efficiency_jlms), 0.99) +}) + +test_that("bc88 matches the textbook formula recomputed by hand", { + d <- .bc88_test_data() + f <- Formula::Formula(y ~ x1 + x2) + + gm <- metafrontier:::.fit_sfa_group(f, d, "hnormal", list()) + + # Half-normal conditional posterior of u given eps + sv <- gm$sigma_v + su <- gm$sigma_u + eps <- gm$residuals + s2 <- sv^2 + su^2 + mu_star <- -eps * su^2 / s2 + sigma_star <- sv * su / sqrt(s2) + + te_hand <- exp(-mu_star + 0.5 * sigma_star^2) * + pnorm(mu_star / sigma_star - sigma_star) / pnorm(mu_star / sigma_star) + + expect_equal(gm$efficiency_bc88, as.numeric(te_hand), tolerance = 1e-10) +}) + +test_that("bc88 and jlms agree with sfaR::sfacross", { + skip_if_not_installed("sfaR") + + d <- .bc88_test_data() + f <- Formula::Formula(y ~ x1 + x2) + + gm <- metafrontier:::.fit_sfa_group(f, d, "hnormal", list()) + sm <- sfaR::sfacross(y ~ x1 + x2, data = d, udist = "hnormal") + effs <- sfaR::efficiencies(sm) + + expect_gt(cor(gm$efficiency_bc88, effs$teBC), 0.999) + expect_lt(max(abs(gm$efficiency_bc88 - effs$teBC)), 1e-3) + expect_gt(cor(gm$efficiency_jlms, effs$teJLMS), 0.999) + expect_lt(max(abs(gm$efficiency_jlms - effs$teJLMS)), 1e-3) +}) + +test_that("estimator argument selects the efficiency vector, default bc88", { + d <- .bc88_test_data() + f <- Formula::Formula(y ~ x1 + x2) + + gm_def <- metafrontier:::.fit_sfa_group(f, d, "hnormal", list()) + gm_jlms <- metafrontier:::.fit_sfa_group(f, d, "hnormal", list(), + estimator = "jlms") + + expect_identical(gm_def$estimator, "bc88") + expect_identical(gm_def$efficiency, gm_def$efficiency_bc88) + expect_identical(gm_jlms$estimator, "jlms") + expect_identical(gm_jlms$efficiency, gm_jlms$efficiency_jlms) +}) + +test_that("bc88 is computed in all distribution branches", { + d <- .bc88_test_data() + set.seed(2) + d$z1 <- rnorm(nrow(d)) + f <- Formula::Formula(y ~ x1 + x2) + fz <- Formula::Formula(y ~ x1 + x2 | z1) + + for (dist in c("hnormal", "tnormal", "exponential")) { + gm <- metafrontier:::.fit_sfa_group(f, d, dist, list()) + expect_true(all(gm$efficiency_bc88 > 0 & gm$efficiency_bc88 <= 1), + label = paste("homoscedastic", dist)) + gmz <- metafrontier:::.fit_sfa_group(fz, d, dist, list()) + expect_true(all(gmz$efficiency_bc88 > 0 & gmz$efficiency_bc88 <= 1), + label = paste("heteroscedastic", dist)) + } +}) diff --git a/tests/testthat/test-dea-features.R b/tests/testthat/test-dea-features.R new file mode 100644 index 0000000..661bd92 --- /dev/null +++ b/tests/testthat/test-dea-features.R @@ -0,0 +1,371 @@ +# Tests for v0.3.0 DEA features: phi-bound fix, FDH, user-defined +# DDF directions, hyperbolic efficiency, and second-stage slacks + +# Small positive-valued production dataset shared across tests +.make_dea_feature_data <- function(n = 40, seed = 123) { + set.seed(seed) + x1 <- runif(n, 1, 10) + x2 <- runif(n, 1, 10) + grp <- rep(c("A", "B"), each = n / 2) + tech <- ifelse(grp == "A", 1, 0.8) + y <- tech * (x1^0.4 * x2^0.4) * runif(n, 0.6, 1) + data.frame(y = y, x1 = x1, x2 = x2, group = grp) +} + + +# ---- (i) phi lower bound fix ---- + +test_that("output-oriented LP solves for super-efficient DMUs (phi < 1)", { + # Reference technology: y = x under CRS. The evaluated point (1, 2) + # lies strictly above it, so phi* = 0.5 < 1. With the old lower + # bound phi >= 1 this LP was infeasible. + X_ref <- matrix(c(1, 2, 3, 4, 5), ncol = 1) + Y_ref <- matrix(c(1, 2, 3, 4, 5), ncol = 1) + + te <- metafrontier:::.dea_solve_lp(1, 2, X_ref, Y_ref, "output", "crs") + expect_false(is.na(te)) + expect_equal(te, 2, tolerance = 1e-8) # Farrell TE = 1/phi + + # The batch solver takes the same cross-period path + te_b <- metafrontier:::.dea_batch_fast(matrix(1, 1, 1), matrix(2, 1, 1), + "output", "crs", + X_ref = X_ref, Y_ref = Y_ref) + expect_equal(te_b, 2, tolerance = 1e-8) +}) + + +test_that("same-period output-oriented scores remain in (0, 1]", { + dat <- .make_dea_feature_data() + X <- cbind(dat$x1, dat$x2) + Y <- matrix(dat$y, ncol = 1) + + te <- metafrontier:::.dea_batch_fast(X, Y, "output", "crs") + expect_true(all(te > 0)) + expect_true(all(te <= 1 + 1e-8)) +}) + + +# ---- (ii) FDH ---- + +test_that("FDH radial efficiency dominates VRS in both orientations", { + skip_on_cran() + dat <- .make_dea_feature_data() + X <- cbind(dat$x1, dat$x2) + Y <- matrix(dat$y, ncol = 1) + + for (orient in c("input", "output")) { + e_fdh <- metafrontier:::.dea_batch_fast(X, Y, orient, "fdh") + e_vrs <- metafrontier:::.dea_batch_fast(X, Y, orient, "vrs") + expect_true(all(e_fdh >= e_vrs - 1e-8)) + expect_true(all(e_fdh > 0)) + expect_true(all(e_fdh <= 1 + 1e-8)) + } +}) + + +test_that("FDH metafrontier TGR lies in (0, 1]", { + skip_on_cran() + dat <- .make_dea_feature_data() + f <- Formula::Formula(y ~ x1 + x2) + gvec <- factor(dat$group) + glev <- levels(gvec) + + gm <- lapply(glev, function(g) { + metafrontier:::.fit_dea_group(f, dat[gvec == g, ], "input", "fdh") + }) + names(gm) <- glev + + res <- metafrontier:::.estimate_dea_metafrontier(f, dat, gvec, glev, gm, + "input", "fdh") + expect_true(all(res$tgr > 0)) + expect_true(all(res$tgr <= 1 + 1e-8)) +}) + + +test_that("FDH returns NA with a warning when no reference point dominates", { + # Cross-period evaluation: no reference point can produce y = 10 + expect_warning( + v <- metafrontier:::.dea_batch_fast(matrix(0.5, 1, 1), matrix(10, 1, 1), + "input", "fdh", + X_ref = matrix(1:3, ncol = 1), + Y_ref = matrix(1:3, ncol = 1)), + "No dominating FDH reference point" + ) + expect_true(is.na(v)) +}) + + +test_that("DDF with FDH (binary MIP) is no larger than the VRS beta", { + skip_on_cran() + dat <- .make_dea_feature_data(n = 20) + X <- cbind(dat$x1, dat$x2) + Y <- matrix(dat$y, ncol = 1) + + b_fdh <- metafrontier:::.dea_solve_ddf(X[1, ], Y[1, ], X, Y, + g_x = X[1, ], g_y = Y[1, ], + rts = "fdh") + b_vrs <- metafrontier:::.dea_solve_ddf(X[1, ], Y[1, ], X, Y, + g_x = X[1, ], g_y = Y[1, ], + rts = "vrs") + expect_true(b_fdh >= -1e-8) + expect_true(b_fdh <= b_vrs + 1e-8) +}) + + +# ---- (iii) user-defined DDF directions ---- + +test_that("numeric direction vector runs and reports additive gaps", { + skip_on_cran() + dat <- .make_dea_feature_data() + f <- Formula::Formula(y ~ x1 + x2) + gvec <- factor(dat$group) + glev <- levels(gvec) + gdir <- c(mean(dat$x1), mean(dat$x2), mean(dat$y)) + + gm <- lapply(glev, function(g) { + metafrontier:::.fit_ddf_group(f, dat[gvec == g, ], "crs", gdir) + }) + names(gm) <- glev + + res <- metafrontier:::.estimate_ddf_metafrontier(f, dat, gvec, glev, gm, + "crs", gdir) + + # te = 1/(1 + beta) is undefined for arbitrary numeric directions + expect_true(all(is.na(res$te_group))) + expect_true(all(is.na(res$te_meta))) + expect_true(all(is.na(res$tgr))) + + # Additive fields carry the results; the gap is non-negative because + # the pooled reference set is a superset + expect_length(res$beta_group, nrow(dat)) + expect_length(res$beta_meta, nrow(dat)) + expect_true(all(res$beta_meta >= res$beta_group - 1e-8)) + expect_equal(res$ddf_gap, res$beta_meta - res$beta_group) +}) + + +test_that("firm-specific direction matrix runs and matches common vector", { + skip_on_cran() + dat <- .make_dea_feature_data() + f <- Formula::Formula(y ~ x1 + x2) + gvec <- factor(dat$group) + glev <- levels(gvec) + n <- nrow(dat) + gdir <- c(mean(dat$x1), mean(dat$x2), mean(dat$y)) + dmat <- matrix(rep(gdir, n), n, 3, byrow = TRUE) + + gm <- lapply(glev, function(g) { + metafrontier:::.fit_ddf_group(f, dat[gvec == g, ], "crs", + dmat[gvec == g, , drop = FALSE]) + }) + names(gm) <- glev + + res_mat <- metafrontier:::.estimate_ddf_metafrontier(f, dat, gvec, glev, + gm, "crs", dmat) + + gm_vec <- lapply(glev, function(g) { + metafrontier:::.fit_ddf_group(f, dat[gvec == g, ], "crs", gdir) + }) + names(gm_vec) <- glev + res_vec <- metafrontier:::.estimate_ddf_metafrontier(f, dat, gvec, glev, + gm_vec, "crs", gdir) + + expect_equal(res_mat$beta_meta, res_vec$beta_meta, tolerance = 1e-10) + expect_equal(res_mat$beta_group, res_vec$beta_group, tolerance = 1e-10) +}) + + +test_that("character presets still populate te and the additive fields", { + dat <- .make_dea_feature_data(n = 20) + f <- Formula::Formula(y ~ x1 + x2) + gvec <- factor(dat$group) + glev <- levels(gvec) + + gm <- lapply(glev, function(g) { + metafrontier:::.fit_ddf_group(f, dat[gvec == g, ], "crs", "proportional") + }) + names(gm) <- glev + + res <- metafrontier:::.estimate_ddf_metafrontier(f, dat, gvec, glev, gm, + "crs", "proportional") + expect_false(any(is.na(res$te_meta))) + expect_false(any(is.na(res$tgr))) + expect_equal(res$ddf_gap, res$beta_meta - res$beta_group) + expect_true(all(res$ddf_gap >= -1e-8)) +}) + + +test_that("malformed numeric directions are rejected", { + X <- matrix(1:6, 3, 2) + Y <- matrix(1:3, 3, 1) + expect_error(metafrontier:::.ddf_direction_mats(c(1, 2), X, Y), + "length m \\+ s") + expect_error(metafrontier:::.ddf_direction_mats(matrix(1, 2, 3), X, Y), + "n x \\(m \\+ s\\)") + expect_error(metafrontier:::.ddf_direction_mats(c(-1, 1, 1), X, Y), + "non-negative") + expect_error(metafrontier:::.ddf_direction_mats(c(0, 0, 0), X, Y), + "at least one positive") +}) + + +# ---- (iv) hyperbolic efficiency ---- + +test_that("hyperbolic CRS closed form agrees with direct bisection", { + skip_on_cran() + dat <- .make_dea_feature_data(n = 20) + X <- cbind(dat$x1, dat$x2) + Y <- matrix(dat$y, ncol = 1) + n <- nrow(X) + + g_closed <- metafrontier:::.hyperbolic_batch(X, Y, "crs") + expect_true(all(g_closed > 0)) + expect_true(all(g_closed <= 1 + 1e-8)) + + # Independent bisection on CRS feasibility LPs: + # gamma feasible iff exists lambda >= 0 with + # X' lambda <= gamma * x_i and Y' lambda >= y_i / gamma + crs_feasible <- function(gamma, i) { + lp <- lpSolveAPI::make.lp(0, n) + lpSolveAPI::set.objfn(lp, rep(0, n)) + lpSolveAPI::lp.control(lp, sense = "min", verbose = "neutral") + for (mm in 1:2) { + lpSolveAPI::add.constraint(lp, X[, mm], "<=", gamma * X[i, mm]) + } + lpSolveAPI::add.constraint(lp, Y[, 1], ">=", Y[i, 1] / gamma) + lpSolveAPI::solve.lpExtPtr(lp) == 0 + } + + for (i in c(1, 5, 10, 15, 20)) { + lo <- 0 + hi <- 1 + for (iter in 1:40) { + if (hi - lo < 1e-9) break + mid <- (lo + hi) / 2 + if (crs_feasible(mid, i)) hi <- mid else lo <- mid + } + expect_equal(g_closed[i], hi, tolerance = 1e-6) + } +}) + + +test_that("hyperbolic metafrontier works under VRS with TGR <= 1", { + skip_on_cran() + dat <- .make_dea_feature_data() + f <- Formula::Formula(y ~ x1 + x2) + gvec <- factor(dat$group) + glev <- levels(gvec) + + gm <- lapply(glev, function(g) { + metafrontier:::.fit_hyperbolic_group(f, dat[gvec == g, ], "vrs") + }) + names(gm) <- glev + + res <- metafrontier:::.estimate_hyperbolic_metafrontier(f, dat, gvec, + glev, gm, "vrs") + + expect_true(all(res$te_group > 0)) + expect_true(all(res$te_group <= 1 + 1e-8)) + expect_true(all(res$te_meta > 0)) + expect_true(all(res$te_meta <= 1 + 1e-8)) + # Bisection tolerance (1e-8 per gamma) is amplified in the ratio + expect_true(all(res$tgr > 0)) + expect_true(all(res$tgr <= 1 + 1e-6)) + expect_identical(res$meta_convergence, 0L) + expect_null(res$meta_coef) +}) + + +test_that("hyperbolic FDH and VRS envelop CRS", { + skip_on_cran() + dat <- .make_dea_feature_data(n = 20) + X <- cbind(dat$x1, dat$x2) + Y <- matrix(dat$y, ncol = 1) + + g_crs <- metafrontier:::.hyperbolic_batch(X, Y, "crs") + g_vrs <- metafrontier:::.hyperbolic_batch(X, Y, "vrs") + g_fdh <- metafrontier:::.hyperbolic_batch(X, Y, "fdh") + + expect_true(all(g_vrs >= g_crs - 1e-6)) + expect_true(all(g_fdh >= g_vrs - 1e-6)) + expect_true(all(g_fdh <= 1 + 1e-8)) +}) + + +# ---- (v) slacks ---- + +test_that("slack = TRUE returns non-negative slack matrices of correct dims", { + skip_on_cran() + dat <- .make_dea_feature_data() + f <- Formula::Formula(y ~ x1 + x2) + gvec <- factor(dat$group) + glev <- levels(gvec) + + gm <- lapply(glev, function(g) { + metafrontier:::.fit_dea_group(f, dat[gvec == g, ], "input", "vrs", + slack = TRUE) + }) + names(gm) <- glev + + for (g in glev) { + n_g <- gm[[g]]$nobs + expect_identical(dim(gm[[g]]$slack_x), c(n_g, 2L)) + expect_identical(dim(gm[[g]]$slack_y), c(n_g, 1L)) + expect_true(all(gm[[g]]$slack_x >= 0)) + expect_true(all(gm[[g]]$slack_y >= 0)) + } + + res <- metafrontier:::.estimate_dea_metafrontier(f, dat, gvec, glev, gm, + "input", "vrs", + slack = TRUE) + expect_identical(dim(res$slack_x_meta), c(nrow(dat), 2L)) + expect_identical(dim(res$slack_y_meta), c(nrow(dat), 1L)) + expect_true(all(res$slack_x_meta >= 0)) + expect_true(all(res$slack_y_meta >= 0)) +}) + + +test_that("a radially efficient but dominated DMU shows positive slack", { + # Classic example: A = (1, 4), B = (4, 1), C = (2, 2) efficient; + # D = (1, 5) has theta* = 1 under VRS but slack of 1 in input 2 + X <- matrix(c(1, 4, 1, 2, + 4, 1, 5, 2), ncol = 2) + Y <- matrix(1, 4, 1) + + theta <- metafrontier:::.dea_batch_fast(X, Y, "input", "vrs") + expect_equal(theta[3], 1, tolerance = 1e-8) + + sl <- metafrontier:::.dea_slacks(X, Y, theta, "input", "vrs", X, Y) + expect_true(all(sl$slack_x >= 0)) + expect_true(all(sl$slack_y >= 0)) + expect_equal(sl$slack_x[3, 2], 1, tolerance = 1e-6) + expect_equal(sl$slack_x[3, 1], 0, tolerance = 1e-6) +}) + + +test_that("FDH slacks are measured against the dominating peer", { + X <- matrix(c(1, 4, 1, 2, + 4, 1, 5, 2), ncol = 2) + Y <- matrix(1, 4, 1) + + theta <- metafrontier:::.dea_batch_fast(X, Y, "input", "fdh") + sl <- metafrontier:::.dea_slacks(X, Y, theta, "input", "fdh", X, Y) + + expect_true(all(sl$slack_x >= 0)) + expect_true(all(sl$slack_y >= 0)) + # D = (1, 5) is dominated by peer A = (1, 4): slack of 1 in input 2 + expect_equal(sl$slack_x[3, 2], 1, tolerance = 1e-8) +}) + + +test_that("output-oriented slacks are non-negative", { + skip_on_cran() + X <- matrix(c(1, 4, 1, 2, + 4, 1, 5, 2), ncol = 2) + Y <- matrix(1, 4, 1) + + te <- metafrontier:::.dea_batch_fast(X, Y, "output", "vrs") + sl <- metafrontier:::.dea_slacks(X, Y, te, "output", "vrs", X, Y) + expect_true(all(sl$slack_x >= 0)) + expect_true(all(sl$slack_y >= 0)) +}) diff --git a/tests/testthat/test-diagnostics.R b/tests/testthat/test-diagnostics.R new file mode 100644 index 0000000..2fec7d8 --- /dev/null +++ b/tests/testthat/test-diagnostics.R @@ -0,0 +1,172 @@ +# Tests for convergence diagnostics and extended coef/vcov/summary methods + +sim_cs <- simulate_metafrontier(n_groups = 2, n_per_group = 40, seed = 42) +fit_det <- metafrontier(log_y ~ log_x1 + log_x2, data = sim_cs$data, + group = "group") +fit_stoch <- metafrontier(log_y ~ log_x1 + log_x2, data = sim_cs$data, + group = "group", meta_type = "stochastic") + +sim_panel <- simulate_panel_metafrontier(n_groups = 2, + n_firms_per_group = 20, + n_periods = 4, + eta = 0.05, seed = 123) +fit_panel <- metafrontier(log_y ~ log_x1 + log_x2, data = sim_panel$data, + group = "group", + panel = list(id = "firm", time = "year"), + panel_dist = "bc92") + + +test_that("check_convergence reports all stages for a deterministic SFA fit", { + conv <- check_convergence(fit_det) + + expect_s3_class(conv, "metafrontier_convergence") + expect_s3_class(conv, "data.frame") + expect_named(conv, c("stage", "method", "code", "converged", "note")) + expect_equal(nrow(conv), length(fit_det$groups) + 1L) + expect_equal(conv$stage, + c(paste0("group: ", fit_det$groups), "metafrontier")) + expect_true(all(conv$converged)) + expect_true(all(conv$method[seq_along(fit_det$groups)] == "MLE")) + expect_true(conv$method[nrow(conv)] %in% c("LP", "QP", "QP (barrier)")) + expect_true(is.integer(conv$code)) +}) + + +test_that("check_convergence on a stochastic fit reports MLE meta stage", { + conv <- check_convergence(fit_stoch) + expect_equal(conv$method[nrow(conv)], "MLE") + expect_true(all(conv$converged)) +}) + + +test_that("check_convergence on a DEA fit uses method DEA and NA codes", { + fit_dea <- metafrontier(log_y ~ log_x1 + log_x2, data = sim_cs$data, + group = "group", method = "dea", rts = "crs") + conv <- check_convergence(fit_dea) + + expect_true(all(conv$method == "DEA")) + expect_true(all(is.na(conv$code))) + expect_true(all(conv$converged)) + expect_true(all(conv$note == "")) +}) + + +test_that("check_convergence flags externally fitted groups", { + mods <- lapply(fit_det$group_models, function(gm) { + list(coefficients = gm$coefficients, efficiency = gm$efficiency, + X = gm$X, y = gm$y) + }) + fit_ext <- metafrontier(models = mods) + conv <- check_convergence(fit_ext) + + n_grp <- length(fit_ext$groups) + expect_true(all(conv$method[seq_len(n_grp)] == "external")) + expect_true(all(is.na(conv$code[seq_len(n_grp)]))) + expect_true(all(is.na(conv$converged[seq_len(n_grp)]))) + expect_match(conv$note[1], "fitted externally") + + expect_warning( + coef(fit_ext, which = "group", extraPar = TRUE), + "externally" + ) +}) + + +test_that("check_convergence.default errors clearly", { + expect_error(check_convergence(lm(dist ~ speed, data = cars)), + "not implemented") +}) + + +test_that("coef(..., extraPar = TRUE) exposes eta for a BC92 panel fit", { + cf <- coef(fit_panel, which = "group", extraPar = TRUE) + + expect_true(is.list(cf)) + expect_named(cf, fit_panel$groups) + for (g in fit_panel$groups) { + expect_true(all(c("sigmaV", "sigmaU", "eta") %in% names(cf[[g]]))) + expect_gt(unname(cf[[g]]["sigmaV"]), 0) + expect_gt(unname(cf[[g]]["sigmaU"]), 0) + } + + # Backward compatibility: default call unchanged + expect_identical(coef(fit_panel, which = "group"), fit_panel$group_coef) +}) + + +test_that("coef(..., extraPar = TRUE) appends stage-2 variances (meta)", { + cf_plain <- coef(fit_stoch) + cf_extra <- coef(fit_stoch, extraPar = TRUE) + + expect_gt(length(cf_extra), length(cf_plain)) + expect_true(all(c("sigmaV", "sigmaU") %in% names(cf_extra))) + expect_equal(cf_extra[seq_along(cf_plain)], cf_plain) + + # Deterministic metafrontier has no auxiliary parameters + expect_message(cf_det <- coef(fit_det, extraPar = TRUE), + "No auxiliary parameters") + expect_identical(cf_det, coef(fit_det)) +}) + + +test_that("vcov(..., extraPar = TRUE) returns the full stage-2 matrix", { + v <- vcov(fit_stoch) + v_full <- vcov(fit_stoch, extraPar = TRUE) + + expect_gt(nrow(v_full), nrow(v)) + expect_equal(unname(v_full[seq_len(nrow(v)), seq_len(ncol(v))]), + unname(v)) +}) + + +test_that("vcov(..., which = 'group') returns per-group matrices", { + v_g <- vcov(fit_stoch, which = "group") + + expect_true(is.list(v_g)) + expect_named(v_g, fit_stoch$groups) + for (g in fit_stoch$groups) { + gm <- fit_stoch$group_models[[g]] + expect_true(is.matrix(v_g[[g]])) + expect_equal(nrow(v_g[[g]]), length(gm$all_params)) + expect_true("log_sigma_v" %in% rownames(v_g[[g]])) + } +}) + + +test_that("summary of a BC92 fit shows eta with a standard error", { + s <- summary(fit_panel) + + for (g in fit_panel$groups) { + tab <- s$group_tables[[g]] + expect_true("eta" %in% rownames(tab)) + expect_true(is.finite(tab["eta", "Std. Error"])) + } + + expect_output(print(s), "eta") +}) + + +test_that("print methods report convergence status", { + expect_output(print(fit_det), "Convergence: OK") + expect_output(print(summary(fit_det)), + "All estimation stages converged") + expect_output(print(check_convergence(fit_det)), + "Convergence of estimation stages") + + # A non-converged stage triggers warnings in all three displays + fit_bad <- fit_det + fit_bad$group_models[[1]]$convergence <- 1L + expect_output(print(fit_bad), "Convergence: WARNING") + expect_output(print(summary(fit_bad)), "See \\?check_convergence") + expect_output(print(check_convergence(fit_bad)), + "did not converge") +}) + + +test_that("print shows Estimator and Objective lines when present", { + fit2 <- fit_det + if (is.null(fit2$estimator)) fit2$estimator <- "bc88" + if (is.null(fit2$objective)) fit2$objective <- "lp" + expect_output(print(fit2), "Estimator:\\s+bc88") + expect_output(print(fit2), "Objective:\\s+lp") +}) diff --git a/tests/testthat/test-frontier-equivalence.R b/tests/testthat/test-frontier-equivalence.R new file mode 100644 index 0000000..7dfdc1d --- /dev/null +++ b/tests/testthat/test-frontier-equivalence.R @@ -0,0 +1,26 @@ +# Equivalence check against frontier::sfa() +# The package's custom MLE should match frontier::sfa() slope and intercept +# estimates to within 1e-4 on the same simulated data. + +skip_if_not_installed("frontier") + +test_that("group SFA coefficients match frontier::sfa() to 1e-4", { + sim <- simulate_metafrontier(n_groups = 3, n_per_group = 150, + tech_gap = c(0, 0.3, 0.5), seed = 42) + dat <- sim$data + + fit_mf <- metafrontier(log_y ~ log_x1 + log_x2, data = dat, + group = "group", method = "sfa", + meta_type = "deterministic") + + for (g in c("G2", "G3")) { + dat_g <- dat[dat$group == g, , drop = FALSE] + f_g <- frontier::sfa(log_y ~ log_x1 + log_x2, data = dat_g) + + b_mf <- coef(fit_mf$group_models[[g]]) + b_fr <- coef(f_g)[names(b_mf)] + + expect_lt(max(abs(b_mf - b_fr)), 1e-4, + label = paste("max |coef diff| for group", g)) + } +}) diff --git a/tests/testthat/test-malmquist-id.R b/tests/testthat/test-malmquist-id.R new file mode 100644 index 0000000..8b76d49 --- /dev/null +++ b/tests/testthat/test-malmquist-id.R @@ -0,0 +1,188 @@ +# Tests for firm matching via `id`, infeasibility accounting, and +# method visibility in malmquist_meta(). + +# Deterministic positive-level panel: 2 groups x n_per_group firms x +# n_periods periods, suitable for DEA on levels. +make_id_panel <- function(n_per_group = 6, n_periods = 3) { + df <- expand.grid(firm = seq_len(2 * n_per_group), + time = seq_len(n_periods)) + df$group <- ifelse(df$firm <= n_per_group, "A", "B") + df$x1 <- 1 + 0.1 * (df$firm %% 5) + 0.05 * df$time + df$x2 <- 1.5 + 0.08 * (df$firm %% 7) + eff <- 0.6 + 0.4 * ((df$firm * 7) %% 10) / 10 + df$y <- eff * df$x1^0.4 * df$x2^0.3 * (1 + 0.1 * df$time) + df +} + +test_that("id gives order-invariant results on a balanced panel", { + pd <- make_id_panel() + + base <- suppressMessages( + malmquist_meta(y ~ x1 + x2, data = pd, + group = "group", time = "time") + ) + set.seed(1) + pd_scr <- pd[sample(nrow(pd)), ] + scr <- malmquist_meta(y ~ x1 + x2, data = pd_scr, + group = "group", time = "time", id = "firm") + + # Positional ids in the baseline are within-group positions; map + # them back to firm numbers (group A = 1..6, group B = 7..12). + base_m <- base$malmquist + base_m$firm <- base_m$id + ifelse(base_m$group == "B", 6L, 0L) + scr_m <- scr$malmquist + + base_m <- base_m[order(base_m$period_from, base_m$group, base_m$firm), ] + scr_m <- scr_m[order(scr_m$period_from, scr_m$group, scr_m$id), ] + + expect_equal(scr_m$id, base_m$firm) + expect_equal(scr_m$MPI, base_m$MPI, tolerance = 1e-8) + expect_equal(scr_m$TEC, base_m$TEC, tolerance = 1e-8) + expect_equal(scr_m$TGC, base_m$TGC, tolerance = 1e-8) + expect_equal(scr_m$TC, base_m$TC, tolerance = 1e-8) + + # Method visibility + expect_identical(scr$method, "dea") + expect_output(print(scr), "Method") +}) + +test_that("unbalanced panels drop unmatched firms with a consolidated warning", { + pd <- make_id_panel() + + # Make firms 3 (group A) and 9 (group B) strictly dominated (same + # inputs as firms 2 and 8, much lower output) so that dropping them + # leaves every DEA reference technology unchanged. + for (tt in unique(pd$time)) { + r3 <- pd$firm == 3 & pd$time == tt + r2 <- pd$firm == 2 & pd$time == tt + pd[r3, c("x1", "x2")] <- pd[r2, c("x1", "x2")] + pd$y[r3] <- 0.4 * pd$y[r2] + r9 <- pd$firm == 9 & pd$time == tt + r8 <- pd$firm == 8 & pd$time == tt + pd[r9, c("x1", "x2")] <- pd[r8, c("x1", "x2")] + pd$y[r9] <- 0.4 * pd$y[r8] + } + + bal <- malmquist_meta(y ~ x1 + x2, data = pd, + group = "group", time = "time", id = "firm") + + pd_unb <- pd[!(pd$firm == 3 & pd$time == 2) & + !(pd$firm == 9 & pd$time == 1), ] + # 3 dropped in total: firm 3 + firm 9 in pair 1 -> 2, firm 3 in 2 -> 3 + expect_warning( + unb <- malmquist_meta(y ~ x1 + x2, data = pd_unb, + group = "group", time = "time", id = "firm"), + "3 observation" + ) + + # The right firms are dropped + expect_false(any(unb$malmquist$id == 3)) + expect_false(any(unb$malmquist$id == 9 & unb$malmquist$period_from == 1)) + expect_true(any(unb$malmquist$id == 9 & unb$malmquist$period_from == 2)) + + # Retained firms match the balanced results exactly + key_b <- paste(bal$malmquist$id, bal$malmquist$period_from) + key_u <- paste(unb$malmquist$id, unb$malmquist$period_from) + expect_true(all(key_u %in% key_b)) + pos <- match(key_u, key_b) + expect_equal(unb$malmquist$MPI, bal$malmquist$MPI[pos], tolerance = 1e-8) + expect_equal(unb$malmquist$TC, bal$malmquist$TC[pos], tolerance = 1e-8) +}) + +test_that("duplicated (id, period) combinations within a group error", { + pd <- make_id_panel() + pd_dup <- rbind(pd, pd[pd$firm == 1 & pd$time == 1, ]) + expect_error( + malmquist_meta(y ~ x1 + x2, data = pd_dup, + group = "group", time = "time", id = "firm"), + "Duplicated" + ) + expect_error( + malmquist_meta(y ~ x1 + x2, data = pd_dup, + group = "group", time = "time", id = "firm"), + "period 1" + ) +}) + +test_that("id = NULL messages, and unequal group sizes warn", { + pd <- make_id_panel() + + expect_message( + malmquist_meta(y ~ x1 + x2, data = pd, + group = "group", time = "time"), + "row position" + ) + + pd_unb <- pd[!(pd$firm == 3 & pd$time == 2), ] + expect_warning( + suppressMessages( + malmquist_meta(y ~ x1 + x2, data = pd_unb, + group = "group", time = "time") + ), + "positional matching" + ) +}) + +test_that("infeasible cross-period programs are counted and warned under vrs", { + pd <- make_id_panel(n_per_group = 5, n_periods = 2) + + # Anchor firms (one per group) define the CRS frontier in both + # periods, so all cross-period CRS programs stay feasible. + pd[pd$firm == 1, c("x1", "x2", "y")] <- list(1, 1, 2) + pd[pd$firm == 6, c("x1", "x2", "y")] <- list(1, 1, 2) + # Tiny firm whose period-1 inputs undercut every period-2 + # observation: its period-1 cross-period VRS programs (group and + # meta) are genuinely infeasible. + pd[pd$firm == 2 & pd$time == 1, c("x1", "x2", "y")] <- + list(0.01, 0.01, 0.001) + pd[pd$firm == 2 & pd$time == 2, c("x1", "x2", "y")] <- + list(0.05, 0.05, 0.002) + + expect_warning( + malm_vrs <- malmquist_meta(y ~ x1 + x2, data = pd, + group = "group", time = "time", + id = "firm", rts = "vrs"), + "cross-period DEA programs were infeasible" + ) + expect_equal(malm_vrs$n_infeasible, 2L) + expect_true(any(is.na(malm_vrs$malmquist$TC))) + expect_true(is.data.frame(malm_vrs$infeasible_by_period)) + expect_equal(sum(malm_vrs$infeasible_by_period$n_infeasible), + malm_vrs$n_infeasible) + expect_output(print(malm_vrs), "Infeasible") + expect_output(print(summary(malm_vrs)), "Infeasible") + + # CRS on the same data has no infeasible programs and no NAs + expect_no_warning( + malm_crs <- malmquist_meta(y ~ x1 + x2, data = pd, + group = "group", time = "time", + id = "firm", rts = "crs") + ) + expect_identical(malm_crs$n_infeasible, 0L) + expect_false(any(is.na(malm_crs$malmquist$TC))) + expect_false(any(is.na(malm_crs$malmquist$MPI))) +}) + +test_that("sfa path notes the approximation and carries id through", { + set.seed(7) + n <- 12 + df <- expand.grid(firm = seq_len(2 * n), time = 1:2) + df$group <- ifelse(df$firm <= n, "A", "B") + df$log_x1 <- rnorm(nrow(df)) + df$log_x2 <- rnorm(nrow(df)) + df$log_y <- 1 + 0.4 * df$log_x1 + 0.3 * df$log_x2 + + 0.2 * (df$group == "B") + + rnorm(nrow(df), sd = 0.1) - abs(rnorm(nrow(df), sd = 0.3)) + + suppressWarnings( + expect_message( + malm <- malmquist_meta(log_y ~ log_x1 + log_x2, data = df, + group = "group", time = "time", + id = "firm", method = "sfa"), + "pointwise maximum" + ) + ) + expect_identical(malm$method, "sfa") + expect_true(all(malm$malmquist$id %in% df$firm)) + expect_output(print(malm), "pointwise-maximum") +}) diff --git a/tests/testthat/test-objective.R b/tests/testthat/test-objective.R new file mode 100644 index 0000000..84f5762 --- /dev/null +++ b/tests/testthat/test-objective.R @@ -0,0 +1,95 @@ +# Tests for the objective = c("lp", "qp") choice in the +# deterministic SFA metafrontier (Battese, Rao and O'Donnell 2004 +# propose both minimum absolute and minimum squared deviations). + +sim_obj <- simulate_metafrontier(seed = 42) +dat_obj <- sim_obj$data + +fit_lp_obj <- metafrontier(log_y ~ log_x1 + log_x2, + data = dat_obj, + group = "group", + method = "sfa", + meta_type = "deterministic") +X_obj <- model.matrix(log_y ~ log_x1 + log_x2, data = dat_obj) + +test_that("QP metafrontier satisfies envelope and TGR bounds", { + skip_if_not_installed("quadprog") + + qp <- metafrontier:::.deterministic_metafrontier_lp( + X_obj, fit_lp_obj$group_frontier, fit_lp_obj$group_vec, + fit_lp_obj$groups, fit_lp_obj$group_coef, ncol(X_obj), + objective = "qp" + ) + + expect_identical(qp$meta_solver, "qp") + + meta_frontier_qp <- as.numeric(X_obj %*% qp$meta_coef) + expect_true(all(meta_frontier_qp >= fit_lp_obj$group_frontier - 1e-6)) + + tgr_qp <- exp(fit_lp_obj$group_frontier - meta_frontier_qp) + expect_true(all(tgr_qp <= 1 + 1e-6)) + + # User-facing argument (added by a concurrent change); skip the + # assertions if metafrontier() does not honour 'objective' yet. + fit_qp <- tryCatch( + metafrontier(log_y ~ log_x1 + log_x2, + data = dat_obj, + group = "group", + method = "sfa", + meta_type = "deterministic", + objective = "qp"), + error = function(e) NULL + ) + if (!is.null(fit_qp) && identical(fit_qp$objective, "qp")) { + expect_identical(fit_qp$meta_solver, "qp") + expect_true(all(fit_qp$meta_frontier >= fit_qp$group_frontier - 1e-6)) + expect_true(all(fit_qp$tgr <= 1 + 1e-6)) + } +}) + +test_that("LP and QP objectives give similar TGR", { + skip_if_not_installed("quadprog") + + qp <- metafrontier:::.deterministic_metafrontier_lp( + X_obj, fit_lp_obj$group_frontier, fit_lp_obj$group_vec, + fit_lp_obj$groups, fit_lp_obj$group_coef, ncol(X_obj), + objective = "qp" + ) + + # Mirror the TGR construction of the deterministic pipeline + tgr_qp <- pmin(exp(fit_lp_obj$group_frontier - + as.numeric(X_obj %*% qp$meta_coef)), 1.0) + + max_diff <- max(abs(fit_lp_obj$tgr - tgr_qp)) + cat("\nMax abs TGR difference (LP vs QP):", format(max_diff), "\n") + + expect_lt(max_diff, 0.05) +}) + +test_that("rank-deficient X falls back to the barrier QP with a message", { + set.seed(99) + n <- 50 + x1 <- runif(n, 1, 2) + X_bad <- cbind("(Intercept)" = 1, x1 = x1, x1_dup = x1) + gf_bad <- 1 + 0.4 * x1 + rnorm(n, sd = 0.05) + gv_bad <- factor(rep(c("a", "b"), length.out = n)) + gc_bad <- list(a = c(1, 0.2, 0.2), b = c(1.1, 0.25, 0.15)) + + res <- NULL + expect_message( + res <- metafrontier:::.deterministic_metafrontier_lp( + X_bad, gf_bad, gv_bad, levels(gv_bad), gc_bad, ncol(X_bad), + objective = "qp" + ), + regexp = "constrOptim" + ) + + expect_identical(res$meta_solver, "qp-barrier") + expect_true(all(is.finite(res$meta_coef))) + expect_true(all(X_bad %*% res$meta_coef >= gf_bad - 1e-4)) +}) + +test_that("default fit records the LP solver and objective", { + expect_identical(fit_lp_obj$meta_solver, "lp") + expect_identical(fit_lp_obj$objective, "lp") +}) diff --git a/tests/testthat/test-panel-fixes.R b/tests/testthat/test-panel-fixes.R new file mode 100644 index 0000000..4d53e51 --- /dev/null +++ b/tests/testthat/test-panel-fixes.R @@ -0,0 +1,136 @@ +# Regression tests for panel SFA fixes: row order, NA handling, +# global-T decay anchoring, and unbalanced BC92 estimation + +.panel_test_group <- function(seed = 42, eta = 0.1, n_firms = 30, + n_periods = 5) { + sim <- simulate_panel_metafrontier( + n_groups = 2, n_firms_per_group = n_firms, + n_periods = n_periods, eta = eta, seed = seed + ) + g1 <- sim$data[sim$data$group == "G1", ] + # Firm IDs F1...F30: string sort order (F1, F10, F11, ...) differs + # from numeric order, which is what exposed the released row-order bug + g1$firm <- sub("^G1_", "", g1$firm) + list(data = g1, params = sim$params) +} + +.panel_formula <- function() Formula::Formula(log_y ~ log_x1 + log_x2) +.panel_info <- function() list(id = "firm", time = "year") + +test_that("panel fit is invariant to row order", { + g1 <- .panel_test_group()$data + set.seed(7) + g1_shuffled <- g1[sample(nrow(g1)), ] + + fit_sorted <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1, "hnormal", "bc92", .panel_info(), list() + ) + fit_shuffled <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1_shuffled, "hnormal", "bc92", .panel_info(), list() + ) + + expect_equal(fit_sorted$coefficients, fit_shuffled$coefficients, + tolerance = 1e-8) + + # Efficiency must agree per (firm, time) key, i.e. returned vectors + # are in input-row order for both fits + key_sorted <- paste(g1$firm, g1$year) + key_shuffled <- paste(g1_shuffled$firm, g1_shuffled$year) + expect_equal(fit_shuffled$efficiency, + fit_sorted$efficiency[match(key_shuffled, key_sorted)], + tolerance = 1e-8) +}) + +test_that("NA rows are dropped consistently across y, X and firm index", { + g1 <- .panel_test_group()$data + g1_na <- g1 + g1_na$log_x1[5] <- NA + + fit_na <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1_na, "hnormal", "bc92", .panel_info(), list() + ) + fit_dropped <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1[-5, ], "hnormal", "bc92", .panel_info(), list() + ) + + expect_length(fit_na$efficiency, nrow(g1) - 1) + expect_identical(as.character(fit_na$firms), as.character(g1$firm[-5])) + expect_equal(fit_na$efficiency, fit_dropped$efficiency, tolerance = 1e-10) + expect_equal(fit_na$coefficients, fit_dropped$coefficients, + tolerance = 1e-10) +}) + +test_that("unbalanced BC92 converges and recovers eta", { + gen <- .panel_test_group(seed = 456, eta = 0.1, n_firms = 50, + n_periods = 6) + g1 <- gen$data + set.seed(3) + g1_unbal <- g1[-sample(nrow(g1), 30), ] + + fit <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1_unbal, "hnormal", "bc92", .panel_info(), list() + ) + + expect_identical(fit$convergence, 0L) + expect_true(all(fit$efficiency > 0 & fit$efficiency <= 1)) + # Loose tolerance: eta is weakly identified in short panels + expect_lt(abs(fit$eta - gen$params$eta), 0.15) +}) + +test_that("BC92 decay is anchored at the global final period", { + g1 <- .panel_test_group(seed = 42, eta = 0.15)$data + # Remove firm F1's final period so its own last period < global T + g1_unbal <- g1[!(g1$firm == "F1" & g1$year == max(g1$year)), ] + + fit <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1_unbal, "hnormal", "bc92", .panel_info(), list() + ) + + # Recompute BC92 efficiency by hand with d_t = exp(-eta * (t - T)), + # T the GLOBAL final period; must match the stored values exactly + T_max <- max(g1_unbal$year) + sv <- fit$sigma_v + su <- fit$sigma_u + eta_hat <- fit$eta + idx_f1 <- which(fit$firms == "F1") + eps_f1 <- fit$residuals[idx_f1] + d_t <- exp(-eta_hat * (g1_unbal$year[g1_unbal$firm == "F1"] - T_max)) + + sigma_star2 <- 1 / (1 / su^2 + sum(d_t^2) / sv^2) + sigma_star <- sqrt(sigma_star2) + mu_star <- -sigma_star2 * sum(eps_f1 * d_t) / sv^2 + ratio <- mu_star / sigma_star + + te_bc88_hand <- pnorm(ratio - d_t * sigma_star) / pnorm(ratio) * + exp(-d_t * mu_star + 0.5 * d_t^2 * sigma_star2) + + expect_equal(fit$efficiency_bc88[idx_f1], te_bc88_hand, tolerance = 1e-10) + + # With eta_hat > 0 inefficiency decays over time, so F1's efficiency + # path must be monotone in the direction implied by the sign of eta + eff_path <- fit$efficiency_bc88[idx_f1][order(g1_unbal$year[g1_unbal$firm == "F1"])] + if (fit$eta > 0) { + expect_true(all(diff(eff_path) > 0)) + } else if (fit$eta < 0) { + expect_true(all(diff(eff_path) < 0)) + } +}) + +test_that("panel estimator argument works like the cross-sectional one", { + g1 <- .panel_test_group()$data + + fit_def <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1, "hnormal", "bc92", .panel_info(), list() + ) + fit_jlms <- metafrontier:::.fit_sfa_panel_group( + .panel_formula(), g1, "hnormal", "bc92", .panel_info(), list(), + estimator = "jlms" + ) + + expect_identical(fit_def$estimator, "bc88") + expect_identical(fit_def$efficiency, fit_def$efficiency_bc88) + expect_identical(fit_jlms$estimator, "jlms") + expect_identical(fit_jlms$efficiency, fit_jlms$efficiency_jlms) + expect_true(all(fit_def$efficiency_bc88 > 0 & fit_def$efficiency_bc88 <= 1)) + expect_gt(cor(fit_def$efficiency_bc88, fit_def$efficiency_jlms), 0.99) +}) diff --git a/tests/testthat/test-plots-interop.R b/tests/testthat/test-plots-interop.R new file mode 100644 index 0000000..d0ed827 --- /dev/null +++ b/tests/testthat/test-plots-interop.R @@ -0,0 +1,111 @@ +# Tests for revised plots (Figures 4 and 5) and model interoperability + +# ---------- Plot revisions ---------- + +test_that("autoplot.boot_tgr distribution shows CI bounds as vlines", { + skip_if_not_installed("ggplot2") + + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 40, seed = 42) + fit <- metafrontier(log_y ~ log_x1 + log_x2, data = sim$data, + group = "group") + boot <- boot_tgr(fit, R = 20, seed = 1, progress = FALSE) + + p <- ggplot2::autoplot(boot, which = "distribution") + expect_s3_class(p, "gg") + + is_vline <- vapply(p$layers, function(l) inherits(l$geom, "GeomVline"), + logical(1)) + expect_true(any(is_vline)) + + # One dashed line per group and per CI bound + vline_data <- ggplot2::layer_data(p, which(is_vline)[1]) + expect_equal(nrow(vline_data), 2L * length(boot$groups)) +}) + +test_that("autoplot.malmquist_meta mpi_trend plots at end periods", { + skip_if_not_installed("ggplot2") + + psim <- simulate_panel_metafrontier(n_groups = 2, n_firms_per_group = 15, + n_periods = 3, seed = 7) + malm <- suppressWarnings(suppressMessages( + malmquist_meta(log_y ~ log_x1 + log_x2, data = psim$data, + group = "group", time = "year", method = "dea") + )) + + p <- ggplot2::autoplot(malm, which = "mpi_trend") + expect_s3_class(p, "gg") + + built <- ggplot2::ggplot_build(p) + xvals <- sort(unique(unlist(lapply(built$data[1:2], + function(d) d$x)))) + expect_equal(as.numeric(xvals), 2:3) + expect_match(p$labels$x, "change from previous period") +}) + +# ---------- as_metafrontier_model interoperability ---------- + +test_that("as_metafrontier_model is idempotent and feeds metafrontier()", { + mk_group <- function(dat) { + X <- model.matrix(~ log_x1 + log_x2, dat) + ols <- lm.fit(X, dat$log_y) + res <- ols$residuals + list(coefficients = ols$coefficients, + efficiency = exp(res - max(res)), + X = X, y = dat$log_y) + } + grps <- levels(test_data$group) + conv1 <- as_metafrontier_model(mk_group(test_data[test_data$group == + grps[1], ])) + conv2 <- as_metafrontier_model(mk_group(test_data[test_data$group == + grps[2], ])) + + expect_s3_class(conv1, "metafrontier_model") + # Converting twice is a no-op + expect_identical(as_metafrontier_model(conv1), conv1) + + fit <- metafrontier(models = list(A = conv1, B = conv2)) + expect_s3_class(fit, "metafrontier") + expect_equal(fit$groups, c("A", "B")) +}) + +test_that("frontier::sfa objects pass directly through metafrontier(models=)", { + skip_if_not_installed("frontier") + + grps <- levels(test_data$group) + fits <- lapply(grps, function(g) { + suppressWarnings( + frontier::sfa(log_y ~ log_x1 + log_x2, + data = test_data[test_data$group == g, ]) + ) + }) + names(fits) <- grps + + expect_s3_class(fits[[1]], "frontier") + + fit <- metafrontier(models = fits) + expect_s3_class(fit, "metafrontier") + expect_equal(unname(fit$nobs["total"]), nrow(test_data)) + expect_true(all(fit$tgr > 0 & fit$tgr <= 1)) +}) + +test_that("Farrell objects warn on conversion and error in metafrontier()", { + skip_if_not_installed("Benchmarking") + + X <- exp(as.matrix(test_data[, c("log_x1", "log_x2")])) + y <- matrix(exp(test_data$log_y), ncol = 1) + d <- Benchmarking::dea(X, y, RTS = "vrs") + + expect_warning(as_metafrontier_model(d), "cannot be used") + + expect_error( + suppressWarnings(metafrontier(models = list(A = d, B = d))), + "does not provide" + ) +}) + +test_that("unsupported classes get an informative error", { + expect_error( + as_metafrontier_model(structure(1, class = "no_such_model")), + "sfacross.*frontier.*Farrell" + ) +}) diff --git a/tests/testthat/test-poolability-dea.R b/tests/testthat/test-poolability-dea.R new file mode 100644 index 0000000..ef00aec --- /dev/null +++ b/tests/testthat/test-poolability-dea.R @@ -0,0 +1,110 @@ +# Test the DEA permutation poolability test and the data.name fix + +test_that("permutation test rejects poolability when groups differ", { + skip_on_cran() + + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 60, + tech_gap = c(0, 0.6), seed = 1) + fit_dea <- metafrontier(log_y ~ log_x1 + log_x2, + data = sim$data, group = "group", + method = "dea", rts = "vrs") + + pt <- poolability_test(fit_dea, B = 99, seed = 123) + + expect_lt(pt$p.value, 0.05) +}) + + +test_that("permutation test does not reject under a pooled technology", { + skip_on_cran() + + # Sanity check under H0 with a fixed seed, not a power study + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 60, + tech_gap = c(0, 0), seed = 2) + fit_dea <- metafrontier(log_y ~ log_x1 + log_x2, + data = sim$data, group = "group", + method = "dea", rts = "vrs") + + pt <- poolability_test(fit_dea, B = 99, seed = 123) + + expect_gt(pt$p.value, 0.05) +}) + + +test_that("permutation test returns a well-formed htest object", { + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 25, + tech_gap = c(0, 0.4), seed = 7) + fit_dea <- metafrontier(log_y ~ log_x1 + log_x2, + data = sim$data, group = "group", + method = "dea", rts = "vrs") + + pt <- poolability_test(fit_dea, B = 49, seed = 9) + + expect_s3_class(pt, "htest") + expect_identical( + pt$method, + "Permutation test for poolability of group frontiers (DEA)" + ) + expect_named(pt$parameter, "B") + expect_lte(pt$parameter[["B"]], 49) + expect_named(pt$statistic, "mean technology gap") + expect_equal(pt$statistic[["mean technology gap"]], + mean(1 - fit_dea$tgr, na.rm = TRUE)) + expect_true(pt$p.value > 0 && pt$p.value <= 1) +}) + + +test_that("data.name is the symbol passed for DEA objects", { + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 25, + tech_gap = c(0, 0.4), seed = 7) + fit_dea <- metafrontier(log_y ~ log_x1 + log_x2, + data = sim$data, group = "group", + method = "dea", rts = "vrs") + + pt <- poolability_test(fit_dea, B = 9, seed = 1) + + expect_identical(pt$data.name, "fit_dea") +}) + + +test_that("SFA objects still take the LR branch with fixed data.name", { + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 60, + tech_gap = c(0, 0.4), seed = 5) + fit_sfa <- metafrontier(log_y ~ log_x1 + log_x2, + data = sim$data, group = "group", + method = "sfa") + + pt <- poolability_test(fit_sfa) + + expect_s3_class(pt, "htest") + expect_identical( + pt$method, + "Likelihood Ratio Test for Poolability of Group Frontiers" + ) + expect_named(pt$statistic, "LR") + expect_identical(pt$data.name, "fit_sfa") +}) + + +test_that("permutation test is reproducible with a fixed seed", { + skip_on_cran() + + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 25, + tech_gap = c(0, 0.4), seed = 7) + fit_dea <- metafrontier(log_y ~ log_x1 + log_x2, + data = sim$data, group = "group", + method = "dea", rts = "vrs") + + pt1 <- poolability_test(fit_dea, B = 49, seed = 42) + pt2 <- poolability_test(fit_dea, B = 49, seed = 42) + + expect_identical(pt1$p.value, pt2$p.value) +}) + + +test_that("non-metafrontier objects give an informative error", { + expect_error( + poolability_test(list()), + "SFA- or DEA-based metafrontier" + ) +}) diff --git a/tests/testthat/test-simulate-dgp.R b/tests/testthat/test-simulate-dgp.R new file mode 100644 index 0000000..cd9b6ec --- /dev/null +++ b/tests/testthat/test-simulate-dgp.R @@ -0,0 +1,147 @@ +# DGP extensions: group-specific slopes, input distributions, attrition +# (v0.3.0). Legacy reference values were generated with +# the pre-0.3.0 implementation. + +test_that("simulate_metafrontier reproduces legacy output for seed 42", { + sim <- simulate_metafrontier(seed = 42) + + expect_equal(nrow(sim$data), 200L) + expect_equal(ncol(sim$data), 9L) + expect_equal(sim$data$log_x1[1], 4.574030217481777, tolerance = 1e-12) + expect_equal(sim$data$log_y[1], 3.553174757064474, tolerance = 1e-12) + expect_equal(sim$data$true_te[57], 0.801290099145637, tolerance = 1e-12) + expect_equal(sim$data$log_y[200], 1.470010450855711, tolerance = 1e-12) + expect_equal(sim$data$true_te_star[150], 0.401389206182210, + tolerance = 1e-12) +}) + +test_that("beta_groups with differing slopes gives per-observation TGR in (0, 1]", { + bg <- rbind(c(1.0, 0.5, 0.2), + c(0.9, 0.6, 0.1)) + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 200, + beta_groups = bg, seed = 1) + + expect_true(all(sim$data$true_tgr > 0 & sim$data$true_tgr <= 1)) + + # TGR varies within each group when slopes differ + sds <- tapply(sim$data$true_tgr, sim$data$group, stats::sd) + expect_true(all(sds > 0)) + + # No single log-linear metafrontier: beta_meta is NULL, + # beta_groups is returned + expect_null(sim$params$beta_meta) + expect_equal(sim$params$beta_groups$G1, c(1.0, 0.5, 0.2)) + expect_equal(sim$params$beta_groups$G2, c(0.9, 0.6, 0.1)) + + # Decomposition identity still holds + expect_equal(sim$data$true_te_star, + sim$data$true_te * sim$data$true_tgr, + tolerance = 1e-10) +}) + +test_that("beta_groups replicating the intercept-only design gives exp(-gap)", { + bg <- rbind(c(1.0, 0.5, 0.2), + c(0.75, 0.5, 0.2)) + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 100, + beta_groups = bg, seed = 7) + + tgr1 <- sim$data$true_tgr[sim$data$group == "G1"] + tgr2 <- sim$data$true_tgr[sim$data$group == "G2"] + expect_equal(tgr1, rep(1, 100), tolerance = 1e-10) + expect_equal(tgr2, rep(exp(-0.25), 100), tolerance = 1e-10) +}) + +test_that("beta_groups accepts a list and warns when tech_gap is also given", { + bg_list <- list(c(1.0, 0.5, 0.2), c(0.9, 0.6, 0.1)) + sim <- simulate_metafrontier(n_groups = 2, beta_groups = bg_list, seed = 2) + expect_equal(sim$params$beta_groups$G2, c(0.9, 0.6, 0.1)) + + expect_warning( + simulate_metafrontier(n_groups = 2, beta_groups = bg_list, + tech_gap = c(0, 0.3), seed = 2), + "tech_gap" + ) + expect_error( + simulate_metafrontier(n_groups = 2, + beta_groups = rbind(c(1, 0.5, 0.2)), + seed = 2), + "beta_groups" + ) +}) + +test_that("input_corr induces the target input correlation", { + R <- matrix(c(1, 0.6, + 0.6, 1), nrow = 2) + sim <- simulate_metafrontier(n_groups = 1, n_per_group = 2000, + input_corr = R, seed = 3) + + realised <- stats::cor(sim$data$log_x1, sim$data$log_x2) + expect_lt(abs(realised - 0.6), 0.15) + + # Spread roughly matches the legacy uniform draws + expect_lt(abs(stats::sd(sim$data$log_x1) - 5 / sqrt(12)), 0.15) +}) + +test_that("input_corr is validated", { + not_psd <- matrix(c(1, 2, + 2, 1), nrow = 2) + expect_error(simulate_metafrontier(input_corr = not_psd, seed = 1), + "positive definite") + not_sym <- matrix(c(1, 0.2, + 0.6, 1), nrow = 2) + expect_error(simulate_metafrontier(input_corr = not_sym, seed = 1), + "symmetric") +}) + +test_that("input_means centres per-group log-inputs", { + m <- rbind(c(1, 2), + c(3, 4)) + sim <- simulate_metafrontier(n_groups = 2, n_per_group = 500, + input_means = m, seed = 4) + + for (g in 1:2) { + d <- sim$data[sim$data$group == paste0("G", g), ] + expect_lt(abs(mean(d$log_x1) - m[g, 1]), 0.15) + expect_lt(abs(mean(d$log_x2) - m[g, 2]), 0.15) + } + + expect_error(simulate_metafrontier(n_groups = 2, + input_means = rbind(c(1, 2)), + seed = 4), + "input_means") +}) + +test_that("attrition = 0 reproduces the legacy balanced panel", { + sim <- simulate_panel_metafrontier(seed = 99) + + expect_equal(nrow(sim$data), 300L) + expect_equal(sim$data$log_y[1], 2.103142882313839, tolerance = 1e-12) + expect_equal(sim$data$log_y[300], 0.928893291616277, tolerance = 1e-12) + expect_equal(sim$data$true_tgr[250], 0.606530659712633, tolerance = 1e-12) + + sim0 <- simulate_panel_metafrontier(seed = 99, attrition = 0) + expect_identical(sim$data, sim0$data) + expect_equal(sim0$params$attrition_share, 0) +}) + +test_that("attrition yields an unbalanced panel with all firms in period 1", { + sim <- simulate_panel_metafrontier(seed = 123, attrition = 0.3) + d <- sim$data + n_firms <- 2 * 30 + + expect_lt(nrow(d), 300) + expect_equal(length(unique(d$firm)), n_firms) + expect_equal(sum(d$year == 1), n_firms) + + # Unbalanced: firms have differing numbers of periods + periods_per_firm <- table(d$firm) + expect_gt(length(unique(as.integer(periods_per_firm))), 1L) + + # Realised share stored and consistent with the rows dropped + expect_equal(sim$params$attrition_share, + (300 - nrow(d)) / (300 - n_firms)) + expect_gt(sim$params$attrition_share, 0) + + expect_error(simulate_panel_metafrontier(seed = 1, attrition = 0.9), + "attrition") +}) diff --git a/vignettes/introduction.Rmd b/vignettes/introduction.Rmd index ea1b0e4..71a8baf 100644 --- a/vignettes/introduction.Rmd +++ b/vignettes/introduction.Rmd @@ -19,16 +19,26 @@ knitr::opts_chunk$set( ## What is a metafrontier? In efficiency analysis, we often study firms that operate under -fundamentally different technologies. Steel producers using electric arc -furnaces (EAF) face a different production possibility set than those using -the blast furnace-basic oxygen furnace (BF-BOF) route. Hospitals in rural -areas face different constraints than urban ones. Banks in developing -economies operate under different regulatory and technological environments -than those in advanced economies. +different production environments. Steel producers using electric arc +furnaces (EAF) face a different feasible set of input-output combinations +than those using the blast furnace-basic oxygen furnace (BF-BOF) route. +Hospitals in rural areas face different constraints than urban ones. +Banks face different regulatory environments across jurisdictions. + +Following Battese, Rao, and O'Donnell (2004) and O'Donnell, Rao, and +Battese (2008), we conceive of a single industry **metatechnology** +$T^*$: the set of all input-output combinations that are technically +feasible in the industry. Each group of firms operates within a +**restricted subset** $T_j \subseteq T^*$ of this metatechnology, where +the restrictions arise from regulation, the physical environment, +resource endowments, or the cost of switching production systems. Groups +do not possess fundamentally different technologies; they face different +restrictions of a common metatechnology. Standard stochastic frontier analysis (SFA) or data envelopment analysis -(DEA) applied to the pooled sample implicitly assumes all firms share the -same technology -- an assumption that may be unrealistic. Estimating +(DEA) applied to the pooled sample implicitly assumes that all firms +have unrestricted access to the same technology set, an assumption +that may be unrealistic. Estimating separate frontiers for each group solves this problem but makes efficiency scores incomparable across groups: a firm that is 90\% efficient relative to a less advanced group frontier may actually be less productive than a @@ -111,10 +121,15 @@ The deterministic metafrontier is estimated in two stages: 1. **Stage 1**: Fit separate SFA models for each group via maximum likelihood. 2. **Stage 2**: Find metafrontier coefficients $\hat\beta^*$ by minimising - $$\sum_i \left[\ln f(x_i; \hat\beta^*) - \ln f(x_i; \hat\beta_j)\right]^2$$ + $$\sum_i \left|\ln f(x_i; \hat\beta^*) - \ln f(x_i; \hat\beta_j)\right|$$ subject to the constraint that the metafrontier envelops all group frontiers: $\ln f(x_i; \hat\beta^*) \ge \ln f(x_i; \hat\beta_j)$ for - all $i$ and $j$. + all $i$ and $j$. Since the envelopment constraints force the + deviations to be non-negative, this reduces to a linear programme + (`objective = "lp"`, the default). The alternative minimum sum of + squared deviations criterion is available via `objective = "qp"`; + both criteria are proposed by Battese, Rao, and O'Donnell (2004), + and the `methods` vignette discusses them in detail. This is the default method: @@ -198,6 +213,16 @@ te_star <- efficiencies(fit_det, type = "meta") all.equal(te_star, te * tgr) ``` +For SFA fits, technical efficiencies are computed with the Battese and +Coelli (1988) conditional expectation estimator by default +(`estimator = "bc88"`). The Jondrow et al. (1982) estimator is computed +and stored alongside it, so you can switch without refitting: + +```{r efficiencies-jlms} +te_jlms <- efficiencies(fit_det, type = "group", estimator = "jlms") +cor(te, te_jlms) +``` + ### Technology gap ratio The `technology_gap_ratio()` function returns TGR values grouped by @@ -277,9 +302,24 @@ statistically different from a single pooled frontier: poolability_test(fit_det) ``` -A significant result (small p-value) indicates that the technology groups -have genuinely different production technologies, justifying the -metafrontier approach. +A significant result (small p-value) indicates that the group frontiers +genuinely differ, that is, the groups face different restrictions of the +common metatechnology, justifying the metafrontier approach. + +### Convergence diagnostics + +Every estimation stage of a metafrontier fit can be inspected with +`check_convergence()`, which returns one row per stage (each group +frontier and the metafrontier itself) with the estimation method, the +optimiser's convergence code, and a logical convergence indicator: + +```{r convergence} +check_convergence(fit_det) +``` + +The `summary()` method also prints a convergence block, so estimation +problems are flagged even if you never call `check_convergence()` +directly. ## Inefficiency distributions @@ -455,6 +495,10 @@ head(data.frame( ## References +- Battese, G.E. and Coelli, T.J. (1988). Prediction of firm-level + technical efficiencies with a generalized frontier production function + and panel data. *Journal of Econometrics*, 38(3), 387--399. + - Battese, G.E., Rao, D.S.P. and O'Donnell, C.J. (2004). A metafrontier production function for estimation of technical efficiencies and technology gaps for firms operating under different technologies. @@ -464,6 +508,10 @@ head(data.frame( estimating the metafrontier production function based on a stochastic frontier framework. *Journal of Productivity Analysis*, 42(3), 241--254. +- Jondrow, J., Lovell, C.A.K., Materov, I.S. and Schmidt, P. (1982). On + the estimation of technical inefficiency in the stochastic frontier + production function model. *Journal of Econometrics*, 19(2--3), 233--238. + - O'Donnell, C.J., Rao, D.S.P. and Battese, G.E. (2008). Metafrontier frameworks for the study of firm-level efficiencies and technology ratios. *Empirical Economics*, 34(2), 231--255. diff --git a/vignettes/malmquist.Rmd b/vignettes/malmquist.Rmd index 4bd6306..f25678a 100644 --- a/vignettes/malmquist.Rmd +++ b/vignettes/malmquist.Rmd @@ -24,9 +24,10 @@ library(metafrontier) Standard Malmquist productivity indices measure productivity change over time by decomposing it into efficiency change and technical change. When -firms operate under different technologies, however, this decomposition -misses an important dimension: changes in the *technology gap* between a -group's frontier and the global best practice. +groups of firms face different restrictions of a common metatechnology, +however, this decomposition misses an important dimension: changes in +the *technology gap* between a group's frontier and the global best +practice. The **metafrontier Malmquist TFP index** of O'Donnell, Rao, and Battese (2008) extends the standard index with a three-way decomposition that @@ -97,6 +98,7 @@ malm <- malmquist_meta( data = panel_data, group = "group", time = "time", + id = "id", orientation = "output", rts = "crs" ) @@ -104,6 +106,10 @@ malm <- malmquist_meta( malm ``` +The `id` argument identifies firms across periods. Firms not observed +in both periods of a consecutive pair are dropped from that comparison, +and `malmquist_meta()` warns about how many observations were dropped. + ## Detailed results @@ -190,6 +196,7 @@ malm_vrs <- malmquist_meta( data = panel_data, group = "group", time = "time", + id = "id", rts = "vrs" ) @@ -218,13 +225,16 @@ malm_us <- malmquist_meta( data = Produc, group = "region", time = "year", + id = "state", rts = "crs" ) summary(malm_us) ``` Similarly, `sfaR::utility` provides electric utility data with a binary -`regu` variable (regulated vs. deregulated) over 1986--1996: +`regu` variable (retail-access vs. traditionally regulated states) over +1986--1996, where `regu = 1` indicates that the state had enacted +retail-access (restructuring) legislation: ```{r utility-example, eval = FALSE} library(sfaR) @@ -235,14 +245,59 @@ malm_util <- malmquist_meta( data = utility, group = "regu", time = "year", + id = "firm", rts = "vrs" ) summary(malm_util) ``` +## Caveats + +A few limitations should be kept in mind when interpreting the results. + +### Interpretation of the index + +The Malmquist index, and by extension its metafrontier version, is not +a proper index in the sense of O'Donnell (2012): it is not +multiplicatively complete, and its binary comparisons are not +transitive, so chained comparisons across several periods depend on the +path taken. Interpretation should therefore be confined to binary +period-to-period comparisons; avoid cumulating the indices into a +multi-period productivity level series. + +### Cross-period infeasibility + +Under `rts = "vrs"`, `"drs"`, `"irs"`, or `"fdh"`, the cross-period +distance functions (evaluating a period-$t$ observation against the +period-$s$ frontier) may have no feasible solution. Infeasible cases +are returned as `NA` (never `Inf`), counted, and reported, so summary +statistics over the remaining observations should be read with the +reduced sample in mind. Under `rts = "crs"` the programmes are always +feasible. + +### Firm matching + +Observations are matched across consecutive periods via the `id` +variable. Firms observed in only one period of a pair are dropped from +that comparison, and `malmquist_meta()` warns about the number of +dropped observations. With substantial attrition or entry, the index is +computed on the surviving matched sample only. + +### SFA-based index + +With `method = "sfa"`, the metafrontier at each period is formed as the +pointwise maximum of the fitted group frontiers, an approximation that +the function announces with a message. For exact decompositions, prefer +the default `method = "dea"`. + + ## References +- O'Donnell, C.J. (2012). An aggregate quantity framework for measuring + and decomposing productivity change. *Journal of Productivity + Analysis*, 38(3), 255--272. + - O'Donnell, C.J., Rao, D.S.P. and Battese, G.E. (2008). Metafrontier frameworks for the study of firm-level efficiencies and technology ratios. *Empirical Economics*, 34(2), 231--255. diff --git a/vignettes/methods.Rmd b/vignettes/methods.Rmd index 6d716ba..2f45be4 100644 --- a/vignettes/methods.Rmd +++ b/vignettes/methods.Rmd @@ -43,7 +43,12 @@ The group-specific technical efficiency is: $$TE_{ij} = \exp(-u_{ij}) \in (0, 1]$$ -estimated via the Jondrow et al. (1982) conditional mean estimator. +estimated by default via the Battese and Coelli (1988) conditional +expectation estimator $E[\exp(-u_{ij}) \mid \varepsilon_{ij}]$ +(`estimator = "bc88"`). The Jondrow et al. (1982) estimator +$\exp(-E[u_{ij} \mid \varepsilon_{ij}])$ is also computed and stored, and +`efficiencies(fit, estimator = "jlms")` switches between the two without +refitting. ### 1.2 The metafrontier @@ -53,9 +58,14 @@ that: $$x'\beta^* \ge x'\beta_j \quad \text{for all } x \text{ and all } j$$ -That is, the metafrontier weakly dominates all group frontiers. It -represents the production technology available to firms with unrestricted -access to all technologies. +That is, the metafrontier weakly dominates all group frontiers. It is +the boundary of the industry **metatechnology** $T^*$: the set of all +input-output combinations that are technically feasible in the industry +(Battese, Rao and O'Donnell, 2004; O'Donnell, Rao and Battese, 2008). +Each group frontier bounds a restricted subset $T_j \subseteq T^*$ of +this metatechnology, where the restrictions reflect regulation, the +physical environment, resource endowments, or the cost of switching +production systems. ### 1.3 The efficiency decomposition @@ -76,15 +86,36 @@ that input mix; values below 1 indicate a technology gap. ### 2.1 Estimation -After obtaining group estimates $\hat\beta_j$ in Stage 1, the metafrontier -parameters $\hat\beta^*$ are estimated by solving: +After obtaining group estimates $\hat\beta_j$ in Stage 1, the +metafrontier parameters $\hat\beta^*$ are chosen so that the metafrontier +lies as close as possible to the group frontiers while enveloping all of +them. Battese, Rao and O'Donnell (2004) propose two optimisation +criteria, both implemented in the package via the `objective` argument. -$$\min_{\beta^*} \sum_{j=1}^{J} \sum_{i=1}^{n_j} \left(x_{ij}'\beta^* - x_{ij}'\hat\beta_j\right)^2$$ +**Minimum sum of absolute deviations** (`objective = "lp"`, the default): + +$$\min_{\beta^*} \sum_{j=1}^{J} \sum_{i=1}^{n_j} \left|x_{ij}'\beta^* - x_{ij}'\hat\beta_j\right|$$ $$\text{subject to: } x_{ij}'\beta^* \ge x_{ij}'\hat\beta_j \quad \forall\, i, j$$ -This is a convex quadratic program. The `metafrontier` package solves it -using `constrOptim()` from base R, which implements an adaptive barrier -algorithm for linearly constrained optimisation. +Because the envelopment constraints force every deviation to be +non-negative, the absolute values can be dropped and the problem reduces +to a linear programme (O'Donnell, Rao and Battese, 2008, Eqs. 23-25). +This LP is the primary presentation in O'Donnell, Rao and Battese (2008) +and is the package default; it is solved with `lpSolveAPI`. + +**Minimum sum of squared deviations** (`objective = "qp"`): + +$$\min_{\beta^*} \sum_{j=1}^{J} \sum_{i=1}^{n_j} \left(x_{ij}'\beta^* - x_{ij}'\hat\beta_j\right)^2$$ + +subject to the same envelopment constraints. This is a convex quadratic +programme, solved exactly with `quadprog` when that package is +available, with an adaptive-barrier fallback via `constrOptim()` from +base R otherwise. + +The two criteria weight deviations differently (the squared criterion +penalises large gaps between the metafrontier and a group frontier more +heavily), but in practice they typically produce very similar +coefficients. ### 2.2 Properties @@ -129,6 +160,25 @@ group_b0 <- sapply(coef(fit_det, which = "group"), `[`, 1) meta_b0 >= group_b0 ``` +To use the minimum sum of squared deviations criterion instead, set +`objective = "qp"`: + +```{r det-qp} +fit_qp <- metafrontier( + log_y ~ log_x1 + log_x2, + data = sim$data, + group = "group", + meta_type = "deterministic", + objective = "qp" +) + +# LP and QP criteria typically give very similar coefficients +cbind( + LP = coef(fit_det, which = "meta"), + QP = coef(fit_qp, which = "meta") +) +``` + ## 3. Stochastic metafrontier (Huang, Huang, and Liu, 2014) @@ -349,8 +399,9 @@ on the pooled sample and $LL_j$ are the group-specific log-likelihoods. poolability_test(fit_det) ``` -A significant test (p < 0.05) confirms that the groups operate under -different technologies and the metafrontier decomposition is warranted. +A significant test (p < 0.05) confirms that the group frontiers differ, +that is, the groups face different restrictions of the common +metatechnology, and the metafrontier decomposition is warranted. ## 8. Simulation for Monte Carlo studies @@ -397,6 +448,10 @@ The `simulate_metafrontier()` function supports: ## References +- Battese, G.E. and Coelli, T.J. (1988). Prediction of firm-level + technical efficiencies with a generalized frontier production function + and panel data. *Journal of Econometrics*, 38(3), 387--399. + - Battese, G.E., Rao, D.S.P. and O'Donnell, C.J. (2004). A metafrontier production function for estimation of technical efficiencies and technology gaps for firms operating under different technologies.