Cálculo de la probabilidad de distribuciones normales bivariadas basadas en coeficientes de regresión bootstrapped

Question

Hola a todos, Me anime a hacer esta pregunta aquí, así como en stackoverflow y estaría muy agradecido de cualquier respuestas...

Debido a hetereoscedasticity estoy haciendo bootstrap de regresión lineal (apela más a mí que la regresión robusta). Me gustaría crear una trama a lo largo de las líneas de lo que he hecho en la secuencia de comandos aquí. Sin embargo el fill=int no es correcto desde int debe (creo) se calcula utilizando una distribución normal bivariante.

• Alguna idea de cómo podría hacer que en esta configuración?
• También hay una manera para bootcov de devolución de corrección de sesgo percentiles?

ejemplo de secuencia de comandos:

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") 

Answer 1

En primer lugar me gustaría señalar que para obtener el objeto ks3 usted no necesita hacer mucho de el libro de mantenimiento de código. Uso de las funciones de 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])
})

También scales=free funciona mucho mejor si quieres ver qué tipo de gráfico que se desea proporcionar.

Respecto a su primera pregunta, si la estimación de los errores estándar de coeficiente de uso de bootstrap, usted está absolutamente en lo correcto en el uso estimado de los percentiles. No hay necesidad de utilizar los percentiles de bivariante normal, ya que, a continuación, usted asume la normalidad y bootstrap se usa exactamente para el propósito de no sujetar a sí mismo para algunos teóricos de la distribución.

Para tu segunda pregunta, primero es necesario definir ¿qué quieres decir con corrección de sesgo de percentiles. Además, su pregunta implica que bootcov devuelve la no corrección de sesgo percentiles, pero no se devuelve ningún percentiles.

Actualización

Si sabemos que los coeficientes son normales podemos utilizar la confianza de los puntos suspensivos. La media y la matriz de covarianza se calcula mediante bootstrap. El siguiente código comprueba que de las observaciones de la bivariante de la muestra cae en confianza 0.95 elipse región:

 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
    }

Aquí está el ejemplo de cómo esta función funciona:

      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)