R/bdf3.mef.R

In bdf3: Efficient Block Designs for 3-Level Factorial Experiments in Block Size 3

#' @importFrom dplyr %>%
#' @importFrom stats na.omit
#'
bdf3.mef <- function(n_factors = 3, show_efficiency = TRUE) {
  get_principal_blocks <- function(n) {
    factors <- LETTERS[1:n]
    all_runs <- expand.grid(replicate(n, 0:2, simplify = FALSE))
    all_runs <- as.matrix(all_runs)
    zero_run <- rep(0, n)
    blocks <- list()
    for (i in 1:(nrow(all_runs) - 1)) {
      for (j in (i + 1):nrow(all_runs)) {
        r1 <- all_runs[i, ]
        r2 <- all_runs[j, ]
        block <- rbind(zero_run, r1, r2)
        if (all(apply(block, 2, function(col) length(unique(col)) == 3))) {
          blocks[[length(blocks) + 1]] <- as.data.frame(block)
          colnames(blocks[[length(blocks)]]) <- factors
        }
      }
    }
    return(blocks)
  }
  label_effect <- function(vec, factors) {
    label <- ""
    superscripts <- c("0" = "", "1" = "", "2" = "\u00B2")
    for (i in seq_along(vec)) {
      if (vec[i] != 0) {
        label <- paste0(label, factors[i], superscripts[as.character(vec[i])])
      }
    }
    if (label == "") "I" else label
  }
  generate_effects <- function(n) {
    total <- 3^n
    effects <- matrix(NA, nrow = total - 1, ncol = n)
    row <- 1
    for (i in 1:(total - 1)) {
      vec <- integer(n)
      num <- i
      for (j in 1:n) {
        vec[j] <- num %% 3
        num <- num %/% 3
      }
      effects[row, ] <- vec
      row <- row + 1
    }
    effects
  }
  is_confounded <- function(block, effect) {
    res <- (as.matrix(block) %*% effect) %% 3
    length(unique(as.vector(res))) == 1
  }
  normalize_effect <- function(eff) {
    first_nonzero <- which(eff != 0)[1]
    if (length(first_nonzero) == 0) return(eff)
    if (eff[first_nonzero] == 2) {
      eff <- (3 - eff) %% 3
    }
    eff
  }
  classify_effect <- function(eff) {
    nonzeros <- sum(eff != 0)
    if (nonzeros == 1) return("Main")
    else if (nonzeros == 2) return("2FI")
    else return("Other")
  }
  compute_efficiency_factors <- function(blocks, normalized_effects, effect_labels) {
    num_blocks <- length(blocks)
    estimable_count <- numeric(length(effect_labels))
    for (blk in blocks) {
      confounded <- apply(normalized_effects, 1, function(eff) is_confounded(blk, eff))
      estimable_count <- estimable_count + as.numeric(!confounded)
    }
    eff_factor <- round(estimable_count / num_blocks, 2)
    names(eff_factor) <- effect_labels
    return(eff_factor)
  }
  generate_unique_blocks <- function(principal_block) {
    n_factors <- ncol(principal_block)
    pb <- as.matrix(principal_block)
    used_treatments_matrix <- matrix(NA, nrow = 0, ncol = n_factors)
    blocks <- list()
    blocks[[1]] <- pb
    used_treatments_matrix <- rbind(used_treatments_matrix, pb)
    shifts <- expand.grid(replicate(n_factors, 0:2, simplify = FALSE))
    shifts <- shifts[!apply(shifts, 1, function(x) all(x == 0)), ]
    for (i in 1:nrow(shifts)) {
      shift <- as.numeric(shifts[i, ])
      derived_block <- (pb + matrix(rep(shift, each = nrow(pb)), nrow = nrow(pb), byrow = FALSE)) %% 3
      already_used <- apply(derived_block, 1, function(row)
        any(apply(used_treatments_matrix, 1, function(used_row) all(used_row == row))))
      if (!any(already_used)) {
        blocks[[length(blocks) + 1]] <- derived_block
        used_treatments_matrix <- rbind(used_treatments_matrix, derived_block)
      }
      if (nrow(used_treatments_matrix) == 3^n_factors) break
    }
    return(blocks)
  }
  factors <- LETTERS[1:n_factors]
  all_blocks <- get_principal_blocks(n_factors)
  all_effects <- generate_effects(n_factors)
  classified <- apply(all_effects, 1, classify_effect)
  target_effects <- all_effects[classified %in% c("Main", "2FI"), , drop = FALSE]
  normalized_effects <- unique(t(apply(target_effects, 1, normalize_effect)))
  valid_rows <- apply(normalized_effects, 1, function(eff) {
    first <- which(eff != 0)[1]
    !is.na(first) && eff[first] == 1
  })
  normalized_effects <- normalized_effects[valid_rows, , drop = FALSE]
  effect_labels <- apply(normalized_effects, 1, function(eff) label_effect(eff, factors))

  block_estimable <- lapply(all_blocks, function(blk) {
    apply(normalized_effects, 1, function(eff) !is_confounded(blk, eff))
  })

  covered <- c()
  chosen_blocks <- list()
  used_indices <- c()
  while (length(covered) < nrow(normalized_effects)) {
    best_block <- NULL
    max_new <- 0
    for (i in seq_along(block_estimable)) {
      if (i %in% used_indices) next
      new_cov <- which(block_estimable[[i]])[!(which(block_estimable[[i]]) %in% covered)]
      if (length(new_cov) > max_new) {
        max_new <- length(new_cov)
        best_block <- i
      }
    }
    if (is.null(best_block)) break
    covered <- union(covered, which(block_estimable[[best_block]]))
    chosen_blocks[[length(chosen_blocks) + 1]] <- as.matrix(all_blocks[[best_block]])
    used_indices <- c(used_indices, best_block)
  }

  eff_factors <- NULL
  if (show_efficiency) {
    eff_factors <- compute_efficiency_factors(chosen_blocks, normalized_effects, effect_labels)
  }

  structure(
    list(
      factors = factors,
      normalized_effects = normalized_effects,
      effect_labels = effect_labels,
      chosen_principal_blocks = chosen_blocks,
      efficiency_factors = eff_factors,
      replications = lapply(chosen_blocks, generate_unique_blocks)
    ),
    class = "bdf3mef"
  )
}
print.bdf3mef <- function(x, ...) {
  is_confounded <- function(block, effect) {
    res <- (as.matrix(block) %*% effect) %% 3
    length(unique(as.vector(res))) == 1
  }
  label_effect <- function(vec, factors) {
    label <- ""
    superscripts <- c("0" = "", "1" = "", "2" = "\u00B2")
    for (i in seq_along(vec)) {
      if (vec[i] != 0) {
        label <- paste0(label, factors[i], superscripts[as.character(vec[i])])
      }
    }
    if (label == "") "I" else label
  }
  classify_effect <- function(eff) {
    nonzeros <- sum(eff != 0)
    if (nonzeros == 1) return("Main")
    else if (nonzeros == 2) return("2FI")
    else return("Other")
  }

  cat("Minimal Principal Blocks Covering All Main Effects and 2FIs\n")
  cat("Factors:", paste(x$factors, collapse = " "), "\n\n")

  for (i in seq_along(x$chosen_principal_blocks)) {
    block <- x$chosen_principal_blocks[[i]]
    colnames(block) <- x$factors
    cat(paste0("Principal Block ", i, ":\n"))
    print(as.data.frame(block), row.names = FALSE)

    confounded <- apply(x$normalized_effects, 1, function(eff) {
      eff <- as.numeric(eff)
      if (is_confounded(block, eff)) label_effect(eff, x$factors) else NA
    })

    classified <- apply(x$normalized_effects, 1, classify_effect)
    confounded_main <- na.omit(confounded[classified == "Main"])
    confounded_2FI  <- na.omit(confounded[classified == "2FI"])

    cat("Confounded Main Effects: ",
        if (length(confounded_main) == 0) "None" else paste(confounded_main, collapse = ", "),
        "\n")
    cat("Confounded 2FI Effects: ",
        if (length(confounded_2FI) == 0) "None" else paste(confounded_2FI, collapse = ", "),
        "\n\n")
  }

  if (!is.null(x$efficiency_factors)) {
    eff_factors <- x$efficiency_factors
    effect_names <- names(eff_factors)
    eff_values <- as.vector(eff_factors)
    col_width <- max(nchar(effect_names), nchar(sprintf("%.2f", eff_values))) + 4
    cols_per_line <- 8
    center_align <- function(text, width) {
      pad_total <- width - nchar(text)
      pad_left <- floor(pad_total / 2)
      pad_right <- ceiling(pad_total / 2)
      paste0(strrep(" ", pad_left), text, strrep(" ", pad_right))
    }
    centered_effects <- sapply(effect_names, center_align, width = col_width)
    centered_values  <- sapply(sprintf("%.2f", eff_values), center_align, width = col_width)
    n <- length(eff_values)
    chunks <- ceiling(n / cols_per_line)

    cat("\nEfficiency Factor:\n\n")
    for (i in seq_len(chunks)) {
      idx <- ((i - 1) * cols_per_line + 1):min(i * cols_per_line, n)
      cat("Effects:           ", paste(centered_effects[idx], collapse = ""), "\n")
      cat("Efficiency Factor: ", paste(centered_values[idx], collapse = ""), "\n\n")
    }
  }

  for (pb_index in seq_along(x$replications)) {
    cat("\n\nReplication", pb_index, "\n")
    full_blocks <- x$replications[[pb_index]]
    n_blocks <- length(full_blocks)
    n_factors <- ncol(x$chosen_principal_blocks[[pb_index]])
    factor_names <- x$factors
    output_matrix <- matrix("", nrow = nrow(x$chosen_principal_blocks[[pb_index]]) + 2,
                            ncol = n_blocks * (n_factors + 1))
    for (i in seq_along(full_blocks)) {
      block <- full_blocks[[i]]
      start_col <- (i - 1) * (n_factors + 1) + 1
      end_col <- start_col + n_factors
      output_matrix[1, start_col] <- sprintf("Block %d", i)
      output_matrix[2, (start_col + 1):end_col] <- factor_names
      output_matrix[3:(nrow(x$chosen_principal_blocks[[pb_index]]) + 2),
                    (start_col + 1):end_col] <- apply(block, 2, as.character)
    }
    output_df <- as.data.frame(output_matrix, stringsAsFactors = FALSE)
    colnames(output_df) <- rep("", ncol(output_df))
    print(output_df, row.names = FALSE, quote = FALSE, right = FALSE)
  }

  invisible(x)
}