Put in one file all S3 generic methods.

R/mc_S3_methods.R

+#' @title ANOVA method for McGLMs.
+#' @name anova.mcglm
+#'
+#' @description ANOVA method for object of class McGLMS.
+#'
+#' @param object an object of class \code{mcglm}, usually, a result of a
+#'     call to \code{mcglm()}.
+#' @param ... additional arguments affecting the summary produced. Note
+#'     that there is no extra options for mcglm object class.
+#'
+#' @return A \code{data.frame} with Chi-square statistic to test the
+#'     null hypothesis of a parameter, or a set of parameters, be
+#'     zero. The Wald test based on the observed covariance matrix of
+#'     the parameters is used.
+#'
+#' @method anova mcglm
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#' @export
+
+anova.mcglm <- function(object, ...) {
+    n_resp <- length(object$mu_list)
+    n_beta <- lapply(object$list_X, ncol)
+    idx.list <- list()
+    for (i in 1:n_resp) {
+        idx.list[[i]] <- rep(i, n_beta[i])
+    }
+    vv <- vcov(object)
+    n_par <- dim(vv)[1]
+    idx.vec <- do.call(c, idx.list)
+    n_cov <- n_par - length(idx.vec)
+    idx.vec <- c(idx.vec, rep(0, n_cov))
+    temp.vcov <- list()
+    temp.beta <- list()
+    for (i in 1:n_resp) {
+        idx.id = idx.vec == i
+        temp.vcov[[i]] <- vv[idx.id, idx.id]
+        temp.beta[[i]] <-
+            coef(object, type = "beta", response = i)$Estimates
+    }
+    saida <- list()
+    for (i in 1:n_resp) {
+        idx <- attr(object$list_X[[i]], "assign")
+        names <- colnames(object$list_X[[i]])
+        if (names[1] == "(Intercept)") {
+            idx <- idx[-1]
+            names <- names[-1]
+            temp.beta[[i]] <- temp.beta[[i]][-1]
+            temp.vcov[[i]] <- temp.vcov[[i]][-1, -1]
+        }
+        n_terms <- length(unique(idx))
+        X2.resp <- list()
+        for (j in 1:n_terms) {
+            idx.TF <- idx == j
+            temp <- as.numeric(
+                t(temp.beta[[i]][idx.TF]) %*%
+                    solve(as.matrix(temp.vcov[[i]])[idx.TF, idx.TF]) %*%
+                    temp.beta[[i]][idx.TF])
+            nbeta.test <- length(temp.beta[[i]][idx.TF])
+            X2.resp[[j]] <- data.frame(
+                Covariate = names[idx.TF][1],
+                Chi.Square = round(temp, 4),
+                Df = nbeta.test,
+                p.value = round(pchisq(temp, nbeta.test,
+                                       lower.tail = FALSE), 4))
+        }
+        saida[[i]] <- do.call(rbind, X2.resp)
+    }
+    cat("Wald test for fixed effects\n")
+    return(saida)
+}
+
+#' @title Extract model coefficients for mcglm class
+#' @name coef.mcglm
+#'
+#' @description \code{coef.mcglm} is a function which extracts model
+#'     coefficients from objects of \code{mcglm} class.
+#'
+#' @param object An object of \code{mcglm} class.
+#' @param std.error Logical. If \code{TRUE} returns the standard errors
+#'     of the estimates. Default is \code{FALSE}.
+#' @param response A numeric vector specyfing for which response
+#'     variables the coefficients should be returned.
+#' @param type A string vector (can be 1 element length) specyfing which
+#'     coefficients should be returned. Options are \code{"beta"},
+#'     \code{"tau"}, \code{"power"}, \code{"tau"} and
+#'     \code{"correlation"}.
+#' @param ... additional arguments affecting the summary produced. Note
+#'     that there is no extra options for mcglm object class.
+#'
+#' @return A \code{data.frame} with parameters names, estimates,
+#'     response number and parameters type.
+#'
+#' @method coef mcglm
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#' @export
+
+coef.mcglm <- function(object, std.error = FALSE,
+                       response = c(NA, 1:length(object$beta_names)),
+                       type = c("beta", "tau", "power", "correlation"),
+                       ...) {
+    n_resp <- length(object$beta_names)
+    cod_beta <- list()
+    cod_power <- list()
+    cod_tau <- list()
+    type_beta <- list()
+    type_power <- list()
+    type_tau <- list()
+    resp_beta <- list()
+    resp_power <- list()
+    resp_tau <- list()
+    response_for <- 1:n_resp
+    for (i in response_for) {
+        cod_beta[[i]] <- paste0(
+            paste0("beta", i), 0:c(object$Information$n_betas[[i]] - 1))
+        type_beta[[i]] <- rep("beta", length(cod_beta[[i]]))
+        resp_beta[[i]] <- rep(response_for[i], length(cod_beta[[i]]))
+        if (object$Information$n_power[[i]] != 0 |
+            object$power_fixed[[i]] == FALSE) {
+            cod_power[[i]] <- paste0(
+                paste0("power", i), 1:object$Information$n_power[[i]])
+            type_power[[i]] <- rep("power",
+                                   length(cod_power[[i]]))
+            resp_power[[i]] <- rep(response_for[i],
+                                   length(cod_power[[i]]))
+        }
+        if (object$Information$n_power[[i]] == 0) {
+            cod_power[[i]] <- rep(1, 0)
+            type_power[[i]] <- rep(1, 0)
+            resp_power[[i]] <- rep(1, 0)
+        }
+        cod_tau[[i]] <- paste0(
+            paste0("tau", i), 1:object$Information$n_tau[[i]])
+        type_tau[[i]] <- rep("tau", length(cod_tau[[i]]))
+        resp_tau[[i]] <- rep(response_for[i], length(cod_tau[[i]]))
+    }
+    rho_names <- c()
+    if (n_resp != 1) {
+        combination <- combn(n_resp, 2)
+        for (i in 1:dim(combination)[2]) {
+            rho_names[i] <- paste0(
+                paste0("rho", combination[1, i]), combination[2, i])
+        }
+    }
+    type_rho <- rep("correlation", length(rho_names))
+    resp_rho <- rep(NA, length(rho_names))
+    cod <- c(do.call(c, cod_beta), rho_names,
+             do.call(c, Map(c, cod_tau)))
+    type_cod <- c(do.call(c, type_beta), type_rho,
+                  do.call(c, Map(c, type_tau)))
+    response_cod <- c(do.call(c, resp_beta), resp_rho,
+                      do.call(c, Map(c, resp_tau)))
+
+    if (length(cod_power) != 0) {
+        cod <- c(do.call(c, cod_beta), rho_names,
+                 do.call(c, Map(c, cod_power, cod_tau)))
+        type_cod <- c(do.call(c, type_beta), type_rho,
+                      do.call(c, Map(c, type_power, type_tau)))
+        response_cod <- c(do.call(c, resp_beta), resp_rho,
+                          do.call(c, Map(c, resp_power, resp_tau)))
+    }
+
+    Estimates <- c(object$Regression, object$Covariance)
+    coef_temp <- data.frame(
+        Estimates = Estimates,
+        Parameters = cod,
+        Type = type_cod,
+        Response = response_cod)
+    if (std.error == TRUE) {
+        coef_temp <- data.frame(
+            Estimates = Estimates,
+            Std.error = sqrt(diag(object$vcov)),
+            Parameters = cod, Type = type_cod,
+            Response = response_cod)
+    }
+    output <- coef_temp[
+        which(coef_temp$Response %in% response &
+                  coef_temp$Type %in% type), ]
+    return(output)
+}
+
+#' @title Confidence Intervals for mcglm
+#' @name confint.mcglm
+#'
+#' @description Computes confidence intervals for parameters in a fitted
+#'     \code{mcglm} model.
+#'
+#' @param object a fitted \code{mcglm} object.
+#' @param parm a specification of which parameters are to be given
+#'     confidence intervals, either a vector of number or a vector of
+#'     strings. If missing, all parameters are considered.
+#' @param level the nominal confidence level.
+#' @param ... additional arguments affecting the confidence interval
+#'     produced. Note that there is no extra options for \code{mcglm}
+#'     object class.
+#'
+#' @return A \code{data.frame} with confidence intervals, parameters
+#'     names, response number and parameters type.
+#'
+#' @method confint mcglm
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#' @export
+
+confint.mcglm <- function(object, parm, level = 0.95, ...) {
+    temp <- coef(object, std.error = TRUE)
+    if (missing(parm)) {
+        parm <- 1:length(temp$Estimates)
+    }
+    a <- (1 - level)/2
+    a <- c(a, 1 - a)
+    fac <- qnorm(a)
+    ci <- temp$Estimates + temp$Std.error %o% fac
+    colnames(ci) <- paste0(format(a, 2), "%")
+    rownames(ci) <- temp$Parameters
+    return(ci[parm])
+}
+#' @title Extract Model Fitted Values of McGLM
+#' @name fitted.mcglm
+#'
+#' @description Extract fitted values for objects of \code{mcglm} class.
+#'
+#' @param object An object of \code{mcglm} class.
+#' @param ... additional arguments affecting the summary produced. Note
+#'     that there is no extra options for \code{mcglm} object class.
+#'
+#' @return Depending on the number of response variable, the function
+#'     \code{fitted.mcglm} returns a vector (univariate models) or a
+#'     matrix (multivariate models) of fitted values.
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#'
+#' @method fitted mcglm
+#' @export
+
+fitted.mcglm <- function(object, ...) {
+    n_resp <- length(object$beta_names)
+    output <- Matrix(object$fitted, ncol = n_resp, nrow = object$n_obs)
+    return(output)
+}
+#' @title Default Multivariate Covariance Generalized Linear models plotting
+#' @name plot.mcglm
+#'
+#' @description takes a fitted \code{mcglm} object and do plots based on
+#'     residuals, influence diagnostic measures and algorithm check.
+#'
+#' @param x a fitted \code{mcglm} object.
+#' @param type Specify which graphical analysis will be performed.
+#'     Options are: \code{"residuals"}, \code{"influence"} and
+#'     \code{"algorithm"}.
+#' @param ... additional arguments affecting the plot produced. Note
+#'     that there is no extra options for mcglm object class.
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#'
+#' @method plot mcglm
+#' @export
+
+plot.mcglm <- function(x, type = "residuals", ...) {
+    object = x
+    n_resp <- length(object$beta_names)
+    if (type == "residuals") {
+        par(mar = c(2.6, 2.5, 0.1, 0.1),
+            mgp = c(1.6, 0.6, 0),
+            mfrow = c(2, n_resp))
+        for (i in 1:n_resp) {
+            res <- residuals(object, type = "pearson")[, i]
+            fit_values <- fitted(object)[, i]
+            plot(res ~ fit_values,
+                 ylab = "Pearson residuals",
+                 xlab = "Fitted values")
+            temp <- loess.smooth(fitted(object)[, i],
+                                 residuals(object, type = "pearson")[, i])
+            lines(temp$x, temp$y)
+            qqnorm(res)
+            qqline(res)
+        }
+    }
+    if (type == "algorithm") {
+        n_iter <- length(na.exclude(object$IterationCovariance[, 1]))
+        par(mar = c(2.6, 2.5, 0.1, 0.1),
+            mgp = c(1.6, 0.6, 0),
+            mfrow = c(2, 2))
+        matplot(object$IterationRegression[1:c(n_iter + 5), ],
+                type = "l", lty = 2,
+                ylab = "Regression", xlab = "Iterations")
+        matplot(object$IterationCovariance[1:c(n_iter + 5), ],
+                type = "l", lty = 2,
+                ylab = "Covariance", xlab = "Iterations")
+        matplot(object$ScoreRegression[1:c(n_iter + 5), ],
+                type = "l", lty = 2,
+                ylab = "Quasi-score Regression", xlab = "Iterations")
+        matplot(object$ScoreCovariance[1:c(n_iter + 5), ],
+                type = "l", lty = 2,
+                ylab = "Quasi-score Covariance", xlab = "Iterations")
+    }
+    if (type == "partial_residuals") {
+        list_beta <- mc_updateBeta(list_initial = object$list_initial,
+                                   betas = object$Regression,
+                                   n_resp = n_resp,
+                                   information = object$Information)
+        comp_X <- list()
+        for (i in 1:n_resp) {
+            comp_X[[i]] <- as.matrix(object$list_X[[i]]) *
+                as.numeric(list_beta$regression[[i]])
+        }
+        for (i in 1:n_resp) {
+            res <- residuals(object, type = "pearson")[, i]
+            dev.new()
+            n_cov <- dim(comp_X[[i]])[2]
+            par(mar = c(2.6, 2.5, 0.5, 0.5),
+                mgp = c(1.6, 0.6, 0),
+                mfrow = c(1, c(n_cov - 1)))
+            for (j in 2:n_cov) {
+                p1 <- comp_X[[i]][, j] + res
+                plot(p1 ~ object$list_X[[i]][, j],
+                     xlab = object$beta_names[[i]][j],
+                     ylab = "Partial residuals ")
+            }
+        }
+    }
+}
+
+#' @title Print method for Multivariate Covariance Generalized Linear
+#'     Model
+#' @name print.mcglm
+#'
+#' @description The default print method for a \code{mcglm} object.
+#'
+#' @param x fitted model objects of class mcglm as produced by mcglm().
+#' @param ... further arguments passed to or from other methods.
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#'
+#' @rdname print.mcglm
+#' @method print mcglm
+#' @export
+
+print.mcglm <- function(x, ...) {
+    object <- x
+    n_resp <- length(object$beta_names)
+    regression <- mc_updateBeta(list_initial = list(),
+                                betas = object$Regression,
+                                information = object$Information,
+                                n_resp = n_resp)
+    for (i in 1:n_resp) {
+        cat("Call: ")
+        print(object$linear_pred[[i]])
+        cat("\n")
+        cat("Link function:", object$link[[i]])
+        cat("\n")
+        cat("Variance function:", object$variance[[i]])
+        cat("\n")
+        cat("Covariance function:", object$covariance[[i]])
+        cat("\n")
+        names(regression[[1]][[i]]) <- object$beta_names[[i]]
+        cat("Regression:\n")
+        print(regression[[1]][[i]])
+        cat("\n")
+        cat("Dispersion:\n")
+        tau_temp <- coef(object, response = i, type = "tau")$Estimate
+        names(tau_temp) <- rep("", length(tau_temp))
+        print(tau_temp)
+        cat("\n")
+        power_temp <- coef(object, response = i, type = "power")$Estimate
+        if (length(power_temp) != 0) {
+            names(power_temp) <- ""
+            cat("Power:\n")
+            print(power_temp)
+            cat("\n")
+        }
+    }
+}
+
+#' @title Residuals for Multivariate Covariance Generalized Linear
+#'     Models (McGLM)
+#' @name residuals.mcglm
+#'
+#' @description Compute residuals based on fitting \code{mcglm} models.
+#'
+#' @param object An of class \code{mcglm}, typically the result of a
+#'     call to \code{mcglm}.
+#' @param type the type of residuals which should be returned. The
+#'     alternatives are: \code{"raw"} (default), \code{"pearson"} and
+#'     \code{"standardized"}.
+#' @param ... additional arguments affecting the residuals
+#'     produced. Note that there is no extra options for mcglm object
+#'     class.
+#'
+#' @return Depending on the number of response variable the function
+#'     \code{residuals.mcglm} returns a vector (univariate models) or a
+#'     matrix (multivariate models) of residuals values.
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#'
+#' @method residuals mcglm
+#' @export
+
+residuals.mcglm <- function(object, type = "raw", ...) {
+    n_resp <- length(object$beta_names)
+    output <- Matrix(object$residuals,
+                     ncol = n_resp, nrow = object$n_obs)
+    if (type == "standardized") {
+        output <- Matrix(
+            as.numeric(object$residuals %*% chol(object$inv_C)),
+            ncol = n_resp, nrow = object$n_obs)
+    }
+    if (type == "pearson") {
+        output <- Matrix(
+            as.numeric(object$residuals/sqrt(diag(object$C))),
+            ncol = n_resp, nrow = object$n_obs)
+    }
+    return(output)
+}
+
+#' @title Summarizing Multivariate Covariance Generalized Linear Models
+#'     fits.
+#' @name summary.mcglm
+#'
+#' @description Summary for McGLMs objects.
+#'
+#' @param object an object of class \code{mcglm}, usually, a result of a
+#'     call to \code{mcglm}.
+#' @param ... additional arguments affecting the summary produced. Note
+#'     the there is no extra options for mcglm object class.
+#'
+#' @return Print an \code{mcglm} object.
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#'
+#' @method summary mcglm
+#' @export
+
+summary.mcglm <- function(object, ...) {
+    n_resp <- length(object$beta_names)
+    output <- list()
+    for (i in 1:n_resp) {
+        cat("Call: ")
+        print(object$linear_pred[[i]])
+        cat("\n")
+        cat("Link function:", object$link[[i]])
+        cat("\n")
+        cat("Variance function:", object$variance[[i]])
+        cat("\n")
+        cat("Covariance function:", object$covariance[[i]])
+        cat("\n")
+        cat("Regression:\n")
+        tab_beta <- coef(object, std.error = TRUE, response = i, type = "beta")[, 1:2]
+        tab_beta$"Z value" <- tab_beta[, 1]/tab_beta[, 2]
+        rownames(tab_beta) <- object$beta_names[[i]]
+        output[i][[1]]$Regression <- tab_beta
+        print(tab_beta)
+        cat("\n")
+        tab_power <- coef(object, std.error = TRUE, response = i, type = "power")[, 1:2]
+        tab_power$"Z value" <- tab_power[, 1]/tab_power[, 2]
+        rownames(tab_power) <- NULL
+        if (dim(tab_power)[1] != 0) {
+            cat("Power:\n")
+            print(tab_power)
+            output[i][[1]]$Power <- tab_power
+            cat("\n")
+        }
+        cat("Dispersion:\n")
+        tab_tau <- coef(object, std.error = TRUE, response = i, type = "tau")[, 1:2]
+        tab_tau$"Z value" <- tab_tau[, 1]/tab_tau[, 2]
+        rownames(tab_tau) <- NULL
+        output[i][[1]]$tau <- tab_tau
+        print(tab_tau)
+        cat("\n")
+    }
+    tab_rho <- coef(object, std.error = TRUE, response = NA, type = "correlation")[, c(3, 1, 2)]
+    tab_rho$"Z value" <- tab_rho[, 2]/tab_rho[, 3]
+    if (dim(tab_rho)[1] != 0) {
+        cat("Correlation matrix:\n")
+        print(tab_rho)
+        cat("\n")
+    }
+    names(object$con$correct) <- ""
+    iteration_cov <- length(na.exclude(object$IterationCovariance[, 1]))
+    names(iteration_cov) <- ""
+    names(object$con$method) <- ""
+    cat("Algorithm:", object$con$method)
+    cat("\n")
+    cat("Correction:", object$con$correct)
+    cat("\n")
+    cat("Number iterations:", iteration_cov)
+}
+
+#' @title Calculate Variance-Covariance matrix for a fitted McGLM
+#'     object.
+#' @name vcov.mcglm
+#'
+#' @description Returns the variance-covariance matrix for all
+#'     parameters of a \code{mcglm} fitted model object.
+#'
+#' @param object a fitted model \code{mcglm} object.
+#' @param ... additional arguments affecting the summary produced. Note
+#'     that there is no extra options for mcglm object class.
+#'
+#' @return A variance-covariance matrix.
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#'
+#' @method vcov mcglm
+#' @export
+
+vcov.mcglm <- function(object, ...) {
+    cod <- coef(object)$Parameters
+    colnames(object$vcov) <- cod
+    rownames(object$vcov) <- cod
+    return(object$vcov)
+}