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<benam...@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}}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
