Вопрос по r, plyr, regression, statistics-bootstrap – Блокировка начальной загрузки из списка тем

11

Я пытаюсь эффективно реализовать технику блочной загрузки, чтобы получить распределение коэффициентов регрессии. Основная схема выглядит следующим образом.

У меня есть набор данных панели, и я говорю, что фирма и год являются показателями. Для каждой итерации начальной загрузки я хочу выбрать n предметов с заменой. Из этого примера мне нужно построить новый фрейм данных, который являетсяrbind() составьте список всех наблюдений для каждого выбранного объекта, запустите регрессию и извлеките коэффициенты. Повторите для нескольких итераций, скажем, 100.

  • Each firm can potentially be selected multiple times, so I need to include it data multiple times in each iteration's data set.
  • Using a loop and subset approach, like below, seems computationally burdensome.
  • Note that for my real data frame, n, and the number iterations is much larger than the example below.

Сначала я хочу разбить существующий фрейм данных на список по предмету, используяsplit() команда. Оттуда используйте

sample(unique(df1$subject),n,replace=TRUE)

чтобы получить новый список, то, возможно, реализоватьquickdf отplyr пакет для создания нового фрейма данных.

Пример медленного кода:

require(plm)
data("Grunfeld", package="plm")

firms = unique(Grunfeld$firm)
n = 10
iterations = 100
mybootresults=list()

for(j in 1:iterations){

  v = sample(length(firms),n,replace=TRUE)
  newdata = NULL

  for(i in 1:n){
    newdata = rbind(newdata,subset(Grunfeld, firm == v[i]))
  }

  reg1 = lm(value ~ inv + capital, data = newdata)
  mybootresults[[j]] = coefficients(reg1)

}

mybootresults = as.data.frame(t(matrix(unlist(mybootresults),ncol=iterations)))
names(mybootresults) = names(reg1$coefficients)
mybootresults

  (Intercept)      inv    capital
1    373.8591 6.981309 -0.9801547
2    370.6743 6.633642 -1.4526338
3    528.8436 6.960226 -1.1597901
4    331.6979 6.239426 -1.0349230
5    507.7339 8.924227 -2.8661479
...
...
Это немного выдумка и, вероятно, излишнее использованиеboot, но я думаю, что ответ, который я отправил, делает работу. seancarmody
Вы смотрели наboot пакет? seancarmody
Да, я посмотрел и использовалboot пакет. Однако я не думаю, что у него есть возможность делать этот тип начальной загрузки блока. baha-kev
Обратите внимание, чтоboot позволяет легко запускать регрессии на нескольких ядрах, что может обеспечить значительное повышение скорости при большом количестве повторов. lmo

Ваш Ответ

5   ответов
0

который обычно должен быть быстрее принятого ответа, возвращает те же результаты и не полагаться на дополнительные пакеты (кромеboot). Ключ здесь заключается в использованииwhich и целочисленная индексация для создания каждой копии данных.split/subset а такжеdo.call/rbind.

# get function for boot
myIndex <- function(x, i) {
  # select the observations to subset. Likely repeated observations
  blockObs <- unlist(lapply(i, function(n) which(x[n] == Grunfeld$firm)))
  # run regression for given replicate, return estimated coefficients
  coefficients(lm(value~ inv + capital, data=Grunfeld[blockObs,]))
}

теперь, начальная загрузка

# get result
library(boot)
set.seed(1234)
b1 <- boot(firms, myIndex, 200)

Запустите принятый ответ

set.seed(1234)
b0 <- boot(firms, myfit, 200)

Давайте сравним глазное яблоко

используя индексацию

b1

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = firms, statistic = myIndex, R = 200)


Bootstrap Statistics :
       original      bias    std. error
t1* 410.8155650 -6.64885086 197.3147581
t2*   5.7598070  0.37922066   2.4966872
t3*  -0.6152727 -0.04468225   0.8351341

Оригинальная версия

b0

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = firms, statistic = myfit, R = 200)


Bootstrap Statistics :
       original      bias    std. error
t1* 410.8155650 -6.64885086 197.3147581
t2*   5.7598070  0.37922066   2.4966872
t3*  -0.6152727 -0.04468225   0.8351341

Они выглядят довольно близко. Теперь немного больше проверки

identical(b0$t, b1$t)
[1] TRUE

а также

identical(summary(b0), summary(b1))
[1] TRUE

Наконец, мы сделаем быстрый тест

library(microbenchmark)
microbenchmark(index={b1 <- boot(firms, myIndex, 200)}, 
              rbind={b0 <- boot(firms, myfit, 200)})

На моем компьютере это возвращается

Unit: milliseconds
  expr      min       lq     mean   median       uq      max neval
 index 292.5770 296.3426 303.5444 298.4836 301.1119 395.1866   100
 rbind 712.1616 720.0428 729.6644 724.0777 731.0697 833.5759   100

Таким образом, прямая индексация более чем в 2 раза быстрее на каждом уровне распределения.

note on missing fixed effects
Как и в большинстве ответов, проблема отсутствия «фиксированных эффектов» может появиться. Обычно фиксированные эффекты используются в качестве контроля, и исследователь интересуется одной или несколькими переменными, которые будут включены в каждое выбранное наблюдение. В этом доминирующем случае нет (или очень мало) вреда в ограничении возвращаемого результатаmyIndex или жеmyfit функция, которая включает только переменные, представляющие интерес в возвращаемом векторе.

1

чтобы управлять фиксированными эффектами.

library(boot)  # for boot
library(plm)   # for Grunfeld
library(dplyr) # for left_join

## Get the Grunfeld firm data (10 firms, each for 20 years, 1935-1954)
data("Grunfeld", package="plm")

## Create dataframe with uniq,ue firm identifier (one line per firm)
firms <- data.frame(firm=unique(Grunfeld$firm),junk=1)

## for boot(), X is the firms dataframe; i index the sampled firms
myfit <- function(X, i) {
    ## join the sampled firms to their firm-year data
    mydata <- left_join(X[i,], Grunfeld, by="firm")
    ## Distinguish between multiple resamples of the same firm
    ## Otherwise they have the same id in the fixed effects regression
    ## And trouble ensues
    mydata  <- mutate(group_by(mydata,firm,year),
                      firm_uniq4boot = paste(firm,"+",row_number())
                      )
    ## Run regression with and without firm fixed effects
    c(coefficients(lm(value ~ inv + capital, data = mydata)),
    coefficients(lm(value ~ inv + capital + factor(firm_uniq4boot), data = mydata)))
    }

set.seed(1)
system.time(b <- boot(firms, myfit, 1000))

summary(b)

summary(lm(value ~ inv + capital, data=Grunfeld))
summary(lm(value ~ inv + capital + factor(firm), data=Grunfeld))
1

dplyr::left_join это немного более кратко, занимает всего около 60% времени и дает те же результаты, что и в ответе Шона. Вот полный автономный пример.

library(boot)  # for boot
library(plm)   # for Grunfeld
library(dplyr) # for left_join

# First get the data
data("Grunfeld", package="plm")

myfit1 <- function(x, i) {
  # x is the vector of firms
  # i are the indexes into x
  mydata <- do.call("rbind", lapply(i, function(n) subset(Grunfeld, firm==x[n])))
  coefficients(lm(value ~ inv + capital, data = mydata))
}

myfit2 <- function(x, i) {
  # x is the vector of firms
  # i are the indexes into x
  mydata <- left_join(data.frame(firm=x[i]), Grunfeld, by="firm")
  coefficients(lm(value ~ inv + capital, data = mydata))
}

# rbind method
set.seed(1)
system.time(b1 <- boot(firms, myfit1, 5000))
  ##  user  system elapsed 
  ## 13.51    0.01   13.62 

# left_join method
set.seed(1)
system.time(b2 <- boot(firms, myfit2, 5000))
   ## user  system elapsed 
   ## 8.16    0.02    8.26 

summary(b1)
##      R  original bootBias    bootSE   bootMed
## 1 5000 410.81557 14.78272 195.62461 413.70175
## 2 5000   5.75981  0.49301   2.42879   6.00692
## 3 5000  -0.61527 -0.13134   0.78854  -0.76452
summary(b2)
##      R  original bootBias    bootSE   bootMed
## 1 5000 410.81557 14.78272 195.62461 413.70175
## 2 5000   5.75981  0.49301   2.42879   6.00692
## 3 5000  -0.61527 -0.13134   0.78854  -0.76452
4

tsboot функция вboot пакет с фиксированной схемой повторной выборки.

require(plm)
require(boot)
data(Grunfeld)

### each firm is of length 20
table(Grunfeld$firm)
##  1  2  3  4  5  6  7  8  9 10 
## 20 20 20 20 20 20 20 20 20 20


blockboot <- function(data) 
{
 coefficients(lm(value ~ inv + capital, data = data))

}

### fixed length (every 20 obs, so for each different firm) block bootstrap
set.seed(321)
boot.1 <- tsboot(Grunfeld, blockboot, R = 99, l = 20, sim = "fixed")

boot.1    
## Bootstrap Statistics :
##      original     bias    std. error
## t1* 410.81557 -25.785972    174.3766
## t2*   5.75981   0.451810      2.0261
## t3*  -0.61527   0.065322      0.6330

dim(boot.1$t)
## [1] 99  3

head(boot.1$t)
##        [,1]   [,2]      [,3]
## [1,] 522.11 7.2342 -1.453204
## [2,] 626.88 4.6283  0.031324
## [3,] 479.74 3.2531  0.637298
## [4,] 557.79 4.5284  0.161462
## [5,] 568.72 5.4613 -0.875126
## [6,] 379.04 7.0707 -1.092860
@ baha-kev Нет ... потому что первым шагом алгоритма является разбиение набора данных (из 200 строк) на 10 блоков по 20 строк (y_1, y_2, ..., y_10) и повторную выборку набора данных по блокам, например (y_1, y_10, y_3, ...., y_1), (y2, y_9, ..., y3) и так далее. А из-за организации набора данных Грюнфельда каждый y_ * соответствует одной фирме (каждые 20 строк). Для получения дополнительной информации посмотрите наBootstrap method and their Application (стр. 396).
l=20 в вашемtsboot() Заявление неизбежно займет блок, скажем, 15 наблюдений одной фирмы и 5 наблюдений другой (= 20), что я не хочу делать. Я либо хочу включить все наблюдения фирмы, либо ни одного, иногда включая наблюдения одной и той же фирмы несколько раз. baha-kev
14

myfit <- function(x, i) {
   mydata <- do.call("rbind", lapply(i, function(n) subset(Grunfeld, firm==x[n])))
   coefficients(lm(value ~ inv + capital, data = mydata))
}

firms <- unique(Grunfeld$firm)

b0 <- boot(firms, myfit, 999)
Спасибо, это было действительно полезно. Вы случайно не решили проблему с фиксированными эффектами? По этой же причине у меня возникают проблемы при сравнении моей целевой статистики по повторениям. Моя текущая стратегия состоит в том, чтобы опустить все коэффициенты FE в векторе, который я возвращаю в «myfit». функция, но это кажется немного дилетантским ..
Вы знаете, как действовать с фиксированными эффектами (т.е.factor(region)) вlm() вboot парадигма. Поскольку некоторые факторы могут не отображаться на каждой итерации, количество коэффициентов отличается, и это приводит к ошибкам. Мысли? baha-kev
Это, кажется, делает именно то, что я хочу - спасибо и дополнительный кредит за использованиеboot пакет. Я не пользоваласьdo.call раньше, так что это полезно и по этой причине. Cheers- baha-kev

Похожие вопросы