Blazars are a subclass of active galactic nuclei (AGNs) with a relativistically boosted jet pointing towards the Earth. They are one of the most violent and variable objects in the Universe, showing extreme variability across the entire electromagnetic spectrum. The jet generally dominates the optical emission of these sources. However, other components such as the stellar emission from the host galaxy or the accretion disk may contribute significantly to this emission. Disentangling the different contributions to the optical emission of these sources is challenging due to the high dominance of the jet, but is crucial to study and understand the variability detected in these objects. In this work we will present the results of the spectral variability analysis performed in a sample of gamma-ray blazars monitored by the Steward Observatory from 2008 to 2018. We make use of the non-negative matrix factorization (NMF), a statistical tool of dimensionality reduction to decompose the spectra of each source in the minimum number of components required to explain its spectral variability. We propose an approach based on a meaningful association of each component with the different expected emitting regions in the AGN, namely the jet, broad line region, accretion disk and stellar population. Using this decomposition we study the contribution of each component, and its corresponding physical association, to the flux evolution and the overall variability of the AGN, together with a possible interpretation of the different variability signatures observed in each blazar.