We consider time-average quantities of chaotic systems and their sensitivity to system parameters. When the parameters are random variables with a prescribed probability density function, the sensitivities are also random. The central aim of the paper is to study and quantify the uncertainty of the sensitivities; this is useful to know in robust design applications. To this end, we couple the nonintrusive polynomial chaos expansion (PCE) with the multiple shooting shadowing (MSS) method, and apply the coupled method to two standard chaotic systems, the Lorenz system and the Kuramoto-Sivashinsky equation. The method leads to accurate results that match well with Monte Carlo simulations (even for low chaos orders, at least for the two systems examined), but it is costly. However, if we apply the concept of shadowing to the system trajectories evaluated at the quadrature integration points of PCE, then the resulting regularization can lead to significant computational savings. We call the new method shadowed PCE (sPCE).