Intervalos de confianza bootstrap

CREACIÓN DE UN INTERVALO DE CONFIANZA BOOSTRAP DE UNOS DATOS.

dades <- read.table('Law School.txt',header=T)

Se sabe que I.C para colas iguales

Siendo ^r=r(x) entonces el I.C:

[] donde es b=b1, ... ,B

recordemos que y que =

¿Cómo estimaremos el valor de error estándar?

Una opción es tomando que los datos son normales entonces

se = y la estimación boostrap será: se* =

Creación de las funciones de la correlación muestral, el error estándar y el estadístico t.

CORRELACIÓN

correl <-function(xy)

{

return(cor(xy[,1],xy[,2]))

}

correl(dades)

[1] 0.7599979

ERROR ESTÁNDAR DE LA CORRELACIÓN MUESTRAL

se.corr <-function(rho.val,rho.factor)

{

return((1- rho.val * rho.val )* rho.factor)

}

#rho.factor es el cálculo de 1/raiz(n-3)

se.corr(correl(dades),82)

[1] 34.63707

ESTADÍSTICO T

t.corr <- function(rho.val,rho.val.boot, se.val)

{

return((rho.val.boot- rho.val)/se.val)

}

A partir de aquí generamos las remuestras y llamamos a la función. En este caso utilizaremos la función for para que no se complique la cosa, pero considerando que debemos tener calculado el error estandar, la correlación sobre la muestra original y el tamaño del valor poblacional.

xy <- read.table('Law School.txt',header=T)

ic.boot<- function(xy,alpha=0.05,B=999)

{

n<-nrow(xy)

rho.factor<-1/sqrt(n-3)

rho.val<-correl(xy)

se.val<-se.corr(rho.val, rho.factor)

t.boot<-numeric(B)

for (b in 1 :B)

{

xy.boot<- xy[sample(1:n,replace=T),]

rho.val.boot <-correl(xy.boot)

se.val.boot <-se.corr(rho.val.boot, rho.factor)

t.boot[b] <- t.corr (rho.val,rho.val.boot, se.val.boot)

}

t.boot<-sort(t.boot)

t.boot.inf<- t.boot [ round((B+1)*(1-alpha/2))]

t.boot.sup<- t.boot [ round((B+1)*(alpha/2))-1]

# como deseamos que el interior de nuestro intervalo contenga el 95% de los datos, hacemos unos interalos conservadores con la aproximación de B+1 para los dos límites, y además para el superior le quitaremos una posición asegurando casi el 0,025% de los datos.

par(c(2,2))

plot(t.boot)

hist(t.boot)

return(list(rho=rho.val ,se= se.val,t.inf= t.boot.inf ,t.sup= t.boot.sup ,ic.inf=rho.val- t.boot.inf *se.val, ic.sup=rho.val- t.boot.sup *se.val))

}

xy.ic <-ic.boot(xy)

#comprobamos que se realiza de forma correcta el cálculo del intervalo de confianza inferior, será idéntico para el superior.

xy.ic$ic.inf

[1] 0.6325478

xy.ic$rho- xy.ic$t.inf* xy.ic$se

[1] 0.6325478

plot(xy.ic$t)

xy.ic

$rho

[1] 0.7599979

$se

[1] 0.04752408

$t.inf

[1] 2.681800

$t.sup

[1] -1.750585

$ic.inf

[1] 0.6325478

$ic.sup

[1] 0.8431928

hist(xy.ic$t)

Generamos una distribución bivariante normal fijando los parámetros de la r

Sabemos por teoría que f=1/2*log((1+ r)/(1-r)) que estimandolo ^f=1/2*log((1+ ^r)/(1-^r))

Obtienes el intervalo ^f+-z(a/2)/raiz(n-3) => [^f1,^f2]

Que quitando la transformación de normalización de la correlación obtenemos [^r1,^r2],

Introducimos la función de la binomial normalizada

rbivnorm <-

function(n, mean1=0, mean2=0, sd1=1, sd2=1, rho=0){

if (n < 0) stop("Numero de replicas de simulación negativo")

if (sd1 < 0 || sd2 < 0) stop("La desviacion estándar no puede ser negativa")

if (rho < -1.0 || rho > 1.0) stop("La correlación tiene que estar entre -1 i +1")

bet <- rho * sd2 / sd1

alf <- mean2 - bet * mean1

sdy <- sd2 * sqrt(1 - rho*rho)

x <- rnorm(n, mean1, sd1)

return (matrix( c(x,rnorm(n, bet*x + alf, sdy)), ncol=2))

}

xy<-rbivnorm(n=2000,rho=0.67)

cor<-correl(xy)

cor

[1] 0.6806052

DEFINIMOS LA INVERSA DE F

fi.invers <-function(fi)

{

e<-exp(2*fi)

return((e-1)/(e+1))

}

CALCULAMOS LA ESTIMACIÓN DE INTERVALOS DE CONFIANZA DE LA CORRELACIÓN HACIENDO LA TRANSFORMACIÓN A NORMAL Y DESHACIÉNDOLA, YA QUE ESTAMOS EN EL CASO DE DISTRIBUCIONES NORMALES BIVARIADAS.

fi<-function(xy,alpha=0.05){

cor<-correl(xy)

fi.val<-0.5*log((1+cor)/(1-cor))

n<-nrow(xy)

rho.factor<-1/sqrt(n-3)

fi.sup <-fi.val+ qnorm(alpha/2)*rho.factor

fi.inf <- fi.val- qnorm(alpha/2)*rho.factor

return(fi.sup, fi.inf,cor.sup=fi.invers(fi.sup), cor.inf=fi.invers(fi.inf))

}

#El intervalo real de confianza es:

fi(xy)

$fi.sup

[1] 0.7863815

$fi.inf

[1] 0.8740996

$cor.sup

[1] 0.6563543

$cor.inf

[1] 0.703451

#El intervalo de confianza Boostrap es:

xy.ic<-ic.boot(xy)

xy.ic

$rho

[1] 0.6806052

$se

[1] 0.01201170

$t.inf

[1] 2.136675

$t.sup

[1] -1.956239

$ic.inf

[1] 0.65494

$ic.sup

[1] 0.7041029

El código R

> source("T:/Alumnes/Román/Clase_20304/bootstrap-t.r")

> source("T:/Alumnes/Román/Clase_20304/rbivnorm.r")

> law.school<-read.table("Law School.txt",header=T)

> n <- nrow(law.school) # n es el tamaño muestral

> indices <- sample(1:n, replace=T) # "indices" define una remuestra

> indices

[1] 75 25 19 42 60 43 15 10 72 6 40 19 51 46 3 16 54 77 47 26 73 56 32 5 57

[26] 77 56 37 54 44 55 67 52 63 68 65 35 4 77 53 26 26 31 66 58 23 7 26 38 59

[51] 68 10 6 32 55 79 73 61 49 53 2 32 9 41 1 5 50 12 20 59 36 14 76 62 52

[76] 81 53 1 79 40 56 35

> correl(law.school, indices) # correlacion de la remuestra definida por "indices", calculada de otra forma

[1] 0.7021442

> correl(law.school, 1:n)

[1] 0.7599979

#hay otra manera alternativa para el bootstrap no parametrico

#es emplear vectores no paramétricos

> presamp <-tabulate(indices)

> indices

[1] 75 25 19 42 60 43 15 10 72 6 40 19 51 46 3 16 54 77 47 26 73 56 32 5 57

[26] 77 56 37 54 44 55 67 52 63 68 65 35 4 77 53 26 26 31 66 58 23 7 26 38 59

[51] 68 10 6 32 55 79 73 61 49 53 2 32 9 41 1 5 50 12 20 59 36 14 76 62 52

[76] 81 53 1 79 40 56 35

> presamp

[1] 2 1 1 1 2 2 1 0 1 2 0 1 0 1 1 1 0 0 2 1 0 0 1 0 1 4 0 0 0 0 1 3 0 0 2 1 1 1

[39] 0 2 1 1 1 1 0 1 1 0 1 1 1 2 3 2 2 3 1 1 2 1 1 1 1 0 1 1 1 2 0 0 0 1 2 0 1 1

[77] 3 0 2 0 1

# Para ordenar los indices de menor a mayor

> sort(indices)

[1] 1 1 2 3 4 5 5 6 6 7 9 10 10 12 14 15 16 19 19 20 23 25 26 26 26

[26] 26 31 32 32 32 35 35 36 37 38 40 40 41 42 43 44 46 47 49 50 51 52 52 53 53

[51] 53 54 54 55 55 56 56 56 57 58 59 59 60 61 62 63 65 66 67 68 68 72 73 73 75

[76] 76 77 77 77 79 79 81

> presamp <- presamp / n

> presamp

[1] 0.02439024 0.01219512 0.01219512 0.01219512 0.02439024 0.02439024

[7] 0.01219512 0.00000000 0.01219512 0.02439024 0.00000000 0.01219512

[13] 0.00000000 0.01219512 0.01219512 0.01219512 0.00000000 0.00000000

[19] 0.02439024 0.01219512 0.00000000 0.00000000 0.01219512 0.00000000

[25] 0.01219512 0.04878049 0.00000000 0.00000000 0.00000000 0.00000000

[31] 0.01219512 0.03658537 0.00000000 0.00000000 0.02439024 0.01219512

[37] 0.01219512 0.01219512 0.00000000 0.02439024 0.01219512 0.01219512

[43] 0.01219512 0.01219512 0.00000000 0.01219512 0.01219512 0.00000000

[49] 0.01219512 0.01219512 0.01219512 0.02439024 0.03658537 0.02439024

[55] 0.02439024 0.03658537 0.01219512 0.01219512 0.02439024 0.01219512

[61] 0.01219512 0.01219512 0.01219512 0.00000000 0.01219512 0.01219512

[67] 0.01219512 0.02439024 0.00000000 0.00000000 0.00000000 0.01219512

[73] 0.02439024 0.00000000 0.01219512 0.01219512 0.03658537 0.00000000

[79] 0.02439024 0.00000000 0.01219512

#en un contexto de remuestreo no paramétrico tiene que ser lo mismo

#la función correl hace el mismo cálculo

> correl(law.school, p=presamp)

[1] 0.7356295

Warning messages:

1: longer object length

is not a multiple of shorter object length in: x * p

2: longer object length

is not a multiple of shorter object length in: y * p

3: longer object length

is not a multiple of shorter object length in: y * p

4: longer object length

is not a multiple of shorter object length in: x * x * p

> correl

function (datos, indices, p, ...) {

if (missing (indices)){

if (missing (p))

return (cor (datos[,1], datos[,2]))

else {

x <- datos[,1]

y <- datos[,2]

x <- x - sum (x * p)

y <- y - sum (y * p)

yp <- y*p

return ((sum(x*yp))/sqrt ((sum(x*x*p))*(sum(y*yp))))

}

else

return (cor (datos[indices,])[1,2])

}

# los puntos suspensivos quieren decir que significa que (ver después del # # return) llamaria a la función y lo pasaría donde los puntos suspensivos

# Definicion de se.correl, calculo del error estándar de la correlacion #muestral

> se.correl

function(datos, rho, se.fact, ...){

if (missing (rho))

rho <- correl(datos, ...)

if (missing (se.fact))

se.fact <- 1 / sqrt (nrow(datos) - 3)

return ((1 - rho * rho) * se.fact)

}

> se.r <- se.correl(law.school) # A partir de los datos originales

> se.correl(law.school, rho=r) # Sin necesidad de calcular de nuevo la correlacion

Error in se.correl(law.school, rho = r) : Object "r" not found

> se.const <- 1 / sqrt(n - 3)

> se.correl(rho=r, se.fact=se.const)

Error in se.correl(rho = r, se.fact = se.const) :

Object "r" not found

#hay una posibilidad de la que hablaremos de aquí a unos dias

#que es una manera reciente de hallar errores estandard y que

#es mucho mas cómodo plantearla empleando vectores de remuestreo.

#puede cortocircuitar, por razones de eficiencia, algunas cosas, si le damos

#rho ya no lo calcula y si le damos el denominador de la expresión, ya no

#lo calculose.correl

> se.correl

function(datos, rho, se.fact, ...){

if (missing (rho))

rho <- correl(datos, ...)

if (missing (se.fact))

se.fact <- 1 / sqrt (nrow(datos) - 3)

return ((1 - rho * rho) * se.fact)

}

#Finalmente tenemos otra función que calcula la correlación studentizada

> stud.correl

function(datos, rho.samp, rho, se, ...){

if (missing (rho.samp))

rho.samp <- correl(datos, ...)

if (missing (se))

se <- se.correl(datos, rho = rho.samp, ...)

return ((rho.samp - rho) / se)

}

> stud.boot.correl

function(datos, B = 999){

if (B < 0) stop("Numero negativo o nulo de replicas bootstrap")

tboot <- numeric(length=B)

n <- nrow(datos)

r <- correl(datos)

se.const <- 1 / sqrt (n - 3)

se.r <- se.correl(rho=r, se.fact=se.const)

indexos <- 1:n

for (b in 1:B) {

tboot[b] <- stud.correl(datos[sample(indexos, replace=T),],rho=r, se.fact=se.const)

}

attr(tboot, "rho") <- r

attr(tboot, "se.rho") <- se.r

return (tboot)

}

#la función attr permite definir atributos de cualquier objeto y de manera #dinámica crea o permite calcular atributos

> tboot <- stud.boot.correl(law.school)

attr(,"rho")

[1] 0.7599979

attr(,"se.rho")

[1] 0.04752408

#este tboot lo emplearemos como un vector, si decimos siete nos

#referimos a la posición siete del vector. Para recuperarlo llamamos

> attr(tboot, "rho")

[1] 0.7599979

> r <- attr(tboot, "rho")

> icboot.t.correl

function(datos, alfa=0.05, simetrizado=F, B=999){

remuestreo <- boot.t.correl( datos, B)

if (simetrizado) {

tU <- quantile (abs (remuestreo), probs = 1 - alfa)

tL <- -tU

}

else {

alfa2 <- 0.5*alfa

t.boot <- quantile (remuestreo, probs = c(alfa2,1-alfa2))

tL <- t.boot[1]

tU <- t.boot[2]

}

rho <- attr(remuestreo,"rho")

se.rho <- attr(remuestreo,"se.rho")

ic <- c(rho - tU*se.rho, rho - tL*se.rho)

names(ic) <- NULL

return(ic)

}

#finalmente esta función nos devuelve el intervalo de confianza bootstrap-t

> icboot.t.correl(law.school)

[1] 0.6272623 0.8409062

> icboot.t.correl(law.school, alfa=0.01)

[1] 0.5902067 0.8628937

> icboot.t.correl(law.school,simetrizado=T)

[1] 0.6572470 0.8627487

> icboot.t.correl(law.school,sim=T)

[1] 0.6463448 0.8736509

#están hechas para la correlación pero a poco que se utilicen pueden

#emplearse para cualquier estadístico, con ciertas limitaciones, etc.

#Para el caso concreto de datos normal bivariante y el coeficiente de correlación de Pearson tenemos el #cambio de Fisher que si tengo rho puedo cambiar de escala mediante la siguiente transformación: