[BUG] Bootstrap avec au moins une strate avec une seule obs.

[Dominique Soudant]

La question ne porte pas sur la pertinence de l'opération mais de son éventuelle occurrence dans le flux d'un grand nombre de bootstrap et des conséquences de cette opération.

R 2.4.1
boot 1.2-27

Soit donc 8 observations se répartissant sur 8 mois et un bootstrap stratifié par mois, mettons 10 répétitions pour fixer les idées. On s'attend à avoir 10 fois la même chose, n'est-ce pas ? Ben non, il y a une variabilité de résultat.

Dans mon cas, j'ai identifié le problème dans la fonction ordinary.array() qui créé l'index de rééchantillonnage:

Code : Tout sélectionner

ordinary.array <- function (n, R, strata)
{
    output <- matrix(0, R, n)
    inds <- as.integer(names(table(strata)))
    for (is in inds) {
        gp <- c(1:n)[strata == is]
        output[, gp] <- matrix(sample(gp, R * length(gp), replace = TRUE),
            nrow = R)
    }
    output
}

Voici l'index de rééchantillonnage produit lorsque l'on a plus d'une observation par strate :

Code : Tout sélectionner

Browse[1]> ordinary.array(8,10,rep(1:4,each=2))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
 [1,]    1    1    4    4    5    6    7    8
 [2,]    1    1    4    4    5    6    7    8
 [3,]    1    1    4    4    5    5    8    7
 [4,]    2    1    4    4    5    5    7    8
 [5,]    1    2    3    4    6    6    8    8
 [6,]    2    1    3    4    6    5    8    7
 [7,]    2    1    3    4    6    5    7    7
 [8,]    1    2    4    4    5    6    8    8
 [9,]    1    2    3    4    6    6    7    8
[10,]    1    2    4    4    6    6    7    7

Nickel, les deux premières colonnes ne concernent que les deux premiers éléments et ainsi de suite, c'est tout bon.

Voici un exemple de résultat avec une observation par strate :

Code : Tout sélectionner

Browse[1]> ordinary.array(8,10,1:8)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
 [1,]    1    1    2    2    1    1    5    2
 [2,]    1    1    2    2    3    3    2    5
 [3,]    1    2    1    4    2    1    1    7
 [4,]    1    1    3    2    5    5    3    6
 [5,]    1    1    2    4    5    3    7    1
 [6,]    1    2    2    4    4    6    1    8
 [7,]    1    1    2    2    4    2    7    5
 [8,]    1    2    3    2    3    4    1    7
 [9,]    1    2    2    1    1    3    1    2
[10,]    1    1    2    4    5    1    1    5

La première colonne est ok, mais la seconde contient des 1 des 2 alors qu'elle ne devrait contenir que des 2 etc. jusqu'a la 8 qui contient des éléments entre 1 et 8.

Le problème vient de sample() qui :

If x has length 1 and x >= 1, sampling takes place from 1:x

lorsqu'on lui fournit un scalaire x en entrée il effectue un échantillonnage sur (1, 2, 3, ... x), et c'est bien ce que l'on observe dans mon exemple ci-dessus.



Bon ...




J'raconte n'importe quoi ou c'est un bug ? Parce que, pour en revenir à mon titre :

Code : Tout sélectionner

> ordinary.array(7,10,c(rep(1:3,each=2),4))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
 [1,]    1    2    3    3    6    5    6
 [2,]    2    2    3    4    5    5    7
 [3,]    2    2    4    4    6    5    5
 [4,]    1    1    4    4    6    5    1
 [5,]    2    2    3    3    6    5    7
 [6,]    1    2    3    4    5    5    1
 [7,]    1    2    4    3    6    6    5
 [8,]    2    1    3    4    5    6    2
 [9,]    2    2    3    4    6    6    4
[10,]    2    2    4    3    6    5    2


dans ce cas j'ai bien une seule strate avec une observation, mais l'index de cette strate avec observation unique va de 1 à 7 ... :(

[Dominique Soudant]

Avec

R 2.8.1

et

Package: boot
Version: 1.2-35
Date: 2009-01-01

ca fait :

Code : Tout sélectionner

> ordinary.array <- function (n, R, strata)
+ {
+     inds <- as.integer(names(table(strata)))
+     if (length(inds) == 1L) {
+         output <- sample(n, n * R, replace = TRUE)
+         dim(output) <- c(R, n)
+     }
+     else {
+         output <- matrix(as.integer(0L), R, n)
+         for (is in inds) {
+             gp <- (1L:n)[strata == is]
+             output[, gp] <- sample(gp, R * length(gp), replace = TRUE)
+         }
+     }
+     output
+ }
> ordinary.array(7,10,c(rep(1:3,each=2),4))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
 [1,]    1    2    4    3    5    5    4
 [2,]    1    2    3    3    6    6    1
 [3,]    2    1    3    4    5    6    3
 [4,]    2    1    4    3    5    5    6
 [5,]    1    2    3    4    5    5    7
 [6,]    2    2    4    3    5    6    5
 [7,]    2    1    3    3    6    6    6
 [8,]    1    2    4    4    5    5    2
 [9,]    2    1    3    3    6    6    5
[10,]    1    2    4    3    5    6    4

Le bug est toujours là.

[Dominique Soudant]

Bon ben j'vais signaler un bug alors.

EDIT :
Hop !
http://bugs.r-project.org/cgi-bin/R/con ... ctid=13584

[Renaud Lancelot]

Et c'est maintenant corrigé !

Code : Tout sélectionner

> boot:::ordinary.array(7, 10, c(rep(1:3, each = 2), 4))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
 [1,]    1    1    3    3    5    6    7
 [2,]    2    2    3    3    6    5    7
 [3,]    2    1    4    4    5    5    7
 [4,]    2    1    3    3    6    5    7
 [5,]    2    1    4    4    6    5    7
 [6,]    1    2    4    3    6    6    7
 [7,]    1    2    3    4    6    6    7
 [8,]    1    2    3    3    6    6    7
 [9,]    1    1    4    3    6    6    7
[10,]    2    2    3    3    5    6    7


Nice shot !

Renaud

[Dominique Soudant]

Merci pour l'info.

Si tu ne me l'avais pas signalé, je n'aurais pas cessé de checker frénétiquement ma notification de bug.