to be continued...

Rや心理統計学の備忘録的な

There are three kinds of lies: lies, damned lies, and statistics.

- 嘘には三種類ある。嘘、大嘘、そして統計だ -
Benjamn Disraeli(19世紀のイギリス首相)

cmdstanrで推定後にWAICを計算

結論

この関数を使おう ※stanファイルの対数尤度の変数名がlog_likとして推定している事を仮定

waic2 <- function(fit){
  # log_lik <- rstan::extract(fit)$log_lik # rtan用
  log_lik <- fit$draws("log_lik") %>% as_draws_df() %>% select(contains("log_lik")) # cmdstanr用
  lppd <- sum(log(colMeans(exp(log_lik))))
  p_waic <- sum(apply(log_lik,2,var))
  waic <- -2*lppd+2*p_waic
  return(list(waic=waic,lppd=lppd,p_waic=p_waic))
}

小噺

Rのcmdstanrパッケージを使ってMCMCをしたときの推定結果は、rstanパッケージのそれとは異なります。rstanで推定した場合はstanfitオブジェクト、cmdstanrで推定した場合はCmdStanFitオブジェクトとなります。

対数尤度からWAIC((Watanabe-Akaike) Widely Applicable Information Criterion)を計算するとき、looパッケージを使うこともできますが、looパッケージのwaic()は引数がstanfitなので、cmdstanrで推定した場合、stanfitに直してあげる必要があります。

なので、自分はこちらの記事にある自作関数で計算します。

ただ、この記事のwaic2()もコードの通り、stanfitを引数としているので、CmdStanFit用に個人的に直しました。

補足

stanfitオブジェクトから特定の推定結果を抜き出すときrstan::extract()でパラメータを指定します。この時の出力は2次元配列(パラメータ×mcmcサンプル総数)です。 同様の処理をcmdstanrの時は、fit$draws("log_lik")で行えます。 しかし、cmdstanrだと3次元配列(パラメータ×チェイン×iteration)で出力されます。

なので、as_draws_ds()でdata.frame型に変換します。この時、2次元データに落とし込めますが、log_likの他に'.chain', '.iteration', '.draw'という変数もついてきます。これが便利かどうかはおいておいて、select()で除外してあげれば、あとはstanfitのときと同じ処理でいける。

以上で終わりだ それだけ

to be continued ...

俺的{dplyr1.0.0}メモ~select(),rename(),relocate()編~

最近は、水曜どうでしょうにハマっています。対決列島面白いですね

今更、{dplyr1.0.0}をキャッチアップしていこうのコーナー

今回はカラム(変数)にまつわるselect(),rename(),relocate()の話

{dplyr1.0.0}select(),rename(),relocate()で追加された俺的に特徴的な機能はこちら

  1. ~if(), ~at()を使わなくてよくなった‼
  2. 関数を使った変数名変更が可能に‼
  3. カラム(変数)界のarrange()...relocate()の登場‼

※カラム = 変数です。呼び方が安定しなくてすみません

準備

今回はdplyr::starwarsデータを使います。映画starwarsに登場するキャラクターの情報が入っているデータです(当方、スターウォーズライトセーバーが赤いと敵くらいしか知りません)。

> library(dplyr,warn.conflicts = F)
> dplyr::starwars
# A tibble: 87 x 14
   name  height  mass hair_color skin_color eye_color birth_year sex  
   <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr>
 1 Luke~    172    77 blond      fair       blue            19   male 
 2 C-3PO    167    75 NA         gold       yellow         112   none 
 3 R2-D2     96    32 NA         white, bl~ red             33   none 
 4 Dart~    202   136 none       white      yellow          41.9 male 
 5 Leia~    150    49 brown      light      brown           19   fema~
 6 Owen~    178   120 brown, gr~ light      blue            52   male 
 7 Beru~    165    75 brown      light      blue            47   fema~
 8 R5-D4     97    32 NA         white, red red             NA   none 
 9 Bigg~    183    84 black      light      brown           24   male 
10 Obi-~    182    77 auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

1. ~if(), ~at()を使わなくてよくなった‼

今まで、select(),rename()のラッパーとしてあったこれらの関数ですが、使わなくてよくなりました(使えなくはなっておらず、非推奨になった感じです)。

dplyr1.0.0ではどうするかというと、where()any_of(),all_of()という関数を使うことで楽に変数指定できるようになりました

select_if()select(where())

例えば、データフレーム内の数値型データだけ抜き出したとき

> starwars %>% 
+   select(where(is.numeric))
# A tibble: 87 x 3
   height  mass birth_year
    <int> <dbl>      <dbl>
 1    172    77       19  
 2    167    75      112  
 3     96    32       33  
 4    202   136       41.9
 5    150    49       19  
 6    178   120       52  
 7    165    75       47  
 8     97    32       NA  
 9    183    84       24  
10    182    77       57  
# ... with 77 more rows

select(where(~~))といった指定をすることで、従来のselect_if()と同じ出力を得ることができます。これが楽かどうかは置いておいて、これにより予期せぬエラーを回避できるそうです(詳しくはコチラのupdate notice読んでね)。

もちろん、他の関数や条件式を使って変数指定することもできます

# factor型以外のデータの抽出
> starwars %>% select(!where(is.factor))
# A tibble: 87 x 14
   name  height  mass hair_color skin_color eye_color birth_year sex  
   <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr>
 1 Luke~    172    77 blond      fair       blue            19   male 
 2 C-3PO    167    75 NA         gold       yellow         112   none 
 3 R2-D2     96    32 NA         white, bl~ red             33   none 
 4 Dart~    202   136 none       white      yellow          41.9 male 
 5 Leia~    150    49 brown      light      brown           19   fema~
 6 Owen~    178   120 brown, gr~ light      blue            52   male 
 7 Beru~    165    75 brown      light      blue            47   fema~
 8 R5-D4     97    32 NA         white, red red             NA   none 
 9 Bigg~    183    84 black      light      brown           24   male 
10 Obi-~    182    77 auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

# 数値型データかつ変数の頭文字が"h"の変数を抽出
> starwars %>% select(where(is.numeric) & starts_with("h"))
# A tibble: 87 x 1
   height
    <int>
 1    172
 2    167
 3     96
 4    202
 5    150
 6    178
 7    165
 8     97
 9    183
10    182
# ... with 77 more rows

select_at()select(any_of()) or select(all_of())

変数名の 文字ベクトルを受け渡して変数指定するときは

# 変数名の文字ベクトル作成
> vars <- c("name","mass","skin_color","sex")
# any_of()は部分一致:文字ベクトルに一致するものがあれば抜き出す
> starwars %>% select(any_of(vars))
# A tibble: 87 x 4
   name                mass skin_color  sex   
   <chr>              <dbl> <chr>       <chr> 
 1 Luke Skywalker        77 fair        male  
 2 C-3PO                 75 gold        none  
 3 R2-D2                 32 white, blue none  
 4 Darth Vader          136 white       male  
 5 Leia Organa           49 light       female
 6 Owen Lars            120 light       male  
 7 Beru Whitesun lars    75 light       female
 8 R5-D4                 32 white, red  none  
 9 Biggs Darklighter     84 light       male  
10 Obi-Wan Kenobi        77 fair        male  
# ... with 77 more rows

このように、select()の中にany_of()を使うことでselect_at()と同じ出力を得ることができます

any_of()は部分一致、all_of()は完全一致です。なので、文字ベクトルの中にデータフレームのカラムと符合しない要素があると

# starwarsの中に"umr"というカラムはない
> vars2 <- c("name","mass","skin_color","sex","umr")
> starwars %>% select(all_of(vars2))
 エラー: Can't subset columns that don't exist.
x Column `umr` doesn't exist.

エラーを返しつつ、「"umr"ってカラムはないよ」と教えてくれます。

基本的にany_of()を使えばいいと思いますが、厳密性を持たせるときとかカラムを探索するときはall_of()がつかえるのかな?

select()の指定方法として当然ですが、文字ベクトルの変数以外を抜き出したいときは、「-」を使えばおけです。

> starwars %>% select(-any_of(vars))
# A tibble: 87 x 10
   height hair_color eye_color birth_year gender homeworld species
    <int> <chr>      <chr>          <dbl> <chr>  <chr>     <chr>  
 1    172 blond      blue            19   mascu~ Tatooine  Human  
 2    167 NA         yellow         112   mascu~ Tatooine  Droid  
 3     96 NA         red             33   mascu~ Naboo     Droid  
 4    202 none       yellow          41.9 mascu~ Tatooine  Human  
 5    150 brown      brown           19   femin~ Alderaan  Human  
 6    178 brown, gr~ blue            52   mascu~ Tatooine  Human  
 7    165 brown      blue            47   femin~ Tatooine  Human  
 8     97 NA         red             NA   mascu~ Tatooine  Droid  
 9    183 black      brown           24   mascu~ Tatooine  Human  
10    182 auburn, w~ blue-gray       57   mascu~ Stewjon   Human  
# ... with 77 more rows, and 3 more variables: films <list>,
#   vehicles <list>, starships <list>

2. 関数を使った変数変更が可能に‼

変数名指定に使われるrename()ですが、新たにrename_with()が追加されました。

rename_with()では変更したい変数の指定・変数名の変更方法を関数を使って記述することが可能です

変数名を大文字にしてみる

> starwars %>% rename_with(toupper)
# A tibble: 87 x 14
   NAME  HEIGHT  MASS HAIR_COLOR SKIN_COLOR EYE_COLOR BIRTH_YEAR SEX  
   <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr>
 1 Luke~    172    77 blond      fair       blue            19   male 
 2 C-3PO    167    75 NA         gold       yellow         112   none 
 3 R2-D2     96    32 NA         white, bl~ red             33   none 
 4 Dart~    202   136 none       white      yellow          41.9 male 
 5 Leia~    150    49 brown      light      brown           19   fema~
 6 Owen~    178   120 brown, gr~ light      blue            52   male 
 7 Beru~    165    75 brown      light      blue            47   fema~
 8 R5-D4     97    32 NA         white, red red             NA   none 
 9 Bigg~    183    84 black      light      brown           24   male 
10 Obi-~    182    77 auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: GENDER <chr>,
#   HOMEWORLD <chr>, SPECIES <chr>, FILMS <list>, VEHICLES <list>,
#   STARSHIPS <list>

rename_with(toupper)で瞬殺です。

{tidyselect}where(),any_of()を使って変数名を変更してみましょう

# 変数名の最初が"n"の変数だけ大文字に
> starwars %>% rename_with(toupper, starts_with("n"))
# A tibble: 87 x 14
   NAME  height  mass hair_color skin_color eye_color birth_year sex  
   <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr>
 1 Luke~    172    77 blond      fair       blue            19   male 
 2 C-3PO    167    75 NA         gold       yellow         112   none 
 3 R2-D2     96    32 NA         white, bl~ red             33   none 
 4 Dart~    202   136 none       white      yellow          41.9 male 
 5 Leia~    150    49 brown      light      brown           19   fema~
 6 Owen~    178   120 brown, gr~ light      blue            52   male 
 7 Beru~    165    75 brown      light      blue            47   fema~
 8 R5-D4     97    32 NA         white, red red             NA   none 
 9 Bigg~    183    84 black      light      brown           24   male 
10 Obi-~    182    77 auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

# 数値型変数名だけ大文字に
> starwars %>% rename_with(toupper, where(is.numeric))
# A tibble: 87 x 14
   name  HEIGHT  MASS hair_color skin_color eye_color BIRTH_YEAR sex  
   <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr>
 1 Luke~    172    77 blond      fair       blue            19   male 
 2 C-3PO    167    75 NA         gold       yellow         112   none 
 3 R2-D2     96    32 NA         white, bl~ red             33   none 
 4 Dart~    202   136 none       white      yellow          41.9 male 
 5 Leia~    150    49 brown      light      brown           19   fema~
 6 Owen~    178   120 brown, gr~ light      blue            52   male 
 7 Beru~    165    75 brown      light      blue            47   fema~
 8 R5-D4     97    32 NA         white, red red             NA   none 
 9 Bigg~    183    84 black      light      brown           24   male 
10 Obi-~    182    77 auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

# 文字ベクトル(vars)の変数だけ大文字に
> starwars %>% rename_with(toupper, all_of(vars))
# A tibble: 87 x 14
   NAME  height  MASS hair_color SKIN_COLOR eye_color birth_year SEX  
   <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr>
 1 Luke~    172    77 blond      fair       blue            19   male 
 2 C-3PO    167    75 NA         gold       yellow         112   none 
 3 R2-D2     96    32 NA         white, bl~ red             33   none 
 4 Dart~    202   136 none       white      yellow          41.9 male 
 5 Leia~    150    49 brown      light      brown           19   fema~
 6 Owen~    178   120 brown, gr~ light      blue            52   male 
 7 Beru~    165    75 brown      light      blue            47   fema~
 8 R5-D4     97    32 NA         white, red red             NA   none 
 9 Bigg~    183    84 black      light      brown           24   male 
10 Obi-~    182    77 auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

このように、rename_with(toupper, 変更変数の指定)といった記述で変更できます。

他にもgsub()を用いた正規表現を使った変数名の変更などもできますがここでは省略


3. カラム(変数)界のarrange()...relocate()の登場‼

変数名の順序を変えたいとき困ってましたよね...そこまで重要な作業ではないのですが、あんなこといいなできたらいいな処理です。従来は、select()を使って変数名の順序を変更することは可能でしたが、もうそんなハックみあふれることをしなくて大丈夫です。

relocate()は指定した変数をデータフレームの左に持ってきます

> starwars %>% relocate(skin_color,birth_year)
# A tibble: 87 x 14
   skin_color birth_year name  height  mass hair_color eye_color sex  
   <chr>           <dbl> <chr>  <int> <dbl> <chr>      <chr>     <chr>
 1 fair             19   Luke~    172    77 blond      blue      male 
 2 gold            112   C-3PO    167    75 NA         yellow    none 
 3 white, bl~       33   R2-D2     96    32 NA         red       none 
 4 white            41.9 Dart~    202   136 none       yellow    male 
 5 light            19   Leia~    150    49 brown      brown     fema~
 6 light            52   Owen~    178   120 brown, gr~ blue      male 
 7 light            47   Beru~    165    75 brown      blue      fema~
 8 white, red       NA   R5-D4     97    32 NA         red       none 
 9 light            24   Bigg~    183    84 black      brown     male 
10 fair             57   Obi-~    182    77 auburn, w~ blue-gray male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

name,height...という並びでしたが、実行後はskin_color,birth_yearが最も左に来ています。

勿論、{tidyselect}where()などを使った指定方法も可能です

> starwars %>% relocate(starts_with("n"))
# A tibble: 87 x 14
   name  height  mass hair_color skin_color eye_color birth_year sex  
   <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr>
 1 Luke~    172    77 blond      fair       blue            19   male 
 2 C-3PO    167    75 NA         gold       yellow         112   none 
 3 R2-D2     96    32 NA         white, bl~ red             33   none 
 4 Dart~    202   136 none       white      yellow          41.9 male 
 5 Leia~    150    49 brown      light      brown           19   fema~
 6 Owen~    178   120 brown, gr~ light      blue            52   male 
 7 Beru~    165    75 brown      light      blue            47   fema~
 8 R5-D4     97    32 NA         white, red red             NA   none 
 9 Bigg~    183    84 black      light      brown           24   male 
10 Obi-~    182    77 auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>


> starwars %>% relocate(where(is.character))
# A tibble: 87 x 14
   name  hair_color skin_color eye_color sex   gender homeworld
   <chr> <chr>      <chr>      <chr>     <chr> <chr>  <chr>    
 1 Luke~ blond      fair       blue      male  mascu~ Tatooine 
 2 C-3PO NA         gold       yellow    none  mascu~ Tatooine 
 3 R2-D2 NA         white, bl~ red       none  mascu~ Naboo    
 4 Dart~ none       white      yellow    male  mascu~ Tatooine 
 5 Leia~ brown      light      brown     fema~ femin~ Alderaan 
 6 Owen~ brown, gr~ light      blue      male  mascu~ Tatooine 
 7 Beru~ brown      light      blue      fema~ femin~ Tatooine 
 8 R5-D4 NA         white, red red       none  mascu~ Tatooine 
 9 Bigg~ black      light      brown     male  mascu~ Tatooine 
10 Obi-~ auburn, w~ fair       blue-gray male  mascu~ Stewjon  
# ... with 77 more rows, and 7 more variables: species <chr>,
#   height <int>, mass <dbl>, birth_year <dbl>, films <list>,
#   vehicles <list>, starships <list>


> starwars %>% relocate(name,where(is.numeric))
# A tibble: 87 x 14
   name  height  mass birth_year hair_color skin_color eye_color sex  
   <chr>  <int> <dbl>      <dbl> <chr>      <chr>      <chr>     <chr>
 1 Luke~    172    77       19   blond      fair       blue      male 
 2 C-3PO    167    75      112   NA         gold       yellow    none 
 3 R2-D2     96    32       33   NA         white, bl~ red       none 
 4 Dart~    202   136       41.9 none       white      yellow    male 
 5 Leia~    150    49       19   brown      light      brown     fema~
 6 Owen~    178   120       52   brown, gr~ light      blue      male 
 7 Beru~    165    75       47   brown      light      blue      fema~
 8 R5-D4     97    32       NA   NA         white, red red       none 
 9 Bigg~    183    84       24   black      light      brown     male 
10 Obi-~    182    77       57   auburn, w~ fair       blue-gray male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

.before .after引数を使った指定

.before,.after引数を使えば、任意の変数を別の任意の変数の前後に持ってくることができます。

IDと併せて、当該の変数情報を見たいときとかに便利ですかね?

# before
> starwars %>% relocate(name, .before = hair_color)
# A tibble: 87 x 14
   height  mass name  hair_color skin_color eye_color birth_year sex  
    <int> <dbl> <chr> <chr>      <chr>      <chr>          <dbl> <chr>
 1    172    77 Luke~ blond      fair       blue            19   male 
 2    167    75 C-3PO NA         gold       yellow         112   none 
 3     96    32 R2-D2 NA         white, bl~ red             33   none 
 4    202   136 Dart~ none       white      yellow          41.9 male 
 5    150    49 Leia~ brown      light      brown           19   fema~
 6    178   120 Owen~ brown, gr~ light      blue            52   male 
 7    165    75 Beru~ brown      light      blue            47   fema~
 8     97    32 R5-D4 NA         white, red red             NA   none 
 9    183    84 Bigg~ black      light      brown           24   male 
10    182    77 Obi-~ auburn, w~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

# after
> starwars %>% relocate(name, .after = hair_color)
# A tibble: 87 x 14
   height  mass hair_color name  skin_color eye_color birth_year sex  
    <int> <dbl> <chr>      <chr> <chr>      <chr>          <dbl> <chr>
 1    172    77 blond      Luke~ fair       blue            19   male 
 2    167    75 NA         C-3PO gold       yellow         112   none 
 3     96    32 NA         R2-D2 white, bl~ red             33   none 
 4    202   136 none       Dart~ white      yellow          41.9 male 
 5    150    49 brown      Leia~ light      brown           19   fema~
 6    178   120 brown, gr~ Owen~ light      blue            52   male 
 7    165    75 brown      Beru~ light      blue            47   fema~
 8     97    32 NA         R5-D4 white, red red             NA   none 
 9    183    84 black      Bigg~ light      brown           24   male 
10    182    77 auburn, w~ Obi-~ fair       blue-gray       57   male 
# ... with 77 more rows, and 6 more variables: gender <chr>,
#   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
#   starships <list>

# last_col()
> starwars %>%
+   select(1:5) %>%
+   print() %>% 
+   relocate(name, .after = last_col())
# A tibble: 87 x 5
   name               height  mass hair_color    skin_color 
   <chr>               <int> <dbl> <chr>         <chr>      
 1 Luke Skywalker        172    77 blond         fair       
 2 C-3PO                 167    75 NA            gold       
 3 R2-D2                  96    32 NA            white, blue
 4 Darth Vader           202   136 none          white      
 5 Leia Organa           150    49 brown         light      
 6 Owen Lars             178   120 brown, grey   light      
 7 Beru Whitesun lars    165    75 brown         light      
 8 R5-D4                  97    32 NA            white, red 
 9 Biggs Darklighter     183    84 black         light      
10 Obi-Wan Kenobi        182    77 auburn, white fair       
# ... with 77 more rows
# A tibble: 87 x 5
   height  mass hair_color    skin_color  name              
    <int> <dbl> <chr>         <chr>       <chr>             
 1    172    77 blond         fair        Luke Skywalker    
 2    167    75 NA            gold        C-3PO             
 3     96    32 NA            white, blue R2-D2             
 4    202   136 none          white       Darth Vader       
 5    150    49 brown         light       Leia Organa       
 6    178   120 brown, grey   light       Owen Lars         
 7    165    75 brown         light       Beru Whitesun lars
 8     97    32 NA            white, red  R5-D4             
 9    183    84 black         light       Biggs Darklighter 
10    182    77 auburn, white fair        Obi-Wan Kenobi    
# ... with 77 more rows

実行環境と参考

環境

> sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

参考

to be continued...

俺的{dplyr1.0.0}メモ~summrise()編~

最近は、fgoキャメロットまで行きました。ガヴェインを令呪1画使って倒しましたとさ(ぎゃてぇ)

今更、{dplyr1.0.0}をキャッチアップしていこうのコーナー

今回はsummarise()、便利な集計関数ですね。

以下を参考にしました

{dplyr1.0.0}summarise()で追加された俺的に特徴的な機能はこちら

  1. ungroup()がいらなくなった
  2. 複数の値やデータフレームを関数を使えるようになった
  3. summarise()を使ってフォルダ内のデータを一気に読み込めるようになった

3は2の機能を使ったものなので、実質1,2ですかね

準備

今回はdplyr::stormsデータを使います。NOAA大西洋ハリケーンデータベースのベストトラックデータのサブセットです。このデータには、198個の熱帯性暴風雨の位置と属性が含まれており、暴風雨が観測されてから6時間ごとに測定されています。

#初めて、またはバージョンが古いときはインストールしなおしましょう
#install.packages("dplyr")
library(dplyr,warn.conflicts = F)
 
storms
# A tibble: 10,010 x 13
   name   year month   day  hour   lat  long status              category  wind pressure ts_diameter hu_diameter
   <chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <chr>               <ord>    <int>    <int>       <dbl>       <dbl>
 1 Amy    1975     6    27     0  27.5 -79   tropical depression -1          25     1013          NA          NA
 2 Amy    1975     6    27     6  28.5 -79   tropical depression -1          25     1013          NA          NA
 3 Amy    1975     6    27    12  29.5 -79   tropical depression -1          25     1013          NA          NA
 4 Amy    1975     6    27    18  30.5 -79   tropical depression -1          25     1013          NA          NA
 5 Amy    1975     6    28     0  31.5 -78.8 tropical depression -1          25     1012          NA          NA
 6 Amy    1975     6    28     6  32.4 -78.7 tropical depression -1          25     1012          NA          NA
 7 Amy    1975     6    28    12  33.3 -78   tropical depression -1          25     1011          NA          NA
 8 Amy    1975     6    28    18  34   -77   tropical depression -1          30     1006          NA          NA
 9 Amy    1975     6    29     0  34.4 -75.8 tropical storm      0           35     1004          NA          NA
10 Amy    1975     6    29     6  34   -74.8 tropical storm      0           40     1002          NA          NA
# ... with 10,000 more rows

ungroup()がいらなくなった

引数に.groupsが追加されました

  • drop_last: group_by()の最後の水準(変数)のグループ化を解除

    バージョン 1.0.0以前のデフォルト機能

    例:group_by(A,B)としたときに、Bのgroup化を解除する

  • drop: 全変数のgroup化の解除

    ungroup()と同じ

  • keep: group化の維持

  • rowwise: 各行をgruop化するのかな?

.groups = "drop"と指定することでungroup()する必要がなくなりました!


複数の値やデータフレームを関数を使えるようになった

summarise()range(),quantile()といった関数が使えるようになりました

> storms %>% 
+   group_by(name) %>% 
+   summarise(rng = range(wind))
`summarise()` regrouping output by 'name' (override with `.groups` argument)
# A tibble: 396 x 2
# Groups:   name [198]
   name       rng
   <chr>    <int>
 1 AL011993    25
 2 AL011993    30
 3 AL012000    25
 4 AL012000    25
 5 AL021992    25
 6 AL021992    30
 7 AL021994    15
 8 AL021994    30
 9 AL021999    25
10 AL021999    30
# ... with 386 more rows

quantile()も

> storms %>% 
+   group_by(name) %>% 
+   summarise(wind = quantile(wind, c(0.25, 0.5, 0.75)),
+           q = c(0.25, 0.5, 0.75))         
`summarise()` regrouping output by 'name' (override with `.groups` argument)
# A tibble: 594 x 3
# Groups:   name [198]
   name      wind     q
   <chr>    <dbl> <dbl>
 1 AL011993  25    0.25
 2 AL011993  27.5  0.5 
 3 AL011993  30    0.75
 4 AL012000  25    0.25
 5 AL012000  25    0.5 
 6 AL012000  25    0.75
 7 AL021992  30    0.25
 8 AL021992  30    0.5 
 9 AL021992  30    0.75
10 AL021994  21.2  0.25
# ... with 584 more rows

summarise()の返り値としてData-frame columnが可能になったからでしょうか。nest()unnest()をいい感じやってくれるような感じのイメージです。あまり、使うことはないかもですができることは広がりますね。


summarise()を使ってフォルダ内のデータを一気に読み込めるようになった

Data-frame columnをとれるようになったので、ファイルのパスを取得できればデータを{tidyverse}を使って一気に読み取れるようになりました。個人的にはこれがかなり便利かもと思っています。フォルダ内のデータを一気に読み込む方法はこれまでもいろんな方法が提案されていますが、これが一番速く読み取れるんじゃないでしょうか。どちらにせよ以下の様に三行で読み込めるのでいいと思います。

tibble(path = dir(pattern = "\\.csv$")) %>% 
  #rowwise()は各行をgroup化する
  rowwise(path) %>%
  #path列にある、各ファイルを読み込む
  summarise(read_csv(path))

summarise()単体で覚えておきたい機能はこんな感じでしょうか。これ以外にもacross()という関数を使えばさらに便利に使えるといわれたりしてますが。それは、また別の記事で言及しようかなと思います。今回はここまでー

Summarise each group to fewer rowsでも、「Data-frame columnは便利だけど、どのように使うのがベストかは一行の余地あるよね」的なことを書いているので、この辺りは変化する可能性があると思います。

to be continued...

R stanで得られた潜在クラスモデルの推定結果から、各個人がどのクラスに属するかを最頻値から求める

最近は、fgoを始めました(初めての星5は葛飾北斎)

今回は、ほぼメモです

stanで潜在クラスモデルを回した後に、個人がどのクラスに振り分けられることが多いかを知るときに困ったので

以下のようなデータがあったときに、個人がどのクラスに属してそうかを知りたいわけです。

いろいろ方法はあると思いますが、各行において最も大きなclの値を抽出して、個人ごとにその値の最頻値を獲れば、いけるかなと思います

データの構造としては

  • 参加者(id)が50人いる
  • 各参加者(id)につき10個のclデータがある(今回clは乱数)
  • clは4つある
> set.seed(42)
> id <- 1:50
> umr <- data.frame(id = rep(id,each=10),
+                   cl1 = rnorm(id,0,1),
+                   cl2 = rnorm(id,0,1),
+                   cl3 = rnorm(id,0,1),
+                   cl4 = rnorm(id,0,1)) 
> head(umr)
  id        cl1         cl2        cl3         cl4
1  1  1.3709584  0.32192527  1.2009654 -0.04069848
2  1 -0.5646982 -0.78383894  1.0447511 -1.55154482
3  1  0.3631284  1.57572752 -1.0032086  1.16716955
4  1  0.6328626  0.64289931  1.8484819 -0.27364570
5  1  0.4042683  0.08976065 -0.6667734 -0.46784532
6  1 -0.1061245  0.27655075  0.1055138 -1.23825233

プロセス

  1. 各行においてclで一番大きな数を抽出
  2. 個人ごとに一番大きなclの最頻値を抽出
# 個人の識別子
id <- 1:50
umr <- data.frame(id = rep(id, each = 10),
                 cl1 = rnorm(id, 0, 1),
                 cl2 = rnorm(id, 0, 1),
                 cl3 = rnorm(id, 0, 1),
                 cl4 = rnorm(id, 0, 1)) %>% 
  tibble::rowid_to_column(var = "subid")
 
umr %>% 
  # clをロング型に
  pivot_longer(cols = starts_with("cl"),
               names_to = "var_name",
               values_to = "value") %>% 
  # idとsubidでグループ処理
  group_by(id, subid) %>% 
  # subid内での最大値フィルタ
  filter(value == max(value)) %>% 
  # countを準備
  group_by(id, var_name) %>% 
  count(name = "cnt") %>% 
  # group切り直し
  group_by(id) %>% 
  # countがidごとに最頻値なレコードを残す
  filter(cnt == max(cnt))

# A tibble: 50 x 3
# Groups:   id [50]
      id var_name   cnt
   <int> <chr>    <int>
 1     1 cl3          4
 2     2 cl1          4
 3     3 cl3          4
 4     4 cl4          6
 5     5 cl1          4
 6     6 cl3          4
 7     7 cl1          4
 8     8 cl3          4
 9     9 cl4          6
10    10 cl1          4
# ... with 40 more rows

idが1の人は、cl3に振り分けられることが最も多いということがわかりました(以下省略)

{dplyr} 1.0.0を使えば、また違う書き方ができると思います。

braveで「はてなブログ」にログインする方法

braveは広告を排除することで高速かつ見やすいブラウザです(お金ももらえるよ)。要は、cookieをブロックすることでネット広告を出さないようにしてるのですが、それゆえに障害のあるページもあるわけです。

braveも初期では、「はてなブログ」もログインできないので、メモっておきます。

braveで「はてなブログ」ログインする方法

  1. braveを開いてsetting画面に行く
  2. setting画面の検索から🔎「cookie」と検索
  3. Site settings → Site Settings → Cookies and site data
  4. Allowの「Add」をクリックして「[*.]hatena.ne.jp」を記入して「Add」

OK!!

「Rで始める心理学Web実験」を書いた話

※本記事はR Advent Calendar 2019の4日目の記事です。


最近は「偽物語」を見てます。僕は忍派です。

10月のTokyo.Rに参加した際、30秒自己紹介タイムで

最近嬉しかったことは、卒論の心理実験をRで書けたことです。

と発言したら、後の懇親会にてある方々から声をかけていただきました。

「実験課題作成方法を教えてほしい」という旨でしたが、リップサービスだと思って話半分でLINEを交換したところ、後日、勉強会をやろうということになりました。

大変、恐縮でありますが若輩なりに資料を作成して、先日、無事勉強会を終えることが出来ました。結構(卒論に手を付けていない程度には)、資料を作成するのに時間をかけたので、ここだけで終わらせるのは勿体ないと思い、Rのアドカレに参加した次第であります。


発表資料

「Rで始める心理学Web実験」というテーマでの資料です。

サイトはこちら↓↓

Rで始める心理学Web実験 - Speaker Deck

結論から言うと、専修大学の国里先生が「jsPsychを用いた認知課題の作成」というチュートリアルをまとめてくださっていたので、その記事を参考にした次第です。

kunisatolab.github.io

2019年9月に、そろそろ卒業研究に着手しようと考えていた頃、自分の研究では

  • Web実験であること
  • 条件分岐する動的な質問紙

という二点が必要でした。qualtricsの様なサービスを使えば容易ですが、お金が無いマンなのでそれは出来ません。

pythonベースのPsychoPy3を使う方法もありました。

しかし、2019年9月時点でPsychoPy3でWeb実験をする際には、pavloviaというサービスを使うわけですが、pavloviaはビルダー(GUI)で作成した課題しか使えず、柔軟な実験プログラムを組むには少し不向きでした。

また、2019年10月よりサービス利用が有料化するという動きがあったためその辺も頭を悩ませていました。

pavlovia.org

そんなとき、某ベイズ塾で国里先生が、Rstudioを使う方法を紹介してくださり、無事に実験課題を作成し、データも取り終えられたということです。

資料内容

スライドの内容としては

  • どんな感じで実験課題を作成していくのか
  • 環境構築
  • 簡単な実験課題の作成(jsPsych)
  • firebaseを使ったweb実験のやりかた
  • Tipsや陥りやすいミス

という感じです。

個人的には、僕も初心者なのでRが良くわからない人目線で、なるべく丁寧に書いたつもりです。


個人的な所感

実験課題の作成には、jsPyschの勉強と自分の実験課題作成合わせて2週間程度で出来ました。実際には、2週間の間に他の作業のあったのでもっと短い時間で習得できると思います。

また、勉強会をした際も4時間程度で実験課題の作成まで行けましたので、特にR,Rstudio,Rmarkdownの知識が少しあればもっと早く学習可能ではないでしょうか。

実験課題作成には基本的には変数作成とプラグインの組み合わせなので、そこまで知識はいらないと思います。ただ、条件分岐や装飾としてhtml,css,jsの知識があるともっと可能性が広がると思います。


補足

RのjsPsychパッケージとしてjspsychrなるものがあります。これも見た感じ非常にいい感じのパッケージだと思いますが、2019年12月現在、バージョンが 0.0.0.9000 なので、バージョンアップがあったときに大きな変更がある可能性があるため、スライドでは紹介しておりません。

to be continued...

R rename_all() で変数名を一気に変更

最近は、 Dr.STONE みてます。そそりますね。

e-statからデータを取ってきたときとか、たまに変数名以外にもいらない要素がついてるときってありますよね。

そんな時、dplyr::rename_all正規表現が少しわかれば人生が豊かになるかもしれません。


環境

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

仮想データ生成

今回は変数名の先頭部分の要素を落とすという設定で書いてみます。

> library(tidyverse)

umr <- data.frame(
  V122_sakana = 1:10,
  V3455_oniku = runif(10),
  V974333_yasai = 1,
  V63_omame = sample(c(T,F),10,replace = T)
)

> umr
   V122_sakana V3455_oniku V974333_yasai V63_omame
1            1   0.6840092             1      TRUE
2            2   0.4950602             1     FALSE
3            3   0.8039745             1      TRUE
4            4   0.2525578             1     FALSE
5            5   0.8230913             1     FALSE
6            6   0.4581368             1     FALSE
7            7   0.9348471             1      TRUE
8            8   0.1878316             1     FALSE
9            9   0.6919405             1      TRUE
10          10   0.8223904             1      TRUE

変数名の先頭部分に「V122_」といったように、V+可変数値+「アンダーバー」がついて主要な変数名が後に続くとしますします。

主要な変数名だけを残して、それ以外のところは除外したいですよね。

そんな時は、dplyr::rename_allを使うといいんじゃないでしょうか。

dolyr::rename()

dplyr::renameは変数名を変更する関数です。

> rename(umr,sakana = V122_sakana)
   sakana V3455_oniku V974333_yasai V63_omame
1       1   0.6840092             1      TRUE
2       2   0.4950602             1     FALSE
3       3   0.8039745             1      TRUE
4       4   0.2525578             1     FALSE
5       5   0.8230913             1     FALSE
6       6   0.4581368             1     FALSE
7       7   0.9348471             1      TRUE
8       8   0.1878316             1     FALSE
9       9   0.6919405             1      TRUE
10     10   0.8223904             1      TRUE

一つ目の変数が「V122_sakanasakana」になってますね。これを一つ一つするのは大変なので、一気に変数名に適用してしまおうというのがdplyr::rename_all()

解答例

では先に、やり方を書きます

> umr %>% 
+   rename_all(vars(umr %>%
                    #変数名抽出
+                     colnames() %>%
                    #一文字目が文字or数字で可変から_まで除外
+                     str_remove(.,"^[:alnum:]+_")))
   sakana     oniku yasai omame
1       1 0.6840092     1  TRUE
2       2 0.4950602     1 FALSE
3       3 0.8039745     1  TRUE
4       4 0.2525578     1 FALSE
5       5 0.8230913     1 FALSE
6       6 0.4581368     1 FALSE
7       7 0.9348471     1  TRUE
8       8 0.1878316     1 FALSE
9       9 0.6919405     1  TRUE
10     10 0.8223904     1  TRUE

うぇい。

説明します。

dplyr::rename_all()

dplyr::rename_all()では、引数vars()を取る必要があります。今回は、中に関数処理を施していますが、変数名とオブジェクト(文字型ベクトル)の要素数が等しければオブジェクトで受け渡し可能です

> umr_name <- c("fish","meat","vege","bean")
> umr %>% rename_all(vars(umr_name))
   fish      meat vege  bean
1     1 0.6840092    1  TRUE
2     2 0.4950602    1 FALSE

上記と同じ処理を施したものをumr_nameというオブジェクトに入れて、rename_all(vars(umr_name))としても受け渡せます。

rename_all()の本質的な情報はここですね。最悪、このルールさえわかれば力業で何とかなります。ただ、変数名が628個あったらどうします?(A. 死にます)

そんなとき、正規表現を少し知っておくと生きることができます。

dplyr::rename_all()の中のお噺

> umr %>%
+     colnames() %>%
+     str_remove(.,"^[:alnum:]+_")
[1] "sakana" "oniku"  "yasai"  "omame"

vars()の中にピックアップしてみてみると、

  1. 変数名を変えたいデータセットオブジェクトを選択(今回はumr)
  2. colnames()で変数名をベクトルとして抽出
  3. stringr::str_remove()で選択した要素を除外

って感じですね。1,2はrename_all()を使うとしたら、大抵使うと思います。

では、最後の正規表現の説明を...

先っちょだけ正規表現のお噺

"^[:alnum:]+_"」こいつですね。

  • ^:先頭を示す。今回は一文字目から取り除くのが確定していたのでつけます。つけなくても動作しますが、予期せぬエラーを防げます。ちなみに、末尾の場合は$
  • [:alnum:]:アルファベットor数字なら何でもOK(記号はダメよ)
  • +:直前の要素が1回以上繰り返す場合に使います。今回は「V+可変数値」だったのでつけます。
  • _:変数名についてるアンダーバー

つまり「先頭が可変のアルファベットor数値でアンダーバーの部分までマッチする要素」を選択していることになります。

この辺の、正規表現の知識は私も詳しくないので、Rチートシート

stringrとregular expression(正規表現)が非常に参考になります。

そそりますね。こいつは

to be continued...