[R] loess: choose span to minimize AIC?

[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