For a given diagrammatic approximation in many-body perturbation theory it is not guaranteed that positive observables, such as the density or the spectral function, retain their positivity. For zero-temperature systems we developed a method [Phys.Rev.B{\bf 90},115134 (2014)] based on so-called cutting rules for Feynman diagrams that enforces these properties diagrammatically, thus solving the problem of negative spectral densities observed for various vertex approximations. In this work we extend this method to systems at finite temperature by formulating the cutting rules in terms of retarded $N$-point functions, thereby simplifying earlier approaches and simultaneously solving the issue of non-vanishing vacuum diagrams that has plagued finite temperature expansions. Our approach is moreover valid for nonequilibrium systems in initial equilibrium and allows us to show that important commonly used approximations, namely the $GW$, second Born and $T$-matrix approximation, retain positive spectral functions at finite temperature. Finally we derive an analytic continuation relation between the spectral forms of retarded $N$-point functions and their Matsubara counterparts and a set of Feynman rules to evaluate them.