I'm surprised that no one has yet mentionned the Christoffersen, Heston and Jacobs paper entitled "The Shape and Term Structure of the Index Option Smirk: Why Multifactor Stochastic Volatility Models Work so Well" published in 2009 in Management Science. You can find it here.
They actually look into the evolution of implied volatility surface of (European) option contracts written on the S\&P500. They filter out "odd" contracts as is routinely done in the literature (basically, they apply the Bakshi, Cao and Chen (1997) filters and other consistency checks). Then, they fit a quadratic polynomial in both moneyness and maturity on the implied volatility surface by ordinary least squares. The idea is that this gives you a statistically sound way of spanning the whole surface. Once this is done, they generate a fitted value surface for a large range of maturity and moneyness that they want to look into. Finally, they apply principal component analysis to extract latent factors from this surface and find that 2 factors explain most of the variance present in the data across time.
If you dig into their paper, you will find that they perform regressions in the moneyness space with these two factors and it seems that one of the factor captures the shifting "levels" of the smile, while the other captures the "slope." So, if your interest lies in thinking about the evolution of the (risk-neutral) volatility, what this tells you is that at least two sources of variations are required in the volatility process of your model because your smile jumps up and down, as well moves between being more and less "smirky" for a given maturity.
Now, going back to your concern about volatility and skewness, your hesitation might stem from missing an important detail. We can produce realtively reliable estimates of stock market volatility using high frequency data -- that is, we can get estimated time series of volatility under the physical measure. Neglecting concerns over jumps, say you believe realized volatility is good enough. Well, that's negatively correlated with returns on the index. The same thing would happen if you used a model to filter out a conditional volatility time series under the physical measure. In other words, under the physical measure, you have a negative correlation. You would expect that for most maturities and moneyness, you'd find the same thing under the risk-neutral measure.