Mean field approximation is a powerful technique which has been used in many settings to study large-scale stochastic systems. In the case of two-timescale systems, the approximation is obtained by a combination of scaling arguments and the use of the averaging principle. This paper analyzes the approximation error of this `average' mean field model for a two-timescale model $(\boldsymbol{X}, \boldsymbol{Y})$, where the slow component $\boldsymbol{X}$ describes a population of interacting particles which is fully coupled with a rapidly changing environment $\boldsymbol{Y}$. The model is parametrized by a scaling factor $N$, e.g. the population size, which as $N$ gets large decreases the jump size of the slow component in contrast to the unchanged dynamics of the fast component. We show that under relatively mild conditions, the `average' mean field approximation has a bias of order $O(1/N)$ compared to $\mathbb{E}[\boldsymbol{X}]$. This holds true under any continuous performance metric in the transient regime, as well as for the steady-state if the model is exponentially stable. To go one step further, we derive a bias correction term for the steady-state, from which we define a new approximation called the refined `average' mean field approximation whose bias is of order $O(1/N^2)$. This refined `average' mean field approximation allows computing an accurate approximation even for small scaling factors, i.e., $N\approx 10 -50$. We illustrate the developed framework and accuracy results through an application to a random access CSMA model.