Primeiro commit

.Rbuildignore

+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

+*.aux
+*.log
+*.out
+.Rhistory
+.Rproj.user
+*~
+.#*
+#*
+*.Rout
+*.RData

DESCRIPTION

+Package: mcglm
+Type: Package
+Title: Fit multivariate covariance generalized linear models
+Version: 0.0.0.9000
+Date: 2015-07-06
+Authors@R: person("Wagner Hugo", "Bonat", email = "wbonat@ufpr.br",
+  role = c("aut", "cre"))
+Description: This package fit multivariate covariance generalized linear models (McGLM).
+Depends: R (>= 3.2.1), Matrix, assertthat
+License: GPL-2
+LazyData: TRUE
+imports: Matrix
+Suggests: testthat

NAMESPACE

+exportPattern("^[[:alpha:]]+")
+import(Matrix)

R/fit_mcglm.R

+#' Chaser and reciprocal likelihood algorithms.
+#'
+#' @description This function implements the two manin algorithms used to fitting McGLMs. The chaser
+#' and the reciprocal likelihood algorithms. In general the chaser algorithm is faster than the
+#' reciprocal likelihood and should be preferred.
+#'
+#' @param list_initial A list of initial values for regression and covariance parameters.
+#' @param list_link A list of link function names.
+#' @param list_variance A list of variance function names.
+#' @param list_covariance A list of covariance function names.
+#' @param list_X A list of design matrices.
+#' @param list_Z A list of knowm matrices to used on the matrix linear predictor.
+#' @param list_offset A list of offset values.
+#' @param list_Ntrial A list of number of trials, useful only when analysis binomial data.
+#' @param list_power_fixed A list of logicals indicating if the power parameters should be estimated
+#' or not.
+#' @param list_sparse A list of logicals indicating if the matrices should be set up as sparse
+#' matrices. This argument is useful only when using exponential-matrix covariance link function.
+#' In the other cases the algorithm detects automatically if the matrix should be sparse or not.
+#' @param y_vec A vector of the response variables stacked.
+#' @param correct A logical indicating if the algorithm will use the correction term or not.
+#' @param max_iter Maximum number of iterations.
+#' @param tol A numeric spcyfing the tolerance.
+#' @param method A string specyfing the method used to fit the models (\code{chaser} or \code{rc}).
+#' @param tunning A numeric value in general near zero for the rc method and near 1 for the chaser
+#' method. This argument control the step-length.
+#' @param verbose A logical if TRUE print the values of the covariance parameters used on each iteration.
+#' @return A list with estimated regression and covariance parameters.
+
+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 = 1e-03, method = "rc", tunning = 0,
+                      verbose){
+  ## Transformation from list to vector
+  parametros <- mc_list2vec(list_initial, list_power_fixed)
+  n_resp <- length(list_initial$regression)
+  if(n_resp == 1){parametros$cov_ini <- parametros$cov_ini[-1]}
+  ## Getting information about the number of parameters
+  inf <- mc_getInformation(list_initial, list_power_fixed, n_resp = n_resp)
+  ## Creating a matrix to sote all values used in the fitting step
+  solucao_beta <- matrix(NA, max_iter,length(parametros$beta_ini))
+  solucao_cov <- matrix(NA, max_iter, length(parametros$cov_ini))
+  ## Setting the initial values
+  solucao_beta[1,] <- parametros$beta_ini
+  solucao_cov[1,] <- parametros$cov_ini
+  beta_ini <- parametros$beta_ini
+  cov_ini <- parametros$cov_ini
+  for(i in 2:max_iter){
+    ## Step 1 - Quasi-score function
+    # Step 1.1 - Computing the mean structure
+  mu_list <- Map(mc_link_function, beta = list_initial$regression, offset = list_offset, X = list_X,
+                 link = list_link)
+  mu_vec <- do.call(c,lapply(mu_list, function(x)x$mu))
+  D <- bdiag(lapply(mu_list, function(x)x$D))
+    # Step 1.2 - Computing the inverse of C matrix. I should improve this step.
+    # I have to write a new function to compute only C or inv_C to be more efficient in this step.
+  Cfeatures <- mc_build_C(list_mu = mu_list, list_Ntrial = list_Ntrial, rho = list_initial$rho,
+                          list_tau = list_initial$tau, list_power = list_initial$power,
+                          list_Z = list_Z, list_sparse = list_sparse, list_variance = list_variance,
+                          list_covariance = list_covariance, list_power_fixed = list_power_fixed,
+                          compute_C = FALSE, compute_derivative_beta = FALSE,
+                          compute_derivative_cov = FALSE)
+    # Step 1.3 - Update the regression parameters
+  beta_temp <- mc_quasi_score(D = D, inv_C = Cfeatures$inv_C, y_vec = y_vec, mu_vec = mu_vec)
+  solucao_beta[i,] <- as.numeric(beta_ini - solve(beta_temp$Sensitivity, beta_temp$Score))
+  list_initial <- updatedBeta(list_initial, solucao_beta[i,], information = inf, n_resp = n_resp)
+    # Step 1.4 - Updated the mean structure to use in the Pearson estimating function step.
+  mu_list <- Map(mc_link_function, beta = list_initial$regression, offset = list_offset, X = list_X,
+                 link = list_link)
+  mu_vec <- do.call(c,lapply(mu_list, function(x)x$mu))
+  D <- bdiag(lapply(mu_list, function(x)x$D))
+    # Step 2 - Updating the covariance parameters
+  Cfeatures <- mc_build_C(list_mu = mu_list, list_Ntrial = list_Ntrial, rho = list_initial$rho,
+                          list_tau = list_initial$tau, list_power = list_initial$power,
+                          list_Z = list_Z,
+                          list_sparse = list_sparse, list_variance = list_variance,
+                          list_covariance = list_covariance,
+                          list_power_fixed = list_power_fixed, compute_C = FALSE,
+                          compute_derivative_beta = FALSE)
+    # Step 2.1 - Using beta(i+1)
+  beta_temp2 <- mc_quasi_score(D = D, inv_C = Cfeatures$inv_C, y_vec = y_vec, mu_vec = mu_vec)
+  inv_J_beta <- solve(beta_temp2$Sensitivity)
+  if(method == "chaser"){
+    cov_temp <- mc_pearson(y_vec = y_vec, mu_vec = mu_vec, Cfeatures = Cfeatures,
+                           inv_J_beta = inv_J_beta, D = D, correct = correct,
+                           compute_variability = TRUE)
+    step <- tunning*solve(cov_temp$Sensitivity, cov_temp$Score)
+  }
+  if(method == "rc"){
+    cov_temp <- mc_pearson(y_vec = y_vec, mu_vec = mu_vec, Cfeatures = Cfeatures,
+                           inv_J_beta = inv_J_beta, D = D, correct = correct,
+                           compute_variability = TRUE)
+    step <- solve(tunning*cov_temp$Score%*%t(cov_temp$Score)%*%
+                    solve(cov_temp$Variability)%*%cov_temp$Sensitivity +
+                    cov_temp$Sensitivity)%*%cov_temp$Score
+  }
+  ## Step 2.2 - Upedating the covariance parameters
+  cov_next <- as.numeric(cov_ini - step)
+  list_initial <- updatedCov(list_initial = list_initial, list_power_fixed = list_power_fixed,
+                             covariance = cov_next,
+                             information = inf, n_resp = n_resp)
+  ## print the parameters values
+  if(verbose == TRUE){print(round(cov_next,4))}
+  ## Step 2.3 - Updating the initial values for the next step
+  beta_ini <- solucao_beta[i,]
+  cov_ini <- cov_next
+  solucao_cov[i,] <- cov_next
+  ## Checking the convergence
+  tolera <- abs(c(solucao_beta[i,],solucao_cov[i,])  - c(solucao_beta[i-1,],solucao_cov[i-1,]))
+  if(all(tolera < tol) == TRUE)break
+  }
+  D_C_beta <- mc_build_C(list_mu = mu_list, list_Ntrial = list_Ntrial, rho = list_initial$rho,
+                         list_tau = list_initial$tau, list_power = list_initial$power,
+                         list_Z = list_Z,
+                         list_sparse = list_sparse, list_variance = list_variance,
+                         list_covariance = list_covariance,
+                         list_power_fixed = list_power_fixed, compute_C = TRUE,
+                         compute_derivative_beta = TRUE,
+                         compute_derivative_cov = TRUE)
+  Product_beta <- lapply(D_C_beta$D_C_beta, mc_multiply, bord2 = D_C_beta$inv_C)
+  S_cov_beta <- mc_cross_sensitivity(Product_cov = cov_temp$Extra,
+                                     Product_beta = Product_beta,
+                                     n_beta_effective = length(beta_temp$Score))
+  res <- y_vec - mu_vec
+  V_cov_beta <- mc_cross_variability(Product_cov = cov_temp$Extra, inv_C = D_C_beta$inv_C,
+                                     res = res, D = D)
+  p1 <- rbind(beta_temp2$Variability, t(V_cov_beta))
+  p2 <- rbind(V_cov_beta, cov_temp$Variability)
+  joint_variability <- cbind(p1,p2)
+  inv_S_beta <- inv_J_beta
+  inv_S_cov <- solve(cov_temp$Sensitivity)
+  mat0 <- Matrix(0, ncol = dim(S_cov_beta)[1], nrow = dim(S_cov_beta)[2])
+  cross_term <- -inv_S_cov%*%S_cov_beta%*%inv_S_beta
+  p1 <- rbind(inv_S_beta,cross_term)
+  p2 <- rbind(mat0, inv_S_cov)
+  joint_inv_sensitivity <- cbind(p1,p2)
+  VarCov <- joint_inv_sensitivity%*%joint_variability%*%t(joint_inv_sensitivity)
+  output <- list("IterationRegression" = solucao_beta, "IterationCovariance" = solucao_cov,
+                 "Regression" = beta_ini, "Covariance" = cov_ini, "vcov" = VarCov,
+                 "fitted" = mu_vec, "residuals" = res, "inv_C" = D_C_beta$inv_C,
+                 "C" = D_C_beta$C, "Information" = inf)
+  return(output)
+}
+
+
+
+
+

R/hello.R

+# Hello, world!
+#
+# This is an example function named 'hello' 
+# which prints 'Hello, world!'.
+#
+# You can learn more about package authoring with RStudio at:
+#
+#   http://r-pkgs.had.co.nz/
+#
+# Some useful keyboard shortcuts for package authoring:
+#
+#   Build and Reload Package:  'Ctrl + Shift + B'
+#   Check Package:             'Ctrl + Shift + E'
+#   Test Package:              'Ctrl + Shift + T'
+
+hello <- function() {
+  print("Hello, world!")
+}

R/mc_auxiliar.R

+#' Matrix product in sandwich form
+#'
+#' @description The function \code{mc_sandwich} is just an auxiliar function to compute product matrix
+#' in the sandwich form bord1*middle*bord2. An special case appears when computing the derivative of
+#' 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.
+#'
+#' @param middle A matrix.
+#' @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)
+## Auxiliar function to multiply matrices
+mc_sandwich <- function(middle, bord1, bord2){bord1%*%middle%*%bord2}
+
+#' @rdname mc_sandwich
+mc_sandwich_negative <- function(middle, bord1, bord2){-bord1%*%middle%*%bord2}
+
+#' @rdname mc_sandwich
+mc_sandwich_power <- function(middle, bord1, bord2){
+  temp1 <- mc_sandwich(middle = middle, bord1 = bord1, bord2 = bord2)
+  return(temp1 + t(temp1))
+}
+
+#' @rdname mc_sandwich
+mc_sandwich_cholesky <- function(bord1, middle, bord2){
+  p1 <- bord1%*%middle%*%bord2
+  return(p1 + t(p1))
+}
+
+#' @rdname mc_sandwich
+mc_multiply <- function(bord1, bord2){
+  return(bord2%*%bord1)
+}
+
+#' @rdname mc_sandwich
+mc_multiply2 <- function(bord1, bord2){
+  return(bord1%*%bord2)
+}

R/mc_build_C.R

+#' Build the joint covariance matrix
+#'
+#' @description This function builds the joint variance-covariance matrix using the Generalized
+#' Kronecker product and its derivatives with respect to rho, power and tau parameters.
+#'
+#'@param list_mu A list with values of the mean.
+#'@param list_Ntrial A list with the number of trials. Usefull only for binomial responses.
+#'@param list_tau A list with values for the tau parameters.
+#'@param list_power A list with values for the power parameters.
+#'@param list_Z A list of matrix to be used in the matrix linear predictor.
+#'@param list_sparse A list with Logical.
+#'@param list_variance A list specifying the variance function to be used for each response variable.
+#'@param list_covariance A list specifying the covariance function to be used for each response variable.
+#'@param list_power_fixed A list of Logical specifying if the power parameters are fixed or not.
+#'
+#'@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.
+
+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]])
+  n_rho <- n_resp*(n_resp - 1)/2
+  if(n_resp != 1){assert_that(n_rho == length(rho))}
+  list_Sigma_within = suppressWarnings(Map(mc_build_sigma, mu =list_mu, Ntrial = list_Ntrial,
+                                           tau = list_tau,
+                                           power = list_power,Z = list_Z, sparse = list_sparse,
+                                           variance = list_variance, covariance = list_covariance,
+                                           power_fixed = list_power_fixed,
+                                           compute_derivative_beta = compute_derivative_beta))
+  list_Sigma_chol <- lapply(list_Sigma_within, function(x)x$Sigma_chol)
+  list_Sigma_inv_chol <- lapply(list_Sigma_within, function(x)x$Sigma_chol_inv)
+  Sigma_between <- mc_build_sigma_between(rho = rho, n_resp = n_resp)
+  II <- Diagonal(n_obs, 1)
+  nucleo <- kronecker(Sigma_between$Sigmab,II)
+  Bdiag_chol_Sigma_within <- bdiag(list_Sigma_chol)
+  t_Bdiag_chol_Sigma_within <- t(Bdiag_chol_Sigma_within)
+  Bdiag_inv_chol_Sigma <- bdiag(list_Sigma_inv_chol)
+  inv_C <- Bdiag_inv_chol_Sigma%*%kronecker(solve(Sigma_between$Sigmab),II)%*%t(Bdiag_inv_chol_Sigma)
+  output <- list("inv_C" = inv_C)
+  if(compute_derivative_cov == TRUE){
+  list_D_Sigma <- lapply(list_Sigma_within, function(x)x$D_Sigma)
+  ## Derivatives of C with respect to power and tau parameters
+  list_D_chol_Sigma <- Map(mc_derivative_cholesky, derivada = list_D_Sigma,
+                           inv_chol_Sigma = list_Sigma_inv_chol, chol_Sigma = list_Sigma_chol)
+  mat_zero <- mc_build_bdiag(n_resp = n_resp, n_obs = n_obs)
+  Bdiag_D_chol_Sigma <- mapply(mc_transform_list_bdiag, list_mat = list_D_chol_Sigma,
+                               response_number = 1:n_resp, MoreArgs = list(mat_zero = mat_zero))
+  Bdiag_D_chol_Sigma <- do.call(c, Bdiag_D_chol_Sigma)
+  D_C = lapply(Bdiag_D_chol_Sigma, mc_sandwich_cholesky, middle = nucleo,
+                                         bord2 = t_Bdiag_chol_Sigma_within)
+  ## Finish the derivatives with respect to power and tau parameters
+  if(n_resp > 1){
+    D_C_rho <- mc_derivative_C_rho(D_Sigmab = Sigma_between$D_Sigmab,
+                                   Bdiag_chol_Sigma_within = Bdiag_chol_Sigma_within,
+                                   t_Bdiag_chol_Sigma_within = t_Bdiag_chol_Sigma_within,
+                                   II = II)
+    D_C <- c(D_C_rho, D_C)
+  }
+  output$D_C <- D_C
+  }
+  if(compute_C == TRUE){
+    C = t_Bdiag_chol_Sigma_within%*%kronecker(Sigma_between$Sigmab,II)%*%Bdiag_chol_Sigma_within
+    output$C <- C
+  }
+  if(compute_derivative_beta == TRUE){
+    list_D_Sigma_beta <- lapply(list_Sigma_within, function(x)x$D_Sigma_beta)
+    list_D_chol_Sigma_beta <- Map(mc_derivative_cholesky, derivada = list_D_Sigma_beta,
+                                  inv_chol_Sigma = list_Sigma_inv_chol, chol_Sigma = list_Sigma_chol)
+    mat_zero <- mc_build_bdiag(n_resp = n_resp, n_obs = n_obs)
+    Bdiag_D_chol_Sigma_beta <- mapply(mc_transform_list_bdiag, list_mat = list_D_chol_Sigma_beta,
+                                      response_number = 1:n_resp, MoreArgs = list(mat_zero = mat_zero))
+    Bdiag_D_chol_Sigma_beta <- do.call(c, Bdiag_D_chol_Sigma_beta)
+    D_C_beta = lapply(Bdiag_D_chol_Sigma_beta, mc_sandwich_cholesky, middle = nucleo,
+                      bord2 = t_Bdiag_chol_Sigma_within)
+    output$D_C_beta <- D_C_beta
+  }
+  return(output)
+}
+
+
+
+
+
+

R/mc_build_bdiag.R

+#' Build a block-diagonal matrix of zeros.
+#'
+#' @description Build a block-diagonal matrix of zeros. Such functions is used when computing
+#' the derivatives of the Cholesky decomposition of C.
+#'
+#' @param n_resp A numeric specifyng the number of response variables.
+#' @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.
+
+mc_build_bdiag <- function(n_resp, n_obs){
+  list_zero <- list()
+  for(i in 1:n_resp){list_zero[[i]] <- Matrix(0, n_obs, n_obs, sparse = TRUE)}
+  return(list_zero)
+}

R/mc_build_omega.R

+#' Build omega matrix
+#'
+#' @description This function build \eqn{\Omega} matrix according the covariance link function.
+#'
+#' @param tau A vector
+#' @param Z A list of matrices.
+#' @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}.
+
+mc_build_omega <- function(tau, Z, covariance_link, sparse = FALSE){
+  if(covariance_link == "identity"){
+    Omega <- mc_matrix_linear_predictor(tau = tau, Z = Z)
+    output <- list("Omega" = Omega, "D_Omega" = Z)
+  }
+  if(covariance_link == "expm"){
+    U <- mc_matrix_linear_predictor(tau = tau, Z = Z)
+    temp <- mc_expm(U = U, inverse = FALSE, sparse = sparse)
+    D_Omega <- lapply(Z, mc_derivative_expm, UU = temp$UU, inv_UU = temp$inv_UU,
+                      Q = temp$Q, sparse = sparse)
+    output <- list("Omega" = forceSymmetric(temp$Omega), "D_Omega" = D_Omega)
+    }
+  if(covariance_link == "inverse"){
+    inv_Omega <- mc_matrix_linear_predictor(tau = tau, Z = Z)
+    output <- list("inv_Omega" = inv_Omega, "D_inv_Omega" = Z)
+  }
+  return(output)
+}
+

R/mc_build_sigma.R

+#' Build variance-covariance 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 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 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.
+#'@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) {
+  if(variance == "constant") {
+    if(covariance == "identity" | covariance == "expm") {
+      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,
+                     "D_Sigma" = Omega$D_Omega)
+      }
+    if(covariance == "inverse") {
+      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)
+    }
+  }
+
+  if(variance == "tweedie" | variance == "binomialP" | variance == "binomialPQ"){
+    if(variance == "tweedie"){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)
+      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)
+      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)
+        }
+        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)
+        }
+      }
+      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)
+        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_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)
+      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)
+        }
+        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)
+        }
+      }
+      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)
+      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)
+      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)
+      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)
+      }
+      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)
+        output$D_Sigma_beta <-D_Sigma_beta
+        }
+    }
+    if(covariance == "inverse"){
+      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)
+      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 <- c(D_Sigma_p, 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)
+      output$D_Sigma_beta <- D_Sigma_beta
+    }
+    }
+    }
+  return(output)
+}
+

R/mc_build_sigmab.R

+#' Build the correlation matrix between response variables
+#'
+#' @description This function builds the correlation matrix between response variable, its inverse and
+#' derivatives.
+#' @param rho A numeric vector.
+#' @param n_resp A numeric.
+#' @param inverse Logical
+#' @return A list with sigmab and its derivatives with respect to rho.
+
+mc_build_sigma_between <- function(rho, n_resp, inverse = FALSE){
+  output <- list("Sigmab" = 1, "D_Sigmab" = 1)
+  if(n_resp > 1){
+  Sigmab <- Diagonal(n_resp, 1)
+  Sigmab[lower.tri(Sigmab)] <- rho
+  Sigmab <- forceSymmetric(t(Sigmab))
+  D_Sigmab <- mc_derivative_sigma_between(n_resp = n_resp)
+  if(inverse == FALSE){
+    output <- list("Sigmab" = Sigmab, "D_Sigmab" = D_Sigmab)
+  }
+  if(inverse == TRUE){
+    inv_Sigmab <- solve(Sigmab)
+    D_inv_Sigmab <- lapply(D_Sigmab, mc_sandwich_negative, bord1 = inv_Sigmab, bord2 = inv_Sigmab)
+    output <- list("inv_Sigmab" = inv_Sigmab, "D_inv_Sigmab" = D_inv_Sigmab)
+  }
+  }
+  return(output)
+}
+
+#' @rdname mc_build_sigma_between
+mc_derivative_sigma_between <- function(n_resp){
+  position <- combn(n_resp,2)
+  list.Derivative <- list()
+  n_par <- n_resp*(n_resp-1)/2
+  for(i in 1:n_par){
+    Derivative <- Matrix(0, ncol = n_resp, nrow = n_resp)
+    Derivative[position[1,i],position[2,i]] <- Derivative[position[2,i],position[1,i]] <- 1
+    list.Derivative[i][[1]] <- Derivative}
+  return(list.Derivative)
+}
+
+

R/mc_coef.R

+#' Extract model coefficients for mcglm class
+#'
+#' coef.mcglm is a function which extracts model coefficients from objects of mcglm class.
+#' @param object An object of mcglm class.
+#' @param response A numeric or vector specyfing for which response variables the coefficients
+#' should be returned.
+#' @param type A string or string vector specyfing which coefficients should be returned.
+#' Options are "beta", "tau", "power", "tau" and "correlation".
+#' @return A data.frame with estimates, parameters names, response number and parameters type.
+#' @exportMethod
+coef.mcglm <- function(object, std.error = FALSE, response = c(NA,1:length(object$beta_names)),
+                       type = c("beta","tau","power","correlation")){
+  n_resp <- length(object$beta_names)
+  cod_beta <- list()
+  cod_power <- list()
+  cod_tau <- list()
+  type_beta <- list()
+  type_power <- list()
+  type_tau <- list()
+  resp_beta <- list()
+  resp_power <- list()
+  resp_tau <- list()
+  for(i in response){
+    cod_beta[[i]] <- paste(paste("beta", i, sep = ""), 0:c(object$Information$n_betas[[i]]-1), sep = "")
+    type_beta[[i]] <- rep("beta", length(cod_beta[[i]]))
+    resp_beta[[i]] <- rep(response[i], length(cod_beta[[i]]))
+    cod_power[[i]] <- paste(paste("power", i, sep = ""), 1:object$Information$n_power[[i]], sep = "")
+    type_power[[i]] <- rep("power", length(cod_power[[i]]))
+    resp_power[[i]] <- rep(response[i], length(cod_power[[i]]))
+    cod_tau[[i]] <- paste(paste("tau", i, sep = ""), 1:object$Information$n_tau[[i]], sep = "")
+    type_tau[[i]] <- rep("tau", length(cod_tau[[i]]))
+    resp_tau[[i]] <- rep(response[i], length(cod_tau[[i]]))
+  }
+  rho_names <- c()
+  if(n_resp != 1){
+    combination <- combn(n_resp,2)
+    for(i in 1:dim(combination)[2]){
+      rho_names[i] <- paste(paste("rho", combination[1,i], sep = ""), combination[2,i], sep = "")
+    }
+  }
+  type_rho <- rep("correlation", length(rho_names))
+  resp_rho <- rep("NA", length(rho_names))
+  cod <- c(do.call(c,cod_beta),rho_names,
+           do.call(c,Map(c,cod_power, cod_tau)))
+  type_cod <- c(do.call(c,type_beta),type_rho,
+            do.call(c,Map(c,type_power, type_tau)))
+  response_cod <- c(do.call(c,resp_beta),resp_rho,
+                    do.call(c,Map(c,resp_power, resp_tau)))
+
+  Estimates <- c(object$Regression, object$Covariance)
+  coef_temp <- data.frame("Estimates" = Estimates, "Parameters" = cod,
+                          "Type" = type, "Response" = response_cod)
+  if(std.error == TRUE){
+    coef_temp <- data.frame("Estimates" = Estimates, "Std.error" = sqrt(diag(object$vcov)),
+                            "Parameters" = cod,
+                            "Type" = type_cod, "Response" = response_cod)
+  }
+
+  output <- coef_temp[which(coef_temp$Response %in% response & coef_temp$Type %in% type),]
+  return(output)
+}
+
+
+

R/mc_confint.mcglm.R

+#' Confidence Intervals for mcglm
+#'
+#' @description Computes confidence intervals for parameters in a fitted mcglm.
+#'
+#' @param object a fitted mcglm object.
+#' @level the confidence level required.
+#' @return A data.frame with confidence intervals, parameters names, response number and parameters type.
+
+confint.mcglm <- function(object, level = 0.95){
+  n_resp <- length(object$beta_names)
+  cod_beta <- list()
+  cod_power <- list()
+  cod_tau <- list()
+  for(i in 1:n_resp){
+    cod_beta[[i]] <- paste(paste("beta", i, sep = ""), 0:c(object$Information$n_betas[[i]]-1), sep = "")
+    cod_power[[i]] <- paste(paste("power", i, sep = ""), 1:object$Information$n_power[[i]], sep = "")
+    cod_tau[[i]] <- paste(paste("tau", i, sep = ""), 1:object$Information$n_tau[[i]], sep = "")
+  }
+  rho_names <- c()
+  if(n_resp != 1){
+    combination <- combn(n_resp,2)
+    for(i in 1:dim(combination)[2]){
+      rho_names[i] <- paste(paste("rho", combination[1,i], sep = ""), combination[2,i], sep = "")
+    }
+  }
+  cod <- c(do.call(c,cod_beta),rho_names,
+           do.call(c,Map(c,cod_power, cod_tau)))
+  parameters <- c(object$Regression,object$Covariance)
+  std.errors <- sqrt(diag(object$vcov))
+  a <- (1 - level)/2
+  a <- c(a, 1 - a)
+  fac <- qnorm(a)
+  ci <- parameters + std.errors%o%fac
+  colnames(ci) <- paste(format(a,2),"%", sep = "")
+  rownames(ci) <- cod
+  return(ci)
+}

R/mc_core_cross_variability.R

+#' Cross variability matrix
+#'
+#' @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)
+  saida <- (res%*%W%*%res)%*%(t(res)%*%A)
+  return(as.numeric(saida))
+}

R/mc_core_pearson.R

+#' Core of the Pearson estimating function.
+#' @description Core of the Pearson estimating function.
+#'
+#' @param product A matrix.
+#' @param inv_C A matrix.
+#' @param res A vector of residuals.
+#' @return A vector
+#' @details It is an internal function.
+
+mc_core_pearson <- function(product, inv_C, res){
+  output <- t(res)%*%product%*%(inv_C%*%res) - sum(diag(product))
+  return(as.numeric(output))
+}

R/mc_correction.R

+#' Pearson correction term
+#'
+#' @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 \code{[mcglm]{mc_link_function}}.
+#' @param inv_C A matrix. In general the output of the \code{[mcglm]{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)
+  return(unlist(output))
+}

R/mc_cross_sensitivity.R

+#' Cross-sensitivity
+#'
+#' @description Compute the cross-sensitivit 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.
+#' @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)){
+  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)
+  }
+  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]])
+    }
+  }
+  }
+  return(cross_sensitivity)
+}
+

R/mc_cross_variability.R

+#' Compute the cross-variability matrix
+#'
+#' @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.
+mc_cross_variability <- function(Product_cov, inv_C, res, D){
+  Wlist <- lapply(Product_cov, mc_multiply2, bord2 = 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]], r = res)
+    }
+  }
+  return(cross_variability)
+}

R/mc_derivative_C_rho.R

+#' Derivative of C with respect to 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.
+
+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)
+}
+
+
+
+