R/nlin_djikstra_correction.R

In plssem: Complex Partial Least Squares Structural Equation Modeling

# Modified code, written by Florian Schubert
# REPO: https://github.com/FloSchuberth/cSEM/
# FILE: R/helper_estimators_paths.R
# COMMENT: I couldn't be bothered to write all this my self
#          it's just too much...

classifyConstructs <- function(.terms) {
  ## Split term
  terms_split <- strsplit(.terms, ":")
  
  ## Count instances of each construct name (used for classifying)
  terms_classified <- lapply(terms_split, function(.x) {
    x <- as.data.frame(table(.x), stringAsFactors = FALSE)
    
    ## To save typing
    a <- sum(x$Freq)
    b <- length(unique(x$.x))
    
    ## Do the actual classification --------------------------------------------
    
    if(a > 3) {
      stop("The nonlinear term(s): ", paste0("`", .terms, "`", collapse = ", "), 
           ifelse(length(.terms == 1), " is", " are"), " currently not supported.\n",
           "Please see ?classifyConstructs for a list of supported terms.",
           call. = FALSE)
    } else {
      switch(a,
             "1" = {
               x <- data.frame("Term_class"     = "Single",
                               "Component"      = x$.x,
                               "Component_type" = "Single",
                               "Component_freq" = x$Freq,
                               stringsAsFactors = FALSE)
             },
             "2" = {
               if(b == 1) {
                 x <- data.frame("Term_class"     = "Quadratic",
                                 "Component"      = x$.x,
                                 "Component_type" = "Quadratic",
                                 "Component_freq" = x$Freq,
                                 stringsAsFactors = FALSE)
               } else if (b == 2){
                 x <- data.frame("Term_class"     = "TwInter",
                                 "Component"      = x$.x,
                                 "Component_type" = c("Single", "Single"),
                                 "Component_freq" = x$Freq,
                                 stringsAsFactors = FALSE)
               }
             },
             "3" = {
               if(b == 1) {
                 x <- data.frame("Term_class"     = "Cubic",
                                 "Component"      = x$.x,
                                 "Component_type" = "Cubic",
                                 "Component_freq" = x$Freq,
                                 stringsAsFactors = FALSE)
               } else if (b == 2) {
                 x <- data.frame("Term_class"     = "QuadTwInter",
                                 "Component"      = x$.x,
                                 "Component_type" =
                                   ifelse(x$Freq == 1, "Single", "Quadratic"),
                                 "Component_freq" = x$Freq,
                                 stringsAsFactors = FALSE)
               } else if (b == 3) {
                 
                 x <- data.frame("Term_class"     = "ThrwInter",
                                 "Component"      = x$.x,
                                 "Component_type" =
                                   c("Single", "Single", "Single"),
                                 "Component_freq" = x$Freq,
                                 stringsAsFactors = FALSE)
               }
             }
      ) # END switch
    } # END else
  }) # END lapply
  
  names(terms_classified) <- unlist(.terms)
  terms_classified
}

# Internal: Utility functions for the estimation of nonlinear models
# 
# 
# @param .i Row index
# @param .j Column index
# @param .select_from matrix to select from
# @inheritParams csem_arguments
# 
# @name nonlinear_estimation_utilities
# @rdname nonlinear_estimation_utilities
# @keywords internal
f1 <- function(.i, .j) {

  tab_i <- classifyConstructs(.i)[[1]]
  tab_j <- classifyConstructs(.j)[[1]]

  class_ij <- paste(unique(tab_i$Term_class),
                    unique(tab_j$Term_class), sep = "_")

  ## Rename/reorder joint class so that "class1_class2" = "class2_class1"
  switch (class_ij,
          # Class "Single_x"
          "Quadratic_Single"      = {class_ij <- "Single_Quadratic"},
          "Cubic_Single"          = {class_ij <- "Single_Cubic"},
          "TwInter_Single"        = {class_ij <- "Single_TwInter"},
          "ThrwInter_Single"      = {class_ij <- "Single_ThrwInter"},
          "QuadTwInter_Single"    = {class_ij <- "Single_QuadTwInter"},
          # Class "Quadratic_x"
          "Cubic_Quadratic"       = {class_ij <- "Quadratic_Cubic"},
          "TwInter_Quadratic"     = {class_ij <- "Quadratic_TwInter"},
          "ThrwInter_Quadrtic"    = {class_ij <- "Quadratic_ThrwInter"},
          "QuadTwInter_Quadratic" = {class_ij <- "Quadratic_QuadTwInter"},
          # Class "Cubic_x"
          "TwInter_Cubic"         = {class_ij <- "Cubic_TwInter"},
          "ThrwInter_Cubic"       = {class_ij <- "Cubic_ThrwInter"},
          "QuadTwInter_Cubic"     = {class_ij <- "Cubic_QuadTwInter"},
          # Class "TwInter_x"
          "ThrwInter_TwInter"     = {class_ij <- "TwInter_ThrwInter"},
          "QuadTwInter_TwInter"   = {class_ij <- "TwInter_QuadTwInter"},
          # Class "ThrwInter_x"
          "QuadTwInter_ThrwInter" = {class_ij <- "ThrwInter_QuadTwInter"},
          # Dont change the object if none of the above apply
          # Note: If none of the alternatives apply,
          # switch evaluates the last alternative...not sure anymore if thats
          # true. Needs to be checked!
          "None_of_the_above"     = {class_ij}
  )
  ## Return
  class_ij
}


# @rdname nonlinear_estimation_utilities
f2 <- function(.i, .j, .select_from, .Q, .H) {

  x <- .select_from[.i, .j, drop = FALSE]

  ## Select estimator according to joint class of the ij'th element of
  ## the VCV matrix of the explanatory variables

  args <- list(.i = .i, .j = .j, .Q = .Q, .H = .H)

  switch (x,
          # Class "Single_x"
          "Single_Single"           = {res <- do.call(SingleSingle, args)},
          "Single_Quadratic"        = {res <- do.call(SingleQuadratic, args)},
          "Single_Cubic"            = {res <- do.call(SingleCubic, args)},
          "Single_TwInter"          = {res <- do.call(SingleTwInter, args)},
          "Single_ThrwInter"        = {res <- do.call(SingleThrwInter, args)},
          "Single_QuadTwInter"      = {res <- do.call(SingleQuadTwInter, args)},
          # Class "Quadratic_x"
          "Quadratic_Quadratic"     = {res <- do.call(QuadraticQuadratic, args)},
          "Quadratic_Cubic"         = {res <- do.call(QuadraticCubic, args)},
          "Quadratic_TwInter"       = {res <- do.call(QuadraticTwInter, args)},
          "Quadratic_ThrwInter"     = {res <- do.call(QuadraticThrwInter, args)},
          "Quadratic_QuadTwInter"   = {res <- do.call(QuadraticQuadTwInter, args)},
          # Class "Cubic_x"
          "Cubic_Cubic"             = {res <- do.call(CubicCubic, args)},
          "Cubic_TwInter"           = {res <- do.call(CubicTwInter, args)},
          "Cubic_ThrwInter"         = {res <- do.call(CubicThrwInter, args)},
          "Cubic_QuadTwInter"       = {res <- do.call(CubicQuadTwInter, args)},
          # Class "TwInter_x"
          "TwInter_TwInter"         = {res <- do.call(TwInterTwInter, args)},
          "TwInter_ThrwInter"       = {res <- do.call(TwInterThrwInter, args)},
          "TwInter_QuadTwInter"     = {res <- do.call(TwInterQuadTwInter, args)},
          # Class "ThrwInter_x"
          "ThrwInter_ThrwInter"     = {res <- do.call(ThrwInterThrwInter, args)},
          "ThrwInter_QuadTwInter"   = {res <- do.call(ThrwInterQuadTwInter, args)},
          # Class "QuadTwInter_x"
          "QuadTwInter_QuadTwInter" = {res <- do.call(QuadTwInercQuadTwInter, args)}
  )
  ## Retun result
  res
}


# @rdname nonlinear_estimation_utilities
# Internal: Calculate consistent moments of a nonlinear model 
# 
# Collection of various moment estimators. See [classifyConstructs] for a list of
# possible moments.
# 
# M is the matrix of the sample counterparts (estimates) of the
# left-hand side terms in Equation (21) - (24) \insertCite{Dijkstra2014}{cSEM}.
# The label "M" did not appear in the paper and is only used in the package.
# Similar is suggested by \insertCite{Wall2000;textual}{cSEM} using classical factor scores.
# 
# @param .i Row index
# @param .j Column index
# @inheritParams csem_arguments
# 
# @references
#   \insertAllCited{} 
#     
# @name moments
# @rdname moments
# @keywords internal

SingleSingle <- function(.i, .j, .Q, .H) {
  
  ## Calculate M
  M <- mean(.H[, .i] * .H[, .j])
  
  ## Calculate denominator D
  D <- prod(.Q[c(.i, .j)])
  
  ## Calculate VCV element
  if(.i != .j) { # Covariance (rowindex != columnindex)
    x <- M / D
  } else { # Variance: (rowindex == columnindex)
    x <- 1
  }
  ## Return
  x
}


# @rdname moments
SingleQuadratic <- function(.i, .j, .Q, .H) {
  
  ## Label
  i_split <- unlist(strsplit(.i, ":"))
  j_split <- unlist(strsplit(.j, ":"))
  ij      <- c(i_split, j_split)
  
  ## Calculate M
  M <- mean(matrixStats::rowProds(.H[, ij]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  x <- M / D
  
  ## Return
  x
}
# @rdname moments
SingleCubic <- function(.i, .j, .Q, .H) {
  
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Single" term and which one contains
  ## the "Cubic" (=other) term?
  single_tab <- if(term_class_i == "Single") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Single") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Single" term
  # j := name of the "Cubic" term
  # j_split   := the "Cubic" term split into its components
  # j_single  := the name of the component(s) of the "Cubic" term
  # ij        := all components
  
  i           <- names(single_tab)
  j           <- names(other_tab)
  j_split     <- unlist(strsplit(j, ":"))
  j_single    <- unique(j_split)
  ij          <- c(i, j_split)
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, c(i, j_single)]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(intersect(i, j_split)) == 0) {
    # 1.) The "Single" term does not equal the components of the "Cubic" term
    
    x <- (M1 - 3 * (1 - .Q[j_single]^2) * M2) / D
    
  } else {
    # 2.) The "Single" term is equal to the components of the "Cubic" term
    
    x <- (M1 - 3 * (1 - .Q[j_single]^4)) / D
    
  }
  ## Return
  x
}
# @rdname moments
SingleTwInter <- function(.i, .j, .Q, .H) {
  
  ## Label
  i_split <- unlist(strsplit(.i, ":"))
  j_split <- unlist(strsplit(.j, ":"))
  ij      <- c(i_split, j_split)
  
  ## Calculate M
  M <- mean(matrixStats::rowProds(.H[, ij]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  x <- M / D
  
  ## Return
  x
}
# @rdname moments
SingleThrwInter <- function(.i, .j, .Q, .H) {
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Single" term and which one contains
  ## the "ThrwInter" (=other) term?
  single_tab <- if(term_class_i == "Single") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Single") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Single" term
  # j := name of the "ThrwInter term
  # j_split     := the "ThrwInter" term split into its components
  # ij          := all components
  # ij_match    := the name of the component of the "ThrwInter" term that matches
  #                the "Single" term
  
  i        <- names(single_tab)
  j        <- names(other_tab)
  j_split  <- unlist(strsplit(j, ":"))
  ij       <- c(i, j_split)
  ij_match <- intersect(i, j_split)
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, c(setdiff(j_split, ij_match))]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) Non of the components match
    
    x <- M1 / D
    
  } else {
    # 2.) The (component of) "Single" matches one of the components of the
    #     "ThrwInter" term
    
    x <- (M1 - (1 - .Q[ij_match]^2) * M2) / D
    
  }
  ## Return
  x
}
# @rdname moments
SingleQuadTwInter <- function(.i, .j, .Q, .H) {
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Single" term and which one contains
  ## the "QuadTwInter" (=other) term?
  single_tab <- if(term_class_i == "Single") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Single") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Single" term
  # j := name of the "QuadTwInter" term
  # ij          := all components
  # j_split     := the "QuadTwInter" term split into its components
  # j_single    := name of the "Single" component of the "QuadTwInter" term
  # j_quadratic := name of the "Quadratic" component of the "QuadTwInter" term
  
  i           <- names(single_tab)
  j           <- names(other_tab)
  j_split     <- unlist(strsplit(j, ":"))
  ij          <- c(i, j_split)
  j_single    <- other_tab[[1]][other_tab[[1]]$Component_type == "Single", ]$Component
  j_quadratic <- other_tab[[1]][other_tab[[1]]$Component_type == "Quadratic", ]$Component
  
  ## Calculation ---------------------------------------------------------------
  ## Calculate M
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, c(i, j_single)]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(intersect(i, j_split)) == 0) {
    # 1.) None of the components match
    
    x <- (M1 - (1 - .Q[j_quadratic]^2) * M2) / D
    
  } else if (i == j_single) {
    # 2.) The "Single" term matches the "Single" component of the
    #     "QuadTwInter" term
    
    x <- (M1 - 1) / D + 1
    
  } else {
    # 3.) The "Single" term matches the "Quadratic" component of the
    #     "QuadTwInter" term
    
    x <- (M1 - 3 * (1 - .Q[j_quadratic]^2) * M2) / D
    
  }
  ## Return
  x
}
# @rdname moments
QuadraticQuadratic <- function(.i, .j, .Q, .H) {
  
  ## Label
  i_split  <- unlist(strsplit(.i, ":"))
  j_split  <- unlist(strsplit(.j, ":"))
  i_single <- unique(i_split)
  ij       <- c(i_split, j_split)
  
  ## Calculate M
  M <- mean(matrixStats::rowProds(.H[, ij]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(.i != .j) {
    # 1.) None of the components match
    
    x <- (M - 1) / D
    
  } else {
    # 2.) All components match (terms are identical)
    
    x <- (M - 3 * (1 - .Q[i_single]^4)) / D - 1
    
  }
  ## Return
  x
}
# @rdname moments
QuadraticCubic <- function(.i, .j, .Q, .H) {
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Quadratic" term and which
  ## one contains the "Cubic" (=other) term?
  quad_tab   <- if(term_class_i == "Quadratic") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Quadratic") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Quadratic" term
  # j := name of the "Cubic" term
  # i_split   := the "Quadratic" term split into its components
  # j_split   := the "Cubic" term split into its components
  # i_single  := the name of the component(s) of the "Quadratic" term
  # j_single  := the name of the component(s) of the "Cubic" term
  # ij        := all components
  
  i         <- names(quad_tab)
  j         <- names(other_tab)
  i_split   <- unlist(strsplit(i, ":"))
  j_split   <- unlist(strsplit(j, ":"))
  i_single  <- unique(i_split)
  j_single  <- unique(j_split)
  ij        <- c(i_split, j_split)
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, c(i_split, j_single)]))
  M3 <- mean(matrixStats::rowProds(.H[, j_split]))
  
  ## Calculate denominator (D)
  D0 <- prod(.Q[j_split])
  D1 <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(i_single != j_single) {
    # 1.) Non of the components match
    
    x <- (M1 - 3*(1 - .Q[j_single]^2) * M2 - (1 - .Q[i_single]^2) * M3) / D1 -
      M3 / D0
    
  } else {
    # 2.) The component(s) of the "Quadratic" term matche those of the
    #    "Cubic" term
    
    x <- (M1 - 10 * (1 - .Q[j_single]^2) * M3) / D1 - M3/D0
    
  }
  ## Return
  x
}
# @rdname moments
QuadraticTwInter <- function(.i, .j, .Q, .H) {
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Quadratic" term and which
  ## one contains the "TwInter" (=other) term?
  quad_tab   <- if(term_class_i == "Quadratic") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Quadratic") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Quadratic" term
  # j := name of the "TwInter" term
  # i_split   := the "Quadratic" term split into its components
  # j_split   := the "TwInter" term split into its components
  # i_single  := the name of the component(s) of the "Quadratic" term
  # ij        := all components
  
  i         <- names(quad_tab)
  j         <- names(other_tab)
  i_split   <- unlist(strsplit(i, ":"))
  j_split   <- unlist(strsplit(j, ":"))
  i_single  <- unique(i_split)
  ij        <- c(i_split, j_split)
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, j_split]))
  
  ## Calculate denominator (D)
  D0 <- prod(.Q[j_split])
  D  <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(intersect(i_split, j_split)) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - (1 - .Q[i_single]^2) * M2) / D - M2 / D0
    
  } else {
    # 2.) The component(s) of the "Quadratic" term match one of the
    #     component(s) of the "TwInter" term
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M2) / D - M2 / D0
    
  }
  ## Return
  x
}
# @rdname moments
QuadraticThrwInter <- function(.i, .j, .Q, .H) {
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Quadratic" term and which one contains
  ## the "ThrwInter" (=other) term?
  quad_tab   <- if(term_class_i == "Quadratic") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Quadratic") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Quadratic" term
  # j := name of the "ThrwInter" term
  # i_split   := the "Quadratic" term split into its components
  # j_split   := the "ThrwInter" term split into its components
  # i_single  := the name of the component(s) of the "Quadratic" term
  # ij        := all components
  
  i         <- names(quad_tab)
  j         <- names(other_tab)
  i_split   <- unlist(strsplit(i, ":"))
  j_split   <- unlist(strsplit(j, ":"))
  i_single  <- unique(i_split)
  ij        <- c(i_split, j_split)
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, j_split]))
  
  ## Calculate denominator (D)
  D0 <- prod(.Q[j_split])
  D  <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(intersect(i_split, j_split)) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - (1 - .Q[i_single]^2) * M2) / D - M2 / D0
    
  } else {
    # 2.) The component(s) of the "Quadratic" match one of the component(s)
    #     of the "ThrwInter" term
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M2) / D - M2 / D0
    
  }
  ## Return
  x
}
# @rdname moments
QuadraticQuadTwInter <- function(.i, .j, .Q, .H) {
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Quadratic" term and which
  ## one contains the "TwInter" (=other) term?
  quad_tab   <- if(term_class_i == "Quadratic") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Quadratic") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Quadratic" term
  # j := name of the "QuadTwInter" term
  # i_split     := the "Quadratic" term split into its components
  # j_split     := the "QuadTwInter" term split into its components
  # i_single    := the name of the component(s) of the "Quadratic" term
  
  # j_single    := name of the "Single" component of the "QuadTwInter" term
  # j_quadratic := name of the "Quadratic" component of the "QuadTwInter" term
  
  i           <- names(quad_tab)
  j           <- names(other_tab)
  i_split     <- unlist(strsplit(i, ":"))
  j_split     <- unlist(strsplit(j, ":"))
  i_single    <- unique(i_split)
  ij          <- c(i_split, j_split)
  j_single    <- other_tab[[1]][other_tab[[1]]$Component_type == "Single", ]$Component
  j_quadratic <- other_tab[[1]][other_tab[[1]]$Component_type == "Quadratic", ]$Component
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, j_split]))
  M3 <- mean(matrixStats::rowProds(.H[, c(i_split, j_single)]))
  
  ## Calculate denominator (D)
  D0 <- prod(.Q[j_split])
  D  <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(intersect(i_split, j_split)) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - (1 - .Q[i_single]^2) * M2 - (1 - .Q[j_quadratic]^2) * M3) / D -
      M2 / D0
    
  } else if (i_single == j_single) {
    # 2.) The component(s) of the "Quadratic" term match the "Single" component
    #     of the "QuadTwInter" term
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M2 - (1 - .Q[j_quadratic]^2) * M3) / D -
      M2 / D0
    
  } else {
    # 3.) The component(s) of the "Quadratic" term match the "Quadratic"
    #     component of the "QuadTwInter" term
    
    x <- (M1 - 6 * (1 - .Q[i_single]^2) * M2) / D -
      M2 / D0
    
  }
  ## Return
  x
}
# @rdname moments
CubicCubic <- function(.i, .j, .Q, .H) {
  
  ## Label
  # i_split  := the term .i split into its components
  # j_split  := the term .j split into its components
  # i_single := the name of the component of the "Cubic" term .i
  # j_single := the name of the component of the "Cubic" term .j
  # ij       := all components
  
  i_split  <- unlist(strsplit(.i, ":"))
  j_split  <- unlist(strsplit(.j, ":"))
  i_single <- unique(i_split)
  j_single <- unique(j_split)
  ij       <- c(i_split, j_split)
  
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, c(i_split, j_single)]))
  M3  <- mean(matrixStats::rowProds(.H[, c(j_split, i_single)]))
  M4  <- mean(matrixStats::rowProds(.H[, c(i_single, j_single)]))
  
  
  ## Calculate denominator (D)
  D0a <- prod(.Q[i_single])
  D0b <- prod(.Q[j_single])
  D  <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(i_single != j_single) {
    # 1.) Non of the components match
    
    x <- (M1 - 3 * (1 - .Q[j_single]^2) * M2 - 3 * (1 - .Q[i_single]^2) * M3 +
            9 * (1 - .Q[i_single]^2) * (1 - .Q[j_single]) * M4) / D -
      (M0a / D0a) * (M0b / D0b)
    
  } else {
    # 2.) All components match (terms are identical)
    
    x <- (M1 - 15 * (1 - .Q[i_single]^2) * M2 + 15 * .Q[i_single]^6 -
            45 * .Q[i_single]^2 + 30) / D -
      (M0a / D0a)^2
    
  }
  ## Return
  x
}
# @rdname moments
CubicTwInter <- function(.i, .j, .Q, .H) {
  ### Prepartion ---------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Cubic" term and which one contains
  ## the "TwInter" (=other) term?
  cubic_tab  <- if(term_class_i == "Cubic") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Cubic") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Cubic" term
  # j := name of the "TwInter" term
  # i_split   := the "Cubic" term split into its components
  # j_split   := the "TwInter" term split into its components
  # i_single  := the name of the component(s) of the "Cubic" term
  # ij        := all components
  
  i         <- names(cubic_tab)
  j         <- names(other_tab)
  i_split   <- unlist(strsplit(i, ":"))
  j_split   <- unlist(strsplit(j, ":"))
  i_single  <- unique(i_split)
  ij        <- c(i_split, j_split)
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, c(i_single, j_split)]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(intersect(i_split, j_split)) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M2) / D - (M0a * M0b) / D
    
  } else {
    # 2.) The component(s) of the "Cubic" term match one of the component(s)
    #     of the "TwInter" term
    
    x <- (M1 - 6 * (1 - .Q[i_single]^2) * M2) / D - (M0a * M0b) / D
    
  }
  ## Return
  x
}
# @rdname moments
CubicThrwInter <- function(.i, .j, .Q, .H) {
  ### Prepartion ---------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Cubic" term and which one contains
  ## the "ThrwInter" (=other) term?
  cubic_tab  <- if(term_class_i == "Cubic") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Cubic") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Cubic" term
  # j := name of the "ThrwInter" term
  # i_split     := the "Cubic" term split into its components
  # j_split     := the "ThrwInter" term split into its components
  # i_single    := the name of the component(s) of the "Cubic" term
  # ij          := all components
  # ij_match    := the name of the component of the "ThrwInter" term that matches
  #                the component of the "Cubic" term
  
  i         <- names(cubic_tab)
  j         <- names(other_tab)
  i_split   <- unlist(strsplit(i, ":"))
  j_split   <- unlist(strsplit(j, ":"))
  i_single  <- unique(i_split)
  ij        <- c(i_split, j_split)
  ij_match  <- intersect(i_single, j_split)
  
  ### Prepartion ---------------------------------------------------------------
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, c(i_single, j_split)]))
  M3  <- mean(matrixStats::rowProds(.H[, setdiff(j_split, ij_match)]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M2) / D - (M0a * M0b) / D
    
  } else {
    # 2.) The component(s) of the "Cubic" term match one of the component(s)
    #     of the "ThrwInter" term
    
    x <- (M1 - 6 * (1 - .Q[i_single]^2) * M2 + 3 *
            (1 - .Q[i_single]^2)^2 * M3) / D - (M0a * M0b) / D
    
  }
  ## Return
  x
}
# @rdname moments
CubicQuadTwInter <- function(.i, .j, .Q, .H) {
  ### Prepartion ---------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "Cubic" term and which one contains
  ## the "QuadTwInter" (=other) term?
  cubic_tab  <- if(term_class_i == "Cubic") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "Cubic") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "Cubic" term
  # j := name of the "QuadTwInter" term
  # i_split     := the "Cubic" term split into its components
  # j_split     := the "QuadTwInter" term split into its components
  # i_single    := the name of the component(s) of the "Cubic" term
  # j_single    := name of the "Single" component of the "QuadTwInter" term
  # j_quadratic := name of the "Quadratic" component of the "QuadTwInter" term
  # ij          := all components
  
  i           <- names(cubic_tab)
  j           <- names(other_tab)
  i_split     <- unlist(strsplit(i, ":"))
  j_split     <- unlist(strsplit(j, ":"))
  i_single    <- unique(i_split)
  ij          <- c(i_split, j_split)
  j_single    <- other_tab[[1]][other_tab[[1]]$Component_type == "Single", ]$Component
  j_quadratic <- other_tab[[1]][other_tab[[1]]$Component_type == "Quadratic", ]$Component
  
  ### Prepartion ---------------------------------------------------------------
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, c(i_single, j_split)]))
  M3  <- mean(matrixStats::rowProds(.H[, c(i_split, j_single)]))
  M4  <- mean(matrixStats::rowProds(.H[, c(i_single, j_single)]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(intersect(i_split, j_split)) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M2 - (1 - .Q[j_quadratic]^2) * M3 +
            3 * (1 -  .Q[i_single]^2) * (1 - .Q[j_single]^2) * M4) / D -
      (M0a * M0b) / D
    
  } else if (i_single == j_single) {
    # 2.) The component(s) of the "Cubic" term match the "Single" component
    #     of the "QuadTwInter" term
    
    x <- (M1 - 6 * (1 - .Q[i_single]^2) * M2 - (1 - .Q[j_quadratic]^2) * M3 -
            3 * (1 - .Q[i_single]^2) * (.Q[i_single]^2 * .Q[j_quadratic]^2 +
                                          .Q[j_quadratic]^2 - 2)) / D -
      (M0a * M0b) / D
    
  } else {
    # 3.) The component(s) of the "Cubic" term match the "Quadratic" component
    #     of the "QuadTwInter" term
    
    x <- (M1 - 10 * (1 - .Q[i_single]^2) * M2 + 15 *
            (1 - .Q[i_single]^2)^2 * M4) / D - (M0a * M0b) / D
    
  }
  ## Return
  x
}
# @rdname moments
TwInterTwInter <- function(.i, .j, .Q, .H) {
  ## Label
  # i_split  := the term .i split into its components
  # j_split  := the term .j split into its components
  # i_single := the names of the components of the "TwInter" term .i
  # j_single := the names of the components of the "TwInter" term .j
  # ij       := all components
  # ij_match := the name of the component of .i that match those of .j
  
  i_split  <- unlist(strsplit(.i, ":"))
  j_split  <- unlist(strsplit(.j, ":"))
  i_single <- unique(i_split)
  j_single <- unique(j_split)
  ij       <- c(i_split, j_split)
  ij_match <- intersect(i_split, j_split)
  
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_single]))
  M0b <- mean(matrixStats::rowProds(.H[, j_single]))
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match)]]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) Non of the components match
    
    x <- M1 / D - M0a * M0b / D
    
  } else if (length(ij_match) == 1){
    # 2.) One of the components of .i matches one of .j
    
    x <- (M1 - (1 - .Q[ij_match]^2) * M2) / D - (M0a * M0b) / D
    
  } else {
    # 3.) All components of .i match those of .j
    
    x <- (M1 - 1) / D + 1 - (M0a * M0b) / D
    
  }
  ## Return
  x
}
# @rdname moments
TwInterThrwInter <- function(.i, .j, .Q, .H) {
  
  ## Label
  # i_split  := the term .i split into its components
  # j_split  := the term .j split into its components
  # ij       := all components
  # ij_match := the name of the component of .i that match those of .j
  
  i_split  <- unlist(strsplit(.i, ":"))
  j_split  <- unlist(strsplit(.j, ":"))
  ij       <- c(i_split, j_split)
  ij_match <- intersect(i_split, j_split)
  
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match)], drop = FALSE]))
  M3  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match[1])], drop = FALSE]))
  M4  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match[2])], drop = FALSE]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) Non of the components match
    
    x <- M1 / D - (M0a * M0b) / D
    
  } else if (length(ij_match) == 1){
    # 2.) One of the components of .i matches one of .j
    
    x <- (M1 - (1 - .Q[ij_match]^2) * M2) / D - (M0a * M0b) / D
    
  } else {
    # 3.) All components of .i match with any two of the components of .j
    
    x <- (M1 - (1 - .Q[ij_match[1]]^2) * M3 - (1 - .Q[ij_match[2]]^2) * M4) / D - 
      (M0a * M0b) / D
    
  }
  ## Return x
  x
}
# @rdname moments
TwInterQuadTwInter <- function(.i, .j, .Q, .H) {
  
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "TwInter" term and which
  ## one contains the "QuadTwInter" (=other) term?
  tw_tab    <- if(term_class_i == "TwInter") {tab_i} else {tab_j}
  other_tab <- if(term_class_i != "TwInter") {tab_i} else {tab_j}
  
  ## Label
  # i := name of the "TwInter" term
  # j := name of the "QuadTwInter" term
  # i_split  := the "TwInter" term split into its components
  # j_split  := the "QuadTwInter" term split into its components
  # ij       := all components
  # ij_match := the names of the components of the "TwInter" term that match
  #             those of the "QuadTwInter" term
  # j_single    := name of the "Single" component of the "QuadTwInter" term
  # j_quadratic := name of the "Quadratic" component of the "QuadTwInter" term
  
  i           <- names(tw_tab)
  j           <- names(other_tab)
  i_split     <- unlist(strsplit(j, ":"))
  j_split     <- unlist(strsplit(j, ":"))
  ij          <- c(i_split, j_split)
  ij_match    <- intersect(i_split, j_split)
  j_single    <- other_tab[[1]][other_tab[[1]]$Component_type == "Single", ]$Component
  j_quadratic <- other_tab[[1]][other_tab[[1]]$Component_type == "Quadratic", ]$Component
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% j_quadratic)]]))
  M3 <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match)]]))
  M4 <- mean(matrixStats::rowProds(.H[, c(i_split, j_single)]))
  M5 <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% j_single)]]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) None of the components of the "TwInter" term match those of the
    #     "QuadTwInter" term
    
    x <- (M1 - (1 - .Q[j_quadratic]^2) * M2) / D - (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && ij_match == j_single) {
    # 2.) One of the components of the "TwInter" term matches the "Single"
    #     component of the "QuadTwInter" term
    
    x <- (M1 - (1 - .Q[j_quadratic]^2) * M2 - (1 - .Q[ij_match]^2) * M3) / D -
      (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && ij_match == j_quadratic) {
    # 3.) One of the components of the "TwInter" term matches the "Quadratic"
    #     component of the "QuadTwInter" term
    
    x <- (M1 - 3 * (1 - .Q[ij_match]^2) * M4) / D -
      (M0a * M0b) / D
    
  } else {
    # 4.) One of the components of the "TwInter" term matches the "Single"
    #     and the other matches the "Quadratic" component  of the
    #     "QuadTwInter" term
    
    x <- (M1 - 3* (1 - .Q[j_quadratic]^2) * M4 - (1 - .Q[j_single]^2) * M5) / D -
      (M0a * M0b) / D
    
  }
  ## Return
  x
}
# @rdname moments
ThrwInterThrwInter <- function(.i, .j, .Q, .H){
  ## Label
  # i_split     := the "ThrwInter" term .i split into its components
  # j_split     := the "ThrwInter" term .j split into its components
  # ij          := all terms
  # ij_match    := the name of the component of the "ThrwInter" terms that match
  
  i_split    <- unlist(strsplit(.i, ":"))
  j_split    <- unlist(strsplit(.j, ":"))
  ij         <- c(i_split, j_split)
  ij_match   <- intersect(i_split, j_split)
  
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match)]]))
  M3  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match[1])]]))
  M4  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match[2])]]))
  M5  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match[3])]]))
  
  ## Calculate denominator
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) Non of the terms matches
    
    x <- M1 / D - (M0a * M0b) / D
    
  } else if (length(ij_match) == 1) {
    # 2.) Exactly one of the components in the first "ThrwInter" term
    #     matches one in the second term
    
    x <- (M1 - (1 - .Q[ij_match]^2) * M2) / D - (M0a * M0b) / D
    
  } else if (length(ij_match) == 2) {
    # 3.) Two components of the first "ThrwInter" term match with two components
    #     of the other "ThrwInter" term.
    
    x <- (M1 - (1 - .Q[ij_match[1]]^2) * M3 - (1 - .Q[ij_match[2]]^2) * M4 -
            (1 - .Q[ij_match[1]]^2) * (1 - .Q[ij_match[2]]^2) * M2) / D -
      (M0a * M0b) / D
  } else {
    # 4.) All components of the first "ThrwInter" term match those of the second
    
    x <- (M1 - (1 - .Q[ij_match[1]]^2) * M3 - (1 - .Q[ij_match[2]]^2) * M4 -
            (1 - .Q[ij_match[3]]^2) * M5) / D +
      (.Q[ij_match[1]]^2 * .Q[ij_match[2]]^2 * .Q[ij_match[3]]^2 -
         .Q[ij_match[1]]^2 - .Q[ij_match[2]]^2 - .Q[ij_match[3]]^2 + 2) / D -
      (M0a * M0b) / D
    
  }
  ## Return x
  x
}
# @rdname moments
ThrwInterQuadTwInter <- function(.i, .j, .Q, .H){
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)
  tab_j <- classifyConstructs(.j)
  
  ## Get the classification for term .i
  term_class_i <- unique(tab_i[[.i]]$Term_class)
  
  ## Which classification table contains the "ThrwInter" term and which
  ## one contains the "QuadTwInter" (=other) term?
  thrw_tab   <- if(term_class_i == "ThrwInter") {tab_i} else {tab_j}
  other_tab  <- if(term_class_i != "ThrwInter") {tab_i} else {tab_j}
  
  ## Label
  # i_split     := the "ThrwInter" term .i split into its components
  # j_split     := the "QuadTwInter" term .j split into its components
  # ij          := all component names
  # ij_match    := the name of the component of the two terms that match
  # ij_nomatch  := the name of the component of the two terms that dont match
  # j_single    := name of the "Single" component of the "QuadTwInter" term
  # j_quadratic := name of the "Quadratic" component of the "QuadTwInter" term
  
  i           <- names(thrw_tab)
  j           <- names(other_tab)
  i_split     <- unlist(strsplit(i, ":"))
  j_split     <- unlist(strsplit(j, ":"))
  j_single    <- other_tab[[1]][other_tab[[1]]$Component_type == "Single", ]$Component
  j_quadratic <- other_tab[[1]][other_tab[[1]]$Component_type == "Quadratic", ]$Component
  
  ij          <- c(i_split, j_split)
  ij_match    <- intersect(i_split, j_split)
  ij_nomatch  <- setdiff(ij, ij_match)
  
  ### Calculation ---------------------------------
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1 <- mean(matrixStats::rowProds(.H[, ij]))
  M2 <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% j_quadratic)]]))
  M3 <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% j_single)]]))
  M4 <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% c(j_quadratic, ij_match))]]))
  M5 <- mean(matrixStats::rowProds(.H[, c(i_split, j_single)]))
  M6 <- mean(matrixStats::rowProds(.H[, c(ij_nomatch, j_quadratic)]))
  
  ## Calculate denominator
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - (1 - .Q[j_quadratic]^2) * M2) / D - (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && ij_match == j_single) {
    # 2.) One of the components of the "ThrwInter" term matches the "Single"
    #     component of the "QuadTwInter" term.
    
    x <- (M1 - (1 - .Q[j_quadratic]^2) * M2 - (1 - .Q[j_single]^2) * M3 -
            (1 - .Q[j_single]^2) * (1 - .Q[j_quadratic]^2) * M4) / D -
      (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && ij_match == j_quadratic) {
    # 3.) One of the components of the "ThrwInter" term matches the "Quadratic"
    #     component of the "QuadTwInter" term.
    
    x <- (M1 - 3 * (1 - .Q[j_quadratic]^2) * M5) / D - (M0a * M0b) / D
    
  } else {
    # 4.) One component of the "ThrwInter" term matches the "Single" and another
    #     the "Quadratic" component of the "QuadTwInter" term.
    
    x <- (M1 - 3 * (1 - .Q[j_quadratic]^2) * M5 - (1 - .Q[j_single]^2) * M3 -
            3 * (1 - .Q[j_quadratic]^2) * ( 1 - .Q[j_single]^2) *  M6) / D -
      (M0a * M0b) / D
    
  }
  ## Return x
  x
}
# @rdname moments
QuadTwInercQuadTwInter <- function(.i, .j, .Q, .H) {
  ### Preparation --------------------------------------------------------------
  ## Get the classification table
  tab_i <- classifyConstructs(.i)[[1]]
  tab_j <- classifyConstructs(.j)[[1]]
  
  ## Label
  # i := name of the "Quadratic" term
  # j := name of the "QuadTwInter" term
  # i_split     := the "Quadratic" term split into its components
  # j_split     := the "QuadTwInter" term split into its components
  # i_single    := the name of the component(s) of the "Quadratic" term
  # j_single    := name of the "Single" component of the "QuadTwInter" term
  # j_quadratic := name of the "Quadratic" component of the "QuadTwInter" term
  
  i_split     <- unlist(strsplit(.i, ":"))
  j_split     <- unlist(strsplit(.j, ":"))
  ij          <- c(i_split, j_split)
  ij_match    <- intersect(i_split, j_split)
  
  i_single    <- tab_i[tab_i$Component_type == "Single", ]$Component
  i_quadratic <- tab_i[tab_i$Component_type == "Quadratic", ]$Component
  j_single    <- tab_j[tab_j$Component_type == "Single", ]$Component
  j_quadratic <- tab_j[tab_j$Component_type == "Quadratic", ]$Component
  
  ### Calculation --------------------------------------------------------------
  ## Calculate M
  M0a <- mean(matrixStats::rowProds(.H[, i_split]))
  M0b <- mean(matrixStats::rowProds(.H[, j_split]))
  M1  <- mean(matrixStats::rowProds(.H[, ij]))
  M2  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% i_quadratic)]]))
  M3  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% j_quadratic)]]))
  M4  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% c(i_quadratic, j_quadratic))]]))
  M5  <- mean(matrixStats::rowProds(.H[, ij[!(ij %in% ij_match)]]))
  M6  <- mean(matrixStats::rowProds(.H[, c(i_split, j_single)]))
  M7  <- mean(matrixStats::rowProds(.H[, c(j_split, i_single)]))
  M8  <- mean(matrixStats::rowProds(.H[, c(i_single, j_single)]))
  M9  <- mean(matrixStats::rowProds(.H[, ij[ij == i_quadratic]]))
  
  ## Calculate denominator (D)
  D <- prod(.Q[ij])
  
  ## Calculate VCV element
  if(length(ij_match) == 0) {
    # 1.) Non of the components match
    
    x <- (M1 - (1 - .Q[i_quadratic]^2) * M2 - (1 - .Q[j_quadratic]^2) * M3 -
            (1 - .Q[i_quadratic]^2) * (1 - .Q[j_quadratic]^2) * M4) / D -
      (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && i_single == j_single) {
    # 2.) The "Single" component of the term .i matches the "Single"
    #     component of the term .j
    
    x <- (M1 - (1 - .Q[i_quadratic]^2) * M2 - (1 - .Q[j_quadratic]^2) * M3 -
            (1 - .Q[ij_match]^2) * M5) / D +
      (.Q[i_quadratic]^2 * .Q[j_quadratic]^2 * .Q[ij_match]^2 -
         .Q[i_quadratic]^2 - .Q[j_quadratic]^2 - .Q[ij_match]^2 + 2) / D -
      (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && i_single == j_quadratic) {
    # 3.) The "Single" component of the term .i matches the "Quadratic"
    #     component of the term .j
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M6 - (1 - .Q[i_quadratic]) * M7 -
            3 * (1 - .Q[i_single]^2) * (1 - .Q[i_quadratic]^2) * M8) / D -
      (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && i_quadratic == j_single) {
    # 4.) The "Single" component of the term .j matches the "Quadratic"
    #     component of the term .i
    
    x <- (M1 - 3 * (1 - .Q[j_single]^2) * M7 - (1 - .Q[j_quadratic]) * M6 -
            3 * (1 - .Q[j_single]^2) * (1 - .Q[j_quadratic]^2) * M8) / D -
      (M0a * M0b) / D
    
  } else if (length(ij_match) == 1 && i_quadratic == j_quadratic) {
    # 5.) The "Quadratic" component of the term .i matches the "Quadratic"
    #     component of the term .j
    
    x <- (M1 - 6 * (1 - .Q[i_quadratic]^2) * M6 +
            3 * (1 - .Q[i_quadratic]^2)^2 * M8) / D -
      (M0a * M0b) / D
    
  } else if (length(ij_match) == 2 &&
             i_single         == j_single &&
             i_quadratic      == j_quadratic) {
    # 6.) Both components match. The "Single" component of .i matches the "Single"
    #     component of .j and the "Quadratic" component of .i matches the
    #     "Quadratic" component of .j
    
    x <- (M1 - (1 - .Q[i_single]^2) * M9 - 6 * (1 - .Q[i_quadratic]^2) * M6 -
            3 * (1 - .Q[i_quadratic]^2) *
            (.Q[i_single]^2 * .Q[i_quadratic]^2 + .Q[i_single]^2 - 2)) / D -
      (M0a * M0b) / D
    
  } else {
    # 7.) Both components match. The "Single" component of .i matches the "Quadratic"
    #     component of .j and the "Quadratic" component of .i matches the
    #     "Single" component of .j
    
    x <- (M1 - 3 * (1 - .Q[i_single]^2) * M6 - 3 * (1 - .Q[i_quadratic]^2) * M7 +
            9 * (1 - .Q[i_quadratic]^2) * (1 - .Q[i_quadratic]^2) * M8) -
      (M0a * M0b) / D
    
  }
  ## Return
  x
}