diff --git a/R/mc_build_sigma.R b/R/mc_build_sigma.R
index 3e858d92821e8605673989574f4cf881220365e7..7d3e282586d6a3d107475ab2aa5471fa2b32a818 100644
--- a/R/mc_build_sigma.R
+++ b/R/mc_build_sigma.R
@@ -1,151 +1,259 @@
-#' Build variance-covariance matrix
+#' @title Build variance-covariance matrix
+#' @author Wagner Hugo Bonat
#'
-#'@description This function builds a variance-covariance matrix, based on the variance function and
-#' omega matrix.
+#' @description This function builds a variance-covariance matrix, based
+#' on the variance function and omega matrix.
#'
-#'@param mu A numeric vector. In general the output from \code{\link{mc_link_function}}.
-#'@param Ntrial A numeric vector, or NULL or a numeric specifing the number of trials in the binomial
-#'experiment. It is usefull only when using variance = binomialP or binomialPQ. In the other cases
-#'it will be ignored.
+#'@param mu A numeric vector. In general the output from
+#' \code{\link{mc_link_function}}.
+#'@param Ntrial A numeric vector, or NULL or a numeric specifing the
+#' number of trials in the binomial experiment. It is usefull only
+#' when using variance = binomialP or binomialPQ. In the other cases
+#' it will be ignored.
#'@param tau A numeric vector.
-#'@param power A numeric or numeric vector. It should be one number for all variance functions except
-#'binomialPQ, in that case the argument specifies both p and q.
+#'@param power A numeric or numeric vector. It should be one number for
+#' all variance functions except binomialPQ, in that case the
+#' argument specifies both p and q.
#'@param Z A list of matrices.
#'@param sparse Logical.
-#'@param variance String specifing the variance function: constant, tweedie, poisson_tweedie,
-#'binomialP or binomialPQ.
-#'@param covariance String specifing the covariance function: identity, inverse or expm.
-#'@param power_fixed Logical if the power parameter is fixed at initial value (TRUE). In the case
-#'power_fixed = FALSE the power parameter will be estimated.
-#'@param compute_derivative_beta Logical. Compute or not the derivative with respect to regression parameters.
-#'@return A list with the Cholesky decomposition of \eqn{\Sigma}, \eqn{\Sigma^{-1}} and the derivative
-#'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}}.
+#'@param variance String specifing the variance function: constant,
+#' tweedie, poisson_tweedie, binomialP or binomialPQ.
+#'@param covariance String specifing the covariance function: identity,
+#' inverse or expm.
+#'@param power_fixed Logical if the power parameter is fixed at initial
+#' value (TRUE). In the case power_fixed = FALSE the power parameter
+#' will be estimated.
+#'@param compute_derivative_beta Logical. Compute or not the derivative
+#' with respect to regression parameters.
+#'@return A list with the Cholesky decomposition of \eqn{\Sigma},
+#' \eqn{\Sigma^{-1}} and the derivative 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}}.
-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)
+ Omega <- mc_build_omega(tau = tau, Z = Z,
+ covariance_link = covariance,
+ sparse = sparse)
chol_Sigma <- chol(Omega$Omega)
inv_chol_Sigma <- solve(chol_Sigma)
- output <- list(Sigma_chol = chol_Sigma,
- Sigma_chol_inv = inv_chol_Sigma,
+ output <- list(Sigma_chol = chol_Sigma,
+ Sigma_chol_inv = inv_chol_Sigma,
D_Sigma = Omega$D_Omega)
}
if (covariance == "inverse") {
- inv_Sigma <- mc_build_omega(tau = tau, Z = Z, covariance_link = "inverse", sparse = sparse)
+ inv_Sigma <- mc_build_omega(tau = tau, Z = Z,
+ covariance_link = "inverse",
+ sparse = sparse)
chol_inv_Sigma <- chol(inv_Sigma$inv_Omega)
chol_Sigma <- solve(chol_inv_Sigma)
- Sigma <- (chol_Sigma) %*% t(chol_Sigma) # Because a compute the inverse of chol_inv_Omega
- D_Sigma <- lapply(inv_Sigma$D_inv_Omega, mc_sandwich_negative, bord1 = Sigma, bord2 = Sigma)
- output <- list(Sigma_chol = t(chol_Sigma), Sigma_chol_inv = t(chol_inv_Sigma), D_Sigma = D_Sigma)
+ ## Because a compute the inverse of chol_inv_Omega
+ Sigma <- (chol_Sigma) %*% t(chol_Sigma)
+ D_Sigma <- lapply(inv_Sigma$D_inv_Omega,
+ mc_sandwich_negative,
+ bord1 = Sigma, bord2 = Sigma)
+ output <- list(Sigma_chol = t(chol_Sigma),
+ Sigma_chol_inv = t(chol_inv_Sigma),
+ D_Sigma = D_Sigma)
}
}
-
- if (variance == "tweedie" | variance == "binomialP" | variance == "binomialPQ") {
+
+ if (variance == "tweedie" | variance == "binomialP" |
+ variance == "binomialPQ") {
if (variance == "tweedie") {
- variance = "power"
+ variance <- "power"
}
if (covariance == "identity" | covariance == "expm") {
- Omega <- mc_build_omega(tau = tau, Z = Z, covariance_link = covariance, sparse = sparse)
- V_sqrt <- mc_variance_function(mu = mu$mu, power = power, Ntrial = Ntrial, variance = variance, inverse = FALSE,
- derivative_power = !power_fixed, derivative_mu = compute_derivative_beta)
- Sigma <- forceSymmetric(V_sqrt$V_sqrt %*% Omega$Omega %*% V_sqrt$V_sqrt)
+ Omega <- mc_build_omega(tau = tau, Z = Z,
+ covariance_link = covariance,
+ sparse = sparse)
+ V_sqrt <- mc_variance_function(
+ mu = mu$mu, power = power,
+ Ntrial = Ntrial,
+ variance = variance,
+ inverse = FALSE,
+ derivative_power = !power_fixed,
+ derivative_mu = compute_derivative_beta)
+ Sigma <- forceSymmetric(V_sqrt$V_sqrt %*% Omega$Omega %*%
+ V_sqrt$V_sqrt)
chol_Sigma <- chol(Sigma)
inv_chol_Sigma <- solve(chol_Sigma)
- D_Sigma <- lapply(Omega$D_Omega, mc_sandwich, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$V_sqrt)
+ D_Sigma <- lapply(Omega$D_Omega, mc_sandwich,
+ bord1 = V_sqrt$V_sqrt,
+ bord2 = V_sqrt$V_sqrt)
if (power_fixed == FALSE) {
if (variance == "power" | variance == "binomialP") {
- D_Sigma_power <- mc_sandwich_power(middle = Omega$Omega, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$D_V_sqrt_p)
- D_Sigma <- c(D_Sigma_power = D_Sigma_power, D_Sigma_tau = D_Sigma)
+ D_Sigma_power <- mc_sandwich_power(
+ middle = Omega$Omega, bord1 = V_sqrt$V_sqrt,
+ bord2 = V_sqrt$D_V_sqrt_p)
+ D_Sigma <- c(D_Sigma_power = D_Sigma_power,
+ D_Sigma_tau = D_Sigma)
}
if (variance == "binomialPQ") {
- D_Sigma_p <- mc_sandwich_power(middle = Omega$Omega, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$D_V_sqrt_p)
- D_Sigma_q <- mc_sandwich_power(middle = Omega$Omega, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$D_V_sqrt_q)
- D_Sigma <- c(D_Sigma_p, D_Sigma_q, D_Sigma)
+ D_Sigma_p <- mc_sandwich_power(
+ middle = Omega$Omega,
+ bord1 = V_sqrt$V_sqrt,
+ bord2 = V_sqrt$D_V_sqrt_p)
+ D_Sigma_q <- mc_sandwich_power(
+ middle = Omega$Omega,
+ bord1 = V_sqrt$V_sqrt,
+ bord2 = V_sqrt$D_V_sqrt_q)
+ D_Sigma <- c(D_Sigma_p, D_Sigma_q, D_Sigma)
}
}
- output <- list(Sigma_chol = chol_Sigma, Sigma_chol_inv = inv_chol_Sigma, D_Sigma = D_Sigma)
+ output <- list(Sigma_chol = chol_Sigma,
+ Sigma_chol_inv = inv_chol_Sigma,
+ D_Sigma = D_Sigma)
if (compute_derivative_beta == TRUE) {
- D_Sigma_beta <- mc_derivative_sigma_beta(D = mu$D, D_V_sqrt_mu = V_sqrt$D_V_sqrt_mu, Omega = Omega$Omega,
- V_sqrt = V_sqrt$V_sqrt, variance = variance)
+ D_Sigma_beta <- mc_derivative_sigma_beta(
+ D = mu$D, D_V_sqrt_mu = V_sqrt$D_V_sqrt_mu,
+ Omega = Omega$Omega, V_sqrt = V_sqrt$V_sqrt,
+ variance = variance)
output$D_Sigma_beta <- D_Sigma_beta
}
}
if (covariance == "inverse") {
- inv_Omega <- mc_build_omega(tau = tau, Z = Z, covariance_link = "inverse", sparse = sparse)
- V_inv_sqrt <- mc_variance_function(mu = mu$mu, power = power, Ntrial = Ntrial, variance = variance, inverse = TRUE,
- derivative_power = !power_fixed, derivative_mu = compute_derivative_beta)
- inv_Sigma <- forceSymmetric(V_inv_sqrt$V_inv_sqrt %*% inv_Omega$inv_Omega %*% V_inv_sqrt$V_inv_sqrt)
+ inv_Omega <- mc_build_omega(tau = tau, Z = Z,
+ covariance_link = "inverse",
+ sparse = sparse)
+ V_inv_sqrt <- mc_variance_function(
+ mu = mu$mu, power = power, Ntrial = Ntrial,
+ variance = variance, inverse = TRUE,
+ derivative_power = !power_fixed,
+ derivative_mu = compute_derivative_beta)
+ inv_Sigma <- forceSymmetric(V_inv_sqrt$V_inv_sqrt %*%
+ inv_Omega$inv_Omega %*%
+ V_inv_sqrt$V_inv_sqrt)
inv_chol_Sigma <- chol(inv_Sigma)
chol_Sigma <- solve(inv_chol_Sigma)
Sigma <- chol_Sigma %*% t(chol_Sigma)
- D_inv_Sigma <- lapply(inv_Omega$D_inv_Omega, mc_sandwich, bord1 = V_inv_sqrt$V_inv_sqrt, bord2 = V_inv_sqrt$V_inv_sqrt)
- D_Sigma <- lapply(D_inv_Sigma, mc_sandwich_negative, bord1 = Sigma, bord2 = Sigma)
+ D_inv_Sigma <- lapply(inv_Omega$D_inv_Omega, mc_sandwich,
+ bord1 = V_inv_sqrt$V_inv_sqrt,
+ bord2 = V_inv_sqrt$V_inv_sqrt)
+ D_Sigma <- lapply(D_inv_Sigma, mc_sandwich_negative,
+ bord1 = Sigma, bord2 = Sigma)
if (power_fixed == FALSE) {
if (variance == "power" | variance == "binomialP") {
- D_Omega_p <- mc_sandwich_power(middle = inv_Omega$inv_Omega, bord1 = V_inv_sqrt$V_inv_sqrt, bord2 = V_inv_sqrt$D_V_inv_sqrt_power)
- D_Sigma_p <- mc_sandwich_negative(middle = D_Omega_p, bord1 = Sigma, bord2 = Sigma)
- D_Sigma <- c(D_Sigma_p, D_Sigma)
+ D_Omega_p <- mc_sandwich_power(
+ middle = inv_Omega$inv_Omega,
+ bord1 = V_inv_sqrt$V_inv_sqrt,
+ bord2 = V_inv_sqrt$D_V_inv_sqrt_power)
+ D_Sigma_p <- mc_sandwich_negative(
+ middle = D_Omega_p,
+ bord1 = Sigma, bord2 = Sigma)
+ D_Sigma <- c(D_Sigma_p, D_Sigma)
}
if (variance == "binomialPQ") {
- D_Omega_p <- mc_sandwich_power(middle = inv_Omega$inv_Omega, bord1 = V_inv_sqrt$V_inv_sqrt, bord2 = V_inv_sqrt$D_V_inv_sqrt_p)
- D_Sigma_p <- mc_sandwich_negative(middle = D_Omega_p, bord1 = Sigma, bord2 = Sigma)
- D_Omega_q <- mc_sandwich_power(middle = inv_Omega$inv_Omega, bord1 = V_inv_sqrt$V_inv_sqrt, bord2 = V_inv_sqrt$D_V_inv_sqrt_q)
- D_Sigma_q <- mc_sandwich_negative(middle = D_Omega_q, bord1 = Sigma, bord2 = Sigma)
- D_Sigma <- c(D_Sigma_p, D_Sigma_q, D_Sigma)
+ D_Omega_p <- mc_sandwich_power(
+ middle = inv_Omega$inv_Omega,
+ bord1 = V_inv_sqrt$V_inv_sqrt,
+ bord2 = V_inv_sqrt$D_V_inv_sqrt_p)
+ D_Sigma_p <- mc_sandwich_negative(
+ middle = D_Omega_p,
+ bord1 = Sigma, bord2 = Sigma)
+ D_Omega_q <- mc_sandwich_power(
+ middle = inv_Omega$inv_Omega,
+ bord1 = V_inv_sqrt$V_inv_sqrt,
+ bord2 = V_inv_sqrt$D_V_inv_sqrt_q)
+ D_Sigma_q <- mc_sandwich_negative(
+ middle = D_Omega_q,
+ bord1 = Sigma, bord2 = Sigma)
+ D_Sigma <- c(D_Sigma_p, D_Sigma_q, D_Sigma)
}
}
- output <- list(Sigma_chol = t(chol_Sigma), Sigma_chol_inv = t(inv_chol_Sigma), D_Sigma = D_Sigma)
+ output <- list(Sigma_chol = t(chol_Sigma),
+ Sigma_chol_inv = t(inv_chol_Sigma),
+ D_Sigma = D_Sigma)
if (compute_derivative_beta == TRUE) {
- D_inv_Sigma_beta <- mc_derivative_sigma_beta(D = mu$D, D_V_sqrt_mu = V_inv_sqrt$D_V_inv_sqrt_mu, Omega = inv_Omega$inv_Omega,
- V_sqrt = V_inv_sqrt$V_inv_sqrt, variance = variance)
- D_Sigma_beta <- lapply(D_inv_Sigma_beta, mc_sandwich_negative, bord1 = Sigma, bord2 = Sigma)
+ D_inv_Sigma_beta <- mc_derivative_sigma_beta(
+ D = mu$D,
+ D_V_sqrt_mu = V_inv_sqrt$D_V_inv_sqrt_mu,
+ Omega = inv_Omega$inv_Omega,
+ V_sqrt = V_inv_sqrt$V_inv_sqrt, variance = variance)
+ D_Sigma_beta <- lapply(D_inv_Sigma_beta,
+ mc_sandwich_negative,
+ bord1 = Sigma, bord2 = Sigma)
output$D_Sigma_beta <- D_Sigma_beta
}
}
}
-
+
if (variance == "poisson_tweedie") {
if (covariance == "identity" | covariance == "expm") {
- Omega <- mc_build_omega(tau = tau, Z = Z, covariance_link = covariance, sparse = sparse)
- V_sqrt <- mc_variance_function(mu = mu$mu, power = power, Ntrial = Ntrial, variance = "power", inverse = FALSE,
- derivative_power = !power_fixed, derivative_mu = compute_derivative_beta)
- Sigma <- forceSymmetric(Diagonal(length(mu$mu), mu$mu) + V_sqrt$V_sqrt %*% Omega$Omega %*% V_sqrt$V_sqrt)
+ Omega <- mc_build_omega(tau = tau, Z = Z,
+ covariance_link = covariance,
+ sparse = sparse)
+ V_sqrt <- mc_variance_function(
+ mu = mu$mu, power = power,
+ Ntrial = Ntrial, variance = "power", inverse = FALSE,
+ derivative_power = !power_fixed,
+ derivative_mu = compute_derivative_beta)
+ Sigma <- forceSymmetric(Diagonal(length(mu$mu), mu$mu) +
+ V_sqrt$V_sqrt %*%
+ Omega$Omega %*% V_sqrt$V_sqrt)
chol_Sigma <- chol(Sigma)
inv_chol_Sigma <- solve(chol_Sigma)
- D_Sigma <- lapply(Omega$D_Omega, mc_sandwich, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$V_sqrt)
+ D_Sigma <- lapply(Omega$D_Omega, mc_sandwich,
+ bord1 = V_sqrt$V_sqrt,
+ bord2 = V_sqrt$V_sqrt)
if (power_fixed == FALSE) {
- D_Sigma_power <- mc_sandwich_power(middle = Omega$Omega, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$D_V_sqrt_p)
- D_Sigma <- c(D_Sigma_power = D_Sigma_power, D_Sigma_tau = D_Sigma)
+ D_Sigma_power <- mc_sandwich_power(
+ middle = Omega$Omega,
+ bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$D_V_sqrt_p)
+ D_Sigma <- c(D_Sigma_power = D_Sigma_power,
+ D_Sigma_tau = D_Sigma)
}
- output <- list(Sigma_chol = chol_Sigma, Sigma_chol_inv = inv_chol_Sigma, D_Sigma = D_Sigma)
+ output <- list(Sigma_chol = chol_Sigma,
+ Sigma_chol_inv = inv_chol_Sigma,
+ D_Sigma = D_Sigma)
if (compute_derivative_beta == TRUE) {
- D_Sigma_beta <- mc_derivative_sigma_beta(D = mu$D, D_V_sqrt_mu = V_sqrt$D_V_sqrt_mu, Omega$Omega, V_sqrt = V_sqrt$V_sqrt,
- variance = variance)
+ D_Sigma_beta <- mc_derivative_sigma_beta(
+ D = mu$D,
+ D_V_sqrt_mu = V_sqrt$D_V_sqrt_mu, Omega$Omega,
+ V_sqrt = V_sqrt$V_sqrt, variance = variance)
output$D_Sigma_beta <- D_Sigma_beta
}
}
if (covariance == "inverse") {
- inv_Omega <- mc_build_omega(tau = tau, Z = Z, covariance_link = "inverse", sparse = sparse)
+ inv_Omega <- mc_build_omega(tau = tau, Z = Z,
+ covariance_link = "inverse",
+ sparse = sparse)
Omega <- chol2inv(chol(inv_Omega$inv_Omega))
- V_sqrt <- mc_variance_function(mu = mu$mu, power = power, Ntrial = Ntrial, variance = "power", inverse = FALSE,
- derivative_power = !power_fixed, derivative_mu = compute_derivative_beta)
- D_Omega <- lapply(inv_Omega$D_inv_Omega, mc_sandwich_negative, bord1 = Omega, bord2 = Omega)
- D_Sigma <- lapply(D_Omega, mc_sandwich, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$V_sqrt)
- Sigma <- forceSymmetric(Diagonal(length(mu$mu), mu$mu) + V_sqrt$V_sqrt %*% Omega %*% V_sqrt$V_sqrt)
+ V_sqrt <- mc_variance_function(
+ mu = mu$mu, power = power,
+ Ntrial = Ntrial, variance = "power", inverse = FALSE,
+ derivative_power = !power_fixed,
+ derivative_mu = compute_derivative_beta)
+ D_Omega <- lapply(inv_Omega$D_inv_Omega,
+ mc_sandwich_negative, bord1 = Omega,
+ bord2 = Omega)
+ D_Sigma <- lapply(D_Omega, mc_sandwich,
+ bord1 = V_sqrt$V_sqrt,
+ bord2 = V_sqrt$V_sqrt)
+ Sigma <- forceSymmetric(Diagonal(length(mu$mu), mu$mu) +
+ V_sqrt$V_sqrt %*% Omega %*%
+ V_sqrt$V_sqrt)
chol_Sigma <- chol(Sigma)
inv_chol_Sigma <- solve(chol_Sigma)
if (power_fixed == FALSE) {
- D_Sigma_p <- mc_sandwich_power(middle = Omega, bord1 = V_sqrt$V_sqrt, bord2 = V_sqrt$D_V_sqrt_power)
+ D_Sigma_p <- mc_sandwich_power(
+ middle = Omega,
+ bord1 = V_sqrt$V_sqrt,
+ bord2 = V_sqrt$D_V_sqrt_power)
D_Sigma <- c(D_Sigma_p, D_Sigma)
}
- output <- list(Sigma_chol = chol_Sigma, Sigma_chol_inv = inv_chol_Sigma, D_Sigma = D_Sigma)
+ output <- list(Sigma_chol = chol_Sigma,
+ Sigma_chol_inv = inv_chol_Sigma,
+ D_Sigma = D_Sigma)
if (compute_derivative_beta == TRUE) {
- D_Sigma_beta <- mc_derivative_sigma_beta(D = mu$D, D_V_sqrt_mu = V_sqrt$D_V_sqrt_mu, Omega = Omega, V_sqrt = V_sqrt$V_sqrt,
- variance = variance)
+ D_Sigma_beta <- mc_derivative_sigma_beta(
+ D = mu$D,
+ D_V_sqrt_mu = V_sqrt$D_V_sqrt_mu, Omega = Omega,
+ V_sqrt = V_sqrt$V_sqrt, variance = variance)
output$D_Sigma_beta <- D_Sigma_beta
}
}