A method for fitting regression models to data that exhibit spatial correlation and Heteroskedasticity is proposed. A combination of parametric and nonparametric regression techniques is used to iteratively estimate the various components of the model. The approach is demonstrated on a large dataset of predicted nitrogen runoff from agricultural lands in the Midwest and Northern Plains regions of the U.S. For this dataset, the model is comprised of three main components: (1) the mean function which includes farming practice variables, local soil and climate characteristics and the nitrogen application treatment, is assumed linear in the parameters and fitted by generalized least squares. (2) the variance function, which contains a local as well as a spatial component whose shapes are left unspecified, is estimated by local as well as a spatial component whose shapes are left unspecified, is estimated by local linear regression, and (3) the spatial correlation function is estimated by fitting parametric variogram model to the standardized residuals and after adjusting the variogram for the presence of Heteroskedasticity. The fitting of these three components is iterated until convergence. The model provides and improved fit to the data compared to a previous model that ignored the Heteroskedasticity and the spatial correlation.