Here is a function for finding all of the modes (i.e. peaks), in a set of data.   The return object is a list of the found peaks by their position (x value), along with their height relative to the highest peak.   Plus a lot of other information.   To be clear, in a simple case you can just do something like:

order(table(x), decreasing=TRUE)[1]

But I often come across data that is composed of large sets of (potentially) multi-modal, real value numbers, so that method won’t work.  Hence the following code:

#
# Examples:
# x = round(runif(1000, min=1, max=23))
# findMode(x)
# findMode(x, smallest.expected=0.5)
# findMode(x, adjust=0.1)
#
findMode <- function (data, smallest.expected=0, ...) {

  d              = density(data, ...)
  maxLocations   = vector()
  minLocations   = vector()
  maxHeights     = vector()
  minDepths      = vector()
  maxValue       = max(d$y)
  exactMax       = order(table(data), decreasing=TRUE)[1]

  # This will not capture the extreme edges as either a min or a max
  wasIncreasing <- TRUE;
  for (j in 2:length(d$y)) {
    
    # Hit a peak
    if (wasIncreasing & (d$y[j] < d$y[j-1])) {
	  maxLocations  = c(maxLocations, d$x[j])
	  maxHeights    = c(maxHeights, d$y[j]/maxValue)
	  wasIncreasing = FALSE

	# Hit a valley
	} else if (! wasIncreasing & (d$y[j] > d$y[j-1])) {
	  minLocations   = c(minLocations, d$x[j]);
	  minDepths      = c(minDepths, d$y[j]/maxValue);
	  wasIncreasing  = TRUE;

	# In a flat spot
	} else if (d$y[j] == d$y[j-1]) {
	  next;
	}
  }

  # In addition to changing the density estimation model (see 'density'),
  # one may want to sweep away tiny peaks that are just noise.  
  if (smallest.expected > 0) {
	if (length(minDepths) > 0) {
	  if (max(minDepths) > 0) {
#            keeperIndex    = ((minDepths / max(minDepths)) > smallest.expected);
		keeperIndex    = (minDepths  > smallest.expected);
		minDepths      = minDepths[keeperIndex];
		minLocations   = minLocations[keeperIndex];
	  }
	}

#        keeperIndex    = ((maxHeights / maxValue) > smallest.expected);
	keeperIndex    = (maxHeights > smallest.expected);
	maxHeights     = maxHeights[keeperIndex];
	maxLocations   = maxLocations[keeperIndex];
  }

  # Collect these features for this channel of the flowFrame
  modeList            = list();
  modeList$exactMax   = exactMax
  modeList$maximum    = d$x[d$y == max(d$y)][1];
  modeList$maxima     = maxLocations;
  modeList$maxHeights = maxHeights;
  modeList$minima     = minLocations;
  modeList$minDepths  = minDepths;
  return(modeList);
}