We propose a formulation to obtain the line shape of a magnetic response with dissipative effects that directly reflects the nature of the environment. Making use of the fact that the time evolution of a response function is described by the same equation as the reduced density operator, we formulate a full description of the complex susceptibility. We describe the dynamics using the equation of motion for the reduced density operator, including the term for the initial correlation between the system and a thermal bath. In this formalism, we treat the full description of non-Markovian dynamics including the initial correlation. We present an explicit and compact formula up to the second order of cumulants, which can be applied in a straightforward way to multiple-spin systems. We also take into account the frequency shift by the system-bath interaction. We study the dependence of the line shape on the type of interaction between the system and the thermal bath. We demonstrate that the present formalism is a powerful tool for investigating various kinds of systems and we show how it is applied to spin systems, including those with up to three spins. We distinguish the contributions of the initial correlation and the frequency shift and make clear the role of each contribution in the Ohmic coupling spectral function. As examples of applications to multispin systems, we obtain the dependence of the line shape on the spatial orientation in relation to the direction of the static field (Nagata-Tazuke effect), including the effects of the thermal environment, in a two-spin system, along with the dependence on the arrangement of a triangle in a three-spin system.