Bootstrap data hierarkis / multilevel (resampling cluster)

9

Saya memproduksi skrip untuk membuat sampel bootstrap dari catsdataset (dari -MASS-paket).

Mengikuti buku teks Davidson dan Hinkley [1] Saya menjalankan regresi linier sederhana dan mengadopsi prosedur dasar non-parametrik untuk bootstrap dari pengamatan awal, yaitu pemasangan kembali pasangan .

Sampel asli dalam bentuk:

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

Melalui model linear univariat, kami ingin menjelaskan berat badan kucing melalui berat otaknya.

Kode tersebut adalah:

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)

Anggaplah sekarang bahwa ada variabel pengelompokan cluster = 1, 2,..., 24(misalnya, masing-masing kucing milik sebuah sampah yang diberikan). Untuk kesederhanaan, anggaplah bahwa data seimbang: kami memiliki 6 pengamatan untuk setiap cluster. Oleh karena itu, masing-masing dari 24 liter terdiri dari 6 kucing (yaitu n_cluster = 6dan n = 144).

Dimungkinkan untuk membuat clustervariabel palsu melalui:

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

Saya punya dua pertanyaan terkait:

Bagaimana mensimulasikan sampel sesuai dengan struktur dataset (berkerumun)? Yaitu, bagaimana melakukan resample pada level cluster? Saya ingin sampel cluster dengan penggantian dan untuk mengatur pengamatan dalam setiap cluster yang dipilih seperti dalam dataset asli (yaitu pengambilan sampel dengan penggantian cluster dan tanpa penggantian pengamatan dalam masing-masing cluster).

Ini adalah strategi yang diusulkan oleh Davidson (hal. 100). Misalkan kita menggambar B = 100sampel. Masing-masing dari mereka harus dikomposisikan oleh 24 cluster yang mungkin berulang (misalnya 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), dan masing-masing cluster harus berisi 6 pengamatan yang sama dari dataset asli. Bagaimana cara melakukannya R? (baik dengan atau tanpa -boot-paket.) Apakah Anda memiliki saran alternatif untuk melanjutkan?

Pertanyaan kedua menyangkut model regresi awal. Misalkan saya mengadopsi model efek tetap , dengan intersep tingkat cluster. Apakah itu mengubah prosedur resampling yang diadopsi?

[1] Davidson, AC, Hinkley, DV (1997). Metode bootstrap dan aplikasinya . Pers Universitas Cambridge.

Stefano Lombardi
sumber

Jawaban:

9

Resampling seluruh cluster telah dikenal dalam statistik survei selama semua metode resampling telah digunakan di sana sama sekali (yaitu, sejak pertengahan 1960-an), sehingga merupakan metode yang mapan. Lihat koleksi tautan saya di http://www.citeulike.org/user/ctacmo/tag/survey_resampling . Apakah bootbisa melakukan ini atau tidak, saya tidak tahu; Saya menggunakan surveypaket ketika saya harus bekerja dengan bootstraps survei, meskipun terakhir kali saya memeriksa, itu tidak memiliki semua fungsi yang saya butuhkan (seperti beberapa koreksi sampel kecil, sejauh yang saya ingat).

Saya tidak berpikir menerapkan model tertentu seperti efek tetap mengubah banyak hal, tetapi, IMO, bootstrap residual membuat banyak asumsi yang kuat (residu iid, model ditentukan dengan benar). Masing-masing dari mereka mudah rusak, dan struktur cluster pasti mematahkan asumsi iid.

Sudah ada beberapa literatur ekonometrik tentang bootstrap cluster liar. Mereka berpura-pura bekerja dalam ruang hampa tanpa semua 50 tahun survei statistik penelitian menjadi topik, jadi saya tidak yakin apa yang harus dilakukan.

Tugas
sumber
Keraguan utama saya tentang membuat efek tetap pada tingkat klaster adalah bahwa dalam beberapa sampel simulasi dapat terjadi bahwa kami belum memilih beberapa kluster awal, sehingga intersep spesifik kluster terkait tidak dapat diidentifikasi. Jika Anda melihat kode yang saya posting, seharusnya tidak menjadi masalah dari sudut pandang "mekanis" (pada setiap iterasi kita dapat menyesuaikan model FE yang berbeda hanya dengan penyadapan cluster sampel '). Saya bertanya-tanya apakah ada masalah "statistik" dalam semua ini
Stefano Lombardi
3

Saya mencoba memecahkan masalah sendiri, dan saya menghasilkan kode berikut.

Meskipun berfungsi, mungkin bisa ditingkatkan dalam hal kecepatan. Juga, jika mungkin saya akan lebih suka untuk menemukan cara untuk menggunakan -boot-paket, karena memungkinkan untuk secara otomatis menghitung sejumlah interval kepercayaan bootstrap melalui boot.ci...

Untuk kesederhanaan, dataset awal terdiri dari 18 kucing (pengamatan "tingkat rendah") bersarang di 6 laboratorium (variabel pengelompokan). Dataset seimbang ( n_cluster = 3untuk setiap kluster). Kami memiliki satu regresi x,, untuk menjelaskan y.

Dataset palsu dan matriks tempat menyimpan hasil adalah:

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

Pada setiap Biterasi, loop berikut sampel 6 cluster dengan penggantian, masing-masing terdiri oleh 3 kucing sampel tanpa penggantian (yaitu komposisi internal cluster dipertahankan tidak berubah). Perkiraan koefisien regressor dan kesalahan standarnya disimpan dalam matriks yang dibuat sebelumnya:

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

  }

Semoga ini bisa membantu, Lando

Stefano Lombardi
sumber
menggunakan for for harus didominasi dengan menggunakan replicate; sebagai bonus itu secara otomatis mengembalikan b.samplearray untuk Anda. Selain itu, dengan semua penggabungan di sini, Anda hampir pasti lebih baik menggunakan data.tabledan mencoba kembali key. Saya dapat berkontribusi jawaban ketika saya sampai di komputer ... Pertanyaan: mengapa Anda melacak kesalahan standar koefisien?
MichaelChirico
Terima kasih @MichaelChirico, saya setuju. Jika saya ingat betul saya menyimpan kesalahan standar untuk merencanakan interval kepercayaan nanti.
Stefano Lombardi
tidakkah seharusnya interval kepercayaan hanya menjadi kuantil dari distribusi koefisien bootstrap? yaitu untuk interval kepercayaan 95%,quantile(b.sample[,2], c(.025, .975))
MichaelChirico
3

Berikut ini cara yang lebih sederhana (dan hampir pasti lebih cepat) untuk melakukan bootstrap menggunakan data.table(pada data @ lando.carlissian):

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

Saya harus melakukan ini baru-baru ini dan digunakan dplyr. Solusinya tidak seanggun dengan data.table, tetapi:

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

The inner_joinmengulangi setiap baris memiliki nilai yang diberikan xdari clusterdengan jumlah kali yang xmuncul dalam cluster_sample.

dash2
sumber
0

Hai solusi yang sangat sederhana berdasarkan split dan lapply, tidak perlu paket khusus kecuali "boot", contoh dengan estimasi ICC berdasarkan prosedur 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)
  }

KETIGA: eksekusi bootstrap

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