A hierarchical mixed-effects model is proposed to account for both individual- and population-level variability in human immunodeficiency virus (HIV) dynamics. This model is implemented by formulating the crucial parameters as random variables in an in-host HIV model. Model reduction is used to guide the choice for a minimal set of parameters, whose distributions are estimated by the global two-stage method. We analyze the system of ordinary differential equations with random coefficients and provide numerical simulations illustrating its asymptotic behaviors.