R/dce_efficiency.R

In JedStephens/ExpertChoice: Design of Discrete Choice and Conjoint Analysis

#' Efficiency Measures for Discrete Choice Experiments
#'
#' @param augmented_full_factorial The level augmented full factorial. See tutorial step 2.
#' @param choice_sets A list of choice sets generated by one of the methods used to convert from fractional factorial designs.
#'
#' @references Street, D.J., Burgess, L. and Louviere, J.J., 2005. Quick and easy choice sets: constructing optimal and nearly optimal stated choice experiments. International Journal of Research in Marketing, 22(4), pp.459-470.
#'
#' @return a list of named output.
#' @export
#'
#' @examples
#' # See Step 8 of the Practical Introduction to ExpertChoice vignette.
#' # Step 1
#' attrshort  = list(condition = c("0", "1", "2"),
#' technical =c("0", "1", "2"),
#' provenance = c("0", "1"))
#'
#' #Step 2
#' # ff stands for "full fatorial"
#'  ff  <-  full_factorial(attrshort)
#'  af  <-  augment_levels(ff)
#' # af stands for "augmented factorial"
#'
#' # Step 3
#' # Choose a design type: Federov or Orthogonal. Here an Orthogonal one is used.
#' nlevels <- unlist(purrr::map(ff, function(x){length(levels(x))}))
#' fractional_factorial <- DoE.base::oa.design(nlevels = nlevels, columns = "min34")
#'
#' # Step 4 & 5
#' # The functional draws out the rows from the original augmented full factorial design.
#' colnames(fractional_factorial) <- colnames(ff)
#' fractional <- search_design(ff, fractional_factorial)
#' # Step 5 (skipped, but important, see vignette)
#'
#' # Step 6
#' # Two modulators c(1,1,1) and c(0,1,1) are specified.
#' dce_modulo <- modulo_method(
#' fractional,
#' list(c(1,1,1),c(0,1,1))
#' )
#'
#' # Step 7 (skipped)
#'
#' # Step 8! -- Inspect the D-efficiency using the Street et. al method of the DCE design.
#' # NOTE: the af is used at this stage not the ff.
#' dce_efficiency(af, dce_modulo)
dce_efficiency <- function(augmented_full_factorial, choice_sets) {
  m <- unique(unlist(purrr::map(choice_sets, function(x){length(x)})))
  p <- ncol(attributes(augmented_full_factorial)$X_full) - 1
  Lamda <- .Lambda_matrix(augmented_full_factorial, choice_sets, m)
  C_mat <- .C_matrix(Lamda, augmented_full_factorial)
  # Determine the theoretical upper bound for main effects using the multinomial regression (See Street.)
  upper_bound_detC <- .detC_optimal(attributes(augmented_full_factorial)$out.attrs$dim, m)
  # Then determine the desings D-optimality.
  dce_d_effiency <- .dce_d_effiency(C_mat, upper_bound_detC, p)
  # Give some user output.
  message("\n\nThe D-efficiency of this discrete choice experiment is", round(dce_d_effiency, 3), "%", "\n")

  # Return the useful output.
  return(list(dce_d_effiency = dce_d_effiency, upper_bound_detC = upper_bound_detC, C_matrix = C_mat, Lamda = Lamda))
}


# Diagnositics about the DCE design.

.Lambda_matrix <- function(augmented_factorial, set_list, m) {
  # Street et al., 2005, pg. 462
  profiles <- augmented_factorial$levels
  lamda <- base::matrix(data = integer(), nrow = length(profiles), ncol = length(profiles))
  colnames(lamda) <- profiles
  rownames(lamda) <- profiles

  check_profiles <- function(x, p) {
    pair1 <- profiles %in% x
    pair2 <- pair1[p]
    pair_present <- pair1 * pair2
    return(pair_present)
  }
  for (r in 1:nrow(lamda)) {
    set_presence <- lapply(set_list, check_profiles, p = r)
    # return(set_presence)
    lamda_row_presence <- colSums(matrix(unlist(set_presence),
      byrow = TRUE,
      nrow = length(set_list), ncol = length(profiles)
    ))
    # return(lamda_row_presence)
    lamda[r, ] <- lamda_row_presence
  }
  diag(lamda) <- 0
  lamda <- -lamda
  diag(lamda) <- -rowSums(lamda)

  return(list(mat = lamda, prefactor = 1 / ((m^2) * (length(set_list)))))
}

.C_matrix <- function(lamda_mat, augmented_factorial) {
  # attributes(augmented_factorial)$B_mat %*% (lamda_mat$prefactor * lamda_mat$mat) %*% t(attributes(augmented_factorial)$B_mat)
  tcrossprod(attributes(augmented_factorial)$B_mat %*% (lamda_mat$prefactor * lamda_mat$mat), attributes(augmented_factorial)$B_mat)
}

# Helper function.
.is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol


.detC_optimal <- function(levels_vector, m) {
  # <- attributes(augmented_full_factorial)$out.attrs$dim
  k <- length(levels_vector)
  # A function.
  determine_s <- function(q, m, l) {
    # return
    if ((((l == 2) & (m %% 2 != 0)))) {
      # CASE 1
      message("  ", "Case 1\n")
      return((m^2 - 1) / 4)
    } else if (((l == 2) & (m %% 2 == 0))) {
      # CASE 2
      message("  ", "Case 2\n")
      return((m^2) / 4)
    } else if (((l > 2) & (l <= m))) {
      # CASE 3
      message("  ", "Case 3\n")
      # First determine y and then x.
      for (y in 0:(l - 1)) {
        x <- (m - y) / l # See page 463 of Street, Burgess, Louviere (Reference above in Function Descp).
        # If x is a whole number then we have a winner!
        if (.is.wholenumber(x)) {
          message("  ", "The implied x is", x, "and y is", y, "\n")
          return((m^2 - (l * x^2 + 2 * x * y + y)) / 2)
        }
      }
    } else if ((l >= m)) {
      # Case 4
      message("Case 4\n")
      return((m * (m - 1)) / 2)
    }
  }

  detC_theroy <- 1
  for (q in 1:k) {
    l <- levels_vector[q]
    message("\nq is", q, "\n")
    message("  ","L is", l, "\n")
    s <- determine_s(q, m, l)
    message("  ", "s is", s, "\n")

    notQ <- (1:k)[-q]
    bracket <- (2 * s) / (m^2 * (l - 1) * prod(levels_vector[notQ]))
    detC_theroy <- detC_theroy * ((bracket)^(l - 1))
  }
  names(detC_theroy) <- "upper_bound_detC"
  return(detC_theroy)
}


.dce_d_effiency <- function(C, C_optimal, p) {
  # message(paste("p is", p))
  dced <- ((det(C)) / C_optimal)^(1 / p) * 100
  names(dced) <- "dce_d_effiency"
  return(dced)
}