The dynamic analysis of nonlinear viscoelastic systems in the frequency domain is not an easy task. In most cases, it is due to the frequency- and temperature-dependent properties of the viscoelastic part. Additionally, due to the inherent uncertainties affecting the viscoelastic efficiency in practical situations, their handling in the nonlinear modeling methodology becomes essential nowadays. However, it is still an issue. Thus, this paper presents a numerical modeling methodology intended to perform dynamic analyses in the frequency domain of thin sandwich plates under large displacements. The uncertainties characterizing the nonlinear dynamics of the viscoelastic system are introduced on the random linear and nonlinear finite element matrices by performing the Karhunen–Loève expansion technique. The Latin hypercube sampling method is used herein as the stochastic solver, and the nonlinear frequency responses are computed using the harmonic balance method combined with the Galerkin bases. To overcome the difficulty in solving the resulting complex nonlinear eigenproblem with a frequency-dependent viscoelastic stiffness, making the stochastic nonlinear analyses in the frequency domain very costly, sometimes unfeasible, an efficient and accurate iterative reduction method is proposed to approximate the complex eigenmodes. The envelopes of nonlinear frequency responses demonstrate clearly the relevance of considering the uncertainties in design variables of viscoelastic systems having nonlinear behavior to deal with more realistic situations.