Abstract An RPA type approach to nuclear rotation and its coupling to other modes is proposed, giving both energies and transition rates in deformed nuclei. It is based on a general expansion of any multipole operator, in terms of earlier introduced elementary transition operators, acting in and between nuclear rotational bands. The expansion is formally obtained by using symmetry properties only. The matrix element of any expansion term is calculated from the algebraic properties of the transition operators. Applied to the Hamiltonian it gives a phenomenological model, accounting for vibrations and rotation with coupling between these modes, together with the prescription for evaluation of the model parameters. Applied to a multipole transition inducing field, it leads to a model-independent derivation of the spin ( I)-dependent corrections to adiabatic unified model transition rates. A microscopic procedure for calculating the model-dependent expansion coefficients is proposed. The latter are expressed in terms of the elementary transition operators, or in terms of an extended density matrix involving these operators, giving their quasiparticle structure. The procedure is based on a generalized equations-of-motion method for the density matrix. The method is illustrated by the example of coupled ground-state (g), beta (β)-state and gamma (γ)-state rotational bands in even-even nuclei. Results for E2 non-adiabatic transition rate corrections inside the g-band and for transitions from the γ- and β-bands to the g-band are obtained. This example is to be developed in the following paper towards numerical applications and an extension of the description to high-spin states.