# Fitting a generalized least squares model with correlated data; use ML or REML?

Reading the Linear Mixed Model (LMM) literature I am aware that fitting a model using REML provides better estimates of variance parameters than fitting via ML. However, we should not compare nested models fitted with REML that have different fixed effects.

Recently, I have been fitting some models using GLS via the gls() function in the nlme package for R. The default fitting method for that function is REML. Do the same principals of REML vs ML for LMM also apply to GLS?

Specifically, I am fitting a model with and without a linear trend with a correlation structure in the residuals:

m1 <- gls(Response ~ Time, data = foo, correlation = corAR1(form ~ Time))
m0 <- gls(Response ~ 1, data = foo, correlation = corAR1(form ~ Time))


In the above, I should fit the models using ML as they have different fixed effects. Is this correct?

Secondly, consider two GLS models that differ in the correlation structure:

m1 <- gls(Response ~ Time, data = foo, correlation = corARMA(form ~ Time, p = 1))
m2 <- gls(Response ~ Time, data = foo, correlation = corARMA(form ~ Time, p = 2))


What fitting method should ideally be used here? REML or ML? Here my intuition would say fit via REML as we are estimating (co)variance parameters. Is my intuition correct or have I got this all mixed up?

-

Your intuition is correct, the same principles apply. I looked in Pinheiro/Bates section 5.4, where gls is introduced, but it doesn't say so explicitly, so you'll just have to trust me, I guess. :)