広告/統計/アニメ/映画 等に関するブログ

広告/統計/アニメ/映画 等に関するブログ

TVCMが昨対比を何%押し上げたか?をGoogle Spread Sheet やEXCELで(敢えて乱暴に)分析する

「テレビCMやったけどさ、どんだけ効果があったのよ?」と質問されたときにどう回答するのが適切だろうか。

データ分析ができる人にとってはデータさえ揃えば難しい話ではない。投下したTVCMのGRPを説明変数として売上なり申込件数を目的変数とした予測モデルを作れば良い。しかし、データ分析の素養がない人に予測モデルの式本体を見せたところでチンプンカンプンだろうし、報告する人が普通の人だった場合、RやPythonモデリングしろ、というわけにもいかない。

例えば、緊急事態宣言下の4月にオンライン教材を売っていた企業が、TVCMを放映してみたエリアと放映してみなかったエリアでA/Bテストしていたとしよう。1 その場合、どのようなデータ分析をすればよいか考えてみたい

モデルの検討

報告すべき情報、わかりたい情報を決める

売上を予測するモデルを作ることも、TVCM効果の有無を検定する方法もあるだろう。しかし管理職以上は統計学もデータ分析の素養はないし、効果は”量的に”把握をしたいだろう。そして今回のケースは、「それって緊急事態宣言の影響もあるんじゃない?」とか「そもそも今年がオンライン学習元年でよそも伸びてない?」と指摘されるだろうから、それとも区別したい。

そこで、 昨対比で申込件数が何%あがったか を求めることにした。

取得できるデータを確認する

今回入手できる情報は、"各県で緊急事態宣言があったかどうか"と"TVCMを放映したかどうか" といういずれもゼロと1の数値。そして、本年度の申込件数と昨年度の前年同期の申込件数だ。2

f:id:yyhhyy:20200815212637j:plain

手法を決める

今回はEXCELGoogle Spread Sheetで可能な作業として、「線形の重回帰分析」とした。

  • 目的変数:申込件数の昨対比(本年度の申込件数/昨年度の申込件数*100:→%にした)
  • 説明変数:「TVCM放映の有無」と「緊急事態宣言の有無」
  • 切片:あり

これはハッキリ言って嘘である。

「y = a1x1 + a2x2 +b」といったお互いに独立の変数が線形で関係している、という式を検討しているわけだが、昨対比などというものは、

  • TVCMも緊急事態宣言もあった場合:a1 * a2倍(%)
  • TVCMを放映し緊急事態宣言がなかった場合:a1倍(%)
  • TVCM放映がなく緊急事態宣言はあった場合:a2倍(%)
  • TVCM放映がなく緊急事態宣言もなかった場合:1倍(100%)

というわけだから、2つの変数が「足し算」になっているはずがない。

とはいえ、掛け算を足し算に戻すには、対数を取らなければならないし、対数変換についてこれる文系役員などそうそう居ない。3 かといって再び「昨対比」を避けてしまうと、

「本年度の申込件数/一定人口辺り=昨年度の申込件数/一定人口辺り + TVCMが加えてくれる件数 + 緊急事態宣言が加えてくれる件数 + その他の要因」

という効果の量的なイメージが判断し辛いものが出来上がる4

だが、この手法であれば、TVCMと緊急事態宣言の効果は一応別々で出てくるし、EXCELGoogle Spread Sheetで結果が出てくるのだから、よくわからない人でも見様見真似で算出することができる。

分析の実際

Google Spread Sheetの作業

  • LINEST(目的変数の列,説明変数の列群,定数項の有無,出力結果の詳細の有無)

をシートに入力する

今回の場合は、

  • 目的変数:昨対比
  • 説明変数の列群:TVCM放映の有無 と 緊急事態宣言の有無

f:id:yyhhyy:20200815222551j:plain

となる。

結果の見方

このような出力が出てくる。重要なところに色付けをした。

f:id:yyhhyy:20200815222648j:plain

99.3 は切片に当たる。昨対比を%表記にしたので切片は本来100であるべきだが少しずれている。

31.3は一つめの変数「TVCM放映の有無」に掛かる係数だ。「TVCM放映は、約31%押し上げる効果があった」と読みとける

11.7は二つ目の変数「緊急事態宣言の有無」に掛かる係数で、「緊急事態宣言があったエリアでは、約12%押し上げることになった」と読める。

もう1箇所網掛けした0.994は、これらの変数2つによって目的変数(の分散)を約99%は説明できている、と解釈する。

なお、Google Spread Sheet のLINEST関数の出力はかなり直感に反している。1つ目の変数の係数→2つ目の変数の係数と左から並んで欲しいものだが、なぜか右から、切片→1つ目の変数→2つ目の変数の係数と並ぶ。わりと面倒だ。

その他の出力の値に興味がある人は公式ドキュメントを参照してほしい

support.google.com

EXCELの場合

EXCELにはもう少しグラフィカルインターフェイスな機能がある。EXCELの場合は丁寧にまとめている人が多いので割愛する。

sigma-eye.com

今回のデータを作った方法

EXCELで乱数を発生させた。

  • TVCM放映の有無と緊急事態宣言の有無:RANDBETWEEN(0,1)
    • ゼロと1の乱数を発生させてくれる。
  • 2019年度の申込件数:RANDBETWEEN(100,150)
    • 使用する関数は同じ。数値の上下を変えた

効果の数値設定と目的変数の作成

TVCMの効果を1.3倍、緊急事態宣言の効果を1.1倍として、

  • TVCM効果:IF(TVCM放映の有無=1,1.3,1)
    • 放映があった場合は1.3、なければ1を掛け算したいため
  • 緊急事態宣言効果:IF(緊急事態宣言の有無=1,1.1,1)
    • 放映があった場合は1.1、なければ1を掛け算したいため

という新たな変数を作り

  • 2020年度の申込件数:ROUND(2019年度申込件数TVCM放映効果緊急事態宣言効果+RANDBETWEEN(-2,2),0)
    • 昨年度の申込件数にTVCM放映効果と緊急事態宣言効果を乗算、更に2件程度は上下にブレるように乱数を加えた上で、最後に小数点以下四捨五入している

を求めている。こんな感じだ。因みにリロードしてしまったので上で貼り付けた数値と乱数がまた更新されている。

f:id:yyhhyy:20200815225134j:plain

先程はかなり乱暴に線形の重回帰分析モデルにしてしまったのだが、31%や12%は上振れしているものの、そこまで真の値からは離れていない。統計学者やデータサイエンティストからすれば余りにも気持ち悪い話かもしれないが、"偉い人に普通の人が効果を説明する" 行為を日常的に行うことがデータドリブン経営なのだということとなると、この程度のいい加減さを我慢しないといつまでも組織として根付くことはないのではないだろうか。この先、データ分析の教養がある学生が社会人になってくれるようになれば話はべつかもしれないが、その間、今の20代~60代の文系社会人が何もしなくて良いということにもならないだろう。5


  1. 極めて当たり前だが、放映しなかったエリアが存在しなければ効果検証できない。そして多くの日本企業のマーケティング・プロモーションは、事前の効果検証もせずに一発で一か八かの勝負をするので、"TVCMを放映しないエリア"なんてものが存在しないことの方が寧ろ多い。この問題の根本原因はデータ分析リテラシーの欠如ではなく、「スケジュールに余裕がなくなるレベルの組織の意思決定の遅さ」「A/Bテストによって発生するエリア別の不公平を許容できないビジネスモデルや取引先との力関係」などの方だ。

  2. 簡単のため各エリアの人口は同じだとしている。実際には人口○○人辺りや、ターゲット○○人辺りなどにすればよい。

  3. 現実には多項ロジスティック回帰分析をしたいことも多いし、対数変換を避けてでも嘘をつくべきかは非常に悩ましい。

  4. 本当を言うとStanを使ってシミュレーションしてしまえば、効果量を設定して推定値を出せば良い。データ分析がわかる自分が毎回計算してあげれば良い、ということであればこれで構わないし、自分が効果の程を把握したいだけならそれで充分だ。だが組織としてそれで良いのか?というとそうではないだろう。

  5. “文系に統計学の教養がない"という前提も実は嘘だ。心理学・社会学では日常的に使うし、最近は経済学部でも因果推論という名称で社会学的アプローチを参照し始めている。更には文学部でも文献をデータ分析して文書の著者を推定するなど統計学をフル活用する分野もある。とはいえ、理系でデータ分析を使わないのは建築学部のデザイン系だとかかなり限られていると思うので、総じて文系がビジネスパーソンとしての教養不足で苦しむことになると思う。

GPUパソコンにtensorflow-gpuを入れたはなし

最近GeforceGPUが入ったPC(もはや名前はワークステーションでした)を買って、暫くずっと`tensorflow-gpuが入らない、と苦労していたのですが、

偶然検索したこちらの記事の通り、Anaconda Prompt に任せたら一瞬で適切なものをインストールしてくれました。

qiita.com

検索で最初にヒットするサイトの多くが自力でNVIDAドライバを入れさせる手法の紹介なのですが、今では上記の方法で充分です。

時代の変化は早いですね!

Pythonでアンケート調査のクラスター分析と決定木分析を行う

アンケート調査の分析をするのはマーケティング担当者で、恐らく大学時代は社会学や心理学といった文系出身だと思います。昔ならSPSS、最近ならRだと思います。

一方で、Pythonはどちらかというと情報学系の人やシステムエンジニアが使うツール(言語)でPythonでアンケート分析を真っ向からしている書籍は存外少ないものです。最近私はRからPythonへの全面的な移行を考えているのですが、備忘録も兼ねて、Pythonでアンケート調査を行ってみました。

事前準備・前処理

先ずは予め読み込んでおいた方が良いLibrary類をインポートしておきます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

seabornのテーマをデフォルトで選ぶようにしておきます。

#sns.set()
sns.set(font="Noto Sans CJK JP")

データの準備

アンケート調査のクラスター分析のためには集計されたクロス表ではなく、その前のローデータ(0,1の回答レコードの行列)が必要なのですが、なかなか個票データを公開しているところがなく、UCIのサイトから拝借することにしました。1

http://archive.ics.uci.edu/ml/datasets.php

ある都市に対する満足度の評価です

http://archive.ics.uci.edu/ml/datasets/Somerville+Happiness+Survey

提供されているCSVファイルがWindows文字コードだったようですのでエンコードの指定をして読み込みます2

df = pd.read_csv("data/SomervilleHappinessSurvey2015_2.csv",
                 encoding="cp932",
                 header=0)
df.head(3)
D X1 X2 X3 X4 X5 X6
0 0 3 3 3 4 2 4
1 0 3 2 3 5 4 3
2 1 5 3 3 3 3 5

各設問についてはMCUのサイトに以下の補足があります。 今のままではわかりにくいので簡単な列名に変えたいと思います3

D = decision attribute (D) with values 0 (unhappy) and 1 (happy) 
X1 = the availability of information about the city services 
X2 = the cost of housing 
X3 = the overall quality of public schools 
X4 = your trust in the local police 
X5 = the maintenance of streets and sidewalks 
X6 = the availability of social community events 

Attributes X1 to X6 have values 1 to 5.
df = df.rename(columns={
    "D":"幸福かどうか",
    "X1":"行政サービス情報へのアクセスしやすさ",
    "X2":"住宅供給の高さ",
    "X3":"公立の学校の全般的な質の良さ",
    "X4":"地域警察への信頼の高さ",
    "X5":"道路・歩道のメンテナンス状況",
    "X6":"地域社会行事の利用のしやすさ"})
df.head(3)
幸福かどうか 行政サービス情報へのアクセスしやすさ 住宅供給の高さ 公立の学校の全般的な質の良さ 地域警察への信頼の高さ 道路・歩道のメンテナンス状況 地域社会行事の利用のしやすさ
0 0 3 3 3 4 2 4
1 0 3 2 3 5 4 3
2 1 5 3 3 3 3 5

何を分析するか?」をまだ決めていませんでしたが、例えば、「幸福度の高い人は市のどこを評価しているのか?」を探ることにしてみましょう。

その場合、「幸福だ」と答えているデータだけに絞る必要があります。

df_happy = df[df["幸福かどうか"]==1]
df_happy.head(3)
幸福かどうか 行政サービス情報へのアクセスしやすさ 住宅供給の高さ 公立の学校の全般的な質の良さ 地域警察への信頼の高さ 道路・歩道のメンテナンス状況 地域社会行事の利用のしやすさ
2 1 5 3 3 3 3 5
5 1 5 5 3 5 5 5
7 1 5 4 4 4 4 5

念のためにこれで何件あるか確認しましょう。あまりデータが少なすぎると分析の信憑性が下がります

len(df_happy)
77

また、幸福ではない人のデータは今回の分析では使わないので、「幸福かどうか」の列を削除します。

更に、インデックスの番号を新しいデータに対して振り直します。

df_happy = df_happy.drop("幸福かどうか", axis=1)
df_happy = df_happy.reset_index(drop=True)
df_happy.head(3)
行政サービス情報へのアクセスしやすさ 住宅供給の高さ 公立の学校の全般的な質の良さ 地域警察への信頼の高さ 道路・歩道のメンテナンス状況 地域社会行事の利用のしやすさ
0 5 3 3 3 3 5
1 5 5 3 5 5 5
2 5 4 4 4 4 5

データの標準化

5段階スケールできいている各設問について、「1」に集中しているものや、「5」に集中しているものがある筈です。4

df_happy.describe()
行政サービス情報へのアクセスしやすさ 住宅供給の高さ 公立の学校の全般的な質の良さ 地域警察への信頼の高さ 道路・歩道のメンテナンス状況 地域社会行事の利用のしやすさ
count 77.000000 77.000000 77.000000 77.000000 77.000000 77.000000
mean 4.545455 2.558442 3.415584 3.792208 3.831169 4.389610
std 0.679502 1.117958 1.004603 0.878660 1.056342 0.763576
min 3.000000 1.000000 1.000000 1.000000 1.000000 1.000000
25% 4.000000 2.000000 3.000000 3.000000 3.000000 4.000000
50% 5.000000 2.000000 3.000000 4.000000 4.000000 5.000000
75% 5.000000 3.000000 4.000000 4.000000 5.000000 5.000000
max 5.000000 5.000000 5.000000 5.000000 5.000000 5.000000

どの列も分散を等しくすることで、情報量を同じにできます。5

from sklearn import preprocessing
ss = preprocessing.StandardScaler()
df_happy_s = pd.DataFrame(ss.fit_transform(df_happy))
df_happy_s.head(3)
0 1 2 3 4 5
0 0.673326 0.397559 -0.416393 -0.907521 -0.791997 0.804625
1 0.673326 2.198267 -0.416393 1.383598 1.113745 0.804625
2 0.673326 1.297913 0.585552 0.238038 0.160874 0.804625

また列名が消えてしまいましたので振り直します。

幸い、先程の「df_happy」の列名をそのままコピペすればいいだけです

print(df_happy.columns)
Index(['行政サービス情報へのアクセスしやすさ', '住宅供給の高さ', '公立の学校の全般的な質の良さ', '地域警察への信頼の高さ',
       '道路・歩道のメンテナンス状況', '地域社会行事の利用のしやすさ'],
      dtype='object')
df_happy_s.columns = df_happy.columns
df_happy_s.head(3)
行政サービス情報へのアクセスしやすさ 住宅供給の高さ 公立の学校の全般的な質の良さ 地域警察への信頼の高さ 道路・歩道のメンテナンス状況 地域社会行事の利用のしやすさ
0 0.673326 0.397559 -0.416393 -0.907521 -0.791997 0.804625
1 0.673326 2.198267 -0.416393 1.383598 1.113745 0.804625
2 0.673326 1.297913 0.585552 0.238038 0.160874 0.804625

分析

因子分析

設問が多岐に渡るときは、設問を縮約すべきで、そういうときは因子分析を行いますが、今回は既に6問と少ないので行いません。

↓自分がやるまでもなく丁寧にまとめられた記事がありました。因子得点は fit_transformで出てきます。また、事前に自分で標準化していないと変な値になるので因子分析だけをするときも注意して下さい。

hk29.hatenablog.jp

scikit-learn.org

クラスター分析

「幸福度の高い人は市のどこを評価しているのか?」を探ると言っても、人によって重視ポイントが異なります。それを幾つかにタイプ分けをして把握するためにクラスター分析を行います

適切なクラスター数を探す

適切なクラスター数を決める手法に完全な決まりはありませんが、見た目で判断するために、階層的クラスタ分析を一度行うという手法がママあります。

from scipy.cluster.hierarchy import linkage, dendrogram

統計に細かい人にとってはどの手法を選ぶか?は重要だと思いますが、比較的一般的なユークリッド距離、ウォード法で階層的クラスター分析を行います

df_happy_s_hclust = linkage(df_happy_s,metric="euclidean",method="ward")
plt.figure(figsize=(12,8))
dendrogram(df_happy_s_hclust)
plt.savefig('figure_1.png')
plt.show()

f:id:yyhhyy:20190707224159p:plain

どの辺りで線をひくか?はかなり恣意的ですが、この場合は4グループぐらいが適切でしょうか?

f:id:yyhhyy:20190707224239p:plain

k平均法によるクラスター分析

実務においてはクラスターに分かれればそれでいいというわけではなく、解釈の容易性が求められます。それぞれのクラスターにはどういう違いがあるのか?そういったことを知るには、kmenasクラスター分析の方が便利です

from sklearn.cluster import KMeans

クラスター数を4つと決めたので引数に入れて関数を作成

km = KMeans(n_clusters=4,random_state=42)

skit-learnは直接pandasデータフレームを読み込まないのでnumpyの行列に変換する必要がある

df_happy_s_ar = df_happy_s.values
display(df_happy_s_ar)
array([[ 0.67332598,  0.39755886, -0.41639284, -0.90752108, -0.79199672,
         0.80462467],
       [ 0.67332598,  2.19826665, -0.41639284,  1.38359771,  1.11374539,
         0.80462467],
       [ 0.67332598,  1.29791276,  0.58555244,  0.23803832,  0.16087433,
         0.80462467],
       [ 0.67332598, -0.50279503,  0.58555244,  1.38359771,  1.11374539,
         0.80462467],
       [-2.28930833, -0.50279503,  0.58555244, -0.90752108,  0.16087433,
        -0.51359022],
       [ 0.67332598, -1.40314893,  0.58555244, -0.90752108,  0.16087433,
         0.80462467],
       [-0.80799118, -0.50279503, -0.41639284, -0.90752108,  0.16087433,
        -0.51359022],
       [-0.80799118, -0.50279503, -0.41639284, -0.90752108,  0.16087433,
        -0.51359022],
       [ 0.67332598, -1.40314893, -1.41833812,  1.38359771, -1.74486778,
        -0.51359022],
       [-0.80799118,  0.39755886, -0.41639284, -0.90752108, -0.79199672,
        -0.51359022],
       [-2.28930833,  0.39755886, -0.41639284,  1.38359771,  1.11374539,
         0.80462467],
       [-2.28930833,  0.39755886, -2.4202834 , -0.90752108, -0.79199672,
        -0.51359022],
       [-2.28930833,  0.39755886, -2.4202834 , -0.90752108, -0.79199672,
        -0.51359022],
       [ 0.67332598,  0.39755886, -0.41639284, -0.90752108,  1.11374539,
        -1.83180511],
       [-2.28930833, -0.50279503,  0.58555244,  0.23803832,  0.16087433,
         0.80462467],
       [-2.28930833, -0.50279503,  0.58555244,  0.23803832,  0.16087433,
         0.80462467],
       [-0.80799118, -1.40314893, -0.41639284, -3.19863987, -2.69773883,
        -0.51359022],
       [ 0.67332598,  0.39755886,  0.58555244, -0.90752108,  0.16087433,
         0.80462467],
       [ 0.67332598,  0.39755886,  0.58555244, -0.90752108,  0.16087433,
         0.80462467],
       [ 0.67332598, -0.50279503, -0.41639284, -0.90752108, -1.74486778,
         0.80462467],
       [-0.80799118,  1.29791276, -0.41639284,  0.23803832, -1.74486778,
        -0.51359022],
       [-0.80799118, -0.50279503,  0.58555244, -0.90752108, -1.74486778,
        -0.51359022],
       [-2.28930833, -1.40314893, -1.41833812,  0.23803832, -0.79199672,
         0.80462467],
       [ 0.67332598,  0.39755886,  0.58555244,  0.23803832,  0.16087433,
         0.80462467],
       [ 0.67332598,  0.39755886, -0.41639284,  0.23803832,  0.16087433,
         0.80462467],
       [ 0.67332598, -0.50279503,  1.58749772,  1.38359771,  1.11374539,
        -1.83180511],
       [ 0.67332598, -1.40314893, -0.41639284, -0.90752108,  0.16087433,
        -0.51359022],
       [ 0.67332598, -1.40314893, -0.41639284, -0.90752108,  0.16087433,
        -0.51359022],
       [ 0.67332598, -1.40314893, -0.41639284, -0.90752108,  0.16087433,
        -0.51359022],
       [ 0.67332598, -0.50279503,  0.58555244, -0.90752108,  0.16087433,
         0.80462467],
       [ 0.67332598, -0.50279503,  0.58555244, -0.90752108,  0.16087433,
         0.80462467],
       [-0.80799118,  0.39755886, -1.41833812,  0.23803832, -0.79199672,
        -0.51359022],
       [-0.80799118,  0.39755886, -1.41833812,  0.23803832, -0.79199672,
        -0.51359022],
       [ 0.67332598,  0.39755886,  1.58749772,  1.38359771,  0.16087433,
         0.80462467],
       [ 0.67332598, -0.50279503,  0.58555244, -2.05308048, -1.74486778,
        -0.51359022],
       [-0.80799118,  1.29791276, -0.41639284, -0.90752108, -1.74486778,
         0.80462467],
       [ 0.67332598, -0.50279503,  0.58555244,  0.23803832,  1.11374539,
         0.80462467],
       [ 0.67332598, -0.50279503,  0.58555244,  0.23803832,  1.11374539,
         0.80462467],
       [-0.80799118, -1.40314893, -0.41639284,  0.23803832,  0.16087433,
        -0.51359022],
       [-0.80799118, -1.40314893, -0.41639284,  0.23803832,  0.16087433,
        -0.51359022],
       [ 0.67332598, -1.40314893,  1.58749772,  1.38359771,  1.11374539,
         0.80462467],
       [ 0.67332598,  1.29791276,  1.58749772,  1.38359771,  1.11374539,
         0.80462467],
       [-0.80799118, -0.50279503, -1.41833812,  0.23803832,  0.16087433,
         0.80462467],
       [-0.80799118,  0.39755886, -0.41639284,  0.23803832, -0.79199672,
        -0.51359022],
       [-0.80799118,  0.39755886, -0.41639284,  0.23803832, -1.74486778,
        -0.51359022],
       [ 0.67332598, -0.50279503,  1.58749772,  1.38359771,  1.11374539,
         0.80462467],
       [ 0.67332598, -0.50279503, -0.41639284,  1.38359771,  1.11374539,
         0.80462467],
       [ 0.67332598,  2.19826665,  1.58749772,  1.38359771,  1.11374539,
         0.80462467],
       [-0.80799118,  0.39755886, -1.41833812,  0.23803832,  0.16087433,
        -0.51359022],
       [-0.80799118, -0.50279503,  0.58555244,  0.23803832,  0.16087433,
        -0.51359022],
       [ 0.67332598,  0.39755886, -1.41833812,  0.23803832,  0.16087433,
        -0.51359022],
       [ 0.67332598, -0.50279503, -0.41639284,  0.23803832,  0.16087433,
         0.80462467],
       [ 0.67332598,  0.39755886, -0.41639284, -0.90752108,  1.11374539,
         0.80462467],
       [ 0.67332598, -0.50279503,  0.58555244,  1.38359771,  0.16087433,
         0.80462467],
       [ 0.67332598, -1.40314893, -0.41639284,  0.23803832,  1.11374539,
         0.80462467],
       [ 0.67332598,  1.29791276,  1.58749772,  1.38359771,  1.11374539,
        -0.51359022],
       [ 0.67332598,  2.19826665, -0.41639284,  0.23803832,  0.16087433,
         0.80462467],
       [ 0.67332598,  1.29791276,  0.58555244, -0.90752108, -0.79199672,
        -0.51359022],
       [ 0.67332598,  1.29791276,  0.58555244,  0.23803832,  0.16087433,
         0.80462467],
       [ 0.67332598,  2.19826665,  1.58749772,  1.38359771,  1.11374539,
         0.80462467],
       [ 0.67332598,  1.29791276,  1.58749772,  0.23803832,  1.11374539,
        -0.51359022],
       [-0.80799118, -0.50279503, -0.41639284,  0.23803832, -0.79199672,
        -1.83180511],
       [ 0.67332598, -0.50279503, -1.41833812,  0.23803832,  1.11374539,
         0.80462467],
       [ 0.67332598,  0.39755886, -1.41833812,  0.23803832,  0.16087433,
         0.80462467],
       [ 0.67332598,  0.39755886,  0.58555244, -0.90752108,  0.16087433,
        -1.83180511],
       [-0.80799118,  0.39755886,  0.58555244,  0.23803832, -0.79199672,
        -0.51359022],
       [ 0.67332598, -1.40314893,  0.58555244, -0.90752108,  1.11374539,
         0.80462467],
       [ 0.67332598, -1.40314893,  1.58749772, -0.90752108,  1.11374539,
         0.80462467],
       [ 0.67332598,  0.39755886,  0.58555244,  0.23803832,  0.16087433,
        -0.51359022],
       [ 0.67332598, -0.50279503,  0.58555244,  0.23803832, -1.74486778,
        -1.83180511],
       [-2.28930833,  1.29791276,  0.58555244,  1.38359771, -2.69773883,
        -1.83180511],
       [ 0.67332598, -1.40314893,  1.58749772,  1.38359771,  1.11374539,
         0.80462467],
       [-0.80799118,  0.39755886, -0.41639284,  0.23803832,  0.16087433,
        -0.51359022],
       [ 0.67332598,  2.19826665, -2.4202834 , -3.19863987,  1.11374539,
        -4.4682349 ],
       [ 0.67332598, -0.50279503, -0.41639284,  0.23803832,  0.16087433,
        -1.83180511],
       [ 0.67332598, -0.50279503, -0.41639284,  0.23803832, -1.74486778,
         0.80462467],
       [ 0.67332598,  0.39755886, -0.41639284,  0.23803832,  0.16087433,
         0.80462467]])

kmeansを適用した結果のグルーピングの配列が出力として渡される

df_happy_s_ar_pred = km.fit_predict(df_happy_s_ar)
display(df_happy_s_ar_pred)
array([3, 3, 3, 0, 1, 0, 2, 2, 1, 1, 1, 1, 1, 2, 1, 1, 2, 3, 3, 2, 1, 2,
       1, 3, 3, 0, 2, 2, 2, 0, 0, 1, 1, 3, 2, 1, 0, 0, 2, 2, 0, 3, 1, 1,
       1, 0, 0, 3, 1, 1, 2, 0, 3, 0, 0, 3, 3, 3, 3, 3, 3, 2, 0, 3, 2, 1,
       0, 0, 3, 2, 1, 0, 1, 2, 2, 1, 3])

割り振られた結果を元の(標準化する前の)データにクラスターIDとして一列追加する

df_happy_clust = df_happy[:]
df_happy_clust["cluster_ID"] = df_happy_s_ar_pred
df_happy_clust.head(3)
行政サービス情報へのアクセスしやすさ 住宅供給の高さ 公立の学校の全般的な質の良さ 地域警察への信頼の高さ 道路・歩道のメンテナンス状況 地域社会行事の利用のしやすさ cluster_ID
0 5 3 3 3 3 5 3
1 5 5 3 5 5 5 3
2 5 4 4 4 4 5 3

因みにこのままではクラスターIDも数値として扱われてしまいます

print(df_happy_clust.dtypes)
行政サービス情報へのアクセスしやすさ    int64
住宅供給の高さ               int64
公立の学校の全般的な質の良さ        int64
地域警察への信頼の高さ           int64
道路・歩道のメンテナンス状況        int64
地域社会行事の利用のしやすさ        int64
cluster_ID            int32
dtype: object

カテゴリカル変数に変えておきましょう

df_happy_clust["cluster_ID"] = df_happy_clust["cluster_ID"].astype("category")
print(df_happy_clust.dtypes)
行政サービス情報へのアクセスしやすさ       int64
住宅供給の高さ                  int64
公立の学校の全般的な質の良さ           int64
地域警察への信頼の高さ              int64
道路・歩道のメンテナンス状況           int64
地域社会行事の利用のしやすさ           int64
cluster_ID            category
dtype: object

クラスタが何人くらいになったのかを把握しておきます

print(df_happy_clust["cluster_ID"].value_counts())
1    22
3    20
2    18
0    17
Name: cluster_ID, dtype: int64

まぁまぁまどのクラスターも同じ人数ですね

クラスターの分析

クラスターでグルーピングし、回答の片寄りがどこに出たのか?を確認します。

今回は、間隔尺度の設問でしたので単純に平均するのは乱暴ですので、踏まえて更にデータを再構築する必要があります。

df_happy_clust = df_happy_clust[:].astype("category")
print(df_happy_clust.dtypes)
行政サービス情報へのアクセスしやすさ    category
住宅供給の高さ               category
公立の学校の全般的な質の良さ        category
地域警察への信頼の高さ           category
道路・歩道のメンテナンス状況        category
地域社会行事の利用のしやすさ        category
cluster_ID            category
dtype: object

ダミー変数化したい列を指定するために列名を取得してリスト化、更にクラスタIDを除きます

dummy_list = list(df_happy_clust.columns)[0:-1]
print(dummy_list)
['行政サービス情報へのアクセスしやすさ', '住宅供給の高さ', '公立の学校の全般的な質の良さ', '地域警察への信頼の高さ', '道路・歩道のメンテナンス状況', '地域社会行事の利用のしやすさ']

ダミー変数化したい列名を指定して全ての設問をダミー変数化します6

df_happy_clust_dmy = pd.get_dummies(df_happy_clust,columns=dummy_list)
df_happy_clust_dmy.head(3)
cluster_ID 行政サービス情報へのアクセスしやすさ_3 行政サービス情報へのアクセスしやすさ_4 行政サービス情報へのアクセスしやすさ_5 住宅供給の高さ_1 住宅供給の高さ_2 住宅供給の高さ_3 住宅供給の高さ_4 住宅供給の高さ_5 公立の学校の全般的な質の良さ_1 ... 地域警察への信頼の高さ_5 道路・歩道のメンテナンス状況_1 道路・歩道のメンテナンス状況_2 道路・歩道のメンテナンス状況_3 道路・歩道のメンテナンス状況_4 道路・歩道のメンテナンス状況_5 地域社会行事の利用のしやすさ_1 地域社会行事の利用のしやすさ_3 地域社会行事の利用のしやすさ_4 地域社会行事の利用のしやすさ_5
0 3 0 0 1 0 0 1 0 0 0 ... 0 0 0 1 0 0 0 0 0 1
1 3 0 0 1 0 0 0 0 1 0 ... 1 0 0 0 0 1 0 0 0 1
2 3 0 0 1 0 0 0 1 0 0 ... 0 0 0 0 1 0 0 0 0 1

3 rows × 28 columns

クラスタIDでグループ化し数値を集約します

df_happy_clust_dmy_gp = df_happy_clust_dmy.groupby("cluster_ID")

グループ別に各設問の回答者数の合計を出します

df_happy_clust_dmy_gp_g = df_happy_clust_dmy_gp.sum().T
display(df_happy_clust_dmy_gp_g)
cluster_ID 0 1 2 3
行政サービス情報へのアクセスしやすさ_3 0 8 0 0
行政サービス情報へのアクセスしやすさ_4 0 12 7 0
行政サービス情報へのアクセスしやすさ_5 17 2 11 20
住宅供給の高さ_1 6 2 6 0
住宅供給の高さ_2 11 6 8 0
住宅供給の高さ_3 0 11 3 10
住宅供給の高さ_4 0 3 0 6
住宅供給の高さ_5 0 0 1 4
公立の学校の全般的な質の良さ_1 0 2 1 0
公立の学校の全般的な質の良さ_2 1 6 1 1
公立の学校の全般的な質の良さ_3 3 8 12 6
公立の学校の全般的な質の良さ_4 8 6 4 7
公立の学校の全般的な質の良さ_5 5 0 0 6
地域警察への信頼の高さ_1 0 0 2 0
地域警察への信頼の高さ_2 0 0 1 0
地域警察への信頼の高さ_3 5 5 9 5
地域警察への信頼の高さ_4 5 14 6 9
地域警察への信頼の高さ_5 7 3 0 6
道路・歩道のメンテナンス状況_1 0 1 1 0
道路・歩道のメンテナンス状況_2 0 5 4 0
道路・歩道のメンテナンス状況_3 0 8 1 2
道路・歩道のメンテナンス状況_4 5 7 10 11
道路・歩道のメンテナンス状況_5 12 1 2 7
地域社会行事の利用のしやすさ_1 0 0 1 0
地域社会行事の利用のしやすさ_3 1 1 5 0
地域社会行事の利用のしやすさ_4 0 14 11 4
地域社会行事の利用のしやすさ_5 16 7 1 16

最近読んだ本で知ったのですが、Jupyter notebookであれば、HTML上で各セルに棒グラフを入れたり、数値によって色をつけたりしてくれるので、是非活用したほうが良いです。

pandas.pydata.org

df_happy_clust_dmy_gp_g.style.bar(color="#4285F4")
cluster_ID 0 1 2 3
行政サービス情報へのアクセスしやすさ_3 0 8 0 0
行政サービス情報へのアクセスしやすさ_4 0 12 7 0
行政サービス情報へのアクセスしやすさ_5 17 2 11 20
住宅供給の高さ_1 6 2 6 0
住宅供給の高さ_2 11 6 8 0
住宅供給の高さ_3 0 11 3 10
住宅供給の高さ_4 0 3 0 6
住宅供給の高さ_5 0 0 1 4
公立の学校の全般的な質の良さ_1 0 2 1 0
公立の学校の全般的な質の良さ_2 1 6 1 1
公立の学校の全般的な質の良さ_3 3 8 12 6
公立の学校の全般的な質の良さ_4 8 6 4 7
公立の学校の全般的な質の良さ_5 5 0 0 6
地域警察への信頼の高さ_1 0 0 2 0
地域警察への信頼の高さ_2 0 0 1 0
地域警察への信頼の高さ_3 5 5 9 5
地域警察への信頼の高さ_4 5 14 6 9
地域警察への信頼の高さ_5 7 3 0 6
道路・歩道のメンテナンス状況_1 0 1 1 0
道路・歩道のメンテナンス状況_2 0 5 4 0
道路・歩道のメンテナンス状況_3 0 8 1 2
道路・歩道のメンテナンス状況_4 5 7 10 11
道路・歩道のメンテナンス状況_5 12 1 2 7
地域社会行事の利用のしやすさ_1 0 0 1 0
地域社会行事の利用のしやすさ_3 1 1 5 0
地域社会行事の利用のしやすさ_4 0 14 11 4
地域社会行事の利用のしやすさ_5 16 7 1 16

このままの表が表示できない環境のであれば、グラフ化する必要があります。

plt.figure(figsize=(12,8))
sns.clustermap(df_happy_clust_dmy_gp_g,cmap="viridis")
plt.savefig('figure_2.png')
plt.show()

f:id:yyhhyy:20190707224604p:plain

ただ、無理にグラフにするより、スプレッドシートに吐き出して条件付き書式で見た方が楽だとは思います

df_happy_clust_dmy_gp_g.to_csv("df_happy_clust_dmy_gp_g.csv")

クラスター分析で各クラスターについての説明を考えるのはいつも恣意的ですが、例えば以下のように分類できます。

  • クラスター「0」はどの項目にも高い評価で満足し、住宅も高くないと感じているようです。
  • クラスター「1」は学校の質は低いと思っていますが、市には満足をしている。
  • クラスター「2」は警察の質は低いと思っていますが、市には満足しています。
  • クラスター「3」は住宅は高いと感じていますが、市には満足している。

それぞれ幸福度は高いわけですから、評価が低い箇所についてはあまり気にしていないということと同義だと判断しました。もちろんこれ以外の着眼点で分けるのも良いでしょう。

決定木分析

最後に、仮に今後新しく取得されたアンケートデータから、各クラスターに回答者を振り分けたい。となれば、様々な機械学習手法が可能なのですが、一般的な人にも理解して貰いやすい、という意味では、決定木分析が良いのではないかと思います。

先ず先程のデータをラベルとデータとにわけ、どちらもnumpyの配列にします

y = np.array(df_happy_clust["cluster_ID"].values)
display(y)
array([3, 3, 3, 0, 1, 0, 2, 2, 1, 1, 1, 1, 1, 2, 1, 1, 2, 3, 3, 2, 1, 2,
       1, 3, 3, 0, 2, 2, 2, 0, 0, 1, 1, 3, 2, 1, 0, 0, 2, 2, 0, 3, 1, 1,
       1, 0, 0, 3, 1, 1, 2, 0, 3, 0, 0, 3, 3, 3, 3, 3, 3, 2, 0, 3, 2, 1,
       0, 0, 3, 2, 1, 0, 1, 2, 2, 1, 3], dtype=int64)
X = df_happy_clust.drop("cluster_ID",axis=1).values
display(X)
array([[5, 3, 3, 3, 3, 5],
       [5, 5, 3, 5, 5, 5],
       [5, 4, 4, 4, 4, 5],
       [5, 2, 4, 5, 5, 5],
       [3, 2, 4, 3, 4, 4],
       [5, 1, 4, 3, 4, 5],
       [4, 2, 3, 3, 4, 4],
       [4, 2, 3, 3, 4, 4],
       [5, 1, 2, 5, 2, 4],
       [4, 3, 3, 3, 3, 4],
       [3, 3, 3, 5, 5, 5],
       [3, 3, 1, 3, 3, 4],
       [3, 3, 1, 3, 3, 4],
       [5, 3, 3, 3, 5, 3],
       [3, 2, 4, 4, 4, 5],
       [3, 2, 4, 4, 4, 5],
       [4, 1, 3, 1, 1, 4],
       [5, 3, 4, 3, 4, 5],
       [5, 3, 4, 3, 4, 5],
       [5, 2, 3, 3, 2, 5],
       [4, 4, 3, 4, 2, 4],
       [4, 2, 4, 3, 2, 4],
       [3, 1, 2, 4, 3, 5],
       [5, 3, 4, 4, 4, 5],
       [5, 3, 3, 4, 4, 5],
       [5, 2, 5, 5, 5, 3],
       [5, 1, 3, 3, 4, 4],
       [5, 1, 3, 3, 4, 4],
       [5, 1, 3, 3, 4, 4],
       [5, 2, 4, 3, 4, 5],
       [5, 2, 4, 3, 4, 5],
       [4, 3, 2, 4, 3, 4],
       [4, 3, 2, 4, 3, 4],
       [5, 3, 5, 5, 4, 5],
       [5, 2, 4, 2, 2, 4],
       [4, 4, 3, 3, 2, 5],
       [5, 2, 4, 4, 5, 5],
       [5, 2, 4, 4, 5, 5],
       [4, 1, 3, 4, 4, 4],
       [4, 1, 3, 4, 4, 4],
       [5, 1, 5, 5, 5, 5],
       [5, 4, 5, 5, 5, 5],
       [4, 2, 2, 4, 4, 5],
       [4, 3, 3, 4, 3, 4],
       [4, 3, 3, 4, 2, 4],
       [5, 2, 5, 5, 5, 5],
       [5, 2, 3, 5, 5, 5],
       [5, 5, 5, 5, 5, 5],
       [4, 3, 2, 4, 4, 4],
       [4, 2, 4, 4, 4, 4],
       [5, 3, 2, 4, 4, 4],
       [5, 2, 3, 4, 4, 5],
       [5, 3, 3, 3, 5, 5],
       [5, 2, 4, 5, 4, 5],
       [5, 1, 3, 4, 5, 5],
       [5, 4, 5, 5, 5, 4],
       [5, 5, 3, 4, 4, 5],
       [5, 4, 4, 3, 3, 4],
       [5, 4, 4, 4, 4, 5],
       [5, 5, 5, 5, 5, 5],
       [5, 4, 5, 4, 5, 4],
       [4, 2, 3, 4, 3, 3],
       [5, 2, 2, 4, 5, 5],
       [5, 3, 2, 4, 4, 5],
       [5, 3, 4, 3, 4, 3],
       [4, 3, 4, 4, 3, 4],
       [5, 1, 4, 3, 5, 5],
       [5, 1, 5, 3, 5, 5],
       [5, 3, 4, 4, 4, 4],
       [5, 2, 4, 4, 2, 3],
       [3, 4, 4, 5, 1, 3],
       [5, 1, 5, 5, 5, 5],
       [4, 3, 3, 4, 4, 4],
       [5, 5, 1, 1, 5, 1],
       [5, 2, 3, 4, 4, 3],
       [5, 2, 3, 4, 2, 5],
       [5, 3, 3, 4, 4, 5]], dtype=object)
from sklearn import tree

決定木は永遠に細かくできます。仮の数字で段階を決めてしまいます。

dtree = tree.DecisionTreeClassifier(max_depth=4)
dtree = dtree.fit(X,y)

作ったモデルがどの程度の精度なのか?確認してみます

dtree_pred = dtree.predict(X)
display(dtree_pred)
array([3, 3, 3, 0, 1, 0, 2, 2, 0, 1, 1, 1, 1, 2, 1, 1, 2, 3, 3, 1, 1, 2,
       1, 3, 3, 0, 2, 2, 2, 0, 0, 1, 1, 3, 2, 1, 0, 0, 2, 2, 0, 3, 1, 1,
       1, 0, 0, 3, 1, 2, 2, 0, 3, 0, 0, 3, 3, 3, 3, 3, 3, 2, 0, 2, 2, 1,
       0, 0, 3, 2, 1, 0, 1, 2, 2, 1, 3], dtype=int64)

ラベルとどれくらいあっているか?正直、過学習してそうですが、今は面倒なのでこのまま進めます

sum(dtree_pred == y) / len(y)
0.948051948051948

jupyter notebookならgraphvizパッケージを使って可視化が可能です7

import pydotplus
from IPython.display import Image
from graphviz import Digraph

配列にした際に特徴量の名称が消えてしまっているので、列名を代入し、またクラスター名も数字のままだとエラーになるため、ラベルデータはastypeを使って文字列(string)に変更しておきます

dot_data = tree.export_graphviz(dtree,out_file=None,
                                feature_names = df_happy_clust.columns[0:-1],
                                class_names = y.astype("str"))

日本語を出力するにあたっては随分苦労しましたがこのブログが正解のようです。

mk-55.hatenablog.com

graph = pydotplus.graph_from_dot_data(dot_data)

graph.set_fontname('Noto Sans CJK JP')
for node in graph.get_nodes():
    node.set_fontname('Noto Sans CJK JP')
for e in graph.get_edges():
    e.set_fontname('Noto Sans CJK JP')

graph.write_png("dtree.png")
Image(graph.create_png())

f:id:yyhhyy:20190707224750p:plain

例えば、「行政サービス情報のアクセスしやすさ」が4.5以上つまり5で、「住宅供給の高さ」が2.5以上だとクラスター3だと判断されるわけですが、実際のデータと比べても特徴をよく表しています。

df_happy_clust_dmy[df_happy_clust_dmy["cluster_ID"]==3].sum().head(10)
cluster_ID              60.0
行政サービス情報へのアクセスしやすさ_3     0.0
行政サービス情報へのアクセスしやすさ_4     0.0
行政サービス情報へのアクセスしやすさ_5    20.0
住宅供給の高さ_1                0.0
住宅供給の高さ_2                0.0
住宅供給の高さ_3               10.0
住宅供給の高さ_4                6.0
住宅供給の高さ_5                4.0
公立の学校の全般的な質の良さ_1         0.0
dtype: float64

f:id:yyhhyy:20190707224827p:plain

分析からわかること

例えば、市であれば自分たちが改善できることできなことがあるはずで、対応できることを対応する場合、どのクラスタの市民に残って欲しいか?増えて欲しいか?を絞り込むことができます。

また、重点ターゲットとなるクラスタを決めたら、決定木のジャッジのポイントを元に訴求項目を作れば、どういう謳い文句が必要か?ということが自ずと決まってきます。

今回はアンケート項目が少なかったので行えませんでしたが、通常はクラスタ分析に使った質問以外にもっと多くの質問をしているはずなので、その質問で決定木分析を行うと

  • どんなメディアに接触しているのか?
  • どんな世代に多いのか?
  • どんなことに興味関心があるのか?

など、広告でターゲティングしやすいセグメントを元に決定木分析を行うこともできるでしょう


決定木以外にも通常のグラフで日本語が文字化けします。seabornを使う場合は、最初にフォントを指定しておけば大丈夫のようで、今のところ私もこれでうまく行っています

qiita.com


最近、Pythonの本がすごくわかりやすいものが増えてきました。今まではエンジニア出身の人の本が多かったですが、徐々にデータサイエンス側に向けて適切な分量の本が増えています

コンパクトにまとまっているけれどそれでも応用範囲が広い

こちらはドリルのような本でも知っていると便利

Pythonではないもの、データサイエンスにとって必要なものは何か?を丁寧に解説していて、この本は最初の頃に読んでおきたかったと思う。


  1. ここのサイト、機械学習に特化しているので、どういう目的で?どういうタイプのデータが欲しいか?をチェックしていくとデータの候補がセレクトされるんですよ!凄いですね。

  2. 因みにこのファイル文字コードが特殊でうまく読み込めず、一度Googleスプレッドシートで読み込んで再保存しています。

  3. 「cost_of_housing」と5段階スケールで聞いてますが、これは5の方が価格が高い、という意味ですかね?

  4. 片寄っていることの確認のため簡便にdescribeを使っていますが、そもそも間隔尺度は通常は平均などは意味がないので出しません。

  5. 学校のテストの偏差値と根本的な考えは同じです。みんなの点数が高かった数学の90点と、みんなの点数が低かった英語の90点とは平等に足すことはできません。

  6. drop_first=True」は今回は指定しない。

  7. Win10ではgraphivzのdot.exeをシステム環境PATHで通す必要があるようです。

機械学習屋と統計屋とデータサイエンティストの違いとは

統計学は最強の学問である」という本がブームになって既に何年も過ぎましたが、未だに「データサイエンティスト」界隈はポジショントークに溢れているようで、データサイエンスについて本を読もうとすると、データサイエンスブームに乗って我田引水しようとする2系統の間で彷徨わされてしまいます。

統計屋の場合

統計学出身の人が先ず進めるのは以下の本ではないでしょうか?

これはこれでとても勉強になる本なのですが、「早くデータサイエンスでビジネスに必要な判断をしてみたい!」というニーズを抱えている人にとっては随分と遠回りに感じます。

「同じことを100回試してみたら95回は”AとBに違いがある”と出るであろう。実際に違いがない可能性もあるし、違いがあるとは断言できないが可能性はとても高い」

学問的なことを言えばその通りかもしれませんが、ビジネスでは「意思決定」をしないといけないので、こんな歯切れの悪い言い方にこだわられても「じゃぁどうすりゃいいんだよ?」と怒られるだけでしょう。1

有意 - Wikipedia

機械学習屋の場合(あるいはAI人材)

一方でプログラミング系の人が勧めるのはこちらでしょう

これもこれで、機械学習(とくにディープラーニング)について知りたい人にとってとても分かり易い素晴らしい本ですが、これを読んだからといって何かの「意思決定」ができるようにはなりません。

「画像から猫か猫以外を分類できました!」

という報告を貰ったところで「で?新しい発見はないの?うちは新たに何に着目したら競合に勝てるの?」と怒られるだけでしょう。

勿論、自社サービスから離反しそうな人をディープラーニングで判断してクーポンを送りつける、といったことはできるでしょうが、ディープラーニングの最大の欠点は ブラックボックス になってしまうことで、自社サービスから離反しそうになった原因をつきとめることができないので、改善策は出てきません。

データを元に意思決定したい人は「データマイニング」から入った方が良いかもしれない

今の「データサイエンティストブーム」のターニングポイントは2つあって、1つは「ビッグデータブーム」で、もう一つは「ディープラーニングブーム」です。この2つのバズワードに乗っかるように「データサイエンス」という言葉が広まり始めました。

しかし「データを元に何かを発見し意思決定する」という行為自体は、データサイエンティストブームの前から存在していて、それは「データマイニング」と呼ばれていました。

どうしてもバズワードが生まれてしまうと乗っかる方が本も売れますし名前も売れますので、戦略的にバズワードに乗っかる人が出てしまいます。 [^1] その結果「データサイエンス」という言葉を見て本を手に取ったときに、「やたらp値にこだわる統計屋」と「何でもディープラーニングで済ませようとする機械学習屋」 の本にあたってしまい、求めていた情報と違うものに出会って期待外れを感じてしまうのです。

データマイニング」目的で「データサイエンス」を学びたい人に

データマイニング目的で学ぶべき統計学機械学習のポイントは、「起きている現象について分析し意思決定できるか?」という点です。その意味では、全ての統計学機械学習を学ぶ必要はなく、以下のようなラインナップがあれば充分です。

要因のうち何が重要かわかる

  • 線形回帰分析
  • 重回帰分析
  • ロジスティック回帰分析
  • 決定木

データの把握をしやくする

  • EDA(探索的データ分析)/或いは、記述統計
  • 因子分析
  • k平均法によるクラスター分析
  • コレスポンデンス分析

A/Bテストのパターンを決める

  • 実験計画法

どの程度の効果があるか?を探る

時系列の取り扱い

  • 状態空間モデル
  • 移動平均・自己回帰(ARIMA)

お勧めの本

上記のようなポイントを抑えて優先順位をつけて学びたい場合、以下のような本が良いでしょう。

ライトな本

先ず「線形回帰分析」ができるだけでもエクセルの棒グラフ地獄から抜け出せます。もちろん、どこかで次の一歩としてロジスティック回帰分析が必要になってきますが、最初はこれだけでも良いでしょう。

一通りしっかり知りたい人はこの本が今のところとてもバランスよくまとまっています。一つ一つの説明は少ないので、細かいことは別の本を読まないといけないかもしれませんが、自分にとって何が重要で何が重要でないか?の指南はこの本に頼れば良いと思います。とても良い本です。


  1. 定義が曖昧な「人工知能」などという言葉を学者は使うべきではないと思いますが、もはやその方が通りが良いので、機械学習の人も人工知能強化学習の人も人工知能と言ってしまいがち、というのも同じ現象です。

Pamdasで特定列が条件を満たすデータをグループ単位で抜き出す

Webサイトのアクセスログを分析していると

コンバージョンに至ったユニークブラウザだけ抜き出してアクセス履歴を抜き出したい

ということがあると思います。

これが意外と面倒だったので備忘録的にまとめておきたいと思います。

マーケティングオートメーションのメールからコンテンツにアクセスした人のリストを抜き出す

「~/hogehoge_1/」というページにアクセスしたデータを抜き出すにはログデータの「URL」列に対して文字列検索をします。例えば、str.containsを使うのが良いでしょう。

hoge_list = df["url"].str.contains("hogehoge_1/')

更に、マーケティングオートメーションツール(例えばマルケト等)で配信したメールのリンクURLに「mkt」というパラメータを振り振っていたとして、Marketo経由で(=条件1)かつ「~/hogehoge_1/」にアクセスした(条件=2)ような行(レコード)を抜き出すには、「&」でつなげます。

hoge_list = df["url"].str.contains("mkt") & df["url"].str.contains("hogehoge_1/")

但し、このままでは「~/hogehoge_1/」にマルケトメール経由でアクセスした履歴が終えるだけで、他にどのページを見たか?という情報が消えてしまいます。

そこで、「~/hogehoge_1/」にアクセスした「IPアドレスの一覧」1を一旦取得し、リスト化します。

hoge_ip_list = df[hoge_list]["ip"]

コンバージョンしたリストを元に元のデータフレームにラベルをはる

普通にこのIPアドレスリストを使ってデータフレームから抜き出しても良いのですが、コンバージョンの種類が1つではないなど、今後の汎用性を考えると、元のデータフレームにラベルを付与する方が良いでしょう。

一度データフレーム化し、インデックスをリセットしています。そうしないとラベルが正しく作成できません。

hoge_ip_list_df = pd.DataFrame(lhoge_ip_list)
hoge_ip_list_df = hoge_ip_list_df.reset_index(drop=True)

コンバージョンしたIPアドレスリストの個数を「len(~.index)」を使って計測し、その数だけ「hogehoge_cv」 という文字列を作成。

hoge_label = pd.DataFrame({"hoge_cv":["hogehoge_cv"] * len(hoge_ip_list_df.index)},columns=["hoge_cv"])

hoge_cv」という列名でコンバージョンしたIPアドレスリストのデータフレームに列を追加します。

hoge_ip_list_df["hoge_cv"] = hoge_label

コンバージョンしたIPアドレスリストのデータフレームと元のデータフレームをleeft_joinすれば、条件を満たすものだけが残ります。

df_mkt = pd.merge(hoge_ip_list_df,df,how="left")

CSVファイル化すればExcelユーザーに共有可能です。

df_mkt.to_csv("df_mkt_all.csv")
同様に別のコンバージョンもラベル化した場合

他のCVのラベルも作った場合はそれらのIPアドレス一覧データフレームをmergeしていけば良いでしょう。

hoge_ip_list_1_2_df = pd.merge(hoge_ip_list_df,hoge_ip_list_2_df,how="outer")

メールで配信したコンテンツの閲覧だけ欲しいとき

全ページではなく、メールで配信したコンテンツの閲覧状況だけに更に絞りたい場合。 先程のデータの中から配信コンテンツをマーケティングオートメーションのメール経由で閲覧したレコードだけ集計してconcatで縦に結合していけば絞れます。

df_mrkt_mail2 = df_mkt[df_mkt["url"].str.contains("mkt") & df_mkt["url"].str.contains("hoge_mail_2/")]
df_mrkt_mail3 = df_mkt[df_mkt["url"].str.contains("mkt") & df_mkt["url"].str.contains("hoge_mail_3/")]

df_mkt_mails = pd.concat([df_mkt_mail2,df_mkt_mail3], ignore_index=False)

パラメータなしのURLが欲しいとき

パラメータ付きURL以外にパラメータなしのURLが欲しいときは、str.splitを使い「?」で分割し、その1つ目の要素(つまりPython的には0番目)を抜き出して、新しい列「path」として插入します。

df_mkt["path"] = df_mkt["url"].str.split("?").str.get(0)

  1. 別にユニークブラウザでも構いません。BtoBを意識してIPアドレスにしてみました。

データ分析で重要なのは「視点」~データサイエンティスト達に本当に必要だったと思われるものについて~

「統計のことがわかっていない自称データサイエンティストが増えて現場は困っている」という”データサイエンティスト”のつぶやきをよく見るようになりました。

これはいったいどういうことなのか?を考えてみたいと思います。

データサイエンティストブーム

データサイエンティストという職業が注目されたのは2010年頃だったでしょうか?

この本が出たのはもう5年前。この頃は「データサイエンティストが不足する」などと言われ、大学でも統計学を教えるべきではないか?などの議論があったように思いますが、その後、データサイエンティストブームは沈静化したように思います。

実態としては、「そこまで”データサイエンティスト”が必要ではなかった」ということではないかと思いますし、必要ではなかったと気がつき始めたからこそ”データサイエンティスト”として生き残ろうと思っていた人は、「方向転換のサンクコストを払うべきか?」、「自分たちの存在価値を喧伝するべきか?」の岐路に立たされているのではないか?と思います。

それでは、実際に求められていて、今も足りない人材は何だったのでしょうか?

データ分析に必要な人材

データ分析に必要なスキルは主に以下の3つです

  • ビジネスの視点
  • 統計学の視点
  • システムの視点

ビジネスの視点は「マーケター」の役割で、システムの視点は「エンジニア」の役割です。

近年の急激なマシンスペック向上により、これまでマーケターが使っていた統計学レベルでは対応できない統計分析ができるようになってきたので、改めて「統計学」の視点が必要となったのでした。

そして、”これからはビッグデータだ”と言っていた頃は、まさに急激に増えた”機械学習”、”ディープラーニング”といった能力が必要とされていたのです。

※なおデータ分析に必要な人材とチームのイメージについてはこの本が最もバランスが良いと思います

必要とされる統計学にはバラつきがあった

一方で、求められる「統計学」についてはバラつきがありました。

Googleの中でアルゴリズムを作りディープラーニングのライブラリを開発しなければならない人、音声チャットサービスの原型を作る人、画像処理で人を認識するエンジンを作らないといけない人、といった専門人材は、数学的な理解が必要で、大学で工学科を出てないと難しいと思います。

一方で、誰かが作ったライブラリを使って自社のサービスに利用する、というレベルになってくると、専門人材である必要はなく、 エンジニアサイドの人が勉強すれば事足りる ということになると思います。

例えば、「自社にプライベートDMPを導入してデジタルマーケティングを強化したい」という目的であれば、マーケティングオートメーションツール、タグ管理ツール、コンテンツマネジメントシステム等々、全て外部のサービスをどう組み合わせるか?であって、トレジャーデータなどになってくると主要なツールとは簡単に連携できるように準備されています。このレベルであれば、情報システム部門の人とマーケティング部の人が連携すれば事足りてしまいます。

転機が明確になったのは、「GoogleがTensorflowをオープンにし、IBMがWatsonの門戸を広げたとき」だったと思います。

全てのアルゴリズムを自社で作り上げる必要があるのであれば、「統計学」に深い造詣のある人材が各企業に必要となりますが、それがサービスとして提供されるのであれば、”何が出来て何が出来ないかわかっていれば良い”ということになります。

一方で求められているのはビジネス視点

寧ろ不足しているのは、”統計学の素養があるマーケター”の方です。

多くの企業でマーケティング部に居る人は文系社員で、統計学についてはせいぜい大学時代にSPSSSASを使っていました、というレベルであり、使うのは専ら因子分析や線形回帰分析などです。勿論、彼等がエンジニアのように引き続き勉強してくれたら物事は解決するのですが、元が文系だと改めて学習するのが辛いこと、普段使い慣れていないコードの入力が必要なこと、などで思ったように勉強が難しいのではないかと思います。

冒頭の「”データサイエンティスト”達から文句を言われている”自称データサイエンティスト”」とは、恐らくこの手の見よう見マネで統計学を復習しはじめた文系マーケターのことではないかと思います。

他方、”データサイエンティスト”側は、エンジニア上がりや学者上がりのことが多く、「ビジネスの視点」というのが不足しています。

「A/Bテストをして良いサイトデザインやバナークリエイティブを残す」などの試みはやって頂いて構わないのですが、「それはビジネスに対してどの程度利益が出るのか?」という視点が抜け落ちています。もし、この辺りの勘所を”データサイエンティスト”たちが鍛えられれば、マーケターという職種に鞍替えできるのではないか?と思います。

”データサイエンティスト”が”自称データサイエンティスト”に仕事を奪われているのは、アサインする側の経営者サイドと会話が成立していないからです。彼等の視点は常に”ビジネスへのインパクト”です。全ての業務をそこから逆算して考えるクセがついていないと会話になりません。

例えば、この本にあるように

「月額6000円でマニアだけが加入するもの」だったADSLを「生産ロットを増やしてモデムが安くできるなら、月額2830円で100万人が入ってペイするもの」にしてしまえばいい!

と発想ができるかどうか? です。

データを取得するのにかかるコストを弾くのはエンジニアに頼むとして、そのデータを取得することでどのようなリターンが見込めるか?を考えて提示する能力が求められています。

幸いにも、工学科出身でビジネスコンサルタントをやっている人は世の中に大量に居ますので、”データサイエンティスト”側からマーケター側への転向は、マーケターが統計学のスキルを強化するよりも早くたどり着けるのではないか?

と文系でマーケティング業界の端っこに居る私は思います。


参考図書

大局的な視点を得る為に良い本

幾つかの企業のビジネスモデルが分析されていますが、特に、 外側からはわからない「コスト構造」までインタビューしている 本は貴重です。

タップスと言えば人工知能や宇宙開発という変わった会社ですが、佐藤さんの先を読む力は飛び抜けていると思います。本書では、仮想通貨ブームを人々の経済の在り方の視点にまで広げて見ています。

視点が大事、という意味ではやはりこの本が一番でしょうか。書いてある通りに発想するのが難しいんですよ、とは思うものの、出発点の視点がとにかく大事です。

Rで時系列データをクール単位で集計する

過去何度か時系列のエントリーで「xts」が便利だと書いて居たのですが、日足データを月次にすることはできてもクール単位に集計する便利な関数がない、ということに気が付き色々調べました。

事例としてJEITAの「民生用電子機器国内出荷統計」を使ってみます。

ダウンロード先はこちら http://www.jeita.or.jp/japanese/stat/shipment/

データのクオーター化

データについて

上記サイトから集めたデータは以下のような形をしています。

df <- read.csv("data.csv",header=T)
head(df)

このように右向きに時系列のデータが並んでいます。

              商品名     種別 X1月 X2月 X3月
1         薄型テレビ     合計  416  579  508
2         薄型テレビ 29型以下  103  147  180
3         薄型テレビ 3036144  207  160
4         薄型テレビ 3749114  157  112
5         薄型テレビ 50型以上   54   69   56
6 ブルーレイディスク     合計  216  272  272
7 ブルーレイディスク レコーダ  169  224  208
8 ブルーレイディスク プレーヤ   47   48   64

データの抜粋

このうち薄型テレビの合計だけを取り出します。 列名でデータを抜き出す為に「dplyr」を使います。

library("dplyr", lib.loc="C:/hogehoge/R-3.4.1/library")
dat <- df %>%
  dplyr::filter(商品名 == "薄型テレビ" & 種別 == "合計")
head(dat)
      商品名 種別 X1月 X2月 X3月 X4月 X5月 X6月
1 薄型テレビ 合計  416  579  508

ts型にするために

時系列にするためts型にしますが、先ず横向きに並んだデータでは上手く認識しませんので転置します。

head(t(dat))
       [,1]        
商品名 "薄型テレビ"
種別   "合計"      
X1月   "416"       
X2月   "579"       
X3月   "508"       
X4月   "309"  

最初の2行は数値ではないので、3行目からのデータを使います。

head(t(dat[1,3:ncol(dat)]))
       1
X1月 416
X2月 579
X3月 508
X4月 309
X5月 345
X6月 456

ts型にする

今回のデータは2014年1月からの月次データでしたので、スタートは2014,1、頻度は12ということになります。

dat_t <- ts(t(dat[1,3:ncol(dat)]),start = c(2014,1),frequency = 12)
print(dat_t)

自動的に月名が割り当てられます。 なおデータフレームではないので「print」でデータ確認します。

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct
2014 416 579 508 309 345 456 395 305 515 337
2015 392 445 627 339 322 421 364 330 353 371
2016 348 381 510 334 350 332 374 325 374 352
2017 305 362 498 337 329 355 290 297 357 

ts型のグラフはts.plotで確認できます

ts.plot(dat_t)

f:id:yyhhyy:20171112175845p:plain

四半期(クオーター)データにする

「aggregate」関数を使います。 四半期データは3ヶ月ごとのため、 12ヶ月単位のデータを4分割することになります。

また、合計値にしたいためファンクションは「sum」を選びます。

dat_t_q <- aggregate(dat_t,nfrequency=4,FUN=sum)
print(dat_t_q)

クオーター単位に合計されています。

     Qtr1 Qtr2 Qtr3 Qtr4
2014 1503 1110 1215 1662
2015 1464 1082 1047 1529
2016 1239 1016 1073 1420
2017 1165 1021  944  

四半期データのグラフになりました。

f:id:yyhhyy:20171112180012p:plain

データフレームにする

データフレームにすると残念ながらts型で割り振ってくれた名称が消えてしまうため、データフレームにするには、時系列パートの列名が必要になる。 1

列名の元となるデータを用意する

dmn_1 <- c("2014年","2015年","2016年","2017年")
dmn_2 <- c("1Q","2Q","3Q","4Q")

組み合わせて結合する

for分を使って上記の組み合わせを作ります。 空のデータフレームを用意してrbindしていくことと、stringrパッケージを使って文字列を結合することがポイントです。

library("stringr", lib.loc="C:/hogehoge/R-3.4.1/library")
dmn <- as.data.frame(NULL)
for(i in 1:4){
  for(j in 1:4){
    dat <- str_c(dmn_1[i],dmn_2[j],sep='', collapse=NULL)
    dat <- as.data.frame(dat)
    dmn <- rbind(dmn,dat)
    dat <- NULL
  }
}
head(dmn)
            dat
1 2014年1Q
2 2014年2Q
3 2014年3Q
4 2014年4Q
5 2015年1Q
6 2015年2Q

ただし、今回のデータは2016年の3Qで終了ですので、数を絞る必要があります。 一度転置させないといけません。

dmn <- t(dmn)
TIME_dmn <- dmn[1:14]
head(TIME_dmn)
[1] "2014年1Q" "2014年2Q"
[3] "2014年3Q" "2014年4Q"
[5] "2015年1Q" "2015年2Q"

列名を使ってデータをデータフレーム化

先程用意した列名を使いながら、時系列データをデータフレーム化します。

dat_q <- data.frame(TIME=TIME_dmn,薄型テレビ_合計 = c(dat_t_q))
head(dat_q)
           TIME 薄型テレビ_合計
1 2014年1Q            1503
2 2014年2Q            1110
3 2014年3Q            1215
4 2014年4Q            1662
5 2015年1Q            1464
6 2015年2Q            1082

同様のことを複数データに行う

クオーター化を関数化

一連のデータフレームへの変化までを関数にします。

※列名は関数外です。

tv_goukei_func <- function(x){
  dat <- df %>%
    dplyr::filter(商品名 == "薄型テレビ" & 種別 == "合計")
  dat_t <- ts(t(dat[1,3:ncol(dat)]),start = c(2014,1),frequency = 12)
  dat_t_q <- aggregate(dat_t,nfrequency=4,FUN=sum)
  dat_q <- data.frame(TIME=TIME_dmn,薄型テレビ_合計 = c(dat_t_q))
  return(dat_q)
}

同様に、サイズ別のデータ用の関数も作ります

tv_goukei_func <- function(x){
  dat <- x %>%
    dplyr::filter(商品名 == "薄型テレビ" & 種別 == "合計")
  dat_t <- ts(t(dat[1,3:ncol(dat)]),start = c(2014,1),frequency = 12)
  dat_t_q <- aggregate(dat_t,nfrequency=4,FUN=sum)
  dat_q <- data.frame(TIME=TIME_dmn,薄型テレビ_合計 = c(dat_t_q))
  return(dat_q)
}

tv_29_func <- function(x){
  dat <- x %>%
    dplyr::filter(商品名 == "薄型テレビ" & 種別 == "29型以下")
  dat_t <- ts(t(dat[1,3:ncol(dat)]),start = c(2014,1),frequency = 12)
  dat_t_q <- aggregate(dat_t,nfrequency=4,FUN=sum)
  dat_q <- data.frame(TIME=TIME_dmn,薄型テレビ_29型以下 = c(dat_t_q))
  return(dat_q)
}

tv_30_func <- function(x){
  dat <- x %>%
    dplyr::filter(商品名 == "薄型テレビ" & 種別 == "30~36型")
  dat_t <- ts(t(dat[1,3:ncol(dat)]),start = c(2014,1),frequency = 12)
  dat_t_q <- aggregate(dat_t,nfrequency=4,FUN=sum)
  dat_q <- data.frame(TIME=TIME_dmn,薄型テレビ_30_36型 = c(dat_t_q))
  return(dat_q)
}

tv_37_func <- function(x){
  dat <- x %>%
    dplyr::filter(商品名 == "薄型テレビ" & 種別 == "37~49型")
  dat_t <- ts(t(dat[1,3:ncol(dat)]),start = c(2014,1),frequency = 12)
  dat_t_q <- aggregate(dat_t,nfrequency=4,FUN=sum)
  dat_q <- data.frame(TIME=TIME_dmn,薄型テレビ_37_49型 = c(dat_t_q))
  return(dat_q)
}

tv_50_func <- function(x){
  dat <- x %>%
    dplyr::filter(商品名 == "薄型テレビ" & 種別 == "50型以上")
  dat_t <- ts(t(dat[1,3:ncol(dat)]),start = c(2014,1),frequency = 12)
  dat_t_q <- aggregate(dat_t,nfrequency=4,FUN=sum)
  dat_q <- data.frame(TIME=TIME_dmn,薄型テレビ_50型以上 = c(dat_t_q))
  return(dat_q)
}

関数を順番に適応

これらの関数をデータに適応しデータフレームに結合します。

df_tv_goukei <- tv_goukei_func(df)
df_tv_29 <- tv_29_func(df)
df_tv_30 <- tv_30_func(df)
df_tv_37 <- tv_37_func(df)
df_tv_50 <- tv_50_func(df)

df_q <- dplyr::full_join(df_tv_goukei,df_tv_29)
df_q <- dplyr::full_join(df_q,df_tv_30)
df_q <- dplyr::full_join(df_q,df_tv_37)
df_q <- dplyr::full_join(df_q,df_tv_50)
head(df_q)
           TIME 薄型テレビ_合計 薄型テレビ_29型以下 薄型テレビ_30_36型
1 2014年1Q            1503                   430                    511
2 2014年2Q            1110                   336                    323
3 2014年3Q            1215                   321                    411
4 2014年4Q            1662                   457                    524
5 2015年1Q            1464                   536                    447
6 2015年2Q            1082                   330                    328

グラフ化

ggplot2を使います

library("reshape2", lib.loc="C:/hogehoge/library")
library("ggplot2", lib.loc="C:/hogehoge/library")

日本語の列名は時にggplotで時系列の順番が狂いますので、reorderを使って並べ替えます。 そのために番号を先ず振ります。

df_q$num <- rep(1:nrow(df_q))

今回はfacet_wrapで種類別にし、凡例の位置をtiopへ、また塗りつぶしの色を手動で指定しています。

df_q_m <- melt(df_q,id.vars = c("TIME","num"))
head(df_q_m)
g <- ggplot(df_q_m,aes(x=reorder(TIME,num),y=value,group=variable,fill=variable))
g <- g + theme_classic()
g <- g + geom_bar(stat = "identity",position = "stack")
g <- g + facet_grid(variable~.)
g <- g + geom_text(aes(x=reorder(TIME,num),y=value,label=value),vjust=-0.5)
g <- g + scale_y_continuous(labels = scales::comma,limits =c(0,2200))
g <- g + ylab("")
g <- g + xlab("")
g <- g + ggtitle("薄型テレビ量販店四半期販売台数数位")
g <- g + scale_fill_manual(values = c("brown4",
                           "deepskyblue4",
                           "deepskyblue3",
                           "deepskyblue2",
                           "deepskyblue1"))
g <- g + theme(legend.position = "top")
g <- g + theme(axis.text.x = element_text(angle = 45, hjust = 1,size=6))
plot(g)

グラフを保存します

ggsave(plot=g,filename = "g-3.png",scale=2,height = 2.5,width = 5)

f:id:yyhhyy:20171112190139p:plain


  1. ほんとはもっとスマートな手法があると良いのですが。。。