## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  eval      = TRUE
)

knitr::knit_hooks$set(output = function(x, options) {
  if (isTRUE(options$fold.output)) {
    paste0(
      "
<details>
  <summary>Toggle to see the output</summary>

```plaintext
",
      x,
      "```

</details>
"
    )
  } else {
    paste0("
```plaintext
", x, "```
")
  }
})

## ----fit, fold.output=TRUE----------------------------------------------------
library(smfa)
data("ricephil", package = "sfaR")

ricephil$group <- cut(
  ricephil$AREA,
  breaks        = quantile(ricephil$AREA, probs = c(0, 1/3, 2/3, 1), na.rm = TRUE),
  labels        = c("small", "medium", "large"),
  include.lowest = TRUE
)

meta_lp <- smfa(
  formula    = log(PROD) ~ log(AREA) + log(LABOR) + log(NPK),
  data       = ricephil,
  group      = "group",
  S          = 1,
  udist      = "hnormal",
  groupType  = "sfacross",
  metaMethod = "lp"
)

## ----summary, fold.output=TRUE------------------------------------------------
summary(meta_lp)

## ----efficiencies, fold.output=TRUE-------------------------------------------
eff <- efficiencies(meta_lp)
head(eff)

## ----subset, fold.output=TRUE-------------------------------------------------
# All small farms
eff_small <- eff[eff$group == "small", ]

# Descriptive statistics
summary(eff_small[, c("TE_group_BC", "TE_meta_BC", "MTR_BC")])

# Group-level means
aggregate(cbind(TE_group_BC, TE_meta_BC, MTR_BC) ~ group, data = eff, FUN = mean)

## ----plot, fold.output=TRUE---------------------------------------------------
# MTR distribution by group (base R)
boxplot(MTR_BC ~ group, data = eff,
        main = "Metatechnology Ratio by Farm Size",
        xlab = "Group", ylab = "MTR (BC)",
        col  = c("#e8f5ef", "#b2dfdb", "#80cbc4"))

# Histogram of TE_meta_BC
hist(eff$TE_meta_BC, breaks = 30,
     main = "Distribution of Metafrontier TE",
     xlab = "TE_meta_BC", col = "#2d8f65", border = "white")

## ----coef, fold.output=TRUE---------------------------------------------------
# First fit a QP model
meta_qp <- smfa(
  formula    = log(PROD) ~ log(AREA) + log(LABOR) + log(NPK),
  data       = ricephil, group = "group", S = 1, udist = "hnormal",
  groupType  = "sfacross", metaMethod = "qp"
)
coef(meta_qp)

## ----vcov, fold.output=TRUE---------------------------------------------------
vcov(meta_qp)

## ----loglik, fold.output=TRUE-------------------------------------------------
logLik(meta_lp)

## ----ic-----------------------------------------------------------------------
ic(meta_lp)
#>        AIC       BIC      HQIC
#> 1  184.579  253.710   212.113

## ----nobs---------------------------------------------------------------------
nobs(meta_lp)  # Total observations across all groups

## ----fitted, fold.output=TRUE-------------------------------------------------
fv <- fitted(meta_lp)
head(fv)

## ----residuals, fold.output=TRUE----------------------------------------------
res <- residuals(meta_lp)
head(res)

## ----compare, fold.output=TRUE------------------------------------------------
meta_lp    <- smfa(log(PROD) ~ log(AREA) + log(LABOR) + log(NPK),
                               data = ricephil, group = "group", S = 1,
                               udist = "hnormal", groupType = "sfacross",
                               metaMethod = "lp")
meta_qp    <- smfa(log(PROD) ~ log(AREA) + log(LABOR) + log(NPK),
                               data = ricephil, group = "group", S = 1,
                               udist = "hnormal", groupType = "sfacross",
                               metaMethod = "qp")
meta_huang <- smfa(log(PROD) ~ log(AREA) + log(LABOR) + log(NPK),
                               data = ricephil, group = "group", S = 1,
                               udist = "hnormal", groupType = "sfacross",
                               metaMethod = "sfa", sfaApproach = "huang")

# Combined information criteria table
models <- list(LP = meta_lp, QP = meta_qp, Huang = meta_huang)
do.call(rbind, lapply(names(models), function(nm) {
  ic_vals <- ic(models[[nm]])
  data.frame(Model = nm, AIC = ic_vals[["AIC"]],
             BIC = ic_vals[["BIC"]], HQIC = ic_vals[["HQIC"]])
}))