diff --git a/DESCRIPTION b/DESCRIPTION
index f5028c06ea556626e75b63aca71c77797793aed2..0fb0b4bab2b576c933ab79fc9e5dbc9e96a06061 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -11,12 +11,13 @@ Authors@R: as.person(c(
 Description: This package fit multivariate covariance generalized linear
     models (McGLM).
 Depends:
-    R (>= 3.2.1),
-    assertthat
-Imports:
-    Matrix
+    R (>= 3.2.1)
 Suggests:
-    testthat
+    testthat,
+    plyr
+Imports:
+    Matrix,
+    asserthat
 License: GPL-2
 LazyData: TRUE
 URL: http://git.leg.ufpr.br/wbonat/mcglm
diff --git a/NAMESPACE b/NAMESPACE
index 42e7eb9fa3bce240ae229ae3c526ffb295ad8ded..ac113d4a61f1ee18f1ba84a70e91bf52faedef2f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,2 +1,28 @@
-exportPattern("^[[:alpha:]]+")
-import(Matrix, assertthat)
+# Generated by roxygen2 (4.1.1): do not edit by hand
+
+S3method(anova,mcglm)
+S3method(coef,mcglm)
+S3method(confint,mcglm)
+S3method(fitted,mcglm)
+S3method(plot,mcglm)
+S3method(print,mcglm)
+S3method(residuals,mcglm)
+S3method(summary,mcglm)
+S3method(vcov,mcglm)
+export(fit_mcglm)
+export(mc_bias_corrected_std)
+export(mc_dexp_gold)
+export(mc_influence)
+export(mc_initial_values)
+export(mc_link_function)
+export(mc_matrix_linear_predictor)
+export(mc_qll)
+export(mc_robust_std)
+export(mc_rw1)
+export(mc_rw2)
+export(mc_unstructured)
+export(mc_variance_function)
+export(mcglm)
+export(qic.mcglm)
+import(Matrix)
+import(assertthat)
diff --git a/R/fit_mcglm.R b/R/fit_mcglm.R
index 3c2d3f8669027234964061ccf1458780c06cf013..f2b9a81cd0f4bd35c5c5737c407db2935e565fbc 100644
--- a/R/fit_mcglm.R
+++ b/R/fit_mcglm.R
@@ -28,8 +28,10 @@
 #' @return A list with estimated regression and covariance parameters.
 #' @export
 
-fit_mcglm <- function(list_initial, list_link, list_variance, list_covariance, list_X, list_Z, list_offset, list_Ntrial,
-    list_power_fixed, list_sparse, y_vec, correct = TRUE, max_iter, tol = 0.001, method = "rc", tunning = 0, verbose) {
+fit_mcglm <- function(list_initial, list_link, list_variance, list_covariance, list_X, list_Z,
+                      list_offset, list_Ntrial, list_power_fixed, list_sparse, y_vec,
+                      correct = TRUE, max_iter, tol = 0.001, method = "rc",
+                      tunning = 0, verbose) {
     ## Transformation from list to vector
     parametros <- mc_list2vec(list_initial, list_power_fixed)
     n_resp <- length(list_initial$regression)
diff --git a/R/mc_anova.R b/R/mc_anova.R
index b23a1d9a0eebd7adabc72bb08e3ea134457d24ad..f8c2e8377dc5a98128d531ea67a1f0ebb8e56589 100644
--- a/R/mc_anova.R
+++ b/R/mc_anova.R
@@ -3,9 +3,11 @@
 #' @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) {
+anova.mcglm <- function(object, ...) {
     n_resp <- length(object$mu_list)
     n_beta <- lapply(object$list_X, ncol)
     idx.list <- list()
diff --git a/R/mc_auxiliar.R b/R/mc_auxiliar.R
index ab99821a27761f5da52f77d9627429572a36288f..f59257a0e4290c6aea20ea8c0720d8e96f1d8346 100644
--- a/R/mc_auxiliar.R
+++ b/R/mc_auxiliar.R
@@ -9,12 +9,6 @@
 #' @param bord1 A matrix.
 #' @param bord2 A matrix.
 #' @return The matrix product bord1*middle*bord2.
-#' @examples
-#' X1 <- Diagonal(5, 2)
-#' M <- Matrix(rep(0.8,5)%*%t(rep(0.8,5)))
-#' mc_sandwich(middle = M, bord1 = X1, bord2 = X1)
-#' mc_sandwich_negative(middle = M, bord1 = X1, bord2 = X1)
-#' @export
 
 ## Auxiliar function to multiply matrices
 mc_sandwich <- function(middle, bord1, bord2) {
diff --git a/R/mc_build_C.R b/R/mc_build_C.R
index d436e3d53c7df1e7b7153154191e7d6464ec56bb..3cf758cb092357b4d95542c9d45652f432519687 100644
--- a/R/mc_build_C.R
+++ b/R/mc_build_C.R
@@ -19,10 +19,11 @@
 #'
 #'@return A list with the inverse of the C matrix and the derivatives of the C matrix with respect to
 #'rho, power and tau parameters.
-#'@export
 
-mc_build_C <- function(list_mu, list_Ntrial, rho, list_tau, list_power, list_Z, list_sparse, list_variance,
-                       list_covariance, list_power_fixed, compute_C = FALSE, compute_derivative_beta = FALSE,
+mc_build_C <- function(list_mu, list_Ntrial, rho, list_tau, list_power,
+                       list_Z, list_sparse, list_variance,
+                       list_covariance, list_power_fixed, compute_C = FALSE,
+                       compute_derivative_beta = FALSE,
                        compute_derivative_cov = TRUE) {
     n_resp <- length(list_mu)
     n_obs <- length(list_mu[[1]][[1]])
diff --git a/R/mc_build_bdiag.R b/R/mc_build_bdiag.R
index c41c93e2c8debc869c02fdf504e453cd510ef4fc..7d9a7290a408d52d2ffcd7e8435b5f143365a5e9 100644
--- a/R/mc_build_bdiag.R
+++ b/R/mc_build_bdiag.R
@@ -7,7 +7,6 @@
 #' @param n_obs A numeric specifying the number of observations in the data set.
 #' @return A list of zero matrices.
 #' @details It is an internal function.
-#' @export
 
 mc_build_bdiag <- function(n_resp, n_obs) {
     list_zero <- list()
@@ -15,4 +14,4 @@ mc_build_bdiag <- function(n_resp, n_obs) {
         list_zero[[i]] <- Matrix(0, n_obs, n_obs, sparse = TRUE)
     }
     return(list_zero)
-} 
+}
diff --git a/R/mc_build_omega.R b/R/mc_build_omega.R
index 80b1b3a57ce4d09ac519dc4c650ccffd053836c2..73ad728a4b0532668b166d293ebe1fa1bb21a946 100644
--- a/R/mc_build_omega.R
+++ b/R/mc_build_omega.R
@@ -7,7 +7,7 @@
 #' @param covariance_link String specifing the covariance link function: identity, inverse, expm.
 #' @param sparse Logical force to use sparse matrix representation 'dsCMatrix'.
 #' @return A list with the \eqn{\Omega} matrix its inverse and derivatives with respect to \eqn{\tau}.
-#' @export
+
 mc_build_omega <- function(tau, Z, covariance_link, sparse = FALSE) {
     if (covariance_link == "identity") {
         Omega <- mc_matrix_linear_predictor(tau = tau, Z = Z)
@@ -24,4 +24,4 @@ mc_build_omega <- function(tau, Z, covariance_link, sparse = FALSE) {
         output <- list(inv_Omega = inv_Omega, D_inv_Omega = Z)
     }
     return(output)
-} 
+}
diff --git a/R/mc_build_sigma.R b/R/mc_build_sigma.R
index bdf47bad24757d81d684bde6c7a420087c19939c..d46aa5e5b8fe853f399c2fd393941ed2eb5c7cde 100644
--- a/R/mc_build_sigma.R
+++ b/R/mc_build_sigma.R
@@ -22,10 +22,9 @@
 #'of \eqn{\Sigma} with respect to the power and tau parameters.
 #'@seealso \code{\link{mc_link_function}}, \code{\link{mc_variance_function}},
 #'\code{\link{mc_build_omega}}.
-#'@export
 
-mc_build_sigma <- function(mu, Ntrial = 1, tau, power, Z, sparse, variance, covariance, power_fixed,
-                           compute_derivative_beta = FALSE) {
+mc_build_sigma <- function(mu, Ntrial = 1, tau, power, Z, sparse, variance,
+                           covariance, power_fixed, compute_derivative_beta = FALSE) {
     if (variance == "constant") {
         if (covariance == "identity" | covariance == "expm") {
             Omega <- mc_build_omega(tau = tau, Z = Z, covariance_link = covariance, sparse = sparse)
diff --git a/R/mc_build_sigmab.R b/R/mc_build_sigmab.R
index 446df48876ad7bbfde490ffebb5dc363bd848037..18c0389e1d2c833657d09d0ac1583a465143e146 100644
--- a/R/mc_build_sigmab.R
+++ b/R/mc_build_sigmab.R
@@ -6,7 +6,6 @@
 #' @param n_resp A numeric.
 #' @param inverse Logical
 #' @return A list with sigmab and its derivatives with respect to rho.
-#' @export
 
 mc_build_sigma_between <- function(rho, n_resp, inverse = FALSE) {
     output <- list(Sigmab = 1, D_Sigmab = 1)
@@ -38,4 +37,4 @@ mc_derivative_sigma_between <- function(n_resp) {
         list.Derivative[i][[1]] <- Derivative
     }
     return(list.Derivative)
-} 
+}
diff --git a/R/mc_coef.R b/R/mc_coef.R
index 6f7a5b47f955772abef852c4503420a348e4f4aa..b118f31634937eb42b5468ff4f419c4e1026f252 100644
--- a/R/mc_coef.R
+++ b/R/mc_coef.R
@@ -7,10 +7,13 @@
 #' 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")) {
+
+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()
diff --git a/R/mc_confint.mcglm.R b/R/mc_confint.mcglm.R
index 457b52afc19c712eba269d617e6ba7424bdd27ff..2542e214368128485d041c3e2d6cee5977e4f3a7 100644
--- a/R/mc_confint.mcglm.R
+++ b/R/mc_confint.mcglm.R
@@ -3,17 +3,24 @@
 #' @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, level = 0.95) {
+
+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)
-} 
+    return(ci[parm])
+}
diff --git a/R/mc_core_pearson.R b/R/mc_core_pearson.R
index 82e169391298484123233a3d84432cd61f9df845..22719d872533dbf7a0f8b6a03d346f001a2123d6 100644
--- a/R/mc_core_pearson.R
+++ b/R/mc_core_pearson.R
@@ -6,9 +6,8 @@
 #' @param res A vector of residuals.
 #' @return A vector
 #' @details It is an internal function.
-#' @export
 
 mc_core_pearson <- function(product, inv_C, res) {
     output <- t(res) %*% product %*% (inv_C %*% res) - sum(diag(product))
     return(as.numeric(output))
-} 
+}
diff --git a/R/mc_correction.R b/R/mc_correction.R
index b57edcec312abfa8d4e96712cc327059af5e99c2..aebd8ef9da4ec17e747fd83b15b9384fd4f6cc4a 100644
--- a/R/mc_correction.R
+++ b/R/mc_correction.R
@@ -5,14 +5,13 @@
 #' @param D_C A list of matrices.
 #' @param inv_J_beta A matrix. In general it is computed based on the output of the
 #' \code{[mcglm]{mc_quasi_score}}.
-#' @param D A matrix. In general it is the output of the \code{[mcglm]{mc_link_function}}.
-#' @param inv_C A matrix. In general the output of the \code{[mcglm]{mc_build_C}}.
+#' @param D A matrix. In general it is the output of the \link{mc_link_function}.
+#' @param inv_C A matrix. In general the output of the \link{mc_build_C}.
 #' @return A vector with the correction terms to be used on the Pearson estimating function.
 #' @details It is an internal function useful inside the fitting algorithm.
-#' @export
 
 mc_correction <- function(D_C, inv_J_beta, D, inv_C) {
     term1 <- lapply(D_C, mc_sandwich, bord1 = t(D) %*% inv_C, bord2 = inv_C %*% D)
     output <- lapply(term1, function(x, inv_J_beta) sum(x * inv_J_beta), inv_J_beta = inv_J_beta)
     return(unlist(output))
-} 
+}
diff --git a/R/mc_cross_sensitivity.R b/R/mc_cross_sensitivity.R
index 23a3bf20b838821d9a4152e9e84be32da338bc65..5df0b1f7ab875cd3ad6e411421e0c5c182472e95 100644
--- a/R/mc_cross_sensitivity.R
+++ b/R/mc_cross_sensitivity.R
@@ -6,7 +6,6 @@
 #' @param Product_beta A list of matrices.
 #' @param n_beta_effective Numeric. Effective number of regression parameters.
 #' @return The cross-sensitivity matrix. Equation (10) of Bonat and Jorgensen (2015).
-#' @export
 
 mc_cross_sensitivity <- function(Product_cov, Product_beta, n_beta_effective = length(Product_beta)) {
     n_beta <- length(Product_beta)
diff --git a/R/mc_cross_variability.R b/R/mc_cross_variability.R
index de7663a8c76732b289ff9eae29f8180e7640cb5e..ae9da6ac0f681f8a2c9c82be893da1f3bdb9bdba 100644
--- a/R/mc_cross_variability.R
+++ b/R/mc_cross_variability.R
@@ -7,7 +7,7 @@
 #' @param res A vector.
 #' @param D A matrix.
 #' @return The cross-variability matrix between regression and covariance parameters.
-#' @export
+
 mc_cross_variability <- function(Product_cov, inv_C, res, D) {
     Wlist <- lapply(Product_cov, mc_multiply2, bord2 = inv_C)
     A = t(D) %*% inv_C
diff --git a/R/mc_derivative_C_rho.R b/R/mc_derivative_C_rho.R
index c31bd0a9505e9dd3285b6672fd800c2426269fd2..bcbbcaaafcb376611912236f97a143a65dec73f1 100644
--- a/R/mc_derivative_C_rho.R
+++ b/R/mc_derivative_C_rho.R
@@ -8,7 +8,7 @@
 #'@param II A diagonal matrix.
 #'@return A matrix.
 #'@details It is an internal function used to build the derivatives of the C matrix.
-#'@export
+
 mc_derivative_C_rho <- function(D_Sigmab, Bdiag_chol_Sigma_within, t_Bdiag_chol_Sigma_within, II) {
     output <- lapply(D_Sigmab, function(x) {
         t_Bdiag_chol_Sigma_within %*% kronecker(x, II) %*% Bdiag_chol_Sigma_within
diff --git a/R/mc_derivative_cholesky.R b/R/mc_derivative_cholesky.R
index 5dd3d3ecc305168af901c8cefaa5ac03e183fe28..5defd4e5ef8cdc78db5ab04cdd82f7e5c494ad94 100644
--- a/R/mc_derivative_cholesky.R
+++ b/R/mc_derivative_cholesky.R
@@ -6,7 +6,7 @@
 #' @param chol_Sigma A matrix.
 #' @return A list of matrix.
 #' @details It is an internal function.
-#' @export
+
 mc_derivative_cholesky <- function(derivada, inv_chol_Sigma, chol_Sigma) {
     faux <- function(derivada, inv_chol_Sigma, chol_Sigma) {
         t1 <- inv_chol_Sigma %*% derivada %*% t(inv_chol_Sigma)
@@ -17,4 +17,4 @@ mc_derivative_cholesky <- function(derivada, inv_chol_Sigma, chol_Sigma) {
     }
     list_D_chol <- lapply(derivada, faux, inv_chol_Sigma = inv_chol_Sigma, chol_Sigma = chol_Sigma)
     return(list_D_chol)
-} 
+}
diff --git a/R/mc_derivative_expm.R b/R/mc_derivative_expm.R
index dad9c665df33723beb4ee86f255f6f02dc7de932..46bde0a5f2b83a524c157b4b9cf345b93eac78a6 100644
--- a/R/mc_derivative_expm.R
+++ b/R/mc_derivative_expm.R
@@ -14,7 +14,7 @@
 #' @details Many arguments required by this function are provide by the \code{link[mcglm]{mc_dexpm}}.
 #' The argument dU is the derivative of the U matrix with respect to the models parameters. It should
 #' be computed by the user.
-#' @export
+
 mc_derivative_expm <- function(dU, UU, inv_UU, Q, n = dim(UU)[1], sparse = FALSE) {
     H = inv_UU %*% dU %*% UU
     P <- Matrix(0, ncol = n, nrow = n)
@@ -24,4 +24,4 @@ mc_derivative_expm <- function(dU, UU, inv_UU, Q, n = dim(UU)[1], sparse = FALSE
     P <- forceSymmetric(P)
     D_Omega = Matrix(UU %*% P %*% inv_UU, sparse = FALSE)
     return(D_Omega)
-} 
+}
diff --git a/R/mc_derivative_sigma_beta.R b/R/mc_derivative_sigma_beta.R
index 5757aa17259180c777669f476e8457c910989b87..fee51bd54f0a834afde05b67d1c45807e757f4a0 100644
--- a/R/mc_derivative_sigma_beta.R
+++ b/R/mc_derivative_sigma_beta.R
@@ -9,7 +9,7 @@
 #' @param variance A string specifying the variance function name.
 #' @return A list of matrices, containg the derivatives of V^{1/2} with respect to the regression
 #' parameters.
-#' @export
+
 mc_derivative_sigma_beta <- function(D, D_V_sqrt_mu, Omega, V_sqrt, variance) {
     n_beta <- dim(D)[2]
     n_obs <- dim(D)[1]
@@ -27,4 +27,4 @@ mc_derivative_sigma_beta <- function(D, D_V_sqrt_mu, Omega, V_sqrt, variance) {
         }
     }
     return(output)
-} 
+}
diff --git a/R/mc_dexp_gold.R b/R/mc_dexp_gold.R
index 4f241e5bf636b94bd51c1b5cce0512180d0e3ff8..77ff1a0b10abc5f7984e0615a882fbc8211b5a38 100644
--- a/R/mc_dexp_gold.R
+++ b/R/mc_dexp_gold.R
@@ -12,8 +12,9 @@
 #' @examples
 #' M <- matrix(c(1,0.8,0.8,1), 2,2)
 #' dM <- matrix(c(0,1,1,0),2,2)
-#' mc_dexp_gold(M = M, dM = dM)
+#' mcglm::mc_dexp_gold(M = M, dM = dM)
 #' @export
+
 mc_dexp_gold <- function(M, dM) {
     N = dim(M)
     AM = rbind(cbind(M, matrix(0, N[1], N[2])), cbind(dM, M))
@@ -21,4 +22,4 @@ mc_dexp_gold <- function(M, dM) {
     F = PSI[1:N[1], 1:N[1]]
     dF = PSI[c(N[1] + 1):c(2 * N[1]), 1:N[1]]
     return(list(F, dF))
-} 
+}
diff --git a/R/mc_dexpm.R b/R/mc_dexpm.R
index 8623496a0d775d41ae790d47fd770ca0f8a886f9..212410d6ba63646de339f6ca7f87e861b7de7b84 100644
--- a/R/mc_dexpm.R
+++ b/R/mc_dexpm.R
@@ -15,7 +15,6 @@
 #'
 #'@seealso \code{\link[Matrix]{expm}}, \code{\link[base]{eigen}},
 #'\code{link[mcglm]{mc_dexp_gold}}.
-#'@export
 
 mc_expm <- function(U, n = dim(U)[1], sparse = FALSE, inverse = FALSE) {
     tt = eigen(U, symmetric = TRUE)
diff --git a/R/mc_fitted.mcglm.R b/R/mc_fitted.mcglm.R
index 2c212c9812ec42a791a526bb0870f9358805b650..ac02339009eae48df18db120a9e4401dd81bf1d7 100644
--- a/R/mc_fitted.mcglm.R
+++ b/R/mc_fitted.mcglm.R
@@ -3,10 +3,13 @@
 #' @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) {
+
+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_getInformation.R b/R/mc_getInformation.R
index 49cf2d0c48beed444fcd105cb8c6c31acf95dfb6..946f0444ebb63139524fc7bdbd0a37d85e24d0cf 100644
--- a/R/mc_getInformation.R
+++ b/R/mc_getInformation.R
@@ -6,7 +6,7 @@
 #' @param list_power_fixed A list of logical specyfing if the power parameters should be estimated or not.
 #' @param n_resp A number specyfing the nmber of response variables.
 #' @return The number of \eqn{\beta}'s, \eqn{\tau}'s, power and correlation parameters.
-#' @export
+
 mc_getInformation <- function(list_initial, list_power_fixed, n_resp) {
     n_betas <- lapply(list_initial$regression, length)
     n_taus <- lapply(list_initial$tau, length)
@@ -25,4 +25,4 @@ mc_getInformation <- function(list_initial, list_power_fixed, n_resp) {
     n_cov <- sum(do.call(c, n_power)) + n_rho + sum(do.call(c, n_taus))
     saida <- list(n_betas = n_betas, n_taus = n_taus, n_power = n_power, n_rho = n_rho, n_cov = n_cov)
     return(saida)
-} 
+}
diff --git a/R/mc_influence.R b/R/mc_influence.R
index 42db2c39e3331a001860d87ace195667be6f8696..51e71ba81b3dd42a6438f44ec2f9eedc29a71f4d 100644
--- a/R/mc_influence.R
+++ b/R/mc_influence.R
@@ -72,9 +72,9 @@ mc_influence <- function(object, id) {
     }
     DFBETA <- lapply(DFBETA_clust, as.matrix)
     DFBETA <- lapply(DFBETA, t)
-    DFBETA <- ldply(DFBETA, data.frame)
+    DFBETA <- plyr::ldply(DFBETA, data.frame)
     names(DFBETA) <- object$beta_names[[1]]
-    DFBETAOij <- ldply(DFBETAOij, data.frame)
+    DFBETAOij <- plyr::ldply(DFBETAOij, data.frame)
     names(DFBETAOij) <- object$beta_names[[1]]
     DCLSi <- apply(as.matrix(DFBETA), MARGIN = 1, FUN = padroniza, M = M)
     DCLOij <- apply(as.matrix(DFBETAOij), MARGIN = 1, FUN = padroniza, M = M)
diff --git a/R/mc_initial_values.R b/R/mc_initial_values.R
index 1e59bc21ea76cec30f122c0baa1157274f9e94ae..53f2dcef07476066962882edf0d1da27aea1b400 100644
--- a/R/mc_initial_values.R
+++ b/R/mc_initial_values.R
@@ -17,7 +17,8 @@
 #' @return Return a list of initial values to be used while fitting McGLMs.
 #' @export
 
-mc_initial_values <- function(linear_pred, matrix_pred, link, variance, covariance, offset, Ntrial, contrasts = NULL, data) {
+mc_initial_values <- function(linear_pred, matrix_pred, link, variance, covariance, offset, Ntrial,
+                              contrasts = NULL, data) {
     n_resp <- length(linear_pred)
     if (!is.null(contrasts)) {
         list_X <- list()
diff --git a/R/mc_link_function.R b/R/mc_link_function.R
index 660b254936957dc8487801db63767c137c5e2d3d..40100b88d5654579b32d660a8d732b3860670a55 100644
--- a/R/mc_link_function.R
+++ b/R/mc_link_function.R
@@ -31,11 +31,13 @@
 #' mc_link_function(beta = c(1,0.5), X = X, offset = NULL, link = 'log')
 #' mc_link_function(beta = c(1,0.5), X = X, offset = rep(10,5), link = 'identity')
 #' @export
+#' @import assertthat
+
 # Generic link function ---------------------------
 mc_link_function <- function(beta, X, offset, link) {
     assert_that(noNA(beta))
     assert_that(noNA(X))
-    if (!is.null(offset)) 
+    if (!is.null(offset))
         assert_that(noNA(offset))
     switch(link, logit = {
         output <- mc_logit(beta = beta, X = X, offset = offset)
@@ -174,4 +176,4 @@ mc_inverse <- function(beta, X, offset) {
     mu = make.link("inverse")$linkinv(eta = eta)
     Deri = make.link("inverse")$mu.eta(eta = eta)
     return(list(mu = mu, D = X * Deri))
-} 
+}
diff --git a/R/mc_list2vec.R b/R/mc_list2vec.R
index 8584e39cb4cd2669be32ae51585d5a3a131715d3..fa06b41697c4ed181c118f86bb6f3f4e33055eef 100644
--- a/R/mc_list2vec.R
+++ b/R/mc_list2vec.R
@@ -7,7 +7,7 @@
 #' @return A vector of model parameters.
 #' @details It is an internal function, in general the users never will use this function.
 #' It will be useful, only if the user wants to implement a different variance-covariance matrix.
-#' @export
+
 mc_list2vec <- function(list_initial, list_power_fixed) {
     cov_ini <- do.call(c, Map(c, list_initial$power, list_initial$tau))
     n_resp <- length(list_initial$regression)
@@ -30,4 +30,4 @@ mc_list2vec <- function(list_initial, list_power_fixed) {
     beta_ini <- do.call(c, list_initial$regression)
     cov_ini <- c(rho = list_initial$rho, cov_ini)
     return(list(beta_ini = as.numeric(beta_ini), cov_ini = as.numeric(cov_ini)))
-} 
+}
diff --git a/R/mc_main_function.R b/R/mc_main_function.R
new file mode 100644
index 0000000000000000000000000000000000000000..eda852cae44d13414286c4859b4bbe23ecb285b9
--- /dev/null
+++ b/R/mc_main_function.R
@@ -0,0 +1,107 @@
+#' Fitting Multivariate Covariance Generalized Linear Models (McGLM)
+#'
+#' @description \code{mcglm} is used to fit multivariate covariance generalized linear models.
+#' The models are specified by a set of lists giving a symbolic description of the linear predictor.
+#' The user can choose between a list of link, variance and covariance functions. The models are
+#' fitted using an estimating function approach, combining quasi-score functions for regression
+#' parameters and Pearson estimating function for covariance parameters. For details see Bonat and
+#' Jorgensen (2015).
+#'
+#' @param linear_pred A list of formula see \code{\link[stats]{formula}} for details.
+#' @param matrix_pred A list of known matrices to be used on the matrix linear predictor. Details
+#' can be obtained on \code{\link[mcglm]{mc_matrix_linear_predictor}}.
+#' @param link A list of link functions names, see \code{\link[mcglm]{mc_link_function}} for details.
+#' @param variance A list of variance functions names, see \code{\link[mcglm]{mc_variance_function}}
+#' for details.
+#' @param covariance A list of covariance link functions names, current options are: identity, inverse
+#' and exponential-matrix (expm).
+#' @param offset A list with values of offset values if any.
+#' @param Ntrial A list with values of the number of trials on Bernoulli experiments. It is useful only
+#' for binomialP and binomialPQ variance functions.
+#' @param power_fixed A list of logicals indicating if the values of the power parameter should be
+#' estimated or not.
+#' @param control_initial A list of initial values for the fitting algorithm. See details below.
+#' @param control_algorithm A list of arguments to be passed for the fitting algorithm. See
+#' \code{\link[mcglm]{fit_mcglm}} for details.
+#' @param contrasts Extra arguments to passed to \code{\link[stats]{model.matrix}}.
+#' @param data A dta frame.
+#' @return mcglm returns an object of class 'mcglm'.
+#' @export
+#' @import Matrix
+mcglm <- function(linear_pred, matrix_pred, link, variance, covariance, offset,
+                  Ntrial, power_fixed, data, control_initial = "automatic",
+                  contrasts = NULL, control_algorithm = list()) {
+  n_resp <- length(linear_pred)
+  linear_pred <- as.list(linear_pred)
+  matrix_pred <- as.list(matrix_pred)
+  if (missing(link)) {
+    link = rep("identity", n_resp)
+  }
+  if (missing(variance)) {
+    variance = rep("constant", n_resp)
+  }
+  if (missing(covariance)) {
+    covariance = rep("identity", n_resp)
+  }
+  if (missing(offset)) {
+    offset = rep(list(NULL), n_resp)
+  }
+  if (missing(Ntrial)) {
+    Ntrial = rep(list(rep(1, dim(data)[1])),n_resp)
+  }
+  if (missing(power_fixed)) {
+    power_fixed <- rep(TRUE, n_resp)
+  }
+  if (missing(contrasts)) {
+    constrasts = NULL
+  }
+  link <- as.list(link)
+  variance <- as.list(variance)
+  covariance <- as.list(covariance)
+  offset <- as.list(offset)
+  Ntrial <- as.list(Ntrial)
+  power_fixed = as.list(power_fixed)
+  if (class(control_initial) != "list") {
+    control_initial <- mc_initial_values(linear_pred = linear_pred, matrix_pred = matrix_pred, link = link, variance = variance,
+                                         covariance = covariance, offset = offset, Ntrial = Ntrial, contrasts = contrasts, data = data)
+    cat("Automatic initial values selected.")
+  }
+  con <- list(correct = TRUE, max_iter = 20, tol = 1e-04, method = "chaser", tunning = 1, verbose = FALSE)
+  con[(namc <- names(control_algorithm))] <- control_algorithm
+  if (!is.null(contrasts)) {
+    list_X <- list()
+    for (i in 1:n_resp) {
+      list_X[[i]] <- model.matrix(linear_pred[[i]], contrasts = contrasts[[i]], data = data)
+    }
+  } else {
+    list_X <- lapply(linear_pred, model.matrix, data = data)
+  }
+
+  list_model_frame <- lapply(linear_pred, model.frame, data = data)
+  list_Y <- lapply(list_model_frame, model.response)
+  y_vec <- as.numeric(do.call(c, list_Y))
+  sparse <- lapply(matrix_pred, function(x) {
+    if (class(x) == "dgeMatrix") {
+      FALSE
+    } else TRUE
+  })
+  model_fit <- try(fit_mcglm(list_initial = control_initial, list_link = link, list_variance = variance, list_covariance = covariance,
+                             list_X = list_X, list_Z = matrix_pred, list_offset = offset, list_Ntrial = Ntrial, list_power_fixed = power_fixed,
+                             list_sparse = sparse, y_vec = y_vec, correct = con$correct, max_iter = con$max_iter, tol = con$tol, method = con$method,
+                             tunning = con$tunning, verbose = con$verbose))
+  if (class(model_fit) != "try-error") {
+    model_fit$beta_names <- lapply(list_X, colnames)
+    model_fit$power_fixed <- power_fixed
+    model_fit$list_initial <- control_initial
+    model_fit$n_obs <- dim(data)[1]
+    model_fit$link <- link
+    model_fit$variance <- variance
+    model_fit$covariance <- covariance
+    model_fit$linear_pred <- linear_pred
+    model_fit$con <- con
+    model_fit$observed <- Matrix(y_vec, ncol = length(list_Y), nrow = dim(data)[1])
+    model_fit$list_X <- list_X
+    class(model_fit) <- "mcglm"
+  }
+  return(model_fit)
+}
diff --git a/R/mc_matrix_linear_predictor.R b/R/mc_matrix_linear_predictor.R
index b50d1577321f33b46be80a123c30216c413b02c8..32016c9587fa8ebe110d5e7b7254e5d94c8feda7 100644
--- a/R/mc_matrix_linear_predictor.R
+++ b/R/mc_matrix_linear_predictor.R
@@ -6,9 +6,11 @@
 #' @param Z   A list of known matrices.
 #' @return A matrix.
 #' @export
+#' @import Matrix
 #' @details Given a list with a set of known matrices (\eqn{Z_0,...,Z_D}) the function
 #' \code{mc_matrix_linear_predictor} returns \eqn{U = \tau_0 Z_0 + ... + \tau_D Z_D}.
 #' @examples
+#' require(Matrix)
 #' Z0 <- Diagonal(5, 1)
 #' Z1 <- Matrix(rep(1,5)%*%t(rep(1,5)))
 #' Z <- list(Z0, Z1)
@@ -22,4 +24,4 @@ mc_matrix_linear_predictor <- function(tau, Z) {
     output <- mapply("*", Z, tau, SIMPLIFY = FALSE)  ## Multiply each matrix by the parameter tau
     output <- Reduce("+", output)  ## Sum all matrices inside the list
     return(output)
-} 
+}
diff --git a/R/mc_pearson.R b/R/mc_pearson.R
index 1c0748353a09046df5d43340441627c88f02d1e0..705c2ae44e4dc18e5c4bdee7680d6d2a7c5d3934 100644
--- a/R/mc_pearson.R
+++ b/R/mc_pearson.R
@@ -13,7 +13,7 @@
 #' and (iii) variability matrices associated with the Pearson estimating function.
 #' @details Compute the Pearson estimating function its sensitivity and variability matrices.
 #' For more details see Bonat and Jorgensen (2015) equations 6, 7 and 8.
-#' @export
+
 mc_pearson <- function(y_vec, mu_vec, Cfeatures, inv_J_beta = NULL, D = NULL, correct = FALSE, compute_variability = FALSE) {
     product <- lapply(Cfeatures$D_C, mc_multiply, bord2 = Cfeatures$inv_C)
     res <- y_vec - mu_vec
@@ -25,9 +25,9 @@ mc_pearson <- function(y_vec, mu_vec, Cfeatures, inv_J_beta = NULL, D = NULL, co
         output <- list(Score = pearson_score + correction, Sensitivity = sensitivity, Extra = product)
     }
     if (compute_variability == TRUE) {
-        variability <- mc_variability(sensitivity = sensitivity, product = product, inv_C = Cfeatures$inv_C, C = Cfeatures$C, 
+        variability <- mc_variability(sensitivity = sensitivity, product = product, inv_C = Cfeatures$inv_C, C = Cfeatures$C,
             res = res)
         output$Variability <- variability
     }
     return(output)
-} 
+}
diff --git a/R/mc_plot.mcglm.R b/R/mc_plot.mcglm.R
index 7fb6442ec4762457a589bd353befc048793e73e9..53dd62f1139eaf558f443d1949f2b6f9b1d3363f 100644
--- a/R/mc_plot.mcglm.R
+++ b/R/mc_plot.mcglm.R
@@ -3,11 +3,15 @@
 #' @description takes a fitted mcglm object by \code{mcglm()} and plots, residuals,
 #' influence diagnostic measures and algorithm check.
 #'
-#' @param object a fitted mcglm object as produced by \code{mcglm()}.
+#' @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(object, type = "residuals") {
+
+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))
diff --git a/R/mc_print.mcglm.R b/R/mc_print.mcglm.R
index 7cc8ae98ffa10cc758a7eb5fa35118c6097572c3..5e54ccfb0eea8932667f4c473c23a069824145fe 100644
--- a/R/mc_print.mcglm.R
+++ b/R/mc_print.mcglm.R
@@ -1,11 +1,14 @@
-#' Print a Multivariate Covariance Generalized Linear Model
+#' Print method for Multivariate Covariance Generalized Linear Model
 #'
 #' @description The default print method for a mcglm object.
-#'
-#' @param object fitted model objects of class mcglm as produced by mcglm().
+#' @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(object) {
+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) {
diff --git a/R/mc_quasi_score.R b/R/mc_quasi_score.R
index 1f398bcfb6f1e6dc6d2395d90f5e25555c167e1c..6bdcee8db1fef74433f60e010519df7dda4f3a96 100644
--- a/R/mc_quasi_score.R
+++ b/R/mc_quasi_score.R
@@ -7,7 +7,7 @@
 #' @param y_vec A vector.
 #' @param mu_vec A vector.
 #' @return The quasi-score vector, the Sensivity and variability matrices.
-#' @export
+
 mc_quasi_score <- function(D, inv_C, y_vec, mu_vec) {
     res <- y_vec - mu_vec
     t_D <- t(D)
diff --git a/R/mc_residuals.mcglm.R b/R/mc_residuals.mcglm.R
index 0a6b1b7e88cb03bcff13f299c91e7509e550d13a..0732c56b41d0e03c29798d3d274a367479815d43 100644
--- a/R/mc_residuals.mcglm.R
+++ b/R/mc_residuals.mcglm.R
@@ -5,10 +5,13 @@
 #' @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") {
+
+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") {
@@ -18,4 +21,4 @@ residuals.mcglm <- function(object, type = "raw") {
         output <- Matrix(as.numeric(object$residuals/sqrt(diag(object$C))), ncol = n_resp, nrow = object$n_obs)
     }
     return(output)
-} 
+}
diff --git a/R/mc_sensitivity.R b/R/mc_sensitivity.R
index 8b2e8e7fbca9a2638e0e88a3149077fcfad34e76..169cc185e31b6cd63b0a30ff8ed6262568b302ae 100644
--- a/R/mc_sensitivity.R
+++ b/R/mc_sensitivity.R
@@ -5,7 +5,7 @@
 #' @param product A list of matrix.
 #' @return The sensitivity matrix associated with the Pearson estimating function.
 #' @details This function implements the equation 7 of Bonat and Jorgensen (2015).
-#' @export
+
 mc_sensitivity <- function(product) {
     n_par <- length(product)
     Sensitivity <- matrix(0, n_par, n_par)
@@ -21,4 +21,4 @@ mc_sensitivity <- function(product) {
     # print(forceSymmetric(Sensitivity)) print(forceSymmetric(Sensitivity_temp)) print(forceSymmetric(Sensitivity1))
     # print(all.equal(Sensitivity1, Sensitivity))
     return(Sensitivity)
-} 
+}
diff --git a/R/mc_summary.mcglm.R b/R/mc_summary.mcglm.R
index acf9ed91d4cd91c843880f36567422c5a37f34ad..015079da29c1596e8d03d0a9cbbd26ab817b3e4f 100644
--- a/R/mc_summary.mcglm.R
+++ b/R/mc_summary.mcglm.R
@@ -1,11 +1,14 @@
 #' Summarizing Multivariate Covariance Generalized Linear Models fits.
 #'
-#' @description Summary for mcglm objects.
+#' @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) {
+summary.mcglm <- function(object, ...) {
     n_resp <- length(object$beta_names)
     output <- list()
     for (i in 1:n_resp) {
diff --git a/R/mc_transform_list_bdiag.R b/R/mc_transform_list_bdiag.R
index d6c82cd5ec6b73e4478b1108fde8eaadec4da8b7..c9765be5838568430c3c8baabf00b1040545dded 100644
--- a/R/mc_transform_list_bdiag.R
+++ b/R/mc_transform_list_bdiag.R
@@ -8,7 +8,7 @@
 #' \code{link[mcglm]{mc_build_bdiag}}.
 #' @param response_number A numeric specifying the response variable number.
 #' @return A list of block-diagonal matrices.
-#' @export
+
 mc_transform_list_bdiag <- function(list_mat, mat_zero, response_number) {
     aux.f <- function(x, mat_zero, response_number) {
         mat_zero[[response_number]] <- x
@@ -16,4 +16,4 @@ mc_transform_list_bdiag <- function(list_mat, mat_zero, response_number) {
     }
     output <- lapply(list_mat, aux.f, mat_zero = mat_zero, response_number = response_number)
     return(output)
-} 
+}
diff --git a/R/mc_updatedBeta.R b/R/mc_updatedBeta.R
index e71f14c425b51a0b57153503d10a63c14953e71e..3f257a3766646163cbdd27afafe2cf8648e29a8e 100644
--- a/R/mc_updatedBeta.R
+++ b/R/mc_updatedBeta.R
@@ -6,10 +6,10 @@
 #' @param list_initial A list of initial values.
 #' @param betas A vector with actual regression parameters values.
 #' @param information A list with information about the number of parameters in the model. In general
-#' the output from \code{[mcglm]{mc_getInformation}}.
+#' the output from \link{mc_getInformation}.
 #' @param n_resp A numeric specyfing the number of response variables.
 #' @return A list with updated values of the regression parameters.
-#' @export
+
 mc_updateBeta <- function(list_initial, betas, information, n_resp) {
     cod <- rep(1:n_resp, information$n_betas)
     temp <- data.frame(beta = betas, cod)
@@ -17,4 +17,4 @@ mc_updateBeta <- function(list_initial, betas, information, n_resp) {
         list_initial$regression[[k]] <- temp[which(temp$cod == k), ]$beta
     }
     return(list_initial)
-} 
+}
diff --git a/R/mc_updatedCov.R b/R/mc_updatedCov.R
index a28ae7f5a1438d152416fa47c0710f2e306c7dfb..1610d7fd8b04ed6b120eccbf90f48fa99118dbfa 100644
--- a/R/mc_updatedCov.R
+++ b/R/mc_updatedCov.R
@@ -6,12 +6,12 @@
 #' @param list_initial A list of initial values.
 #' @param covariance A vector with actual covariance parameters values.
 #' @param information A list with information about the number of parameters in the model. In general
-#' the output from \code{[mcglm]{mc_getInformation}}.
+#' the output from \link{mc_getInformation}.
 #' @param list_power_fixed A list of logicals indicating if the power parameter should be
 #' estimated or not.
 #' @param n_resp A numeric specyfing the number of response variables.
 #' @return A list with updated values of the covariance parameters.
-#' @export
+
 mc_updateCov <- function(list_initial, covariance, list_power_fixed, information, n_resp) {
     rho_cod <- rep("rho", information$n_rho)
     tau_cod <- list()
@@ -35,4 +35,4 @@ mc_updateCov <- function(list_initial, covariance, list_power_fixed, information
         list_initial$rho <- temp[which(temp$cod == "rho"), ]$values
     }
     return(list_initial)
-} 
+}
diff --git a/R/mc_variability.R b/R/mc_variability.R
index 6775da353d16de03e9ac56fa9f121e6a94b99103..8bce489091d8a28228d5d718e3d37b6419798d4c 100644
--- a/R/mc_variability.R
+++ b/R/mc_variability.R
@@ -2,14 +2,14 @@
 #'
 #' @description Compute the variability matrix associated with the Pearson estimating function.
 #'
-#' @param sensitivity A matrix. In general the output from \code{[mcglm]{mc_sensitivity}}.
+#' @param sensitivity A matrix. In general the output from \code{mc_sensitivity}.
 #' @param product A list of matrix.
-#' @param inv_C A matrix. In general the output from \code{[mcglm]{mc_build_C}}.
-#' @param C A matrix. In general the output from \code{[mcglm]{mc_build_C}}.
+#' @param inv_C A matrix. In general the output from \code{mc_build_C}.
+#' @param C A matrix. In general the output from \code{mc_build_C}.
 #' @param res A vector. The residuals vector, i.e. (y_vec - mu_vec).
 #' @return The variability matrix associated witht the Pearson estimating function.
 #' @details This function implements the equation 8 of Bonat and Jorgensen (2015).
-#' @export
+
 mc_variability <- function(sensitivity, product, inv_C, C, res) {
     W <- lapply(product, mc_multiply2, bord2 = inv_C)
     n_par <- length(product)
@@ -22,4 +22,4 @@ mc_variability <- function(sensitivity, product, inv_C, C, res) {
         }
     }
     return(Variability)
-} 
+}
diff --git a/R/mc_vcov.R b/R/mc_vcov.R
index f3ea538b4c60fc84c81722b707936227b646c6d1..9448b399a6e907fd159dae9e41294d4dc33761e9 100644
--- a/R/mc_vcov.R
+++ b/R/mc_vcov.R
@@ -3,12 +3,14 @@
 #' @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) {
+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 8af85309457af8e163e69a577be2cf0b9f1e3541..8497592ed3ef67143e16cc7340eb7a8ccac2c285 100644
--- a/R/mcglm.R
+++ b/R/mcglm.R
@@ -1,106 +1,67 @@
-#' Fitting Multivariate Covariance Generalized Linear Models (McGLM)
+##' @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.
+##'
+##'
+##' @docType package
+##' @name mcglm
+NULL
+
+#' @name ahs
 #'
-#' @description \code{mcglm} is used to fit multivariate covariance generalized linear models.
-#' The models are specified by a set of lists giving a symbolic description of the linear predictor.
-#' The user can choose between a list of link, variance and covariance functions. The models are
-#' fitted using an estimating function approach, combining quasi-score functions for regression
-#' parameters and Pearson estimating function for covariance parameters. For details see Bonat and
-#' Jorgensen (2015).
+#' @title Australian health survey
 #'
-#' @param linear_pred A list of formula see \code{\link[stats]{formula}} for details.
-#' @param matrix_pred A list of known matrices to be used on the matrix linear predictor. Details
-#' can be obtained on \code{\link[mcglm]{mc_matrix_linear_predictor}}.
-#' @param link A list of link functions names, see \code{\link[mcglm]{mc_link_function}} for details.
-#' @param variance A list of variance functions names, see \code{\link[mcglm]{mc_variance_function}}
-#' for details.
-#' @param covariance A list of covariance link functions names, current options are: identity, inverse
-#' and exponential-matrix (expm).
-#' @param offset A list with values of offset values if any.
-#' @param Ntrial A list with values of the number of trials on Bernoulli experiments. It is useful only
-#' for binomialP and binomialPQ variance functions.
-#' @param power_fixed A list of logicals indicating if the values of the power parameter should be
-#' estimated or not.
-#' @param control_initial A list of initial values for the fitting algorithm. See details below.
-#' @param control_algorithm A list of arguments to be passed for the fitting algorithm. See
-#' \code{\link[mcglm]{fit_mcglm}} for details.
-#' @param contrasts Extra arguments to passed to \code{\link[stats]{model.matrix}}.
-#' @param data A dta frame.
-#' @return mcglm returns an object of class 'mcglm'.
-#' @export
-mcglm <- function(linear_pred, matrix_pred, link, variance, covariance, offset, Ntrial,
-                  power_fixed, data, control_initial = "automatic",
-                  contrasts = NULL, control_algorithm = list()) {
-    n_resp <- length(linear_pred)
-    linear_pred <- as.list(linear_pred)
-    matrix_pred <- as.list(matrix_pred)
-    if (missing(link)) {
-        link = rep("identity", n_resp)
-    }
-    if (missing(variance)) {
-        variance = rep("constant", n_resp)
-    }
-    if (missing(covariance)) {
-        covariance = rep("identity", n_resp)
-    }
-    if (missing(offset)) {
-        offset = rep(list(NULL), n_resp)
-    }
-    if (missing(Ntrial)) {
-        Ntrial = rep(list(rep(1, dim(data)[1])),n_resp)
-    }
-    if (missing(power_fixed)) {
-        power_fixed <- rep(TRUE, n_resp)
-    }
-    if (missing(contrasts)) {
-        constrasts = NULL
-    }
-    link <- as.list(link)
-    variance <- as.list(variance)
-    covariance <- as.list(covariance)
-    offset <- as.list(offset)
-    Ntrial <- as.list(Ntrial)
-    power_fixed = as.list(power_fixed)
-    if (class(control_initial) != "list") {
-        control_initial <- mc_initial_values(linear_pred = linear_pred, matrix_pred = matrix_pred, link = link, variance = variance,
-            covariance = covariance, offset = offset, Ntrial = Ntrial, contrasts = contrasts, data = data)
-        cat("Automatic initial values selected.")
-    }
-    con <- list(correct = TRUE, max_iter = 20, tol = 1e-04, method = "chaser", tunning = 1, verbose = FALSE)
-    con[(namc <- names(control_algorithm))] <- control_algorithm
-    if (!is.null(contrasts)) {
-        list_X <- list()
-        for (i in 1:n_resp) {
-            list_X[[i]] <- model.matrix(linear_pred[[i]], contrasts = contrasts[[i]], data = data)
-        }
-    } else {
-        list_X <- lapply(linear_pred, model.matrix, data = data)
-    }
-
-    list_model_frame <- lapply(linear_pred, model.frame, data = data)
-    list_Y <- lapply(list_model_frame, model.response)
-    y_vec <- as.numeric(do.call(c, list_Y))
-    sparse <- lapply(matrix_pred, function(x) {
-        if (class(x) == "dgeMatrix") {
-            FALSE
-        } else TRUE
-    })
-    model_fit <- try(fit_mcglm(list_initial = control_initial, list_link = link, list_variance = variance, list_covariance = covariance,
-        list_X = list_X, list_Z = matrix_pred, list_offset = offset, list_Ntrial = Ntrial, list_power_fixed = power_fixed,
-        list_sparse = sparse, y_vec = y_vec, correct = con$correct, max_iter = con$max_iter, tol = con$tol, method = con$method,
-        tunning = con$tunning, verbose = con$verbose))
-    if (class(model_fit) != "try-error") {
-        model_fit$beta_names <- lapply(list_X, colnames)
-        model_fit$power_fixed <- power_fixed
-        model_fit$list_initial <- control_initial
-        model_fit$n_obs <- dim(data)[1]
-        model_fit$link <- link
-        model_fit$variance <- variance
-        model_fit$covariance <- covariance
-        model_fit$linear_pred <- linear_pred
-        model_fit$con <- con
-        model_fit$observed <- Matrix(y_vec, ncol = length(list_Y), nrow = dim(data)[1])
-        model_fit$list_X <- list_X
-        class(model_fit) <- "mcglm"
-    }
-    return(model_fit)
-}
+#' @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.
+#' }
+#'
+#' @docType data
+#'
+#' @keywords datasets
+#'
+#' @usage data(ahs)
+#'
+#' @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.
+#'
+NULL
diff --git a/data/data.rda b/data/data.rda
deleted file mode 100644
index 45849863b8a4166508f3bdd6a41d2467afdbc2a8..0000000000000000000000000000000000000000
Binary files a/data/data.rda and /dev/null differ
diff --git a/data/mcglm.rda b/data/mcglm.rda
deleted file mode 100644
index 30f7ffcc63641564f873dcb6308b2ea1f6303138..0000000000000000000000000000000000000000
Binary files a/data/mcglm.rda and /dev/null differ
diff --git a/man/ahs.Rd b/man/ahs.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..4569dd4eca3aaa17f45a2a28db2ebe070e5d7685
--- /dev/null
+++ b/man/ahs.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/mcglm.R
+\docType{data}
+\name{ahs}
+\alias{ahs}
+\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.
+}
+\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.
+
+\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.
+}
+}
+\keyword{datasets}
+
diff --git a/man/anova.mcglm.Rd b/man/anova.mcglm.Rd
index 3ba5b6be04e9892ad93f21d74fd3e09946b963b6..d367ebbf685bcb454cb9bf677f598f741b2458d2 100644
--- a/man/anova.mcglm.Rd
+++ b/man/anova.mcglm.Rd
@@ -4,10 +4,13 @@
 \alias{anova.mcglm}
 \title{ANOVA method for McGLMs.}
 \usage{
-\method{anova}{mcglm}(object)
+\method{anova}{mcglm}(object, ...)
 }
 \arguments{
 \item{object}{an object of class 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.}
 }
 \description{
 ANOVA method for McGLMS.
diff --git a/man/coef.mcglm.Rd b/man/coef.mcglm.Rd
index ad9b83da33440191eb5d449925ba4e4461d29dee..741388b2237b81fe9ebe23671621ddc3c8ece256 100644
--- a/man/coef.mcglm.Rd
+++ b/man/coef.mcglm.Rd
@@ -6,7 +6,7 @@
 \usage{
 \method{coef}{mcglm}(object, std.error = FALSE, response = c(NA,
   1:length(object$beta_names)), type = c("beta", "tau", "power",
-  "correlation"))
+  "correlation"), ...)
 }
 \arguments{
 \item{object}{An object of mcglm class.}
@@ -18,6 +18,9 @@ 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{...}{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.
diff --git a/man/confint.mcglm.Rd b/man/confint.mcglm.Rd
index 37dbe1db0e4ed3950f4d40f0f9b52d52a51d5c65..efdb6bb31b69cae0964019622e95df923c4ffd36 100644
--- a/man/confint.mcglm.Rd
+++ b/man/confint.mcglm.Rd
@@ -4,12 +4,19 @@
 \alias{confint.mcglm}
 \title{Confidence Intervals for mcglm}
 \usage{
-\method{confint}{mcglm}(object, level = 0.95)
+\method{confint}{mcglm}(object, parm, level = 0.95, ...)
 }
 \arguments{
 \item{object}{a fitted 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{level}{the confidence level required.}
+
+\item{...}{additional arguments affecting the confidence interval produced.
+Note that there is no extra options for
+mcglm object class.}
 }
 \value{
 A data.frame with confidence intervals, parameters names,
diff --git a/man/fitted.mcglm.Rd b/man/fitted.mcglm.Rd
index f0024801a480f5f72fa9af69f85cba950abd0a47..115d5810f201bc4391de693565a8e49413d55eb2 100644
--- a/man/fitted.mcglm.Rd
+++ b/man/fitted.mcglm.Rd
@@ -4,10 +4,13 @@
 \alias{fitted.mcglm}
 \title{Extract Model Fitted Values of McGLM}
 \usage{
-\method{fitted}{mcglm}(object)
+\method{fitted}{mcglm}(object, ...)
 }
 \arguments{
 \item{object}{An object of mcglm class.}
+
+\item{...}{additional arguments affecting the summary produced. Note that there is no extra options for
+mcglm object class.}
 }
 \value{
 Depending on the number of response variable the function \code{fitted.mcglm} returns
diff --git a/man/mc_correction.Rd b/man/mc_correction.Rd
index 39c10b02ef83c2a864a402b001a771a3b8730940..a24379d92bdde9919386b3324992026b9f905ceb 100644
--- a/man/mc_correction.Rd
+++ b/man/mc_correction.Rd
@@ -12,9 +12,9 @@ mc_correction(D_C, inv_J_beta, D, inv_C)
 \item{inv_J_beta}{A matrix. In general it is computed based on the output of the
 \code{[mcglm]{mc_quasi_score}}.}
 
-\item{D}{A matrix. In general it is the output of the \code{[mcglm]{mc_link_function}}.}
+\item{D}{A matrix. In general it is the output of the \link{mc_link_function}.}
 
-\item{inv_C}{A matrix. In general the output of the \code{[mcglm]{mc_build_C}}.}
+\item{inv_C}{A matrix. In general the output of the \link{mc_build_C}.}
 }
 \value{
 A vector with the correction terms to be used on the Pearson estimating function.
diff --git a/man/mc_dexp_gold.Rd b/man/mc_dexp_gold.Rd
index a51631ad77119c6a6dd075011a355aaffd928e4d..174642287aad7168da3793aea84f59b47173466e 100644
--- a/man/mc_dexp_gold.Rd
+++ b/man/mc_dexp_gold.Rd
@@ -23,7 +23,7 @@ function to test my own implementation based on eigen values decomposition.
 \examples{
 M <- matrix(c(1,0.8,0.8,1), 2,2)
 dM <- matrix(c(0,1,1,0),2,2)
-mc_dexp_gold(M = M, dM = dM)
+mcglm::mc_dexp_gold(M = M, dM = dM)
 }
 \seealso{
 \code{\link[Matrix]{expm}}, \code{\link[base]{eigen}}.
diff --git a/man/mc_matrix_linear_predictor.Rd b/man/mc_matrix_linear_predictor.Rd
index 0e56a8cf4ecf82ef7f3c5eee59979ebf065ef088..07ab6c50b7b04566d9894d2655d87e4da30c5c1c 100644
--- a/man/mc_matrix_linear_predictor.Rd
+++ b/man/mc_matrix_linear_predictor.Rd
@@ -22,6 +22,7 @@ Given a list with a set of known matrices (\eqn{Z_0,...,Z_D}) the function
 \code{mc_matrix_linear_predictor} returns \eqn{U = \tau_0 Z_0 + ... + \tau_D Z_D}.
 }
 \examples{
+require(Matrix)
 Z0 <- Diagonal(5, 1)
 Z1 <- Matrix(rep(1,5)\%*\%t(rep(1,5)))
 Z <- list(Z0, Z1)
diff --git a/man/mc_sandwich.Rd b/man/mc_sandwich.Rd
index 40e79a545df30326230bd32a5dd7f587bb91e6b8..f59c67b792a82f58723117028ac37c8b95b4d333 100644
--- a/man/mc_sandwich.Rd
+++ b/man/mc_sandwich.Rd
@@ -37,10 +37,4 @@ in the sandwich form bord1*middle*bord2. An special case appears when computing
 the covariance matrix with respect to the power parameter. Always the bord1 and bord2 should be
 diagonal matrix. If it is not true, this product is too slow.
 }
-\examples{
-X1 <- Diagonal(5, 2)
-M <- Matrix(rep(0.8,5)\%*\%t(rep(0.8,5)))
-mc_sandwich(middle = M, bord1 = X1, bord2 = X1)
-mc_sandwich_negative(middle = M, bord1 = X1, bord2 = X1)
-}
 
diff --git a/man/mc_updateBeta.Rd b/man/mc_updateBeta.Rd
index cd647bd58147a996b510eb22d8ca7f1f0ca276be..21e7027ccc47f9145c0492923eb134ff8f4f24e7 100644
--- a/man/mc_updateBeta.Rd
+++ b/man/mc_updateBeta.Rd
@@ -12,7 +12,7 @@ mc_updateBeta(list_initial, betas, information, n_resp)
 \item{betas}{A vector with actual regression parameters values.}
 
 \item{information}{A list with information about the number of parameters in the model. In general
-the output from \code{[mcglm]{mc_getInformation}}.}
+the output from \link{mc_getInformation}.}
 
 \item{n_resp}{A numeric specyfing the number of response variables.}
 }
diff --git a/man/mc_updateCov.Rd b/man/mc_updateCov.Rd
index 4cf72785692b0b4f9ef2dea2e15f58002a617d2d..c286f7b2fc02432b73a6faf7756fa34f22784ea0 100644
--- a/man/mc_updateCov.Rd
+++ b/man/mc_updateCov.Rd
@@ -15,7 +15,7 @@ mc_updateCov(list_initial, covariance, list_power_fixed, information, n_resp)
 estimated or not.}
 
 \item{information}{A list with information about the number of parameters in the model. In general
-the output from \code{[mcglm]{mc_getInformation}}.}
+the output from \link{mc_getInformation}.}
 
 \item{n_resp}{A numeric specyfing the number of response variables.}
 }
diff --git a/man/mc_variability.Rd b/man/mc_variability.Rd
index d52228ec6763e64097a8321e2b145b70536450b2..357dd3e1113ecad1c25dde441b7114a43441eac2 100644
--- a/man/mc_variability.Rd
+++ b/man/mc_variability.Rd
@@ -7,13 +7,13 @@
 mc_variability(sensitivity, product, inv_C, C, res)
 }
 \arguments{
-\item{sensitivity}{A matrix. In general the output from \code{[mcglm]{mc_sensitivity}}.}
+\item{sensitivity}{A matrix. In general the output from \code{mc_sensitivity}.}
 
 \item{product}{A list of matrix.}
 
-\item{inv_C}{A matrix. In general the output from \code{[mcglm]{mc_build_C}}.}
+\item{inv_C}{A matrix. In general the output from \code{mc_build_C}.}
 
-\item{C}{A matrix. In general the output from \code{[mcglm]{mc_build_C}}.}
+\item{C}{A matrix. In general the output from \code{mc_build_C}.}
 
 \item{res}{A vector. The residuals vector, i.e. (y_vec - mu_vec).}
 }
diff --git a/man/mcglm.Rd b/man/mcglm.Rd
index c61917fe7c94c64ccf63844802128a610608fe69..63c324fce7fbcaf6a8667b8cc84e3b626b68aa91 100644
--- a/man/mcglm.Rd
+++ b/man/mcglm.Rd
@@ -1,7 +1,9 @@
 % Generated by roxygen2 (4.1.1): do not edit by hand
-% Please edit documentation in R/mcglm.R
+% Please edit documentation in R/mc_main_function.R, R/mcglm.R
+\docType{package}
 \name{mcglm}
 \alias{mcglm}
+\alias{mcglm-package}
 \title{Fitting Multivariate Covariance Generalized Linear Models (McGLM)}
 \usage{
 mcglm(linear_pred, matrix_pred, link, variance, covariance, offset, Ntrial,
@@ -49,5 +51,17 @@ The user can choose between a list of link, variance and covariance functions. T
 fitted using an estimating function approach, combining quasi-score functions for regression
 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.
 }
 
diff --git a/man/plot.mcglm.Rd b/man/plot.mcglm.Rd
index 8ec40c94ecfe0dfc00e3c56dcb2f099ab2a0feb2..a1b852cc05cd1bff264060e921f0b7c96864cdb6 100644
--- a/man/plot.mcglm.Rd
+++ b/man/plot.mcglm.Rd
@@ -4,13 +4,16 @@
 \alias{plot.mcglm}
 \title{Default Multivariate Covariance Generalized Linear models plotting}
 \usage{
-\method{plot}{mcglm}(object, type = "residuals")
+\method{plot}{mcglm}(x, type = "residuals", ...)
 }
 \arguments{
-\item{object}{a fitted mcglm object as produced by \code{mcglm()}.}
+\item{x}{a fitted mcglm object as produced by \code{mcglm()}.}
 
 \item{type}{Specify which graphical analysis will be performed, options are: residuals, influence
 and algorithm.}
+
+\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,
diff --git a/man/print.mcglm.Rd b/man/print.mcglm.Rd
index 46218fa8eb3da33a0ec63a10d6de8f70e9e5f3e1..008141d86c9bb9cac15fb6e7b5427670d0d3eff6 100644
--- a/man/print.mcglm.Rd
+++ b/man/print.mcglm.Rd
@@ -2,12 +2,14 @@
 % Please edit documentation in R/mc_print.mcglm.R
 \name{print.mcglm}
 \alias{print.mcglm}
-\title{Print a Multivariate Covariance Generalized Linear Model}
+\title{Print method for Multivariate Covariance Generalized Linear Model}
 \usage{
-print.mcglm(object, ...)
+\method{print}{mcglm}(x, ...)
 }
 \arguments{
-\item{object}{fitted model objects of class mcglm as produced by mcglm().}
+\item{x}{fitted model objects of class mcglm as produced by mcglm().}
+
+\item{...}{further arguments passed to or from other methods.}
 }
 \description{
 The default print method for a mcglm object.
diff --git a/man/residuals.mcglm.Rd b/man/residuals.mcglm.Rd
index 1545596fef7247d7e0b3fa2340766934dd66bb82..0e5cc27de00e03907f2492c1e420db57674f0fa4 100644
--- a/man/residuals.mcglm.Rd
+++ b/man/residuals.mcglm.Rd
@@ -4,13 +4,16 @@
 \alias{residuals.mcglm}
 \title{Residuals for Multivariate Covariance Generalized Linear Models (McGLM)}
 \usage{
-\method{residuals}{mcglm}(object, type = "raw")
+\method{residuals}{mcglm}(object, type = "raw", ...)
 }
 \arguments{
 \item{object}{An of class 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{...}{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
diff --git a/man/summary.mcglm.Rd b/man/summary.mcglm.Rd
index 664168dd541c70417d334c789b931597fdbfec10..8e94d1de27313318a23c0653e52c4da6ec0cbb7b 100644
--- a/man/summary.mcglm.Rd
+++ b/man/summary.mcglm.Rd
@@ -4,12 +4,18 @@
 \alias{summary.mcglm}
 \title{Summarizing Multivariate Covariance Generalized Linear Models fits.}
 \usage{
-summary.mcglm(object)
+\method{summary}{mcglm}(object, ...)
 }
 \arguments{
 \item{object}{an object of class 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.}
+}
+\value{
+Print an mcglm object.
 }
 \description{
-Summary for mcglm objects.
+Summary for McGLMs objects.
 }
 
diff --git a/man/vcov.mcglm.Rd b/man/vcov.mcglm.Rd
index 68dcae6b18112d7926353389d7bdda64183c3b4a..31f20be578f2a7d118ded6e55b24a24c585511c4 100644
--- a/man/vcov.mcglm.Rd
+++ b/man/vcov.mcglm.Rd
@@ -4,10 +4,13 @@
 \alias{vcov.mcglm}
 \title{Calculate Variance-Covariance matrix for a fitted McGLM object.}
 \usage{
-\method{vcov}{mcglm}(object)
+\method{vcov}{mcglm}(object, ...)
 }
 \arguments{
 \item{object}{a fitted model mcglm object.}
+
+\item{...}{additional arguments affecting the summary produced. Note that there is no extra options for
+mcglm object class.}
 }
 \value{
 A variance-covariance matrix.