I am trying to (semi) calculate the confidence intervals for a loess smoother
(function: loess()), but have been thus far unsuccessful.
The CI for the loess predicted values, yhat, are apparently
yhat +- t*s * sqrt(w^2), where s is the residual sum of squares and w is the
weight function
Correct me of I'm wrong, but R uses the tricubic function (1-abs(z)^3)^3,
where z = (x-xi)/h, where h is the the number of neighboring values within
the bandwidth. Assuming that's correct, here is my code (for the first
observation in x1:
x1 <- 1:156 # x
set.seed(123)
y <- arima.sim(list(order=c(2,0,0), ar=c(1,-.1)), n=156) # randomly
generated y
lo1 <- loess(y ~ x1, span=0.4) # loess smoother for y ~ x
df <- lo1$one.delta # estimated df from r
wconstant <- summary(lo1)[17]$weights # set weight; otherwise 1 (as in
this case)
res <- residuals(lo1) # y-yhat
ss <- sum(wconstant *res^2)
s <- sqrt(ss/df) # r terms this residual standard error
x0 <- x1[1] # focal x for first observation
dist1 <- abs(x1 - x0) # distance from focal x
h <- sort(dist1)[.4*length(x1)] # bandwidth for span: bandwidth
span*length(x1)
inwindow <- dist1 <= h # observations within window
d <- dist1[dist1 <= h]
z <- d/h
w1 <- (1-abs(z)^3)^3
w2 <- sum(w1^2)
s*sqrt(w2) # calculated sefit
> s*sqrt(w2)
[1] 8.052887
this is quite high, and r gives me
> predict(lo1, se=T)$se.fit[1]
1
0.6055714
although I did calculate s correctly
> lo1$s
[1] 1.484307> s
[1] 1.484307
So clearly something is wrong with the weight part of my equation. I'm not
sure whether my CI equation doesn't match R's, my weight function is
wrong,
R does things quite differently than above, or simply my code is off. Any
help would be great.
--
View this message in context:
http://r.789695.n4.nabble.com/Loess-CI-tp4501215p4501215.html
Sent from the R help mailing list archive at Nabble.com.
