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.