#mode function
getmode <- function(v) {
 uniqv <- unique(v)
 uniqv[which.max(tabulate(match(v, uniqv)))]
}

incomeR_mode_income <- incomeR %>% 
 group_by(INCOME) %>%
 summarise(mode = getmode(OCCUPATION))
 
incomeR_mode_income