metafrontier v0.3.0

.Rbuildignore

@@ -13,9 +13,16 @@
 ^QA_ROUND2_EDITORIAL_REPORT\.md$
 ^ROADMAP\.md$
 ^mc_review_r10\.R$
+^mc_murphy_topel\.R$
+^review_malmquist_r14\.R$
 ^reviewer2_dso_test\.R$
+^reviewer4_dea_test\.R$
+^reviewer8_test\.R$
 ^test_edge_cases\.R$
 ^test_panel_review\.R$
 ^test_reviewer16_edgecases\.R$
 ^test_reviewer9\.R$
 ^test_steel_review\.R$
+^test_vignette\.R$
+^submit_cran2\.R$
+^metafrontier_0\.2\.1\.tar\.gz$

DESCRIPTION

@@ -1,25 +1,32 @@
 Package: metafrontier
 Type: Package
 Title: Analysis of Metafrontier Models for Efficiency and Productivity
-Version: 0.2.2
+Version: 0.3.0
 Authors@R: c(
     person("Erik", "Enstad", email = "erik.enstad@nhh.no",
            role = c("aut", "cre"),
            comment = c(ORCID = "0009-0006-1053-3849")))
 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
+    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>. 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.
+    <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.
 License: GPL (>= 3)
 Encoding: UTF-8
 RoxygenNote: 7.3.3
@@ -39,6 +46,7 @@ Suggests:
     sfaR,
     frontier,
     Benchmarking,
+    quadprog,
     plm,
     ggplot2,
     parallel

NAMESPACE

@@ -2,9 +2,14 @@
 
 S3method(as_metafrontier_model,Farrell)
 S3method(as_metafrontier_model,default)
+S3method(as_metafrontier_model,frontier)
 S3method(as_metafrontier_model,list)
-S3method(as_metafrontier_model,sfa)
+S3method(as_metafrontier_model,metafrontier_model)
 S3method(as_metafrontier_model,sfacross)
+S3method(check_convergence,default)
+S3method(check_convergence,lc_metafrontier)
+S3method(check_convergence,malmquist_meta)
+S3method(check_convergence,metafrontier)
 S3method(coef,lc_metafrontier)
 S3method(coef,metafrontier)
 S3method(confint,boot_tgr)
@@ -24,6 +29,7 @@ S3method(print,boot_tgr)
 S3method(print,lc_metafrontier)
 S3method(print,malmquist_meta)
 S3method(print,metafrontier)
+S3method(print,metafrontier_convergence)
 S3method(print,summary.lc_metafrontier)
 S3method(print,summary.malmquist_meta)
 S3method(print,summary.metafrontier)
@@ -34,6 +40,7 @@ S3method(summary,metafrontier)
 S3method(vcov,metafrontier)
 export(as_metafrontier_model)
 export(boot_tgr)
+export(check_convergence)
 export(efficiencies)
 export(latent_class_metafrontier)
 export(malmquist_meta)

NEWS.md

@@ -1,3 +1,109 @@
+# metafrontier 0.3.0
+
+## Breaking changes
+
+- The default technical efficiency estimator for SFA models is now the
+  Battese-Coelli (1988) conditional expectation `E[exp(-u)|eps]`
+  (`estimator = "bc88"`), which provides consistent efficiency
+  estimates. The previous JLMS estimator remains available via
+  `estimator = "jlms"`, and both are stored on every fit;
+  `efficiencies(fit, estimator = )` switches without refitting.
+  Deterministic and stochastic TGRs are unaffected (they do not depend
+  on the estimator), but TE and TE* values change slightly relative to
+  0.2.x.
+- `malmquist_meta()` output `id` column now carries the user-supplied
+  firm identifier (previously a within-group loop index).
+- The `as_metafrontier_model()` method for `frontier::sfa()` fits is
+  now registered for the correct class `"frontier"`; the previous
+  registration (`"sfa"`) never dispatched.
+
+## Bug fixes
+
+- Panel SFA (BC92/BC95): fixed a row-alignment bug where
+  observation-level efficiencies were returned in internally sorted
+  (id, time) order but assigned back in input-row order, scrambling
+  `te_group`/`te_meta` whenever the input rows were not already in
+  string-sorted order. Coefficients, TGRs, and group means were
+  unaffected. Results are now row-order invariant (regression-tested).
+- BC92 decay is now anchored at the panel-wide final period `T`
+  (previously the firm-specific last period), matching Battese and
+  Coelli's (1992) unbalanced formulation and the package's own
+  simulator. Balanced panels are numerically unchanged.
+- Panel SFA now aligns the firm/time index with rows dropped by
+  `na.action`, and uses the same BFGS-to-Nelder-Mead fallback cascade
+  as the cross-sectional path.
+- DEA: removed an erroneous lower bound (`phi >= 1`) in the
+  output-oriented LPs that made cross-period evaluation of
+  super-efficient DMUs infeasible even under CRS, silently biasing
+  Malmquist TC and MPI means. Cross-period scores with `phi < 1` now
+  solve correctly.
+- `as_metafrontier_model()` is now idempotent (converting an already
+  converted object is a no-op), so the previously documented
+  pre-conversion workflow works.
+- `poolability_test()` now derives `data.name` from the passed
+  expression instead of deparsing the stored call.
+- `autoplot(boot)` (and the base `plot()` method) now actually draw
+  the dashed CI bound lines promised by the documentation.
+- `autoplot(malm, which = "mpi_trend")` now plots each transition at
+  its end period (axis 2..T), consistent with the caption convention
+  "change relative to the previous period".
+- Removed a redundant marginal likelihood computation in the latent
+  class EM loop.
+
+## New features
+
+- `estimator = c("bc88", "jlms")` on `metafrontier()` and
+  `malmquist_meta()` (see Breaking changes).
+- `objective = c("lp", "qp")` on `metafrontier()`: both identification
+  criteria of Battese, Rao and O'Donnell (2004) for the deterministic
+  metafrontier. The default LP minimises the sum of absolute
+  deviations (O'Donnell, Rao and Battese, 2008, Eqs. 23-25); the QP
+  minimises the sum of squared deviations, solved exactly via
+  `quadprog` (new in Suggests) with a `constrOptim()` barrier
+  fallback. The bootstrap respects the choice.
+- `engine = c("internal", "sfaR", "frontier", "Benchmarking")` on
+  `metafrontier()`: delegate group-frontier estimation to external
+  packages; `engine = "Benchmarking"` also delegates the pooled DEA
+  metafrontier via `XREF`/`YREF`.
+- `check_convergence()`: new exported diagnostic reporting one row per
+  estimation stage (group frontiers and metafrontier). `print()` and
+  `summary()` methods now include convergence status; every stage
+  warns on non-zero optimiser codes.
+- `malmquist_meta(id = )`: explicit firm matching across periods, with
+  errors on duplicated (id, period) pairs, warnings counting dropped
+  observations on unbalanced panels, and a message when falling back
+  to positional matching. Cross-period infeasible programmes (possible
+  under vrs/drs/irs/fdh) are now counted and reported in a
+  consolidated warning, stored as `n_infeasible`, and shown by
+  `print()`/`summary()`; the SFA path announces its pointwise-maximum
+  approximation in a message and in the documentation.
+- Unbalanced panels are fully supported in BC92 (firm-specific period
+  sets), and `simulate_panel_metafrontier()` gains an `attrition`
+  argument for generating unbalanced test panels.
+- DEA: `rts = "fdh"` (free disposable hull; exact enumeration for
+  radial measures, binary MIP for DDF), `type = "hyperbolic"` (graph
+  efficiency; closed form under CRS, bisection otherwise; always
+  feasible cross-period), user-supplied direction vectors or
+  firm-specific direction matrices for DDF (reported via the additive
+  `ddf_gap`), and `slack = TRUE` two-stage slack maximisation.
+- `poolability_test()` now dispatches a permutation test for DEA fits
+  (group labels exchangeable under the pooled-technology null), with
+  `B` and `seed` arguments.
+- `coef()`, `vcov()`, and `summary()` expose all estimated parameters:
+  `extraPar = TRUE` returns variance parameters (and `eta` for BC92)
+  with back-transformed values; `vcov(which = "group")` returns full
+  per-group covariance matrices; group summary tables now include the
+  variance parameters and `eta` with standard errors.
+- `simulate_metafrontier()` gains `beta_groups` (group-specific slope
+  coefficients; the true TGR is then computed against the pointwise
+  maximum over group frontiers and varies within groups),
+  `input_means` (group-specific input distributions), and
+  `input_corr` (correlated log inputs).
+- Latent class: the EM now records whether the best start met the
+  convergence tolerance (`em_converged`) and warns when it did not.
+
+---
+
 # metafrontier 0.2.2
 
 ## Bug fixes

R/bootstrap.R

@@ -60,6 +60,10 @@ boot_tgr <- function(object, R = 999,
     stop("'object' must be a fitted metafrontier model.", call. = FALSE)
   }
 
+  # The parametric bootstrap draws new noise and inefficiency terms
+  # from the estimated error distributions, which requires a
+  # distributional model. DEA is nonparametric and provides no such
+  # model, so only case resampling is valid for DEA fits.
   if (type == "parametric" && object$method == "dea") {
     stop("Parametric bootstrap is not available for DEA metafrontiers. ",
          "Use type = 'nonparametric'.", call. = FALSE)
@@ -190,10 +194,13 @@ boot_tgr <- function(object, R = 999,
   groups <- object$groups
 
   if (type == "parametric") {
-    # Parametric: resample residuals from estimated distributions
+    # Parametric: keep the design fixed and regenerate the response by
+    # drawing new noise (v) and inefficiency (u) terms from the fitted
+    # group-specific error distributions
     boot_data <- .parametric_resample(object)
   } else {
-    # Nonparametric: resample rows within groups
+    # Nonparametric: case resampling, i.e. resample whole rows with
+    # replacement within each group so group sizes are preserved
     boot_data <- .nonparametric_resample(data, group_vec, groups)
   }
 
@@ -218,6 +225,19 @@ boot_tgr <- function(object, R = 999,
               !is.null(object$rts)) {
       object$rts
     } else "crs",
+    # Forward the estimation choices of the original fit so the
+    # bootstrap distribution reflects the reported point estimates.
+    type = if (!is.null(object$type)) object$type else "radial",
+    direction = if (!is.null(object$direction)) {
+      object$direction
+    } else "proportional",
+    estimator = if (!is.null(object$estimator)) {
+      object$estimator
+    } else "bc88",
+    objective = if (!is.null(object$objective)) {
+      object$objective
+    } else "lp",
+    engine = if (!is.null(object$engine)) object$engine else "internal",
     ...
   )
 
@@ -433,6 +453,12 @@ plot.boot_tgr <- function(x, which = c("distribution", "ci"),
       graphics::hist(group_means, main = paste("TGR:", g),
                      xlab = "Mean TGR", col = "lightblue", border = "white")
       graphics::abline(v = mean(x$tgr_original[idx]), col = "red", lwd = 2)
+      # Dashed lines at the group-level CI bounds (columns 3:4 of
+      # ci_group; names are level-dependent, e.g. "2.5%"/"97.5%")
+      if (!is.null(x$ci_group)) {
+        bounds <- unlist(x$ci_group[x$ci_group$Group == g, 3:4])
+        graphics::abline(v = bounds, col = "red", lwd = 1, lty = 2)
+      }
     }
   } else {
     # CI plot per group