The anisotropic shear relaxation functions of confined thin films of unentangled polymer melts are measured via nonequilibrium step strain simulations of in-plane and out-of-plane shear using the finitely extensible, nonlinear-elastic (FENE) model. We show that the classical Rouse model unsurprisingly fails to predict the thin-film relaxation functions in response to out-of-plane shear, due in part to non-Gaussian conformation statistics in the dimension perpendicular to the sub/superstrate. Using an alternate empirical model for the out-of-plane response, we quantify decreases in the plateau modulus G(perpendicular to)(P) relaxation time lambda(perpendicular to), and viscosity eta(perpendicular to) and an increase in the logarithmic relaxation rate r(perpendicular to) as functions of film thickness, and we discuss these anisotropic changes in stress-relaxation properties in terms of structural/conformation changes on the microscopic level, namely the relative contraction and non-Gaussian quality of polymer conformations in the dimension normal to the substrate and the resulting phenomenon of cooperative relaxation. We then incorporate these into a semiempirical extension to the Rouse model which closely predicts our computational results and which will be useful for further study of polymer thin films.