Подтвердить что ты не робот

Загрузочные иерархические/многоуровневые данные (передискретизирующие кластеры)

Я создаю script для создания выборок начальной загрузки из набора данных cats (из пакета -MASS-).

Следуя учебнику Дэвидсона и Хинкли [1], я провел простую линейную регрессию и принял фундаментальную непараметрическую процедуру для начальной загрузки из наблюдений iid, а именно: par resampling.

Оригинальный образец находится в форме:

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

Через одномерную линейную модель мы хотим объяснить вес кошачьих очагов через их вес мозга.

Код:

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)

Предположим теперь, что существует кластерная переменная cluster = 1, 2,..., 24 (например, каждая кошка принадлежит данному подстилку). Для простоты предположим, что данные сбалансированы: у нас есть 6 наблюдений для каждого кластера. Следовательно, каждый из 24 пометов состоит из 6 кошек (т.е. n_cluster = 6 и n = 144).

Можно создать поддельную переменную cluster через:

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

У меня есть два связанных вопроса:

Как имитировать образцы в соответствии с структурой (кластеризованным) набором данных? То есть как выполнить повторную выборку на уровне кластера? Я хотел бы пробовать кластеры с заменой и устанавливать наблюдения в каждом выбранном кластере, как в исходном наборе данных (т.е. Выборка с заменой кластеров и без замена наблюдений в каждом кластере).

Это стратегия, предложенная Дэвидсоном (стр. 100). Предположим, что мы рисуем B = 100 выборки. Каждый из них должен состоять из 24 возможных рекуррентных кластеров (например, 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), и каждый кластер должен содержать те же 6 наблюдений исходного набора данных. Как это сделать в R? (либо с пакетом -boot-, либо без него.) Есть ли у вас альтернативные предложения для продолжения?

Второй вопрос касается начальной модели регрессии. Предположим, что я использую модель фиксированных эффектов с перехватами на уровне кластера. Изменяется ли процедура повторной дискретизации?

[1] Davidson, A. C., Hinkley, D. V. (1997). Методы Bootstrap и их приложения. Пресса Кембриджского университета.

4b9b3361

Ответ 1

Если я правильно вас понимаю, это то, что вы пытаетесь сделать с c.data в качестве ввода:

  • Классы Resample с заменой
  • Поддерживать связь между каждым кластером в случайной выборке и ее точками из исходного набора данных (т.е. c.data)
  • Создайте бутстрап с выбранными кластерами

Вот script, которые достигают этого, который вы можете обернуть в функцию, чтобы повторить его R раз, где R - количество реплик bootstrap

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

# get a vector with all clusters
c <- sort(unique(c.data$cluster))

# group the data points per cluster
clust.group <- function(c) {
    c.data[c.data$cluster==c,]
}

clust.list <- lapply(c,clust.group)

# resample clusters with replacement
c.sample <- sample(c, replace=T)

clust.sample <- clust.list[c.sample]

clust.size <- 6

# combine the cluster list back to a single data matrix
clust.bind <- function(c) {
    matrix(unlist(c),nrow=clust.size)
}

c.boot <- do.call(rbind,lapply(clust.sample,clust.bind))

# Just to maintain columns name
colnames(c.boot) <- names(c.data)

# the new data set (single bootstrap replicate)
c.boot

Ответ 2

Я попытался решить проблему следующим образом. Хотя он работает, его можно, вероятно, улучшить с точки зрения скорости и "элегантности". Кроме того, если возможно, я предпочел бы найти способ использования пакета -boot-, поскольку он позволяет автоматически вычислять количество загрузочных доверительных интервалов через boot.ci...

Для простоты исходный набор данных состоит из 18 кошек ( "наблюдения нижнего уровня" ), вложенных в 6 лабораторий (переменная кластеризации). Набор данных сбалансирован (n_cluster = 3 для каждого кластера). У нас есть один регресс, x, для объяснения y.

Поддельный набор данных и матрица для хранения результатов:

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

На каждой итерации B следующие петлевые выборки 6 кластеров с заменой, каждая из которых состоит из 3 кошек, взятых без замены (т.е. внутренняя композиция кластеров поддерживается неизменной). Оценки коэффициента регрессора и его стандартной ошибки сохраняются в ранее созданной матрице:

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

  }

Конечная стандартная ошибка начальной загрузки:

 boot_se_x <- sum(b.sample[,3])/(B-1) 
 boot_se_x

Есть ли способ принять пакет -boot- для этого?

Кроме того, что касается принятия фиксированных эффектов на уровне кластера вместо простой линейной регрессии, то мое основное сомнение заключается в том, что в некоторых симулированных выборках может случиться так, что мы не выбрали некоторые из исходных кластеров, связанные с кластером перехваты не могут быть идентифицированы. Если вы посмотрите на код, который я опубликовал, это не должно быть проблемой с "механической" точки зрения (на каждой итерации мы можем поместиться в другую FE-модель с перехватами выборочных кластеров).

Мне было интересно, есть ли статистическая проблема во всех этих