BOOTSTRAP

Si tratta di una tecnica introdotta a fine anni ’70 con cui calcolare il valore atteso, varianza o altre proprietà (PARAMETRO THETA) di uno stimatore o di una procedura statistica, quando la risoluzione analitica non è agevole.

L’idea generale è quella di fare inferenza con simulazioni in ambito frequentista (RICAMPIONAMENTO dall’unico campione x a disposizione; campionamento casuale semplice, i.e. estrazioni con reinserimento).

Sui B-campioni ottenuti si calcola la stima di theta (se theta=varianza, si calcola la varianza campionaria ecc.) e si prende la media sui B valori come stima finale (Bootstrap Estimate).

x <- c(28,-44, 29, 30, 24, 28, 37, 32, 36, 27, 26, 28, 29, 26, 27, 22, 23, 20, 25, 25, 36, 23, 31, 32, 24, 27, 33, 16, 24, 29, 36, 21, 28, 26, 27, 27, 32, 25, 28, 24, 40, 21, 31, 32, 28, 26, 30, 27, 26, 24, 32, 29, 34, -2, 25, 19, 36, 29, 30, 22, 28, 33, 39, 25, 16, 23
)

bootstrap <- function(data, B, theta, ...) {
    z <- list()
  # genero x1,...,xB campioni a partire dal campione x
  datas <- matrix(sample(data, size=length(data)*B, replace=T),nrow=B)
  # calcolo theta(x1), ..., theta(xB)
  B.estimate <- apply(datas, 1, theta, ...)
  # calcolo theta(x) stima classica
  z$est <- theta(data, ...)
  z$distn <- B.estimate
  # stima del Bias (distorsione delle stime bootstrap)
  z$bias <- mean(B.estimate) - z$est
  z$sd <- sd(B.estimate)
  z
}

# vogliamo stimare la tendenza centrale
# considera 2 stimatori: media e media troncata con alpha=.05

mean.b <- bootstrap(data=x, B=1000, theta=mean)
df1 <- data.frame(media=mean.b$est, bias=mean.b$bias, sd=mean.b$sd)
df1
##      media        bias       sd
## 1 26.21212 -0.03842424 1.339121
trigger.mean.b <- bootstrap(data=x, B=1000, theta=mean, trim=.05)
df2 <- data.frame(media.trigger=trigger.mean.b$est, bias=trigger.mean.b$bias, sd=trigger.mean.b$sd)
df2
##   media.trigger        bias        sd
## 1          27.4 -0.08301667 0.8175357
hist(mean.b$distn, freq=F, xlab='bootstrap estimates', main='')
abline(v=mean.b$est, col='red', lty=2)

hist(trigger.mean.b$distn, freq=F, xlab='bootstrap estimates', main='')
abline(v=trigger.mean.b$est, col='red', lty=2)

Tema d’esame: generare un campione di 20 valori casuali da una Gamma(4,10) con seme dell’estrazione 127. Calcola per il parametro media (=alpha/beta) una stima bootstrap e costruisci un intervallo di confidenza bootstrap base di livello 95%.

set.seed(127)
n = 20
x <- rgamma(n, 4, 10)
m <- 4/10
# genero B=200 campioni a partire da x
B = 200
datas <- matrix(data=sample(x, size=length(x)*B, replace=T), nrow=B)
# stime
m.est <- mean(x)                  # classica
m.distn <- apply(datas, 1, mean)  # B stime 
m.bias <- mean(m.distn) - m.est   # bias
m.sd <- sd(m.distn)               # sqrt(varianza campionaria B)

# se la distorsione/(standard error) è < 0.25 non è necessario usare metodi per ridurre la distorsione
m.bias/m.sd <= 0.25
## [1] TRUE
# costruisco IC di livello 95% per theta=media di una Gamma(4,10)
# 1 - 2*alpha = 0.95 -> alpha = 0.05/2 = 0.025
# estraggo da Gamma(4,10) i quantili di livello 0.975,0.025
quantiles <- c(qgamma(p=0.975, 4, 10), qgamma(p=0.025, 4, 10))
IC <- c(m.est - m.sd*quantiles[1], m.est - m.sd*quantiles[2])

# l'IC include la vera media 
df <- data.frame(media.vera=4/10, media.stimata=m.est, lower.bound=IC[1], upper.bound=IC[2]); df
##   media.vera media.stimata lower.bound upper.bound
## 1        0.4     0.4373786   0.3951094   0.4321241

Il Bootstrap può essere utilizzato anche in ambito di regressione e modelli lineari generalizzati. Si usa il pacchetto boot con l’omonima funzione, che prevede 3 argomenti: il primo è il dataframe su cui lavoriamo, il secondo è il modello lineare (generalizzato) che vogliamo implementare con indici per le unità ricampionate, il terzo è il numero di ricampionamenti B. Per fornire il secondo argomento dobbiamo quindi fare uno step intermedio e creare una funzione bootstrap customizzata, funzione dei dati e degli indici.

Esempio 1 : Un certo dispositivo medico atto a rilasciare un certo ormone antinfiammatorio nel corpo in modo continuo nel tempo è stato testato su 27 soggetti. Ai soggetti è stato casualmente assegnato un dispositivo proveniente da uno di tre lotti (Lot).E’ stato poi misurato l’ammontare di ormone rimanente (amount) nel dispositivo dopo averlo indossato per un determinato periodo di tempo (hours).

library(bootstrap)
## 
## Attaching package: 'bootstrap'
## The following object is masked _by_ '.GlobalEnv':
## 
##     bootstrap
library(boot)
library(MASS)
data(hormone)
set.seed(123); B = 1000

# inferenza classica 
hormone.lm <- lm(amount ~ 0 + Lot + hrs, data=hormone)
coef(hormone.lm)
##        LotA        LotB        LotC         hrs 
## 32.13159495 36.10509461 35.59732433 -0.06013605
# bootstrap sulle unità statistiche (x,y) , x=c(Lot, hrs)
# bootstrap statistic : function of data 
boot.lm <- function(data, index) {
  data <- data[index,]
  mod <- lm(amount ~ 0 + Lot + hrs, data=data)
  coef(mod)
}

bootstrap.hormone <- boot(data=hormone, statistic=boot.lm, R=B)
bootstrap.hormone
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = hormone, statistic = boot.lm, R = B)
## 
## 
## Bootstrap Statistics :
##        original        bias    std. error
## t1* 32.13159495  0.0435835807  0.76888509
## t2* 36.10509461  0.1834510547  1.21962760
## t3* 35.59732433  0.0604831462  0.62038195
## t4* -0.06013605 -0.0004758288  0.00395957
# per ridurre il bias posso applicare la correzione
# theta.BIAS.CORRECTED = 2theta - bootstrap.estimate

BC <- 2*coef(hormone.lm)[2] - (1/B)*sum(bootstrap.hormone$t[1001:2000])

new.bias <- coef(hormone.lm)[2] - BC
new.bias
##      LotB 
## 0.1834511
#plot(bootstrap.hormone, index=1)
#plot(bootstrap.hormone, index=2)
#plot(bootstrap.hormone, index=3)
#plot(bootstrap.hormone, index=4)

#head(boot.array(bootstrap.hormone), n=5)

# bootstrap sui residui (bias minore se modello correttamente specificato; meno robusto della variante sulle unità statistiche)

fit <- fitted(hormone.lm)

e <- residuals(hormone.lm)

X <- model.matrix(hormone.lm)
boot.lm.res <- function(data, index) {
  y = fit + e[index]
  mod = lm(y ~ X - 1)
  coef(mod)
}

boot(data=hormone, statistic=boot.lm.res, R=B)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = hormone, statistic = boot.lm.res, R = B)
## 
## 
## Bootstrap Statistics :
##        original        bias    std. error
## t1* 32.13159495 -0.0076871101 0.701040512
## t2* 36.10509461  0.0069272759 0.909694907
## t3* 35.59732433 -0.0072456100 0.564658786
## t4* -0.06013605  0.0000594167 0.003212912

Esempio 2 : E’ stato condotto un esperimento in cui un diverso numero di scarafaggi (m) è stato esposto a diverse dosi del veleno (ldose, in scala logaritmica), registrando poi il numero di scarafaggi morti a ciascun dosaggio (r). Il numero di decessi è funzione delle dosi del veleno e si distribuisce come una binomiale (family binomial) con size m e probabilità p=G(a+b*ldose), dove G è una funzione di ripartizione continua.
1.Normale probit
2.Logistica logit
3.a valori estremi complementary log-log

# obiettivo: trovare il valore di dose (in scala logaritmica)
#             t.c. il 50% della popolazione di scarafaggi muore.
X <- matrix(data=c(59, 6, 1.6907, 60, 13, 1.7242, 62, 18, 1.7552, 56, 28, 1.7842, 63, 52, 1.8113, 59, 53, 1.8369, 62, 61, 1.8610, 60, 60, 1.8839), ncol=3, byrow=T)

X <- data.frame(m=X[,1], r=X[,2], ldose=X[,3])

# inferenza classica
# G(a + b*x50) = 0.5; G probit (cdf N)
# x50 = (G^-1(0.5) - aGLS)/bGLS
mod <- glm(cbind(r, m-r) ~ ldose, data=X, family=binomial(link='probit'))
xp <- dose.p(mod, p=0.5) #p=0.5 default value
xp
##              Dose          SE
## p = 0.5: 1.770852 0.003803312
# bootstrap
boot.glm <- function(data, index) {
data <- data[index,]
mod <- glm(cbind(r, m-r) ~ ldose, data=data,family=binomial(link='probit'))
dose.p(mod)
}

boot.estimate <- boot(data=X, statistic=boot.glm, R=1000)
boot.estimate
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = X, statistic = boot.glm, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original       bias    std. error
## t1* 1.770852 0.0003508581 0.005340713
# il bias è 4*10^-4 (praticamente nullo) 
# lo SE è superiore rispetto alla stima classica

JACKKNIFE

Dato un campione x d’ampiezza n, definiamo l’i-esimo campione JackKnife come x[-i] d’ampiezza n-1 (rimuovo semplicemente l’i-esima unità di x).
Su ciascun campione JackKnife viene calcolato lo stimatore d’interesse e, al solito, si prende la media degli n valori (non più B, user-specified) così calcolati.

Le stime della distorsione e della varianza dello stimatore vengono corrette per un fattore moltiplicativo (n-1) che tiene conto del fatto che i campioni JackKnife sono molto simili tra loro.

# esempio: stima MEDIA e DEV.STD dei tempi di sopravvivenza dopo un intervento

x <- c(94,197,16,38,99,141,23)
media <- mean(x)
dev.std <- sd(x)/sqrt(length(x))

# jacknife 
require(bootstrap)

theta <- function(x) {mean(x)}
results <- jackknife(x, theta)
c(media, mean(results$jack.values))
## [1] 86.85714 86.85714
c(dev.std, results$jack.se)
## [1] 25.23549 25.23549
results$jack.bias
## [1] 0

Per concludere il JackKnife si può vedere come un’approssimazione del BootStrap, computazionalmente molto meno onerosa (se n è inferiore a B), però meno efficiente (soprattutto in caso di stimatori non lineari nei dati).