Estou produzindo um script para criar exemplos de inicialização do catsconjunto de dados (do -MASS-pacote).

Seguindo o livro de Davidson e Hinkley [1], executei uma regressão linear simples e adotei um procedimento não paramétrico fundamental para o bootstrapping a partir de observações da iid, ou seja, pares de reamostragem .

A amostra original está no formato:

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

Através de um modelo linear univariado, queremos explicar o peso da lareira dos gatos através do seu peso cerebral.

O código é:

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)

Suponha agora que exista uma variável de agrupamento cluster = 1, 2,..., 24(por exemplo, cada gato pertence a uma ninhada). Para simplificar, suponha que os dados sejam equilibrados: temos 6 observações para cada cluster. Assim, cada uma das 24 ninhadas é composta por 6 gatos (ie n_cluster = 6e n = 144).

É possível criar uma clustervariável falsa através de:

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

Eu tenho duas perguntas relacionadas:

Como simular amostras de acordo com a estrutura do conjunto de dados (em cluster)? Ou seja, como reamostrar no nível do cluster? Gostaria de amostrar os clusters com substituição e definir as observações em cada cluster selecionado, como no conjunto de dados original (ou seja, amostragem com reposição dos clusters e sem substituir as observações em cada cluster).

Essa é a estratégia proposta por Davidson (p. 100). Suponha que desenhemos B = 100amostras. Cada um deles deve ser composto por 24 grupos possivelmente recorrentes (por exemplo 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), e cada grupo deve conter as mesmas 6 observações do conjunto de dados original. Como fazer isso R? (com ou sem o -boot-pacote.) Você tem sugestões alternativas para prosseguir?

A segunda pergunta diz respeito ao modelo de regressão inicial. Suponha que eu adote um modelo de efeitos fixos , com interceptações em nível de cluster. Isso altera o procedimento de reamostragem adotado?

[1] Davidson, AC, Hinkley, DV (1997). Métodos de inicialização e seus aplicativos . Cambridge University Press.

A reamostragem de todos os clusters é conhecida nas estatísticas da pesquisa desde que todos os métodos de reamostragem tenham sido utilizados (ou seja, desde meados da década de 1960), por isso é um método bem estabelecido. Veja minha coleção de links em http://www.citeulike.org/user/ctacmo/tag/survey_resampling . Se bootposso fazer isso ou não, eu não sei; Uso o surveypacote quando preciso trabalhar com autoinicializações de pesquisa, embora a última vez que verifiquei, ele não tivesse toda a funcionalidade necessária (como algumas pequenas correções de amostra, até onde me lembro).

Não acho que a aplicação de um modelo específico, como efeitos fixos, mude muito, mas, na IMO, a inicialização residual faz muitas suposições fortes (os resíduos são iid, o modelo está especificado corretamente). Cada um deles é facilmente quebrado, e a estrutura do cluster certamente quebra a suposição iid.

Tem havido alguma literatura econométrica sobre o bootstrap de cluster selvagem. Eles fingiram que trabalharam no vácuo sem todos aqueles cinquenta anos de pesquisa estatística sobre o tópico, então não tenho certeza do que fazer com ele.

Tentei resolver o problema sozinho e produzi o seguinte código.

Embora funcione, provavelmente poderia ser melhorado em termos de velocidade. Além disso, se possível, eu preferiria encontrar uma maneira de usar o -boot-pacote, pois ele permite calcular automaticamente vários intervalos de confiança de inicialização através de boot.ci...

Para simplificar, o conjunto de dados inicial consiste em 18 gatos (as observações de "nível inferior") aninhados em 6 laboratórios (a variável de agrupamento). O conjunto de dados é equilibrado ( n_cluster = 3para cada cluster). Temos um regressor x, para explicar y.

O conjunto de dados falso e a matriz onde armazenar resultados são:

  # 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)

Em cada uma das Biterações, o loop seguinte mostra 6 clusters com substituição, cada um composto por 3 gatos amostrados sem substituição (ou seja, a composição interna dos clusters é mantida inalterada). As estimativas do coeficiente de regressor e de seu erro padrão são armazenadas na matriz criada anteriormente:

  ####################################
  #   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

  }

Espero que isso ajude, Lando

Aqui está uma maneira muito mais simples (e quase indubitavelmente mais rápida) de data.tableexecutar o bootstrap usando (nos dados de @ lando.carlissian):

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

Eu tive que fazer isso recentemente e usado dplyr. A solução não é tão elegante quanto com data.table, mas:

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))
})

A inner_joinrepete cada linha tendo um determinado valor xde clusterpelo número de vezes que xaparece no cluster_sample.

Oi, uma solução muito simples, baseada em split e lapply, sem necessidade de pacote específico, exceto "boot", por exemplo, com uma estimativa do ICC com base no procedimento 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)
  }

TERCEIRO: execução de auto-inicialização

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