Подмножество данных по многим факторам одновременно в R

В моем решении используются функции broom с функциями tidy.

Чтение данных:

library(readr)

data <- read_table("region   season   year   altitud   response
IT       wint     2013   800       45
IT       wint     2013   815       47
IT       wint     2013   840       54
IT       wint     2014   800       49
IT       wint     2014   815       59")

Фактическое решение:

library(dplyr)
library(broom)
data_fit <- data %>%
    group_by(region, season, year) %>%
    do(fit = lm(altitud ~ response, data = .))

dfCoefs <- tidy(data_fit, fit)
dfCoefs

Что дает следующие коэффициенты регрессии для данных примера:

# A tibble: 4 x 8
# Groups:   region, season, year [2]
  region season  year term        estimate std.error statistic  p.value
  <chr>  <chr>  <dbl> <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 IT     wint    2013 (Intercept)   613.      34.7       17.7    0.0360
2 IT     wint    2013 response        4.22     0.711      5.93   0.106 
3 IT     wint    2014 (Intercept)   726.     NaN        NaN    NaN     
4 IT     wint    2014 response        1.5    NaN        NaN    NaN    

Однако вы хотите altitud ~ response (т.е. прогнозировать высоту по ответу) или response ~ altitud (предсказывать ответ по высоте?)

Надеюсь, я вас правильно понял, вот решение для мурлыки:

library(purrr)
library(dplyr)
nested<-df %>% 
  mutate_if(is.character,as.factor) %>% 
  group_by(year,season,region) %>% 
  nest()
my_model<-function(df){
  lm(altitud~response,data=df)
}

nested %>% 
  mutate(Mod=map(data,my_model)) 

Результат: Частично изменены данные для получения коэффициентов.

 A tibble: 3 x 5
   year season region data             Mod     
  <int> <fct>  <fct>  <list>           <list>  
1  2013 wint   IT     <tibble [3 x 2]> <S3: lm>
2  2014 wint   IT     <tibble [1 x 2]> <S3: lm>
3  2014 Summer IF     <tibble [1 x 2]> <S3: lm>

Прогнозирование с помощью modelr. Вы можете получить статистику, используя broom, как показано в другом ответе.

require(modelr)
nested %>% 
  mutate(Mod=map(data,my_model)) %>% 
  mutate(Preds=map2(data,Mod,add_predictions)) %>% 
  unnest(Preds)
# A tibble: 5 x 6
   year season region altitud response  pred
  <int> <fct>  <fct>    <int>    <int> <dbl>
1  2013 wint   IT         800       45  44.4
2  2013 wint   IT         815       47  47.9
3  2013 wint   IT         840       54  53.7
4  2014 wint   IT         800       49  49  
5  2014 Summer IF         815       59  59  

Получение сводной статистики с broom и purrr:

# A tibble: 4 x 8
   year season region term        estimate std.error statistic p.value
  <int> <fct>  <fct>  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1  2013 wint   IT     (Intercept) -140.      31.8        -4.40   0.142
2  2013 wint   IT     altitud        0.231    0.0389      5.93   0.106
3  2014 wint   IT     (Intercept)   49      NaN         NaN    NaN    
4  2014 Summer IF     (Intercept)   59      NaN         NaN    NaN

nested %>% 
  mutate(Mod=map(data,my_model)) %>% 
  mutate(Preds=map2(data,Mod,add_predictions),Tidy=map(Mod,tidy)) %>% 
  unnest(Tidy)

Данные:

df<-read.table(text="region   season   year   altitud   response
IT       wint     2013   800       45
               IT       wint     2013   815       47
               IT       wint     2013   840       54
               IT       wint     2014   800       49
               IF       Summer     2014   815       59",header=T)

Для полноты здесь также представлены база R и data.table решения.

База R

Одним из возможных базовых подходов R с использованием split() и lapply() был предложено Джого:

result <- lapply(split(DT, list(DT$region, DT$season, DT$year)), 
                 lm, formula = response ~ altitud)
print(result)

$IT.wint.2013

Call:
FUN(formula = ..1, data = X[[i]])

Coefficients:
(Intercept)      altitud  
  -140.0510       0.2306  


$IT.wint.2014

Call:
FUN(formula = ..1, data = X[[i]])

Coefficients:
(Intercept)      altitud  
  -484.3333       0.6667

Или, используя трубопровод для улучшения читаемости

library(magrittr)
result <- split(DT, list(DT$region, DT$season, DT$year)) %>% 
  lapply(lm, formula = response ~ altitud)

Таблица данных

С некоторой помощью broom:

library(data.table)
library(magrittr)
setDT(DT)[, lm(response ~ altitud, .SD) %>% broom::tidy(), by = .(region, season, year)]

   region season year        term     estimate   std.error statistic   p.value
1:     IT   wint 2013 (Intercept) -140.0510204 31.82553603 -4.400586 0.1422513
2:     IT   wint 2013     altitud    0.2306122  0.03888277  5.930962 0.1063382
3:     IT   wint 2014 (Intercept) -484.3333333         NaN       NaN       NaN
4:     IT   wint 2014     altitud    0.6666667         NaN       NaN       NaN

setDT(DT)[, lm(response ~ altitud, .SD) %>% broom::glance(), by = .(region, season, year)]

   region season year r.squared adj.r.squared    sigma statistic   p.value df    logLik      AIC    BIC deviance df.residual
1:     IT   wint 2013 0.9723576     0.9447152 1.111168  35.17631 0.1063382  2 -2.925132 11.85026 9.1461 1.234694           1
2:     IT   wint 2014 1.0000000           NaN      NaN       NaN       NaN  2       Inf     -Inf   -Inf 0.000000           0

Если вычисление lm() для разных групп занимает много времени, возможно, стоит сохранить результат и использовать его для последующих шагов обработки:

mod <- setDT(DT)[, .(model = .(lm(response ~ altitud, .SD))), by = .(region, season, year)]
mod

   region season year models
1:     IT   wint 2013   <lm>
2:     IT   wint 2014   <lm>

mod$models — это список моделей, эквивалентных result.

Теперь мы можем извлечь нужную информацию из вычисленных моделей, например,

mod[, models[[1]] %>% broom::tidy(), by = .(region, season, year)]

   region season year        term     estimate   std.error statistic   p.value
1:     IT   wint 2013 (Intercept) -140.0510204 31.82553603 -4.400586 0.1422513
2:     IT   wint 2013     altitud    0.2306122  0.03888277  5.930962 0.1063382
3:     IT   wint 2014 (Intercept) -484.3333333         NaN       NaN       NaN
4:     IT   wint 2014     altitud    0.6666667         NaN       NaN       NaN

Данные

library(data.table)
DT <- fread("
region   season   year   altitud   response
IT       wint     2013   800       45
IT       wint     2013   815       47
IT       wint     2013   840       54
IT       wint     2014   800       49
IT       wint     2014   815       59")