R/gen_momentsmatrix_ideal.R

In tylermorganwall/skpr: Design of Experiments Suite: Generate and Evaluate Optimal Designs

#'@title Generates the moment matrix, assuming no constraints
#'
#'@description Returns number of levels prior to each parameter
#'
#'@keywords internal
#'@return Returns a vector consisting of the number
#'of levels preceding each parameter (including the intercept)
gen_momentsmatrix_ideal = function(modelfactors, classvector) {

  #parse factor types
  isintercept = modelfactors == "(Intercept)"
  ishigherorder = lapply(
    strsplit(modelfactors, split = "^", fixed = TRUE),
    length
  ) >
    1
  isinteraction = lapply(
    strsplit(modelfactors, split = ":", fixed = TRUE),
    length
  ) >
    1
  isnumeric = modelfactors %in% names(classvector)[!classvector]
  islinear = rep(TRUE, length(modelfactors)) &
    !isintercept &
    !ishigherorder &
    !isinteraction
  isfullfactor = !isintercept & !ishigherorder & !isinteraction & !isnumeric

  linearfactors = modelfactors[isfullfactor]
  linearterms = modelfactors[islinear]

  if (!all(islinear) && all(isinteraction[-1])) {
    stop(
      "skpr: No main effects in model. skpr only supports interactions between terms when their main effects are present."
    )
  }

  ordermatrix = matrix(
    0,
    nrow = length(linearterms),
    ncol = length(modelfactors)
  )
  momentsmatrix = matrix(
    list(2),
    nrow = length(modelfactors),
    ncol = length(modelfactors)
  )
  momentsmatrixresults = matrix(
    0,
    nrow = length(modelfactors),
    ncol = length(modelfactors)
  )

  for (i in 1:length(which(islinear))) {
    ordermatrix[i, which(islinear)[i]] = 1
  }

  for (i in 1:length(linearterms)) {
    if (any(ishigherorder)) {
      for (j in 2:ncol(ordermatrix)) {
        if (
          length(strsplit(
            modelfactors[j],
            split = paste(linearterms[i], "^", sep = ""),
            fixed = TRUE
          )[[1]]) >
            1
        ) {
          ordermatrix[i, j] = as.numeric(gsub(
            "\\D",
            "",
            strsplit(modelfactors[j], split = "^", fixed = TRUE)[[1]][2]
          ))
        }
      }
    }
  }
  for (i in 1:length(linearterms)) {
    for (j in 2:ncol(ordermatrix)) {
      if (
        isinteraction[j] &&
          (linearterms[i] %in%
            strsplit(modelfactors[j], split = ":", fixed = TRUE)[[1]])
      ) {
        ordermatrix[i, j] = ordermatrix[i, j] + 1
      }
    }
  }
  ordermatrix[, 1] = 0

  for (i in 1:length(modelfactors)) {
    for (j in 1:length(modelfactors)) {
      momentsmatrix[i, j] = list(ordermatrix[, i] + ordermatrix[, j])

      if (any((ordermatrix[, i] + ordermatrix[, j]) %% 2 == 1)) {
        momentsmatrix[i, j] = 0
      } else {
        momentsmatrixresults[i, j] = prod(1 / (momentsmatrix[i, j][[1]] + 1))
      }
    }
  }

  for (i in 1:length(modelfactors)) {
    for (j in 1:length(modelfactors)) {
      if (
        (isfullfactor[j] || isfullfactor[i]) &&
          !(isinteraction[i] || isinteraction[j])
      ) {
        momentsmatrixresults[i, j] = 3 * momentsmatrixresults[i, j]
      } else {
        if (
          (isfullfactor[j] || isfullfactor[i]) ||
            (isinteraction[i] || isinteraction[j])
        ) {
          term = strsplit(modelfactors[i], split = ":", fixed = TRUE)[[1]]
          catterms = sum(term %in% linearfactors)
          momentsmatrixresults[i, j] = 3^catterms * momentsmatrixresults[i, j]
        }
      }
    }
  }
  return(momentsmatrixresults)
}