(可能な限り)
baseRで統計学を
楽しむ

TokyoR #105
Kien Y. Knot

誰？

• きぬいと
  • statistician-ja管理人。
  • 労働者6年目になり転職、管理職となってしまった。
  • 一緒に仕事したい人を探しています
    (ダブルミーニング)

この報告でやること

• あえてtidyverseを使わないRで統計学を楽しむ。

なぜbaseRなのか？

Rによる実践的マーケティングリサーチと分析より

base Rの能力はすべてのRユーザーにとって不可欠なので学ぶ必要があります

base Rによるコードはtidyverseによるものほど簡潔でないかもしれませんが、それは常に機能します(p9)

御託

• すべてのRユーザーにとって不可欠ならやるしかない(？)
• 多くの便利なパッケージはbaseな文法で記述
  • たとえばynakahashiさんのブログなどを見るとたのしい
• 自分で関数を作る場合はbaseな文法がオススメ
  • tidyverseで関数を作ろうとすると、
    いくつか応用的な理解が必要。
• PythonやJuliaなどにも興味があるなら割と馴染む書き方
  (個人の見解)

御託はここまで

• 「統計学」つってるけど多変量解析も検定もないです。
• とにかくこの辺書ければなんとかなるみたいなのを見ていく
• 雑な示唆出し・分析はゆるして

使用するデータ

• この本の4.1節のデータを使います
• 小売店での顧客取引を想定した模擬データです

データの読み込みread.csv

• stringAsFactorsをTRUEにするのは、本当はよくない
  • でもsummaryでは良い動きをするのであえて使う

usedata <- read.csv(
    #csvファイルならURL突っ込んでも持ってこれるよ
    "http://r-marketing.r-forge.r-project.org/data/rintro-chapter4.csv",
    stringsAsFactors = TRUE)
head(usedata, n = 5) # 先頭5行だけ見る

  cust.id      age credit.score email distance.to.store online.visits
1       1 22.89437     630.6089   yes          2.582494            20
2       2 28.04994     748.5746   yes         48.175989           121
3       3 35.87942     732.5459   yes          1.285712            39
4       4 30.52740     829.5889   yes          5.253992             1
5       5 38.73575     733.7968    no         25.044693            35
  online.trans online.spend store.trans store.spend sat.service sat.selection
1            3     58.42999           4   140.32321           3             3
2           39    756.88008           0     0.00000           3             3
3           14    250.32801           0     0.00000          NA            NA
4            0      0.00000           2    95.91194           4             2
5           11    204.69331           0     0.00000           1             1

データの確認1: データの規模

• 基本的なデータの規模(行数、列数)の確認
• 変数名もあらかじめ確認しておくと良い
  • 日本語の場合は文字化けや、文字列の"."への置換が
    起きていることもあるので、確認したほうが安全

dim(usedata) # dimentionの略

[1] 1000   12

colnames(usedata) # columns names

 [1] "cust.id"           "age"               "credit.score"     
 [4] "email"             "distance.to.store" "online.visits"    
 [7] "online.trans"      "online.spend"      "store.trans"      
[10] "store.spend"       "sat.service"       "sat.selection"    

データの確認2: 合計金額

• このデータにはonline.spendとstore.spendがある。
  • ECと店舗それぞれでの支払金額
• だいたいこれらは「合計」を見て意思決定をしたくなる
  • 平均でも良いが、平均の定義上「合計」があれば計算できる
• 後々年齢とかDMの有無とかで分解したくなってくる
  • データを絞り込まない限りは必ず合計は同じ値になる

データの確認2: 合計金額

• というわけで合計金額は控えておくと良い。

online_spend_sum <- sum(usedata$online.spend)
store_spend_sum <- sum(usedata$store.spend)
cat( # 
  paste0(
    "EC購入金額合計: ", online_spend_sum, "\n",
    "店舗購入金額合計: ", store_spend_sum 
    )
)

EC購入金額合計: 170318.241439745
店舗購入金額合計: 47582.0987232549

• ECは合計170,318円、店舗は合計47,582円。
  ECが多めという示唆も得られる。こういうことの積み重ね。

全体的な傾向を見るsummary

• summaryでデータの要約が一覧できる
  • stringsAsFactors=TRUEのご利益はここ

summary(usedata)

    cust.id            age         credit.score   email     distance.to.store 
 Min.   :   1.0   Min.   :19.34   Min.   :543.0   no :186   Min.   :  0.2136  
 1st Qu.: 250.8   1st Qu.:31.43   1st Qu.:691.7   yes:814   1st Qu.:  3.3383  
 Median : 500.5   Median :35.10   Median :725.5             Median :  7.1317  
 Mean   : 500.5   Mean   :34.92   Mean   :725.5             Mean   : 14.6553  
 3rd Qu.: 750.2   3rd Qu.:38.20   3rd Qu.:757.2             3rd Qu.: 16.6589  
 Max.   :1000.0   Max.   :51.86   Max.   :880.8             Max.   :267.0864  
                                                                              
 online.visits     online.trans      online.spend      store.trans    
 Min.   :  0.00   Min.   :  0.000   Min.   :   0.00   Min.   : 0.000  
 1st Qu.:  0.00   1st Qu.:  0.000   1st Qu.:   0.00   1st Qu.: 0.000  
 Median :  6.00   Median :  2.000   Median :  37.03   Median : 1.000  
 Mean   : 28.29   Mean   :  8.385   Mean   : 170.32   Mean   : 1.323  
 3rd Qu.: 31.00   3rd Qu.:  9.000   3rd Qu.: 177.89   3rd Qu.: 2.000  
 Max.   :606.00   Max.   :169.000   Max.   :3593.03   Max.   :12.000  
                                                                      
  store.spend      sat.service   sat.selection  
 Min.   :  0.00   Min.   :1.00   Min.   :1.000  
 1st Qu.:  0.00   1st Qu.:3.00   1st Qu.:2.000  
 Median : 30.05   Median :3.00   Median :2.000  
 Mean   : 47.58   Mean   :3.07   Mean   :2.401  
 3rd Qu.: 66.49   3rd Qu.:4.00   3rd Qu.:3.000  
 Max.   :705.66   Max.   :5.00   Max.   :5.000  
                  NA's   :341    NA's   :341    

データの分布を見るhist

• par(mfrow=c(縦,横))でグラフを並べられます

par(mfrow = c(1,2))
hist(usedata$store.spend, breaks = 30, main = "店舗支払い")
hist(usedata$online.spend, breaks = 30, main = "EC支払い")

連続値を離散値にする

• 例えば年齢とかは「○歳代」に置き換える事が多い
• ちょっとした算数を知っていると、簡単に変換できる。
• 関数を探すよりも、思いつくほうが速いこともある

# %/%は割り算の結果の整数部分を返す演算子。
usedata$ages <- paste0((usedata$age %/% 10) * 10, "代")

# 実は$だけじゃなくてこういう形でも変数を表記できる
# pandasっぽくてわかりやすいこともある
head(usedata[, c("age", "ages")])

       age ages
1 22.89437 20代
2 28.04994 20代
3 35.87942 30代
4 30.52740 30代
5 38.73575 30代
6 42.41277 40代

グループ別集計aggregate

• aggregate関数でグループ別集計が可能。
• 先に作った年代別にonline.spendの合計を出してみる

# byの中に複数の変数を入れてもOK
online_spend_by_ages <- aggregate(usedata$online.spend,
by = list(usedata$ages), sum)
print(online_spend_by_ages)

  Group.1         x
1    10代      0.00
2    20代  33380.28
3    30代 119398.47
4    40代  17539.49
5    50代      0.00

sum(online_spend_by_ages[, 2])

[1] 170318.2

• さっき見た金額合計と一致するのを確認できる。

データの関連を見るplot

• 2つの変数の関係性を見たいならplot
  • ふんわり「山なり」な関係……？

par(mfrow = c(1,2))
plot(x = usedata$age, y = usedata$store.spend, main = "年齢×店舗支払い")
plot(x = usedata$age, y = usedata$online.spend, main = "年齢×EC支払い")

離散値同士の関連を見るtable

離散値の場合tableでクロス表にしてしまうのも手。
applyをうまく使うと縦横それぞれの比率も作れる。

cross_table <- table(usedata$ages, usedata$email)
# roundで丸めができる(四捨五入ではないので注意)。
cross_table_p_col <- round(apply(cross_table, 2, function(x){x/sum(x)}), 2)
cross_table_p_row <- round(t(apply(cross_table, 1, function(x){x/sum(x)})), 2)

cbind(cross_table, cross_table_p_col, cross_table_p_row)

      no yes   no  yes   no  yes
10代   0   2 0.00 0.00 0.00 1.00
20代  27 131 0.15 0.16 0.17 0.83
30代 135 555 0.73 0.68 0.20 0.80
40代  23 126 0.12 0.15 0.15 0.85
50代   1   0 0.01 0.00 1.00 0.00

相関

• 2つの変数の関係性を見るにはcorがある
• 無相関の検定までやりたいならcor.test

print(cor(usedata$age, usedata$credit.score))

[1] 0.2545045

cor.test(usedata$age, usedata$credit.score)

    Pearson's product-moment correlation

data:  usedata$age and usedata$credit.score
t = 8.3138, df = 998, p-value = 3.008e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1955974 0.3115816
sample estimates:
      cor 
0.2545045 

相関

• 複数の変数それぞれの相関を見たい場合もcor。

usecols <- colnames(usedata[, c(2, 3, 5:10)])
cor(usedata[, usecols])

                           age credit.score distance.to.store online.visits
age                1.000000000  0.254504457        0.00198741   -0.06138107
credit.score       0.254504457  1.000000000       -0.02326418   -0.01081827
distance.to.store  0.001987410 -0.023264183        1.00000000   -0.01460036
online.visits     -0.061381070 -0.010818272       -0.01460036    1.00000000
online.trans      -0.063019935 -0.005018400       -0.01955166    0.98732805
online.spend      -0.060685729 -0.006079881       -0.02040533    0.98240684
store.trans        0.024229708  0.040424158       -0.27673229   -0.03666932
store.spend        0.003841953  0.042298123       -0.24149487   -0.05068554
                  online.trans online.spend store.trans  store.spend
age                -0.06301994 -0.060685729  0.02422971  0.003841953
credit.score       -0.00501840 -0.006079881  0.04042416  0.042298123
distance.to.store  -0.01955166 -0.020405326 -0.27673229 -0.241494870
online.visits       0.98732805  0.982406842 -0.03666932 -0.050685537
online.trans        1.00000000  0.993346657 -0.04024588 -0.052244650
online.spend        0.99334666  1.000000000 -0.04089133 -0.051690053
store.trans        -0.04024588 -0.040891332  1.00000000  0.892756851
store.spend        -0.05224465 -0.051690053  0.89275685  1.000000000

実務に連動した変換

• 店舗までの距離と、店舗での取引回数との相関は-0.2414949。
  • 店舗までの距離が遠いほど、取引回数は減るのが自然ではある
  • なんならもっと強い相関があっても良さそうだが？
• 本当にこの相関は現実を捉えられているのか？？

実務に連動した変換

• 相関だけで安易に判断しないことがいいかも
• 散布図をみると「反比例」している……？

plot(usedata$distance.to.store, usedata$store.spend)

実務に連動した変換

• 距離の逆数との相関を取ってみる

cor(x = 1 / sqrt(usedata$distance.to.store), y = usedata$store.trans)

[1] 0.5340744

plot(x=1 / sqrt(usedata$distance.to.store), y = usedata$store.spend,
xlab = "1 / √店舗までの距離", ylab = "店舗での支出額")

実務に連動した変換

• 相関はあくまで「直線的な連関」しか示唆しない
  • そうでない連関の場合相関係数は注意
• 散布図や実務・経験をもとに変換を行うことで
  示唆を得られることもある
• (経験的に)有効な変換がある
  • 売上や所得など: log(対数)変換する
  • 距離：逆数への変換をとる　など

回帰分析lm

• 「統計学」っぽいことも試して終わろう。
• DMによる施策でどれだけオンライン訪問が増えたか確かめたい
  • 施策の効果検証でよくある枠組み。
• formulaという記法y ~ xでモデル式を書く。

linear_model <- lm(
  online.visits ~ age + email + credit.score + sat.selection,
  data = usedata)

回帰分析lm

• DMが届いた人は、届かない人より平均15回程度
  オンライン訪問回数が増えそう

summary(linear_model)$coefficients

                 Estimate  Std. Error    t value    Pr(>|t|)
(Intercept)   59.15012203 26.85283195  2.2027517 0.027960321
age           -0.53306729  0.38247112 -1.3937452 0.163867965
emailyes      15.08167808  4.59455006  3.2825147 0.001083616
credit.score  -0.03003998  0.03638709 -0.8255672 0.409350715
sat.selection -1.45660311  2.00107850 -0.7279090 0.466929795

• でも、結果を安易に信じてはいけない
  • 統計的な因果関係の検証は実際難しい
  • オンラインを使う人にDM送ってない？(選択バイアス)

まとめ1

• baseな書き方で基本的な集計を楽しもう
  • 簡単な確認くらいなら慣れたらこっちのほうが速い
• ただデータを眺めるだけでも示唆が得られるが、
  現実に合わせて数値の変換を考えられるとより洞察が深まることもある
  • 割とこのあたり簡単な算数で表現できがち

まとめ2

• 手法は高度になるほど考えることが増える
  • 安易に使うと意思決定を誤る
• とはいえ使わないことには何も覚えない
  • いっぱい使っていっぱい失敗しよう。
  • わかる人は失敗を許し、守ってあげよう(自戒)

Enjoy!

• 発表では言わないけど、次スライド以降で
  仕事で使えるパッケージを列挙してます。

お仕事でお世話になるシリーズ①

• tidyverse
  • 操作の直感性が高い。ググれば使い方がいっぱい
• data.table
  • クセがあるけど大規模なデータの処理が得意
  • Uryuさんのスライドを見るといい
• summarytools
  • summaryをよりきれいにわかりやすく出してくれる
  • pythonでいうpandas_profilingに近い

お仕事でお世話になるシリーズ②

• polars
  • 使ってみるとtidyverseの100倍速いときがある
• tidymodels
  • 今日も発表があると思う。
  • 特に機械学習モデル構築を便利にする
• psych
  • いぶし銀の因子分析パッケージ
  • その他データ要約も得意

このスライドは

• Quartoで記述しています。
  • なんかいいドキュメンテーションライブラリ
  • Pythonとかも書いてやっていけるので楽しい。