diff --git a/R/mc_S3_methods.R b/R/mc_S3_methods.R
new file mode 100644
index 0000000000000000000000000000000000000000..2219d5f6116c91890b22806bf4461ed691ec0b40
--- /dev/null
+++ b/R/mc_S3_methods.R
@@ -0,0 +1,512 @@
+#' @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)
+}