#!/usr/bin/env python
# coding: utf-8

# # 食べログ3.8問題に終止符を打つ

# 「評価3.8以上は年会費を払わなければ3.6に下げられる」という問題が最近バズっていますが、  
# 
# * 食べログのスコアがどのように決まっているのか
# * 食べログは年会費を払わないとスコアを下げているのか
# 
# の2点が混同されているなぁと感じました。  
# そこで、因果推論の手法を用いて  
# 
# * 年会費を払う=有料会員になることは、お店のスコアをどれくらいあげるのか
# 
# を明らかにします。

# 参考(本家): https://clean-copy-of-onenote.hatenablog.com/entry/tabelog38_problem  
# 参考(nardtreeさんのEDA): https://gist.github.com/GINK03/8826e84c2c7b8ccb4c26158f17d2eb4c  
# 参考(考察): https://news.yahoo.co.jp/byline/inoueakito/20191012-00146614/

# 食べログは有料会員の中でもランクがあったり、非会員と無料会員などがありますが、説明の簡便性のため有料会員/無料会員に2分します

# # スクレイピング周り

# 東京都のお店をエリア別に計15476店舗取得しました。  
# コードは一番下に載せてあります。

# # データの読み込みと前処理

# In[178]:


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import json
import glob
import math
from pathlib import Path
from collections import Counter
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import StratifiedKFold
import lightgbm
import japanize_matplotlib
from scipy import stats

sns.set_style('darkgrid')
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', 500)
plt.rcParams['font.family'] = 'IPAexGothic'

get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('config', "InlineBackend.figure_formats = {'png', 'retina'}")


# In[117]:


def force_to_int(x):
    try:
        return int(x)
    except:
        return -1
objs = []
for fn in glob.glob('../input/*'):
    try:
        arr = json.loads(open(fn).read())
    except Exception as exc:
        print(exc)
        continue
    for obj in arr:
        objs.append(obj)

df = pd.DataFrame(objs)
df.drop_duplicates(subset=['mise'], keep='last', inplace=True)
df['score'] = df.score.apply(float)
df['save_num'] = df.save_num.apply(force_to_int)
df['review_num'] = df.review_num.apply(force_to_int)
df['idv_score'] = df.idv_score.apply(lambda x: [float(i) for i in x if i != '-'])
tmp = df.drop(['comments', 'idv_score', 'genres'], axis=1)
# レビュー毎のスコアが3.5を超えたか,3.0を下回ったかを特徴量に加える
tmp['idv_over35'] = df['idv_score'].apply(lambda x: sum([1 for i in x if i >= 3.5])/len(x) if len(x)>=1 else 0)
tmp['idv_under30'] = df['idv_score'].apply(lambda x: sum([1 for i in x if i < 3.0])/len(x) if len(x) >= 1 else 0)


# # 基礎統計量の確認

# ## データ数など

# In[118]:


tmp.shape


# In[119]:


tmp.head(3)


# ## スコアの分布

# 全ての店のスコアの分布を確認します。  
# 今までの方々の調査では、一定以上レビューが集まった店の分布において3.6や3.8の手前に山がありました。

# まず全てのお店のスコアの分布を確認します。

# In[120]:


# 全てのお店のスコア
fig, axes = plt.subplots(nrows=2, figsize=(16,8))
sns.distplot(tmp['score'], kde=False, bins=np.arange(3, 4.5, 0.02), ax=axes[0])
sns.distplot(tmp['score'], kde=False, bins=np.arange(3, 4.5, 0.02), ax=axes[1])
axes[1].set_ylim(0, 600)
axes[0].set_title('score hist')
axes[1].set_title('score hist set_ylim(0, 600)')
axes[0].set_xticks(np.arange(3, 5, 0.1))
axes[1].set_xticks(np.arange(3, 5, 0.1))
axes[0].set_title('全てのお店のスコアの分布(0.02刻み)')
axes[1].set_title('全てのお店のスコアの分布(0.02刻み, 高さの上限を600店舗に設定)')
plt.show()


# 次に、一定以上の評価数に絞ります。

# nardtreeさんの分析を参考に、レビュー数が30以上、60以上、90以上でデータを絞ります。

# また、OsciiArtさんを参考に、A以上B未満でbinningするのではなく、Aより大きくB以下でbinningします。  
# これにより、scoreの上限が3.6や3.8に設定されているのかを明らかにします。  
# 
# 参考: https://gist.github.com/GINK03/8826e84c2c7b8ccb4c26158f17d2eb4c

# In[121]:


size30 = tmp.query('review_num>=30').shape[0]
size60 = tmp.query('review_num>=60').shape[0]
size90 = tmp.query('review_num>=90').shape[0]


# In[122]:


fig, axes = plt.subplots(nrows=3, figsize=(16,12))
sns.distplot(tmp.query('review_num>=30')['score'].values+0.00001, kde=False, bins=np.arange(3, 4.5, 0.02), ax=axes[0])
sns.distplot(tmp.query('review_num>=60')['score'].values+0.00001, kde=False, bins=np.arange(3, 4.5, 0.02), ax=axes[1])
sns.distplot(tmp.query('review_num>=90')['score'].values+0.00001, kde=False, bins=np.arange(3, 4.5, 0.02), ax=axes[2])
axes[0].set_title(f'score hist review_num >= 30, size: {size30}')
axes[1].set_title(f'score hist review_num >= 60, size: {size60}')
axes[2].set_title(f'score hist review_num >= 90, size: {size90}')
axes[0].set_xlim(3,5)
axes[1].set_xlim(3,5)
axes[2].set_xlim(3,5)
axes[0].set_xticks(np.arange(3, 5, 0.1))
axes[1].set_xticks(np.arange(3, 5, 0.1))
axes[2].set_xticks(np.arange(3, 5, 0.1))
plt.show()


# スコアを0.02刻みでプロットしました。  
# 今までの議論にあったように、3.6, 3.8の前に壁があります。  
# それだけでなく、[3.68, 3.70)が少なく[3.70, 3.72)が多いなどの偏りが見られます。

# ## その他の基礎統計量

# せっかくなので、スコア以外にも色々データを眺めてみましょう。

# ### 主要店ジャンル

# In[123]:


tmp['genre'].value_counts().head(10)


# ### 主要店ジャンル別有料会員/無料会員別店舗数

# In[124]:


major_genres = set(tmp['genre'].value_counts().head(10).index) - {'その他'}
tmp[tmp['genre'].isin(major_genres)].groupby(['genre', 'is_paid'])[['mise']].count()


# イタリアンや居酒屋、ダイニングバー焼肉は有料店舗が無料店舗の2倍以上あります。    
# 中華料理、寿司、焼き鳥は同じくらいです。  
# カフェとラーメンは無料店舗の方が多いです。

# ### 有料会員/無料会員数で分けたスコア分布

# In[125]:


fig, ax = plt.subplots(figsize=(16,6))
sns.distplot(tmp.query('is_paid==1')['score'], kde=False, bins=np.arange(3, 5, 0.02), label='is_paid=1', ax=ax)
sns.distplot(tmp.query('is_paid==0')['score'], kde=False, bins=np.arange(3, 5, 0.02), label='is_paid=0', ax=ax)
ax.set_title('有料会員/無料会員別スコアの分布(is_paid=1:有料会員, is_paid=0:無料会員)')
ax.set_xticks(np.arange(3, 5, 0.1))
ax.set_xlabel('スコア')
ax.set_ylabel('店舗数')
plt.legend()
plt.show()


# 点数別で有料会員の割合を計算してみます  
# 0.05刻みで、有料会員店舗が100店舗以上あった(3.75,3.8]までをプロットします

# In[126]:


s1 = pd.cut(tmp.query('is_paid==1')['score'], bins=np.arange(3, 3.8, 0.05)).value_counts().sort_index()
s = pd.cut(tmp['score'], bins=np.arange(3, 3.8, 0.05)).value_counts().sort_index()
fig, ax = plt.subplots(figsize=(16, 4))
sns.lineplot(x=np.arange(3, 3.75, 0.05), y=s1/s, ax=ax)
ax.set_title('スコア別で有料会員が全体に占める割合の分布((3.0, 3.05]から(3.65, 3.7]まで0.05刻み)')
ax.set_ylabel('有料会員の割合')
ax.set_xlabel('スコア')
plt.show()


# 店全体にしめる有料会員の割合はスコアが上がるほど高まりますが、(3.55, 3.6]を境に逆に下がっていることがわかります。    
# 有料会員になると露出が増えてレビューが集まりやすくなるので、点数が初期値から真の値に近づきやすくなるのかもしれません。  
# 3.6を超えるには露出を増やすだけでない何か別の要因が必要そうです。

# # 本題: 有料会員になると、お店のスコアはどれくらい上がるのか

# 食べログ3.8問題でみんなが気になっているのは「食べログは店のスコアを有料会員かどうかで操作しているのか」です。  
# なので、お店が有料会員である事によるスコアの上昇、因果効果がどれくらいあるのかを調べます。

# ## 因果効果調査対象のデータ選別

# 食べログは初期スコア(3.04あたり)からスタートし、レビューが集まるとより真のスコアに近づいていきます。  
# スコアが安定しているであろう、レビューが30件以上集まっている店に絞ります。

# In[127]:


tmp2 = tmp


# In[128]:


tmp = tmp2.query('review_num>30')


# In[129]:


cov_cols = ['genre', 'lunch', 'dinner', 'save_num', 'review_num', 'idv_over35', 'idv_under30']


# 出現数が少ないカテゴリ変数の処理

# In[135]:


over_10000 = ['￥10,000～￥14,999', '￥15,000～￥19,999', '￥20,000～￥29,999',
              '￥30,000～￥39,999', '￥40,000～￥49,999', '￥50,000～￥59,999']


# In[141]:


tmp.loc[tmp['lunch'].isin(over_10000), 'lunch'] = '¥10,000~'
tmp.loc[tmp['dinner'].isin(over_10000), 'dinner'] = '¥10,000~'


# In[142]:


minor_genres = list(pd.DataFrame(tmp['genre'].value_counts()).query('genre < 50').index)


# In[143]:


tmp.loc[tmp['genre'].isin(minor_genres), 'genre'] = 'その他'


# ## 因果効果の定義

# お店が有料会員になることによって、スコアがどう変化するのかを推定します。  
# あるお店が有料会員であるときのスコアを$y_1$, そうでない場合のスコアを$y_0$とすると、有料会員になる事によるスコアの変化量は以下のように表せます。  
# $
# E[y_1 - y_0] = E[y_1] - E[y_0]
# $
# 
# これを、因果効果と呼びます。  

# 実際には、あるお店$i$に対してそのお店が有料会員だった時のスコア$y_{1,i}$と無料会員だった時のスコア$y_{0,i}$両方を観測することは出来ません。  
# ある有料お店のスコア$y_{1,i}$を観測した後に、もしもボックスで因果を捻じ曲げてそのお店が無料会員だったことにして、タイムマシンで観測時点に戻れば$y_{0,i}$を観測出来ます。

# # 手法１: そのまま比較

# まずはざっくり有料店舗会員と無料店舗会員でスコアに差があるのかを調べてみます

# お店$i$が有料会員なら$z_i = 1$, 無料会員なら$z_i = 0$とすると、そのまま比較した式は以下のようになります。

# $$
# E[z y_1 - (1 - z) y_0] = E[y | z=1] - E[y | z=0]
# $$

# 有料会員になるか無料会員かが完全にランダムな場合、この推定量は因果効果$E[y_1 - y_0]$と一致します。

# In[145]:


fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(18, 12))
sns.barplot(x='is_paid', y='score', data=tmp, capsize=.2, ax=axes[0,0])
sns.barplot(x='score', y='genre', hue='is_paid', capsize=.2, data=tmp[tmp['genre'].isin(major_genres)], ax=axes[0,1])
sns.barplot(x='score', y='dinner', hue='is_paid', capsize=.2, data=tmp, ax=axes[1,0])
sns.barplot(x='score', y='lunch', hue='is_paid', capsize=.2, data=tmp, ax=axes[1,1])
axes[0,1].legend(loc='lower left')
axes[0,0].set_title('有料会員と無料会員のスコア比較')
axes[0,1].set_title('主要ジャンル別有料会員と無料会員のスコア比較')
axes[1,0].set_title('ランチ価格帯別有料会員と無料会員のスコア比較')
axes[1,1].set_title('ディナー価格帯有料会員と無料会員のスコア比較')
plt.show()


# 黒い線は信頼区間です。  
# 有料会員と無料会員でスコアの平均値に**差はほぼないことがわかります**  
# ジャンル別や昼/夜金額別にみても、有料会員と無料会員で目立った差は見られません。

# In[151]:


tmp.groupby(['is_paid'])[['score']].agg(['count', 'mean', 'median'])


# In[152]:


naive = tmp.query("is_paid==1")['score'].mean() - tmp.query("is_paid==0")['score'].mean()


# In[153]:


naive


# 全体の平均値の差は-0.007です。

# ## 母平均の差の検定

# 参考URL: https://ohke.hateblo.jp/entry/2018/05/19/230000

# 母平均の差-0.007が統計的に有意な差なのか偶然なのかを調べるために、母平均の差の検定を行います。  

# 店が有料会員なのか無料会員なのかが完全にランダムで、スコアの平均がそれぞれ独立した$\mu_1, \mu_0$の正規分布に従うと仮定します。  
# $$
# y_1 \sim \mathcal{N}(\mu_1, \sigma_1^2) \\
# y_2 \sim ,\mathcal{N}(\mu_0, \sigma_0^2) 
# $$

# (かなり無理のある仮定です。実際の分布の形はこのように歪で、特に無料会員は多峰です。)

# In[154]:


fig, axes = plt.subplots(ncols=2, figsize=(16,4))
sns.distplot(tmp.query('is_paid==1')['score'], kde=False, bins=np.arange(3, 5, 0.02), ax=axes[0])
sns.distplot(tmp.query('is_paid==0')['score'], kde=False, bins=np.arange(3, 5, 0.02), ax=axes[1])
axes[0].set_title('有料会員のスコア分布')
axes[1].set_title('無料会員のスコア分布')
axes[0].set_xticks(np.arange(3, 5, 0.1))
axes[1].set_xticks(np.arange(3, 5, 0.1))
axes[0].set_xlabel('スコア')
axes[0].set_ylabel('店舗数')
axes[1].set_xlabel('スコア')
axes[1].set_ylabel('店舗数')
plt.show()


# ### 1. 等分散が成り立つのか確認

# まず2群の母分散$\sigma_1^2, \sigma_0^2$が等しいとみなして良いのかを確認します。

# In[155]:


sigma2_1 = tmp.query('is_paid==1')['score'].var()
sigma2_0 = tmp.query('is_paid==0')['score'].var()
F = sigma2_1/sigma2_0
n = tmp.query('is_paid==1').shape[0]
m = tmp.query('is_paid==0').shape[0]
if stats.f.cdf(F, n-1, m-1) < 0.05:
    print('reject')
else:
    print('accept')


# 等分散性は仮定出来ませんでした。

# ### 2. 等分散性を仮定出来ない場合の検定

# In[156]:


stats.ttest_ind(tmp.query('is_paid==1')['score'], tmp.query('is_paid==0')['score'], equal_var=False)


# pvalue=0.33なので、有意水準10%で2群の差は無いと言えます。

# 実際は、店が有料会員なのか無料会員なのかはランダムとは限りません。  
# そこで、傾向スコアを用いてより正確な因果効果を推定します。

# # 手法2: IPW推定量を用いた因果効果の推定

# ### 概要

# 傾向スコアを用いてバイアスを取り除いた比較を行います。

# お店が有料会員になることによって、スコアがどう変化するのかを知りたいです。  
# あるお店が有料会員であるときのスコアを$y_1$, そうでない場合のスコアを$y_0$とすると、有料会員になる事によるスコアの変化量は以下のように表せます。  
# $
# E[y_1 - y_0] = E[y_1] - E[y_0]
# $
# 
# これを、因果効果と呼びます。  
# 
# 実際には、あるお店が有料会員だった場合のスコアとそうでなかった場合のスコアを両方観測することはできません。  
# そこで、有料会員になるかどうかが完全にランダムであれば、観測された有料会員と無料会員の平均値を比較することで因果効果を推定出来ます。  
# 
# しかし、実際はランダムではないです。例えば、居酒屋は有料会員の割合が高いのに対し、ラーメンは無料会員の割合が高いです。  
# 居酒屋の平均的なスコアが低く、ラーメンの平均的なスコアが高い場合、そのまま平均値の差分をとると因果効果はマイナスの値になってしまいます。  
# 
# そこで、有料お店会員/無料お店会員の平均値を計算するときに、あるお店が有料お店会員である確率を求めて、その確率の逆で重み付けします。  
# 先程の例だと、居酒屋は有料会員が多いので、有料会員の平均をとる際に小さい値で重み付けします。  
# 代わりに、ラーメンは有料会員が少ないので、有料会員の平均をとる際に大きい値で重み付けします。  
# 無料会員の場合も同様に、無料会員の居酒屋はレアなので大きく重み付けし、無料会員のラーメンはありふれているので小さく重み付けします。  
# これによって、有料会員と無料会員における居酒屋とラーメンの割合が均等になり、ランダムな状況を擬似的に作り出せます。

# <img src="rebalance.png" width=500>

# 理論の詳細についてはめっちゃくちゃわかりやすいスライドがあるので、参考にしてみてください。  
# 上の図はこのスライドの22ページから抜粋いたしました。  
# https://speakerdeck.com/tomoshige_n/causal-inference-and-data-analysis

# ### (興味がある人向け)数式

# お店$i$が有料会員なら$z_i=1$, 無料会員なら$z_i=0$, 有料会員である確率(傾向スコア)を$e_i$, スコアを$y_i$とします。  
# IPW推定量は以下のようになります。  

# $$
# \hat{E}[y_1] = \frac{\sum_{i=1}^N \frac{z_i y_i}{e_i}}{\sum_{i=1}^N \frac{z_i}{e_i}} \\
# \hat{E}[y_0] = \frac{\sum_{i=1}^N \frac{(1 - z_i) y_i}{1 - e_i}}{\sum_{i=1}^N \frac{(1 - z_i)}{1 - e_i}} \\
# $$

# ちなみに$E$にハットが付いているのは、真の値ではなく推定値だからです。  
# $E[y]$の時はyの個数に制約はなく、母集団と見なしていましたが、今回の$\hat{E}[y_1]$はN個のサンプル$y_1 \sim y_N$を用いて推定しています。

# ここで、真の傾向スコアの値$z_i$がわかっていて、「強く無視できる割り当て」条件(説明は省きます)が成立するなら、$\hat{E}[y_1]$は$E[y_1]$と一致します。

# ### どうやるのか

# 1. お店が有料会員なのか無料会員なのかを分類するモデルを作ります。  
# 2. お店が有料/無料会員である確率を求めて、その逆数をかけた値でスコアの平均値を重み付け平均します。  
# 3. 分類モデルが正確なら、重み付け平均の差分はバイアスが除去された純粋な因果効果となります。

# ### コード

# にのぴらさんの実装などが参考になります。  
# URL: https://pira-nino.hatenablog.com/entry/causal_inference_implement

# ### 1. 有料会員である確率を求めるモデルの作成

# ロジスティック回帰モデルを使います

# In[157]:


X = pd.get_dummies(tmp[cov_cols], drop_first=True).values


# In[159]:


y = tmp['is_paid'].values


# In[160]:


y_pred = np.zeros(y.shape[0])


# In[161]:


# StratifiedKFoldで予測値を計算する
sss = StratifiedKFold(n_splits=4, random_state=1234, shuffle=True)
for i, (train_index, test_index) in enumerate(sss.split(X, y)):
    model = LogisticRegression()
    model.fit(X[train_index, :], y[train_index])
    y_pred[test_index] = model.predict(X[test_index])


# モデルの正答率も一応確認します  

# In[162]:


y_pred = model.predict(X)
y_true = y


# In[163]:


confusion_matrix(y_true, y_pred)


# In[164]:


accuracy_score(y_true, y_pred)


# In[179]:


roc_auc_score(y_true, y_pred)


# AUCが0.7を超えたのでギリギリOKです。  

# ### 2. 傾向スコアの計算

# In[165]:


ps = model.predict_proba(X)


# In[166]:


tmp['ps_lr'] = ps[:,1]


# ### 3. 傾向スコアの分布を確認

# 傾向スコアの分布が有料会員と無料会員で左右対称になっているか、所属確率が満遍なく分布しているのかを確かめます。

# In[167]:


tmp.groupby(['is_paid'])['mise'].count()


# In[168]:


fig, ax = plt.subplots(figsize=(12, 4))
sns.distplot(tmp.query('is_paid==1')['ps_lr'], kde=False, norm_hist=True, bins=np.arange(0,1,0.05), ax=ax, label='有料会員')
sns.distplot(tmp.query('is_paid==0')['ps_lr'], kde=False, norm_hist=True, bins=np.arange(0,1,0.05), ax=ax, label='無料会員')
ax.legend()
plt.show()


# この分布が有料会員と無料会員で左右対称になるのがベストです。  
# どちらの分布でも、所属確率がきちんと[0,1]で分布していることを確認します。

# 有料会員になる確率が0.0~0.2となる場所の分布がイマイチです。  
# また、無料会員における分布が右に寄っています。  
# 教師データのバランスをアンダーサンプリング等で揃えてあげた方が良いのでしょうか...？  
# 
# ちょっとここら辺はわからないので、誰かご存知の方いらっしゃれば教えて頂けると嬉しいです。

# 次に、主要な説明変数に対して傾向スコアの逆数をかけたら平らに補正されるかを確かめます。  
# 雑にlightGBMのfeature importanceを確認します。

# In[191]:


# lightgbmのfeatureimportanceはバイナリより連続な変数を重要視する傾向がある
model = lightgbm.LGBMClassifier()
model.fit(X, y)
importance = pd.DataFrame(model.feature_importances_,
                          index=pd.get_dummies(tmp[cov_cols], drop_first=True).columns, 
                          columns=['importance'])
importance.sort_values(['importance'], ascending=False).head(3)


# In[200]:


fig, axes = plt.subplots(nrows=2, figsize=(16, 10))
sns.distplot(tmp.query('is_paid==1')['save_num'], kde=False, bins=np.arange(0, 20000, 200), label='有料会員', ax=axes[0])
sns.distplot(tmp.query('is_paid==0')['save_num'], kde=False, bins=np.arange(0, 20000, 200), label='無料会員', ax=axes[0])
sns.distplot(tmp.query('is_paid==1')['save_num']/tmp.query('is_paid==1')['ps_lr'], kde=False, bins=np.arange(0, 20000, 200), label='有料会員', ax=axes[1])
sns.distplot(tmp.query('is_paid==0')['save_num']/tmp.query('is_paid==0')['ps_lr'], kde=False, bins=np.arange(0, 20000, 200), label='無料会員', ax=axes[1])
axes[0].legend()
axes[1].legend()
axes[0].set_title('IPS補正前のsave_num分布')
axes[1].set_title('IPS補正後のsave_num分布')
plt.show()


# 傾向スコア(PS)の逆数(Inverse)をかけると分布の偏りが若干補正されて、よりランダム化された平らな分布になっているのがわかります。

# 統計量を取って確認してみます。

# In[221]:


tmp.groupby(['is_paid'])['save_num'].agg(['count', 'mean', 'median'])


# In[222]:


pd.DataFrame({
    'is_paid': tmp['is_paid'],
    'save_num_with_ips': tmp['save_num']/tmp['ps_lr']
}).groupby(['is_paid'])['save_num_with_ips'].agg(['count', 'mean', 'median'])


# In[212]:


# 極端に小さい傾向スコアの値がないかをチェック→まあ問題なさそう
tmp['ps_lr'].min()


# 傾向スコアの逆数をかけた方が、save_numの中央値や平均値がより同じ値に近づいているのがわかります。  
# 無料会員の平均値が若干爆発気味なのが気になります...実データ解析であれば、値が爆発するようなら傾向スコアの下限をクリッピングすることも検討した方が良さそうです。

# ### 4. 因果効果の計算

# In[169]:


e_y1 = np.sum(tmp['is_paid']*tmp['score']/tmp['ps_lr']) / np.sum(tmp['is_paid']/tmp['ps_lr'])


# In[170]:


e_y0 = np.sum((1 - tmp['is_paid']) * tmp['score'] / tmp['ps_lr']) / np.sum((1 - tmp['is_paid']) / tmp['ps_lr'])


# In[171]:


e_y1, e_y0


# In[172]:


ate = e_y1 - e_y0


# In[173]:


ate


# 直接有料会員と無料会員のスコアの平均値を比較した際は-0.007でしたが、IPW推定量は0.005になりました。  

# # 結論

# 有料会員になることは、お店のスコアをどの程度上げるのかという因果効果を調査しました。  
# レビュー数が30以上のお店がもし無料会員でなく有料会員だった場合、平均して約0.005スコアが上昇するという結果が出ました。  
# 
# 有料会員の平均と無料会員の平均の差を直接比較した-0.007と比べると、因果効果は正の値を取るようになりました。  
# 0.005が大きいのか小さいのかは議論の余地がありますが、  
# 
# * 有料会員になる事による因果効果は無いとは言い切れない  
# * しかし、お店全体で見ると、年会費を払わないとスコアを0.2下げる、という大きな因果効果は全体として見られなかった
# 
# といった歯切れの悪い解釈を述べて終わりにしたいと思います。

# ## 感想など

# IPW推定量は割り当てに対するバイアスを取り除いて純粋な効果を推定できるところは便利ですが、有意な差があるのかの検定を行うのが難しいなど、取り扱いがとても難しいです。   
# 傾向スコアマッチングやDRなど他の手法と比較して、様々な角度から検証する必要がありそうです。  
# 割り当て確率を予測するモデルがどの程度正確なら良いのかなどまだまだ勉強不足な点はありますが、**バイアスを全く除去しないで2群を比較するよりはより真の因果効果に近い値が推定出来る**はずなので、使いこなせるようになりたいです。  
# スコアの分布自体や共変量が割り当てに与える影響のメカニズムなど興味のない所は何も仮定を置かなくて良いのは、Rubinの因果モデルの良い所ですね。  
# 「評価が3.8以上の店に置いて年会費を払わないと3.6にスコアを下げる」という課題設定そのものに対する答えを提示できないのは歯がゆいです。もう少し上手いやり方をご存知でしたら教えていただけると嬉しいです。

# # (参考)スクレイピングコード

# In[ ]:


import os
import sys
import bs4
import requests
import time
import json
from hashlib import sha256
from pathlib import Path
import random
import re
import logging

# ログの保存
logging.basicConfig(filename="../log/taberogu_scraping.log",
                    format="%(asctime)s %(levelname)s: %(message)s",
                    level=logging.INFO)

def scan(u, page, chi):
    hs = sha256(bytes(u, 'utf8')).hexdigest()[:16]
    print('try', u)
    logging.info('try', u)
    if Path(f'comps/{hs}').exists():
        return
    headers = {'User-agent': 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_13_3) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/77.0.3865.90 Safari/537.36', 'referer': u}
    r = requests.get(u, headers=headers)
    soup = bs4.BeautifulSoup(r.text)

    lsts = soup.find_all('li', {'class': 'list-rst'})
    objs = []
    for lst in lsts:
        t = lst.find('a', {'class': 'cpy-rst-name'})
        s = lst.find('span', {'class': 'list-rst__rating-val'})
        rn = lst.find('em', {'class': 'list-rst__rvw-count-num cpy-review-count'})
        mise = t.text
        '''
        点が入っていないことがあるためハンドル
        '''
        if s is None:
            continue
        
        # 昼夜の価格帯
        l = lst.find('span', {'class': 'c-rating__val list-rst__budget-val cpy-lunch-budget-val'})
        d = lst.find('span', {'class': 'c-rating__val list-rst__budget-val cpy-dinner-budget-val'})
        
        # 店詳細に飛ぶ
        u2 = t['href']
        r2 = requests.get(u2, headers=headers)
        soup2 = bs4.BeautifulSoup(r2.text)
        # 有料会員か見分ける、有料会員なら店舗詳細ページに広告が載っていない
        p = soup2.find('div', {'id': '_popIn_recommend'})
        # 主要ジャンルを取得する
        g = soup2.find('meta', {'property': 'og:title'})
        # 全てのジャンルを取得する
        gs = soup2.find_all('script', {'type': 'application/ld+json'})
        # 保存済み数
        save_num = soup2.find('span', {'class': 'rdheader-rating__hozon-target'}).find('em', {'class': 'num'}).text
        
        time.sleep(random.uniform(1, 2))
        
        # 主要コメントに飛ぶ
        u3 = u2 + '/dtlrvwlst/'
        r3 = requests.get(u3, headers=headers)
        soup3 = bs4.BeautifulSoup(r3.text)
        idv_score = [r.text for r in soup3.find_all('b', {'class': 'c-rating__val--strong', 'rel': None})]
        comments = [t.text for t in soup3.find_all('a', {'class': 'rvw-item__title-target'})]
           
        is_paid = 1 if p is None else 0
        genre = re.findall('/.*\)', g['content'])[-1][1:-1]
        genres = re.findall(f'"servesCuisine".*', gs[0].text)[0].split('"')[3].split('、')
        lunch = l.text
        dinner = d.text
        score = s.text
        review_num = rn.text
        obj = {'mise':mise, 'is_paid': is_paid, 'genre': genre,
               'lunch': lunch, 'dinner': dinner, 'save_num': save_num,
               'score':score, 'page':page, 'chi':chi, 'genres': genres,
               'review_num':review_num, 'comments': comments, 'idv_score': idv_score}
        objs.append(obj)
        print(f'mise: {mise}, url: {u2}')
        logging.info(f'mise: {mise}, url: {u2}')
        time.sleep(random.uniform(1, 2))
    json.dump(objs, fp=open(f'../input/{hs}', 'w'), ensure_ascii=False, indent=2)
    print('complete', u)
    logging.info('complete', u)
    time.sleep(random.uniform(2, 4))

def main():
    chis = random.sample(list(range(1, 31)), 30)
    while True:
        chi = chis.pop()
        print(chi)
        logging.info(chi)
        pages = random.sample(list(range(1, 61)) , 60)
        for i in range(30):
            page = pages.pop()
            scan(f'https://tabelog.com/tokyo/A13{chi:02d}/rstLst/{page}/', page, chi)

if __name__ == '__main__':
    try:
        main()
    except KeyboardInterrupt as err:
        logging.error(err)
        print(err)
        sys.exit()