We consider statistical inference for a class of mixed-effects models with system noise described by a non-Gaussian integrated Ornstein-Uhlenbeck process. Under the asymptotics where the number of individuals goes to infinity with possibly unbalanced sampling frequency across individuals, we prove some theoretical properties of the Gaussian quasi-likelihood function, followed by the asymptotic normality and the tail-probability estimate of the associated estimator. In addition to the joint inference, we propose and investigate the three-stage inference strategy, revealing that they are first-order equivalent while quantitatively different in the second-order terms. Numerical experiments are given to illustrate the theoretical results.