The purpose of the present work is to construct estimators for the random effects in a fractional diffusion model using a hybrid estimation method where we combine parametric and nonparametric techniques. We precisely consider n stochastic processes \(\left\{ X_t^j,\ 0\le t\le T\right\} \) , \(j=1,\ldots , n\) continuously observed over the time interval [0, T], where the dynamics of each process are described by fractional stochastic differential equations with drifts depending on random effects. We first construct a parametric estimator for random effects using maximum likelihood estimation techniques and study its asymptotic properties when the time horizon T is sufficiently large. Then, on the basis of the obtained estimator for the random effects, we build a nonparametric estimator for their common unknown density function using Bernstein polynomials approximation. Some asymptotic properties of the density estimator, such as its asymptotic bias, variance, and mean integrated squared error, are studied for an infinite time horizon T and a fixed sample size n. The asymptotic normality of the estimator is established for a fixed T, a high frequency, and as long as the order of Bernstein polynomials is sufficiently large. We also investigate a non-asymptotic bound for the expected uniform error between the density function and its estimator. A numerical study is then presented in order to evaluate both qualitative and quantitative performance of the Bernstein estimator compared with the standard kernel estimator within and at boundaries of the support of the density function.