SUMMARY The solution of the remote reference method, a frequently used technique in magnetotelluric (MT) data processing, can be viewed as a product of the two-input–multiple-output relationship between the local electromagnetic field and the reference field at a remote station. By applying a robust estimator to the two-input–multiple-output system, one can suppress the influence of outliers in the local magnetic field as well as those in the local electric field based on regression residuals. Therefore, this study develops a new robust remote reference estimator with the aid of robust multivariate linear regression. By applying the robust multivariate regression S-estimator to the multiple-output system, the present work derives a set of equations for robust estimates of the transfer function, noise variances, and the scale of the Mahalanobis distance simultaneously. The noise variances are necessary for the multivariate analysis to normalize the residuals of dependent variables. The Mahalanobis distance, a distance measure for multivariate data, is a commonly used indicator of outliers in multivariate statistics. By updating those robust estimates iteratively, the new robust remote reference estimator seeks the transfer function that minimizes the robust scale estimate of the Mahalanobis distance. The developed estimator can avoid bias in the MT transfer function even if there are significant noises in the reference magnetic field and handle outlying data more robustly than previously proposed robust remote reference estimators. The authors applied the developed method to a synthetic data set and real-world data. The test results demonstrate that the developed method downweights outliers in the local electric and magnetic fields and gives an unbiased transfer function.