Does anyone know how I can find/calculate an overall p-value for a GLS multiple regression model (made with nlme)?
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
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:
- The fixed effects terms
- 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)
> 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)  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
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