I have a data frame in R that looks like above
Id ln_W Year Exp
1 2.5 2010 15
1 2.3 2011 16
2 2.1 2010 20
3 2.5 2012 17
3 2.5 2013 18
I want to regress ln_W~Exp by each year in my dataset and save the summary of the results in a list format.
Does anyone know how to do this?
In base R, we split by the Year, loop over the list with lapply, create the model with lm and store the output as a list
out <- lapply(split(df1, df1$Year), function(x)
lm(ln_W ~ Exp, data = x))
NOTE: This doesn't require any packages
Or another option is lmList from lme4
library(lme4)
lmList(ln_W ~Exp | Year, data = df1)
#Call: lmList(formula = ln_W ~ Exp | Year, data = df1)
#Coefficients:
# (Intercept) Exp
#2010 3.7 -0.08
#2011 2.3 NA
#2012 2.5 NA
#2013 2.5 NA
#Degrees of freedom: 5 total; -3 residual
#Residual standard error: 0
df1 <- structure(list(Id = c(1L, 1L, 2L, 3L, 3L), ln_W = c(2.5, 2.3,
2.1, 2.5, 2.5), Year = c(2010L, 2011L, 2010L, 2012L, 2013L),
Exp = c(15L, 16L, 20L, 17L, 18L)), class = "data.frame",
row.names = c(NA,
-5L))
You can group by Year and then save lm summary as a list-column:
library(tidyverse)
df %>%
group_by(Year) %>%
summarise(fit = list(lm(ln_W ~ Exp, data = cur_data()) %>% summary))
Output:
# A tibble: 4 x 2
Year fit
<int> <list>
1 2010 <smmry.lm>
2 2011 <smmry.lm>
3 2012 <smmry.lm>
4 2013 <smmry.lm>
Get just the list of summaries by adding %>% pull(fit) to the chain.
(Note that, for the data provided, these summaries won't show much, just the intercept, as there are insufficient observations to fit.)
Why not using by.
res <- lapply(by(d, d$Year, lm, formula=ln_W ~ Exp), summary)
res
# $`2010`
#
# Call:
# FUN(formula = ..1, data = data[x, , drop = FALSE])
#
# Residuals:
# ALL 2 residuals are 0: no residual degrees of freedom!
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 3.70 NA NA NA
# Exp -0.08 NA NA NA
#
# Residual standard error: NaN on 0 degrees of freedom
# Multiple R-squared: 1, Adjusted R-squared: NaN
# F-statistic: NaN on 1 and 0 DF, p-value: NA
#
#
# $`2011`
#
# Call:
# FUN(formula = ..1, data = data[x, , drop = FALSE])
#
# Residuals:
# ALL 1 residuals are 0: no residual degrees of freedom!
#
# Coefficients: (1 not defined because of singularities)
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 2.3 NA NA NA
# Exp NA NA NA NA
#
# Residual standard error: NaN on 0 degrees of freedom
# [...]
Data:
d <- structure(list(Id = c(1L, 1L, 2L, 3L, 3L), ln_W = c(2.5, 2.3,
2.1, 2.5, 2.5), Year = c(2010L, 2011L, 2010L, 2012L, 2013L),
Exp = c(15L, 16L, 20L, 17L, 18L)), class = "data.frame", row.names = c(NA,
-5L))