P-value for ML estimate of error-autoregressive parameters in R using GLS with ARMA(p=1) correlation

Question

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")

summary(test.reg)

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):
   Phi 
0.01887445 

Coefficients:
                 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

 Correlation: 
  (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
   Phi 
0.01887445 

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

Thanks so much, John

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:

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

require(nlme)
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.