[6354] in s-news-athena
Summary: confidence intervals of derivatives
daemon@ATHENA.MIT.EDU (Bill Shipley)
Wed Feb 1 16:16:32 1995
Date: Wed, 1 Feb 1995 15:54:08 -0500
To: S-news@utstat.toronto.edu
From: bshipley@courrier.usherb.ca (Bill Shipley)
Estimating confidence intervals of the 1st derivatives obtained from
smoother splines.
Although one can obtain confidence intervals of the function values
(predicted y values) of smoother splines via gam(), the only way to obtain
estimates of the 1st derivatives of the function values is via
smooth.spline() followed by predict.smooth.spline() and these S-PLUS
functions do not output confidence intervals. There are two ways of
obtaining these confidence intervals.
The first way way provided by Trevor Hastie. The key reference is: Hastie,
T.J. & R.J. Tibshirani. 1990. Generalized Additive Models. Chapman and Hall,
especially chapter 3 and p. 127).
Hastie writes: "One simple approach is to acknowledge that the derivative
operator is linear (for a given value of the smoothing parameter) in the
response y, hence it can be written as d=Gy for some matrix G. Given G, and
appropriate iid error assumptions, one can compute the covariance matrix:
Cov(d)=GG^T*sigma. If you prefer the Bayesian style posterior intervals,
these will also involve G and perhaps some function of the smoothing matrix
S itself.
How to get G? and S?
Well, G=GI where I is the identity matrix (also S=SI), so apply
smooth.spline() and predict.smooth.spline(,deriv=1) repeatedly to the
columns of I, and then paste the results together."
Here is how to get G and S (you must specify df or span to be the same as
used to fit the data):
G<-matrix(0,n,n)
I<-diag(n)
for(i in 1:n){
fit<-smooth.spline(x, I[,i], df=value)
S[,i]<-fit$y
G[,i]<-predict.smooth.spline(fit,x,deriv=1)$y
}
The standard error of the fitted values (sigma.2 is the residual variance) is:
sqrt(diag(S%*%t(S))*sigma.2)
The standard error of the 1st derivative is:
sqrt(diag(G%*%t(G))*sigma.2)
Example: given x,y and using cross-validation to choose the smoother value:
y.fit<-smooth.spline(x,y,cv=T)
deriv.fit<-predict.smooth.spline(y.fit,y.fit$x,deriv=1)
approx.95.confidence interval<-1.95*sqrt(diag(G%*%t(G))*y.fit$cv.crit)
upper.95<-deriv.fit+approx.95.confidence.interval
lower.95<-deriv.fit-approx.95.confidence.interval
Another method of estimating these confidence intervals is to use bootstrap
techniques. Repeatedly sample, with replacement, from the data (using
sample()), and save the calculated 1st derivatives. Do this many times and
extract the desired quantiles.
I have tryed this, and it gives results very similar to Hastie's method
above. Unfortunately, S-PLUS is terrible for doing bootstrap studies
because of its speed and memory problems.
Note two things:
1) These derivatives are systematically biased at the extremes of the data.
This is because the second derivatives are 0 beyond the extremes and the
1st derivatives are therefore constant. Thus, if the 1st derivatives are,
in reality, rapidly changing beyond the extremes, the estimates will be
biased. Even on very "wiggly" true values, I find that the bias disappears
quickly.
2) the cross-validation proceedure may occasionally over- or undersmooth the
data. These proceedures are only asymtotically correct.
If anyone has any ideas on how to correct for the bias (above) please let me
know.
Bill Shipley
Departement de Biologie
Universite de Sherbrooke
Sherbrooke (Quebec)
CANADA J1K 2R1
bshipley@courrier.USherb.ca
