CODEWITHSUNDEEP
rloopsregressionlinear-regressionlm
Lernst
Lernst

Reputation: 13

record linear regression results repeatly

As shown in the following example, what I want to achieve is to run the regression many times, each time R records the estimates of did in one data.frame.

Each time, I changed the year condition in "ifelse", ie., ifelse(mydata$year >= 1993, 1, 0), thus each time I run a different regression. 

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

Can anyone help it? My basic code is as below (the data can be downloaded through browser if R returned errors):

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)

Upvotes: 1

Views: 213

Answers (2)

Zheyuan Li
Zheyuan Li

Reputation: 73385

Generally if you want to repeat a process many times with some different input each time, you need a function. The following function takes a scalar value year_value as its input, creates local variables for regression and exports estimates for model term 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

Note there are several places where your original code can be improved. Say, using "%in%" instead of multiple "|", and using as.numeric instead of ifelse to coerce boolean to numeric.

Now you need something like a loop to iterate this function over several different year_value. I would use 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

Note, don't include year 1990 (the minimum of variable year) as a choice, otherwise time will be a vector of 1, as same as the intercept. The resulting model is rank-deficient and you will get "subscript out of bounds" error. R version since 3.5.0 has a new complete argument to generic function coef. So for stability we may use

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

But you should see all NA or NaN for year 1990.

Upvotes: 1

AndS.
AndS.

Reputation: 8120

Here is another option, here we create a matrix for all the years, join it to mydata, gather to long, nest by grouping, then run regression to extract the estimates. Note that "gt_et_**" stands for "greater than or equal to.."

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

Created on 2018-09-01 by the reprex package (v0.2.0).

Upvotes: 0

Related Questions