R/summary.snreg.R

In snreg: Regression with Skew-Normally Distributed Error Term

#' Summary Method for snreg Objects
#'
#' @title Summary for Skew-Normal Regression Models
#'
#' @description
#' Produces a summary object for objects of class \code{"snreg"}.
#' The function assigns the class \code{"summary.snreg"} to the fitted model
#' object, enabling a dedicated print method (\code{print.summary.snreg}) to
#' display results in a structured format.
#'
#' @param object
#' an object of class \code{"snreg"}, typically returned by \code{\link{snreg}}.
#'
#' @param ...
#' additional arguments (currently not used).
#'
#' @details
#' This method expects a fitted \code{"snreg"} object.
#' 
#' \code{summary.snreg} does not modify the contents of the object; it only
#' updates the class attribute to \code{"summary.snreg"}. The corresponding
#' print method (\code{\link{print.summary.snreg}}) is responsible for
#' formatting and displaying estimation details, such as convergence criteria,
#' log-likelihood, coefficient tables, and (if present) heteroskedastic and
#' skewness components.
#'
#' @return
#' An object of class \code{"summary.snreg"}, identical to the input \code{object}
#' except for its class attribute.
#'
#' @seealso
#' \code{\link{snreg}}, \code{\link{print.summary.snreg}}
#'
#' @examples
#' library(snreg)
#' 
#' data("banks07")
#' head(banks07)
#' 
#' # Translog cost function specification
#' 
#' spe.tl <- log(TC) ~ (log(Y1) + log(Y2) + log(W1) + log(W2))^2 +
#'   I(0.5 * log(Y1)^2) + I(0.5 * log(Y2)^2) +
#'   I(0.5 * log(W1)^2) + I(0.5 * log(W2)^2)
#' 
#' # Specification 1: homoskedastic noise and skewness
#' 
#' formSV <- NULL   # variance equation; constant variance
#' formSK <- NULL   # skewness equation; constant skewness
#' 
#' m1 <- snreg(
#'   formula  = spe.tl,
#'   data     = banks07,
#'   ln.var.v = formSV,
#'   skew.v   = formSK
#' )
#' 
#' summary(m1)
#'
#' @export
# summary.snreg <- function( object, ...) {
#   class(object) <- c("summary.snreg")
#   return(object)
# }

summary.snreg <- function(object, ...) {
  s <- list(
    # all your precomputed summary bits as above...
    original_class = class(object)
  )
  class(s) <- c("summary.snreg","summary")
  attr(s, "object") <- object  # attached but not printed by default
  s
}


#' Print Method for Summary of snreg Objects
#'
#' @title Print Summary of snreg Results
#'
#' @description
#' Prints the contents of a \code{"summary.snreg"} object in a structured
#' format. The method reports convergence status (based on gradient-Hessian
#' scaling), log-likelihood, estimation results, and—when present—summaries
#' for technical/cost efficiencies and marginal effects.
#'
#' @param x
#' an object of class \code{"summary.snreg"} (produced by \code{\link{summary.snreg}}).
#'
#' @param digits
#' integer indicating the number of digits to print; default \code{NULL}
#' (internally set to 4).
#'
#' @param ...
#' additional arguments (currently unused).
#'
#' @details
#'
#' This method expects a fitted \code{"snreg"} object.
#'
#' @return
#' The input \code{obj} is returned (invisibly) after printing.
#'
#' @seealso
#' \code{\link{summary.snreg}}
#'
#' @export
print.summary.snreg <- function( x, digits = NULL, ... ) {

  model <- attr(x, "object")       # the original snreg/snsf/lm.mle object
  cls   <- x$original_class         # character vector of original classes

  # Defensive: handle empty results gracefully when computing max.name.length
  rn <- tryCatch(row.names(model$results), error = function(e) NULL)
  max.name.length <- if (!is.null(rn) && length(rn)) max(nchar(rn)) else 12
  
  if( is.null( digits ) ) {
    digits <- 4
  }
  

  # Compare only if both are present, numeric, length-1, and finite
  if (!is.null(model[["gHg"]]) && !is.null(model[["lmtol"]]) &&
      is.numeric(model[["gHg"]]) && length(model[["gHg"]]) == 1 && is.finite(model[["gHg"]]) &&
      is.numeric(model[["lmtol"]]) && length(model[["lmtol"]]) == 1 && is.finite(model[["lmtol"]])) {
    
    if(model$gHg < model$lmtol){
      cat("\nConvergence given g*inv(H)*g' = ",
          formatC(model$gHg, digits = 1, format = "e"),
          " < lmtol(", model$lmtol, ")\n", sep = "")
      # cat(" was reached in ",model$counts[2]," iteration(s)\n", sep = "")
    }
  else {
    cat("\nCriterion g*inv(H)*g' = ",
        formatC(model$LM, digits = 1, format = "e"),
        " > lmtol(", model$lmtol, ")\n", sep = "")
    
    if (!is.null(model[["bhhh"]]) && 
        is.numeric(model[["bhhh"]]) && 
        length(model[["bhhh"]]) == 1 && 
        is.finite(model[["bhhh"]])) {
      
    if(model$bhhh){
      cat("Note that Hessian is computed as outer product (BHHH)\n", sep = "")
      cat("Criterion g'g = ", model$gg, "\n", sep = "")
    }
    }
    warning("Convergence given g*inv(H)*g' is still not reached; one of optim's convergence criteria is used",
            call. = FALSE)
  }
  }
  .timing(model$esttime, "Log likelihood maximization completed in ")
  cat("Log likelihood = ", formatC(model$ll, digits = 4, format = "f"), "\n",  sep = "")
  # cat("____________________________________________________\n")
  


  
  if ("snsf" %in% cls) {
    # snsf --------------------------------------------------------------------
    
    output <- model$ctab
    colnames(output) <- c("Coef.", "SE ", "z ",  "P>|z|")
    max.name.length <- max(nchar(names(output[,1]) ))
    
    distribution <- model$distribution
    k  <- model$K
    ksv <- model$Ksv
    ksk <- model$Ksk
    ksu <- model$Ksu
    kmu <- model$Kmu
    prod <- model$prod
    n <- model$n
    
    
    cat("",rep("_", max.name.length+42-1),"", "\n", sep = "")
    if(prod){
      cat("\nCross-sectional stochastic (production) frontier model\n",  sep = "")}
    else {
      cat("\nCross-sectional stochastic (cost) frontier model\n",  sep = "")
    }
    cat("\nDistributional assumptions\n\n", sep = "")
    Assumptions <- rep("heteroskedastic",2)
    if(ksv==1){
      Assumptions[1] <- "homoskedastic"
    }
    if(ksu==1){
      Assumptions[2] <- "homoskedastic"
    }
    Distribution = c("skew normal ", "half-normal ")
    if(distribution == "t"){
      Distribution[2] <- "truncated-normal "
    }   
    if(distribution == "e"){
      Distribution[2] <- "exponential "
    }
    a1 <- data.frame(
      Component = c("Random noise: ","Inefficiency: "),
      Distribution = Distribution,
      Assumption = Assumptions
    )
    print(a1, quote = FALSE, right = FALSE)
    cat("\nNumber of observations = ", n, "\n", sep = "")
    # max.name.length <- max(nchar(row.names(output)))
    est.rez.left <- floor( (max.name.length+42-22) / 2 )
    est.rez.right <- max.name.length+42-22 - est.rez.left
    cat("\n",rep("-", est.rez.left)," Estimation results: ",rep("-", est.rez.right),"\n\n", sep ="")
    # cat("\n--------------- Estimation results: --------------\n\n", sep = "")
    .printoutcs(output, digits = digits, k = k, ksv = ksv, ksk = ksk, ksu = ksu, kmu = kmu, na.print = "NA", dist = distribution, max.name.length = max.name.length)
    
    
    
    
    if(model$prod){
      myeff <- "technical"
      ndots <- "...."
    }
    else {
      myeff <- "cost"
      ndots <- "........."
    }
    
    sum.te <- paste0(" Summary of ", myeff, " efficiencies: ")
    est.cle.left  <- floor( (max.name.length + 42 - nchar(sum.te) - 1) / 2 )
    est.cle.right <- max.name.length + 42 - nchar(sum.te) - 1 - est.cle.left
    if(est.cle.left  <= 0) est.cle.left  <- 1
    if(est.cle.right <= 0) est.cle.right <- 1
    cat("\n", rep("-", est.cle.left), sum.te, rep("-", est.cle.right), "\n\n", sep ="")
    .su(model$ eff, transpose = TRUE, print = TRUE, names = "TE_JLMS")
  } 
  else if  ("snreg" %in% cls){
    # snreg ------------------------------------------------------------------
    output <- model$ctab
    colnames(output) <- c("Coef.", "SE ", "z ",  "P>|z|")
    max.name.length <- max(nchar(names(output[,1]) ))
    
    k  <- model$K
    ksv <- model$Ksv
    ksk <- model$Ksk
    n <- model$n
    
    cat("\nNumber of observations = ", n, "\n", sep = "")
    # max.name.length <- max(nchar(row.names(output)))
    est.rez.left <- floor( (max.name.length+42-22) / 2 )
    est.rez.right <- max.name.length+42-22 - est.rez.left
    cat("\n",rep("-", est.rez.left)," Estimation results: ",rep("-", est.rez.right),"\n\n", sep ="")
    # cat("\n--------------- Estimation results: --------------\n\n", sep = "")
    .printoutcs(output, digits = digits, k = k, ksv = ksv, ksk = ksk, ksu = 0, kmu = 0, na.print = "NA", dist = distribution, max.name.length = max.name.length)
    
  } else if  ("lm.mle" %in% cls){
    # lm.mle ------------------------------------------------------------------
    output <- model$ctab
    colnames(output) <- c("Coef.", "SE ", "z ",  "P>|z|")
    max.name.length <- max(nchar(names(output[,1]) ))
    
    
    ksv <- model$Ksv
    k  <- model$K - ksv # (since theta contains them all)
    n <- model$n
    ksk <- 0
    
    cat("\nNumber of observations = ", n, "\n", sep = "")
    # max.name.length <- max(nchar(row.names(output)))
    est.rez.left <- floor( (max.name.length+42-22) / 2 )
    est.rez.right <- max.name.length+42-22 - est.rez.left
    cat("\n",rep("-", est.rez.left)," Estimation results: ",rep("-", est.rez.right),"\n\n", sep ="")
    # cat("\n--------------- Estimation results: --------------\n\n", sep = "")
    .printoutcs(output, digits = digits, k = k, ksv = ksv, ksk = ksk, ksu = 0, kmu = 0, na.print = "NA", dist = distribution, max.name.length = max.name.length)
  } else {
    stop("Unknown model\n")
  }
  
  invisible( x )
}