[R] GAM and AIC: How can I do??? please

Martin Maechler maechler at stat.math.ethz.ch
Mon Oct 24 09:49:00 CEST 2005


>>>>> "ronggui" == ronggui  <042045003 at fudan.edu.cn>
>>>>>     on Mon, 24 Oct 2005 10:09:30 +0800 writes:

    ronggui> ======= 2005-10-24 09:55:32
    ronggui> úÚ´ÅÐ´Àº=======

    >>  Hello, I'm a Korean researcher who have been started to
    >> learn the "R" package.
    >> 
    >> I want to make gam model and AIC value of the model to
    >> compare several models.
    >> 
    >> I did the GAM model, but there were error for AIC.
    >> SO, how can I do? pleas help me!!!
    >> 
    >> 
    >> I did like below;

    >> 
    >> > a.fit <- gam(pi~ s(t1r), family = gaussian(link="log"))
    >> > summary(a.fit)

    >> >   Family: gaussian
    >> >   Link function: log
    >> >
    >> >   Formula:
    >> >   pi ~ s(t1r)
    >> >
    >> >   Parametric coefficients:
    >> >              Estimate  std. err.    t ratio    Pr(>|t|)
    >> >   constant   0.093105   0.005238      17.77    < 2.22e-16
    >> >
    >> >   Approximate significance of smooth terms:
    >> >                 edf       chi.sq     p-value
    >> >   s(t1r)      1.833       24.153     0.00014213
    >> >
    >> >   R-sq.(adj) =  0.435   Deviance explained = 47.1%
    >> >   GCV score = 0.0010938   Scale est. = 0.00099053  n = 30

    ronggui> are you using the mgcv package?  if you are,just
    ronggui> use a.fit$aic to get the aic.

hmm, yes, and no:

It's true what you say,  
BUT is not at all recommended in general:

You should use the generic AIC() function
rather than extracting components yourself.

This is a general priniciple:  If possible use  'extractor functions'
to work on objects rather then relying on internal
representations.

This is particularly relevant for fitted models:

Do use  residuals(.), fitted(.), LogLik(.), AIC(.), vcov(.) 
etc etc!


Now back to this problem:

    >> AIC(a.fit) 
    >> Error in logLik(object) : no applicable method for "logLik"

I can't reproduce this; Eun definitely needs to give more
details, since the following works fine:

> library(mgcv)
> x <- 1:50
> set.seed(1)
> y <- 2^(sin(x/10) + rnorm(50))
> a.fit <- gam(y ~ s(x), family = gaussian(link="log"))
> summary(a.fit)

Family: gaussian 
Link function: log 

Formula:
y ~ s(x)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.3171     0.1251   2.535   0.0147 *
---
Signif. codes:  ..........{UTF-8 code}

Approximate significance of smooth terms:
       edf Est.rank    F p-value   
s(x) 2.858    9.000 3.07 0.00576 **
---
Signif. codes: ............

R-sq.(adj) =    0.4   Deviance explained = 43.5%
GCV score = 0.94391   Scale est. = 0.87107   n = 50

> AIC(a.fit)
[1] 140.6937
>




More information about the R-help mailing list