[6499] in s-news-athena
loess derivatives
daemon@ATHENA.MIT.EDU (clive@research.att.com)
Thu Feb 23 10:52:00 1995
From: clive@research.att.com
Date: Thu, 23 Feb 95 10:22 EST
To: s-news@utstat.toronto.edu
In response to Mary Lindstrom's questions on loess derivatives...
Following is a function lo.co() that reconstructs the kd tree from a loess
object, and also the fitted values and local slopes (first derivative
estimates). By default, loess estimates the local coefficients using
weighted least squares so is non-robust (I hope Trevor and I didn't claim
otherwise). Using family="symmetric" should be more robust to heavy tailed
residuals, although I doubt this has been investigated with regard to
derivative estimation.
Notes:
1) lo.co is slow and ugly, and may crash in cases with duplicate points
in the kd tree; e.g. if the predictors lie on a grid.
2) For interpolation/plotting of derivatives, one could interface
lo.co() with interp(). Since this is crude compared to the methods
used in predict.loess(), it may be necessary to reduce the cell
parameter in loess().
3) Derivative estimates are of course highly variable in sparse neighborhoods
and boundary regions. Beware!
> x1 <- rnorm(100)
> x2 <- rnorm(100)
> y <- 10*x1+5*x2+rnorm(100)
> fit <- loess(y~x1+x2)
> lo.co(fit)
x1 x2 fhat d.x1 d.x2
[1,] -1.76390007 -1.83632316 -29.477514 10.944953 7.2758035
[2,] 3.04106305 -1.83632316 21.271929 10.084532 5.1323782
...
[17,] 0.47300889 0.05036947 5.129923 10.192049 5.3252511
[18,] 0.47300889 2.66211397 17.122967 8.181994 3.9465924
-------------------------------------------------------------------------------
lo.co <- function(fit)
{
aln <- fit$surface$interpolator$all.numeric
d <- aln$parameter["d"]
coef <- matrix(aln$vval, ncol = d + 1, byrow = T) %*% diag(c(1, 1/fit$
surface$divisor))
x <- matrix(0, ncol = d, nrow = nrow(coef))
vc <- nv <- 2^d
for(i in 1:d)
x[1:vc, i] <- rep(rep(aln$vert[c(i, i + d)], rep(2^(i - 1), 2)),
2^(d - i))
cell <- matrix(1:vc, nrow = 1)
i <- nc <- 1
while(i <= nc) {
s <- aln$a[i]
if(s > 0) {
cell <- rbind(cell, cell[i, ], cell[i, ])
z <- x[cell[i, ], s]
u1 <- (1:vc)[z == z[1]]
u2 <- (1:vc)[z == z[vc]]
newv <- (nv + 1):(nv + vc/2)
cell[nc + 1, u2] <- cell[nc + 2, u1] <- newv
x[newv, ] <- x[cell[i, u1], ]
x[newv, s] <- aln$xi[i]
nc <- nc + 2
nv <- nv + vc/2
}
i <- i + 1
}
x <- x %*% diag(fit$surface$divisor)
ret <- cbind(x, coef)
dimnames(ret) <- list(NULL, c(fit$predictors$names.predictors, "fhat",
paste("d.", fit$predictors$names.predictors, sep = "")))
ret
}
