How to report results based on Likert item group differences? I collected Likert Item responses to a series of questions. My goal is to demonstrate that some of the groups are different from the others.
The dataset has the following format:
 UserID    Group      Q1      Q2      Q3 ...
 ---------------------------------
  user1    Group_A     5       4       5
  user2    Group_B     3       1       5 
  ...

Each of Q1, Q2 and Q3 is a Likert item, answered using the scale:
1: strongly disagree
2: disagree 
3: neither agree nor disagree
4: agree
5: Strongly agree 

The Group column contains 6 different groups (Group_A, Group_B, Group_C, Group_D, Group_E, Group_F).
In order to determine if there was a significant difference between any of the groups for each of the Likert Items, I used:
> kruskal.test(Q1 ~ Group, data=dt) 
Kruskal-Wallis rank sum test
data:  Q1 by Group
Kruskal-Wallis chi-squared = 148.14, df = 5, p-value < 2.2e-16

This highlights a significant difference between the groups, which may be further explored using:
> posthoc.kruskal.dunn.test(Q1 ~ group, data=dt)
Pairwise comparisons using Dunn's-test for multiple 
                     comparisons of independent samples 
data:  Q1 by Group 

                 Group_F       Group_A     Group_B       Group_C Group_D
Group_A          0.235           -           -             -                  -                   
Group_B        1.7e-14         < 2e-16     -             -                  -                   
Group_C        7.1e-05         1.2e-09     0.032         -                  -                   
Group_D        1.0e-05         8.6e-10     0.906         0.852              -                   
Group_E        9.5e-08         2.3e-13     0.906         0.726              0.990               
P value adjustment method: holm 

This helps to identify similarities/differences between groups.
Now it comes to reporting this data!
Based on this question, some suggestions are made around how to highlight a difference between two groups. However, how should the differences between multiple groups be published? Would it make sense to just state that a significant difference exists (and give the p-value from the Kruskal-Wallis test), and then somehow visually demonstrate the differences, such as a table with medians or a graph of some kind?
Any references you could recommend on this type of problem or other papers/reports that write up results similar to this?
 A: When presenting multiple comparisons among groups, it is helpful to use a compact letter display.  In this format, each group is assigned one or more letters.  Groups sharing a letter are not significantly different.
Creating a compact letter display can be done by visual inspection of your table of p-values, by an algorithm by hand, or with computer code. 
The package multcompLetters has a function to create the information for a compact letter display from a table of p-values.  However, this function requires a full square matrix of p-values, not the compact matrix that is outputted by several R functions.
The code below uses the rcompanion package to make creating the compact letter display easier.  This package relies on several other packages for things, so it may take a while for initial installation.
In all cases, when creating a compact letter display, be sure to double check the results with the original output.  It is easy along the way to make a mistake manipulating the matrix of p-values or assigning letters.  In general, it is best if you can arrange your matrix or output table by e.g. increasing median before creating a compact letter display to avoid spurious or weird results.
### Install packages

if(!require(PMCMR)){install.packages("PMCMR")}
if(!require(multcompView)){install.packages("multcompView")}
if(!require(rcompanion)){install.packages("rcompanion")}
if(!require(FSA)){install.packages("FSA")}
if(!require(ggplot2)){install.packages("ggplot2")}

### Create data

Group = c(rep("Group.A", 5),rep("Group.B", 5),rep("Group.C", 5),rep("Group.D", 5))
Value = c(5,6,7,7,8,6,7,8,9,9,15,16,17,18,19,10,11,13,14,15)
Data = data.frame(Group, Value) 

### Dunn test with PMCMR

library(PMCMR)

Post = posthoc.kruskal.dunn.test(Value ~ Group, data=Data)

PT = Post$p.value

PT

### Convert matrix to full matrix

library(rcompanion)

PT1 = fullPTable(PT)

PT1

### Create compact letter display information

library(multcompView)

multcompLetters(PT1,  
                compare="<",  
                threshold=0.05,
                Letters=letters,  
                reversed = FALSE)

In FSA there is a function for Dunn test whose output I like better, and I think is easier to convert to a compact letter display.
### Dunn test with FSA

library(FSA)

PT2 = dunnTest(Value ~ Group, data=Data,
              method="holm")

PT3= PT2$res

PT3

### Create compact letter display

library(rcompanion)

cldList(P.adj ~ Comparison, data=PT3,
        threshold  = 0.05)

In my opinion, it is appropriate to present Likert data with medians or other quantile-based measures. We usually do this by assigning numbers to the categories, but this is strictly for convenience.  R can determine quantiles for ordered category variables.
A barplot of counts is also appropriate for ordered category data.
Answers = c("A", "A", "A", "N", "N", "SA", "SA", "D", "SD")

Q1 = factor(Answers, ordered=TRUE,
            levels=c("SD", "D", "N", "A", "SA"))

min(Q1)

quantile(Q1, 0.25, type=1)

median(Q1)

quantile(Q1, 0.75, type=1)

max(Q1)

barplot(table(Q1))


However, if we are using numbers to represent the categories, results could be presented as box plots.
The code below also manually adds the letters from the compact letter display information above.
library(ggplot2)

  A = max(Data$Value[Data$Group=="Group.A"]) + 1
  B = max(Data$Value[Data$Group=="Group.B"]) + 1
  C = max(Data$Value[Data$Group=="Group.C"]) + 1
  D = max(Data$Value[Data$Group=="Group.D"]) + 1
  G = 0
  H = 22

  DF = data.frame(x = c(1,2,3,4), 
                  y = c(A,B,C,D), 
                  g = c("a", "a", "b","ab"))
  ggplot(Data, aes(x = Group, y = Value)) +
  geom_boxplot() +
    ylim(G,H) +

  geom_text(data = DF,
            aes(x=x,y=y,label=g)) +

  ylab("Value") +
  xlab("Group")


A: I simply would like to add this summary to Sal's answer. It provides a bit of additional info and alternative code, such as e.g.
library(emmeans)
library(multcomp)
library(multcompView)

model <- lm(weight ~ group, data = PlantGrowth)

# get weight means per group
model_means <- emmeans::emmeans(object = model,
                                pairwise ~ "group",
                                adjust = "tukey")

# add letters to each mean
model_means_cld <- multcomp::cld(object = model_means$emmeans,
                                 Letters = letters, 
                                 alpha = 0.05)

# show output
model_means_cld
#>  group emmean    SE df lower.CL upper.CL .group
#>  trt1    4.66 0.197 27     4.26     5.07  a    
#>  ctrl    5.03 0.197 27     4.63     5.44  ab   
#>  trt2    5.53 0.197 27     5.12     5.93   b   
#> 
#> Confidence level used: 0.95 
#> P value adjustment: tukey method for comparing a family of 3 estimates 
#> significance level used: alpha = 0.05 
#> NOTE: Compact letter displays can be misleading
#>       because they show NON-findings rather than findings.
#>       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

Created on 2021-08-13 by the reprex package (v2.0.0)

