There are mainly three kinds of regression models: Ordinary Least Squares (OLS) which ignores any spatial effect, Spatial AutoRegression (SAR) and Geographically Weighted Regression (GWR) which take into consideration the spatial effect. OLS is based on certain assumptions, among them the independence of residuals. So before applying OLS for prediction we should check if residuals are independent or not using either graphical tools or statistical tests (for example Moran's I) and for spatial stationarity (for example Breusch-Pagan test). If the independence of residuals and spatial stationarity are accepted then OLS can be used confidently. Otherwise, SAR should be used. In this case, residuals must be independent; however, ther may be spatial non-stationarity. In this case, GWR should be considered for taking into account both dependence and nonstationarity. GWR uses distance matrices for considering explicitly the dependence and these matrices are local and change in space for considering nonstationarity. So, I think in GWR we dont have to consider SAC since it is already included in it.

Do you mean a spatially autocorrelated error or a spatial lag? You've posted this as geostatistical, so I'm assuming you want an autocorrelated error. A GWR model with a spatial error is overspecified. A GWR model is a model with random intercept and random slopes. A spatial error model is essentially a random intercept model. So if you combine them, you'd be adding two random intercepts. If they are identified, it would be solely through whatever you assumed about the functional forms of the kernel and the autocovariance function, rather than identified by the data. So, while you could mathematically do this, you should be concerned about getting robust results (which you should be worried about with GWR already). Paez, Uchida and Miyamoto, Environment and Planning A, 34(5) discuss this model you've asked about and show how it can be used to choose between either GWR and Spatial autocorrelation.

I am not familiar with the terms you are using but drawing a wild analogy with what we do in geostatistics applied to reservoir inversion or optimisation I would say the following. If your first regression process assume a single variogram structure and allows you to fit the data to this single structure, recursively applying a linear fit can make sense because there could be several variogram structures. In that case analyzing the residual of the first treatment might reveal the second structure. On the other hand, it might be better to devise an analysis scheme that allows you to access the 2 structures directly.

Theoretically, because it's a local regression applied in moving windows of variable sizes, the GWR should account for spatial "anomalies". Therefore the residuals shouldn’t be spatially autocorrelated. If this is not the case, you could interpolate the residuals and add them to the GWR model, as in regression-kriging approach, in order to improve you model. It is probably better to interpolate the residuals achieved for each local regression model (within the bandwidth), but I don't see how you could do this. So if a single model of residuals can be achieved for the entire region, you could use it if it improves the accuracy of results.

There are mainly three kinds of regression models: Ordinary Least Squares (OLS) which ignores any spatial effect, Spatial AutoRegression (SAR) and Geographically Weighted Regression (GWR) which take into consideration the spatial effect. OLS is based on certain assumptions, among them the independence of residuals. So before applying OLS for prediction we should check if residuals are independent or not using either graphical tools or statistical tests (for example Moran's I) and for spatial stationarity (for example Breusch-Pagan test). If the independence of residuals and spatial stationarity are accepted then OLS can be used confidently. Otherwise, SAR should be used. In this case, residuals must be independent; however, ther may be spatial non-stationarity. In this case, GWR should be considered for taking into account both dependence and nonstationarity. GWR uses distance matrices for considering explicitly the dependence and these matrices are local and change in space for considering nonstationarity. So, I think in GWR we dont have to consider SAC since it is already included in it.

agree with Nicholas and Ahmed, SAC is built into the GWR model.

GWR is an excellent tool for testing model specification (where does model fit look good, where are you missing something?), identifying nonstationarity in relationships, and residuals are generally lower and not spatially autocorrelated.

I'm preparing to run similar models and have found the following book to be helpful.

Geographically Weighted Regression: the analysis of spatially varying relationships. 2002. Fotheringham, A.S. Brunsdon, C. Charlton, M.

Late to the party but I agree with the above.

Further more, if you see spatially correlated differences in residuals from a GWR, this would suggest to me that some additional term is missing from the model which is only relevant locally to a few places. Spatial autocorrelation should be accounted for by the GWR.

That said, the precise choice of kernel is critical to how well autocorrelation is managed not just in terms of band width (i.e. meters or N neighbours) but also as to whether you want to consider distance itself as continuous and Euclidean, as a value surface or perhaps containing barriers. I noted from GISRUK this year that network based regression models are also emerging now which might be more relevant to disease transmission?

Geographically weighted regression (as well as ordinary regression) essentially estimate the non-constant mean term, thus it would make more sense to apply ordinary kriging to the residuals. This and similar ideas have been treated by a number of authors, do a search on Google for "geographically weighted regression kriging" and you will find a number of references

GWR constructs a separate equation for every feature in the dataset incorporating the dependent and explanatory variables of features falling within the bandwidth of each target feature. The shape and extent of the bandwidth is dependent on user input for the Kernel type, Bandwidth method, Distance, and Number of neighbors parameters with one restriction: when the number of neighboring features would exceed 1000, only the closest 1000 are incorporated into each local equation.

GWR should be applied to datasets with several hundred features for best results. It is not an appropriate method for small datasets. The tool does not work with multipoint data.

@Mohammad

what do you mean by multipoint data? Does that mean one or more observations at the same location?

Why Not? You can model spatial dependency even in the presence of Spatial Heterogeneity!! 

Read Fotheringham's ‘‘The Problem of Spatial Autocorrelation’’

and Local Spatial Statistics" published in Geographical Analysis

Could you please tell me how you added this SAC term to your linear regression models? I mean what exactly was this SAC term, how was it defined?

You should look at a paper that appeared in Mathematical Geosciences several years ago, I think one of the authors wrote the R packag for GWR

If you consider the data as a non-random sample from a random function of the form

Z(s) = m(s) + Y(s)

where  m(s) is the non-constant expected value of Z(s) and Y(s) is an intrinsic or second order stationary random function. Then the GWR is an estimator of m(s) hence you wouldn't add a spatial autocorrelation term in the GWR but rather in the kriging equations to estimate/predict the values of Y(s) at non data locations.

As you know GWR is a somewhat more sophisticated version of a Trend surface

I suggest to look at Ruben Fernandez Casal PhD thesis, where he gives a brilliant account on this problem.