запись результатов линейной регрессии повторно

Question

Как показано в следующем примере, я хочу добиться многократного запуска регрессии, каждый раз, когда R записывает оценки, сделанные в одном data.frame.

Каждый раз я менял условие года в «ifelse», т.е. ifelse (mydata $ year> = 1993, 1, 0), таким образом, каждый раз, когда я запускаю другую регрессию.

mydata$time = ifelse(mydata$year >= 1994, 1, 0)

Кто-нибудь может помочь? Мой основной код такой, как показано ниже (данные могут быть загружены через браузер, если R вернул ошибки):

library(foreign)
mydata = read.dta("http://dss.princeton.edu/training/Panel101.dta")
mydata$time = ifelse(mydata$year >= 1994, 1, 0)
mydata$did = mydata$time * mydata$treated
mydata$treated = ifelse(mydata$country == "E" | mydata$country == "F" | mydata$country == "G", 1, 0)
didreg = lm(y ~ treated + time + did, data = mydata)
summary(didreg)

Answer 1

Вот еще один вариант, здесь мы создаем матрицу для всех лет, присоединяем ее к mydata, собираем в длинные, объединяем в группы, затем запускаем регрессию для извлечения оценок. Обратите внимание, что «gt_et _ **» означает «больше или равно ..»

library(foreign)
library(dplyr)
library(tidyr)
library(purrr)
mydata = read.dta("http://dss.princeton.edu/training/Panel101.dta")

mtrx <- matrix(0, length(min(mydata$year):max(mydata$year)), length(min(mydata$year):max(mydata$year)))
mtrx[lower.tri(mtrx, diag = TRUE)] <- 1
df <- mtrx %>% as.data.frame() %>% mutate(year = min(mydata$year):max(mydata$year))
colnames(df) <- c(paste0("gt_et_", df$year), "year")
models <- df %>% 
  full_join(., mydata, by = "year") %>% 
  gather(mod, time, gt_et_1990:gt_et_1999) %>%
  nest(-mod) %>%
  mutate(data = map(data, ~mutate(.x, treated = ifelse(country == "E"|country == "F"|country == "G", 1, 0),
                                  did = time * treated)),
         mods = map(data, ~lm(y ~ treated + time + did, data = .x) %>% summary() %>% coef())) %>%
  unnest(mods %>% map(broom::tidy)) %>%
  filter(.rownames == "did") %>%
  select(-.rownames)

models
#>          mod    Estimate Std..Error     t.value   Pr...t..
#> 1 gt_et_1991 -2309823993 2410140350 -0.95837738 0.34137018
#> 2 gt_et_1992 -2036098728 1780081308 -1.14382344 0.25682856
#> 3 gt_et_1993 -2784221881 1504348732 -1.85078222 0.06867661
#> 4 gt_et_1994 -2519511630 1455676087 -1.73081886 0.08815711
#> 5 gt_et_1995 -2357323806 1455203186 -1.61992760 0.11001662
#> 6 gt_et_1996   250180589 1511322882  0.16553749 0.86902697
#> 7 gt_et_1997   405842197 1619653548  0.25057346 0.80292231
#> 8 gt_et_1998   -75683039 1852314277 -0.04085864 0.96753194
#> 9 gt_et_1999  2951694230 2452126428  1.20372840 0.23299421

Создано в 2018-09-01 по представ пакет (v0.2.0).

Answer 2

Как правило, если вы хотите повторять процесс много раз с разными входными данными каждый раз, вам нужна функция. Следующая функция принимает скалярное значение year_value в качестве входных данных, создает локальные переменные для регрессии и экспортирует оценки для модельного термина did.

foo <- function (year_value) {
  ## create local variables from `mydata`
  y <- mydata$y
  treated <- as.numeric(mydata$country %in% c("E", "F", "G"))  ## use `%in%`
  time <- as.numeric(mydata$year >= year_value)  ## use `year_value`
  did <- time * treated
  ## run regression using local variables
  didreg <- lm(y ~ treated + time + did)
  ## return estimate for model term `did`
  coef(summary(didreg))["did", ]
  }

foo(1993)
#     Estimate    Std. Error       t value      Pr(>|t|) 
#-2.784222e+09  1.504349e+09 -1.850782e+00  6.867661e-02 

Обратите внимание, что есть несколько мест, где ваш оригинальный код может быть улучшен. Скажем, с использованием "%in%" вместо нескольких "|" и использованием as.numeric вместо ifelse для приведения логического значения к числовому.

Теперь вам нужно что-то вроде цикла для итерации этой функции по нескольким year_value. Я бы использовал lappy.

## raw list of result from `lapply`
year_of_choice <- 1993:1994  ## taken for example
result <- lapply(year_of_choice, foo)

## rbind them into a matrix
data.frame(year = year_of_choice, do.call("rbind", result), check.names = FALSE)
#  year    Estimate Std. Error   t value   Pr(>|t|)
#1 1993 -2784221881 1504348732 -1.850782 0.06867661
#2 1994 -2519511630 1455676087 -1.730819 0.08815711

Обратите внимание, не включайте в качестве выбора год 1990 (минимум переменной year), в противном случае time будет вектором 1, таким же, как и перехват. Получившаяся модель имеет недостаток ранга, и вы получите ошибку «индекс за пределами». Версия R, начиная с 3.5.0, имеет новый аргумент complete для универсальной функции coef. Так что для стабильности мы можем использовать

coef(summary(didreg), complete = TRUE)["did", ]

Но вы должны увидеть все NA или NaN за 1990 год.