# How to compare GLS models?

I have spatial data set for 35 studies. In each study, there are variables y, x1, x2, latitude, and longitude. I want to know whether adding x2 to model y~x1 will improve the simple regression model y~x1.

I use m_study1 = gls(y ~ x1 + x2, data=study1, correlation=corExp(form=~Lat+Long)) for each study. I also use m_study1_x1 = gls(y ~ x1, data=study1, correlation=corExp(form=~Lat+Long)). I used AIC to compare two models and got an error "models are not all fitted to the same number of observations". I think that is because of different df in two models. What should I do? How can I tell whether adding x2 will improve predicting y? Thx.

Why the number of observations change? You are using the same dataset (study1).

It may because the use of REML...

In addition to the ANOVA cited above, you can try the model.sel function from MuMIn package, which rank by AICc(an AIC corrected for smaller samples, it's safer). You could also change the correlation structures, gaussian, spherical, and check which one gives the lower AICc.

Using a preset database from R programming, you could compare the fitted models using ANOVA:

df <- mtcars

fit1 <- gls(mpg ~ cyl + disp, data = df)
fit2 <- gls(mpg ~ cyl + disp + hp, data = df)

fit1ml <- update(fit1, . ~ ., method = "ML")
fit2ml <- update(fit2, . ~ ., method = "ML")
anova(fit1ml, fit2ml)


Which would return comparisons and p-value:

> anova(fit1ml, fit2ml)

Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit1ml     1  4 167.1456 173.0086 -79.57282
fit2ml     2  5 168.0184 175.3471 -79.00921 1 vs 2 1.127214  0.2884

• Thank you. I used default method="REML". It showed the warning that REML comparisons are not meaningful. I tried "ML", it works. However, since there are 35 studies and I did this anova comparison for each study, do I need to adjust "significance level" by Bonferroni? @Ryan Morton Jun 29 '17 at 14:17