Create p-values from bootstrapped multivariate normal distribution

Question

Dear all,

Due to hetereoscedasticity I'm doing bootstrapped linear regression (appeals more to me than robust regression). I'd like to create a plot along the lines of what I've done in the script here. However the fill=int is not right since int should (I believe) be calculated using a bivariate normal distribution.

-Any idea how I could do that in this setting?

-Also is there a way for bootcov to return bias-corrected percentiles?

sample script:

library(ggplot2)
library(Hmisc)

o<-data.frame(value=rnorm(10,20,5),bc=rnorm(1000,60,50),age=rnorm(1000,50,20),ai=as.factor(round(runif(1000,0,4),0)),Gs=as.factor(round(runif(1000,0,6),0)))



reg.s<-function(x){

    ols(value~as.numeric(bc)+as.numeric(age),data=x,x=T,y=T)->temp
    bootcov(temp,B=1000,coef.reps=T)->t2

    return(t2)
    }

dlply(o,.(ai,Gs),function(x) reg.s(x))->b.list
llply(b.list,function(x) x[["boot.Coef"]])->b2

ks<-llply(names(b2),function(x){
    s<-data.frame(b2[[x]])
    s$ai<-x
    return(s)
    })


ks3<-do.call(rbind,ks)
ks3$ai2<-with(ks3,substring(ai,1,1))

ks3$gc2<-sapply(strsplit(as.character(ks3$ai), "\\."), "[[", 2)


k<-ks3
j<-dlply(k,.(ai2,gc2),function(x){
    i1<-quantile(x$Intercept,probs=c(0.025,0.975))[1]
i2<-quantile(x$Intercept,probs=c(0.025,0.975))[2]

j1<-quantile(x$bc,probs=c(0.025,0.975))[1]
j2<-quantile(x$bc,probs=c(0.025,0.975))[2]

o<-x$Intercept>i1 & x$Intercept<i2

p<-x$bc>j1 & x$bc<j2

h<-o & p
    return(h)

    })
m<-melt(j)
ks3$int<-m[,1]  

ggplot(ks3,aes(x=bc,y=Intercept,fill=int))+geom_point(,alpha=0.3,size=1,shape=21)+facet_grid(gc2~ai2,scales = "free_y")+theme_bw()->plott
plott<-plott+opts(panel.grid.minor=theme_blank(),panel.grid.major=theme_blank())
plott<-plott+geom_vline(x=0,color="red")
plott+xlab("BC coefficient")+ylab("Intercept")