I would like to know if there is a way to get p-values when using the GLS function (part of the nlme package) for the ML estimate of the error-autoregressive parameters. I looked through the package information on nlme and have not found a clear example of how to do this.

Here is an example of some code that illustrates my question:

x <- rnorm(100)
y <- rnorm(100)
z <- rnorm(100)
df <- data.frame(x,y,z)

test.reg <- gls(x ~ y+z,data=df,correlation=corARMA(p=1),method="ML")


Generalized least squares fit by maximum likelihood
  Model: x ~ y + z 
  Data: df 
   AIC      BIC    logLik
  289.3276 302.3534 -139.6638

Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):

                 Value  Std.Error   t-value p-value
(Intercept) 0.08418063 0.10214246 0.8241492  0.4119
y           0.00236772 0.09410912 0.0251593  0.9800
z           0.07689263 0.09480738 0.8110405  0.4193

  (Intr) y     
y  0.135       
z -0.013  0.074

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-1.94366502 -0.66310726  0.09546043  0.64881478  2.78384648 

Residual standard error: 0.9781187 
Degrees of freedom: 100 total; 97 residual

What I would like to know more about is the Parameter estimate Phi. I have searched through the nlme package to get more information about p-values for Phi but have come up short. I am unsure if the gls function computes the p-values for Phi.

I do know this gives the coefficients of phi:

coeff <- test.reg$modelStruct$corStruct

Correlation structure of class corAR1 representing

But that is all I have really found. Any help would be great!

Thanks so much, John


1 Answer 1


A very general and often easy way to test whether particular parameter(s) are zero is to use a likelihood ratio test (LRT). In many models that are estimated using maximum likelihood, LRTs are a valid way to compare nested models.

The basic steps are:

1) Fit your full model 2) Fit a reduced model constraining the parameter(s) of interest to zero, often achieved by simply removing the term. 3) Compare the likelihoods compared to the change in number of parameters

The equation is:

$$ -2 * (LL_F - LL_R) \sim \chi^2_{(DF_F - DF_R)} $$

where $LL_F$ = log likelihood for the full model, $DF_F$ = degrees of freedom for the full model, and the $R$ subscript indicates the reduced model.

Taking this to R:

d <- data.frame(
  x = rnorm(100),
  y = rnorm(100),
  z = rnorm(100))

reg.f <- gls(x ~ y + z, data = d, correlation = corARMA(p=1), method = "ML")
reg.r <- gls(x ~ y + z, data = d, method = "ML")

anova(reg.f, reg.r)

##       Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## reg.f     1  5 291.3681 304.3940 -140.6841                         
## reg.r     2  4 289.9794 300.4001 -140.9897 1 vs 2 0.6113316  0.4343

the anova function on two models will compare the log likelihoods and compute the LRT.

  • 1
    $\begingroup$ Thanks for the answer! I would upvote it if I had enough reputation. $\endgroup$ Feb 17, 2013 at 14:50
  • $\begingroup$ @JohnLombardi Quite welcome! Also, I caught a typo in my formula, which is now fixed. $\endgroup$
    – Joshua
    Feb 19, 2013 at 2:02

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.