We investigate model order reduction (MOR) of random parametric linear systems via the regression method. By sampling the random parameters contained in the coefficient matrices of the systems, the iterative rational Krylov algorithm (IRKA) is used to generate sample reduced models corresponding to the sample data.We assemble the resulting reduced models by interpolating the coefficient matrices of reduced sample models with the regression technique, where the generalized polynomial chaos (gPC) is adopted to characterize the random dependence coming from the original systems. Noting the invariance of the transfer function with respect to restricted equivalence transformations, the regression method is conducted based on the controllable canonical form of reduced sample models in such a way to improve the accuracy of reduced models greatly.We also provide a posteriori error bound for the projection reduction method in the stochastic setting. We showcase the efficiency of the proposed approach by two large-scale systems along with random parameters: a synthetic model and a mass-spring-damper system.