Ejercicio Bootstrap

Ejemplo 1 (diapositiva 5 de la Clase 12)

Tengo datos de sobrevivencia de 16 ratones luego de una cirugía de prueba: hay 9 ratones en el grupo control y 7 ratones en el grupo de tratamiento.

Group	Survival time (in days)	Mean
Treatment	94,197,16,38,99,141,23	86.86
Control	52,104,146,10,51,30,40,27,46	56.22

¿Podemos decir que el tratamiento es efectivo?

En estadística, resolvemos esa pregunta estimando \(\bar{X}- \bar{Y} = 30.63\). El problema es cómo calcular la variabilidad, ¿podemos suponer lo mismo de siempre?

Datos

X<-c(94,197,16,38,99,141,23)
Y<-c(52,104,146,10,51,30,40,27,46)

Estimación de \(\bar{X}\).

1. Estimación por intervalo del tiempo promedio de sobrevivencia del tratamiento.

(estimacion<-mean(X))

## [1] 86.85714

sd(X)

## [1] 66.76683

n<-length(X)

Asumiendo normalidad de la población, tenemos que el intervalo de confianza de 95% es dado por:

estimacion+c(-1,1)*qt(p=0.975,df=(n-1))*sd(X)/sqrt(n)

## [1]  25.10812 148.60616

2. IC estándar de Bootstrap

B <- 1000
Tboot_b <- NULL
for(b in 1:B) {
  xb <- sample(X, size = n, replace = TRUE)
  Tboot_b[b] <- mean(xb)
}
Tboot_b[1:10]

##  [1]  97.28571 100.85714  91.85714  75.00000  91.00000  73.85714  90.42857
##  [8]  45.28571 123.57143 107.85714

hist(Tboot_b)

2.1. Sesgo del bootstrap

(sesgo_bootstrap <- mean(Tboot_b))

## [1] 86.322

estimacion

## [1] 86.85714

sesgo_bootstrap-estimacion

## [1] -0.5351429

2.2. Intervalos de confianza de 95%

(cuantil_z <- qnorm(1 - 0.05 / 2))

## [1] 1.959964

(sdboot <- sd(Tboot_b))

## [1] 23.6704

estimacion + c(-1,1)*cuantil_z * sdboot

## [1]  40.46402 133.25027

3. IC bootstrap t

B <- 1000
Tboot_b <- NULL
Tboot_bm <- NULL
sdboot_b <- NULL
for (b in 1:B) {
  xb <- sample(X, size = n, replace = TRUE)
  Tboot_b[b] <- mean(xb)
  for (m in 1:B) {
    xbm <- sample(xb, size = n, replace = TRUE)
    Tboot_bm[m] <- mean(xbm)
  }
  sdboot_b[b] <- sd(Tboot_bm)
}
z_star <- (Tboot_b - estimacion) / sdboot_b

hist(z_star)

summary(z_star)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -5.5573 -0.8638  0.0000 -0.1644  0.6457  7.3303

cuantil_t_empirico <- quantile(z_star, 1 - 0.05 / 2)
estimacion +c(-1,1) *cuantil_t_empirico* sdboot

## [1]  32.87144 140.84285

4. IC percentil de bootstrap

B <- 1000
Tboot_b <- NULL
for(b in 1:B) {
  xb <- sample(X, size = n, replace = TRUE)
  Tboot_b[b] <- mean(xb)
}
hist(Tboot_b)

quantile(Tboot_b,0.05 / 2)

##     2.5% 
## 44.14286

quantile(Tboot_b, 1 - 0.05 / 2)

##    97.5% 
## 135.4321

5. Con el paquete boots de R.

library(boot) 
mean_fun=function(datos,indice){ 
  m=mean(datos[indice]) 
  v=var(datos[indice]) 
  c(m,v)
}

X.boot<-boot(X,mean_fun,R=1000) 
(results<-boot.ci(X.boot,type=c("basic","stud","perc")))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = X.boot, type = c("basic", "stud", "perc"))
## 
## Intervals : 
## Level      Basic             Studentized          Percentile     
## 95%   ( 38.61, 130.57 )   ( 32.62, 167.18 )   ( 43.14, 135.11 )  
## Calculations and Intervals on Original Scale

Repitan el ejercicio pero para la diferencia de medias.