The One with ...

思いついたことや作業メモなどを公開しています

相対的剥奪のオンライン実験用コード

はじめに

とある原稿を書いているうちに,昔作ったztreeトリートメントをアップデートしたい気持ちがふつふつと高まってきたので,otreeで新たに実験用コードを作成しました.以下忘れた時用のメモです.

otree準備

  • アナコンダの最新版をインストール
  • pip install -U otreeでotreeをインストールする

こちらのサイトがたいへん参考になりました.

akrgt.gitbook.io

作業用ディレクトリの設定

異なるマシンで作業ができるようにdropbox上にフォルダを作ります.

たとえばC:\Users\Dropboxにドロップボックスのローカルアドレスがあれば

otree startproject C:\Users\Dropbox\oTree

という感じでホームディレクトリをつくってあげるとよいでしょう.

はじめてディレクトリをつくるとサンプルコードもdlするかと聞かれるので,これはdlしておきましょう.

今回作ったようなトリートメントはcournotゲームをベースにして作ると簡単でした.

サンプルを探せば自分がやりたい実験に似たコードがあるはずなので,まずはそれをベースに少しずつ修正を加えると作業が楽でしょう (フルスクラッチでいきなり書くのは難しそう)

以降はマシンを変える度にdb.sqlite3を消してから,テスト用のローカルサーバーを otree devserverでたてればOK.

すでに作成済みのトリートメントの実行

アナコンダのパワーシェルでトリートメントを含むディレクトリ移動する. たとえばdropbox上に保存した場合,研究室等のマシンから cd C:\Users\hogehoge\Dropbox\oTree などと入力し,ディレクトリを移動する.移動したディレクトリにはsettings.pyがあり,この中に目的のトリートメントを記入しておく. (新しく作った場合はあとでここで記入する)

移動したらテスト用のローカルサーバーを otree devserverでたてる.

ブラウザでhttp://localhost:8000/にアクセス

settingに登録したトリートメントが表示されるので実行(登録してないとメニューに出てこない) .各クライアントを別ウィンドウで開くためのリンクは,参加人数で設定した分だけ表示される.

アプリの準備

サンプルコードにcournotゲームのフォルダがあるのでコピーしてリネームしましょう. 相対剥奪なので,relative_deprivationという名前にします.

このフォルダには主に

  • __init__.py
  • Contribute.html
  • Results.html

の三つのファイルがあります.

この中身を適当に修正して実験コードを作成します.一つ上の階層

C:\Users\Dropbox\oTreeには

  • settings.py

というファイルがあります.この中の

SESSION_CONFIGS = [
    dict(
        name='guess_two_thirds',
        display_name="Guess 2/3 of the Average",
        app_sequence=['guess_two_thirds', 'payment_info'],
        num_demo_participants=3,
    ),
    dict(
        name='survey', app_sequence=['survey', 'payment_info'], num_demo_participants=1
    )]

に自分が作成したrelative_deprivation'を登録します.guess_two_thirdssurveyはデフォルトのdemoです.

結果は次のようなものになるでしょう.

SESSION_CONFIGS = [
    dict(
        name='guess_two_thirds',
        display_name="Guess 2/3 of the Average",
        app_sequence=['guess_two_thirds', 'payment_info'],
        num_demo_participants=3,
    ),
    dict(
        name='survey', app_sequence=['survey', 'payment_info'], num_demo_participants=1
    ),
    dict(
        name='relative_deprivation', app_sequence=['relative_deprivation', 'payment_info'], num_demo_participants=5
    )
]

settings.pyを保存したらotree devserverでテスト用サーバを立ち上げます.demoの中に自分が登録した実験(relative_deprivation)があれば成功です.

コードの中身

アプリを構成するコードの中身はこんな感じです

__init__.py

クラスの定義や計算用関数を定義するpythonコードです.

pythonになじみのある人は,なんとなく何をやっているか分かると思います.詳細は後ほど解説します.

from otree.api import *
import random


class C(BaseConstants):
    NAME_IN_URL = 'relative_deprivation'
    PLAYERS_PER_GROUP = 5
    NUM_ROUNDS = 1
    RETURN_RATE = random.randint(3,5)
    NUM_SUCCESS = random.randint(1,5)
    PROFIT = RETURN_RATE*100-100 #COST=100


class Subsession(BaseSubsession):
    pass


class Group(BaseGroup):
    x = models.IntegerField()

class Player(BasePlayer):
    action = models.IntegerField(
        choices=[[0, 'no'],[1,'yes (invest 100)']],
        label="Do you invest? ",initial=(0)
    )
    outcome = models.IntegerField(
        choices=[[5, 'very satisfied'],[4, 'satisfied'],[3, 'neither'],[2, 'dissatisfied'], [1, 'very dissatisfied']],
        label='How do you feel about the result?',
        widget=widgets.RadioSelect,
    )
    profit = models.IntegerField()


# FUNCTIONS
def set_profit(group: Group):#グループオブジェクトを引数にした関数
    players = group.get_players()#playerを取得
    group.x = sum([p.action for p in players])#投資者数の総数xを計算
    xids=[]#投資者のplayeridを格納するリスト
    for p in players:
        if p.action==1:
            xids.append(p.id_in_group)#投資者のplayeridを全てリストに追加
    #投資者数が多い場合の抽選
    if C.NUM_SUCCESS < group.x:
        win=random.sample(xids,C.NUM_SUCCESS)#NUM_SUCCESSだけxidsから勝者を選択
    for p in players:
        if p.action == 0:
            p.profit=0 # 非投資者の利得を計算
        else:
            if C.NUM_SUCCESS >= group.x:
                p.profit= 100*(C.RETURN_RATE*p.action-p.action)
            else:
                if p.id_in_group in win:
                    p.profit= 100*(C.RETURN_RATE*p.action-p.action)
                else:
                    p.profit= -100

# PAGES
class Contribute(Page):
    form_model = 'player'
    form_fields = ['action']


class ResultsWaitPage(WaitPage):
    after_all_players_arrive = set_profit


class Results(Page):
    form_model = 'player'
    form_fields = ['outcome']


page_sequence = [Contribute, ResultsWaitPage, Results]

Contribute.html

つぎは入力画面のhtmlです.

htmlタグ直書き世代の自分にとっては懐かしいコードです(わしらの若い頃は,これを手打ちして「ホームページ」をつくっておったのじゃよ).

ちなみに{{C.PROFIT}}のように変数名C.PROFITを囲むと,pythonで定義した変数をhtml上に表示させることができます.

{{ extends "global/Page.html" }}
{{ block title }}Contribute{{ endblock }}
{{ block content }}

    <p>
        You can choose "invest" or "not invest".
    </p>

    <p>
        {{C.NUM_SUCCESS}} players among those who invest will get a profit.
    </p>

 <p>
     You will get 
 <ul>
     <strong>{{C.PROFIT}}, when you win.</strong>
 </ul>
 <ul>
     <strong>- 100 , when you lose.</strong>
 </ul>
 <ul>
     <strong>0 , when you do not invest.</strong>
 </ul>


    {{ formfields }}

    {{ next_button }}

{{ endblock }}

Results.html

最後は結果表示と結果に対する満足度を聞くフォームを表示するhtmlです.

{{ block title }}Results{{ endblock }}
{{ block content }}

    <p>
       Your profit is {{ player.profit }}.
    </p>

    <p>
        Please answer the following questions.
    </p>

    {{ formfields }}


    {{ next_button }}

{{ endblock }}

実験の概要

相対的剥奪実験の具体例は以下の通りです.

  • 被験者は100円投資or投資しないを選択できます
  • 投資をえらんだ人は成功すれば300円もらえますが,失敗すると何ももらえません.
  • 成功人数は2名です.投資した人の数が2名を超えると,成功者の2名以外は100円失います

セッションによって利益率Rと成功者数nが変化します.一般的には

  • 被験者は100円投資or投資しないを選択できます
  • 投資をえらんだ人は成功すれば100*R円もらえますが,失敗すると何ももらえません.
  • 成功人数はn名です.投資した人の数がn名を超えると,成功者のn名以外は100円失います

というル-ルです

実験の流れ

直感的には実験は次のような順番で進行します.

  1. 実験用URLにアクセスしたプレイヤーにContribute.htmlが表示され,プレイヤーは行動(投資するか否か)を選択します
  2. 全プレイヤーが選択を終えると __init__.pyで定義した関数(set_profit)がプレイヤーの利得を計算します
  3. Result.htmlが表示され,各プレイヤーが受け取った利得がその中に表示されます.
  4. Result.htmlの中に満足度を入力するフォームが表示されます.プレイヤーが満足度を選択すると実験は終了です
  5. 結果が自動でデータベースに記録されます.

コードの概要

otreeでは基本変数と関数の定義を全て,__init__.py内に記述します.htmlはプレイヤーからの入力の受け取りや,計算結果の表示を担当します.

otreeでは実験用に基本的なクラスである

  • class C(BaseConstants): グローバル変数を定義するクラス
  • class Subsession(BaseSubsession): セッションを定義するクラス
  • class Group(BaseGroup): グループ属性を定義するクラス
  • class Player(BasePlayer): プレイヤー属性を定義するクラス

などがあらかじめ定義されています.

それぞれ簡単に内容を確認します.

class C(BaseConstants):
    NAME_IN_URL = 'relative_deprivation' #実験の名称です
    PLAYERS_PER_GROUP = 5 #被験者数です
    NUM_ROUNDS = 1 #ラウンド数です
    RETURN_RATE = random.randint(3,5) #利益率です
    NUM_SUCCESS = random.randint(1,5) #成功者数です
    PROFIT = RETURN_RATE*100-100 #COST=100 #成功時の利得です

以上がグローバル変数です.ラウンドごとに水準を変化させる場合は, 利益率などをsessionクラスで定義します.

class Group(BaseGroup):
    x = models.IntegerField()

グループの属性です.被験者集団全体で決まる値をここで属性として定義します. 属性xは投資者数を記録します.この値は被験者の行動によって変化するので,グローバル変数としては定義できないことに注意します.

class Player(BasePlayer):
    action = models.IntegerField(
        choices=[[0, 'no'],[1,'yes (invest 100)']],
        label="Do you invest? ",initial=(0)
    )
    outcome = models.IntegerField(
        choices=[[5, 'very satisfied'],[4, 'satisfied'],[3, 'neither'],[2, 'dissatisfied'], [1, 'very dissatisfied']],
        label='How do you feel about the result?',
        widget=widgets.RadioSelect,
    )
    profit = models.IntegerField()

プレイヤーのクラスです.ここでは

  • action #行動
  • outcome #結果への評価(満足度)
  • profit #受け取る利得

の3属性を定義します.idのような基本属性はもともと定義されているので,呼び出すだけで使えます.

def set_profit(group: Group):#グループオブジェクトを引数にした関数
    players = group.get_players()#playerを取得
    group.x = sum([p.action for p in players])#投資者数の総数xを計算
    xids=[]#投資者のplayeridを格納するリスト
    for p in players:
        if p.action==1:
            xids.append(p.id_in_group)#投資者のplayeridを全てリストに追加
    #投資者数が多い場合の抽選
    if C.NUM_SUCCESS < group.x:
        win=random.sample(xids,C.NUM_SUCCESS)#NUM_SUCCESSだけxidsから勝者を選択
    for p in players:
        if p.action == 0:
            p.profit=0 # 非投資者の利得を計算
        else:
            if C.NUM_SUCCESS > group.x:
                p.profit= 100*(C.RETURN_RATE*p.action-p.action)
            else:
                if p.id_in_group in win:
                    p.profit= 100*(C.RETURN_RATE*p.action-p.action)
                else:
                    p.profit= -100

利得計算関数です.本実験の肝の部分です.投資者数xnで設定した成功者数を下回る場合は,投資者全員が成功し, n<xの場合はランダムに成功者が決まるようにrandom.sample(xids,C.NUM_SUCCESS)で計算しています. そのため,投資選択者のidリストをxids.append(p.id_in_group)で作ります.このときotreeのプレイヤーのデフォルト属性であるid_in_group を利用します.

このようにotreeにはクラス毎によく利用する属性が準備されています.

# PAGES
class Contribute(Page):
    form_model = 'player'
    form_fields = ['action']


class ResultsWaitPage(WaitPage):
    after_all_players_arrive = set_profit


class Results(Page):
    form_model = 'player'
    form_fields = ['outcome']


page_sequence = [Contribute, ResultsWaitPage, Results]

htmlのクラスです.フォームからの入力などを定義します. 入力の受け取りにはプレイヤークラスで定義した属性をform_fieldsで指定します. 回答選択肢などは,プレイヤークラスの方で定義しています.




まとめ

otreeで相対的剥奪オンライン実験用コードを作成しました.

ポイントは

  • クールノーゲームのようにベースとなる実験を見つけ,それを上書き修正する
  • クラス毎に定義すべき属性を理解する

です.ztreeに比べるとかなりコードが書きやすくなった印象です. pythonの基本とhtmlのタグを知っていれば,簡単な実験なら2~3時間くらいで書けると思います.

Win10からUbuntu経由でPythonを使う

Win10からUbuntu経由でPythonを使おうとしたら思いのほか手間取ったので作業メモを残しておく.

問題:WindowsPowerShellからWSLをインストールできない

原因:AMDのCPUが仮想環境を使える状態になっていない

解決法:BIOS→アドバンス設定→OC→CPUの機能→SVM mode→状態を「無効」から「有効」に切り替え.AMDのコアのデフォルト設定ではSVMが無効になっているようだ. インテルの場合も似たようなCPUの設定を探すと仮想環境対応にできる. CPUの状態が整ったかどうかは,タスクマネージャのCPUのプロパティで確認できる.仮想環境が有効になっていればWindowsPowerShellからWSLをインストールできる. コアの設定後にマイクロソフトストアから好きなディストリビューション(たとえばUbuntu 20.04LTS)を選択してインストールしてもよい. WSLからインストールする場合はディストリビューションまで指定すればうまくいく.

参考:

intint000.blog96.fc2.com

問題:WSLからgitにSSH接続できない

原因:SSH接続用の公開鍵をGithubで設定していない

解決法:鍵を作成して,Githubアカウントで設定する.まずWSLにて

cd ~/.ssh
ssh-keygen

で鍵を作成.clipコマンドだとうまくいかないので

cat ~/.ssh/id_rsa.pub

で中身を表示してコピーする.その後Win上のブラウザからgithubにログインし,マイプロフィールにて公開鍵を設定. WSLに戻ってgitでクローンすると,アクセスが許可される.

なお別のローカルクライアントからgitにSSH接続する場合は,鍵を追加で作成後,githubにログインしてssh keyを追加する.githubにログインせずにシェルから直接公開できるようだが,うまくいかなかった.

参考:

tonariho.com

問題:WSLからjupyter notebookが起動しない

解決法:エクスプローラでwslのファイルにアクセスし,jupyter_notebook_config.pyを直接書き換える

\\wsl$

でWSL上のファイルにWindowsから直にアクセスして

jupyter_notebook_config.py

を探す.メモ帳などで開けたらファイル冒頭に

c.NotebookApp.use_redirect_file = False

を書き込む.

その後WSLからjupyter notebookを入力.するとブラウザから自動でnotebookが開く. 編集したファイルはWSL上に保存される.

なおディストリビューションが異なる別のローカルクライアントは,この方法ではうまくいかなかった. jupyter notebook listローカルアドレストークンを表示して,直接ブラウザに打ち込んでアクセスした. configファイルの指定が違うのかもしれない.

参考:

cuxe.livedoor.blog

問題:Win10上のpythonでGinzaを使おうとしたらSudachiPyのインストールで失敗する.

解決策:未解決.Ubuntu(WSL)でCコンパイラを先にインストールするとうまくいった.

感想

自分だけが使うぶんには便利だけど,ちょっと授業では使えないかな.

OSごとSSDを別のマザーボードに移植する作業の手順

OSがインストールされたM.2SSDドライブを,異なるマザーボードにCPUと一緒に移植するとOSはそのまま起動できるのでしょうか? またチップセットとCPUが異なる場合はどうでしょうか?

実験した結果を記録を残しておきます.

実験1 (CPUは同じだがチップセットのグレードだけ異なるマザボへの移植)

元の構成(メーカーPC)
移植先
手順

あとでWindowsのPINコードを再設定する必要があるために,WindowsアカウントのログインIDとパスワードを準備しておきます.

またマイクロソフトからの確証コードを受け取るために携帯電話等も必要です.

パスワードを控えたら,もとのマザーボードからCPU,SSD,メモリを外して,B365マザーボードに載せ替えます *1. 構成としてはマザボ,CPU,SSD,メモリのみで映像出力はCPU内蔵グラフィックスを使います. キーボードとマウスはUSB2.0接続の無線タイプをつなげました.

実験1の結果

電源スイッチを押すとWindowsは起動しました.しかしログインしようとすると,PINを再設定しろというメッセージが出て,Windowsのアカウント入力を要求されます. ここであらかじめ控えておいた情報が役立ちます(PCが一台しかない場合は,ここで困ると推測). Windowsアカウントを入力するとSMSで確証コードが携帯電話に送られてきます.コードを入力してPINを再設定すると,Windowsのデスクトップが立ち上がります. 結局必要なのはWindowsアカウントの確証のみでした.確証コードの受け取りに携帯電話かタブレットが必要です. 確証後の動作は特に問題ありません.

実験2 (CPUもチップセットも異なるマザボへの移植)

元の構成
移植先
実験2の結果
  • 世代が異なるチップセット(初代から第9世代)のマザボSSDを載せ替えてもOSは起動した(ちなみにP55は10年以上前に組んだマシンだが現役)
  • win10のためか,PINの再設定やWindowsアカウントの入力も必要なし.ただし起動後の動作は明らかに遅く,なにか不具合が生じている可能性がある.
  • 数日後にログインすると,Windowsアカウントの確証を求められた.マザボMACアドレス等で識別しているのかもしれない.
分かったこと
  • メーカー製のマザボが独自のバックプレートを装着している場合は,トルクスドライバ(T20)で外す.ただしCPUソケットを固定するには専用のステーかナットが必要.このステーはチップセットが異なっても共通なのでジャンクマザボから使い回すことが可能.
  • 高さ164mmのクーラが入るならAS500(PLUS)が最善の選択.他のパーツとの干渉なし,冷える,静か,低価格で最高.OCしないのであればアサシンとかNH-D15は必要なさそう.
  • CPU用ATX12V8ピンコネクタに,4ピンだけ挿しても起動する.今回使用した電源ユニットで使える補助電源はATX12V4ピンのみ.TDP65WくらいのCPUなら問題なさそう.
  • PCケースのフロントパネルから電源スイッチだけ外してマザボ(Front_Panel)につなげば,バラックでOSを起動できる(この方法だとCPUクーラの取り付けとOS起動確認が簡単).
  • M.2 Wifiコネクタにカードを挿せば,PCIスロットを節約できる.アンテナ線はカードにスナップで取り付け可能(ハンダ付けじゃない).アンテナ線の先はケース内部に貼りつけでOK.っていうかM.2にWifi接続できるとか,おっさんは知らんし.
極私的結論
  • ストレージはM.2-MVRe SSD2TBをマザボに2枚挿せば十分.というか起動ドライブもSATA-SSDで十分.
  • もはやPCケースに5インチベイとHDDベイは不要.つまり165mmの空冷クーラー用ハイトさえ確保できればよい.
  • マザボの拡張スロットは2個で十分.wifiもM.2接続で十分.つまりマイクロATXで十分.
  • クーラはAS500かAK400で十分.簡易水冷は不要. 
  • 静音化最後の壁は電源ユニット.ハイブリッドで無音化可能.
  • BIOSでCPUの電圧下げて静音化しても,パフォーマンスはそれほど低下しない.
おまけ CPUの設定

RyzenTDPを下げる方法の一つは,BIOSでEco modeを指定する. MSIの場合,delキーでBIOS設定.Setting→AMD overclocking→Eco mode.Ryzen9-3900xの場合は65Wに下がる.

クロックスピードの調整

  • BIOSをアドバンスモードに切り替える
  • CPU Ratio 43→33くらいに倍率を下げる
  • CPU Core voltage 1.4→1.0(Autoのままでもよい)

なお,DragonCenter(MSIのPCパーツ管理ソフト)でサイレントモードを選ぶとクロックスピードが定格の半分になるようです.

*1:マザボをケースに入れると作業しにくいので,ケースの外に置いたマザボにフロントパネルコネクタから電源スイッチのケーブルだけを挿しておくと便利です.高級なマザボはテスト用にスイッチが基盤に付属してますが,なくても電源スイッチだけ外付けできます.

ggplot2で複数の折れ線をまとめてplot

すぐ忘れるのでメモ.

library("ggplot2")

y1 <- rnorm(10,5,1)
y2 <- rnorm(10,2,1)
y3 <- rnorm(10,3,1)
years <- 2001:2010

data <-data.frame(years,y1,y2,y3)
ggplot(data, aes(years)) +
  geom_line(aes(y = y1, colour = "y1")) +
  geom_line(aes(y = y2, colour = "y2"))+
  geom_line(aes(y = y3, colour = "y3"))

f:id:hamada7418:20211216173840p:plain

確率変数の組の独立性について

確率変数の組の独立性

とあるエコノメのテキストに確率変数の組 $$(X_i,Y_i), i=1,2,...,n$$ が独立ならば $$X_j, j=1,2,...,n$$も独立である(大意),と書かれていた.

まあようするに確率変数のペア$(X_i,Y_i)$と$(X_j,Y_j)$が独立ならば,ペアの片方だけを取り出した $X_i$と$X_j$も独立である,ということなのだが,どうも釈然としない.

そもそも $$(X_i,Y_i), i=1,2,...,n$$ が独立であることの定義も書かれていない.

定義して確かめる

テキストに書いてないのだから,自分で定義してもよかろう,ということで1変数の$n$次元独立のアナロジーで,確率変数の順序対の独立を以下のように定義する.

定義(確率変数の組の独立性)

$$(X_i,Y_i), i=1,2,...,n$$ が独立であるとは順序対$(x_1,y_1), (x_2,y_2), ..., (x_n, y_n)$の同時分布について $$ f( (x_1,y_1), (x_2,y_2), ..., (x_n, y_n) )=f(x_1,y_1)f(x_2,y_2)\cdots f(x_n, y_n) $$ が成り立つことである.

なぜこのような定義が必要かと言えば,次のような理由である.上記の定義において$x$と$y$との間には依存関係があってもよい. よって$x$と$y$との間に推測したい関数関係を念頭に置いているが,データとして個人$i$の情報と個人$j$の情報とは無関連であるとき上記の定義が利用できる.

実際,説明変数が確率変数であるような回帰モデルを考える場合には,この仮定が必要である.

さて知りたいことは「 $(X_i,Y_i), i=1,2,...,n$が独立である$\Longrightarrow X_j, j=1,2,...,n$が独立である」 が成立するかどうかである.

証明

仮定より $$ f( (x_1,y_1), (x_2,y_2), ..., (x_n, y_n) )=f(x_1,y_1)f(x_2,y_2)\cdots f(x_n, y_n) $$ が成立している.

左辺を$y_1, y_2, ..., y_n$で積分する. $$ \iint \cdots \int f( (x_1,y_1), (x_2,y_2), ..., (x_n, y_n) )dy_1 d y_2 \cdots d y_n=f( x_1,x_2, ..., x_n) $$ 次に右辺も$y_1, y_2, ..., y_n$で積分する. $$ \iint \cdots \int f(x_1,y_1)f(x_2,y_2)\cdots f(x_n, y_n)dy_1 d y_2 \cdots d y_n=f(x_1)f(x_2)\cdots f(x_n) $$

この二つは等しいから $$ f( x_1,x_2, ..., x_n)=f(x_1)f(x_2)\cdots f(x_n) $$ である.

このことはつまり, $$ X_j, j=1,2,...,n $$ が独立であることを意味する.

ふう.

これで安心して $$(X_i,Y_i), i=1,2,...,n$$ が独立であるとき $$ E[u_i|X_1,X_2,...,X_n]=E[u_i|X_i] $$ っていう仮定が使える.

傾向スコアによるバイアス補正についてのメモ

はじめに

介入の割り当て確率が共変量に依存して異なるとき,統制群と介入群の応答変数の平均の差によって介入効果を推定するとバイアスが生じます. もし,共変量の値毎に割り当て確率を決定する関数が分かっていると,推定量のバイアスを取り除くことができます.

ただし,確率を決定する関数が分かっている点がポイントで,テキトーに関連がありそうな共変量でロジスティック回帰をすればいいわけではありません.

以下にテキトーなロジスティック回帰だと失敗する例を示します.

バイアス割り当て実験

分析に使うパッケージの読み込み

library("MatchIt")
library("Matching")
library("cobalt")
library("WeightIt")
library("broom")

共変量xによって介入(トリートメント)の割り当て確率が異なる状況を考えます.

x-5から5までの範囲に値をとる一様分布とします. 潜在的結果y0は,この共変量と相関しています.

真の介入効果はtau=3とします.よって潜在的結果y1y1=y0+3で定義できます.

n <- 1000;# sample size   
tau <- 3;# 真の介入効果    
x <- runif(n,-5,5)# 共変量
y0 <-5+2*x+ rnorm(n,mean=0,sd=1)#潜在的結果yo, xの影響を受ける
y1 <- y0+tau #潜在的結果y1

介入の割り当て確率をロジスティック関数で定義します.

p <- 1/(1+exp(-x))# make probabilities by logistic function
d <- rbinom(n,size=1,p)# treatment with bias

pは割り当て確率です.共変量xの値が大きいほど割り当て確率が高くなっています.

dは介入を表すダミー変数です. dを確率pのベルヌーイ分布で定義するため,xの値が大きい個体ほど介入を受けやすくなります.

以上のコードで割り当てにバイアスがある状況を表現できました.

観測データyを以下のように定義します*1

y <- c()
  for(i in 1:n){# 観測値yの生成
    if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}
data <- data.frame(y,d,p,x)

介入があった個体については潜在的結果y1を,介入がなかった個体についてはy0を応答変数として観測します. 各個体は,y0もしくはy1のどちらかしか観測できないことに注意します(因果推論の根本問題)

共変量xと割り当て確率pの関係を図で確認しておきましょう.

以下,割り当て確率のことを「傾向スコア」と呼びます*2

f:id:hamada7418:20211202225451p:plain
共変量xと傾向スコアpの関係

共変量xと介入dの関係も図で確認しておきましょう.

f:id:hamada7418:20211202225737p:plain
共変量xと介入dの関係

共変量xの値が小さいほど介入の割り当てがなく,値が高いほど介入が割り当てられています.

このようなランダムでない割り当ては,介入効果の推定にバイアスをもたらします.

計算して確かめてみましょう.

> summary(lm(data=data,y~d))

Call:
lm(formula = y ~ d, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.3599  -2.7284   0.0375   2.8791  11.4583 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8242     0.1699   4.852 1.42e-06 ***
d            11.8107     0.2384  49.549  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.768 on 998 degrees of freedom
Multiple R-squared:  0.711,  Adjusted R-squared:  0.7107 
F-statistic:  2455 on 1 and 998 DF,  p-value: < 2.2e-16

介入効果の推定値は11.8107でした.真の効果は3ですから,過大推定になっています.

回帰ではなく,条件付き期待値の差をとっても同じ結果を得ます.

mean(y[d==1])-mean(y[d==0])

傾向スコアマッチング

傾向スコアはすでに計算済みですから,それをつかってマッチングデータを作り,ATEを推定します.

> matched.data <- Match(Y = data$y, Tr=data$d, X =p,estimand="ATE")
> summary(matched.data)#Matched number of observations...............  500 

Estimate...  3.4254 
AI SE......  0.23556 
T-stat.....  14.541 
p.val......  < 2.22e-16 

Original number of observations..............  1000 
Original number of treated obs...............  508 
Matched number of observations...............  1000 
Matched number of observations  (unweighted).  1152 

直感的に言えば,関数Matchは「傾向スコアが同じだけど介入の有無が異なる個体をペアにしたデータ」を作ります.

結果,バイアスが減り,真の効果に近いATE推定値を得ることができました. 念のためcobaltパッケージで共変量の調整結果を可視化します

data2 <- matchit(d~x,data=data,method = "nearest",replace = T)
love.plot(data2, thresholds = 0.1)#SMD

f:id:hamada7418:20211202231936p:plain
調整前後での共変量平均値の比較

傾向スコアによるマッチングで共変量のバランスがとれたことが分かります.

真の傾向スコアが手に入らない場合

上記の実験では,真の傾向スコアが分かっている条件下でマッチングデータを作りました.しかし現実の観察データを使って分析する場合,傾向スコアが事前に分かっていることなどありません.そもそも傾向スコアを決定する共変量の全てがデータとして入手できる状況など,極めて希でしょう.

そこで次に,必要な共変量が手に入らない場合にGLMで傾向スコアを推定してマッチングに利用した場合,どのくらい真の効果量の推定にバイアスが生じるのかを確かめます. まず分析用のデータの生成です.

n <- 3000;# sample size   
tau <- 3;# 真の介入効果    
x1 <- runif(n,-5,5)# 共変量
x2 <- runif(n,-5,5)# 共変量
y0 <-x1+x2#潜在的結果y0, x1,x2の影響を受ける
y1 <- y0+tau #潜在的結果y1
p <- 1/(1+exp(-x1-x2))# make probabilities by logistic function
# pは真の傾向スコア.共変量x1とx2の値に依存する.
d <- rbinom(n,size=1,p)# treatment with bias
y <- c()
for(i in 1:n){# 観測値yの生成
  if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}

真の傾向スコアは共変量x1x2によって定まると仮定しています. 次に,必要な共変量が入手できない場合とできた場合を比較します.

out.glm <- glm(d~x1,family = "binomial")#傾向スコアを過小定式化モデルで推定
out.glm2 <- glm(d~x1+x2,family = "binomial")#フルモデルで推定

p.predict <- out.glm$fitted.values
p.predict2 <- out.glm2$fitted.values

p.predictx2が足りない場合の傾向スコアの推定値で,p.predict2x1x2が正しくGLMに含まれる場合の推定値です. 結果を比較してみます.

まず共変量が足りない場合

matched.data <- Match(Y = y, Tr=d, X =p.predict,estimand="ATE")
summary(matched.data)

Estimate...  7.2681 
AI SE......  0.10906 
T-stat.....  66.645 
p.val......  < 2.22e-16 

Original number of observations..............  3000 
Original number of treated obs...............  1509 
Matched number of observations...............  3000 
Matched number of observations  (unweighted).  7892 

真の効果3に対して推定値は7.2681ですから過大推定しています(GLMのMcFadden R2 は 0.2209495でした).

このとき共変量のバランスをチェックすると,共変量x1に関しては調整後にバランスはとれています.

dat <- data.frame(d,x1,x2,y)
match.out <- matchit(data=dat,d~x1,method = "nearest",replace = T)
love.plot(match.out, thresholds = 0.1)

f:id:hamada7418:20211215125227p:plain
調整後の共変量

つまり「投入した共変量のバランスがとれたとしても,効果量の推定には失敗する可能性がある」ので注意が必要です.

次に必要な共変量が入手できる場合のマッチング推定の結果です.

 matched.data2 <- Match(Y = y, Tr=d, X =p.predict2,estimand="ATE")
 summary(matched.data2)

Estimate...  3.4448 
AI SE......  0.1594 
T-stat.....  21.611 
p.val......  < 2.22e-16 

Original number of observations..............  3000 
Original number of treated obs...............  1509 
Matched number of observations...............  3000 
Matched number of observations  (unweighted).  4921 

バイアスが補正されました(GLMのMcFadden R2 は 0.5866154でした).

実際の観測データを用いた分析では,傾向スコアの正しい定義に必要な共変量は不明です. データに含まれる共変量のバランスがとれたとしても,処置への割り当て確率を計算するモデルに欠落変数があると,処置効果の推定に失敗します. つまり処置に影響しそうな共変量を分析者が(勝手に)考えて傾向スコアを予測しても,うまくいく保証はありません*3

補足

必要な共変量が揃っている場合の,誤差項の分散の影響も確かめてみました. 複数回マッチングさせるために関数化します.

f1 <- function(s){
  n <- 1000;# sample size    
  tau <- 3;# 真の介入効果  
  x <- runif(n,-5,5)# 共変量
  y0 <-x+ rnorm(n,mean=0,sd=3)#潜在的結果yo, xの影響を受ける
  y1 <- y0+tau #潜在的結果y1
  p <- 1/(1+exp(-x-rnorm(n,mean=0,sd=s)))# make probabilities by logistic function
  # pは真の傾向スコア.共変量xの値が大きいほど割り当て確率が高い
  d <- rbinom(n,size=1,p)# treatment with bias
  y <- c()
  for(i in 1:n){# 観測値yの生成
    if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}
  out.glm <- glm(d~x,family = binomial)
  Lm <- out.glm$deviance/(-2)
  L0 <- out.glm$null.deviance/(-2)
  p.predict <- pout.glm$fitted.values
  #予測した傾向スコアを使いMatchでマッチさせる
  matched.data <- Match(Y = y, Tr=d, X =p.predict,estimand="ATE")
 c(matched.data$est[1], 1-Lm[1]/L0[1])# McFadden R^2
}

7行目p <- 1/(1+exp(-x-rnorm(n,mean=0,sd=s)))の部分で,共変量xから傾向スコアを計算する過程で誤差を導入します.誤差項の標準偏差sが大きくなると,共変量xによって予測できる傾向スコアの精度が低下します.

f:id:hamada7418:20211212001754p:plain
共変量xと傾向スコアpの関係.誤差項の標準偏差s=3

トリートメントの割り当て確率は,この傾向スコアによって決まるので,割り当てdを共変量xに回帰すると傾向スコアの予測値を得ます.

最後にMatch関数に予測した傾向スコアを与えてマッチングデータを作成します.計算の手続きをまとめ関数f1の引数は,誤差の標準偏差sのみです.アウトプットは「マッチングデータから推定したATE」と「傾向スコアを予測したロジスティック回帰モデルのMcFadden 疑似決定係数R2」です.

計算例を確認してみましょう.

 Mac <- c()
 Est <- c()
 for(i in 1:50){
   Est <- append(Est,f1(s=3)[1])
   Mac <- append(Mac,f1(s=3)[2])
 }
 
> mean(Mac)
[1] 0.2395028
> mean(Est)
[1] 2.998437

McFadden疑似決定係数が0.24程度でも,真の効果をうまく推定できています.

 Mac <- c()
 Est <- c()
 for(i in 1:50){
   Est <- append(Est,f1(s=10)[1])
   Mac <- append(Mac,f1(s=10)[2])
 }
 
> mean(Mac)
[1] 0.03893244
> mean(Est)
[1] 3.002857

McFadden疑似決定係数が0.039程度でも,真の効果をうまく推定できています.この計算例から,疑似決定係数の大きさはバイアス補正にとって本質的ではないことが分かります.たとえ疑似決定係数が高くても,必要な共変量が足りなければ推定には失敗します.

逆確率による重み付け

傾向スコアを使ったIPWによる調整も試しておきます.ウェイトの計算にはWeightItパッケージを使いました. weightit関数がdataを引数にとるのでdata.frameを定義しておきます.

データセットの定義

n <- 3000;# sample size   
tau <- 3;# 真の介入効果    
x1 <- runif(n,-5,5)# 共変量
x2 <- runif(n,-5,5)# 共変量
y0 <-x1+x2#潜在的結果y0, x1,x2の影響を受ける
y1 <- y0+tau #潜在的結果y1
p <- 1/(1+exp(-x1-x2))# make probabilities by logistic function
d <- rbinom(n,size=1,p)# treatment with bias
y <- c()
for(i in 1:n){# 観測値yの生成
  if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}
dat <- data.frame(d,x1,x2,y)

過小定式化モデル(共変量x2が足りないモデル)と適切なモデルによって傾向スコアを計算し,それぞれ重みに変換して回帰に使用します.

## 重み付きデータでの効果の推定 傾向スコアを過小定式化モデルで推定
weighting <- weightit(data=dat, d ~ x1,method = "ps",estimand = "ATE")
IPW_result <- lm(y ~ d,weights = weighting$weights) %>% tidy()
IPW_result


## 重み付きデータでの効果の推定 フルモデルで推定
weighting2 <- weightit(data=dat, d ~ x1+x2,method = "ps",estimand = "ATE")
IPW_result2 <- lm(y ~ d,weights = weighting2$weights) %>% tidy()
IPW_result2

結果の比較です.

> IPW_result
# A tibble: 2 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    -2.16    0.0679     -31.9 4.56e-192
2 d               7.45    0.0959      77.7 0.       


> IPW_result2
# A tibble: 2 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)   -0.269    0.0923     -2.92 3.58e-  3
2 d              4.02     0.134      30.1  6.63e-174

共変量が足りないモデルではうまくバイアスが補正できませんでした. 共変量がそろったモデルではバイアスが(相対的にはうまく)補正できました.ただし傾向スコアが0に近い値がデータセットに多く含まれるせいで,推定結果にはまだバイアスが残っています.

 参考文献

  • 安井 翔太,2020『効果検証入門』技術評論社
  • ローゼンバウム,2017=2021『統計的因果推論入門』共立出版
  • 星野崇宏,2009『調査観察データの統計科学 因果推論・選択バイアス・データ融合』岩波書店

*1:データフレームはなくても構いませんが,ここでは他のコードで利用するために定義しています

*2:傾向スコアとは「共変量Xの値がxであるすべての個体iに関する処置群への割り当て確率の平均値」です

*3:もちろん,偶然成功する場合もあります

数理モデル解説本の内容紹介

はじめに

2018年12月に『その問題,数理モデルが解決します』をベレ出版さんから刊行しました.

当blogではこの本の内容紹介をしていなかったので,以下にざっと解説します*1

本書は数学,統計学,経済学(行動経済学ゲーム理論),社会学(数理社会学)から,筆者が特におもしろいと思う《人や社会に関する数理モデルアルゴリズム》を紹介しています.2人のキャラクターが対話しながら「ある問題」を解決するための「数理モデル」や「アルゴリズム」を考えます.基本的には1章で1話が完結します.

続編『その問題,やっぱり数理モデルが解決します』と世界観は共通ですが,内容は独立しているので,どちらから読み始めていただいても問題ありません.概ね高校・大学初年度レベルの数学を使っています.

特徴

単にこんなモデルがあるよ,と紹介するだけでなく,その数学的な原理をフォローできるように証明や計算を省略せずに書きました.この本で紹介するモデルやアルゴリズムは,「一度はどこかで聞いたことがあるもの」であり,目新しさはないかもしれません.本書は「こういう命題は知ってるけど,どうして成立するんだろう」とか「こういうアルゴリズムは知ってるけど,どうして有効なんだろう」という疑問に答えます.

たとえば,5章の「秘書問題」の解決アルゴリズムが「全体の36.8%を観察してから,次にベストを超えた対象を選ぶ」であることは,よく知られています.しかしどうして36.8%(より正確に書くと1/e)を見送る必要があるのか,その《しくみ》を説明することは簡単ではありません.本書は,強力なアルゴリムが成立する《しくみ》をなるべく丁寧に解説します.

対象読者

本書の読者として以下のような方を想定しています.

  • 数学って何の役に立つの?と考えている中高生
  • 数理モデル*2に興味のある大学生
  • 統計学を使うけど,仕組みがよく分からないという研究者
  • データ分析を勉強しようかなと考えている社会人

各章のタイトル・難易度・トピック

章タイトル  難易度 トピック     
序  モデルとはなにか 
第1 章 隠された事実を知る方法 ☆☆ 回答のランダム化、集合、確率変数、期待値、分散
第2 章 卒業までに彼氏ができる確率 ☆☆ ベルヌーイ分布、確率変数の合成、2 項分布
第3 章 内定をもらう方法 ☆☆☆ 2 項分布の期待値、確率変数の和の期待値、インプリケーション、ベータ分布、ベータ2 項分布
第4 章 先延ばしをしない方法 ☆☆ 先延ばしのメカニズム、時間割引、準双曲型割引、先延ばしの防止
第5 章 理想の部屋を探す方法 ☆☆☆ 秘書問題、グーゴル・ゲーム、最適停止問題、アルゴリズム、全体の36.8% を見送る理由、生活満足度と通勤時間
第6 章 アルバイトの配属方法 選好、安定マッチング、DA アルゴリズム、パレート効率性
第7 章 売り上げをのばす方法 ☆☆☆ ランダム化比較試験、条件付き期待値、潜在的結果、不偏推定量、統計的検定
第8 章 その差は偶然でないと言えるのか? ☆☆☆☆☆ 検定の仕組み、棄却域、対立仮説、正規分布の性質、サンプルサイズの設計
第9 章 ネットレビューは信頼できるのか? ☆☆☆☆ 陪審定理、チェビシェフの不等式、大数の弱法則陪審定理の一般化
第10 章 なぜ0 円が好きなのか? ☆☆☆ ゼロ価格アノマリー、効用関数と導関数プロスペクト理論、価値関数、ゼロ価格効果の一般化
第11 章 取引相手の真意を知る方法 ☆☆☆ ゲーム理論と支配戦略、第2 価格封印入札、メカニズムデザイン、ナッシュ均衡
第12 章 お金持ちになる方法 ☆☆☆ ギャンブルでお金持ちになる方法、所得分布のカタチ、累積効果、対数正規分布の生成

以下,各章の概要です.

各章の内容

第1 章 隠された事実を知る方法

規則やモラルに反する行動の調査では,回答者が正直に答えないことが知られています.そこで「禁煙の場所と知りながら煙草を吸ったことがある」などの逸脱行動を把握するために《回答のランダム化》という手法を使って解決を試みます(この手法のもとで母比率を推定するためには,実は《大数の法則》という命題が必要です.その解説は9章で詳しくかきました)・

第2 章 卒業までに彼氏ができる確率

100人と出会ったとき,1人以上から好かれる確率はどのくらいだろうか?この確率は2項分布をつかったモデルで計算できます. 2項分布は統計学でよく使う基礎的な部分の1つであり,より単純なベルヌーイ分布の合成から《つくる》ことができます.その過程を詳しく解説しました.

第3 章 内定をもらう方法

第2章の内容を展開します.モデルからインプリケーションを引き出したり、ベータ分布を使ったモデルの拡張について考えます.

第4 章 先延ばしをしない方法

先延ばし行動に関する心理学や行動経済学の知見を紹介しつつ,時間割引モデルについて解説します.先延ばし防止に役立つ知識もまとめています.

第5 章 理想の部屋を探す方法

いわゆる「秘書問題」モデルを解説しています.ファーガソンのレビュー論文内の証明を,できるだけ丁寧に再構成しました. 他の章と比べて少しだけ難易度が高くなっていますが,1/eの謎が解けるはずです.コンピュータで計算するための簡単なコードも記載しています. こちらを参照

第6 章 アルバイトの配属方法

ゲールとシャプレーが提唱したマッチングアルゴリズムを紹介します.アルゴリズム自体は簡単なのですが,なぜこれが安定マッチングを生むのか,その理屈(証明)を詳しく書きました.

第7 章 売り上げをのばす方法

ランダム化比較試験をルービンの因果推定モデルをベースにして紹介しました.平均処置効果の推定では条件付き期待値の計算がとても重要なので,その基礎部分を詳しく書きました.

第8 章 その差は偶然でないと言えるのか?

検定の仕組みを解説しました.サンプルサイズを計算する式を見たことのある人は多いと思いますが,あの式がどこから出てくるのか疑問に思うことでしょう. 私自身もどういう理屈なのか気になったので,いろいろ調べて理解した内容をここで解説しています.

第9 章 ネットレビューは信頼できるのか?

陪審定理について解説しています.人文社会系では有名な定理ですが、証明は書かれていないことが多いように思いました. 第1章で使った大数の法則を利用して、陪審定理の証明を説明します.また少しだけ定理の一般化についても考えています.

第10 章 なぜ0 円が好きなのか?

《100円の商品の10円の値引き》と《10円の商品の無料化》はどちらも値引き額は10円ですが,後者の方が人を引きつけることが知られています. どうして同じ値引き額なのに,われわれは無料を選ぶのでしょうか?行動経済学の実験を紹介しつつ,効用関数やプロスペクト理論の価値関数について学びます.

第11 章 取引相手の真意を知る方法

第2 価格封印入札というオークションのルールがあります.これは「入札額が一番高い人が勝者だが,支払い額は2番目に高い入札額となる」というルールです.このルールを使うと,理論上は,入札者は評価額を正直に表明するインセンティブを持つことが知られています. なぜそうなるのか?という理屈をゲーム理論のモデルを使って説明します.

第12 章 お金持ちになる方法

所得の分布はパレート分布や対数正規分布という確率分布で近似できることが知られています.しかし世の中の人々は,所得分布の形など意識することなく,生活していることでしょう.なぜ所得分布が特定の確率分布で近似できるのか,そのしくみを簡単な数理モデルで説明します.このモデルの原型は私が大学院生のときに作りました.

メッセージ

私は大学院に入るまでは,「自分は文系だから数学は人生に無用だ」と考えていました.

しかし師の影響で数理モデルの勉強をはじめたことがきっかけとなり,モデルの楽しさに目覚めました. これ以上の娯楽は,なかなかない,と今では思います.

本書をきっかけに,数理モデル(数学)っておもしろいかも,と思っていただけたら望外の喜びです.

*1:続編(『その問題,やっぱり数理モデルが解決します』2020年9月刊行)との違いはなんですか?という質問を読者から時々いただきます.この紹介を参考にしてください.

*2:数理モデルは,数学を使って現象を表現・説明・予測する方法の総称であり,観察した現象から本質を見抜く洞察やデータを分析するための指針(理論)を与えてくれます.現象が成立するメカニズムを数理モデルで特定できれば,現象の理解や未来の予測に役立ちます.