dplyrに新たに実装されたらしいdo関数を使います.
なお現在のdplyrではdo()はsuperseded(非推奨だが動作する)扱いで,
reframe()やnest_by() +
mutate(list(...))を使うのが今風です.
この文書はdo()の練習なのでそのままdo()を使いますが,
新しく書くコードではreframe()などを検討してください.
group_byしてできたグループごとに関数適用した結果を見たいときに使えばいいのでしょうか.
…でも, dplyr::summariseでもいいですよね. 結果がsingle valueではない時に使うといいんでしょうかね?
とりあえず, vignetteをなぞってみます.
車のデータで遊びます.
mtcars %>%
group_by(cyl) %>%
do(head(., 2))## # A tibble: 6 × 11
## # Groups: cyl [3]
## mpg cyl disp hp drat wt qsec vs am gear carb
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 22.8 4 108 93 3.85 2.32 18.6 1 1 4 1
## 2 24.4 4 147. 62 3.69 3.19 20 1 0 4 2
## 3 21 6 160 110 3.9 2.62 16.5 0 1 4 4
## 4 21 6 160 110 3.9 2.88 17.0 0 1 4 4
## 5 18.7 8 360 175 3.15 3.44 17.0 0 0 3 2
## 6 14.3 8 360 245 3.21 3.57 15.8 0 0 3 4cyl数で分けられたdata.tableにhead(2)した結果がつながってますね.
ピリオドは自分自身.
次. 各cylごとに, データに線形回帰を施した結果を見たい時. group_by(cyl)してからlmをすると, lmの結果が入ったS3クラスのデータ(?)が返ってきます.
mtcars %>%
group_by(cyl) %>%
do(mod = lm(mpg ~ disp, data = .)) ## # A tibble: 3 × 2
## # Rowwise:
## cyl mod
## <dbl> <list>
## 1 4 <lm>
## 2 6 <lm>
## 3 8 <lm>data.tableでは, S3クラスのデータ(?)を持てません. 代わりに, list的なものになっているのですよね多分.
mtcars %>%
group_by(cyl) %>%
do(mod = lm(mpg ~ disp, data = .)) %>%
class## [1] "rowwise_df" "tbl_df" "tbl" "data.frame"tbl_dfってなんだろう. 分かりません. listみたいなものでしょう.
modにはcyl数で分けたデータごとにlmした結果(S3クラス)が入っているはず.
例えば決定係数を取り出してみます.
mtcars %>%
group_by(cyl) %>%
do(mod = lm(mpg ~ disp, data = .)) %>%
dplyr::summarise(cyl = cyl, rsq = summary(mod)$r.squared)## # A tibble: 3 × 2
## cyl rsq
## <dbl> <dbl>
## 1 4 0.648
## 2 6 0.0106
## 3 8 0.270できてますね. 次に, 係数を取り出してみます.
mtcars %>%
group_by(cyl) %>%
do(mod = lm(mpg ~ disp, data = .)) %>%
do(data.frame(cyl = .$cyl, var = names(coef(.$mod)), coef = coef(.$mod))) ## # A tibble: 6 × 3
## # Rowwise:
## cyl var coef
## <dbl> <chr> <dbl>
## 1 4 (Intercept) 40.9
## 2 4 disp -0.135
## 3 6 (Intercept) 19.1
## 4 6 disp 0.00361
## 5 8 (Intercept) 22.0
## 6 8 disp -0.0196なるほど.
2013年4月のメジャーリーグの打席結果データを使って遊びます. コードとデータはここにあります. do.Rmdを実行します. https://github.com/gghatano/analyze_mlbdata_with_R/tree/main/articles/clutch-hitting
dat = fread("../../data/derived/events-2013-04.csv")
head(dat)## GAME_ID AWAY_TEAM_ID INN_CT BAT_HOME_ID OUTS_CT BALLS_CT STRIKES_CT
## <char> <char> <int> <int> <int> <int> <int>
## 1: ANA201304090 OAK 1 0 0 0 1
## 2: ANA201304090 OAK 1 0 1 2 2
## 3: ANA201304090 OAK 1 0 2 3 1
## 4: ANA201304090 OAK 1 0 2 3 1
## 5: ANA201304090 OAK 1 0 2 3 1
## 6: ANA201304090 OAK 1 0 2 0 0
## PITCH_SEQ_TX AWAY_SCORE_CT HOME_SCORE_CT BAT_ID BAT_HAND_CD RESP_BAT_ID
## <char> <int> <int> <char> <char> <char>
## 1: CX 0 0 crisc001 R crisc001
## 2: CBCFBX 0 0 younc004 R younc004
## 3: BBCBB 0 0 lowrj001 R lowrj001
## 4: BBCBB 0 0 cespy001 R cespy001
## 5: BCBBX 0 0 norrd001 R norrd001
## 6: X 1 0 donaj001 R donaj001
## RESP_BAT_HAND_CD PIT_ID PIT_HAND_CD RESP_PIT_ID RESP_PIT_HAND_CD
## <char> <char> <char> <char> <char>
## 1: R wilsc004 L wilsc004 L
## 2: R wilsc004 L wilsc004 L
## 3: R wilsc004 L wilsc004 L
## 4: R wilsc004 L wilsc004 L
## 5: R wilsc004 L wilsc004 L
## 6: R wilsc004 L wilsc004 L
## POS2_FLD_ID POS3_FLD_ID POS4_FLD_ID POS5_FLD_ID POS6_FLD_ID POS7_FLD_ID
## <char> <char> <char> <char> <char> <char>
## 1: iannc001 pujoa001 kendh001 calla001 aybae001 troum001
## 2: iannc001 pujoa001 kendh001 calla001 aybae001 troum001
## 3: iannc001 pujoa001 kendh001 calla001 aybae001 troum001
## 4: iannc001 pujoa001 kendh001 calla001 aybae001 troum001
## 5: iannc001 pujoa001 kendh001 calla001 aybae001 troum001
## 6: iannc001 pujoa001 kendh001 calla001 aybae001 troum001
## POS8_FLD_ID POS9_FLD_ID BASE1_RUN_ID BASE2_RUN_ID BASE3_RUN_ID
## <char> <char> <char> <char> <char>
## 1: bourp001 hamij003
## 2: bourp001 hamij003
## 3: bourp001 hamij003
## 4: bourp001 hamij003 lowrj001
## 5: bourp001 hamij003 cespy001 lowrj001
## 6: bourp001 hamij003 norrd001 cespy001
## EVENT_TX LEADOFF_FL PH_FL BAT_FLD_CD BAT_LINEUP_ID EVENT_CD
## <char> <char> <char> <int> <int> <int>
## 1: 53/G56S T F 8 1 2
## 2: 63/G6 F F 9 2 2
## 3: W F F 6 3 14
## 4: W.1-2 F F 7 4 14
## 5: S8/G6M+.2-H;1-2 F F 2 5 20
## 6: S56/L5+.2-3;1-2 F F 5 6 20
## BAT_EVENT_FL AB_FL H_FL SH_FL SF_FL EVENT_OUTS_CT DP_FL TP_FL RBI_CT
## <char> <char> <int> <char> <char> <int> <char> <char> <int>
## 1: T T 0 F F 1 F F 0
## 2: T T 0 F F 1 F F 0
## 3: T F 0 F F 0 F F 0
## 4: T F 0 F F 0 F F 0
## 5: T T 1 F F 0 F F 1
## 6: T T 1 F F 0 F F 0
## WP_FL PB_FL FLD_CD BATTEDBALL_CD BUNT_FL FOUL_FL BATTEDBALL_LOC_TX ERR_CT
## <char> <char> <int> <char> <char> <char> <char> <int>
## 1: F F 5 G F F 56S 0
## 2: F F 6 G F F 6 0
## 3: F F 0 F F 0
## 4: F F 0 F F 0
## 5: F F 8 G F F 6M 0
## 6: F F 5 L F F 5 0
## ERR1_FLD_CD ERR1_CD ERR2_FLD_CD ERR2_CD ERR3_FLD_CD ERR3_CD BAT_DEST_ID
## <int> <char> <int> <char> <int> <char> <int>
## 1: 0 N 0 N 0 N 0
## 2: 0 N 0 N 0 N 0
## 3: 0 N 0 N 0 N 1
## 4: 0 N 0 N 0 N 1
## 5: 0 N 0 N 0 N 1
## 6: 0 N 0 N 0 N 1
## RUN1_DEST_ID RUN2_DEST_ID RUN3_DEST_ID BAT_PLAY_TX RUN1_PLAY_TX RUN2_PLAY_TX
## <int> <int> <int> <num> <char> <char>
## 1: 0 0 0 53
## 2: 0 0 0 63
## 3: 0 0 0 NA
## 4: 2 0 0 NA
## 5: 2 4 0 NA
## 6: 2 3 0 NA
## RUN3_PLAY_TX RUN1_SB_FL RUN2_SB_FL RUN3_SB_FL RUN1_CS_FL RUN2_CS_FL
## <char> <char> <char> <char> <char> <char>
## 1: F F F F F
## 2: F F F F F
## 3: F F F F F
## 4: F F F F F
## 5: F F F F F
## 6: F F F F F
## RUN3_CS_FL RUN1_PK_FL RUN2_PK_FL RUN3_PK_FL RUN1_RESP_PIT_ID
## <char> <char> <char> <char> <char>
## 1: F F F F
## 2: F F F F
## 3: F F F F
## 4: F F F F wilsc004
## 5: F F F F wilsc004
## 6: F F F F wilsc004
## RUN2_RESP_PIT_ID RUN3_RESP_PIT_ID GAME_NEW_FL GAME_END_FL PR_RUN1_FL
## <char> <char> <char> <char> <char>
## 1: T F F
## 2: F F F
## 3: F F F
## 4: F F F
## 5: wilsc004 F F F
## 6: wilsc004 F F F
## PR_RUN2_FL PR_RUN3_FL REMOVED_FOR_PR_RUN1_ID REMOVED_FOR_PR_RUN2_ID
## <char> <char> <char> <char>
## 1: F F
## 2: F F
## 3: F F
## 4: F F
## 5: F F
## 6: F F
## REMOVED_FOR_PR_RUN3_ID REMOVED_FOR_PH_BAT_ID REMOVED_FOR_PH_BAT_FLD_CD
## <char> <char> <int>
## 1: 0
## 2: 0
## 3: 0
## 4: 0
## 5: 0
## 6: 0
## PO1_FLD_CD PO2_FLD_CD PO3_FLD_CD ASS1_FLD_CD ASS2_FLD_CD ASS3_FLD_CD
## <int> <int> <int> <int> <int> <int>
## 1: 3 0 0 5 0 0
## 2: 3 0 0 6 0 0
## 3: 0 0 0 0 0 0
## 4: 0 0 0 0 0 0
## 5: 0 0 0 0 0 0
## 6: 0 0 0 0 0 0
## ASS4_FLD_CD ASS5_FLD_CD EVENT_ID
## <int> <int> <int>
## 1: 0 0 1
## 2: 0 0 2
## 3: 0 0 3
## 4: 0 0 4
## 5: 0 0 5
## 6: 0 0 6打者ごとに, ヒットを打つor打たないの系列に対して連検定(tseries::runs.test)を実行.
各打席にヒットを打つかどうかについて, ランダム性を検定してみます.
library(tseries)
source("../../R/retrosheet.R")
## 昔のfreadは"T"/"F"をlogicalとして読んでいたので AB_FL == "TRUE" で拾えました.
## 現在のfreadはcharacterとして読むので as_flag() で変換します.
## また, 打席数が少ない打者はruns.testが計算できないので除きます.
seqs =
dat %>%
dplyr::filter(as_flag(AB_FL)) %>% ## 四死球は除いて
mutate(HIT = as.factor(ifelse(H_FL > 0, "HIT", "NOHIT"))) %>%
group_by(BAT_ID) %>% ## 各打者ごとに
dplyr::filter(n() >= 50, dplyr::n_distinct(HIT) == 2)
runstest_res =
seqs %>%
## 現在のdo()はグループの列を裸の名前では参照できないので .$HIT と書きます
do(res = runs.test(.$HIT))
runstest_res## # A tibble: 255 × 2
## # Rowwise:
## BAT_ID res
## <chr> <list>
## 1 ackld001 <htest>
## 2 alony001 <htest>
## 3 altuj001 <htest>
## 4 alvap001 <htest>
## 5 amara001 <htest>
## 6 andre001 <htest>
## 7 ankir001 <htest>
## 8 aokin001 <htest>
## 9 arenj001 <htest>
## 10 avila001 <htest>
## # ℹ 245 more rows検定の結果が入ったS3クラスから, p値を取り出します.
## do()の結果はrowwise data frameなので, 1行ずつsummariseできます.
## もとのコードは summarise(runstest_res, pval = res$p.value) と
## 第1引数にデータそのものを渡してしまっていて動きませんでした.
runstest_res %>%
summarise(BAT_ID = BAT_ID, pval = res$p.value, .groups = "drop") %>%
arrange(pval) %>%
head(5)## # A tibble: 5 × 2
## BAT_ID pval
## <chr> <dbl>
## 1 gyorj001 0.00839
## 2 denoc001 0.0151
## 3 ellsj001 0.0182
## 4 pujoa001 0.0195
## 5 uptob001 0.0213do() は superseded です. 結果が1個で済むなら
summarise() で足ります.
## .by= は group_by 済みのデータには使えないので ungroup してから渡す
seqs |>
ungroup() |>
summarise(p = runs.test(HIT)$p.value, .by = BAT_ID) |>
arrange(p) |>
head(5)## # A tibble: 5 × 2
## BAT_ID p
## <chr> <dbl>
## 1 gyorj001 0.00839
## 2 denoc001 0.0151
## 3 ellsj001 0.0182
## 4 pujoa001 0.0195
## 5 uptob001 0.0213モデルのような「1個の値ではないもの」を持ち回りたいときは,
nest_by() とリスト列を使うのが今風です.
mtcars |>
nest_by(cyl) |>
mutate(mod = list(lm(mpg ~ disp, data = data))) |>
summarise(rsq = summary(mod)$r.squared, .groups = "drop")## # A tibble: 3 × 2
## cyl rsq
## <dbl> <dbl>
## 1 4 0.648
## 2 6 0.0106
## 3 8 0.270p値が小さい打者を並べて「この打者には流れがある」と読んではいけません. 250人ほどに同じ検定を繰り返しているので, 全員がランダムでも 10人強はp < 0.05になります.
この点をきちんと扱った分析は 打席結果に流れはあるか にあります.