[Rd] odd (erroneous?) results from gls

[Uwe Ligges]

Timothy_Handley at nps.gov wrote:
> A couple weeks ago I posted a message on this topic to r-help, the response
> was that this seemed like odd behavior, and that I ought to post it to one
> of the developer lists. I posted to r-sig-mixed-models, but didn't get any
> response. So, with good intentions, I decided to try posting once more, but
> to this more general list.
> 
> The goal is (1) FYI, to make you aware of this issue, in case it is an
> error in gls (2) For my information, in case I have made an error, in the
> hope that one of you folks might be able to correct me. Thanks in advance
> for your time.
> 
> The issue is in 2 parts.
> 
> (A) I've used gls to fit a model with two fixed effects and a corExp
> object. By my count, this fitting process estimates 5 parameters:
> (Intercept), l10area, newx, range, and nugget. With 118 total df, there
> should be 118 - 5 = 113 residual df. However, the output from summary.gls
> reports 115 residual degrees of freedom. Is this an error in summary or
> gls, or is there an error in my count?


Those for the correlation structure are not counted for these residuals 
as you can find in

  nlme:::print.summary.gls

that contains the line

  cat("Degrees of freedom:", dd[["N"]], "total;", dd[["N"]] -
         dd[["p"]], "residual\n")


> (B) Summary.gls reports logLik=-273.6. Using my count of 5 estimated
> parameters, the AIC should be -2*(-273.6) + 2*5 = 557.2. However,
> summary.gls reports an AIC of 559.2. If one works backwards from the
> reported AIC of 559.2, it seems that gls believes it has estimated 6
> parameters in the fitting process. Again, is this an error in gls, or an
> error on my part?


1 additional was used for the estimation of the variance. Accordingly

nlme:::logLik.gls

contains the line

   attr(val, "df") <- p + length(coef(object[["modelStruct"]])) + 1

The AICs should be comparable at least within R and if others think 5 
rather than 6 should be used its still fine since the difference will be 
there in all models to be compared.


Best wishes,
Uwe Ligges



> Copied from R terminal:
> 
>> summary(sppl.i.ex)
> Generalized least squares fit by maximum likelihood
>   Model: all.all.rch ~ l10area + newx
>   Data: gtemp
>       AIC      BIC    logLik
>   559.167 575.7911 -273.5835
> 
> Correlation Structure: Exponential spatial correlation
>  Formula: ~x + y | area
>  Parameter estimate(s):
>      range     nugget
> 15.4448835  0.3741476
> 
> Coefficients:
>                Value Std.Error   t-value p-value
> (Intercept) 7.621306 0.7648135  9.964921  0.0000
> l10area     6.332931 0.5589199 11.330659  0.0000
> newx        0.066535 0.0204417  3.254857  0.0015
> 
>  Correlation:
>         (Intr) l10are
> l10area -0.605
> newx     0.358 -0.024
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max
> -3.0035983 -0.5990432 -0.2226852  0.5113270  2.4444263
> 
> Residual standard error: 2.820337
> Degrees of freedom: 118 total; 115 residual
> 
> Tim Handley
> Fire Effects Monitor
> Santa Monica Mountains National Recreation Area
> 401 W. Hillcrest Dr.
> Thousand Oaks, CA 91360
> 805-370-2347
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel