In this paper we show how to apply classical probabilistic tools for partial sums $\sum_{j=0}^{n-1}\varphi\circτ^j$ generated by a skew product $τ$, built over a sufficiently well mixing base map and a random expanding dynamical system. Under certain regularity assumptions on the observable $\varphi$, we obtain a central limit theorem (CLT) with rates, a functional CLT, an almost sure invariance principle (ASIP), a moderate deviations principle, several exponential concentration inequalities and Rosenthal type moment estimates for skew products with $α, φ$ or $ψ$ mixing base maps and expanding on the average random fiber maps. All of the results are new even in the uniformly expanding case. The main novelty here (contrary to \cite{ANV}) is that the random maps are not independent, they do not preserve the same measure and the observable $\varphi$ depends also on the base space. For stretched exponentially $\al$-mixing base maps our proofs are based on multiple correlation estimates, which make the classical method of cumulants applicable. For $φ$ or $ψ$ mixing base maps, we obtain an ASIP and maximal and concentration inequalities by establishing an $L^\infty$ convergence of the iterates $\cK^n$ of a certain transfer operator $\cK$ with respect to a certain sub-$\sig$-algebra, which yields an appropriate (reverse) martingale-coboundary decomposition.