A hybrid model reduction for use with the automated multilevel substructuring (AMLS) method is presented for frequency response analysis of complex structures. Structure responses to harmonic excitations and quasi-static responses to dominant damping forces are included in a reduced approximation subspace. Both types of responses greatly increase the efficiency of the subspace for solving the frequency response problem (FRP) for systems with high modal density and structural damping, and provide a good preparation for future frequency-dependent problems. A distilled subspace assumed to provide accurate frequency responses is generated from the finite element (FE) models by using the AMLS method. Then, the hybrid model reduction method is used to reduce the distilled subspace into a small new subspace. Three types of vectors are used to construct this subspace. The first type is distilled subspace dynamic response vectors (DRVs), which are exact solutions in the distilled subspace at certain chosen frequencies, called the DRV frequencies. The second type is modal DRVs, which are inexpensive approximate solutions calculated in an eigenspace. The third type is damping deformation vectors (DDVs), which provide information about response of the structure to damping effects. As exact responses, the distilled subspace DRVs eliminate frequency response errors at the DRV frequencies, and improve the accuracy at nearby frequencies as well. A small number of DRV frequencies are chosen carefully to offer maximum benefit with minimal computational cost. The modal DRVs are approximated very inexpensively from a suitable eigenspace. Only the diagonal entries in the modal coefficient matrices are used, along with low-rank updates that improve the accuracy of the modal DRVs and are applied using the Sherman-Morrison-Woodbury formula. Because of their low cost, a large number of modal DRVs constitute the major part of the reduced subspace. A small number of DDVs represent response to provide damping with minimal computational cost. The dimension of the final subspace is minimized by removing any redundancy through a special implementation of the QR factorization. This method results in a much smaller new subspace than the one from traditional modal truncation while achieving the same FRP accuracy. Such an efficiency also establishes a good foundation for future application in frequency-dependent problems. / text