#-----------------------------------------------------------------------
#' @title Summary of the flexible count models (individual t-tests)
#' @param object an object of class \code{flexcm}, a result of call
#' \code{\link{flexcm}(...)}.
#' @param correlation logical; if \code{TRUE}, the correlation matrix of
#' the estimated parameters is returned and printed. Default is
#' \code{FALSE}.
#' @param ... Currently not used.
#' @return an object of class \code{"summary.flexcm"}, a list with
#' components.
#' @author Eduardo Jr <edujrrib@gmail.com>
#' @importFrom stats cov2cor pnorm printCoefmat
#' @export
#'
summary.flexcm <- function(object, correlation = FALSE, ...) {
if (is.null(object$vcov))
stop(paste("Refit the model with `hessian=TRUE` to compute",
"the standard errors."))
power <- object$power
powse <- attr(power, "std_error")
#------------------------------------------
estimates <- object$coefficients
stderror <- sqrt(diag(object$vcov))
zvalue <- estimates/stderror
pvalue <- 2 * pnorm(-abs(zvalue))
ctable <- cbind("Estimate" = estimates,
"Std. Error" = stderror,
"Z value" = zvalue,
"Pr(>|z|)" = pvalue)
if (is.numeric(powse)) {
pline <- c(power, powse, power/powse,
2 * pnorm(-abs(power/powse)))
ctable <- rbind("power" = pline, ctable)
}
#------------------------------------------
if (correlation) {
correlation <- cov2cor(object$vcov)
}
#--------------------------------------------
out <- list(model = object$model,
coeftable = ctable,
loglik = object$loglik,
df.residual = object$df.residual,
correlation = correlation,
power = object$power,
call = object$call)
class(out) <- "summary.flexcm"
return(out)
}
#-----------------------------------------------------------------------
# Print method for summary flexible count models models
#' @rdname flexcm-methods
#' @export
#'
print.summary.flexcm <- function(x,
digits = max(3L, getOption("digits") - 3L),
...) {
name <- `_get_model_name`(x$model)
index <- is.numeric(attr(x$power, "std_error"))
cat(sprintf("\nIndividual Wald-tests for %s regression models",
name), sep = "")
cat("\nCall: ",
paste(deparse(x$call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
start_mean <- ifelse(index, 3, 2)
cat("Mean coefficients:", "\n", sep = "")
printCoefmat(x$coeftable[start_mean:nrow(x$coeftable), ,
drop = FALSE],
digits = digits,
has.Pvalue = TRUE)
cat("\n")
start_disp <- ifelse(index, 2, 1)
cat("Dispersion coefficient:", "\n", sep = "")
printCoefmat(x$coeftable[start_disp, , drop = FALSE],
digits = digits,
has.Pvalue = TRUE)
cat("\n")
if (name == "Poisson-Tweedie") {
if (index) {
cat("Power coefficient:", "\n", sep = "")
printCoefmat(x$coeftable[1, , drop = FALSE],
digits = digits,
has.Pvalue = TRUE)
cat("\n")
} else {
cat("Power coefficient (fixed): power = ", x$power,
"\n\n", sep = "")
}
}
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("Minus twice the log-likelihood: ", -2 * x$loglik,
"\n", sep = "")
invisible(x)
}
#-----------------------------------------------------------------------
#' @name lrtests-flexcm
#' @title Likelihood ratio tests for nested flexible count models
#' @param object an object of class \code{flexcm}, a result of call
#' \code{\link{flexcm}(...)}.
#' @param heading logical; if \code{TRUE}, the model formulas are
#' printed as heading. Default is \code{TRUE}.
#' @param ... Currently not used.
#' @return an object of class \code{"anova"}, a table with compenents
#' for LR test.
#' @author Eduardo Jr <edujrrib@gmail.com>
#' @importFrom stats pchisq
#' @export
#'
anova.flexcm <- function(object, ..., heading) {
dots <- list(...)
isflexcm <- class(object) == "flexcm"
#------------------------------------------
# Organize the list of models
if (isflexcm & length(dots) == 0)
stop("Currently, this method only compare two (or more) models.")
if (!isflexcm & length(dots) == 0) mlist <- object
if (!isflexcm & length(dots) > 0) mlist <- c(object, dots)
if (isflexcm & length(dots) > 0) mlist <- c(list(object), dots)
if (missing(heading)) heading <- is.null(names(mlist))
#------------------------------------------
# Handle the class of the models
inherits_flexcm <- vapply(mlist, inherits, what = "flexcm", 0L)
if (any(!inherits_flexcm))
stop("Not all objects are of class \"flexcm\".")
mnames <- vapply(mlist, function(m) m$model, "")
if (!all(mnames[1L] == mnames)) {
stop("The models were fitted under different likelihoods.")
}
obs <- vapply(mlist, function(x) x$nobs, 0L)
if (diff(range(obs)) != 0) {
warning("The models were fitted for different sample sizes.")
}
forms <- vapply(mlist, function(x) deparse(x$formula), "")
#--------------------------------------------
# Compute stats
rds <- vapply(mlist, function(x) x$df.residual, 0L)
lls <- vapply(mlist, function(x) x$loglik, 0)
aic <- -2 * lls + 2 * (obs - rds)
bic <- -2 * lls + log(obs) * (obs - rds)
lrs <- c(NA, abs(diff(-2 * lls)))
dfs <- c(NA, abs(diff(rds)))
pvs <- pchisq(q = lrs, df = dfs, lower.tail = FALSE)
nam <- if (is.null(names(mlist))) {
sprintf("Model %i", seq_along(mlist))
} else names(mlist)
tab <- data.frame("Model" = nam,
"AIC" = aic,
"BIC" = bic,
"Resid.df" = rds,
"loglik" = lls,
"Chisq.df" = dfs,
"Chisq" = lrs,
"Pr(>Chisq)" = pvs,
check.names = FALSE,
stringsAsFactors = FALSE)
rownames(tab) <- NULL
#--------------------------------------------
# Heading
if (heading) {
diffpower <- vapply(mlist, function(x)
typeof(attr(x$power, "std_error")), "")
if (!all(diffpower[1] == diffpower)) {
powers <- vapply(mlist, function(x) x$power, 0)
heading <- paste(sprintf(" %s (power=%.2f): %s",
format(nam),
powers,
forms),
collapse = "\n")
} else {
heading <- paste(sprintf(" %s: %s",
format(nam),
forms),
collapse = "\n")
}
} else {
heading <- NULL
}
name <- `_get_model_name`(mnames[1L])
attr(tab, "name") <- name
attr(tab, "heading") <- heading
class(tab) <- c("anova.flexcm", "data.frame")
return(tab)
}
#' @rdname lrtests-flexcm
#' @export
anova.list <- function(object, ..., heading) {
out <- anova.flexcm(object, ..., heading = heading)
return(out)
}
#' @export
print.anova.flexcm <- function(x,
digits = max(getOption("digits") - 2L, 3L),
signif.stars = getOption("show.signif.stars"),
...) {
name <- attr(x, "name")
heading <- attr(x, "heading")
cat(sprintf("\nLikelihood ratio test for %s", name),
"regression models", "\n\n")
if (!is.null(heading)) {
cat(heading, "\n\n")
}
xx <- x
rownames(xx) <- x$Model
printCoefmat(xx[, -1],
cs.ind = NA,
tst.ind = c(3L, 5L),
zap.ind = 7L,
na.print = "",
digits = digits,
signif.stars = signif.stars)
invisible(x)
}
#' @export
#' @rdname lrtests-flexcm
lrtest <- function (object, ..., heading) {
UseMethod("anova", object)
}
#-----------------------------------------------------------------------
#' @title Likelihood ratio test for equidispersion
#' @param object an object (or a list of objects) of class
#' \code{"flexcm"} with a model nested to Poisson (i.e. Com-Poisson,
#' Gamma-count, Generalized-Poisson, Double Poisson and
#' Poissson-Tweedie).
#' @param ... possibly, others fitted models (objects of class
#' \code{"flexcm"}).
#' @param heading logical; if \code{TRUE}, the model formulas are
#' printed as heading.
#' @details Note that, when more than one model is passed the models
#' \strong{are not compared among themselves}. Here, each model is
#' compared to Poisson and them organized in a table.
#' @author Eduardo Jr <edujrrib@gmail.com>
#' @export
#'
equitest <- function (object, ..., heading) {
UseMethod("equitest", object)
}
#' @export
equitest.flexcm <- function(object, ..., heading) {
dots <- list(...)
isflexcm <- class(object) == "flexcm"
#------------------------------------------
# Organize the list of models
if (!isflexcm & is.null(dots)) mlist <- object
if ( isflexcm & is.null(dots)) mlist <- list(object)
if (!isflexcm & !is.null(dots)) mlist <- c(object, dots)
if ( isflexcm & !is.null(dots)) mlist <- c(list(object), dots)
if (any(!vapply(mlist, inherits, what = "flexcm", 0L)))
stop("Not all objects are of class \"flexcm\".")
if (missing(heading))
heading <- ifelse(length(mlist) == 1,
FALSE,
is.null(names(mlist)))
fam <- vapply(mlist, function(x) `_get_model_name`(x$model), "")
if ("Discrete Weibull" %in% fam) {
stop("Discrete Weibull models are not nested to Poisson.")
}
#--------------------------------------------
lpo <- vapply(mlist, function(x) x$poissonfit$loglik, 0)
npo <- vapply(mlist, function(x) length(x$poissonfit$coef), 0L)
rds <- vapply(mlist, function(x) x$df.residual, 0L)
lls <- vapply(mlist, function(x) x$loglik, 0)
obs <- vapply(mlist, function(x) x$nobs, 0L)
lrs <- 2 * (lls - lpo)
dfs <- obs - rds - npo
pvs <- pchisq(q = lrs, df = dfs, lower.tail = FALSE)
nam <- if (is.null(names(mlist))) {
sprintf("Model %i", seq_along(mlist))
} else names(mlist)
tab <- data.frame("Family" = fam,
"Model" = nam,
"Resid.df" = rds,
"Loglik" = lls,
"LRT_stat" = lrs,
"LRT_df" = dfs,
"Pr(>LRT_stat)" = pvs,
check.names = FALSE,
stringsAsFactors = FALSE)
rownames(tab) <- NULL
#--------------------------------------------
if (heading) {
calls <- vapply(mlist, function(x) deparse(x$call), "")
cnames <- paste0(" ", format(nam), ": ")
calls <- paste0(cnames, calls)
calls <- paste(calls, collapse = "\n")
} else {
calls <- NULL
}
attr(tab, "heading") <- calls
class(tab) <- c("equitest.flexcm", "data.frame")
return(tab)
}
#' @export
equitest.list <- function(object, ..., heading) {
out <- equitest.flexcm(object, ..., heading = heading)
return(out)
}
#' @export
print.equitest.flexcm <- function(x,
digits = max(getOption("digits") - 2L, 3L),
signif.stars = getOption("show.signif.stars"),
...) {
xx <- x
heading <- attr(x, "heading")
cat("\nLikelihood ratio test for equidispersion", "\n\n")
if (!is.null(heading)) {
if (nrow(x) == 1L) {
heading <- gsub("Model 1:", "", heading)
}
cat(heading, "\n\n")
}
if (nrow(x) > 1L) {
rownames(xx) <- paste(format(x$Family), format(x$Model))
} else {
rownames(xx) <- format(x$Family)
}
printCoefmat(xx[, -(1:2)],
cs.ind = NA,
tst.ind = c(2L, 3L),
zap.ind = 5L,
na.print = "",
digits = digits,
signif.stars = signif.stars,
...)
invisible(x)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.