In this paper, we propose a data-driven model reduction method to solve parabolic inverse source problems with uncertain data efficiently. Our method consists of offline and online stages. In the offline stage, we explore the low-dimensional structures in the solution space of parabolic partial differential equations (PDEs) in the forward problems with a given class of source functions and construct a small number of proper orthogonal decomposition (POD) basis functions to achieve significant dimension reduction. Equipped with the POD basis functions, we can solve the forward problems extremely fast in the online stage. Thus, we develop a fast algorithm to solve the optimization problem in parabolic inverse source problems, which is referred to as the POD method. Moreover, we design an a posteriori algorithm to find the optimal regularization parameter in the optimization problem using the proposed POD method without knowing the noise level. Under a weak regularity assumption on the solution of the parabolic PDEs, we prove the convergence of the POD method in solving the forward parabolic PDEs. In addition, we obtain the error estimate of the POD method for parabolic inverse source problems. Finally, we present numerical examples to demonstrate the accuracy and efficiency of the proposed method. Numerical results show that the POD method provides considerable computational savings over the finite element method while maintaining the same accuracy.