Theoretical calculations of thermally activated decay of skyrmions in systems comprising several magnetic monolayers are presented, with a special focus on bilayer systems. Mechanisms of skyrmion collapse are identified and corresponding energy barriers and thermal collapse rates are evaluated as functions of the interlayer exchange coupling and mutual stacking of the monolayers using transition state theory and an atomistic spin Hamiltonian. In order to contrast the results to monolayer systems, the magnetic interactions within each layer are chosen so as to mimic the well-established Pd/Fe/Ir(111) system. Even bilayer systems demonstrate a rich diversity of skyrmion collapse mechanisms that sometimes coexist. For very weakly coupled layers, the skyrmions in each layer decay successively via radially symmetric shrinking. Slightly larger coupling leads to an asymmetric chimera collapse stabilized by the interlayer exchange. When the interlayer exchange coupling reaches a certain critical value, the skyrmions collapse simultaneously. Interestingly, the overall energy barrier for the skyrmion collapse does not always converge to a multiple of that for a monolayer system in the strongly coupled regime. For a certain stacking of the magnetic layers, the energy barrier as a function of the interlayer exchange coupling features a maximum and then decreases with the coupling strength in the strong coupling regime. Calculated mechanisms of skyrmion collapse are used to ultimately predict the skyrmion lifetime. Our results reveal a comprehensive picture of the thermal stability of skyrmions in magnetic multilayers and provide a perspective for realizing skyrmions with controlled properties.