Metafor rma.mv function for multivariate mixed models

Full link to source code: https://github.com/cran/metafor/blob/master/R/rma.mv.r

# fixed/random/mixed-effects multivariate/multilevel model with:
#    - possibly one or multiple random intercepts (sigma2) with potentially known correlation matrices
#    - possibly correlated random effects for arms/groups/levels within studies (tau2 and rho for 1st term, gamma2 and phi for 2nd term)
# model also allows for correlated sampling errors via non-diagonal V matrix

# V      = variance-covariance matrix of the sampling errors
# sigma2 = (preset) value(s) for the variance of the random intercept(s)
# tau2   = (preset) value(s) for the variance of the random effects
# rho    = (preset) value(s) for the correlation(s) between random effects
# gamma2 = (preset) value(s) for the variance of the random effects
# phi    = (preset) value(s) for the correlation(s) between random effects

# structures when there is an '~ inner | outer' term in the random argument:
# - CS   (compound symmetry)
# - HCS  (heteroscedastic compound symmetry)
# - UN   (general positive-definite matrix with no structure)
# - UNR  (general positive-definite correlation matrix with a single tau2/gamma2 value)
# - AR   (AR1 structure with a single tau2/gamma2 value and autocorrelation rho/phi)
# - HAR  (heteroscedastic AR1 structure with multiple tau2/gamma2 values and autocorrelation rho/phi)
# - CAR  (continuous time AR1 structure)
# - ID   (same as CS but with rho/phi=0)
# - DIAG (same as HCS but with rho/phi=0)
# - SPEXP/SPGAU/SPLIN/SPRAT/SPSPH (spatial structures: exponential, gaussian, linear, rational quadratic, spherical)
# - GEN (general positive-definite matrix for an arbitrary number of predictors)
# - PHYBM/PHYPL/PHYPD (phylogenetic structures: Brownian motion, Pagel's lambda, Pagel's delta)`

`   ### in case user specifies v (instead of V), verbose is set to v, which is non-sensical
   ### - if v is set to the name of a variable in 'data', it won't be found; can check for
   ###   this with try() and inherits(verbose, "try-error")
   ### - if v is set to vi or var (or anything else that might be interpreted as a function),
   ###   then can catch this by checking if verbose is a function

   verbose <- try(verbose, silent=TRUE)

   if (inherits(verbose, "try-error") || is.function(verbose) || length(verbose) != 1L || !(is.logical(verbose) || is.numeric(verbose)))
      stop(mstyle$stop("Argument 'verbose' must be a scalar (logical or numeric/integer)."))

   ### set options(warn=1) if verbose > 2

   if (verbose > 2) {
      opwarn <- options(warn=1)
      on.exit(options(warn=opwarn$warn), add=TRUE)
   }

[ ......... ]

Specifying V and W:


  ### some checks on V (and turn V into a diagonal matrix if it is a column/row vector)

   if (is.null(V))
      stop(mstyle$stop("Must specify 'V' argument."))

   ### catch cases where 'V' is the utils::vi() function

   if (identical(V, utils::vi))
      stop(mstyle$stop("Variable specified for 'V' argument cannot be found."))

   if (is.list(V) && !is.data.frame(V)) {

      ### list elements may be data frames (or scalars), so coerce to matrices

      V <- lapply(V, as.matrix)

      ### check that all elements are square

      if (any(!sapply(V, .is.square)))
         stop(mstyle$stop("All list elements in 'V' must be square matrices."))

      ### turn list into block-diagonal (sparse) matrix

      if (sparse) {
         V <- bdiag(V)
      } else {
         V <- bldiag(V)
      }

   }

   ### check if user constrained V to 0 (can skip a lot of the steps below then)

   if ((.is.vector(V) && length(V) == 1L && V == 0) || (.is.vector(V) && length(V) == k && !anyNA(V) && all(V == 0))) {
      V0 <- TRUE
   } else {
      V0 <- FALSE
   }

   ### turn V into a diagonal matrix if it is a column/row vector
   ### note: if V is a scalar (e.g., V=0), then this will turn V into a kxk
   ### matrix with the value of V along the diagonal

   if (V0 || .is.vector(V) || nrow(V) == 1L || ncol(V) == 1L) {
      if (sparse) {
         V <- Diagonal(k, as.vector(V))
      } else {
         V <- diag(as.vector(V), nrow=k, ncol=k)
      }
   }

   ### turn V into a matrix if it is a data frame

   if (is.data.frame(V))
      V <- as.matrix(V)

   ### remove row and column names (important for isSymmetric() function)
   ### (but only do this if V has row/column names to avoid making an unnecessary copy)

   if (!is.null(dimnames(V)))
      V <- unname(V)

   ### check whether V is square and symmetric (can skip when V0)

   if (!V0 && !.is.square(V))
      stop(mstyle$stop("'V' must be a square matrix."))

   if (!V0 && !isSymmetric(V)) ### note: copy of V is made when doing this
      stop(mstyle$stop("'V' must be a symmetric matrix."))

   ### check length of yi and V

   if (nrow(V) != k)
      stop(mstyle$stop(paste0("Length of 'yi' (", k, ") and length/dimensions of 'V' (", nrow(V), ") is not the same.")))

   ### force V to be sparse when sparse=TRUE (and V is not yet sparse)

   if (sparse && inherits(V, "matrix"))
      V <- Matrix(V, sparse=TRUE)

   ### check if V is numeric (but only for 'regular' matrices, since this is always FALSE for sparse matrices)

   if (inherits(V, "matrix") && !is.numeric(V))
      stop(mstyle$stop("The object/variable specified for the 'V' argument is not numeric."))

### process W if it was specified

   if (!is.null(W)) {

      ### turn W into a diagonal matrix if it is a column/row vector
      ### in general, turn W into A (arbitrary weight matrix)

      if (.is.vector(W) || nrow(W) == 1L || ncol(W) == 1L) {

         W <- as.vector(W)

         ### allow easy setting of W to a single value

         if (length(W) == 1L)
            W <- rep(W, k)

         A <- diag(W, nrow=length(W), ncol=length(W))

      } else {

         A <- W

      }

coraliewilliams / glmmTMB

Metafor rma.mv function for multivariate mixed models #3

Background of rma.mv function:

Specifying V and W: