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 thechniques. We precisely consider $n$ stochastic processes $\left\{X_t^j,\ 0\leq t\leq 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 the random effects using the techniques of maximum likelihood estimation and we study its asymptotic properties when the time horizon $T$ is sufficiently large. Then by taking into account 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 and the uniform convergence of the estimator are investigated for an infinite time horizon $T$, a high frequency and as the order of Bernstein polynomials is sufficiently large. Some numerical simulations are also presented to illustrate the performance of the Bernstein polynomials based estimator compared to standard Kernel estimator for the random effects density function.