雑記

お気付きの点がありましたらご指摘いただけますと幸いです。

  1. 宮川 雅巳. 統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学). 朝倉書店. 2004.

💜

f:id:cookie-box:20211229165400p:plain:w70

  • 何度も挙げている例ですが、「ピアノを習っているか」X から「東大に合格したか」Y への効果がどれくらいあるかを知りたいとしましょう。このとき、「ご家庭の年収」U という交絡因子があるために正しい効果の大きさを測れないことが予想されます。ご家庭の年収で層別するなどの対策は考えられるでしょう。しかし、ご家庭の年収が観測できるとは限りません。
  • そこで [1] の 89 ページでは図 5.5 (b) でいう Z なる変数を観測できていたとしたらどうか、という議論をしています。上の例であれば「ピアノのレッスンに払っていた月謝」などでしょうか。現実にはお金持ちほど高額な月謝のピアノのレッスンに通うのかもしれませんが、ここでは都合上「お金持ちほどピアノを習う傾向があるが、ピアノを習うならばご家庭の年収と月謝は独立である」と仮定しましょう。実際そこまで世の中のピアノのレッスンに月謝の格差があるわけではないのではないでしょうか。全然わかりませんが。ご家庭の年収は観測できずピアノの月謝が観測できる状況というのが考えにくいのも問題ですが致し方ありません。
  • ともかくその独立性を仮定すると、
     「ご家庭の年収とピアノの月謝の相関係数
     =「ご家庭の年収とピアノを習っているかの相関係数
      ・「ピアノを習っているかとピアノの月謝の相関係数
    です。適当な数値例をつくってみましょう。ピアノを習っている下でご家庭の年収とピアノの月謝が独立になるように注意します……この例だとピアノを習っていることの東大合格への効果 \alpha_{yx} はゼロになってしまいますね。

    import statistics
    
    def standardization(l):
        l_mean = statistics.mean(l)
        l_pstdev = statistics.pstdev(l)
        return [(i - l_mean) / l_pstdev for i in l]
    
    x = standardization([0, 0, 0, 0, 0, 0, 1, 1, 1, 1])  # ピアノを習っているか
    y = standardization([0, 0, 0, 1, 0, 0, 1, 0, 0, 1])  # 東大に合格したか
    u = standardization([1, 2, 3, 4, 5, 5, 4, 4, 5, 5])  # ご家庭の年収
    z = standardization([0, 0, 0, 0, 0, 0, 1, 2, 1, 2])  # ピアノの月謝
    
    rho_zu = statistics.correlation(z, u)
    rho_xu = statistics.correlation(x, u)
    rho_xz = statistics.correlation(x, z)
    print(f'rho_zu = {rho_zu}')
    print(f'rho_xu = {rho_xu}')
    print(f'rho_xz = {rho_xz}')
    print(f'rho_xu * rho_xz = {rho_xu * rho_xz}\n')
    # rho_zu = rho_xu rho_xz
    
    # Y = alpha_yz Z + alpha_yu U + epsilon_y ... (5.21.3)
    # => X Y = alpha_yz X Z + alpha_yu X U + epsilon_y X
    # => rho_xy = alpha_yz rho_xz + alpha_yu rho_xu
    # => alpha_yu = (rho_xy - alpha_yz rho_xz) / rho_xu ... (a)
    
    # Y = alpha_yz Z + alpha_yu U + epsilon_y ... (5.21.3)
    # => Y Z = alpha_yz Z Z + alpha_yu Z U + epsilon_y Z
    # => rho_yz = alpha_yz + alpha_yu rho_zu
    # => rho_yz = alpha_yz + alpha_yu rho_xu rho_xz ... (b)
    
    # (a) (b) から alpha_yz を解く.
    # rho_yz = alpha_yz + (rho_xy - alpha_yz rho_xz) rho_xz
    # => alpha_yz = (rho_yz - rho_xy rho_xz) / (1.0 - rho_xz rho_xz)
    
    rho_xy = statistics.correlation(x, y)
    rho_yz = statistics.correlation(y, z)
    print(f'rho_xy = {rho_xy}')
    print(f'rho_yz = {rho_yz}')
    print(f'rho_xy * rho_xz = {rho_xy * rho_xz}')
    
    alpha_yz = (rho_yz - rho_xy * rho_xz) / (1.0 - rho_xz * rho_xz)
    print(f'alpha_yz = {alpha_yz}')
    print(f'alpha_yx = rho_xz * alpha_yz = {rho_xz * alpha_yz}\n')
    
    # Y = alpha_yx X + alpha_yu U + epsilon_y ... (5.20.2)
    # => X Y = alpha_yx X X + alpha_yu X U + epsilon_y X
    # => rho_xy = alpha_yx + alpha_yu rho_xu
    # => alpha_yx = rho_xy - alpha_yu rho_xu
    alpha_yu = (rho_xy - alpha_yz * rho_xz) / rho_xu
    print(f'alpha_yu = {alpha_yu}')
    print(f'alpha_yu * rho_xu = {alpha_yu * rho_xu}')
    print(f'alpha_yx = rho_xy - alpha_yu * rho_xu = {rho_xy - alpha_yu * rho_xu}')
    
    rho_zu = 0.3957336397583147
    rho_xu = 0.4308202184276646
    rho_xz = 0.9185586535436918
    rho_xu * rho_xz = 0.3957336397583148
    
    rho_xy = 0.3563483225498992
    rho_yz = 0.3273268353539886
    rho_xy * rho_xz = 0.3273268353539886
    alpha_yz = 0.0
    alpha_yx = rho_xz * alpha_yz = 0.0
    
    alpha_yu = 0.8271392736637095
    alpha_yu * rho_xu = 0.3563483225498992
    alpha_yx = rho_xy - alpha_yu * rho_xu = 0.0