R/calculate_degrees_of_freedom.R

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

#'@title Calculate Degrees of Freedom
#'
#'@description Calculates the degrees of freedom adjustment for models with random effects from Penheiro and Bates, pg. 91
#'@param run_matrix_processed The design
#'@param nointercept Whether the intercept term is present
#'@return Vector of numbers indicating split plot layers
#'@keywords internal
calculate_degrees_of_freedom = function(
  run_matrix_processed,
  nointercept,
  model,
  contrasts
) {
  blocklist = strsplit(rownames(run_matrix_processed), ".", fixed = TRUE)
  existingBlockStructure = do.call(rbind, blocklist)
  blockgroups = get_block_groups(existingBlockStructure)
  # Check block nesting structure
  block_nesting = c(calculate_block_nesting(
    blockgroups,
    existingBlockStructure
  ))
  blockgroups = blockgroups[-length(blockgroups)]
  #Check if all blocks are uniquely nested--if not, then it does not have a split-plot structure
  if (length(unique(block_nesting)) == length(block_nesting)) {
    split_plot_blocking = TRUE
  } else {
    split_plot_blocking = FALSE
    stop(
      "All blocks are not uniquely nested: design has a random block structure. ",
      "Parametric evaluation may not accurately assess the denominator degrees of freedom: ",
      "use `eval_design_mc(..., blocking = TRUE, adjust_alpha_inflation = TRUE)`",
      "for an accurate estimate of power."
    )
  }
  splitlayerlist = calculate_split_columns(
    run_matrix_processed,
    blockgroups,
    existingBlockStructure
  )

  if (!nointercept) {
    numberblocks = c(
      1,
      unlist(lapply(blockgroups, length)),
      nrow(run_matrix_processed)
    )
  } else {
    numberblocks = c(
      unlist(lapply(blockgroups, length)),
      nrow(run_matrix_processed)
    )
  }
  model = as.formula(paste0(
    "~",
    paste(attr(terms.formula(model), "term.labels"), collapse = " + ")
  ))
  splitlayer = splitlayerlist$splitlayer
  #`allsplit` determines if it's a split-plot or blocking structure by seeing if all the blocks are nested
  split_plot_structure = splitlayerlist$allsplit
  #And here is the level which each layer is nested under
  nested_level = c(0, splitlayerlist$nested_level + 1)
  splitterms = unlist(strsplit(
    as.character(model)[-1],
    split = " + ",
    fixed = TRUE
  ))
  splitterms = gsub("\\n", "", splitterms, perl = TRUE) #Remove newlines inserted when model character string too long
  splitterms = trimws(splitterms)
  interactions = is_intralayer_interaction(
    run_matrix_processed,
    model,
    splitlayer
  )

  if (!nointercept) {
    m = 1
    currentlayer = 1
    for (i in 2:(max(splitlayer) + 1)) {
      interaction_cols_layer = interactions[[currentlayer]]
      names(interaction_cols_layer) = splitterms
      layercols = which(currentlayer == splitlayer)
      #Total number of potential degrees of freedom for this layer is the number of blocks
      #minus the number of blocks from the previous layer.
      mtemp = numberblocks[i] - numberblocks[i - 1]

      #Subtract out degrees of freedom for terms nested within this split-plot layer
      if (split_plot_structure) {
        if (length(layercols) > 0) {
          for (j in seq_len(length(layercols))) {
            if (
              is.numeric(run_matrix_processed[, layercols[j], drop = FALSE])
            ) {
              mtemp = mtemp - 1
            } else {
              cat_degrees = length(unique(run_matrix_processed[,
                layercols[j],
                drop = FALSE
              ])) -
                1
              mtemp = mtemp - cat_degrees
            }
          }
        }
        #Now we have to deal with interaction terms. These can either be interactions between
        #terms fully contained within a layer (and thus will account for degrees of freedom
        #subtracted within that layer) or they will be an interaction with a term nested in a
        #higher layer (for which they will also count) or with a deeper layer (for which they
        #will not count--those degrees of freedom will be accounted for when we calculate the
        #df for that layer).
        for (j in seq_len(length(interaction_cols_layer))) {
          if (interaction_cols_layer[j]) {
            interaction_cols = unlist(strsplit(splitterms[j], ":"))
            higherorder_terms = unlist(grepl(
              "I\\((.+)\\^.+\\)",
              interaction_cols,
              perl = TRUE
            ))
            if (any(higherorder_terms)) {
              for (term in 1:length(higherorder_terms)) {
                if (higherorder_terms[term]) {
                  interaction_cols[term] = gsub(
                    "I\\(",
                    "",
                    interaction_cols[term],
                    perl = TRUE
                  )
                  interaction_cols[term] = gsub(
                    "\\^.+\\)",
                    "",
                    interaction_cols[term],
                    perl = TRUE
                  )
                }
              }
            }
            #Numeric factors get 1 df automatcally, otherwise count number of levels of
            #categorical factors minus 1 for the df. Multiply these numbers together
            #e.g. two 3 level factors interacting = 1*(3-1)*(3-1) = 4 df total.
            temp_degrees = 1
            for (col_name in interaction_cols) {
              if (!is.numeric(run_matrix_processed[, col_name, drop = FALSE])) {
                temp_degrees = temp_degrees *
                  (length(unique(run_matrix_processed[,
                    col_name,
                    drop = FALSE
                  ])) -
                    1)
              }
            }
            #Subtract from mtemp for this layer.
            mtemp = mtemp - temp_degrees
          }
        }
      }
      #Record current degrees of freedom and move to the next layer
      m = c(m, mtemp)
      currentlayer = currentlayer + 1
    }
  } else {
    #Subtract out degrees of freedom for terms nested within the first layer,
    #Which is slightly different in the no intercept case.

    currentlayer = 1
    layercols = which(currentlayer == splitlayer)
    mtemp = numberblocks[1]
    if (split_plot_structure) {
      if (length(layercols) > 0) {
        for (j in seq_len(length(layercols))) {
          if (is.numeric(run_matrix_processed[, layercols[j], drop = FALSE])) {
            mtemp = mtemp - 1
          } else {
            cat_degrees = length(unique(run_matrix_processed[,
              layercols[j],
              drop = FALSE
            ])) -
              1
            mtemp = mtemp - cat_degrees
          }
        }
      }
    }
    m = mtemp
    for (i in 2:(max(splitlayer))) {
      interaction_cols_layer = interactions[[currentlayer]]

      layercols = which(currentlayer == splitlayer)
      #Total number of potential degrees of freedom for this layer is the number of blocks
      #minus the number of blocks from the previous layer.

      mtemp = numberblocks[i] - numberblocks[i - 1]
      #Subtract out degrees of freedom for terms nested within split-plot layers
      if (split_plot_structure) {
        if (length(layercols) > 0) {
          for (j in seq_len(length(layercols))) {
            if (
              is.numeric(run_matrix_processed[, layercols[j], drop = FALSE])
            ) {
              mtemp = mtemp - 1
            } else {
              cat_degrees = length(unique(run_matrix_processed[,
                layercols[j],
                drop = FALSE
              ])) -
                1
              mtemp = mtemp - cat_degrees
            }
          }
        }
      }
      #Now we have to deal with interaction terms. These can either be interactions between
      #terms fully contained within a layer (and thus will account for degrees of freedom
      #subtracted within that layer) or they will be an interaction with a term nested in a
      #higher layer (for which they will also count) or with a deeper layer (for which they
      #will not count--those degrees of freedom will be accounted for when we calculate the
      #df for that layer).
      for (j in seq_len(length(interaction_cols_layer))) {
        if (interaction_cols_layer[j]) {
          interaction_cols = unlist(strsplit(splitterms[j], ":"))
          higherorder_terms = unlist(grepl(
            "I\\((.+)\\^.+\\)",
            interaction_cols,
            perl = TRUE
          ))
          if (any(higherorder_terms)) {
            for (term in 1:length(higherorder_terms)) {
              if (higherorder_terms[term]) {
                interaction_cols[term] = gsub(
                  "I\\(",
                  "",
                  interaction_cols[term],
                  perl = TRUE
                )
                interaction_cols[term] = gsub(
                  "\\^.+\\)",
                  "",
                  interaction_cols[term],
                  perl = TRUE
                )
              }
            }
          }
          #Numeric factors get 1 df automatcally, otherwise count number of levels of
          #categorical factors minus 1 for the df. Multiply these numbers together
          #e.g. two 3 level factors interacting = 1*(3-1)*(3-1) = 4 df total.
          temp_degrees = 1
          for (col_name in interaction_cols) {
            if (!is.numeric(run_matrix_processed[, col_name, drop = FALSE])) {
              temp_degrees = temp_degrees *
                (length(unique(run_matrix_processed[,
                  col_name,
                  drop = FALSE
                ])) -
                  1)
            }
          }
          #Subtract from mtemp for this layer.
          mtemp = mtemp - temp_degrees
        }
      }
      m = c(m, mtemp)
      currentlayer = currentlayer + 1
    }
  }
  degrees_of_freedom = calc_interaction_degrees(
    run_matrix_processed,
    model,
    contrasts,
    nointercept,
    splitlayer,
    m,
    split_plot_structure
  )
  if (!nointercept) {
    # From Pinheiro and Bates, pp. 91:
    # "The intercept, which is the parameter corresponding to the column of 1’s
    # in the model matrices Xi, is treated differently from all the other param-
    # eters, when it is present. As a parameter it is regarded as being estimated
    # at level 0 because it is outer to all the grouping factors. However, its de-
    # nominator degrees of freedom are calculated as if it were estimated at level
    # Q + 1. This is because the intercept is the one parameter that pools
    # information from all the observations at a level even when the corresponding
    # column in Xi doesn’t change with the level."
    degrees_of_freedom[1] = min(degrees_of_freedom[-1])
  }
  return(degrees_of_freedom)
}