Comparing the goodness of fit for multiple GLS models applied to different data sets

I fitted simple linear models to several (6) data sets of different sizes, using GLS. All models have the same predictors (but the coefficients are different). I would like to compare the goodness of fit between these models. In a simple OLS regression, I would have simply compared the R-squared scores for each model to see which has the best fit. However, since data points in each data set are related, I have to use GLS.
I believe comparing AIC (or AICc or BIC) scores is not valid here, since the data sets are different. Maybe some measure of the percentage of variance explained by each model? But how exactly is that calculated?
I am using the gls() function from the nlme package in R, so I'd be happy to get relevant code, but general ideas are also welcome.

Some pseudo $$R^2$$ measures are discussed here. These require one to define the effective sample size $$m$$ to use in expressions such as $$1 - \exp(LR / m)$$ where $$LR$$ is the overall likelihood ratio $$\chi^2$$ statistic for the fitted model. At least for the AR1 correlation structure in GLS the effective sample size has been worked out here. I hope that someone can develop this further.
The effective sample size $$m$$ is a very useful concept with longitudinal data. One way to define it is to say that for the standard error (SE) of a parameter estimate of interest other than time, the SE from the GLS fit on the $$n$$ subjects is equivalent to the standard error from $$m$$ subjects each having only one measurement.