The traditional approaches of modeling and estimation of highly skewed deposits have led to incorrect evaluations, creating challenges and risks in resource management. The low concentration of the rare earth element (REE) deposits, on one hand, and their strategic importance, on the other, enhances the necessity of multivariate modeling of these deposits. The wide variations of the grades and their relation with different rock units increase the complexities of the modeling of REEs. In this work, the Gazestan Magnetite-Apatite deposit was investigated and modeled using the statistical and geostatistical methods. Light and heavy REEs in apatite minerals are concentrated in the form of fine monazite inclusions. Using 908 assayed samples, 64 elements including light and heavy REEs from drill cores were analyzed. By performing the necessary pre-processing and stepwise factor analysis, and taking into account the threshold of 0.6 in six stages, a mineralization factor including phosphorus with the highest correlation was obtained. Then using a concentration-number fractal analysis on the mineralization factor, REEs were investigated in various rock units such as magnetite-apatite units. Next, using the sequential Gaussian simulation, the distribution of light, heavy, and total REEs and the mineralization factor in various realizations were obtained. Finally, based on the realizations, the analysis of uncertainty in the deposit was performed. All multivariate studies confirm the spatial structure analysis, simulation and analysis of rock units, and relationship of phosphorus with mineralization.