R/revealAtte.R

In ramchoice: Revealed Preference and Attention Analysis in Random Limited Attention Models

################################################################################
#' @title Revealed Attention Analysis in Random Limited Attention Models
#'
#' @description Given a random sample of choice problems and choices, \code{revealAtte}
#'   returns the upper and lower bounds on the attention frequency following the construction
#'   of \href{https://arxiv.org/abs/2110.10650}{Cattaneo, Cheung, Ma, and Masatlioglu (2022)}.
#'
#' \code{\link{sumData}} is a low-level function that generates summary statistics. For
#'   revealed preference analysis, see \code{\link{revealPref}}.
#'
#' @param menu Numeric matrix of 0s and 1s, the collection of choice problems.
#' @param choice Numeric matrix of 0s and 1s, the collection of choices.
#' @param alternative Numeric vector, the alternatives for which to compute bounds on the
#'   attention frequency. For example, \code{c(1, 2, 4)} means the first, second, and fourth alternatives.
#' @param S Numeric matrix of 0s and 1s, the collection of choice problems to compute bounds on the
#'   attention frequency.
#' @param lower Boolean, whether lower bounds should be computed (default is \code{TRUE}).
#' @param upper Boolean, whether upper bounds should be computed (default is \code{TRUE}).
#' @param pref Numeric vector, corresponding to the preference. For example, \code{c(2, 3, 1)} means
#'   2 is preferred to 3 and to 1. When set to \code{NULL}, the default, \code{c(1, 2, 3, ...)},
#'   will be used. This option only applies to the upper bounds (i.e., when \code{upper} is set to \code{TRUE}).
#' @param nCritSimu Integer, number of simulations used to construct the critical value. Default is \code{2000}.
#' @param level Numeric, the significance level (default is \code{0.95}).
#'
#' @return
#' \item{sumStats}{Summary statistics, generated by \code{\link{sumData}}.}
#' \item{lowerBound}{Matrix containing the lower bounds.}
#' \item{upperBound}{Matrix containing the upper bounds.}
#' \item{critVal}{The simulated critical value.}
#' \item{opt}{Options used in the function call.}
#'
#' @references
#' M. D. Cattaneo, X. Ma, Y. Masatlioglu, and E. Suleymanov (2020). \href{https://arxiv.org/abs/1712.03448}{A Random Attention Model}. \emph{Journal of Political Economy} 128(7): 2796-2836. \doi{10.1086/706861}
#'
#' M. D. Cattaneo, P. Cheung, X. Ma, and Y. Masatlioglu (2022). \href{https://arxiv.org/abs/2110.10650}{Attention Overload}. Working paper.
#'
#' @author
#' Matias D. Cattaneo, Princeton University. \email{cattaneo@princeton.edu}.
#'
#' Paul Cheung, University of Maryland. \email{hycheung@umd.edu}
#'
#' Xinwei Ma (maintainer), University of California San Diego. \email{x1ma@ucsd.edu}
#'
#' Yusufcan Masatlioglu, University of Maryland. \email{yusufcan@umd.edu}
#'
#' Elchin Suleymanov, Purdue University. \email{esuleyma@purdue.edu}
#'
#' @examples
#' # Load data
#' data(ramdata)
#'
#' # Set seed, to replicate simulated critical values
#' set.seed(42)
#'
#' # preference
#' pref <- matrix(c(1, 2, 3, 4, 5), ncol=5, byrow=TRUE)
#' # list of choice problems
#' S <- matrix(c(1, 1, 0, 0, 0,
#'               1, 1, 1, 0, 0,
#'               1, 1, 1, 0, 1,
#'               1, 1, 1, 1, 1), ncol=5, byrow=TRUE)
#' result <- revealAtte(menu = ramdata$menu, choice = ramdata$choice,
#'   alternative = c(1,2), S = S,
#'   lower = TRUE, upper = TRUE,
#'   pref = pref)
#' summary(result)
#'
#' @export
revealAtte <- function(menu, choice,
                  alternative = NULL, S = NULL,
                  lower = TRUE, upper = TRUE,
                  pref = NULL,
                  nCritSimu = 2000,
                  level = 0.95) {

  ################################################################################
  # Error Check
  ################################################################################

  # menu
  if (!is.matrix(menu)) {
    stop("Input 'menu' has to be a matrix.\n")
  } else if (min(dim(menu)) == 0) {
    stop("Input 'menu' has at least one dimension 0.\n")
  } else if (!all((menu == 0) | (menu == 1))) {
    stop("Input 'menu' can only contain 0 (FALSE) and 1 (TRUE).\n")
  } else {
    # nothing
  }

  # choice
  if (!is.matrix(choice)) {
    stop("Input 'choice' has to be a matrix.\n")
  } else if (min(dim(choice)) == 0) {
    stop("Input 'choice' has at least one dimension 0.\n")
  } else if (!all((choice == 0) | (choice == 1))) {
    stop("Input 'choice' can only contain 0 (FALSE) and 1 (TRUE).\n")
  } else if (!all(rowSums(choice) == 1)) {
    stop("Input 'choice' should contain one and only one 1 (TRUE) in each row.\n")
  } else {
    # nothing
  }

  # menu and choice
  if (!all(dim(menu) == dim(choice))) {
    stop("Input 'menu' and 'choice' have to have the same dimensions.\n")
  }
  else if (!all(as.integer(rowSums((menu == 1) & (choice == 1))))) {
    stop("Input 'menu' and 'choice' should contain one and only one common 1 (TRUE) in each row.\n")
  } else {
    # nothing
  }

  # alternative
  if (length(alternative) == 0) {
    alternative <- 1
  } else if (!all(alternative %in% (1:ncol(menu)))) {
    stop("Option 'alternative' incorrectly specified.\n Should be integer(s) from 1, 2, ...\n")
  } else {
    # nothing
  }

  # S
  if (!is.matrix(S)) {
    stop("Input 'S' has to be a matrix.\n")
  } else if (min(dim(S)) == 0) {
    stop("Input 'S' has at least one dimension 0.\n")
  } else if (!all((S == 0) | (S == 1))) {
    stop("Input 'S' can only contain 0 (FALSE) and 1 (TRUE).\n")
  } else {
    # nothing
  }

  # alternative and S
  if (!all(apply(S, MARGIN=1, FUN=function(x) all(x[alternative] == 1)))) {
    stop("Some 'alternative' do not belong to the choice problems specified in 'S'.\n")
  } else {
    # nothing
  }

  # option lower
  if (length(lower) == 0) {
    lower <- TRUE
  } else if (length(lower) > 1 | !lower[1]%in%c(TRUE, FALSE)) {
    stop("Option 'lower' incorrectly specified.\n")
  }

  # option upper
  if (length(upper) == 0) {
    upper <- TRUE
  } else if (length(upper) > 1 | !upper[1]%in%c(TRUE, FALSE)) {
    stop("Option 'upper' incorrectly specified.\n")
  }

  # check if both lower and upper are FALSE
  if (!lower & !upper) {
    stop("At least one option, lower or upper, has to be used.\n")
  }

  # preference
  if (length(as.vector(pref)) == 0) {
    pref <- matrix(1:ncol(menu), nrow=1)
  } else if (!is.matrix(pref)) {
    stop("Input 'pref' has to be a row vector/matrix.\n")
  } else if (min(dim(pref)) == 0) {
    stop("Input 'pref' has at least one dimension 0.\n")
  } else if (ncol(pref) != ncol(menu)) {
    stop("Input 'pref' has to have the same number of columns as 'menu'.\n")
  } else if (!all(apply(pref, MARGIN=1, FUN=function(x) all(x %in% (1:ncol(menu)))))) {
    stop("Input 'pref' incorrectly specified.\n")
  } else {
    # nothing
  }

  ################################################################################
  # summary statistics
  ################################################################################
  sumStats <- sumData(menu, choice)
  n <- sum(sumStats$sumN) # sample size

  ################################################################################
  # critical value generation
  ################################################################################
  Rlower <- c()
  if (lower) {
    # find all supersets
    SupersetsIndex <- c()
    for (sIndex in 1:nrow(S)) {
      tempSupersetsIndex <- (1:nrow(sumStats$sumMenu))[apply(sumStats$sumMenu, MARGIN=1, FUN=function(x) all(x - S[sIndex, ] >= 0))]
      if (length(tempSupersetsIndex) == 0) {
        stop(paste("The ", sIndex, "th choice problem specified in 'S' does not have a superset in the data.\n Lower bounds cannot be computed.\n", sep=""))
      }
      SupersetsIndex <- c(SupersetsIndex, tempSupersetsIndex)
    }
    SupersetsIndex <- sort(unique(SupersetsIndex))

    # construct the constraint matrix
    for (i in SupersetsIndex) for (j in alternative) {
      temp <- rep(0, length(sumStats$sumProbVec))
      temp[sum(sumStats$sumMsize[1:i]) - sumStats$sumMsize[i] + sum(sumStats$sumMenu[i, 1:j])] <- 1
      Rlower <- rbind(Rlower, temp)
    }
    nRlower <- nrow(Rlower)
  }

  Rupper <- c()
  if (upper) {
    for (j in alternative) {
      SubsetsIndex <- c()
      for (sIndex in 1:nrow(S)) {
        tempSubsetsIndex <- (1:nrow(sumStats$sumMenu))[apply(sumStats$sumMenu, MARGIN=1, FUN=function(x) all(x - S[sIndex, ] <= 0) & (x[j] == 1))]
        if (length(tempSubsetsIndex) == 0) {
          stop(paste("The ", sIndex, "th choice problem specified in 'S' does not have a subset in the data.\n Upper bounds cannot be computed.\n", sep=""))
        }
        SubsetsIndex <- c(SubsetsIndex, tempSubsetsIndex)
      }
      SubsetsIndex <- sort(unique(SubsetsIndex))

      for (i in SubsetsIndex) {
        upperContourIndex <- (pref[1:which.max(pref == j)])[sumStats$sumMenu[i, pref[1:which.max(pref == j)]] == 1]
        # determine if the alternative is least preferred
        if (length(upperContourIndex) < sumStats$sumMsize[i]) {
          temp <- rep(0, length(sumStats$sumProbVec))
          for (k in upperContourIndex) {
            temp[sum(sumStats$sumMsize[1:i]) - sumStats$sumMsize[i] + sum(sumStats$sumMenu[i, 1:k])] <- 1
          }
          Rupper <- rbind(Rupper, temp)
        }
      }
    }
    nRupper <- nrow(Rupper)
  }
  R <- rbind(Rlower, Rupper)

  RSigma <- sqrt(abs(diag(R %*% sumStats$Sigma %*% t(R))))
  normSimu <- t(mvrnorm(n=nCritSimu, mu=rep(0, nrow(sumStats$Sigma)), Sigma=sumStats$Sigma))
  normSimu[diag(sumStats$Sigma) == 0, ] <- 0
  normSimu <- R %*% normSimu
  if (lower & (!upper)) {
    # nothing
  } else if (lower & upper) {
    normSimu[(nRlower+1):(nRlower+nRupper), ] <- -1 * normSimu[(nRlower+1):(nRlower+nRupper), ]
  } else {
    normSimu <- -1 * normSimu
  }
  normSimu <- apply(normSimu, MARGIN=2, FUN=function(x) max(x / RSigma))
  normSimu[is.na(normSimu)] <- 0
  critVal <- quantile(normSimu, level)

  ################################################################################
  # compute the bounds
  ################################################################################
  lowerBound <- upperBound <- matrix(NA, ncol=nrow(S), nrow=length(alternative))

  # lower bounds
  if (lower) {
    for (j in 1:length(alternative)) for (i in 1:nrow(S)) {
      SupersetsIndex <- c()
      tempSupersetsIndex <- (1:nrow(sumStats$sumMenu))[apply(sumStats$sumMenu, MARGIN=1, FUN=function(x) all(x - S[i, ] >= 0))]
      for (k in tempSupersetsIndex) {
        temp <- rep(0, length(sumStats$sumProbVec))
        temp[sum(sumStats$sumMsize[1:k]) - sumStats$sumMsize[k] + sum(sumStats$sumMenu[k, 1:alternative[j]])] <- 1
        SupersetsIndex <- rbind(SupersetsIndex, temp)
      }
      lowerBound[j, i] <- max(SupersetsIndex %*% sumStats$sumProbVec - critVal * sqrt(abs(diag(SupersetsIndex %*% sumStats$Sigma %*% t(SupersetsIndex))) / n))
    }
  }

  # upper bounds
  if (upper) {
    for (j in 1:length(alternative)) for (i in 1:nrow(S)) {
      SubsetsIndex <- c()
      tempSubsetsIndex <- (1:nrow(sumStats$sumMenu))[apply(sumStats$sumMenu, MARGIN=1, FUN=function(x) all(x - S[i, ] <= 0) & (x[alternative[j]] == 1))]
      for (k in tempSubsetsIndex) {
        upperContourIndex <- (pref[1:which.max(pref == alternative[j])])[sumStats$sumMenu[k, pref[1:which.max(pref == alternative[j])]] == 1]
        temp <- rep(0, length(sumStats$sumProbVec))
        for (ii in upperContourIndex) {
          temp[sum(sumStats$sumMsize[1:k]) - sumStats$sumMsize[k] + sum(sumStats$sumMenu[k, 1:ii])] <- 1
        }
        SubsetsIndex <- rbind(SubsetsIndex, temp)
      }
      upperBound[j, i] <- min(SubsetsIndex %*% sumStats$sumProbVec + critVal * sqrt(abs(diag(SubsetsIndex %*% sumStats$Sigma %*% t(SubsetsIndex))) / n))
    }
  }

    Result <- list(sumStats=sumStats,
                   lowerBound = lowerBound,
                   upperBound = upperBound,
                   critVal = critVal,
                   opt=list(alternative = alternative,
                            S = S,
                            lower = lower,
                            upper = upper,
                            pref = pref,
                            nCritSimu = nCritSimu,
                            level = level)
                   )

    class(Result) <- "ramchoiceRevealAtte"
    return(Result)
}

################################################################################
#' Internal function.
#'
#' @param object Class \code{ramchoiceRevealAtte} objects.
#'
#' @keywords internal
#' @export
summary.ramchoiceRevealAtte <- function(object, ...) {
  x <- object
  cat("\n Revealed Attention Analysis in Random Limited Attention Models.\n")
  cat("\n")

  cat(paste(format("# of observations", width=25), toString(sum(x$sumStats$sumN)), sep="")); cat("\n")
  cat(paste(format("# of alternatives", width=25), toString(ncol(x$sumStats$sumMenu)), sep="")); cat("\n")
  cat(paste(format("# of choice problems", width=25), toString(nrow(x$sumStats$sumMenu)), sep="")); cat("\n")
  cat("\n")

  cat(paste(format("Min # of alternatives", width=25), toString(min(x$sumStats$sumMsize)), sep="")); cat("\n")
  cat(paste(format("Max # of alternatives", width=25), toString(max(x$sumStats$sumMsize)), sep="")); cat("\n")
  cat("\n")

  cat(paste(format("Min eff. observations", width=25), toString(min(x$sumStats$sumN)), sep="")); cat("\n")
  cat(paste(format("Max eff. observations", width=25), toString(max(x$sumStats$sumN)), sep="")); cat("\n")
  cat("\n")

  if (x$opt$upper) {
    cat(paste(format("Preference", width=12), gsub(",", " >", toString(x$opt$pref)), sep="")); cat("\n")
    cat("\n")
  }

  cat(paste(format("Critical value", width=25), toString(round(x$critVal, 3)), sep="")); cat("\n")
  cat("\n")

  for (i in 1:length(x$opt$alternative)) {
    cat(paste(rep("=", 60), collapse="")); cat("\n")
    cat(paste("Alternative[", toString(i), "]",": ", sep=""))
    cat(x$opt$alternative[i]); cat("\n")

    cat(paste(format("lower", width=10), sep=""))
    cat(paste(format("upper", width=10), sep=""))
    cat(paste(format("choice problem", width=10), sep="")); cat("\n")

    for (j in 1:nrow(x$opt$S)) {
      if (x$opt$lower) {
        cat(paste(format(toString(round(x$lowerBound[i, j], 4)), width=10), sep=""))
      } else {
        cat(paste(format(" ", width=10), sep=""))
      }

      if (x$opt$upper) {
        cat(paste(format(toString(round(x$upperBound[i, j], 4)), width=10), sep=""))
      } else {
        cat(paste(format(" ", width=10), sep=""))
      }
      cat(paste("{ ", toString((1:ncol(x$opt$S))[x$opt$S[j, ] == 1]), " }", sep="")); cat("\n")
    }
  }
  cat(paste(rep("=", 60), collapse="")); cat("\n")
  cat("\n")
}

################################################################################
#' Internal function.
#'
#' @param object Class \code{ramchoiceRevealAtte} objects.
#'
#' @keywords internal
#' @export
print.ramchoiceRevealAtte <- function(x, ...) {
  cat("\n Revealed Attention Analysis in Random Limited Attention Models.\n")
  cat("\n")

  cat(paste(format("# of observations", width=25), toString(sum(x$sumStats$sumN)), sep="")); cat("\n")
  cat(paste(format("# of alternatives", width=25), toString(ncol(x$sumStats$sumMenu)), sep="")); cat("\n")
  cat(paste(format("# of choice problems", width=25), toString(nrow(x$sumStats$sumMenu)), sep="")); cat("\n")
  cat("\n")

  cat(paste(format("Min # of alternatives", width=25), toString(min(x$sumStats$sumMsize)), sep="")); cat("\n")
  cat(paste(format("Max # of alternatives", width=25), toString(max(x$sumStats$sumMsize)), sep="")); cat("\n")
  cat("\n")

  cat(paste(format("Min eff. observations", width=25), toString(min(x$sumStats$sumN)), sep="")); cat("\n")
  cat(paste(format("Max eff. observations", width=25), toString(max(x$sumStats$sumN)), sep="")); cat("\n")
  cat("\n")

  if (x$opt$upper) {
    cat(paste(format("Preference", width=12), gsub(",", " >", toString(x$opt$pref)), sep="")); cat("\n")
    cat("\n")
  }

  cat(paste(format("Critical value", width=25), toString(round(x$critVal, 3)), sep="")); cat("\n")
  cat("\n")

  for (i in 1:length(x$opt$alternative)) {
    cat(paste(rep("=", 60), collapse="")); cat("\n")
    cat(paste("Alternative[", toString(i), "]",": ", sep=""))
    cat(x$opt$alternative[i]); cat("\n")

    cat(paste(format("lower", width=10), sep=""))
    cat(paste(format("upper", width=10), sep=""))
    cat(paste(format("choice problem", width=10), sep="")); cat("\n")

    for (j in 1:nrow(x$opt$S)) {
      if (x$opt$lower) {
        cat(paste(format(toString(round(x$lowerBound[i, j], 4)), width=10), sep=""))
      } else {
        cat(paste(format(" ", width=10), sep=""))
      }

      if (x$opt$upper) {
        cat(paste(format(toString(round(x$upperBound[i, j], 4)), width=10), sep=""))
      } else {
        cat(paste(format(" ", width=10), sep=""))
      }
      cat(paste("{ ", toString((1:ncol(x$opt$S))[x$opt$S[j, ] == 1]), " }", sep="")); cat("\n")
    }
  }
  cat(paste(rep("=", 60), collapse="")); cat("\n")
  cat("\n")
}