Comparing multiple gls' objects using anova.gls I got a warning message when I was trying to do anova for two nlme::gls objects. Here is an example:
require(nlme)
set.seed(123)
y<-rnorm(100,10,2)
x1<-rnorm(100)
x2<-sample(1:5,100,T)
x3<-rt(100,20)
x4<-rbinom(100,1,0.3)
fit1<-gls(y~x1+x2+x3+x4,correlation=corAR1(form=~1|x4))
fit2<-gls(y~x3+x4,correlation=corAR1(form=~1|x4))
anova(fit1,fit2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit1     1  7 421.5587 439.4359 -203.7794                        
fit2     2  5 413.9855 426.8590 -201.9928 1 vs 2 3.573242  0.1675
Warning message:
In nlme::anova.lme(object = fit1, fit2) :
  fitted objects with different fixed effects. REML comparisons are not meaningful.

Dose anybody know what the message means? Should I do the following rather than above?
> anova(fit1, L=c(-1, -1, 1, 1))
Denom. DF: 95 
 F-test for linear combination(s)
(Intercept)          x1          x2          x3 
         -1          -1           1           1 
  numDF  F-value p-value
1     1 281.9844  <.0001

Why the numDF is one here? It seems 2?
 A: The two GLS models are fitted by maximising a restricted (or residual) maximum likelihood function, the REML in the warning. This approach gives better estimates of the variance components of the model; the correlation=corAR1(form=~1|x4) part. However, REML estimates are not appropriate for comparing models where the fixed effects parts differ. The fixed effects parts of your two models are
x1+x2+x3+x4

and
x3+x4

You should re-estimate the two models using maximum likelihood, and then refit the chosen model using the default REML to get the best estimate for the covariance matrix of the model. This is quite easily done, for example:
anova(update(fit1, . ~ ., method = "ML"),
      update(fit2, . ~ ., method = "ML"))

The update() calls will refit the models in-place and throw away the objects once anova() is done with them. The . ~ . part just indicates that no changes should be made to the model response or fixed effects. I forget if this is even need. If you want the ML models for later, do
fit1ml <- update(fit1, . ~ ., method = "ML")
fit2ml <- update(fit2, . ~ ., method = "ML")
anova(fit1ml, fit2ml)

A: The L= c(-1, -1, 1, 1) specifies the contrast you would like to test; as a result, you are testing for $\beta_0 + \beta_1= \beta_2 + \beta_3$ by specifying the command.
You can test $\beta_1= 0$ and $\beta_2= 0$ using F-tests by issuing the following command:
L<- matrix(0, nrow= 2, ncol= 5, dimnames= list(NULL, names(coef(fit1))))
L[1, "x1"]<- 1
L[2, "x2"]<- 1
L
anova(fit1, L= L)

You will obtain a similar result as in 
anova(fit1, fit2)

Hope it helps.
