diff --git a/.Rbuildignore b/.Rbuildignore
index 8c50e02d2f63c91bc439832d2594aecc2cdde67f..c7b544af4620fe2f6b43e12da6dc0d6ba159b2f3 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -6,3 +6,5 @@
 .Rd2pdf5516
 .Rd2pdf* 
 buildPkg.R
+README.Rmd
+Examples/
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)
+}
diff --git a/R/mc_anova.R b/R/mc_anova.R
deleted file mode 100644
index f8c2e8377dc5a98128d531ea67a1f0ebb8e56589..0000000000000000000000000000000000000000
--- a/R/mc_anova.R
+++ /dev/null
@@ -1,52 +0,0 @@
-#' ANOVA method for McGLMs.
-#'
-#' @description ANOVA method for McGLMS.
-#'
-#' @param object an object of class 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.
-#' @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)
-}
diff --git a/R/mc_coef.R b/R/mc_coef.R
deleted file mode 100644
index b118f31634937eb42b5468ff4f419c4e1026f252..0000000000000000000000000000000000000000
--- a/R/mc_coef.R
+++ /dev/null
@@ -1,74 +0,0 @@
-#' Extract model coefficients for mcglm class
-#'
-#' coef.mcglm is a function which extracts model coefficients from objects of mcglm class.
-#' @param object An object of mcglm class.
-#' @param std.error Logical. Returns or not the standard errors.
-#' @param response A numeric or vector specyfing for which response variables the coefficients
-#' should be returned.
-#' @param type A string or string vector specyfing which coefficients should be returned.
-#' Options are 'beta', 'tau', 'power', 'tau' and 'correlation'.
-#' @param ... additional arguments affecting the summary produced. Note that there is no extra options for
-#' mcglm object class.
-#' @return A data.frame with estimates, parameters names, response number and parameters type.
-#' @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]] <- paste(paste("beta", i, sep = ""), 0:c(object$Information$n_betas[[i]] - 1), sep = "")
-        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]] <- paste(paste("power", i, sep = ""), 1:object$Information$n_power[[i]], sep = "")
-            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]] <- paste(paste("tau", i, sep = ""), 1:object$Information$n_tau[[i]], sep = "")
-        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] <- paste(paste("rho", combination[1, i], sep = ""), combination[2, i], sep = "")
-        }
-    }
-    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)
-}
diff --git a/R/mc_confint.mcglm.R b/R/mc_confint.mcglm.R
deleted file mode 100644
index 2542e214368128485d041c3e2d6cee5977e4f3a7..0000000000000000000000000000000000000000
--- a/R/mc_confint.mcglm.R
+++ /dev/null
@@ -1,26 +0,0 @@
-#' Confidence Intervals for mcglm
-#'
-#' @description Computes confidence intervals for parameters in a fitted mcglm.
-#'
-#' @param object a fitted mcglm object.
-#' @param parm a specification of which parameters are to be given confidence intervals, either a vector of
-#' number or a vector of names. If missing, all parameters are considered.
-#' @param level the confidence level required.
-#' @param ... additional arguments affecting the confidence interval produced.
-#' Note that there is no extra options for
-#' mcglm object class.
-#' @return A data.frame with confidence intervals, parameters names,
-#' response number and parameters type.
-#' @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) <- paste(format(a, 2), "%", sep = "")
-    rownames(ci) <- temp$Parameters
-    return(ci[parm])
-}
diff --git a/R/mc_fitted.mcglm.R b/R/mc_fitted.mcglm.R
deleted file mode 100644
index ac02339009eae48df18db120a9e4401dd81bf1d7..0000000000000000000000000000000000000000
--- a/R/mc_fitted.mcglm.R
+++ /dev/null
@@ -1,16 +0,0 @@
-#' Extract Model Fitted Values of McGLM
-#'
-#' @description Extract fitted values for objects of mcglm class.
-#'
-#' @param object An object of mcglm class.
-#' @param ... additional arguments affecting the summary produced. Note that there is no extra options for
-#' 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.
-#' @export
-
-fitted.mcglm <- function(object, ...) {
-    n_resp <- length(object$beta_names)
-    output <- Matrix(object$fitted, ncol = n_resp, nrow = object$n_obs)
-    return(output)
-}
diff --git a/R/mc_plot.mcglm.R b/R/mc_plot.mcglm.R
deleted file mode 100644
index 53dd62f1139eaf558f443d1949f2b6f9b1d3363f..0000000000000000000000000000000000000000
--- a/R/mc_plot.mcglm.R
+++ /dev/null
@@ -1,53 +0,0 @@
-#' Default Multivariate Covariance Generalized Linear models plotting
-#'
-#' @description takes a fitted mcglm object by \code{mcglm()} and plots, residuals,
-#' influence diagnostic measures and algorithm check.
-#'
-#' @param x a fitted mcglm object as produced by \code{mcglm()}.
-#' @param type Specify which graphical analysis will be performed, options are: residuals, influence
-#' and algorithm.
-#' @param ... additional arguments affecting the plot produced. Note that there is no extra options for
-#' mcglm object class.
-#' @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 ")
-            }
-        }
-    }
-}
diff --git a/R/mc_print.mcglm.R b/R/mc_print.mcglm.R
deleted file mode 100644
index 5e54ccfb0eea8932667f4c473c23a069824145fe..0000000000000000000000000000000000000000
--- a/R/mc_print.mcglm.R
+++ /dev/null
@@ -1,41 +0,0 @@
-#' Print method for Multivariate Covariance Generalized Linear Model
-#'
-#' @description The default print method for a mcglm object.
-#' @method print mcglm
-#' @param x fitted model objects of class mcglm as produced by mcglm().
-#' @param ... further arguments passed to or from other methods.
-#' @export
-#' @rdname print.mcglm
-
-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")
-        }
-    }
-}
diff --git a/R/mc_qic.R b/R/mc_qic.R
index b20d58bca39a6192808f0ea458a10077fee5ae8e..e51bd66ef74c6313cf7f620dc18922ad4b39cc31 100644
--- a/R/mc_qic.R
+++ b/R/mc_qic.R
@@ -1,10 +1,17 @@
-#' Compute Quasi Information Criterion (QIC) for McGLMs.
+#' @title Compute Quasi Information Criterion (QIC) for McGLMs.
+#' @name qic.mcglm
+#'
+#' @description \code{qic.mcglm} is a function which computes the QIC
+#'     for McGLMs.
+#'
+#' @param object An object of \code{mcglm} class.
+#' @param object.iid An object of \code{mcglm} class contained the model
+#'     fitted using independent covariance structure.
 #'
-#' qic.mcglm is a function which computes the QIC for McGLMs.
-#' @param object An object of mcglm class.
-#' @param object.iid An object of mcglm class contained the model fitted using independent
-#' covariance structure.
 #' @return The QIC value.
+#'
+#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
+#'
 #' @export
 
 qic.mcglm <- function(object, object.iid) {
@@ -19,7 +26,10 @@ qic.mcglm <- function(object, object.iid) {
         if (object$power_fixed[[i]] == TRUE) {
             power = object$list_initial$power[[i]]
         }
-        Q[, i] <- mc_qll(y = obs[, i], mu = mu[, i], variance = object$variance[[i]], power = power)
+        Q[, i] <- mc_qll(y = obs[, i],
+                         mu = mu[, i],
+                         variance = object$variance[[i]],
+                         power = power)
     }
     Vbeta <- -object$inv_S_beta
     Vnull <- solve(-object.iid$inv_S_beta)
diff --git a/R/mc_residuals.mcglm.R b/R/mc_residuals.mcglm.R
deleted file mode 100644
index 0732c56b41d0e03c29798d3d274a367479815d43..0000000000000000000000000000000000000000
--- a/R/mc_residuals.mcglm.R
+++ /dev/null
@@ -1,24 +0,0 @@
-#' Residuals for Multivariate Covariance Generalized Linear Models (McGLM)
-#'
-#' @description Compute residuals based on fitting mcglm models.
-#'
-#' @param object An of class mcglm, typically the result of a call to \code{mcglm}.
-#' @param type the type of residuals which should be returned. The alternatives are: 'raw'
-#' (default), 'pearson' and '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.
-#' @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)
-}
diff --git a/R/mc_summary.mcglm.R b/R/mc_summary.mcglm.R
deleted file mode 100644
index 015079da29c1596e8d03d0a9cbbd26ab817b3e4f..0000000000000000000000000000000000000000
--- a/R/mc_summary.mcglm.R
+++ /dev/null
@@ -1,64 +0,0 @@
-#' Summarizing Multivariate Covariance Generalized Linear Models fits.
-#'
-#' @description Summary for McGLMs objects.
-#'
-#' @param object an object of class 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 mcglm object.
-#' @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)
-}
diff --git a/R/mc_vcov.R b/R/mc_vcov.R
deleted file mode 100644
index 9448b399a6e907fd159dae9e41294d4dc33761e9..0000000000000000000000000000000000000000
--- a/R/mc_vcov.R
+++ /dev/null
@@ -1,16 +0,0 @@
-#' Calculate Variance-Covariance matrix for a fitted McGLM object.
-#'
-#' @description Returns the variance-covariance matrix for all parameters of a mcglm fitted model object.
-#'
-#' @param object a fitted model 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.
-#' @export
-
-vcov.mcglm <- function(object, ...) {
-    cod <- coef(object)$Parameters
-    colnames(object$vcov) <- cod
-    rownames(object$vcov) <- cod
-    return(object$vcov)
-}
diff --git a/R/mcglm.R b/R/mcglm.R
index 8497592ed3ef67143e16cc7340eb7a8ccac2c285..84a11369430944b88a905ff6170a00760ed89b36 100644
--- a/R/mcglm.R
+++ b/R/mcglm.R
@@ -1,56 +1,92 @@
 ##' @title Multivariate covariance generalized linear models (McGLMs)
 ##'
-##'
-##' @description Fits a multivariate covariance generalized linear models (McGLMs) to data.
-##' McGLM is a general framework for non-normal multivariate data analysis, designed to handle
-##' multivariate response variables, along with a wide range of temporal and spatial correlation
-##' structures defined in terms of a covariance link function combined with a matrix linear predictor
-##' involving known matrices. The models take non-normality into account in the conventional way by means
-##' of a variance function, and the mean structure is modelled by means of a link function and a linear predictor.
-##' The models are fitted using an efficient Newton scoring algorithm based on quasi-likelihood and
-##' Pearson estimating functions, using only second-moment assumptions.
-##' This provides a unified approach to a wide variety of different types of response variables and
-##' covariance structures, including multivariate extensions of repeated measures, time series, longitudinal,
-##' spatial and spatio-temporal structures.
-##'
+##' @description Fits a multivariate covariance generalized linear
+##'     models (McGLMs) to data.  McGLM is a general framework for
+##'     non-normal multivariate data analysis, designed to handle
+##'     multivariate response variables, along with a wide range of
+##'     temporal and spatial correlation structures defined in terms of
+##'     a covariance link function combined with a matrix linear
+##'     predictor involving known matrices. The models take
+##'     non-normality into account in the conventional way by means of a
+##'     variance function, and the mean structure is modelled by means
+##'     of a link function and a linear predictor.  The models are
+##'     fitted using an efficient Newton scoring algorithm based on
+##'     quasi-likelihood and Pearson estimating functions, using only
+##'     second-moment assumptions.  This provides a unified approach to
+##'     a wide variety of different types of response variables and
+##'     covariance structures, including multivariate extensions of
+##'     repeated measures, time series, longitudinal, spatial and
+##'     spatio-temporal structures.
 ##'
 ##' @docType package
 ##' @name mcglm
 NULL
 
-#' @name ahs
-#'
 #' @title Australian health survey
+#' @name ahs
 #'
-#' @description The Australian health survey was used by Bonat and Jorgensen (2015) as an example of multivariate
-#' count regression model. The data consists of five count response variables concerning health system access
-#' measures and nine covariates concerning  social conditions in Australian for 1987-88.
+#' @description The Australian health survey was used by Bonat and
+#'     Jorgensen (2015) as an example of multivariate count regression
+#'     model. The data consists of five count response variables
+#'     concerning health system access measures and nine covariates
+#'     concerning social conditions in Australian for 1987-88.
 #'
 #' \itemize{
-#'     \item \code{sex} - Factor, two levels (0-Male; 1-Female).
-#'     \item \code{age} - Respondent's age in years divided by 100.
-#'     \item \code{income} - Respondent's annual income in Australian dollars divided by 1000.
-#'     \item \code{levyplus} - Factor, two levels (1- if respondent is covered by private health
-#'     insurance fund for private patients in public hospital (with doctor of choice); 0 - otherwise).
-#'     \item \code{freepoor} - Factor, two levels (1 - if respondent is covered by government because low income,
-#'     recent immigrant, unemployed; 0 - otherwise).
-#'     \item \code{freerepa} - Factor, two levels (1 - if respondent is covered free by government because of
-#'     old-age or disability pension, or because invalid veteran or family of deceased veteran; 0 - otherwise).
-#'     \item \code{illnes} - Number of illnesses in past 2 weeks, with 5 or more weeks coded as 5.
-#'     \item \code{actdays} - Number of days of reduced activity in the past two weeks due to illness or injury.
-#'     \item \code{hscore} - Respondent's general health questionnaire score using Goldberg's method;
-#'     high score indicates poor health.
-#'     \item \code{chcond1} - Factor, two levels (1 - if respondent has chronic condition(s) but is not limited
-#'     in activity; 0 - otherwise).
-#'     \item \code{chcond2} - Factor, two levels (1 if respondent has chronic condition(s) and is limited in
-#'     activity; 0 - otherwise).
-#'     \item \code{Ndoc} - Number of consultations with a doctor or specialist (response variable).
-#'     \item \code{Nndoc} - Number of consultations with health professionals (response variable).
-#'     \item \code{Nadm} - Number of admissions to a hospital, psychiatric hospital, nursing or
-#'     convalescence home in the past 12 months (response variable).
-#'     \item \code{Nhosp} - Number of nights in a hospital during the most recent admission.
-#'     \item \code{Nmed} - Total number of prescribed and non prescribed medications used in the past two days.
-#'     \item \code{id} - Respondent's index.
+#'
+#' \item \code{sex} - Factor, two levels (0-Male; 1-Female).
+#'
+#' \item \code{age} - Respondent's age in years divided by 100.
+#'
+#' \item \code{income} - Respondent's annual income in Australian
+#'     dollars divided by 1000.
+#'
+#' \item \code{levyplus} - Factor, two levels (1- if respondent is
+#'     covered by private health insurance fund for private patients in
+#'     public hospital (with doctor of choice); 0 - otherwise).
+#'
+#' \item \code{freepoor} - Factor, two levels (1 - if respondent is
+#'     covered by government because low income, recent immigrant,
+#'     unemployed; 0 - otherwise).
+#'
+#' \item \code{freerepa} - Factor, two levels (1 - if respondent is
+#'     covered free by government because of old-age or disability
+#'     pension, or because invalid veteran or family of deceased
+#'     veteran; 0 - otherwise).
+#'
+#' \item \code{illnes} - Number of illnesses in past 2 weeks, with 5 or
+#'     more weeks coded as 5.
+#'
+#' \item \code{actdays} - Number of days of reduced activity in the past
+#'     two weeks due to illness or injury.
+#'
+#' \item \code{hscore} - Respondent's general health questionnaire score
+#'     using Goldberg's method; high score indicates poor health.
+#'
+#' \item \code{chcond1} - Factor, two levels (1 - if respondent has
+#'     chronic condition(s) but is not limited in activity; 0 -
+#'     otherwise).
+#'
+#' \item \code{chcond2} - Factor, two levels (1 if respondent has
+#'     chronic condition(s) and is limited in activity; 0 - otherwise).
+#'
+#' \item \code{Ndoc} - Number of consultations with a doctor or
+#'     specialist (response variable).
+#'
+#' \item \code{Nndoc} - Number of consultations with health
+#'     professionals (response variable).
+#'
+#' \item \code{Nadm} - Number of admissions to a hospital, psychiatric
+#'     hospital, nursing or convalescence home in the past 12 months
+#'     (response variable).
+#'
+#' \item \code{Nhosp} - Number of nights in a hospital during the most
+#'     recent admission.
+#'
+#' \item \code{Nmed} - Total number of prescribed and non prescribed
+#'     medications used in the past two days.
+#'
+#' \item \code{id} - Respondent's index.
+#'
 #' }
 #'
 #' @docType data
@@ -61,7 +97,8 @@ NULL
 #'
 #' @format a \code{data.frame} with 5190 records and 17 variables.
 #'
-#' @source Deb, P. and Trivedi, P. K. (1997). Demand for medical care by the elderly: A finite mixture approach,
-#' Journal of Applied Econometrics 12(3):313--336.
+#' @source Deb, P. and Trivedi, P. K. (1997). Demand for medical care by
+#'     the elderly: A finite mixture approach, Journal of Applied
+#'     Econometrics 12(3):313--336.
 #'
 NULL
diff --git a/man/ahs.Rd b/man/ahs.Rd
index 4569dd4eca3aaa17f45a2a28db2ebe070e5d7685..7cb32800493c4d577077f654dad50863d8004854 100644
--- a/man/ahs.Rd
+++ b/man/ahs.Rd
@@ -6,42 +6,76 @@
 \title{Australian health survey}
 \format{a \code{data.frame} with 5190 records and 17 variables.}
 \source{
-Deb, P. and Trivedi, P. K. (1997). Demand for medical care by the elderly: A finite mixture approach,
-Journal of Applied Econometrics 12(3):313--336.
+Deb, P. and Trivedi, P. K. (1997). Demand for medical care by
+    the elderly: A finite mixture approach, Journal of Applied
+    Econometrics 12(3):313--336.
 }
 \usage{
 data(ahs)
 }
 \description{
-The Australian health survey was used by Bonat and Jorgensen (2015) as an example of multivariate
-count regression model. The data consists of five count response variables concerning health system access
-measures and nine covariates concerning  social conditions in Australian for 1987-88.
+The Australian health survey was used by Bonat and
+    Jorgensen (2015) as an example of multivariate count regression
+    model. The data consists of five count response variables
+    concerning health system access measures and nine covariates
+    concerning social conditions in Australian for 1987-88.
 
 \itemize{
-    \item \code{sex} - Factor, two levels (0-Male; 1-Female).
-    \item \code{age} - Respondent's age in years divided by 100.
-    \item \code{income} - Respondent's annual income in Australian dollars divided by 1000.
-    \item \code{levyplus} - Factor, two levels (1- if respondent is covered by private health
-    insurance fund for private patients in public hospital (with doctor of choice); 0 - otherwise).
-    \item \code{freepoor} - Factor, two levels (1 - if respondent is covered by government because low income,
-    recent immigrant, unemployed; 0 - otherwise).
-    \item \code{freerepa} - Factor, two levels (1 - if respondent is covered free by government because of
-    old-age or disability pension, or because invalid veteran or family of deceased veteran; 0 - otherwise).
-    \item \code{illnes} - Number of illnesses in past 2 weeks, with 5 or more weeks coded as 5.
-    \item \code{actdays} - Number of days of reduced activity in the past two weeks due to illness or injury.
-    \item \code{hscore} - Respondent's general health questionnaire score using Goldberg's method;
-    high score indicates poor health.
-    \item \code{chcond1} - Factor, two levels (1 - if respondent has chronic condition(s) but is not limited
-    in activity; 0 - otherwise).
-    \item \code{chcond2} - Factor, two levels (1 if respondent has chronic condition(s) and is limited in
-    activity; 0 - otherwise).
-    \item \code{Ndoc} - Number of consultations with a doctor or specialist (response variable).
-    \item \code{Nndoc} - Number of consultations with health professionals (response variable).
-    \item \code{Nadm} - Number of admissions to a hospital, psychiatric hospital, nursing or
-    convalescence home in the past 12 months (response variable).
-    \item \code{Nhosp} - Number of nights in a hospital during the most recent admission.
-    \item \code{Nmed} - Total number of prescribed and non prescribed medications used in the past two days.
-    \item \code{id} - Respondent's index.
+
+\item \code{sex} - Factor, two levels (0-Male; 1-Female).
+
+\item \code{age} - Respondent's age in years divided by 100.
+
+\item \code{income} - Respondent's annual income in Australian
+    dollars divided by 1000.
+
+\item \code{levyplus} - Factor, two levels (1- if respondent is
+    covered by private health insurance fund for private patients in
+    public hospital (with doctor of choice); 0 - otherwise).
+
+\item \code{freepoor} - Factor, two levels (1 - if respondent is
+    covered by government because low income, recent immigrant,
+    unemployed; 0 - otherwise).
+
+\item \code{freerepa} - Factor, two levels (1 - if respondent is
+    covered free by government because of old-age or disability
+    pension, or because invalid veteran or family of deceased
+    veteran; 0 - otherwise).
+
+\item \code{illnes} - Number of illnesses in past 2 weeks, with 5 or
+    more weeks coded as 5.
+
+\item \code{actdays} - Number of days of reduced activity in the past
+    two weeks due to illness or injury.
+
+\item \code{hscore} - Respondent's general health questionnaire score
+    using Goldberg's method; high score indicates poor health.
+
+\item \code{chcond1} - Factor, two levels (1 - if respondent has
+    chronic condition(s) but is not limited in activity; 0 -
+    otherwise).
+
+\item \code{chcond2} - Factor, two levels (1 if respondent has
+    chronic condition(s) and is limited in activity; 0 - otherwise).
+
+\item \code{Ndoc} - Number of consultations with a doctor or
+    specialist (response variable).
+
+\item \code{Nndoc} - Number of consultations with health
+    professionals (response variable).
+
+\item \code{Nadm} - Number of admissions to a hospital, psychiatric
+    hospital, nursing or convalescence home in the past 12 months
+    (response variable).
+
+\item \code{Nhosp} - Number of nights in a hospital during the most
+    recent admission.
+
+\item \code{Nmed} - Total number of prescribed and non prescribed
+    medications used in the past two days.
+
+\item \code{id} - Respondent's index.
+
 }
 }
 \keyword{datasets}
diff --git a/man/anova.mcglm.Rd b/man/anova.mcglm.Rd
index d367ebbf685bcb454cb9bf677f598f741b2458d2..5f4e4d9b241f5553306a2cbfa7728b513db9c0ee 100644
--- a/man/anova.mcglm.Rd
+++ b/man/anova.mcglm.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_anova.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{anova.mcglm}
 \alias{anova.mcglm}
 \title{ANOVA method for McGLMs.}
@@ -7,12 +7,22 @@
 \method{anova}{mcglm}(object, ...)
 }
 \arguments{
-\item{object}{an object of class mcglm, usually, a result of a call to \code{mcglm}.}
+\item{object}{an object of class \code{mcglm}, usually, a result of a
+call to \code{mcglm()}.}
 
-\item{...}{additional arguments affecting the summary produced. Note that there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the summary produced. Note
+    that there is no extra options for mcglm object class.}
+}
+\value{
+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.
 }
 \description{
-ANOVA method for McGLMS.
+ANOVA method for object of class McGLMS.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/coef.mcglm.Rd b/man/coef.mcglm.Rd
index 741388b2237b81fe9ebe23671621ddc3c8ece256..127d214b5152113798f2404f7aafc658515a13d8 100644
--- a/man/coef.mcglm.Rd
+++ b/man/coef.mcglm.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_coef.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{coef.mcglm}
 \alias{coef.mcglm}
 \title{Extract model coefficients for mcglm class}
@@ -9,23 +9,31 @@
   "correlation"), ...)
 }
 \arguments{
-\item{object}{An object of mcglm class.}
+\item{object}{An object of \code{mcglm} class.}
 
-\item{std.error}{Logical. Returns or not the standard errors.}
+\item{std.error}{Logical. If \code{TRUE} returns the standard errors
+of the estimates. Default is \code{FALSE}.}
 
-\item{response}{A numeric or vector specyfing for which response variables the coefficients
-should be returned.}
+\item{response}{A numeric vector specyfing for which response
+variables the coefficients should be returned.}
 
-\item{type}{A string or string vector specyfing which coefficients should be returned.
-Options are 'beta', 'tau', 'power', 'tau' and 'correlation'.}
+\item{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"}.}
 
-\item{...}{additional arguments affecting the summary produced. Note that there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the summary produced. Note
+    that there is no extra options for mcglm object class.}
 }
 \value{
-A data.frame with estimates, parameters names, response number and parameters type.
+A \code{data.frame} with parameters names, estimates,
+    response number and parameters type.
 }
 \description{
-coef.mcglm is a function which extracts model coefficients from objects of mcglm class.
+\code{coef.mcglm} is a function which extracts model
+    coefficients from objects of \code{mcglm} class.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/confint.mcglm.Rd b/man/confint.mcglm.Rd
index efdb6bb31b69cae0964019622e95df923c4ffd36..677abb619e7afb01b4ed28cffcfec59cb32f5797 100644
--- a/man/confint.mcglm.Rd
+++ b/man/confint.mcglm.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_confint.mcglm.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{confint.mcglm}
 \alias{confint.mcglm}
 \title{Confidence Intervals for mcglm}
@@ -7,22 +7,27 @@
 \method{confint}{mcglm}(object, parm, level = 0.95, ...)
 }
 \arguments{
-\item{object}{a fitted mcglm object.}
+\item{object}{a fitted \code{mcglm} object.}
 
-\item{parm}{a specification of which parameters are to be given confidence intervals, either a vector of
-number or a vector of names. If missing, all parameters are considered.}
+\item{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.}
 
-\item{level}{the confidence level required.}
+\item{level}{the nominal confidence level.}
 
-\item{...}{additional arguments affecting the confidence interval produced.
-Note that there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the confidence interval
+    produced. Note that there is no extra options for \code{mcglm}
+    object class.}
 }
 \value{
-A data.frame with confidence intervals, parameters names,
-response number and parameters type.
+A \code{data.frame} with confidence intervals, parameters
+    names, response number and parameters type.
 }
 \description{
-Computes confidence intervals for parameters in a fitted mcglm.
+Computes confidence intervals for parameters in a fitted
+    \code{mcglm} model.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/fitted.mcglm.Rd b/man/fitted.mcglm.Rd
index 115d5810f201bc4391de693565a8e49413d55eb2..64bb7b17a4a68982952da97b25a98e1d863dbbd3 100644
--- a/man/fitted.mcglm.Rd
+++ b/man/fitted.mcglm.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_fitted.mcglm.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{fitted.mcglm}
 \alias{fitted.mcglm}
 \title{Extract Model Fitted Values of McGLM}
@@ -7,16 +7,20 @@
 \method{fitted}{mcglm}(object, ...)
 }
 \arguments{
-\item{object}{An object of mcglm class.}
+\item{object}{An object of \code{mcglm} class.}
 
-\item{...}{additional arguments affecting the summary produced. Note that there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the summary produced. Note
+    that there is no extra options for \code{mcglm} object class.}
 }
 \value{
-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.
+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.
 }
 \description{
-Extract fitted values for objects of mcglm class.
+Extract fitted values for objects of \code{mcglm} class.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/mcglm.Rd b/man/mcglm.Rd
index 63c324fce7fbcaf6a8667b8cc84e3b626b68aa91..402d73c5cd3c00fe42f62ed64ae1c8e1006f5fd1 100644
--- a/man/mcglm.Rd
+++ b/man/mcglm.Rd
@@ -52,16 +52,22 @@ fitted using an estimating function approach, combining quasi-score functions fo
 parameters and Pearson estimating function for covariance parameters. For details see Bonat and
 Jorgensen (2015).
 
-Fits a multivariate covariance generalized linear models (McGLMs) to data.
-McGLM is a general framework for non-normal multivariate data analysis, designed to handle
-multivariate response variables, along with a wide range of temporal and spatial correlation
-structures defined in terms of a covariance link function combined with a matrix linear predictor
-involving known matrices. The models take non-normality into account in the conventional way by means
-of a variance function, and the mean structure is modelled by means of a link function and a linear predictor.
-The models are fitted using an efficient Newton scoring algorithm based on quasi-likelihood and
-Pearson estimating functions, using only second-moment assumptions.
-This provides a unified approach to a wide variety of different types of response variables and
-covariance structures, including multivariate extensions of repeated measures, time series, longitudinal,
-spatial and spatio-temporal structures.
+Fits a multivariate covariance generalized linear
+    models (McGLMs) to data.  McGLM is a general framework for
+    non-normal multivariate data analysis, designed to handle
+    multivariate response variables, along with a wide range of
+    temporal and spatial correlation structures defined in terms of
+    a covariance link function combined with a matrix linear
+    predictor involving known matrices. The models take
+    non-normality into account in the conventional way by means of a
+    variance function, and the mean structure is modelled by means
+    of a link function and a linear predictor.  The models are
+    fitted using an efficient Newton scoring algorithm based on
+    quasi-likelihood and Pearson estimating functions, using only
+    second-moment assumptions.  This provides a unified approach to
+    a wide variety of different types of response variables and
+    covariance structures, including multivariate extensions of
+    repeated measures, time series, longitudinal, spatial and
+    spatio-temporal structures.
 }
 
diff --git a/man/plot.mcglm.Rd b/man/plot.mcglm.Rd
index a1b852cc05cd1bff264060e921f0b7c96864cdb6..ebb64e8fd5d3dc830eeba23dd354d56390fe49d7 100644
--- a/man/plot.mcglm.Rd
+++ b/man/plot.mcglm.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_plot.mcglm.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{plot.mcglm}
 \alias{plot.mcglm}
 \title{Default Multivariate Covariance Generalized Linear models plotting}
@@ -7,16 +7,20 @@
 \method{plot}{mcglm}(x, type = "residuals", ...)
 }
 \arguments{
-\item{x}{a fitted mcglm object as produced by \code{mcglm()}.}
+\item{x}{a fitted \code{mcglm} object.}
 
-\item{type}{Specify which graphical analysis will be performed, options are: residuals, influence
-and algorithm.}
+\item{type}{Specify which graphical analysis will be performed.
+Options are: \code{"residuals"}, \code{"influence"} and
+\code{"algorithm"}.}
 
-\item{...}{additional arguments affecting the plot produced. Note that there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the plot produced. Note
+that there is no extra options for mcglm object class.}
 }
 \description{
-takes a fitted mcglm object by \code{mcglm()} and plots, residuals,
-influence diagnostic measures and algorithm check.
+takes a fitted \code{mcglm} object and do plots based on
+    residuals, influence diagnostic measures and algorithm check.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/print.mcglm.Rd b/man/print.mcglm.Rd
index 008141d86c9bb9cac15fb6e7b5427670d0d3eff6..c712e79bcad170d71bc11813a36011e95e1532cc 100644
--- a/man/print.mcglm.Rd
+++ b/man/print.mcglm.Rd
@@ -1,8 +1,9 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_print.mcglm.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{print.mcglm}
 \alias{print.mcglm}
-\title{Print method for Multivariate Covariance Generalized Linear Model}
+\title{Print method for Multivariate Covariance Generalized Linear
+    Model}
 \usage{
 \method{print}{mcglm}(x, ...)
 }
@@ -12,6 +13,9 @@
 \item{...}{further arguments passed to or from other methods.}
 }
 \description{
-The default print method for a mcglm object.
+The default print method for a \code{mcglm} object.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/qic.mcglm.Rd b/man/qic.mcglm.Rd
index a98434607d65cc254cbb68a96873b28c16760361..02235c1afd4b475b93b3badb8ff5f0aed4dd594f 100644
--- a/man/qic.mcglm.Rd
+++ b/man/qic.mcglm.Rd
@@ -7,15 +7,19 @@
 qic.mcglm(object, object.iid)
 }
 \arguments{
-\item{object}{An object of mcglm class.}
+\item{object}{An object of \code{mcglm} class.}
 
-\item{object.iid}{An object of mcglm class contained the model fitted using independent
-covariance structure.}
+\item{object.iid}{An object of \code{mcglm} class contained the model
+    fitted using independent covariance structure.}
 }
 \value{
 The QIC value.
 }
 \description{
-qic.mcglm is a function which computes the QIC for McGLMs.
+\code{qic.mcglm} is a function which computes the QIC
+    for McGLMs.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/residuals.mcglm.Rd b/man/residuals.mcglm.Rd
index 0e5cc27de00e03907f2492c1e420db57674f0fa4..495705b9b5e2b99cff28a0cdb11a806cc4146bcc 100644
--- a/man/residuals.mcglm.Rd
+++ b/man/residuals.mcglm.Rd
@@ -1,25 +1,33 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_residuals.mcglm.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{residuals.mcglm}
 \alias{residuals.mcglm}
-\title{Residuals for Multivariate Covariance Generalized Linear Models (McGLM)}
+\title{Residuals for Multivariate Covariance Generalized Linear
+    Models (McGLM)}
 \usage{
 \method{residuals}{mcglm}(object, type = "raw", ...)
 }
 \arguments{
-\item{object}{An of class mcglm, typically the result of a call to \code{mcglm}.}
+\item{object}{An of class \code{mcglm}, typically the result of a
+call to \code{mcglm}.}
 
-\item{type}{the type of residuals which should be returned. The alternatives are: 'raw'
-(default), 'pearson' and 'standardized'.}
+\item{type}{the type of residuals which should be returned. The
+alternatives are: \code{"raw"} (default), \code{"pearson"} and
+\code{"standardized"}.}
 
-\item{...}{additional arguments affecting the residuals produced. Note that there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the residuals
+    produced. Note that there is no extra options for mcglm object
+    class.}
 }
 \value{
-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.
+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.
 }
 \description{
-Compute residuals based on fitting mcglm models.
+Compute residuals based on fitting \code{mcglm} models.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }
 
diff --git a/man/summary.mcglm.Rd b/man/summary.mcglm.Rd
index 8e94d1de27313318a23c0653e52c4da6ec0cbb7b..8061e7980152f9e306756a91f5d40889625ba597 100644
--- a/man/summary.mcglm.Rd
+++ b/man/summary.mcglm.Rd
@@ -1,21 +1,26 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_summary.mcglm.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{summary.mcglm}
 \alias{summary.mcglm}
-\title{Summarizing Multivariate Covariance Generalized Linear Models fits.}
+\title{Summarizing Multivariate Covariance Generalized Linear Models
+    fits.}
 \usage{
 \method{summary}{mcglm}(object, ...)
 }
 \arguments{
-\item{object}{an object of class mcglm, usually, a result of a call to \code{mcglm}.}
+\item{object}{an object of class \code{mcglm}, usually, a result of a
+call to \code{mcglm}.}
 
-\item{...}{additional arguments affecting the summary produced. Note the there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the summary produced. Note
+    the there is no extra options for mcglm object class.}
 }
 \value{
-Print an mcglm object.
+Print an \code{mcglm} object.
 }
 \description{
 Summary for McGLMs objects.
 }
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
+}
 
diff --git a/man/vcov.mcglm.Rd b/man/vcov.mcglm.Rd
index 31f20be578f2a7d118ded6e55b24a24c585511c4..3ba9ffafda81d6fe7a5452222b20d5e94e83507d 100644
--- a/man/vcov.mcglm.Rd
+++ b/man/vcov.mcglm.Rd
@@ -1,21 +1,26 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mc_vcov.R
+% Please edit documentation in R/mc_S3_methods.R
 \name{vcov.mcglm}
 \alias{vcov.mcglm}
-\title{Calculate Variance-Covariance matrix for a fitted McGLM object.}
+\title{Calculate Variance-Covariance matrix for a fitted McGLM
+    object.}
 \usage{
 \method{vcov}{mcglm}(object, ...)
 }
 \arguments{
-\item{object}{a fitted model mcglm object.}
+\item{object}{a fitted model \code{mcglm} object.}
 
-\item{...}{additional arguments affecting the summary produced. Note that there is no extra options for
-mcglm object class.}
+\item{...}{additional arguments affecting the summary produced. Note
+    that there is no extra options for mcglm object class.}
 }
 \value{
 A variance-covariance matrix.
 }
 \description{
-Returns the variance-covariance matrix for all parameters of a mcglm fitted model object.
+Returns the variance-covariance matrix for all
+    parameters of a \code{mcglm} fitted model object.
+}
+\author{
+Wagner Hugo Bonat, \email{wbonat@ufpr.br}
 }