[R] Loess CI

[sw1]

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.