#' @name methods-unitquantreg
#'
#' @title Methods for \code{unitquantreg} and \code{unitquantregs} objects
#'
#' @author André F. B. Menezes
#'
#' @param object,x fitted model object of class \code{\link{unitquantreg}}.
#' @param digits minimal number of _significant_ digits.
#' @param correlation logical; if \code{TRUE}, the correlation matrix of the estimated parameters is returned and printed. Default is \code{FALSE}.
#' @param parm a specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names. If missing, all parameters are considered.
#' @param level the confidence level required.
#' @param type character indicating type of fitted values to return.
#' @param formula. Changes to the formula see \code{\link{update.formula}} for details.
#' @param formula an R formula.
#' @param evaluate If true evaluate the new call else return the call.
#' @param ... additional argument(s) for methods. Currently not used.
#'
#' @description Methods for extracting information from fitted regression models
#' objects of class \code{unitquantreg} and \code{unitquantregs}.
#'
#' @return
#' The \code{summary} method gives Wald tests for the regressions coefficients
#' based on the observed Fisher information matrix. As usual the \code{summary}
#' method returns a list with relevant model statistics and estimates, which
#' can be printed using the \code{print} method.
#'
#' The \code{coef}, \code{vcov}, \code{confint} and \code{fitted} methods can
#' be use to extract, respectively, the estimated coefficients, the
#' estimated covariance matrix, the Wald confidence intervals, and fitted
#' values.
#'
#' A \code{\link{logLik}} method is also provide, then the \code{\link{AIC}}
#' function can be use to calculated the Akaike Information Criterion.
#'
#' The generic methods \code{\link{terms}}, \code{\link{model.frame}},
#' \code{\link{model.matrix}}, \code{\link{update}} and are also provided.
#'
#' @examples
#' data(sim_bounded, package = "unitquantreg")
#' sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ]
#' fit_1 <- unitquantreg(formula = y1 ~ x + z + I(x^2) | z + x,
#' data = sim_bounded_curr,
#' family = "uweibull",
#' tau = 0.5, link.theta = "log")
#' fit_1
#' summary(fit_1)
#' vcov(fit_1)
#' coef(fit_1)
#' confint(fit_1)
#' terms(fit_1)
#' model.frame(fit_1)[1:5, ]
#' model.matrix(fit_1)[1:5, ]
#' update(fit_1, . ~ . -x)
#' update(fit_1, . ~ . -z)
#' update(fit_1, . ~ . -I(x^2))
#' update(fit_1, . ~ . | . -z)
#' update(fit_1, . ~ . | . -x)
#'
#' @importFrom stats pnorm cov2cor coef vcov printCoefmat update formula
#'
#' @rdname methods-unitquantreg
#' @export
print.unitquantreg <- function(x, digits = max(4, getOption("digits") - 3), ...) {
cat("\n", x$family, " quantile regression model \n", sep = "")
cat("\nCall: ", paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n", sep = "")
cat("Mu coefficients (quantile model with ", x$link$mu$name, " link and tau = ", x$tau, "): \n", sep = "")
print.default(.FF(x$coefficients$mu, Digits = digits), print.gap = 2, quote = FALSE)
cat("\n")
if (x$theta_const) {
cat("Model with constant shape parameter:", "\n", sep = "")
print.default(.FF(x$coefficients$theta, Digits = digits), print.gap = 2,
quote = FALSE)
cat("\n")
} else {
names(x$coefficients$theta) <- gsub("(theta)_", "", names(x$coefficients$theta), fixed = TRUE)
cat("Theta coefficients (shape model with ", x$link$theta$name, " link):", "\n", sep = "")
print.default(.FF(x$coefficients$theta, Digits = digits), print.gap = 2, quote = FALSE)
cat("\n")
}
invisible(x)
}
# Summary -----------------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
summary.unitquantreg <- function(object, correlation = FALSE, ...) {
cf <- object$coefficients
names(cf) <- NULL
estimates <- unlist(cf)
stderror <- sqrt(diag(object$vcov))
zvalue <- estimates/stderror
pvalue <- 2 * pnorm(-abs(zvalue))
table <- cbind("Estimate" = estimates, "Std. Error" = stderror,
"Z value" = zvalue, "Pr(>|z|)" = pvalue)
if (correlation) {
correlation <- cov2cor(object$vcov)
}
out <- list(coeftable = table,
loglik = object$loglik,
df.residual = object$df.residual,
correlation = correlation,
call = object$call,
iterations = object$iterations,
tau = object$tau,
family = object$family,
link = object$link,
dims = c(length(cf[[1L]]), length(cf[[2L]])),
theta_const = object$theta_const)
class(out) <- "summary.unitquantreg"
out
}
# Print output summary ----------------------------------------------------
#' @export
print.summary.unitquantreg <- function(x, digits = max(4, getOption("digits") - 3), ...) {
p <- x$dims[[1L]]
q <- x$dims[[2L]]
cat("\n Wald-tests for ", x$family, " quantile regression model", "\n",
sep = "")
cat("\nCall: ", paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n",
sep = "")
cat("Mu coefficients: (quantile model with ", x$link$mu$name, " link and tau = ",
x$tau, "): \n", sep = "")
printCoefmat(x$coeftable[1:p, , drop = FALSE], digits = digits,
has.Pvalue = TRUE)
cat("\n")
if (x$theta_const) {
cat("Model with constant shape:", "\n", sep = "")
printCoefmat(x$coeftable[-(1:p), , drop = FALSE], digits = digits, has.Pvalue = TRUE)
cat("\n")
} else {
rownames(x$coeftable)[-(1:p)] <- gsub("(theta)_", "", rownames(x$coeftable)[-(1:p)], fixed = TRUE)
cat("Theta coefficients (shape model with ", x$link$theta$name, " link):", "\n",
sep = "")
printCoefmat(x$coeftable[-(1:p), , drop = FALSE], digits = digits,
has.Pvalue = TRUE)
cat("\n")
}
if (is.matrix(x$correlation)) {
cat("Correlation of coefficients:", "\n", sep = "")
corr <- x$correlation
corr <- format(round(corr, 2L), nsmall = 2L, digits = digits)
corr[!lower.tri(corr)] <- ""
print(corr[-1, -ncol(corr), drop = FALSE], quote = FALSE)
cat("\n")
}
cat("Residual degrees of freedom: ", x$df.residual, "\n", sep = "")
cat("Log-likelihood: ", x$loglik, "\n", sep = "")
cat("Number of iterations: ", x$iterations, "\n", sep = "")
invisible(x)
}
# coef function -----------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
coef.unitquantreg <- function(object, type = c("full", "quantile", "shape"), ...) {
if (!missing(...)) {
warning("Extra arguments discarded")
}
cf <- object$coefficients
type <- match.arg(type)
out <- switch (type,
"full" = {
names(cf) <- NULL
unlist(cf)
},
"quantile" = cf$mu,
"shape" = cf$theta
)
out
}
# vcov function -----------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
vcov.unitquantreg <- function(object, ...) {
if (!missing(...)) {
warning("Extra arguments discarded")
}
object$vcov
}
# logLik function ---------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
logLik.unitquantreg <- function(object, ...) {
if (!missing(...)) {
warning("Extra arguments discarded")
}
ll <- object$loglik
attr(ll, "df") <- object$npar
attr(ll, "nobs") <- object$nobs
class(ll) <- "logLik"
ll
}
# confint function --------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
confint.unitquantreg <- function(object, parm, level = 0.95, ...)
{
cf <- coef(object)
ses <- sqrt(diag(vcov(object)))
pnames <- names(ses)
if (missing(parm)) {
parm <- pnames
}
else if (is.numeric(parm)) {
parm <- pnames[parm]
}
a <- (1 - level) * 0.5
a <- c(a, 1 - a)
fac <- qnorm(a)
pct <- .format_perc(a, 3)
ci <- array(NA_real_, dim = c(length(parm), 2L),
dimnames = list(parm, pct))
ci[] <- cf[parm] + ses[parm] %o% fac
ci
}
# Fitted values -----------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
fitted.unitquantreg <- function(object, type = c("quantile", "shape", "full"), ...) {
if (!missing(...)) {
warning("Extra arguments discarded")
}
type <- match.arg(type)
switch (type,
"all" = object$fitted.values,
"quantile" = object$fitted.values$mu,
"shape" = object$fitted.values$theta
)
}
# Terms -------------------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
terms.unitquantreg <- function(x, type = c("quantile", "shape"), ...) {
x$terms[[match.arg(type)]]
}
#' @rdname methods-unitquantreg
#' @export
model.frame.unitquantreg <- function(formula, ...) {
# if (!is.null(formula$model)) return(formula$model)
dots <- list(...)
nargs <- dots[match(c("data", "na.action", "subset"), names(dots), 0L)]
if (length(nargs) || is.null(formula$model)) {
fcall <- formula$call
fcall$method <- "model.frame"
## need unitquantreg:: for non-standard evaluation
fcall[[1L]] <- quote(unitquantreg::unitquantreg)
fcall[names(nargs)] <- nargs
env <- environment(formula$terms)
if (is.null(env)) env <- parent.frame()
eval(fcall, env)
}
else formula$model
}
#' @rdname methods-unitquantreg
#' @export
model.matrix.unitquantreg <- function(object, type = c("quantile", "shape"), ...) {
type <- match.arg(type)
if (is.null(object$x[[type]])) {
model.matrix(object$terms[[type]], model.frame(object))
} else {
object$x[[type]]
}
}
# Update ------------------------------------------------------------------
#' @rdname methods-unitquantreg
#' @export
update.unitquantreg <- function (object, formula., ..., evaluate = TRUE)
{
call <- object$call
if (is.null(call)) stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$...
if (!missing(formula.)) {
call$formula <- formula(update(Formula::Formula(formula(object)), formula.))
}
if (length(extras)) {
existing <- !is.na(match(names(extras), names(call)))
for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
if (any(!existing)) {
call <- c(as.list(call), extras[!existing])
call <- as.call(call)
}
}
if (evaluate) eval(call, parent.frame())
else call
}
# Methods unitquantregs (tau vectorized) ----------------------------------
#' @rdname methods-unitquantreg
#' @export
print.unitquantregs <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("\n", x[[1L]]$family, " quantile regression model \n", sep = "")
cat("\nCall: ", paste(deparse(x[[1L]]$call), sep = "\n", collapse = "\n"), "\n\n", sep = "")
quant_coefs <- sapply(x, coef, type = "quantile")
shape_coefs <- sapply(x, coef, type = "shape")
taus <- .FF(sapply(x, "[[", "tau"), Digits = digits)
colnames(quant_coefs) <- paste0("tau = ", taus)
cat("Mu coefficients (quantile model with ", x[[1L]]$link$mu$name, " link): \n", sep = "")
print.default(.FF(quant_coefs, Digits = digits), print.gap = 2, quote = FALSE)
cat("\n")
cat("\n")
if (x[[1L]]$theta_const) {
shape_coefs <- matrix(shape_coefs, nrow = 1)
colnames(shape_coefs) <- paste0("tau = ", taus)
rownames(shape_coefs) <- names(coef(x[[1L]], type = "shape"))
cat("Model with constant shape parameter:", "\n", sep = "")
print.default(.FF(shape_coefs, Digits = digits), print.gap = 2,
quote = FALSE)
cat("\n")
} else {
colnames(shape_coefs) <- paste0("tau = ", taus)
cat("Theta coefficients (shape model with ", x[[1L]]$link$theta$name, " link):", "\n", sep = "")
print.default(.FF(shape_coefs, Digits = digits), print.gap = 2, quote = FALSE)
cat("\n")
}
invisible(x)
}
#' @rdname methods-unitquantreg
#' @export
summary.unitquantregs <- function(object, digits = max(3, getOption("digits") - 3), ...) {
out <- lapply(object, summary)
names(out) <- .FF(sapply(object, "[[", "tau"), Digits = digits)
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.