R help - Dec 2012 - bootstrap based confidence band

I'm trying to find a bootstrap based confidence band for a linear model.
I have created a data set with X and Y
X=runif(n,-1.25,1.25)
e=rnorm(n,0,1)
Y=exp(3*X)+5*sin((30*X)/(2*pi))+2*e
fit=lm(Y~X)
summary(fit)
 
I define a bootstrap function named PairedBootstrap which is not listed here.
Than I try many ways to find the confidence band. One way is to predict Y using
the model I get above for 1000 times.
Data.Orig=data.frame(cbind(X,Y))
B=1000
Boot.Result=matrix(nrow=B,ncol=n)
for ( b in 1:B){
 Data.Orig.Boot=PairedBootstrap(Data.Orig)
      fit.Boot=predict(fit,newdata=Data.Orig.Boot,type="response")
 Boot.Result[b,]=fit.Boot
 }
 
And try to find 95% confidence interval for 1000 copies of y corresponding to
each x.
ConfidenceSet.Pointwise=function(Boot.Result,alpha){
  n=ncol(Boot.Result)
  B=nrow(Boot.Result)
  SetBounds=matrix(ncol=2,nrow=n)
  for(j in 1:n){
            Result.Sort=sort(Boot.Result[,j])
  SetBounds[j,1]=Result.Sort[floor(0.5*B*alpha)]
       SetBounds[j,2]=Result.Sort[ceiling(B*(1-0.5*alpha))]
  }
  return(SetBounds)
  }

And then  try to line them up. But the result is not what I want.
alpha=0.05
 Boot.Pointwise=ConfidenceSet.Pointwise(Boot.Result,alpha)
 lines(Data.Orig[,1],Boot.Pointwise[,1], col=2, lwd=3)
 
Could someone tell where I'm wrong, or is there another better way to do it?
Thank you so much!

	[[alternative HTML version deleted]]

Hello,

Inline.
Em 06-12-2012 03:52, Rebecca escreveu:
> I'm trying to find a bootstrap based confidence band for a linear
model.
> I have created a data set with X and Y
> X=runif(n,-1.25,1.25)
> e=rnorm(n,0,1)
> Y=exp(3*X)+5*sin((30*X)/(2*pi))+2*e
> fit=lm(Y~X)
You cannot expect this model, a first degree polynomial in X,  to fit 
your data.
> summary(fit)
>   
> I define a bootstrap function named PairedBootstrap which is not listed
here. Than I try many ways to find the confidence band. One way is to predict Y
using the model I get above for 1000 times.
> Data.Orig=data.frame(cbind(X,Y))
Get rid of cbind(), please. It's doing nothing and the construct 
data.frame(cbind(...)) is potentially harmful.
> B=1000
> Boot.Result=matrix(nrow=B,ncol=n)
> for ( b in 1:B){
>   Data.Orig.Boot=PairedBootstrap(Data.Orig)
>       
fit.Boot=predict(fit,newdata=Data.Orig.Boot,type="response")
>   Boot.Result[b,]=fit.Boot
>   }
>   
> And try to find 95% confidence interval for 1000 copies of y corresponding
to each x.
> ConfidenceSet.Pointwise=function(Boot.Result,alpha){
>    n=ncol(Boot.Result)
>    B=nrow(Boot.Result)
>    SetBounds=matrix(ncol=2,nrow=n)
>    for(j in 1:n){
>              Result.Sort=sort(Boot.Result[,j])
>    SetBounds[j,1]=Result.Sort[floor(0.5*B*alpha)]
>         SetBounds[j,2]=Result.Sort[ceiling(B*(1-0.5*alpha))]
>    }
>    return(SetBounds)
>    }
Better could be


ConfidenceSet.Pointwise2 <- function(Boot.Result,alpha){
     SetBounds <- matrix(ncol=2, nrow=n)
     for(j in seq_len(ncol(Boot.Result)))
         SetBounds[j, ] <- quantile(Result.Sort, probs = c(alpha/2, 1 - 
alpha/2))
     SetBounds
}

Simpler, no?


Hope this helps,

Rui Barradas
>
> And then  try to line them up. But the result is not what I want.
> alpha=0.05
>   Boot.Pointwise=ConfidenceSet.Pointwise(Boot.Result,alpha)
>   lines(Data.Orig[,1],Boot.Pointwise[,1], col=2, lwd=3)
>   
> Could someone tell where I'm wrong, or is there another better way to
do it?
> Thank you so much!
>
> 	[[alternative HTML version deleted]]
>
>
>
> ______________________________________________
> 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.

	[[alternative HTML version deleted]]