Introdução

Seja uma amostra \(y = (y_1, y_2, \ldots, y_n)\) e quantidade de interesse (parâmetro) \(\theta\) para o qual se tem um estimador \[\hat{\theta} = \hat{\theta}(y).\]

O objetivo do método jackknife é estimar o viés e o erro padrão de \(\hat{\theta}\).

As amostras jackknife \(y(i)\) de tamanho \((n-1)\) são obtidas retirando-se cada um dos dados por vez \[y(i) = (y_1, y_2, \ldots, y_{i-1}, y_{i+1}, \ldots, y_n)\] Há portanto \(n\) amostras jackknife e para cada uma delas calcula-se uma estimativa \[\hat{\theta}_{(i)} = \hat{\theta}(y(i)).\] Simplificando notação denotamos \(\overline{\hat{\theta}_{(i)}} = \bar{\theta}\)

Viés

Estimativa de viés é obtida por \[ viés_J(\hat{\theta}) = (n-1)(\bar{\theta} - \hat{\theta}). \]

A expressão do estimador corrigido pelo viés é: \[ \begin{align} \hat{\theta}_J &= \hat{\theta} - viés_J(\hat{\theta})\\ &= n \hat{\theta} - (n-1)\bar{\theta}. \end{align} \]

Erro padrão

\[\text{se}(\hat{\theta}_J) = \sqrt{\frac{n-1}{n} \sum_{i=1}^n(\hat{\theta}_{(i)} - \bar{\theta})^2}\]

Pseudo-valores

Os pseudo-valores (\(pv_i\)) são valores calculados para cada amostra jackknife. \[pv_{i} = n \hat{\theta} - (n-1)\hat{\theta}_{(i)}.\]

Expressões do viés e erro padrão do estimador corrigido podem ser diretamente obtidas como função dos pseudo-valores. \[ \begin{align} \hat{\theta}_J &= \overline{pv_i} = \frac{1}{n} \sum_{i=1}^n pv_i\\ \text{se}(\hat{\theta}_J) &= \frac{S_{pv}}{\sqrt{n}} = \sqrt{\frac{1}{n}\frac{1}{n-1} \sum_{i=1}^n (pv_i - \overline{pv_i})^2}\\ \end{align} \]

Exemplo 1: Índice de Gini de diversidade

O artigo de (Dixon 2001) fornece dados com os quais ilustra os procedimentos de Jackknife e Bootstrap.

A estatística de interesse é o índice de Gini, \[G = \frac{\sum_{i=1}^n (2i-n-1) y_{(i)}}{(n-1) \sum_{i=1}^n y_i}, \] em que os valores de \(y\) são ordenados e \(y_{(i)}\) denota o valor ordenado na \(i\)-ésima posição.

Seguem código com o estimador na forma de uma função, os dados e a estimativa para dados fornecidos.

gini <- function(x){
    n <- length(x)
    xs <- sort(x)
    Num <- sum((2*(1:n) - n-1)*xs)
    Den <- (n-1)*sum(xs)
    return(Num/Den)
}

plantsI <- c(18, 18, 20, 23, 25, 28)
n <- length(plantsI)
(plantsI.G <- gini(plantsI))    ## estimativa
## [1] 0.1121212

Estimativas individuais para as amostras jackknife e a média destas,

(plantsI.Gi <- sapply(1:n, function(i) gini(plantsI[-i])))
## [1] 0.10964912 0.10964912 0.12053571 0.12385321 0.11682243 0.09134615
(plantsI.G.bar <- mean(plantsI.Gi))  
## [1] 0.111976

permitem estimar o viés (estimativa de) e estimativa corrigida.

(vies <- (n-1) * (mean(plantsI.Gi) - plantsI.G))
## [1] -0.0007262651
plantsI.G - vies
## [1] 0.1128475

O erro padrão étambém obtido a partir das estimativas das amostras jackknife.

sqrt(((n-1)/n) * sum((plantsI.Gi - mean(plantsI.Gi))^2))
## [1] 0.0237038
sqrt(((n-1)^2/n) * var(plantsI.Gi))
## [1] 0.0237038

Os pseudo valores permitem obter os mesmos resultados de forma direta.

(plantsI.pv <- n*plantsI.G - (n-1)*plantsI.Gi)
## [1] 0.12448166 0.12448166 0.07004870 0.05346122 0.08861512 0.21599650
##
(est.jack <- mean(plantsI.pv))
## [1] 0.1128475
(vies <- plantsI.G - est.jack)
## [1] -0.0007262651
##
(se.jack <- sd(plantsI.pv)/sqrt(n))
## [1] 0.0237038

Escrevendo função em R

O procedimento pode ser implementado em uma função para ser usado em outros dados e estimadores.

jack <- function(x, estimador, ...){
    ## x é vetor de dados
    ## estimador é uma função que recebe dados e calcula quantidade de interesse
    th.hat <- estimador(x, ...)
    n <- length(x)
    th.hat.i <- sapply(1:n, function(i) estimador(x[-i], ...))
    th.bar <- mean(th.hat.i)
    pv <- n*th.hat - (n-1)*th.hat.i 
    th.jack <- mean(pv)
    th.se <- sd(pv)/sqrt(n)
    vies <- (n-1) * (mean(th.hat.i) - th.hat)
    return(list(est = th.hat, est.jack = th.jack, mean.jack = th.bar, vies = vies,
                se = th.se, est.i = th.hat.i, pseudovals = pv))
}    
jack(plantsI, gini)
## $est
## [1] 0.1121212
## 
## $est.jack
## [1] 0.1128475
## 
## $mean.jack
## [1] 0.111976
## 
## $vies
## [1] -0.0007262651
## 
## $se
## [1] 0.0237038
## 
## $est.i
## [1] 0.10964912 0.10964912 0.12053571 0.12385321 0.11682243 0.09134615
## 
## $pseudovals
## [1] 0.12448166 0.12448166 0.07004870 0.05346122 0.08861512 0.21599650

Há diversas implementações em pacotes do R. Entre elas aunção do pacote bootstrap.

bootstrap::jackknife(x = plantsI, theta = gini)
## $jack.se
## [1] 0.0237038
## 
## $jack.bias
## [1] -0.0007262651
## 
## $jack.values
## [1] 0.10964912 0.10964912 0.12053571 0.12385321 0.11682243 0.09134615
## 
## $call
## bootstrap::jackknife(x = plantsI, theta = gini)

jackknifeR::jackknife(statistic = gini, d = 1, data = data.frame(x = plantsI))
resample::jackknife(data = plantsI, statistic = gini)

Exemplo 2

Dados de (Manly 2006). Estimar \(\theta = \sigma\) com estimador \[\hat{\sigma} = \sqrt{\frac{\sum_{i=1}^n (y_i - \bar{y})^2}{n}}\]

dat <- c(3.56, 0.69, 0.10, 1.84, 3.93, 1.25, 0.18, 1.13, 0.27, 0.50,
         0.67, 0.01, 0.61, 0.82, 1.70, 0.39, 0.11, 1.20, 1.21, 0.72)
mean(dat)
## [1] 1.0445
var(dat)
## [1] 1.122921
sd(dat)
## [1] 1.05968
n <- length(dat)
(th.hat <- sqrt(((n-1)/n)*var(dat)))
## [1] 1.032848

sigma <- function(x){
    n <- length(x)
    sigmaMLE <- sqrt(((n-1)/n)*var(x))
    return(sigmaMLE)
}
sigma(dat)
## [1] 1.032848

n <- length(dat)
th.i <- numeric(n)
for(i in 1:n){
    yjack <- dat[-i]
    th.i[i] <- sigma(yjack)
}
th.i
##  [1] 0.8788364 1.0563893 1.0360975 1.0430060 0.8134128 1.0585751 1.0399595
##  [8] 1.0594885 1.0438813 1.0519008 1.0560070 1.0313246 1.0547329 1.0583612
## [15] 1.0483872 1.0484218 1.0365998 1.0590473 1.0589633 1.0569234
sapply(1:n, function(i) sigma(dat[-i]))
##  [1] 0.8788364 1.0563893 1.0360975 1.0430060 0.8134128 1.0585751 1.0399595
##  [8] 1.0594885 1.0438813 1.0519008 1.0560070 1.0313246 1.0547329 1.0583612
## [15] 1.0483872 1.0484218 1.0365998 1.0590473 1.0589633 1.0569234

Viés, estimativa corrigida e erro padrão.

(th.bar <- mean(th.i))
## [1] 1.029516
(vies <- (n-1)*(th.bar - th.hat))
## [1] -0.06330983
(th.J <- th.hat - vies)
## [1] 1.096158
n * th.hat - (n-1)*th.bar
## [1] 1.096158
c(th.hat, th.bar, th.J, sd(dat))
## [1] 1.032848 1.029516 1.096158 1.059680
(th.se <- sqrt(((n-1)/n)* sum((th.i - th.bar)^2)))
## [1] 0.2728037

Pseudo valores e cálculo alternativo de viés e erro padrão.

(th.pv <- n*th.hat - (n-1)*th.i)
##  [1] 3.9590656 0.5855601 0.9711049 0.8398438 5.2021138 0.5440315 0.8977269
##  [8] 0.5266770 0.8232137 0.6708425 0.5928255 1.0617899 0.6170325 0.5480939
## [15] 0.7376002 0.7369426 0.9615623 0.5350591 0.5366545 0.5754139
mean(th.pv)
## [1] 1.096158
sd(th.pv)/sqrt(n)
## [1] 0.2728037

Repetir para \(\bar{y}\), \(s\), \(s^2\).

Exemplo 3

Jackknife no logarítmo da variância \[\theta = \log(\sigma^2)\]

dados <- c(17.23, 18.71, 13.93, 18.81, 15.78, 11.29, 14.91, 13.39,
           18.21, 11.57, 14.28, 10.94, 18.83, 15.52, 13.45, 15.25)

jFun <- function(x) log(var(x))
bootstrap::jackknife(x = dados, theta = jFun)
## $jack.se
## [1] 0.2611095
## 
## $jack.bias
## [1] -0.03379347
## 
## $jack.values
##  [1] 1.994428 1.903270 2.024676 1.894992 2.034906 1.880898 2.038602 2.008562
##  [9] 1.940382 1.904688 2.031877 1.847697 1.893301 2.037588 2.010659 2.038948
## 
## $call
## bootstrap::jackknife(x = dados, theta = jFun)
jack(dados, jFun)
## $est
## [1] 1.970095
## 
## $est.jack
## [1] 2.003889
## 
## $mean.jack
## [1] 1.967842
## 
## $vies
## [1] -0.03379347
## 
## $se
## [1] 0.2611095
## 
## $est.i
##  [1] 1.994428 1.903270 2.024676 1.894992 2.034906 1.880898 2.038602 2.008562
##  [9] 1.940382 1.904688 2.031877 1.847697 1.893301 2.037588 2.010659 2.038948
## 
## $pseudovals
##  [1] 1.6051027 2.9724729 1.1513842 3.0966369 0.9979341 3.3080517 0.9424848
##  [8] 1.3930922 2.4157877 2.9511994 1.0433719 3.8060600 3.1220065 0.9576974
## [15] 1.3616347 0.9372996

Exemplo 4

Um exemplo interessante é o caso com \(Y \sim U[0, \theta]\).
Para as amostras jackknife:

o maior dado será a estimativa (n-1) vezes
o segundo maior dado será a estimativa 1 vez (quando excluir o máximo)

theta = 5
n = 10
set.seed(30)
(y <- round(runif(n, min = 0, max = theta), dig = 3))
##  [1] 0.494 2.441 1.820 2.103 1.505 0.738 4.493 1.118 4.830 0.705
sort(y)
##  [1] 0.494 0.705 0.738 1.118 1.505 1.820 2.103 2.441 4.493 4.830

(th.hat <- max(y))
## [1] 4.83
(th.bar <- ((n-1)*max(y) + 1 * sort(y)[n-1])/n)
## [1] 4.7963
(th.jack <- n * th.hat - (n-1)*th.bar)        ## estimativa corrigido para viés
## [1] 5.1333

Usando funções.

bootstrap::jackknife(x = y, theta = max)
## $jack.se
## [1] 0.3033
## 
## $jack.bias
## [1] -0.3033
## 
## $jack.values
##  [1] 4.830 4.830 4.830 4.830 4.830 4.830 4.830 4.830 4.493 4.830
## 
## $call
## bootstrap::jackknife(x = y, theta = max)
jack(y, max)
## $est
## [1] 4.83
## 
## $est.jack
## [1] 5.1333
## 
## $mean.jack
## [1] 4.7963
## 
## $vies
## [1] -0.3033
## 
## $se
## [1] 0.3033
## 
## $est.i
##  [1] 4.830 4.830 4.830 4.830 4.830 4.830 4.830 4.830 4.493 4.830
## 
## $pseudovals
##  [1] 4.830 4.830 4.830 4.830 4.830 4.830 4.830 4.830 7.863 4.830

Exemplo 5

Um caso de avaliação bioequivalência de adesivos médicos apresentado em (Efron and Tibshirani 1993) é frequentemente citado como exemplo para método de jackknife. O objetivo do estudo é determinar se adesivos fabricados em um novo local são (bio-)equivalentes a anteriores. O estudo utilizou oito pacientes que tiveram níveis de certo hormônio mendido após utilizar três adesivos (patches) diferentes:
1. sem hormônio (placebo),
2. adesivo fabricado no local original,
3. adesivo fabricado no novo local.

Os dados são:

individuo	placebo	original	novo
1	9243	17649	16449
2	9671	12013	14614
3	11792	19979	17274
4	13357	21816	23798
5	9055	13850	12560
6	6290	9806	10157
7	12412	17208	16570
8	18806	29044	26325

O critério para considerar que há bioequivalência é \[ \frac{|\mathbb{E}[\text{novo}] - \mathbb{E}[\text{original}]|}{\mathbb{E}[\text{original}] - \mathbb{E}[\text{placebo}]} \leq 0,20 \]

Neste caso o estimador de interesse é uma razão \[ \theta = \frac{\mathbb{E}[\text{novo}] - \mathbb{E}[\text{original}]}{\mathbb{E}[\text{original}] - \mathbb{E}[\text{placebo}]} \] e não há expressão fechada para o erro padrão.

Denotanto as variáveis medidas por \(Y_p\) para placebo, \(Y_o\) para original e \(Y_n\) para novo, definimos \[ Y_1 = Y_o - Y_p \;\mbox{ e }\; Y_2 = Y_n - Y_o\] O estimador de interesse é definido por \[\hat{\theta} = \frac{\overline{Y}_2}{\overline{Y}_1}\] que para os dados em questão produz \[\hat{\theta} = \frac{-452.2}{6342.4} = -0.071.\]

Daqui, pode-se seguir para obter amostras jackknife para estimar o vício, obter o estimador corrigido e o respectivo erro padrão. A função jack() vista anteriormente foi definida para dados em um único vetor e precisa ser adaptada para receber a estrutura de dados e estimador deste problema.

thBioeq <- function(dados){mean(dados[,"y2"])/mean(dados[,"y1"]) }
thBioeq(df)
## [1] -0.0713061
jackBioeq <- function(dados, estatistica, ...){
    ## dados é data.frame (ou matriz) incluindo colunas nomeadas y1 e y2
    ## estatística é função que recebe dados e calcula quantidade de interesse (estimador)
    th.hat <- estatistica(dados)
    n <- nrow(dados)
    th.hat.i <- sapply(1:n, function(i) estatistica(dados[-i,], ...))
    th.bar <- mean(th.hat.i)
    bias <- (n-1)*(th.bar - th.hat)
    th.jack <- th.hat - bias
    th.pv <- n*th.hat - (n-1)*th.hat.i 
    ##bias <- mean(th.pv) - th.hat
    th.se <- sqrt(sum((th.pv-mean(th.pv))^2)/(n*(n-1)))
    return(list(est = th.hat, est.jack = th.bar, est.cor = th.jack, se = th.se,
                bias = bias, est.i = th.hat.i, pseudovalores = th.pv))
}    
jackBioeq(df, thBioeq)
## $est
## [1] -0.0713061
## 
## $est.jack
## [1] -0.07016288
## 
## $est.cor
## [1] -0.07930858
## 
## $se
## [1] 0.1055278
## 
## $bias
## [1] 0.008002488
## 
## $est.i
## [1] -0.05711856 -0.12849970 -0.02145610 -0.13245033 -0.05067038 -0.08404803
## [7] -0.06486298 -0.02219698
## 
## $pseudovalores
## [1] -0.17061885  0.32904914 -0.42025606  0.35670355 -0.21575610  0.01788742
## [7] -0.11640789 -0.41506989

Alternativamente pode-se utilizar alguma função “pronta” de algum pacote como bootstrab::jackknife(), o que também requer organizar dados e função do estimador de maneira adequada.

thBio <- function(x, data){mean(data[x,"y2"])/mean(data[x,"y1"]) }
thBio(1:8, df)
## [1] -0.0713061
bootstrap::jackknife(1:8, thBio, data = df)
## $jack.se
## [1] 0.1055278
## 
## $jack.bias
## [1] 0.008002488
## 
## $jack.values
## [1] -0.05711856 -0.12849970 -0.02145610 -0.13245033 -0.05067038 -0.08404803
## [7] -0.06486298 -0.02219698
## 
## $call
## bootstrap::jackknife(x = 1:8, theta = thBio, data = df)

Comentários

Método pode falhar se a estatística/estimador não é suave, no sentido que que pequenas variações nos dados devem provocar apenas pequenas alterações nas estimativas.
Por exemlo, a mediana não é suave. De forma geral quantis amostrais não são estatísticas suaves.
No exemplo com dados de (Manly 2006) \(\hat{\theta}_{(i)}\) só tem dois possíves valores, 069 e 0,72.
Uma alteração tentando contornar este problema é o delete-d jackknife qe remova \(d\) amostras em cada passo. Neste caso tem-se \(\binom{n}{d}\) amostras de tamanho \(n-d\). Para remover problemas de não suavidade recomenda-se \(\sqrt{n} < d < n\) e a amostra deve ser suficientemente grande.
O método bootstrap não tem problemas de suavidade e é, em geral, mais atrativo.

Jackknife