Skip to content
Snippets Groups Projects
Commit 567bf462 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Add tags and tidy code an text.

parent 31c47c43
Branches
No related tags found
No related merge requests found
#' Cross variability matrix #' @title Cross variability matrix
#' @author Wagner Hugo Bonat
#' #'
#' @description Compute the cross-covariance matrix between covariance and regression parameters. #' @description Compute the cross-covariance matrix between covariance
#' Equation (11) of Bonat and Jorgensen (2015). #' and regression parameters. Equation (11) of Bonat and Jorgensen
#' (2015).
#' #'
#' @param A A matrix. #' @param A A matrix.
#' @param res A vector of residuals. #' @param res A vector of residuals.
#' @param W A matrix of weights. #' @param W A matrix of weights.
covprod <- function(A, res, W) { covprod <- function(A, res, W) {
res = as.numeric(res) res <- as.numeric(res)
saida <- (res %*% W %*% res) %*% (t(res) %*% A) saida <- (res %*% W %*% res) %*% (t(res) %*% A)
return(as.numeric(saida)) return(as.numeric(saida))
} }
#' Core of the Pearson estimating function. #' @title Core of the Pearson estimating function.
#' @author Wagner Hugo Bonat
#'
#' @description Core of the Pearson estimating function. #' @description Core of the Pearson estimating function.
#' #'
#' @param product A matrix. #' @param product A matrix.
...@@ -8,6 +10,7 @@ ...@@ -8,6 +10,7 @@
#' @details It is an internal function. #' @details It is an internal function.
mc_core_pearson <- function(product, inv_C, res) { 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)) return(as.numeric(output))
} }
#' 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 D_C A list of matrices.
#' @param inv_J_beta A matrix. In general it is computed based on the output of the #' @param inv_J_beta A matrix. In general it is computed based on the
#' \code{[mcglm]{mc_quasi_score}}. #' 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 D A matrix. In general it is the output of the
#' @param inv_C A matrix. In general the output of the \link{mc_build_C}. #' \link{mc_link_function}.
#' @return A vector with the correction terms to be used on the Pearson estimating function. #' @param inv_C A matrix. In general the output of the
#' @details It is an internal function useful inside the fitting algorithm. #' \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) { 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) term1 <- lapply(D_C, mc_sandwich, bord1 = t(D) %*% inv_C,
output <- lapply(term1, function(x, inv_J_beta) sum(x * inv_J_beta), inv_J_beta = inv_J_beta) 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)) return(unlist(output))
} }
#' 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_cov A list of matrices.
#' @param Product_beta A list of matrices. #' @param Product_beta A list of matrices.
#' @param n_beta_effective Numeric. Effective number of regression parameters. #' @param n_beta_effective Numeric. Effective number of regression
#' @return The cross-sensitivity matrix. Equation (10) of Bonat and Jorgensen (2015). #' 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_beta <- length(Product_beta)
n_cov <- length(Product_cov) n_cov <- length(Product_cov)
if (n_beta == 0) { 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) { if (n_beta != 0) {
cross_sensitivity <- Matrix(NA, nrow = n_cov, ncol = n_beta) cross_sensitivity <- Matrix(NA, nrow = n_cov, ncol = n_beta)
for (i in 1:n_cov) { for (i in 1:n_cov) {
for (j in 1:n_beta) { 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]])
} }
} }
} }
......
#' 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 Product_cov A list of matrices.
#' @param inv_C A matrix. #' @param inv_C A matrix.
#' @param res A vector. #' @param res A vector.
#' @param D A matrix. #' @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) { mc_cross_variability <- function(Product_cov, inv_C, res, D) {
Wlist <- lapply(Product_cov, mc_multiply2, bord2 = inv_C) 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_beta <- dim(A)[1]
n_cov <- length(Product_cov) n_cov <- length(Product_cov)
cross_variability <- Matrix(NA, ncol = n_cov, nrow = n_beta) cross_variability <- Matrix(NA, ncol = n_cov, nrow = n_beta)
for (j in 1:n_beta) { for (j in 1:n_beta) {
for (i in 1:n_cov) { 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) return(cross_variability)
......
#' 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 D_Sigmab A matrix.
#'@param Bdiag_chol_Sigma_within A block-diagonal matrix. #'@param Bdiag_chol_Sigma_within A block-diagonal matrix.
#'@param t_Bdiag_chol_Sigma_within A block-diagonal matrix. #'@param t_Bdiag_chol_Sigma_within A block-diagonal matrix.
#'@param II A diagonal matrix. #'@param II A diagonal matrix.
#'@return A 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) { mc_derivative_C_rho <- function(D_Sigmab, Bdiag_chol_Sigma_within,
output <- lapply(D_Sigmab, function(x) { t_Bdiag_chol_Sigma_within, II) {
t_Bdiag_chol_Sigma_within %*% kronecker(x, II) %*% Bdiag_chol_Sigma_within output <- lapply(D_Sigmab,
function(x) {
t_Bdiag_chol_Sigma_within %*%
kronecker(x, II) %*%
Bdiag_chol_Sigma_within
}) })
return(output) return(output)
} }
#' Derivatives of the Cholesky decomposition #' @title Derivatives of the Cholesky decomposition
#' @description This function compute the derivative of the Cholesky decomposition. #' @author Wagner Hugo Bonat
#'
#' @description This function compute the derivative of the Cholesky
#' decomposition.
#' #'
#' @param derivada A matrix. #' @param derivada A matrix.
#' @param inv_chol_Sigma A matrix. #' @param inv_chol_Sigma A matrix.
...@@ -7,7 +10,8 @@ ...@@ -7,7 +10,8 @@
#' @return A list of matrix. #' @return A list of matrix.
#' @details It is an internal function. #' @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) { faux <- function(derivada, inv_chol_Sigma, chol_Sigma) {
t1 <- inv_chol_Sigma %*% derivada %*% t(inv_chol_Sigma) t1 <- inv_chol_Sigma %*% derivada %*% t(inv_chol_Sigma)
t1 <- tril(t1) t1 <- tril(t1)
...@@ -15,6 +19,8 @@ mc_derivative_cholesky <- function(derivada, inv_chol_Sigma, chol_Sigma) { ...@@ -15,6 +19,8 @@ mc_derivative_cholesky <- function(derivada, inv_chol_Sigma, chol_Sigma) {
output <- chol_Sigma %*% t1 output <- chol_Sigma %*% t1
return(output) 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) return(list_D_chol)
} }
#' 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 UU A matrix.
#' @param inv_UU A matrix. #' @param inv_UU A matrix.
...@@ -9,19 +11,21 @@ ...@@ -9,19 +11,21 @@
#' @param n A numeric. #' @param n A numeric.
#' @param sparse Logical. #' @param sparse Logical.
#' @return A matrix. #' @return A matrix.
#' @seealso \code{\link[Matrix]{expm}}, \code{link[mcglm]{mc_dexp_gold}} and #' @seealso \code{\link[Matrix]{expm}}, \code{link[mcglm]{mc_dexp_gold}}
#' \code{link[mcglm]{mc_dexpm}}. #' and \code{link[mcglm]{mc_dexpm}}.
#' @details Many arguments required by this function are provide by the \code{link[mcglm]{mc_dexpm}}. #' @details Many arguments required by this function are provide by the
#' The argument dU is the derivative of the U matrix with respect to the models parameters. It should #' \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. #' be computed by the user.
mc_derivative_expm <- function(dU, UU, inv_UU, Q, n = dim(UU)[1], sparse = FALSE) { mc_derivative_expm <- function(dU, UU, inv_UU, Q, n = dim(UU)[1],
H = inv_UU %*% dU %*% UU sparse = FALSE) {
H <- inv_UU %*% dU %*% UU
P <- Matrix(0, ncol = n, nrow = n) P <- Matrix(0, ncol = n, nrow = n)
diag(P) <- diag(H) * exp(Q) diag(P) <- diag(H) * exp(Q)
P[upper.tri(P)] <- H[upper.tri(H)] * c(dist(exp(Q))/dist(Q)) P[upper.tri(P)] <- H[upper.tri(H)] * c(dist(exp(Q))/dist(Q))
P[is.na(P)] <- 0 P[is.na(P)] <- 0
P <- forceSymmetric(P) P <- forceSymmetric(P)
D_Omega = Matrix(UU %*% P %*% inv_UU, sparse = FALSE) D_Omega <- Matrix(UU %*% P %*% inv_UU, sparse = FALSE)
return(D_Omega) return(D_Omega)
} }
#' 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 \eqn{V^{1/2}} matrix with
#' respect to the regression parameters beta.
#' #'
#' @description Compute the derivatives of V^{1/2} matrix with respect to the regression
#' parameters beta.
#' @param D A matrix. #' @param D A matrix.
#' @param D_V_sqrt_mu A matrix. #' @param D_V_sqrt_mu A matrix.
#' @param Omega A matrix. #' @param Omega A matrix.
#' @param V_sqrt A matrix. #' @param V_sqrt A matrix.
#' @param variance A string specifying the variance function name. #' @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 #' @return A list of matrices, containg the derivatives of \eqn{V^{1/2}}
#' parameters. #' 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_beta <- dim(D)[2]
n_obs <- dim(D)[1] n_obs <- dim(D)[1]
output <- list() output <- list()
if (variance == "power" | variance == "binomialP" | variance == "binomialPQ") { if (variance == "power" | variance == "binomialP" |
variance == "binomialPQ") {
for (i in 1:n_beta) { for (i in 1:n_beta) {
D_V_sqrt_beta <- Diagonal(n_obs, D_V_sqrt_mu * D[, i]) 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") { if (variance == "poisson_tweedie") {
for (i in 1:n_beta) { for (i in 1:n_beta) {
D_V_sqrt_beta <- Diagonal(n_obs, D_V_sqrt_mu * D[, i]) 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) return(output)
......
#' 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} #' @description Given a matrix \eqn{M} and its derivative \eqn{dM} the
#' returns the exponential-matrix \eqn{expm(M)} and its derivative. This function is based on #' function \code{dexp_gold} returns the exponential-matrix
#' the \code{\link[Matrix]{expm}} function. It is not really used in the package, but I keep this #' \eqn{expm(M)} and its derivative. This function is based on the
#' function to test my own implementation based on eigen values decomposition. #' \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. #' @param dM A matrix.
...@@ -16,10 +19,10 @@ ...@@ -16,10 +19,10 @@
#' @export #' @export
mc_dexp_gold <- function(M, dM) { mc_dexp_gold <- function(M, dM) {
N = dim(M) N <- dim(M)
AM = rbind(cbind(M, matrix(0, N[1], N[2])), cbind(dM, M)) AM <- rbind(cbind(M, matrix(0, N[1], N[2])), cbind(dM, M))
PSI = expm(AM) PSI <- expm(AM)
F = PSI[1:N[1], 1:N[1]] F <- PSI[1:N[1], 1:N[1]]
dF = PSI[c(N[1] + 1):c(2 * N[1]), 1:N[1]] dF <- PSI[c(N[1] + 1):c(2 * N[1]), 1:N[1]]
return(list(F, dF)) return(list(F, dF))
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment