  We explicitly consider both phonon energy shifts and broadening arising from
both cubic and quartic anharmonicities, as well as diagonal/non-diagonal terms
of heat flux operators in thermal conductivity. Our findings show that the
strong anharmonicity of CsCu$_3$S$_2$ primarily arises from the presence of
$p$-$d$ anti-bonding hybridization between Cu and S atoms, coupled with the
random oscillations of Cs atoms. Notably, the competition between phonon
hardening described by the loop diagram and softening induced by the bubble
diagram significantly influences particle-like propagation, predominantly
reflected in group velocity and energy-conservation rule. Additionally, the
electrical transport properties are determined by employing the precise
momentum relaxation-time approximation (MRTA). At high temperatures, the
thermoelectric performance of $p$-type CsCu$_3$S$_2$ reaches its optimum
theoretical value of 0.94 along the in-plane direction based on advanced phonon
renormalization theory. In striking contrast, the harmonic approximation theory
significantly overestimates the thermoelectric efficiency at the same
temperatures, rendering it an impractical expectation. Conversely, the
first-order renormalization approach leads to a serious underestimation of the
thermoelectric properties due to the over-correction of phonon energy. Our
study not only reveals the pivotal role of anharmonic lattice dynamics in
accurately assessing thermoelectric properties but also underscores the
potential thermoelectric applications for novel copper-based chalcogenides.
