Je produis un script pour créer des échantillons de bootstrap à partir de l' catsensemble de données (à partir du -MASS-package).

En suivant le manuel de Davidson et Hinkley [1], j'ai effectué une régression linéaire simple et adopté une procédure fondamentale non paramétrique pour le bootstrap à partir des observations iid, à savoir le rééchantillonnage des paires .

L'échantillon d'origine se présente sous la forme:

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

Grâce à un modèle linéaire univarié, nous voulons expliquer le poids du foyer du chat à travers le poids de son cerveau.

Le code est:

library(MASS)
library(boot)


##################
#   CATS MODEL   #
##################

cats.lm <- glm(Hwt ~ Bwt, data=cats)
cats.diag <- glm.diag.plots(cats.lm, ret=T)


#######################
#   CASE resampling   #
#######################

cats.fit <- function(data) coef(glm(data$Hwt ~ data$Bwt)) 
statistic.coef <- function(data, i) cats.fit(data[i,]) 

bootl <- boot(data=cats, statistic=statistic.coef, R=999)

Supposons maintenant qu'il existe une variable de regroupement cluster = 1, 2,..., 24(par exemple, chaque chat appartient à une portée donnée). Pour simplifier, supposons que les données soient équilibrées: nous avons 6 observations pour chaque cluster. Ainsi, chacune des 24 portées est composée de 6 chats (ie n_cluster = 6et n = 144).

Il est possible de créer une fausse clustervariable via:

q <- rep(1:24, times=6)
cluster <- sample(q)
c.data <- cbind(cats, cluster)

J'ai deux questions connexes:

Comment simuler des échantillons conformément à la structure (groupée) de l'ensemble de données? Autrement dit, comment rééchantillonner au niveau du cluster? Je voudrais échantillonner les grappes avec remplacement et définir les observations dans chaque grappe sélectionnée comme dans le jeu de données d'origine (c'est-à-dire échantillonner avec remplacement des grappes et sans remplacer les observations dans chaque grappe).

C'est la stratégie proposée par Davidson (p. 100). Supposons que nous prélevions des B = 100échantillons. Chacun d'eux devrait être composé de 24 grappes éventuellement récurrentes (par exemple cluster = 3, 3, 1, 4, 12, 11, 12, 5, 6, 8, 17, 19, 10, 9, 7, 7, 16, 18, 24, 23, 11, 15, 20, 1), et chaque grappe devrait contenir les mêmes 6 observations de l'ensemble de données d'origine. Comment faire ça en R? (avec ou sans le -boot-paquet.) Avez-vous d'autres suggestions pour continuer?

La deuxième question concerne le modèle de régression initial. Supposons que j'adopte un modèle à effets fixes , avec des interceptions au niveau du cluster. Modifie-t-il la procédure de rééchantillonnage adoptée?

[1] Davidson, AC, Hinkley, DV (1997). Méthodes d'amorçage et leurs applications . La presse de l'Universite de Cambridge.

Le rééchantillonnage de l'ensemble des grappes est connu dans les statistiques d'enquête depuis aussi longtemps que des méthodes de rééchantillonnage y ont été utilisées (c'est-à-dire depuis le milieu des années 1960), c'est donc une méthode bien établie. Voir ma collection de liens à http://www.citeulike.org/user/ctacmo/tag/survey_resampling . Que bootce soit possible ou non, je ne sais pas; J'utilise le surveypackage lorsque j'ai besoin de travailler avec des bootstraps d'enquête, bien que la dernière fois que j'ai vérifié, il n'avait pas toutes les fonctionnalités dont j'avais besoin (comme quelques petites corrections d'échantillon, pour autant que je m'en souvienne).

Je ne pense pas que l'application d'un modèle particulier tel que les effets fixes change beaucoup les choses, mais, à mon avis, le bootstrap résiduel fait beaucoup d'hypothèses solides (les résidus sont iid, le modèle est correctement spécifié). Chacun d'entre eux est facilement cassé, et la structure du cluster rompt sûrement l'hypothèse iid.

Il existe de la littérature économétrique sur le bootstrap de cluster sauvage. Ils ont prétendu avoir travaillé dans le vide sans toutes ces cinquante années de recherches statistiques sur le sujet, donc je ne sais pas trop quoi en penser.

J'ai essayé de résoudre le problème moi-même et j'ai produit le code suivant.

Bien que cela fonctionne, il pourrait probablement être amélioré en termes de vitesse. Aussi, si possible, j'aurais préféré trouver un moyen d'utiliser le -boot-package, car il permet de calculer automatiquement un certain nombre d'intervalles de confiance bootstrap à travers boot.ci...

Par souci de simplicité, l'ensemble de données de départ se compose de 18 chats (les observations de «niveau inférieur») imbriqués dans 6 laboratoires (la variable de regroupement). L'ensemble de données est équilibré ( n_cluster = 3pour chaque cluster). Nous avons un régresseur x, pour expliquer y.

Le faux ensemble de données et la matrice où stocker les résultats sont:

  # fake sample 
  dat <- expand.grid(cat=factor(1:3), lab=factor(1:6))
  dat <- cbind(dat, x=runif(18), y=runif(18, 2, 5))

  # empty matrix for storing coefficients estimates and standard errors of x
  B <- 50 # number of bootstrap samples
  b.sample <- matrix(nrow=B, ncol=3, dimnames=list(c(), c("sim", "b_x", "se_x")))
  b.sample[,1] <- rep(1:B)

À chacune des Bitérations, la boucle suivante échantillonne 6 grappes avec remplacement, chacune composée de 3 chats échantillonnés sans remplacement (c'est-à-dire que la composition interne des grappes est maintenue inchangée). Les estimations du coefficient du régresseur et de son erreur standard sont stockées dans la matrice précédemment créée:

  ####################################
  #   loop through "b.sample" rows   #
  ####################################

  for (i in seq(1:B)) {

  ###   sampling with replacement from the clustering variable   

    # sampling with replacement from "cluster" 
    cls <- sample(unique(dat$lab), replace=TRUE)
    cls.col <- data.frame(lab=cls)

    # reconstructing the overall simulated sample
    cls.resample <- merge(cls.col, dat, by="lab")


  ###   fitting linear model to simulated data    

    # model fit
    mod.fit <- function(data) glm(data$y ~ data$x)

    # estimated coefficients and standard errors
    b_x <- summary(mod.fit(data=cls.resample))$coefficients[2,1]
    	se_x <- summary(mod.fit(data=cls.resample))$coefficients[2,2]

    b.sample[i,2] <- b_x
    b.sample[i,3] <- se_x

  }

J'espère que cela vous aide, Lando

Voici une manière beaucoup plus simple (et presque sans doute plus rapide) de faire le bootstrap en utilisant data.table(sur les données de @ lando.carlissian):

library(data.table)
setDT(dat, key = "lab")
b.sample <- 
  replicate(B, dat[.(sample(unique(lab), replace = T)),
                   glm(y ~ x)$coefficients])

J'ai dû le faire récemment et utilisé dplyr. La solution n'est pas aussi élégante qu'avec data.table, mais:

library(dplyr)
replicate(B, {
  cluster_sample <- data.frame(cluster = sample(dat$cluster, replace = TRUE))
  dat_sample <- dat %>% inner_join(cluster_sample, by = 'cluster')
  coef(lm(y ~ x, data = dat_sample))
})

Le inner_joinrépète chaque ligne ayant une valeur donnée xde clusterpar le nombre de fois qui xapparaît dans cluster_sample.

Salut une solution très simple basée sur split et lapply, pas besoin de package spécifique sauf "boot", exemple avec une estimation de ICC basée sur la procédure nagakawa:

# FIRST FUNCTION : "parameter assesment"
nagakawa <- function(dataICC){
    #dataICC <- dbICC
    modele <- lmer(indicateur.L ~ 1 + (1|sujet.L) + (1|injection.L) + experience.L, data = dataICC)
    variance <- get_variance(modele)
    var.fixed <- variance$var.fixed
var.random <- variance$var.random
    var.sujet <- variance$var.intercept[1]
var.resid <- variance$var.residual
    icc.juge1 <- var.random / (var.random + var.fixed + var.resid)

    modele <- lmer(indicateur.L ~ 1 + (1 + injection.L|sujet.L) + experience.L, data = dataICC)
    variance <- VarCorr(modele)
    var.fixed <- get_variance_fixed(modele)
    var.random <- (attributes(variance$sujet.L)$stddev[1])^2 + (attributes(variance$sujet.L)$stddev[2])^2
    var.sujet <- (attributes(variance$sujet.L)$stddev[1])^2
    var.resid <- (attributes(variance)$sc)^2
icc.juge2 <- var.random / (var.random + var.fixed + var.resid)
return(c(as.numeric(icc.juge1),as.numeric(icc.juge2)))
  }
```
#SECOND FONCTION : bootstrap function, split on the hirarchical level as you want
```
  nagakawa.boot <- function(data,x){
list.ICC <- split(x = data, f = paste(data$juge.L,data$injection.L,sep = "_"))
    list.BOOT <- lapply(X = list.ICC, FUN = function(y){
      y[x,]
    })
    db.BOOT <- do.call(what = "rbind", args = list.BOOT)
    nagakawa(dataICC = db.BOOT)
  }

TROISIÈME: exécution du bootstrap

ICC.BOOT <- boot(data = dbICC, statistic = nagakawa.boot, R = 1000)