1
$\begingroup$

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

$\endgroup$
  • 5
    $\begingroup$ Welcome to SO. Your question is actually better suited for crossvalidated.com . Probably one of the mods will move it there soon enough. Next to that, I'd advise you to add some toy dataset and a reproducible example so people can play around a bit as well. It shows you did some effort, and people will be more keen to help you. On how to make a reproducible example, see this question $\endgroup$ – Joris Meys Aug 4 '11 at 15:36
6
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Gavin - very nice explanation and examples, with code. At first I was confused by the need to fix the correlation parameter $\rho$, but then I thought of an analogy with getting a likelihood ratio test for one of the regression parameters in the OLS univariate linear model case. Do you fix $\sigma^2$ when comparing the two models, or do you let if float? [I'm ignoring the fact that we have an exact $F$ test in this case.] $\endgroup$ – Frank Harrell Aug 5 '11 at 12:50
  • $\begingroup$ @Frank Thanks. My first stab allowed $\rho$ to be re-estimated, but then I realised that the LRT would be comparing models with not just different fixed effects but different correlation structures. $\sigma^2$ was re-estimated. I'm not aware of a canned way with gls() to allow fixing that. I came at this from a practical point of view - re your comment in [...], what is the exact $F$? Showing my limited statistical knowledge now I'm afraid. $\endgroup$ – Reinstate Monica - G. Simpson Aug 5 '11 at 13:06
  • $\begingroup$ You're doing fine - this is a tricky issue that I wish I had figured out years ago. Someone can probably point us to a good reference about what parameters to "freeze" when doing likelihood ratio tests, in general. $\endgroup$ – Frank Harrell Aug 5 '11 at 18:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy