2017-09-12 4 views
0

J'ai un problème pour lequel je ne pouvais pas trouver la solution après des heures de recherche, alors peut-être vous pouvez me aider sur celui-ci:comment écrire une boucle dans laquelle un ensemble de données est divisé et la pente de la ligne de tendance pour chaque division est donnée dans R

mon dataframe se présente comme suit:

stations_id phase_id refyear day 
140 10 1992 260 
140 10 1993 263 
140 10 1995 260 
140 10 1995 257 
140 12 1993 286 
140 12 1994 289 
140 12 1997 290 
150 10 1992 260 
150 10 1993 270 
150 10 1994 274 
165 15 1992 310 

le dataframe a environ 600 000 lignes et donc j'essaie désespérément de créer une boucle for qui met la pente de la ligne de régression avec "refyear" comme variable indépendante et "day" comme variable dépendante pour chaque combinaison de "stations_id" et "phase_id"; par conséquent, la division dépend de deux variables. Cependant, je ne trouve vraiment pas la solution et apprécierais vraiment si quelqu'un pouvait m'aider!

Meilleures salutations

Répondre

2

Utilisation dplyr et broom, vous pouvez modéliser refyear vs day par groupe sans avoir recours à une boucle et renvoyer une trame de données avec les coefficients du modèle. Dans le code ci-dessous, les coefficients de régression sont dans la colonne estimate. Les pentes de régression sont dans les rangées où term est égal à "jour".

library(tidyverse) 
library(broom) 

models = dat %>% group_by(stations_id, phase_id) %>% 
    do(tidy(lm(refyear ~ day, data=.))) 
stations_id phase_id  term  estimate std.error statistic  p.value 
     <int> <int>  <chr>  <dbl>  <dbl>  <dbl>  <dbl> 
1   140  10 (Intercept) 2080.4166667 94.44595383 22.0275891 0.002054594 
2   140  10   day -0.3333333 0.36324158 -0.9176629 0.455668946 
3   140  12 (Intercept) 1750.6923077 153.66666453 11.3927917 0.055736327 
4   140  12   day 0.8461538 0.53293871 1.5877132 0.357824750 
5   150  10 (Intercept) 1956.9230769 8.92887743 219.1678734 0.002904693 
6   150  10   day 0.1346154 0.03330867 4.0414519 0.154420958 
7   165  15 (Intercept) 1992.0000000   NaN   NaN   NaN 
0

Voici une tidyverse/solution purrr que je pense est plus propre que pour une boucle la version.

library(tidyverse) 
library(purrr) 
d <- read_csv("stations_id, phase_id, refyear, day 
140, 10, 1992, 260 
140, 10, 1993, 263 
140, 10, 1995, 260 
140, 10, 1995, 257 
140, 12, 1993, 286 
140, 12, 1994, 289 
140, 12, 1997, 290 
150, 10, 1992, 260 
150, 10, 1993, 270 
150, 10, 1994, 274 
165, 15, 1992, 310") 

nested <- d %>% 
    group_by(stations_id, phase_id) %>% 
    nest() 

nested <- nested %>% 
    mutate(mod = map(data, ~lm(day ~ refyear, data = .))) 

map(nested$mod, coef) 

[[1]] 
(Intercept)  refyear 
2032.2222222 -0.8888889 

[[2]] 
    (Intercept)  refyear 
-1399.4615385  0.8461538 

[[3]] 
(Intercept)  refyear 
    -13683   7 

[[4]] 
(Intercept)  refyear 
     310   NA 
0

Vous pouvez utiliser tidyverse pour cela.

Premier groupe par les variables, puis tidyr::nest les données de ce groupe. Vous avez maintenant une liste-colonne qui contient toutes les données des variables non-groupantes pour chaque combinaison des variables de regroupement.

Ensuite, vous pouvez utiliser purrr::map au sein de dplyr::mutate pour parcourir la nouvelle colonne de liste correspondant à votre modèle sur chaque daraframe distincte dans la liste-col. Vous avez maintenant un list-col supplémentaire contenant les modèles. Vous pouvez ensuite itérer sur ceux-ci, en prenant le coefficient désiré de chaque modèle.

Enfin, vous pouvez simplement sélectionner la pente et vous obtiendrez une seule ligne pour chaque combinaison de vars de regroupement avec la pente de vos modèles. Ou vous pouvez ajouter les données et ajouter la pente en tant que nouvelle colonne qui se répète pour toutes les valeurs des variables de regroupement.

Pour guide plus détaillé à ces sortes de flux de travail vérifier le chapitre sur many models de R for Data Science


library(tidyverse) 

nested <- mtcars %>% 
    select(cyl, mpg, wt) %>% 
    group_by(cyl) %>% 
    nest() 

#> # A tibble: 3 x 2 
#>  cyl    data 
#> <dbl>   <list> 
#> 1  6 <tibble [7 x 2]> 
#> 2  4 <tibble [11 x 2]> 
#> 3  8 <tibble [14 x 2]> 

models <- nested %>% 
    mutate(
    model = map(data, ~lm(mpg ~ wt, data = .x)), 
    slope = map_dbl(model, c("coefficients", "wt")) 
) 

#> # A tibble: 3 x 4 
#>  cyl    data model  slope 
#> <dbl>   <list> <list>  <dbl> 
#> 1  6 <tibble [7 x 2]> <S3: lm> -2.780106 
#> 2  4 <tibble [11 x 2]> <S3: lm> -5.647025 
#> 3  8 <tibble [14 x 2]> <S3: lm> -2.192438 

models %>% select(cyl, slope) 

#> # A tibble: 3 x 2 
#>  cyl  slope 
#> <dbl>  <dbl> 
#> 1  6 -2.780106 
#> 2  4 -5.647025 
#> 3  8 -2.192438 

models %>% select(-model) %>% unnest() 

#> # A tibble: 32 x 4 
#>  cyl  slope mpg wt 
#> <dbl>  <dbl> <dbl> <dbl> 
#> 1  6 -2.780106 21.0 2.620 
#> 2  6 -2.780106 21.0 2.875 
#> 3  6 -2.780106 21.4 3.215 
#> 4  6 -2.780106 18.1 3.460 
#> 5  6 -2.780106 19.2 3.440 
#> 6  6 -2.780106 17.8 3.440 
#> 7  6 -2.780106 19.7 2.770 
#> 8  4 -5.647025 22.8 2.320 
#> 9  4 -5.647025 24.4 3.190 
#> 10  4 -5.647025 22.8 3.150 
#> # ... with 22 more rows