論理の流刑地

地獄の底を、爆笑しながら闊歩する

【備忘】交互作用のparameterizationについて考えたこと

実務上の分析をするなかでふと思ったこと。
いちいち書き留めるほどのことでもないけど、そうしないと忘れちゃうので(老兵はつらい)。

何が問題か

たとえば(例はほんとになんでもいいんだけど)
テストの成績Yに対しての通塾X=0,1の効果が、朝ごはんを食べているかどうかZ=0,1によって変わるか否か*1、ということに関して調べたいとする。

そういう関心をもっているとき、我々は往々にして以下のようなモデルを推定する。

Y_i = \alpha X_i + \beta Z_i+  \gamma X_i Z_i +\epsilon_i

このとき、\gammaの係数推定値(やその統計的有意性)をみて、前述の問いを検証することが多い。
\gammaが正でかつ有意であれば、「朝ごはんを食べることにより、通塾の成績への効果が増幅される」といった風に解釈するわけだ。

ここまではまぁとても基礎的なハナシなんだけど、実はひとつ盲点を抱えている
それはXの効果がZによって高められるのは、通塾勢の成績が上がることによってなのか、非通塾勢の成績が下がることによってなのかが、\gammaをみただけではわからない、ということである。
ことばを替えると、天井を引き上げる効果なのか床を押し下げる効果なのかがわからない

図で表すと、以下の場合両方ともXの効果がZ=1のときに20増えているけれど、上の図は天井への効果によるもので、下の図は床への効果によるものである。

f:id:ronri_rukeichi:20201027220535p:plain
交互作用1(天井への効果)
f:id:ronri_rukeichi:20201027220711p:plain
交互作用2(床への効果)

しかし上の定式化であると、一様に\gamma = 20ということになってしまう。
パラメタライゼーションを工夫することにより、天井への効果と床への効果を識別して推定する術を探す。

じゃあどうすればいいか

まずZが離散変数(というかダミー変数)の時を考えて、次に

Zが離散変数:ダミー変数のcodingを工夫

以下のような仮想データを生成する。天井への効果が+15,底への効果が-10で、\gamma=25となるような場合を想定している。

x0z0 <- rnorm(100 ,  30 ,10)
x0z1 <- rnorm( 100 , 20 ,10)
x1z0 <- rnorm(100 , 50 , 10)
x1z1 <- rnorm(100 , 65 , 10)
dta1 <- data.frame(  Y = c( x0z0, x0z1 , x1z0 , x1z1) , X = c(rep(0, 200), rep(1 , 200 )), Z = c(rep(0, 100), rep(1 , 100 ),rep(0, 100), rep(1 , 100 ) ))
lm_res <- lm( Y ~ X *Z , data = dta1)
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  30.0541     0.9335  32.194  < 2e-16 ***
#   X            20.0450     1.3202  15.183  < 2e-16 ***
#   Z            -9.5118     1.3202  -7.205 2.96e-12 ***
#   X:Z        24.3167     1.8671  13.024  < 2e-16 ***

\gamma = 24.3167である
われわれが目指すのは、\gamma_1\fallingdotseq15 , \gamma_2 \fallingdotseq -10でかつ \gamma_1 -  \gamma_2 =  \gammaという推定値を与えるものだ

Xが2値 × Zが2値=4通りがあるので、「Z=0のときのXの効果(\beta_0)」のほかに
「Zの天井への効果(X=1のときの効果)」「Zの天井への効果(X=0のときの効果)」という三つを識別すると、自由度を消費しきってしまう。
なので、上述の通常モデルで把握できる「Z」の主効果(というか実はX=0のときの条件付き効果なのだがはあきらめるよりない。

結局、求める推定値\beta_0 , \gamma_1, \gamma_2と切片\alpha使って、(X, Z)の四通りの状況を以下のようにあらわしたいと考えるのである。
E( Y|X=0, Z =0) = \alpha
E( Y|X=1, Z =0) = \alpha + \beta_0
E( Y|X=0, Z =1) = \alpha + \gamma_2
E( Y|X=1, Z =1) = \alpha +  \beta_0+ \gamma_1

よって、 (X=0, Z =1)のときは1,それ以外は0をとる変数U_1と、 (X=1, Z =1)のときは1,それ以外の場合は0をとる変数U_2を作って、Xとともに回帰式に投入すればよいです。

dta1 <- dplyr::mutate( dta1 , gamma1 = if_else( X==1 & Z==1 , 1 ,0 ), gamma2 = if_else( X==0 & Z==1 , 1,0))
lm_res2 <- lm( Y ~ X+ gamma1 + gamma2 , data = dta1)

# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  30.0541     0.9335  32.194  < 2e-16 ***
#   X            20.0450     1.3202  15.183  < 2e-16 ***
#   gamma1       14.8049     1.3202  11.214  < 2e-16 ***
#   gamma2       -9.5118     1.3202  -7.205 2.96e-12 ***

14.8049  -(-9.5118 )= 24.3167なので、想定通り\gamma_1 - \gamma_2 = \gammaになっている。

Zが連続変数の場合はどう考えるとよいか

状況をより複雑にしよう。さて、Zが連続変数の場合はどう考えるとよいだろうか。
さっきはZを朝ごはんをたべていたかどうか(二値変数)、としたが、今回は平日の自宅での平均勉強時間、にしよう(表記はHとする)。

通塾の効果は+20点だが、相乗効果があって勉強時間が1hあがるごとに、効果は+3.5点されると考える。
勉強自体の効果は1hごとに+3.0点としよう。

仮想データをつくり、そして通常の交互作用モデルのOLSを推定する。

cpt <- rnorm( 200, mean=30 , sd = 0.5) #切片(30点とする)
juku <-  sapply( rnorm( 200) , function( x){
  return( if_else( x >=0 , 1 , 0))
}) #塾にいくか否か
std_h <-  rpois( 200 , 3)/2 #勉強時間(Poisson分布に従うと仮定)
errT <- rnorm( 200, sd= 2.5)#誤差項 
#点数
tensu <- cpt + 20 * juku +  3 *std_h +  3.5 * juku * std_h + errT
test_dta <- data.frame( tensu , juku , std_h )

test_res1 <- lm( tensu ~ juku + std_h + std_h:juku, data = test_dta)
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  31.0673     0.5222  59.491  < 2e-16 ***
#   juku         19.3359     0.7091  27.270  < 2e-16 ***
#   std_h         2.7675     0.3151   8.784 7.90e-16 ***
#   juku:std_h    3.4089     0.4247   8.026 9.07e-14 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

\gamma =  3.4089である
さて、ここで第三の変数が連続変数である本ケースにおいても、「天井」と「床」それぞれの効果を識別したうえで推定値が得たい。

さきほどの離散変数の場合と同じく、
Y_i = \beta_0 X_i + \gamma_1 X_i H_i + \gamma_2 (1-X_i) H_i
といった形での定式化をする。

test_dta <- dplyr::mutate( test_dta , U0 = std_h *juku , U1 = std_h*(1- juku))
test_res2 <- lm( tensu ~ juku + U0 + U1, data = test_dta)

# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  31.0673     0.5222  59.491  < 2e-16 ***
#   juku         19.3359     0.7091  27.270  < 2e-16 ***
#   U0            6.1765     0.2848  21.685  < 2e-16 ***
#   U1            2.7675     0.3151   8.784  7.9e-16 ***

\gamma_1 -\gamma_0 =  6.1765 -2.7675  = 3.4089 が成立している。
勉強時間1h増加につき、天井(通塾してる子)の得点は6.2点のび、底(非通塾勢)2.8点のびていることがわかる。

Conclusion

特に難しいことをみようとしているわけではないが、意外と普段のモデルの定式化だと見えていないものがあることに気づかない、という話でした。


あっかんべーだ / 瀬名航 (covered by 霧島伸行)


Enjoy!!

*1:どんな問いやねん

グラぽで掲載されていた出場時間集中係数の計算プログラム&分析結果再現をこころみる

仕事に疲れたので、休憩用に遊んでいたものをそのまま記事にするだけのやつ*1
車輪の再発明ってやつですね。

Motivation:コナカ先生の知的な試み

先日、グラぽでとても興味深い記事が出ていた。
grapo.net

出場時間がどれだけ特定の選手に集中しているのかを指標化していて、非常に面白い試みであった。

元ネタはジニ係数とのことであるので、結果の再現や計算プログラムは簡単に作れるのではないかと考えた。


また、以下のTwitterのように、計算方法が分からないとつぶやいているかたもおられた。

そこで、小中先生は(本人もおっしゃっている通り紙幅の都合上で)計算方法は直接語っていないので、計算プログラムや計算結果を再現できないか考えてみる。

ジニ係数(オリジナル)の計算方法について

ジニ係数の計算の方法は色々あるが、スマート*2なのは、Jacques Silberが紹介しているG-Matrixを用いた方法である。
G-Matirxとは、対角成分が0, 非対角成分の右上部分が1, 非対角成分の左下部分が-1の正方行列である。

ジニ係数を計算するには、所得シェアを昇順にならべたベクトルをI 、対応する人数シェアをPとし、G-MatrixをGとすると
Gini= P' GI
という計算を行えばよい('は転置)。

それ自体は、以下のような簡単なプログラムで実装できる。

#Gini係数の計算
gini_calc <- function(x){
  x_s <- sort(x) 
  p <-  rep( 1, length(x_s)) / length(x_s)#人数シェア
  incm <- x_s / sum(x_s) #所得シェア
  g_mat <- make_gmat(length(x)) #G-Matrix
  g_index <- t(matrix( p)) %*% g_mat %*% (matrix( incm))
  return(g_index)
} #
#G-matrixの作成関数
make_gmat <- function(n){
  n_mat <- matrix( 0 , nrow=n , ncol=n)
  n_mat[col(n_mat) > row(n_mat)] <- 1
  n_mat[col(n_mat) < row(n_mat)] <- -1
  return(n_mat)
} #G-Matrix

集中係数へのアレンジ

さて、件の指標はもちろんジニ係数に着想を得ただけであって、そのものではない。
なぜなら不平等指標としてのジニ係数の完全不平等状態は「1人」が金や資産を独占した場合を想定しているが
サッカーにおける独占は「11人」がフル出場している場合だからだ。

そのため、1人が独占している場合でなく11人が独占している場合を完全不平等状態=1としなければならない

よって、出場時間を元にジニ係数の計算式を適用して出た数値を、11人が独占した場合のジニ係数で除するという補正を行っているのではないかと推測した。

そのようなアレンジ(推定)を加えて、以下のような関数をつくる

gini_konaka <- function(x){
  #30人に足りない場合は30人にする 
  add_n <- 30- length( x)
  x <- c( x , rep( 0, add_n))
  #11人による独占状態のベクトルをつくる
  total_min <- sum(x)
  dokusen_min <- c( rep(0,19), rep(total_min/11 , 11))
  
  #独占状態のGiniを計算する
  dokusen_gini <- gini_calc(dokusen_min)
  club_gini <- gini_calc(x)
  return(round( as.numeric( club_gini/dokusen_gini) ,3))
}

再現できたかの検証

さて、グラポに掲載されていた出場時間集中係数は、名古屋=0.891, 川崎=0.763であった。

手元のデータ(取得方法は過去記事参照)から、10/15時点までの名古屋や川崎のデータをとってくる。

f:id:ronri_rukeichi:20201022225210p:plain
川崎の出場時間データ(10/15時点)

試した結果は以下の通りである。無事グラぽの記事の値と一致している。

f:id:ronri_rukeichi:20201022225358p:plain
自作関数による集中係数の計算結果

アレンジに関する推定も正しかったようだ。
小中係数は、出場時間を元にジニ係数の計算式を適用して出た数値を、11人が独占した場合のジニ係数で除することによって求められていたようだ。

めでたし、めでたし。

Enjoy!!

追記:Excelでやるには

この記事を書いたあと、上で引用させていただいたTWUさんから直接Feedbackをいただいた。

冷静に考えて、いや確かにこれ読んでも分かる人にしか分からんし不親切かもしれないと思ったのでExcelでやる方法を追記する。
TWUさんの例(おそらく8人制の少年サッカー?)をもとに解説する。

24分×8人×3試合=576分をどう15人に分配するか?という状況を考えよう。
仮想的に、15人の出場時間が、降順に72分×3、58分、55分、52分、51分、32分、31分、20分、15分、14分、13分、10分、9分である状況を考える

Step1:実際の出場時間(降順)と極限集中状態の出場時間を横に並べて記載する

f:id:ronri_rukeichi:20201023060005p:plain
実出場時間と、完全独占の場合(仮想)の出場時間

まず上の画像のように、実際の出場時間(一行目)と完全独占を仮定した場合の出場時間(二行目)をExcelに書き込んでいく
ここではExcelのF列を左端とし、2行目を見出し、3行目を実際の出場時間、4行目を独占状態の場合の出場時間としている。

Step2:G-Matrixを作成する。

f:id:ronri_rukeichi:20201023061913p:plain
G-Matrix@Excel

次に、上の画像のようにG-Matrixを作成する。

このとき、行列の一番左上がちょうどシートの対角線上にくるようにする(上の画像だとF列6行目)と、
以下のような数式を書いてコピペするだけでG-Matrixができる

f:id:ronri_rukeichi:20201023060907p:plain

Step3:行列計算をおこなって集中係数を算出。

実はExcelではmmult()関数を使って行列計算ができる(自分も今初めて知りました)。
これを使うと、非常に簡単に行列計算ができるので、それを使って小中係数を計算する

まず、さきほど書いたG-Matrixの下側で同じ列数だけセルを選択する

f:id:ronri_rukeichi:20201023063129p:plain
(画像の下側、緑のワクのところ)


次に、数式バーに=MMULT( 実際の選手の出場時間を書いた範囲、G-Matrixの範囲)を入力する

f:id:ronri_rukeichi:20201023063352p:plain
入力する数式
f:id:ronri_rukeichi:20201023065139p:plain
数式詳細

数式が入力できたら、(Enterではなく)Ctrlキー+Shiftキー+Enterキーを押下すると、選択した範囲に行列の計算結果が反映されている。

f:id:ronri_rukeichi:20201023063737p:plain

そして、完全独占状態についても同じように行列計算を行う(Ctrl + Shift + Enterキーの同時押しを忘れずに)

f:id:ronri_rukeichi:20201023063928p:plain

すると、以下のように完全状態についての行列計算の結果が出力される
f:id:ronri_rukeichi:20201023064132p:plain

あとは、(一行目の和)/(二行目の和)を計算すると小中係数になっている*3

f:id:ronri_rukeichi:20201023064624p:plain

f:id:ronri_rukeichi:20201023064658p:plain


今回の仮想ケースの場合は、0.738095が小中係数であるようだ。川崎フロンターレよりもさらに出場機会がバラけている。


以上、Excelでも小中係数を計算する方法を示した。

Enjoy!

*1:こういう遊びに走るから仕事ができない

*2:というか他の分離指標などとも結び付きがあって美しい

*3:実は上述のGini-Indexの計算方法から割愛してもよい部分を抜いているが、そこまでは説明しない

【R備忘】既存のテキストファイルの文字コードをUTF8からShift-JISに変換する関数

備忘。たぶんすごい基礎的なことだと思うんだけど、あんまり日本語情報がない気がしたんで。
まぁテキストエディタ文字コード変更して保存しなおせばええやん、って話ではあるんだが、Rだけでなるべく処理は完結したい。

以下みたいな関数をつくっておくとよい。

utf_to_sjis <- function( dir_chr , file_chr ){
  old_path <- paste0( dir_chr , file_chr) #旧ファイル名を取得する
  new_path <- gsub( ".txt","_s.txt",old_path)  #新しいファイル名をつくる
  cont_utf <- readLines( old_path, encoding = "UTF8") #UTF-8でファイル内容を
  f_link <- file( new_path, encoding="SJIS") #ファイルリンクをSJISで作成
  write( cont_utf , file=f_link) #新ファイルに内容を書き込む
  close(f_link) #connectionをclose
  file.remove( old_path) #旧ファイルをdelete
  file.rename( from = new_path , to = old_path) #新ファイルを旧ファイル名に替える
} #dir

#実行例
utf_to_sjis( "./Document/",  "test1.txt")

参考にしたURLは以下。

  1. R: 文字コードを指定してのファイル出力 - Research and Study Notes of Rabbi Zakki
  2. r - RMeCabでUTF-8で書かれたファイルを扱う方法 - スタック・オーバーフロー



yonige -アボカド-【Official Video】

Enjoy!!

【分析準備】Jリーグの会見コメントからコーパスをつくる by rvest/RSelenium/RMeCab

Motivation

最近、マッシモの会見が面白くなってきたともっぱらの評判である。
記者への警戒心がやや解けてきたのか、わりとざっくばらんに色々と話してくれるようになっている。


前任者である風間監督も、なかなか言葉選びに独特のセンスやこだわりを感じさせる監督して定評があった。
そこで、「言葉」という視角から、マッシモさんと風間さんの特徴比較をできたら良いな、と思う。

しかし、その二人のデータを単純に比べるだけではいまいちJリーグの監督全体の言説空間における位置づけがわかりにくいし、なにより面白くない。
したがって、2019-2020のJ1の監督全員の会見コメントをとってきたうえで、LDAなどにかけていきたい。

<参考URL>

実装(データ取得)

具体的には、以下のような工程を踏んで処理することとなる。

  1. 取得対象のURLの取得
  2. URLから原テキストの取得
  3. RMeCab::docDFにかけやすいように、監督×試合ごとの単位でテキストファイルに書き出す

①:取得対象のURLの取得

J1の日程・結果検索のところで、カテゴリを「J1」かつ対象月を1~12月すべてにして、得られたURLから、「監督コメント」のリンクを取得していく。
あとは適当にHadley神のつくりたもうた、SelectorGadgetをつかいつつ得られたCSSセレクタによりアクセスしていく。

一応関数化する(年度やJ1 or J2を指定できるように)

getMatchDay_All <- function( cat ="j1",  year = 2020){
  md_url <- paste0("https://www.jleague.jp/match/search/?category%5B%5D=" ,cat,"&year=",year,"&month=01&month%5B%5D=02&month%5B%5D=03&month%5B%5D=04&month%5B%5D=05&month%5B%5D=06&month%5B%5D=07&month%5B%5D=08&month%5B%5D=09&month%5B%5D=10&month%5B%5D=11&month%5B%5D=12&section=")
  
  md_page <- read_html(md_url)#読み込み
  a_nodes <- html_nodes( md_page , css=".matchTable li:nth-child(3) a")
  
  href_arr <- sapply( a_nodes , function( a_node){
    return( html_attr(a_node , name="href"))
  }) #
  return( href_arr)
} #function 

#実行
href_j1m <- getMatchDay_All(cat = "j1", year = 2020)

#結果
 head(href_j1m)
# [1] "/match/j1/2020/022101/coach/" "/match/j1/2020/022201/coach/"
# [3] "/match/j1/2020/022202/coach/" "/match/j1/2020/022203/coach/"
# [5] "/match/j1/2020/022204/coach/" "/match/j1/2020/022301/coach/"

②:監督コメントの取得

こちらは動的生成されるコンテンツ*1であるので、rvestは力及ばずでRSeleniumの出番となる。

◆対象ノードの情報取得
まずは、cssセレクタを利用し、基礎的な情報をとってくる(list形式、1個目Home , 2個目Away)

get_ComeContents <- function(clnt,url){
  clnt$navigate( url)
  Sys.sleep(0.4)
  
  m_doc <- (clnt$getPageSource())[[1]]
  m_page <- read_html(m_doc)
  come_info <- html_nodes( m_page , ".commentArea .commentBox")
  return(come_info)
} #get_Mng

#実行
come_NV <- get_ComeContents(clnt , "https://www.jleague.jp/match/j1/2020/022201/coach/") #第一節 仙台vs名古屋戦


◆ノードからの情報の構造化
さて、最終的な目的はコーパスをつくることなので、それがしやすいようにデータ化する。
監督会見のコメントのテキストには、記者からの質問も含まれるので、それを除く*2
また、あとで整理がしやすいようにクラブ名と監督名も取得する

comeToData <- function(c_node){
  #内部関数  
  gsub_arr <- function(chr , frm_str, to_str){
      for( i in seq_along( frm_str)){
        chr <- gsub( frm_str[i], to_str[i],chr)
      }
      return(chr)
    } #gsub
    
  h4_node <-  read_html(as.character(html_node( c_node, "h4")))
  p_node <-  read_html(as.character(html_node( c_node, "p")))
  club_name <- html_text(html_node( h4_node , "span"))
  comment_txt_arr <-  sapply(html_nodes(p_node,xpath="//p/text()"), function(x){
     return( html_text(x))
   }) ##
  comment_txt <- paste0( comment_txt_arr, collapse="\n")
  manager_name <- trimws( gsub_arr(  html_text( html_node( p_node , css="span")) , c("\\]","\\[", "監督") , rep("",3)))
   
  return( list(Club= club_name , Manager =manager_name, Comment = comment_txt))
} #func
##実行
xtN <- comeToData( c_node = come_NV[[2]])
txtN$Club
#[1] "名古屋グランパス"
txtN$Manager
#[1] "マッシモ フィッカデンティ"
txtN$Comment
#[1] "\nわれわれのやりたいサッカー、つまり相手の陣地...

あとはファイルに書き込むだけある。
ただし、後述するようにdocDF()はファイル名を列名とするので、監督名、クラブ名、日付あたりの分類に使いやすい情報をファイル名にいれておくのがよい。

#関数、第一引数はget_ComeContents()の戻り値
comeToFile <- function( come_info , url ,dir ="C:/Users/XXX/Documents/01_data/Manager_txt/"){
  ## 日付をとってくる ##
  gsub_arr <- function(chr , frm_str, to_str){
    for( i in seq_along( frm_str)){
      chr <- gsub( frm_str[i], to_str[i],chr)
    }
    return(chr)
  } #gsub_arr
  date_v <-  substr( gsub_arr(url,c(".*/2020/","/coach/") ,rep("",2) )  , 1,4)
  c_long <- c('湘南ベルマーレ','浦和レッズ','ベガルタ仙台','名古屋グランパス','柏レイソル','北海道コンサドーレ札幌','川崎フロンターレ','サガン鳥栖','セレッソ大阪','大分トリニータ','清水エスパルス','FC東京','横浜F・マリノス','ガンバ大阪','サンフレッチェ広島','鹿島アントラーズ','ヴィッセル神戸','横浜FC')
  c_alph <- c("shon","uraw","send", "nago", "kasw","sapp","ka-f","tosu","c-os","oita","shim","fctk","y-fm","g-os","hiro","kasm","kobe","y-fc")
  club_df <- data.frame( Name = c_long, ALPH = c_alph)
    lapply(list( comeToData(come_info[[1]]) , comeToData(come_info[[2]])) ,function(club_info){
    mng_name <- gsub(" ", "",club_info$Manager) #監督名
    club_name <- as.character((dplyr::left_join( data.frame( Name = club_info$Club) ,  club_df , by ="Name"))[["ALPH"]] )
    fname <- paste0(c(date_v, club_name ,mng_name,".txt"), collapse="_")
    full_path <- paste0(dir , fname) #フルパスを指定
    come_txt <- gsub("\n" ,"" , club_info$Comment) #改行だけ消しとく
    write(come_txt , file=full_path)
  }) #lapply  
}#func

#実行する
comeToFile(come_NV , "https://www.jleague.jp/match/j1/2020/022201/coach/")

実行すると、以下のような感じでtxtファイルができあがっている。

あとは、先ほど取得したURL一覧から適当にループで回していけばよい。

実装(取得したtextからのコーパス化)

テキストを監督ごとに結合する

20/10/14時点までの204試合(J1)の監督会見のデータを、上の処理で取得できたわけである。
代行*3を除外すると、19名(神戸がフィンク、三浦の2名)の監督がいるので、監督ごとに結合を行う。

daikou_idx <- grep("コーチ" , mng_names)
kantoku_names <- mng_names[setdiff( 1:length(mng_names) , daikou_idx)]

##--- 関数:テキストファイル結合用 ---##
txtFileMerge <- function( fpath , new_path){
  file_tbl_l <-sapply( fpath, function(fp){
    rt <- read.table(fp)
    rt_chr <-paste0(lapply( rt[1,] ,  as.character),collapse="")
    return( rt_chr)
  })# lapply
  merge_txt <- paste0( file_tbl_l, sep="\n")
  write( merge_txt ,file=new_path)
}# f_name

old_fp_l <- list()
new_fp_l <- list()
for( i in seq_along( kantoku_names)){
  mng <- kantoku_names[i] #監督名
  mng_idx <- grep( mng, list.files(mng_dir))
  tgt_files <- list.files(mng_dir)[ mng_idx]
  tgt_fp <- paste0(mng_dir , "00_merged/",  mng, "_all.txt")  
  txtFileMerge( fpath = paste0(mng_dir, tgt_files), tgt_fp)
  old_fp_l[[i]] <-   paste0(mng_dir, tgt_files)
  new_fp_l[[i]] <-  tgt_fp
} #for

結合できた。
監督名ごとに結合された会見のテキスト

RMeCab::docDF()でコーパスをつくる

さて、監督ごとの結合ファイルもつくることができたのでコーパスをつくる
RMeCabのdocDFで読み込んでいく

library(RMeCab)
mng_j1_df20 <- RMeCab::docDF(target = paste0(mng_dir, "00_merged"),pos=c("名詞","形容詞","動詞","副詞"), type=1)

# > head( mng_j1_df20[,1:7],3)
# TERM POS1     POS2 アンジェポステコグルー_all.txt ザーゴ_all.txt
# 1    , 名詞 サ変接続                              0              0
# 2    0 名詞       数                              2              2
# 3   00 名詞       数                              0              0
# トルステンフィンク_all.txt ネルシーニョ_all.txt
# 1                          0                    0
# 2                          2                    0
# 3                          1                    0

うまく読み込めている。

次回予告

これで準備は完了した、次回記事ではこのコーパスをもとにして、参考URL①を参考にトピックモデルを適用していく。


Spontania feat. AZU - 同じ空みつめてるあなたに

Enjoy!!

*1:最近はこういうサイトが多いので、むしろrvestでゴリゴリいけるFootball LABのほうがめずらしいのかも

*2:厳密に考えると、自分から話す内容と質問への応答では意味合いはかわってくるが、ここでそれは区別しない

*3:四方田@札幌、ビベス@神戸

【R小ネタ】最終引数以外でも、複数の変数を非標準評価(NSE)形式で与えたい場合のやりかた。

なんか簡単そうで意外と思いつかなかったので簡易な備忘

Probrem

関数内でdplyrの諸関数をつかうときに、NSEな感じでquotationをつけずに引数を複数指定したいときがある。
そんなとき、最終の引数であれば 「...」をquos()で受け取ってから、!!!演算子をつかえばよい(下の記事参照)

ronri-rukeichi.hatenablog.com

でもそうじゃなくて、最終以外の引数で

func1( data , list( X, Y , Z) , c1, c2 ,c3)

みたいな感じでの指定がしたいときどうするのかっていう話です。

Solution

おなじことをやりたい人が海外にいたので、StackOverFlowにやりかたがあった。
(というかたいていこういう細かいこと調べたい時英語情報しか出てこない*1
stackoverflow.com

以下のようにするとよい。
例として第二引数で欲しい変数の種類を、第三引数でグルーピングに使う変数を指定する、みたいな関数を考える。

test_fun <- function( dta , vname , ...){
grp_v <- quos(...)

slct_v <-  (as.list( rlang::enexpr(vname)))
if( length(slct_v) > 1){slct_v[[1]] <- NULL }

return( dplyr::select(dplyr::group_by(dta , !!!grp_v) , !!!slct_v))
} #func

#つかってるのはpanelrパッケージのWageData
wg1 <- test_fun(WageData , list(wks , fem , lwage), union , occ )

# wg1
# # A tibble: 4,165 x 5
# # Groups:   union, occ [4]
# union   occ   wks   fem lwage
# <dbl> <dbl> <dbl> <dbl> <dbl>
#   1     0     0    32     0  5.56
# 2     0     0    43     0  5.72

途中listの要素数が1を超えてる場合に、1<- NULLしてるのは、要素が2つ以上ある場合ひとつめに余分な要素*2が入ってしまうからである。
とりあえず、as.list(rlange::enexpr() )でexpressionのlist化ができることがわかった。

めでたし、めでたし。


【MV】バックドロップシンデレラ/さらば青春のパンク【11/2ベストアルバム発売】


Enjoy!!

*1:もっと日本のR界隈頑張れって思ったけど頑張ってる人は英語で発信してるんだろう

*2:list()ってやつ

選手の個人スタッツのデータベースをつくる by rvest&RSelenium

鉄は熱いうちに打て、ということで分析基盤をやる気のあるうちにつくる。
テーマや領域にかぎらず、大抵の分析は「分析可能なデータをつくる」工程に大部分の時間・労力が費やされる説あります。

Introduction

Jリーグをクラブ単位でデータ化することはあらかたできた。
ronri-rukeichi.hatenablog.com


ので、次にすることは、個人単位でのStatsをデータ化することである。

例によって、rvestによるデータ化をこころみていく。。。。。予定だったが、
対象のサイトがコンテンツを(Javascriptなどにより)動的生成している系のやつだったので、RSeleniumを使う必要がでてきた。
その学習記録もふくめた、備忘。

最終目標は、sofascoreの試合詳細ページ(例:名古屋vs清水)から個人の各種Statsをとってくることである。

◆参考URL

  1. vignettes
  2. reference
  3. [R]RSeleniumの環境構築(Windows) - Qiita
  4. [R][RSelenium/rvest]スクレイピングから某夢の国のクチコミを可視化してみた - Qiita ⇐良記事
  5. [翻訳] RSelenium vignette: RSeleniumの基本 - Qiita ⇐最重要

目標の設定

「何を実現したいか」、を先に決めることが学習の効率性を高めるためには肝要である。
ここでは、試合詳細のページ(例:ガンバvs名古屋 20/9/23)から出場選手の各統計指標を取ってくることを目標としよう。

なぜRSeleniumが必要かといえば、まずページを読み込んだ時の初期表示においては、以下のようにフォメ図が表示されているが、

f:id:ronri_rukeichi:20201006134121p:plain
出典:Sofascore

ここで、「Player Statistics」をクリックしないと、指標に関するtableが生成されないからである。
また、テーブルは「Summary」「Attack」など6種類に分かれているので、これも見出しをクリックしていかないと各種統計指標を取得できない。

したがって、以下のような手順でのScrapingのプログラムをかく必要がある。

  1. 試合詳細のURLにアクセス
  2. 「Player Statistics」をclick
  3. 各テーブルの見出しをclick
  4. 表示されたテーブルの内容をrvestで取得&データ化
  5. 3~4を繰り返す

下線部が、ブラウザの自動制御が必要なところとなる。ここにおいて、RSeleniumが必要となってくる。

実装(Rseleniumのとこだけ)

Seleniumサーバーは立ち上がっているとした前提でね。

初期化

library(rvest)
library(Rselenium)

clnt <- RSelenium::remoteDriver(port=4444L,browserName = "chrome") #Chromeで立ち上げる
clnt$open() #接続

これでChromeが立ち上がる。

クリックして表示をかえていく

##--ガンバvs名古屋の試合に移動--##
url1 <-  "https://www.sofascore.com/gamba-osaka-nagoya-grampus-eight/LmbsNmb"
clnt$navigate(url)

##--clickして「Player Statisticsを表示させる」--##
pstat_tag <-  clnt$findElements(using="css" , ".u-mB12 a+a.bfqsCw")[[1]]
pstat_tag$clickElement() #Clickして「Player Statistics」に移動
Sys.sleep(0.5) #0.5秒待機

##--6種類ある見出しを順番にclickしていく--##

click_tgts <- clnt$findElements(using="css" , "div.dtirhQ a.bfqsCw") #見出しのタグをまとめて取得
  for( i in seq_along( click_tgts)){
    click_tgt <- click_tgts[[i]]
    click_tgt$clickElement()
    Sys.sleep(0.3) #待機
 #---ここにrvestでの取得処理を入れる---#
  } #for

一応読み込みの時間を考慮して、0.3秒ずつの待機時間を挟んでいる。
例によってcssセレクタはselectorGadget(あるいはChromeの「検証」機能)で確認しとります。

実装(rvestで実際にテーブルを取得していく)

クリック等でお目当ての状態に遷移させたうえで、ページの情報を取得していくには、
getPageSource()でドキュメントを取得してから、rvestに渡すのが一番ラクだ。

今回の場合、上のクリックで表示されるテーブルを変えていくfooループのところに、以下の関数で書かれているような、tableをデータフレーム化して取得する処理を書いていくのがよい。

#お目あてのテーブルを表示させたうえで、remoteServerのオブジェクトをわたす
pTbl_get <- function(clnt){
  pdoc <- clnt$getPageSource()
  phtml <- read_html( pdoc[[1]])

  ### 対象のテーブルタグを探索sur
  tbl_node <- html_nodes( phtml , css ="table.dygnyR")
  tbl <- html_table( tbl_node[[1]], fill=NA, header=T)
  return(tbl)
} #function

すると、以下のようにrvestだけでは取得できないデータが取得している

> head(ptbl1$Attack,2)
#               + Shots on target Shots off target Shots blocked
#1 NA         Leandro               3                1             0
#2 NA Akihiro Hayashi               0                0             0
#  Dribble attempts (succ.) Notes Position Rating
#1                    2 (2)     -        F    7.9
#2                    0 (0)     -        G    7.4

あとは、分析しやすいように各指標を整形・変換したうえでdata frame化すればよい。

Conclusion

上の基本挙動に加え、データの整形処理*1等も加えた関数をつくったうえで実行してみる。

stat_df1 <-  get_PStat(clnt , url1)
head( stat_df1, 2)


# PlayerName Goals Assists Tackles Passes_All Passes_Suc Duels_All
# 1 Yosuke Ideguchi     0       1       2         70         65         3
# 2   Takashi Usami     1       0       0         14         13         3
# Duels_Suc gDuels_All gDuels_Suc aDuels_All aDuels_Suc Position
# 1         2          2          2          1          0        M
# 2         1          3          1          0          0        M
# Shots_onTgt Shots_offTgt Shots_blocked Dribbles_All Dribbles_Suc
# 1           1            1             0            0            0
# 2           2            0             0            0            0
# Clearances Interceptions Dribbles_Past Block_Shots Touches KeyPasses
# 1          0             3             0           0      84         3
# 2          0             0             0           0      21         0
# Crosses_All Crosses_Suc Longballs_All Longballs_Suc BoalLost
# 1           1           1             3             0        9
# 2           2           0             0             0        5
# Fouls_For Fouls_Against
# 1         1             0
# 2         1             1

出場した個人のStatsがうまく取れている。

まだ、

  1. 各ページのURLの取得・走査ルールをつくる
  2. 差分取得プログラムをつくる
  3. Football Labの方のデータと統合できるようにKey tableをつくる*2

などやるべきことはあるが、ひとまず目標は達成したといってよいだろう。



ACIDMAN - 「Rebirth」MV

Enjoy!!

*1:実際はここがめんどくさくて意外と工数かかったりする。が、丁寧にやらないとあとで分析しようと思った時に泣きをみる

*2:上のプログラムだとアルファベットの選手名しか分からないので、背番号や所属クラブの名前と照合できるようにする

【データ検証】マッシモ名古屋は本当に「先行逃げ切り特化型」なのか?それはなぜか?

データから見るマッシモ名古屋シリーズ*1

◆Outline

Introduction: グラぽにおける問題提起

マッシモ名古屋=「先行逃げ切り特化型」説の提起

清水戦(3-1)のレビュー記事@グラぽにおいて、現在のグランパスは「先行逃げ切り特化型」である、という説が提起された。
grapo.net

この記事では、リードされた状態でグランパスがHTに突入したときには、100%負けるというデータがラグ氏から提示され、多くの名古屋サポが*2衝撃をうけた。

引かれた相手を崩すことができないマッシモ名古屋は、やはり「先行逃げ切り特化型」なのであろうか。

「先行逃げ切り特化型」説の検証の必要性

ラグ氏の提起された説は、確かに、われわれの実感と見事に符合する。
今のグランパスは、引かれた相手を配置の妙やパスワークで崩す仕組みも、空中戦で押し切るパワープレー力も備えていない、と皆が感じているのではないのだろうか。

ただしここで冷静になって考えてみると、「先行逃げ切りは果たしてマッシモ名古屋だけの特徴なのか?」という疑問も浮かびあがってくる。
サッカーはそもそも点が入りにくいスポーツだし、相手が出てこないと苦労するのは名古屋だけではないだろう。

少数であるが、下のように「先行逃げ切り特化型」仮説の根拠となるデータの解釈に慎重な声も上がっていた。


確かに厳密に考えれば、「先行逃げ切り特化型」であることが、マッシモ名古屋だけの特徴なのか否かは、「他チームがHT時点でビハインドの場合の勝ち点取得率」と比較したうえで判断される必要があるだろう。


そこで本稿においては、Football Lab*3のデータを用いたうえで、
マッシモ名古屋が本当に「先行特化逃げ切り型」であるかを判断したい。

分析:マッシモ名古屋は「他チームと比較しても」先行逃げ切り特化型なのか?

上記記事で、ラグ氏は以下のように述べておられる。

お分かりでしょうか。
今年、これまでのグランパス、ハーフタイム時点でリードしてたらほぼ勝つ。同点でもまあ勝つ。リードされてたら負けます。


支配率を見ても、リードしてたら自分達はボールを持たなくなる。一方、リードされてたらボールを持って攻撃しようとするんだけど、そこまでの反発力を持たずそのまま負ける、とデータからわかりますね(分母は少ないですが………)

最終的な勝敗だけで判断するにはケース数が少ないかも..という留保もラグ氏はつけていることからも、いくつかの指標からの検討を試みる。
リードされてたらボールを持って攻撃しようとするんだけど、そこまでの反発力を持たずそのまま負ける」という点に関しても、検証してみたい。

データ引用元はFootball Lab, 対象期間は2020/9/30までに開催された2020シーズンのJ1全176試合である。

データの取得方法に関しては以下記事を参考
ronri-rukeichi.hatenablog.com

メイン指標の確認:追いつけないのは名古屋だけか?

20/10/2時点での、ビハインドでHTを迎えたときの平均勝ち点(Y軸)とそもそものビハインドでHTを迎える確率(X軸)をプロットしたものが下図である。

f:id:ronri_rukeichi:20201003042015p:plain
データ出典:Football Lab

ちなみに、HTをビハインドで迎える割合(X軸)の全チーム平均値=0.298 、HTでビハインドの場合の平均勝ち点(Y軸)の全チーム平均値=0.314であった。
上図から確認できるのは、以下のようなことである。

  • マッシモ名古屋は確かにビハインドでHTを迎えると勝ち点をとれない。しかし、HTビハインド時の平均勝ち点がゼロのクラブは名古屋以外にも5クラブ(セレッソ鳥栖、湘南、横浜FC、清水)ある。1/3(=6/18)のクラブはHT時のビハインドだと勝ち点を獲得できずに終わるのである。
  • 逆に川崎、柏はHT時に劣勢でも平均で1以上の勝ち点を獲得しており「先行されても逃げきられないチーム」だといえる。次いでガンバ、横浜FM、鹿島が続いている。
  • ただし名古屋はそもそもビハインドでHTを迎える確率がとても低いチーム(3/19試合=15.8%)であり、リーグ全体でみてもセレッソ鳥栖、川崎についで4位タイである。

マッシモ名古屋は先行を許すと弱いのは確かであり、川崎*4・柏といった「追いつく力」があるクラブとはその意味で差があるが、HT時にビハインドであると勝ち点を獲得できないクラブは名古屋以外にも5クラブ存在している。

少なくとも勝ち点ベースで判断したとき、マッシモ名古屋を「先行逃げ切り特化型」とするならば、ロティーセレッソも「先行逃げ切り特化型」としなければならないだろう。

関連指標の検討:リードされた後半そもそもシュートが打てているか?

ラグ氏自身も指摘している通り、そもそも名古屋(やセレッソ)はリードされて後半を迎える試合数自体が少ないので、ただ単に「運が悪い」可能性も否めない。
いわゆる「下振れ」が、ビハインドでHTを迎えた場合の「平均勝ち点ゼロ」にあらわれている、という見方だ。

シュートは打てているのに、相手GKの攻守に遭ったり、枠に嫌われたり...という可能性もある。

また、ラグ氏が指摘するところの”リードされてたらボールを持って攻撃しようとするんだけど、そこまでの反発力を持たずそのまま負ける”傾向を確かめるためには、ビハインドで後半を迎えたときのボール支配率もみなければならない。

そこで、ビハインドで後半を迎えたときの各チームのボール支配率およびシュート数を示したのが、下図である。

f:id:ronri_rukeichi:20201003170459p:plain
データ出所:Football LAB

ここから、以下のことがわかる。

  • リーグのなかでみれば、名古屋はビハインドで後半を迎えたときに極端にボールをもたされているわけでもなく*5、極端にシュートが少なくなっているわけではない
  • ただし、「ビハインドでHTを迎えている勝ち点をとれるクラブ」の代表格である川崎と柏は、確かに名古屋より少ない支配率で高いシュートを打っている
  • しかし、「ビハインドでHTを迎えている勝ち点をとれるクラブ」の第二グループであるガンバ、鹿島、横浜FMは名古屋よりも高い支配率である。よって、後半追いつくうえで、必ずしもボールを持たされること自体が致命的であるというわけではない

限定された指標ではあるが、以上からうかがえるのは、名古屋は、ビハインドで後半に突入したときの戦い(方)がそこまでリーグのなかで特異なわけではない、ということである。
しかし、追いつかれたら厳しいという現実を変えるためには、戦い方を変えねばならないかもしれない。


上図でいえば、「あまりボールを持たされ過ぎずにシュートで終わる」方向性(柏のような)に舵をきるのがよいのかな、と個人的には思った。

補足的な分析:「先行されるとズルズル型」チームの共通点

先の分析で、名古屋、セレッソ鳥栖横浜FC、湘南、清水がHTでビハインドの場合に勝ち点がゼロのクラブであると確認された。
彼らの共通点を何か見出すことができるだろうか。

ビハインド時にサポーターがフラストレーションを溜めるのは、「ボールは回る、されど進まず」というという状況である。
この点に関して、Football Labでは、試合ごとの各クラブの「敵陣ペナルティエリアへの侵入回数」のデータが揃っている。

残念ながら、時間帯別ではこの指標は得られないが、すべての試合を対象としてデータを確認してみよう。
X軸にボール支配率を、Y軸にPA侵入回数をとってプロットする。

f:id:ronri_rukeichi:20201003203303p:plain
データ出所: Football LAB

上図から、「ポゼッションの高いチームほど敵陣PAへの侵入率が高い」という傾向がみてとれる(相関係数は0.6306)。
図中の回帰直線より上にプロットされているクラブは、ボール支配率の割にPA侵入率の高いクラブであり、
回帰直線より下に位置しているクラブは、ボールを長くもっていてもあまりPA侵入できない「ボール保持の効用が低い」チームである、といえる。


ここで、HT時点でのビハインドを跳ね返せない6クラブ(名古屋、セレッソ鳥栖横浜FC、湘南、清水)の図中の位置を確認していると、全クラブが見事に回帰直線より下方に位置している*6

すなわち、HT時点でのビハインドを跳ね返せないクラブは、ボール保持を敵陣深くへの侵入につなげないという共通点があり、マッシモ名古屋もそれに該当している、ということである。

Conclusion

総括と希望:サスパサー成瀬とレシーバーシャビエル

本記事の分析では、マッシモ名古屋は「先行逃げ切り特化型」ではあるが、同様の問題を抱えているクラブは少なくない(2位のセレッソもそこに含まれる)こと、そしてその背後にはボール保持を前進に結び付ける効率の悪さが存在していること、がわかった。

他チームと比べてグランパスは極端に「ボールを持たされている」わけではないのだが、ボールを保持している単位時間あたりの前進度合いが少ない、ということである。

さて、ここまで述べてきたことだけではあまりに希望がない記事になってしまうので、最後に「何が希望となりうるか」ということについて述べて終わりたい。

先日のグラぽの記事で、フルゐ氏が「刺すパス」の重要性を指摘していた。

grapo.net


この記事では刺すパス=「相手が引いて守備ブロックを作っているときに、名古屋のディフェンダーボランチから守備ブロック中央の密集地帯に入れるグラウンダーの鋭いパス」の出し手としての成瀬の重要性が力説されている。
成瀬への期待感の本質が「刺すパス」の出し手としての資質にあるというフルゐ氏の指摘は、我々に大きな気づきを与えてくれているといってよい。

ただし「刺すパス」を試合でよく見るためには、「刺すパス」を受けるのが得意な受け手も重要となるだろう。
個人的に「刺すパス」の受け手として注目しているのは、我らがNo.10、シャビエルである。

技術の高さに目が奪われがちだが、彼の本当の強みはパス&ムーブを高速で繰り返せることにこそあると思っている。
いわゆるハーフスペースで受ける動きに関しては、他の二列目の選手にくらべて、一日の長があるのではないだろうか。

そこで、サスパサー成瀬と、レシーバーシャビエルの効用をデータからみるために、以下の二つの「攻撃の効率性」の指標をみていきたい

  • ゴール期待値(参考:Football LABによる定義
  • 枠内シュート到達率=枠内シュート/攻撃回数

ゴール期待値は、「どれだけ可能性の高いシュートを打てているか」に関しての指標であり、
枠内シュート到達率は、「どれだけ攻撃をシュートに結び付けられているか」に関しての指標である。

これを、以下の三つの場合にわけて比較する

  1. 成瀬、シャビエルとも46分以上の出場時間(5試合)
  2. 成瀬、シャビエルのどちらかが46分以上の出場時間(11試合)
  3. 成瀬、シャビエルともに45分以下の出場時間(3試合)

ケース数が少ないというデータ上の制約はある(統計的有意性云々なんて言えない))ものの、簡便な一次近似としてこの区分による比較を行う。

まず、ゴール期待値に関しては下図のとおりである。

f:id:ronri_rukeichi:20201004064628p:plain
データ出所: Football LAB

次に、枠内シュート到達率に関しては下図の通りである。

f:id:ronri_rukeichi:20201004064816p:plain
データ出所: Football LAB

これらの比較から、成瀬やシャビエルの出場が攻撃の効率性を高めていることが分かる。

成瀬/シャビエルの双方が45分以下の出場機会しか与えられない時に比べ、双方46分以上出場したときには一回の攻撃あたりの枠内シュート率は二倍以上になる

むろん攻撃だけでなく(特にマッシモにとっては)守備もあってのサッカーなのでどの試合でも成瀬やシャビエルを起用すべき!というわけではないが、「先行逃げ切り特化型」からの脱却をはかるにあたっては、興味深いデータだといえるのではないだろうか。

今後への課題

  • マッシモ名古屋の「プロビンチャ性」をデータ面から検討する
  • マッシモ名古屋がアウェーで勝てない要因の分析
  • 風間前体制との各指標/Styleの比較

一点目については、みぎ氏のブログの下記記事において問題提起がなされている。
migiright8.hatenablog.com

「そもそもマッシモは〝リスク〟の三文字が大嫌いだ。」
「彼の辞書にリスクの文字はないだろうし、あればきっと轟々と燃やし尽くすはずだ。何故危険な目にあってまでゴールを目指す?ノーリスクで勝利を目指せ。これがマッシモの哲学であり、自らリスクを冒すなど愚の誇張。リスクは冒すものではなく、狙うものだ」(みぎブログより)

また、他のサポーターの方からも以下のような指摘がみられる

要するに、骨の髄まで"プロビンチャ根性"的なものが染み込んでいるのが、マッシモさんの本質だといっていい。


この点に関し、主体的にアクションを起こすクラブを(それなりに)カモりながらも、逆にアクションを促されると右往左往してしまう今のマッシモ名古屋の特性をデータ面からとらえたい。

2点目に関しては、「たまたま」という可能性もありえる(というかそう解釈したほうが自然である)が、
イタリアのサッカー文化や思考法が染みついているマッシモが意識的/無意識的に、Home or Awayで戦い方を変えている可能性はゼロではないので、一応データをみてみたい。

風間さんとの比較に関しては、単に興味というか好奇心です。


cero / Summer Soul【OFFICIAL MUSIC VIDEO】

Enjoy!

*1:シリーズ化するかは未定

*2:ほんとうに多くはともかくとして少なくとも私は衝撃をうけた

*3:利用規約で、個人ブログにおけるデータの利用をみとめている

*4:なんとビハインドでHTを迎えても平均勝ち点1.333である。えげつない。

*5:ちなみにHT時点でビハインドでない場合とのポゼッションの変化率もリーグ平均くらいである。名古屋だけ特に「ボールを持たせる」ような対策がされているわけではない

*6:神戸も下方に位置しているが、先の図でみたように神戸はビハインド時でもシュート数を稼げている。おそらくミドルを打てる選手の存在が大きい