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

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

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.

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