diff --git a/R/mc_core_cross_variability.R b/R/mc_core_cross_variability.R index 610bff7920a2d45c9a1ff354c5eec52bc5c291d3..2bff4076255cebb5c1fc2ee1f8ebb21f06ba7636 100644 --- a/R/mc_core_cross_variability.R +++ b/R/mc_core_cross_variability.R @@ -1,13 +1,16 @@ -#' Cross variability matrix +#' @title Cross variability matrix +#' @author Wagner Hugo Bonat #' -#' @description Compute the cross-covariance matrix between covariance and regression parameters. -#' Equation (11) of Bonat and Jorgensen (2015). +#' @description Compute the cross-covariance matrix between covariance +#' and regression parameters. Equation (11) of Bonat and Jorgensen +#' (2015). #' #' @param A A matrix. #' @param res A vector of residuals. #' @param W A matrix of weights. + covprod <- function(A, res, W) { - res = as.numeric(res) + res <- as.numeric(res) saida <- (res %*% W %*% res) %*% (t(res) %*% A) return(as.numeric(saida)) -} +} diff --git a/R/mc_core_pearson.R b/R/mc_core_pearson.R index 22719d872533dbf7a0f8b6a03d346f001a2123d6..cbe5e0d759713337abd2ea55c1657b32629d0c9f 100644 --- a/R/mc_core_pearson.R +++ b/R/mc_core_pearson.R @@ -1,4 +1,6 @@ -#' Core of the Pearson estimating function. +#' @title Core of the Pearson estimating function. +#' @author Wagner Hugo Bonat +#' #' @description Core of the Pearson estimating function. #' #' @param product A matrix. @@ -8,6 +10,7 @@ #' @details It is an internal function. mc_core_pearson <- function(product, inv_C, res) { - output <- t(res) %*% product %*% (inv_C %*% res) - sum(diag(product)) + 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 aebd8ef9da4ec17e747fd83b15b9384fd4f6cc4a..4154774865190bfb9bbd773529c8e2f5e87d6337 100644 --- a/R/mc_correction.R +++ b/R/mc_correction.R @@ -1,17 +1,26 @@ -#' Pearson correction term +#' @title Pearson correction term +#' @author Wagner Hugo Bonat #' -#' @description Compute the correction term associated with the Pearson estimating function. +#' @description Compute the correction term associated with the Pearson +#' estimating function. #' #' @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 \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. +#' @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 +#' \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. 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) + 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 5df0b1f7ab875cd3ad6e411421e0c5c182472e95..ba87b88875b5af6da5a482e1e66e869c8d51cf9e 100644 --- a/R/mc_cross_sensitivity.R +++ b/R/mc_cross_sensitivity.R @@ -1,23 +1,32 @@ -#' Cross-sensitivity +#' @title Cross-sensitivity +#' @author Wagner Hugo Bonat +#' +#' @description Compute the cross-sensitivity matrix between regression +#' and covariance parameters. Equation 10 of Bonat and Jorgensen +#' (2015). #' -#' @description Compute the cross-sensitivity matrix between regression and covariance parameters. -#' Equation 10 of Bonat and Jorgensen (2015). #' @param Product_cov A list of matrices. #' @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). +#' @param n_beta_effective Numeric. Effective number of regression +#' parameters. +#' @return The cross-sensitivity matrix. Equation (10) of Bonat and +#' Jorgensen (2015). -mc_cross_sensitivity <- function(Product_cov, Product_beta, n_beta_effective = length(Product_beta)) { +mc_cross_sensitivity <- function(Product_cov, Product_beta, + n_beta_effective = + length(Product_beta)) { n_beta <- length(Product_beta) n_cov <- length(Product_cov) if (n_beta == 0) { - cross_sensitivity <- Matrix(0, ncol = n_beta_effective, nrow = n_cov) + cross_sensitivity <- Matrix(0, ncol = n_beta_effective, + nrow = n_cov) } if (n_beta != 0) { cross_sensitivity <- Matrix(NA, nrow = n_cov, ncol = n_beta) for (i in 1:n_cov) { for (j in 1:n_beta) { - cross_sensitivity[i, j] <- -sum(Product_cov[[i]] * Product_beta[[j]]) + cross_sensitivity[i, j] <- + -sum(Product_cov[[i]] * Product_beta[[j]]) } } } diff --git a/R/mc_cross_variability.R b/R/mc_cross_variability.R index ae9da6ac0f681f8a2c9c82be893da1f3bdb9bdba..d6dc3de235e718c3119c4256c9512fbd77837950 100644 --- a/R/mc_cross_variability.R +++ b/R/mc_cross_variability.R @@ -1,22 +1,26 @@ -#' Compute the cross-variability matrix +#' @title Compute the cross-variability matrix +#' @author Wagner Hugo Bonat #' -#' @description Compute the cross-variability matrix between covariance and regression parameters. +#' @description Compute the cross-variability matrix between covariance +#' and regression parameters. #' #' @param Product_cov A list of matrices. #' @param inv_C A matrix. #' @param res A vector. #' @param D A matrix. -#' @return The cross-variability matrix between regression and covariance parameters. +#' @return The cross-variability matrix between regression and +#' covariance parameters. mc_cross_variability <- function(Product_cov, inv_C, res, D) { Wlist <- lapply(Product_cov, mc_multiply2, bord2 = inv_C) - A = t(D) %*% inv_C + A <- t(D) %*% inv_C n_beta <- dim(A)[1] n_cov <- length(Product_cov) cross_variability <- Matrix(NA, ncol = n_cov, nrow = n_beta) for (j in 1:n_beta) { for (i in 1:n_cov) { - cross_variability[j, i] <- covprod(A[j, ], Wlist[[i]], res = res) + cross_variability[j, i] <- + covprod(A[j, ], Wlist[[i]], res = res) } } return(cross_variability) diff --git a/R/mc_derivative_C_rho.R b/R/mc_derivative_C_rho.R index bcbbcaaafcb376611912236f97a143a65dec73f1..913920649a163c9667f04a73e3fa1674b0512d5d 100644 --- a/R/mc_derivative_C_rho.R +++ b/R/mc_derivative_C_rho.R @@ -1,17 +1,24 @@ -#' Derivative of C with respect to rho. +#' @title Derivative of C with respect to rho. +#' @author Wagner Hugo Bonat #' -#'@description Compute the derivative of the C matrix with respect to the correlation parameters rho. +#'@description Compute the derivative of the C matrix with respect to +#' the correlation parameters rho. #' #'@param D_Sigmab A matrix. #'@param Bdiag_chol_Sigma_within A block-diagonal matrix. #'@param t_Bdiag_chol_Sigma_within A block-diagonal matrix. #'@param II A diagonal matrix. #'@return A matrix. -#'@details It is an internal function used to build the derivatives of the C matrix. +#'@details It is an internal function used to build the derivatives of +#' the C matrix. -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 - }) +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 + }) return(output) } diff --git a/R/mc_derivative_cholesky.R b/R/mc_derivative_cholesky.R index 5defd4e5ef8cdc78db5ab04cdd82f7e5c494ad94..d48b1c3e8a7f79678fe1a9c09d1f540c9dc6899b 100644 --- a/R/mc_derivative_cholesky.R +++ b/R/mc_derivative_cholesky.R @@ -1,5 +1,8 @@ -#' Derivatives of the Cholesky decomposition -#' @description This function compute the derivative of the Cholesky decomposition. +#' @title Derivatives of the Cholesky decomposition +#' @author Wagner Hugo Bonat +#' +#' @description This function compute the derivative of the Cholesky +#' decomposition. #' #' @param derivada A matrix. #' @param inv_chol_Sigma A matrix. @@ -7,7 +10,8 @@ #' @return A list of matrix. #' @details It is an internal function. -mc_derivative_cholesky <- function(derivada, inv_chol_Sigma, chol_Sigma) { +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) t1 <- tril(t1) @@ -15,6 +19,8 @@ mc_derivative_cholesky <- function(derivada, inv_chol_Sigma, chol_Sigma) { output <- chol_Sigma %*% t1 return(output) } - list_D_chol <- lapply(derivada, faux, inv_chol_Sigma = inv_chol_Sigma, 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 46bde0a5f2b83a524c157b4b9cf345b93eac78a6..a40876285d0db27ac8bb9c06f0889f520d40506a 100644 --- a/R/mc_derivative_expm.R +++ b/R/mc_derivative_expm.R @@ -1,6 +1,8 @@ -#' Derivative of exponential-matrix function +#' @title Derivative of exponential-matrix function +#' @author Wagner Hugo Bonat #' -#' @description Compute the derivative of the exponential-matrix covariance link function. +#' @description Compute the derivative of the exponential-matrix +#' covariance link function. #' #' @param UU A matrix. #' @param inv_UU A matrix. @@ -9,19 +11,21 @@ #' @param n A numeric. #' @param sparse Logical. #' @return A matrix. -#' @seealso \code{\link[Matrix]{expm}}, \code{link[mcglm]{mc_dexp_gold}} and -#' \code{link[mcglm]{mc_dexpm}}. -#' @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. +#' @seealso \code{\link[Matrix]{expm}}, \code{link[mcglm]{mc_dexp_gold}} +#' and \code{link[mcglm]{mc_dexpm}}. +#' @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. -mc_derivative_expm <- function(dU, UU, inv_UU, Q, n = dim(UU)[1], sparse = FALSE) { - H = inv_UU %*% dU %*% UU +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) diag(P) <- diag(H) * exp(Q) P[upper.tri(P)] <- H[upper.tri(H)] * c(dist(exp(Q))/dist(Q)) P[is.na(P)] <- 0 P <- forceSymmetric(P) - D_Omega = Matrix(UU %*% P %*% inv_UU, sparse = FALSE) + 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 fee51bd54f0a834afde05b67d1c45807e757f4a0..70f3e0f1a4e0f779dd85c4ae6d0b13e54548aad7 100644 --- a/R/mc_derivative_sigma_beta.R +++ b/R/mc_derivative_sigma_beta.R @@ -1,29 +1,37 @@ -#' Derivatives of V^{1/2} with respect to beta. +#' @title Derivatives of V^{1/2} with respect to beta. +#' @author Wagner Hugo Bonat #' -#' @description Compute the derivatives of V^{1/2} matrix with respect to the regression -#' parameters beta. +#' @description Compute the derivatives of \eqn{V^{1/2}} matrix with +#' respect to the regression parameters beta. +#' #' @param D A matrix. #' @param D_V_sqrt_mu A matrix. #' @param Omega A matrix. #' @param V_sqrt A matrix. #' @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. +#' @return A list of matrices, containg the derivatives of \eqn{V^{1/2}} +#' with respect to the regression parameters. -mc_derivative_sigma_beta <- function(D, D_V_sqrt_mu, Omega, V_sqrt, variance) { +mc_derivative_sigma_beta <- function(D, D_V_sqrt_mu, Omega, V_sqrt, + variance) { n_beta <- dim(D)[2] n_obs <- dim(D)[1] output <- list() - if (variance == "power" | variance == "binomialP" | variance == "binomialPQ") { + if (variance == "power" | variance == "binomialP" | + variance == "binomialPQ") { for (i in 1:n_beta) { D_V_sqrt_beta <- Diagonal(n_obs, D_V_sqrt_mu * D[, i]) - output[[i]] <- mc_sandwich_power(middle = Omega, bord1 = V_sqrt, bord2 = D_V_sqrt_beta) + output[[i]] <- + mc_sandwich_power(middle = Omega, + bord1 = V_sqrt, bord2 = D_V_sqrt_beta) } } if (variance == "poisson_tweedie") { for (i in 1:n_beta) { D_V_sqrt_beta <- Diagonal(n_obs, D_V_sqrt_mu * D[, i]) - output[[i]] <- Diagonal(n_obs, D[, i]) + mc_sandwich_power(middle = Omega, bord1 = V_sqrt, bord2 = D_V_sqrt_beta) + output[[i]] <- Diagonal(n_obs, D[, i]) + + mc_sandwich_power(middle = Omega, bord1 = V_sqrt, + bord2 = D_V_sqrt_beta) } } return(output) diff --git a/R/mc_dexp_gold.R b/R/mc_dexp_gold.R index 77ff1a0b10abc5f7984e0615a882fbc8211b5a38..9c0d8a49ac7c6be5c1caadf38f4ddf8296943dde 100644 --- a/R/mc_dexp_gold.R +++ b/R/mc_dexp_gold.R @@ -1,11 +1,14 @@ -#' Exponential-matrix and its derivatives +#' @title Exponential-matrix and its derivatives +#' @author Wagner Hugo Bonat #' -#' @description Given a matrix \eqn{M} and its derivative \eqn{dM} the function \code{dexp_gold} -#' returns the exponential-matrix \eqn{expm(M)} and its derivative. This function is based on -#' the \code{\link[Matrix]{expm}} function. It is not really used in the package, but I keep this -#' function to test my own implementation based on eigen values decomposition. +#' @description Given a matrix \eqn{M} and its derivative \eqn{dM} the +#' function \code{dexp_gold} returns the exponential-matrix +#' \eqn{expm(M)} and its derivative. This function is based on the +#' \code{\link[Matrix]{expm}} function. It is not really used in the +#' package, but I keep this function to test my own implementation +#' based on eigen values decomposition. #' -#' @param M A matrix. +#' @param M A matrix. #' @param dM A matrix. #' @return A list with two elements: \eqn{expm(M)} and its derivatives. #' @seealso \code{\link[Matrix]{expm}}, \code{\link[base]{eigen}}. @@ -16,10 +19,10 @@ #' @export mc_dexp_gold <- function(M, dM) { - N = dim(M) - AM = rbind(cbind(M, matrix(0, N[1], N[2])), cbind(dM, M)) - PSI = expm(AM) - F = PSI[1:N[1], 1:N[1]] - dF = PSI[c(N[1] + 1):c(2 * N[1]), 1:N[1]] + N <- dim(M) + AM <- rbind(cbind(M, matrix(0, N[1], N[2])), cbind(dM, M)) + PSI <- expm(AM) + 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)) }