Bootstrapping jerárquica/multinivel de datos (remuestreo clusters)

Question

Estoy produciendo un script para la creación de muestras bootstrap de la cats conjunto de datos (de la -MASS- paquete).

Siguiendo las Davidson y Hinkley libro de texto [1] me encontré con una regresión lineal simple y adoptó una fundamental no paramétrico de procedimiento para el arranque del sistema desde el iid observaciones, es decir, pares de remuestreo.

La muestra original es de la forma:

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

A través de un modelo lineal univariante queremos explicar los gatos hogar de peso a través de su cerebro de peso.

El código es:

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)

Supongamos ahora que existe una agrupación de la variable cluster = 1, 2,..., 24 (por ejemplo, cada gato pertenece a una camada). Por simplicidad, supongamos que los datos son equilibrada: tenemos 6 observaciones para cada clúster. Por lo tanto, cada uno de los 24 camadas se compone de 6 gatos (es decir, n_cluster = 6 y n = 144).

Es posible crear un falso cluster variable a través de:

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

Tengo dos preguntas:

Cómo simular muestras de acuerdo con la (agrupado) conjunto de datos strucure? Es decir, cómo volver a muestrear en el nivel de clúster? Me gustaría muestra de los grupos con reemplazo y para establecer las observaciones dentro de cada grupo seleccionado, como en el conjunto de datos original (es decir, el muestreo con replacenment los clusters y sin el reemplazo de las observaciones dentro de cada grupo).

Esta es la estrategia propuesta por Davidson (p. 100). Supongamos que dibujar B = 100 de las muestras. Cada uno de ellos debe estar compuesto por 24 posiblemente recurrente grupos (por ejemplo 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), y cada grupo debe contener el mismo 6 de observaciones del conjunto de datos original. Cómo hacer que en R? (ya sea con o sin el -boot- paquete.) ¿Tienes sugerencias alternativas para el procedimiento?

La segunda cuestión se refiere a la inicial del modelo de regresión. Supongamos que adoptar un modelo de efectos fijos, con el grupo de nivel intercepta. ¿Cambia el remuestreo procedimiento adoptado?

[1] Davidson, A. C., Hinkley, D. V. (1997). Bootstrap métodos y sus aplicaciones. Cambridge University press.

Answer 1

Remuestreo de los racimos enteros ha sido conocida en la encuesta de estadísticas para como cualquier métodos de remuestreo se han utilizado en todo (que es, desde mediados de la década de 1960), por lo que es un método bien establecido. Ver mi colección de enlaces en http://www.citeulike.org/user/ctacmo/tag/survey_resampling. Si boot puede hacer esto o no, no sé; yo uso survey paquete cuando tengo que trabajar con la encuesta se levanta, aunque la última vez que lo comprobé, no tiene toda la funcionalidad que necesitaba (como algunos de muestras pequeñas correcciones, tan lejos como puedo recordar).

No creo que la aplicación de un modelo en particular, tales como efectos fijos cambia mucho las cosas, pero, en mi opinión, el residual de bootstrap hace un montón de fuertes supuestos (los residuos son yo.yo.d., el modelo está correctamente especificado). Cada uno de ellos se rompe con facilidad, y la estructura de cluster seguramente rompe el yo.yo.d. de la asunción.

Ha habido algunos econometría de la literatura sobre el clúster de bootstrap. Ellos fingieron que trabajó en el vacío sin todos esos cincuenta años de estadísticas encuesta de investigación en el tema, así que no estoy seguro de qué hacer con ella.

Answer 2

Traté de resolver el problema, y me produjo el siguiente código.

Aunque funciona, es probable que pudiera ser mejorado en términos de velocidad. También, si es posible, yo hubiera preferido encontrar un camino para el uso de la -boot- paquete, ya que permite calcular automáticamente un número de intervalos de confianza bootstrap a través de la boot.ci...

Para simplificar el proceso, a partir del conjunto de datos consiste en 18 de los gatos (el "nivel inferior" observaciones) anidado en 6 laboratorios (la agrupación de la variable). El conjunto de datos es equilibrada (n_cluster = 3 para cada grupo). Tenemos un regresor, x, para explicar y.

El falso conjunto de datos y la matriz donde almacenar los resultados son:

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

En cada una de las B de iteraciones, el siguiente bucle de muestras 6 grupos con reemplazo, cada uno compuesto por 3 gatos muestreados sin reemplazo (es decir, las agrupaciones internas de la composición se mantiene inalterada). Las estimaciones del coeficiente del regresor y de su error estándar se almacenan en los ya creados de la matriz:

  ####################################
  #   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 esto ayude, Lando

Answer 3

Aquí es mucho más fácil (y casi, sin duda, más rápido) manera de hacer el bootstrap usando data.table (en @lando.carlissian 's datos:

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