Abstract Multilayer fibre reinforced composites exhibit significant spatial variabilities in their response due to the variations in the individual lamina properties. Uncertainty quantification of these structures can be carried out using the polynomial chaos based stochastic finite element method (SFEM). The variations in the individual laminae properties are modelled as independent random fields. The SFEM method involves two stages of discretisation. First, the displacement fields are discretised into variables along the spatial dimension using traditional FEM. Next, the random fields are discretised into a vector of correlated random variables using polynomial chaos expansions. As the random field discretisation is carried out for each individual laminae, the number of random variables entering the formulation becomes large, making uncertainty quantification computationally prohibitive. This study focusses on the development of model order reduction strategies that enable reducing the stochastic dimensionality of the problem and enable faster computations. This is achieved by developing a probabilistically equivalent structure scale model using two approaches – one based on probabilistic equivalence of a nodal response and the other seeking equivalence on the probabilistic characteristics of the constitutive matrix. A set of numerical examples are presented to highlight the salient features of the proposed developments.