Is there an R implementation of a scheme for automatic smoothing
parameter selection with loess, e.g., by minimizing one of the AIC/GCV
statistics discussed by Hurvich, Simonoff & Tsai (1998)?
Below is a function that calculates the relevant values of AICC,
AICC1 and GCV--- I think, because I to guess from the names of the
components returned in a loess object.
I guess I could use optimize(), or do a simple line search on span=,
but I'm not sure how to use loess.aic to write a function
that would act as a wrapper for loess() and return the mimimizing
loess fit for a specified criterion.
loess.aic <- function (x) {
# extract values from loess object
if (!(inherits(x,"loess"))) stop("Error: argument must be a
loess object")
span <- x$pars$span
n <- x$n
traceL <- x$trace.hat
sigma2 <- sum( x$residuals^2 ) / (n-1)
delta1 <- x$one.delta
delta2 <- x$two.delta
enp <- x$enp
aicc <- log(sigma2) + 1 + 2* (2*(traceL+1)) / (n-traceL-2)
aicc1<- n*log(sigma2) + n* (
(delta1/(delta2*(n+enp)))/(delta1^2/delta2)-2 )
gcv <- n*sigma2 / (n-traceL)^2
result <- list(span=span, aicc=aicc, aicc1=aicc1, gcv=gcv)
return(result)
}
> cars.lo <- loess(dist ~ speed, cars)
>
> (values <- loess.aic(cars.lo))
$span
[1] 0.75
$aicc
[1] 6.93678
$aicc1
[1] 167.7267
$gcv
[1] 5.275487
>
--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT M3J 1P3 CANADA
On Thu, 17 Nov 2005, Michael Friendly wrote:> Is there an R implementation of a scheme for automatic smoothing > parameter selection with loess, e.g., by minimizing one of the AIC/GCV > statistics discussed by Hurvich, Simonoff & Tsai (1998)?If you particularly want loess smoothing then I don't know, but if penalised spline smoothing will do then in gam() in the mgcv package does minimize GCV. -thomas> Below is a function that calculates the relevant values of AICC, > AICC1 and GCV--- I think, because I to guess from the names of the > components returned in a loess object. > > I guess I could use optimize(), or do a simple line search on span=, > but I'm not sure how to use loess.aic to write a function > that would act as a wrapper for loess() and return the mimimizing > loess fit for a specified criterion. > > loess.aic <- function (x) { > # extract values from loess object > if (!(inherits(x,"loess"))) stop("Error: argument must be a loess object") > span <- x$pars$span > n <- x$n > traceL <- x$trace.hat > sigma2 <- sum( x$residuals^2 ) / (n-1) > delta1 <- x$one.delta > delta2 <- x$two.delta > enp <- x$enp > > aicc <- log(sigma2) + 1 + 2* (2*(traceL+1)) / (n-traceL-2) > aicc1<- n*log(sigma2) + n* ( > (delta1/(delta2*(n+enp)))/(delta1^2/delta2)-2 ) > gcv <- n*sigma2 / (n-traceL)^2 > > result <- list(span=span, aicc=aicc, aicc1=aicc1, gcv=gcv) > return(result) > } > > > > cars.lo <- loess(dist ~ speed, cars) > > > > (values <- loess.aic(cars.lo)) > $span > [1] 0.75 > > $aicc > [1] 6.93678 > > $aicc1 > [1] 167.7267 > > $gcv > [1] 5.275487 > > > > > > -- > Michael Friendly Email: friendly AT yorku DOT ca > Professor, Psychology Dept. > York University Voice: 416 736-5115 x66249 Fax: 416 736-5814 > 4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html > Toronto, ONT M3J 1P3 CANADA > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html >Thomas Lumley Assoc. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle
Dear Mike,
You could try
bestLoess <- function(model, criterion=c("aicc", "aicc1",
"gcv"),
spans=c(.05, .95)){
criterion <- match.arg(criterion)
f <- function(span) {
mod <- update(model, span=span)
loess.aic(mod)[[criterion]]
}
result <- optimize(f, spans)
list(span=result$minimum, criterion=result$objective)
}
A little experimentation suggests that aicc1 doesn't seem to behave
reasonably.
Regards,
John
--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
--------------------------------
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of
> Michael Friendly
> Sent: Thursday, November 17, 2005 9:58 AM
> To: R-help at stat.math.ethz.ch
> Subject: [R] loess: choose span to minimize AIC?
>
> Is there an R implementation of a scheme for automatic
> smoothing parameter selection with loess, e.g., by minimizing
> one of the AIC/GCV statistics discussed by Hurvich, Simonoff
> & Tsai (1998)?
>
> Below is a function that calculates the relevant values of AICC,
> AICC1 and GCV--- I think, because I to guess from the names
> of the components returned in a loess object.
>
> I guess I could use optimize(), or do a simple line search on
> span=, but I'm not sure how to use loess.aic to write a
> function that would act as a wrapper for loess() and return
> the mimimizing loess fit for a specified criterion.
>
> loess.aic <- function (x) {
> # extract values from loess object
> if (!(inherits(x,"loess"))) stop("Error: argument must
> be a loess object")
> span <- x$pars$span
> n <- x$n
> traceL <- x$trace.hat
> sigma2 <- sum( x$residuals^2 ) / (n-1)
> delta1 <- x$one.delta
> delta2 <- x$two.delta
> enp <- x$enp
>
> aicc <- log(sigma2) + 1 + 2* (2*(traceL+1)) / (n-traceL-2)
> aicc1<- n*log(sigma2) + n* (
> (delta1/(delta2*(n+enp)))/(delta1^2/delta2)-2 )
> gcv <- n*sigma2 / (n-traceL)^2
>
> result <- list(span=span, aicc=aicc, aicc1=aicc1, gcv=gcv)
> return(result)
> }
>
>
> > cars.lo <- loess(dist ~ speed, cars)
> >
> > (values <- loess.aic(cars.lo))
> $span
> [1] 0.75
>
> $aicc
> [1] 6.93678
>
> $aicc1
> [1] 167.7267
>
> $gcv
> [1] 5.275487
>
> >
>
>
> --
> Michael Friendly Email: friendly AT yorku DOT ca
> Professor, Psychology Dept.
> York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
> 4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html
> Toronto, ONT M3J 1P3 CANADA
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
The locfit package (which, I believe, contains an independent implementation of loess, plus more) contains the gcvplot() and aicplot() functions that I think can do this. Best, Andy> From: Michael Friendly > > Thanks very much, John > > The formula for AICC1 was transscribed from an ambiguously > rendered version (in the SAS documentation). This is a > corrected version. > > loess.aic <- function (x) { > if (!(inherits(x,"loess"))) stop("Error: argument must > be a loess object") > # extract values from loess object > span <- x$pars$span > n <- x$n > traceL <- x$trace.hat > sigma2 <- sum( x$residuals^2 ) / (n-1) > delta1 <- x$one.delta > delta2 <- x$two.delta > enp <- x$enp > > aicc <- log(sigma2) + 1 + 2* (2*(traceL+1)) / (n-traceL-2) > # aicc1<- n*log(sigma2) + n* ( > (delta1/(delta2*(n+enp)))/(delta1^2/delta2)-2 ) > aicc1<- n*log(sigma2) + n* ( > (delta1/delta2)*(n+enp)/(delta1^2/delta2)-2 ) > gcv <- n*sigma2 / (n-traceL)^2 > > result <- list(span=span, aicc=aicc, aicc1=aicc1, gcv=gcv) > return(result) > } > > > John Fox wrote: > > > Dear Mike, > > > > You could try > > > > bestLoess <- function(model, criterion=c("aicc", "aicc1", "gcv"), > > spans=c(.05, .95)){ > > criterion <- match.arg(criterion) > > f <- function(span) { > > mod <- update(model, span=span) > > loess.aic(mod)[[criterion]] > > } > > result <- optimize(f, spans) > > list(span=result$minimum, criterion=result$objective) > > } > > > > A little experimentation suggests that aicc1 doesn't seem to behave > > reasonably. > > > > Regards, > > John > > > > -------------------------------- > > John Fox > > Department of Sociology > > McMaster University > > Hamilton, Ontario > > Canada L8S 4M4 > > 905-525-9140x23604 > > http://socserv.mcmaster.ca/jfox > > -------------------------------- > > > > > >>-----Original Message----- > >>From: r-help-bounces at stat.math.ethz.ch > >>[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of > >>Michael Friendly > >>Sent: Thursday, November 17, 2005 9:58 AM > >>To: R-help at stat.math.ethz.ch > >>Subject: [R] loess: choose span to minimize AIC? > >> > >>Is there an R implementation of a scheme for automatic > >>smoothing parameter selection with loess, e.g., by minimizing > >>one of the AIC/GCV statistics discussed by Hurvich, Simonoff > >>& Tsai (1998)? > >> > >>Below is a function that calculates the relevant values of AICC, > >>AICC1 and GCV--- I think, because I to guess from the names > >>of the components returned in a loess object. > >> > >>I guess I could use optimize(), or do a simple line search on > >>span=, but I'm not sure how to use loess.aic to write a > >>function that would act as a wrapper for loess() and return > >>the mimimizing loess fit for a specified criterion. > >> > >>loess.aic <- function (x) { > >> # extract values from loess object > >> if (!(inherits(x,"loess"))) stop("Error: argument must > >>be a loess object") > >> span <- x$pars$span > >> n <- x$n > >> traceL <- x$trace.hat > >> sigma2 <- sum( x$residuals^2 ) / (n-1) > >> delta1 <- x$one.delta > >> delta2 <- x$two.delta > >> enp <- x$enp > >> > >> aicc <- log(sigma2) + 1 + 2* (2*(traceL+1)) / (n-traceL-2) > >> aicc1<- n*log(sigma2) + n* ( > >>(delta1/(delta2*(n+enp)))/(delta1^2/delta2)-2 ) > >> gcv <- n*sigma2 / (n-traceL)^2 > >> > >> result <- list(span=span, aicc=aicc, aicc1=aicc1, gcv=gcv) > >> return(result) > >>} > >> > >> > >> > cars.lo <- loess(dist ~ speed, cars) > >> > > >> > (values <- loess.aic(cars.lo)) > >>$span > >>[1] 0.75 > >> > >>$aicc > >>[1] 6.93678 > >> > >>$aicc1 > >>[1] 167.7267 > >> > >>$gcv > >>[1] 5.275487 > >> > >> > > >> > >> > >>-- > >>Michael Friendly Email: friendly AT yorku DOT ca > >>Professor, Psychology Dept. > >>York University Voice: 416 736-5115 x66249 Fax: 416 736-5814 > >>4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html > >>Toronto, ONT M3J 1P3 CANADA > >> > >>______________________________________________ > >>R-help at stat.math.ethz.ch mailing list > >>https://stat.ethz.ch/mailman/listinfo/r-help > >>PLEASE do read the posting guide! > >>http://www.R-project.org/posting-guide.html > > -- > Michael Friendly Email: friendly AT yorku DOT ca > Professor, Psychology Dept. > York University Voice: 416 736-5115 x66249 Fax: 416 736-5814 > 4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html > Toronto, ONT M3J 1P3 CANADA > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > >