It is not difficult to write a function in R to calculate the mode (as a central tendency). The function I have posted a while ago did not work with missing values. Here’s a tweak that does work with NA in the data.
Mode = function(x) {
ux = na.omit(unique(x))
return(ux[which.max(tabulate(match(x, ux)))])
}