[R] Confidence Bands in Polynomial Regression

[Rolf Turner]

Perhaps I'm just being obtuse, but I don't see what ols() has to do with
the question that was asked.

Often with R it is easier to roll your own rather than struggling  
with the
arcana of someone else's software.

Here is a function that seems to do what Ben wants:

pecb <- function(fit,tt,alpha=0.05,fill.col,...) {
#
# pecb <--> ``partial effect confidence band''.
#
ind <- c(2,5,6)
ccc <- coef(fit)[ind]
Sig <- summary(fit)$cov.unscaled[ind,ind]
M   <- cbind(1,tt,tt^2)
pe  <- M%*%ccc
sig <- sqrt(apply(M*t(Sig%*%t(M)),1,sum))
tv  <- -qt(alpha/2,fit$df.residual)
hi  <- pe + tv*sig
lo  <- pe - tv*sig
plot(tt,pe,type="n",ylim=range(hi,lo))
polygon(c(tt,rev(tt)),c(lo,rev(hi)),col=fill.col)
lines(tt,pe,...)
}

And here is a test script to try it out using simulated data:

# 1. Simulate some data.
set.seed(42)
tt <- seq(0,10,length=301)
x  <- sort(runif(301))
y  <- 0.25 + 1.3*x - 0.7*tt + 0.3*tt^2 + 0.4*x*tt - 0.6*x*tt^2 + rnorm 
(301)

# 2. Fit the model.
fit <- lm(y ~ x + tt + I(tt^2) + I(x*tt) + I(x*tt^2))

# 3. Plot the desired result.
pecb(fit,tt,fill.col="red")

# 4. Add in the ``true'' curve.
yt <- 1.3 + 0.4*tt - 0.6*tt^2
lines(tt,yt,lty=5,col="blue")

	cheers,

		Rolf Turner

On 16/06/2009, at 2:11 PM, Dylan Beaudette wrote:

> On Mon, Jun 15, 2009 at 6:57 PM, Ben Amsel<benamsel at gmail.com> wrote:
>>  Hello R users,
>>
>> Given a linear (in the parameters) regression model where one  
>> predictor x
>> interacts with time and time*time (ie, a quadratic effect of time t):
>> y = b0 + b1(x) + b2(t) + b3(t^2) + b4(x*t) + b5(x*t^2) + e,
>>
>> I would like to construct 95% confidence bands (optimally, shaded)  
>> around
>> this function:
>>
>> *dy* = b1 + b4(t) + b5(t^2)
>> *dx*
>>
>> That is, the partial effect of x on y changing over time t
>>
>> Is this possible with predict() or perhaps another function?
>>
>>
>> Thank you very much
>> Ben
>
> Hi Ben,
>
> Check out the 'Design' package, it has all kinds of convenience
> functions for plotting these type of things.
>
> install.packages('Design', dep=TRUE)
> ?ols

######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}