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

Question

У меня есть фрейм данных с несколькими переменными: регион, сезон, год, высота над уровнем моря и отклик (здесь пример):

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

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

subset(region&season&year)   and get  altitud~response
IT&wint&2013
IT&wint&2014
IT&spring&2013
IT&spring&2014

и так далее.Поэтому 24 комбинации.Есть идеи?

Заранее большое спасибо

Дэвид

Answer 1

Мое решение использует 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 (прогнозировать ответ по высоте?)

Answer 2

Для полноты изложения здесь приведены также решения с базой 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")

Answer 3

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

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)