# Calculating probability for bivariate normal distributions based on bootstrapped regression coefficients

Dear all, I was encouraged to ask this question here as well as on stackoverflow and would be very appreciative of any answers...

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)
library(Design) # for ols()

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")  - add comment ## 1 Answer First I would like to point out that to get the object ks3 you do not need to do that much of the book-keeping code. Use the features of ddply: ks3 <- ddply(o,.(ai,Gs),function(d){ temp <- ols(value~as.numeric(bc)+as.numeric(age),data=d,x=T,y=T) t2 <- bootcov(temp,B=1000,coef.reps=T) data.frame(t2$boot.Coef,ai=d$ai[1],gc2=d$Gs[1])
})


Also scales=free works much better if you want to see what kind of graph you intend to provide.

Concerning your first question, if you estimate the standard errors of coefficient using bootstrap, you are absolutely correct in using estimated percentiles. There is no need to use the percentiles of bivariate normal, since then you assume normality and bootstrap is used exactly for the purpose of not attaching oneself to some theoretical distribution.

For your second question, first you need to define what do you mean by bias-corrected percentiles. Furthermore your question implies that bootcov returns non bias-corrected percentiles, but it does not return any percentiles at all.

Update

If we know that coefficients are normal we can use confidence ellipses. The mean and covariance matrix are estimated using bootstrap. The following code checks which of the observations of the bivariate sample falls into 0.95 confidence ellipse region:

 confeps <- function(x, level=0.95, t=qchisq(level,2)) {
cv <- cov(x)
ch <- t(chol(x))
x <- sweep(x, 2, apply(x, 2, mean), "-")
y <- solve(ch,t(x))
tt <- y[1, ]^2 + y[2, ]^2
tt < t
}


Here is the example how this function works:

      bs <- cbind(a <- rnorm(1000), 2-a/3+rnorm(1000)/5)
mn <- apply(bs, 2, mean)
require(ellipse)
eps <- ellipse(cov(bs), centre=mn, npoints=1000)
plot(eps, type="l", xlim=range(bs[,1]), ylim=range(bs[, 2]))
col <- rep(1, nrow(x))
col[!confeps(bs)] <- 2
points(bs[, 1], bs[, 2], col=col)

-
Thx. But both bootstrap estimates have a normal distribution and I'd really like to give a joint probaility for the two occuring at the same time. –  Misha Dec 22 '10 at 13:18
@Misha, why do bootstrap estimates should have normal distribution? Usually bootstrap is used when the distribution is unknown. Or you know that the distribution is normal and you use bootstrap to estimate the covariance matrix? –  mpiktas Dec 23 '10 at 10:38
I was uncertain of the distribution, but as it happens both A & B turn out normally distributed. However what I really would like to know is to estimate the probability of both A & B occuring at the same time. –  Misha Dec 23 '10 at 11:09