Finding overall p-value for GLS model

Question

Does anyone know how I can find/calculate an overall p-value for a GLS multiple regression model (made with nlme)?

Answer 1

The way one is expected to do this in gls() is to use a likelihood ratio test between two models via the anova() methods for "gls" objects. This is a general way of comparing two nested models with one another, but because of the way these models are fitted (REML estimates are used by default as one is generally fitting variance parameters for the correlation or weights arguments, but to compare models with different fixed effects we need ML estimates), you can't really get the significance of the overall model - it depends what you are wanting to test - the fixed effects or the variance parameters for the variance-covariance matrix?

Here is an example using data from the nlme package.

require(nlme)
fm1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), Ovary,
           correlation = corAR1(form = ~ 1 | Mare))

There are two things to test in that model:

1. The fixed effects terms
2. The need for the correlation term

Lets assume the interest is in the fixed effects terms. We need an appropriate Null model that does not include these terms but does include the intercept and the same correlation estimate from the candidate model above.

fm0 <- gls(follicles ~ 1, Ovary,
           correlation = corAR1(value = 0.7532079, form = ~ 1 | Mare, 
                                fixed = TRUE))

(I got the value by looking at intervals(fm1). Both these models have been fitted using REML, which is not useful in testing models with different fixed effects, however, REML gives unbiased estimates of variance parameters (for the corAR1() bit, hence we used that first.)

To check if the fixed effect terms are significant, we use anova() on the updated models that are fitted using ML:

fm0.ml <- update(fm0, . ~ ., method = "ML")
fm1.ml <- update(fm1, . ~ ., method = "ML")
anova(fm0.ml, fm1.ml)

which gives

> anova(fm0.ml, fm1.ml)
       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fm0.ml     1  2 1588.169 1595.629 -792.0844                        
fm1.ml     2  5 1574.387 1593.037 -782.1934 1 vs 2 19.78195   2e-04

As I used the same estimated correlation structure in both models (based on REML estimate from fm1) the difference between the two is isolated to the two fixed effects. However, the difference in degrees of freedom is not right here; because I fixed the value of the AR(1) parameter, the software doesn't treat that as a degree of freedom used, hence it thinks that the models differ by 3 instead of 2 degrees of freedom.

To address this we could just compute the p-value for the L.Ratio shown using the correct d.f.:

> pchisq(19.78195, 2, lower.tail = FALSE)
[1] 5.062956e-05

Or we could update fm1.ml to use a fixed estimate of the AR(1) parameter and redo the LRT:

fm1.ml <- update(fm1.ml, . ~ ., 
                 correlation = corAR1(value = 0.7532079, form = ~ 1 | Mare, 
                                      fixed = TRUE))

> anova(fm0.ml, fm1.ml)
       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fm0.ml     1  2 1588.169 1595.629 -792.0844                        
fm1.ml     2  4 1572.449 1587.369 -782.2243 1 vs 2 19.72008   1e-04

The reason for the difference in p-values between the two is the slightly lower L.Ratio for fm1.ml when fixing the AR(1) parameter.

If you want to check the need for the correlation structure, we need to use the REML estimated models:

fm2 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), Ovary)
anova(fm2, fm1)

> anova(fm2, fm1)
    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fm2     1  4 1804.868 1819.749 -898.4340                        
fm1     2  5 1571.455 1590.056 -780.7273 1 vs 2 235.4135  <.0001