regression with Moran eigenvectors for multiple years of data

Hello,

I think this is mostly a statistics question with possibly some R details.
Any feedback is appreciated.

I have several years of spatial biological sampling data in the same region
but the number and locations of sites vary across year.  Very strong
spatial autocorrelation is present in the data.

I want to construct a regression model using Moran' eigenvectors as
explanatory variables to account for SAC. For example,

y_ijk=intercept+x1_ijk+x2_ijk+ EV_k

where x1,x2 are environmental covariates and EV are Moran eigenvectors; i,j
are location and k is year.
Environmental covariate relationships with response variable are assumed
constant across years.

My plan was to first estimate using all years of data:
y=intercept+x1+x2
then use function ME in spdep to find identify Moran eigenvectors to reduce
residual SAC using a year specific (index k) spatial weights and
year-specific residuals using function ME from spdep package:
EV_k=  ME(residuals_k~1, listw=weights_k),
then linearly combine resulting eigenvectors for a given year into a single
vector and then concatenate each year's vector such that the final Moran
eignevector used in the regression is
EV= c(EV_2014,EV_2015,EV_2016)

and add EV as an offset or covariate as in the first equation shown.

This approach seems to work quite well (eliminates residual SAC, doesn't
shift regression coefficients substantially, improves model fit), I just
don't know if it is statistically sound?

thanks!