2017年6月15日(金)-16日(土)理論合宿@日本橋

統計モデルと一般化推定方程式
推定量の一致性と漸近正規性の証明
・Wedderburn, R. W., Quasi-likelihood functions, generalized linear models and theGauss—Newton method. Biometrika, Vol.61, No.3, pp.439-447, 1974.
・Liang, K. Y., Zeger, S. L., Longitudinal data analysis using generalized linear
models. Biometrika, Vol.73, No.1, pp.13-22, 1986.

2017年9月16日(金) 14:00〜14:40
M1 井澤 佳那子

本日の発表内容
1. 最近の私の研究について：選択行動における量子観測の情報系とその理論化

2. 統計モデルと推定
3. 一般化線形モデルと最尤推定
4. 擬似尤度スコア方程式（Wedderburn, 1974）
5. 一般化推定方程式（Zeger, S. L., Liang, K. Y. , 1986）

Introduction2

モチベーション
都市内の交通行動や都市の形成に対して自然な表現を与えることで評価・分析を行い，それに基づいた都市
交通計画立案に繋げたい．

特に交通行動について，観測量の代数的構造に自然な表現を与えることが，現象（都市内の交通行動原理）
へのより良い理解に繋がる．
cf. 最大エントロピー原理

私の研究について3

背景：観測と行動理論
1. 観測
・データの多様化（quantization, sampling）：情報通信・観測技術の進展やセンサネットワークの普及

・情報量の増大：点単位での行動把握が可能となり，交通行動に関する情報量が増大

・点データの特徴量：点データの特徴量に合わせた空間の集計単位とスケールに焦点を当てた整理（+それら正規化する行動判別手法）が必要．

・情報量的解釈の必要性：点ベースの交通行動データについて，情報量を定義し評価することで，観測量によるモデル性能を直感的に理解することができ，他の観測データとの比較や異なるモビリティ間の比較が容易になる．
私の研究について4背景：観測と行動理論
2. 行動理論・任意の選択行動に対する評価尺度：既出の行動モデルは，実質的に特定のモビリティを対象としているものや実際に真のパラメータを推定することが困難なもの．

・情報量の異質性と相関：パラメータ推定では観測データから得られる情報量が十分に大きいことが前提．交通行動の相関性により当該情報量が減少する可能性．

・推定上の特性，計算性向上：交通行動判別においても行動モデルのパラメータ推定においても，情報量を定義することによる計算性の向上に対する寄与．

交通行動分析における量子観測の情報系を再整理．
私の研究について5

フレームの検討
選択行動の情報量を選択結果チェーンが選択行動に関して持つ情報量と定義し，選択モデルのパラメータの全体分布とチェーンから推定されるパラメータの部分分布の相対性を利用して，
当該情報量を実験的に評価する．
e.g. 交通分析の情報系アプローチ(Sopasakis, A. et al., 2016）𝑐ℎ𝑎𝑖𝑛1𝑐ℎ𝑎𝑖𝑛2𝑐ℎ𝑎𝑖𝑛𝑘𝑑+,-./0𝑐ℎ𝑎𝑖𝑛𝐾

𝑑+,-./+,

観測遷移数 𝐿 = ∑65 𝑥5
任意の空間単位に紐づけられる
選択結果の１遷移

私の研究について

1. 序論
2. 行動モデルの定式化
3. 情報量的解釈
4. 情報量規準（モデル選択）
5. 推定，パラメータの特性
6. 実データ推定
7. 結論
6情報理論（Shannon, 1948）[基礎]
情報：不確実性を減少させてくれるもの．
性質① 発生確率の低いことを知れた方が情報量が多い．
性質② 情報量は足し合わされる．

・エントロピー𝐻 𝑦 ≔ − ; 𝑞 𝑦 log 𝑞 𝑦 𝑑𝑦

・Kullback-Leibler ダイバージェンス𝑞(𝑦)𝐷[𝑞 𝑦 : 𝑝(𝑦)] ≔ ; 𝑞 𝑦 log𝑑𝑦𝑝(𝑦)= ; 𝑞 𝑦 log 𝑞(𝑦) 𝑑𝑦 − ; 𝑞 𝑦 log 𝑝(𝑦) 𝑑𝑦
𝑝 𝑦 , 𝑞(𝑦)：確率変数𝑦の確率分布
私の研究について7

選択行動における情報量的解釈
今，𝒒 𝒚 が真の確率分布，𝒑 𝒚; 𝜷 が推定した確率分布とすると
𝐷[𝑞 𝑦 : 𝑝(𝑦; 𝛽)] = ; 𝑞 𝑦 log 𝑞(𝑦) 𝑑𝑦 − ; 𝑞 𝑦 log 𝑝(𝑦; 𝛽) 𝑑𝑦
一定; 𝒒 𝒚 𝒍𝒐𝒈 𝒑(𝒚; 𝜷) 𝒅𝒚 が大きいと𝑫[𝒒 𝒚 : 𝒑(𝒚; 𝜷)]は小さくなる．

𝟏; 𝒍𝒐𝒈 𝒑(𝒚; 𝜷) 𝒅𝒚 → ; 𝒒(𝒚) 𝐥𝐨𝐠 𝒑(𝒚; 𝜷) 𝒅𝒚
𝑲(データ数 𝑲 → ∞); 𝐥𝐨𝐠 𝒑(𝒚; 𝜷) 𝒅𝒚 をパラメタ𝜷について最大化すれば良い[最尤推定法]．

・選択行動(e.g. MNL)ではこれが選択確率になっている．
・配分(e.g. Dial, 1971)は，ローディングの問題なので単にネットワーク上のエントロピー最小化問題．
私の研究について8情報量的解釈の意義
・大数の法則の仮定の妥当性．
𝟏; 𝒍𝒐𝒈 𝒑(𝒚; 𝜷) 𝒅𝒚 → ; 𝒒(𝒚) 𝐥𝐨𝐠 𝒑(𝒚; 𝜷) 𝒅𝒚
𝑲(データ数 𝑲 → ∞)・本来最小化すべきは，KLダイバージェンス．
𝐷[𝑞 𝑦 : 𝑝(𝑦; 𝛽)] = ; 𝑞 𝑦 log 𝑞(𝑦) 𝑑𝑦 − ; 𝑞 𝑦 log 𝑝(𝑦; 𝛽) 𝑑𝑦
観測量の多様化に伴う真の分布に対する仮定

モデル発展に伴う母数の高次元化

・任意の選択行動に対する評価尺度の一般化．

私の研究について9

・Wedderburn, R. W., Quasi-likelihood functions, generalized linear models and theGauss—Newton method. Biometrika, Vol.61, No.3, pp.439-447, 1974.
・Liang, K. Y., Zeger, S. L., Longitudinal data analysis using generalized linear
models. Biometrika, Vol.73, No.1, pp.13-22, 1986.

10

統計モデルと推定
推定とは確率変数𝒚がある分布族𝒑(𝒚, 𝜷)に従うとして，観測したデータ𝒚をもとに分布を求めること．
線形回帰モデル（最小線形二乗法）

一般化線形モデル（最尤推定法）
母数の次元高度化

GEE，GMM（Newton法など）

ベイズモデル（MCMCなど）

久保拓弥. (2012). データ解析のための統計モデリング入門 (pp. 106-108). 岩波書店.

統計モデルと推定

11

推定と幾何学的理解
確率変数の一部は観測されないことが多い…．
確率変数は，𝒚 = (𝒙, 𝒛)のよう分割され，𝒙は観測できて𝒛は観測できない隠れ変数であるとする．（𝑺：確率分布空間）
隠れ変数のない場合

𝑺𝒙e.g. GEE

隠れ変数のある場合

𝑺𝒚

𝑞b(𝑥)
経験分布

𝛽a

𝑞b 𝑥 𝑟(𝑧|𝑥)𝑞b d/0 𝑥 𝑝(𝑧|𝑥, 𝛽)𝑫[𝒒 𝒙 : 𝒑(𝒙)]MDe.g. emM𝑝(𝑥|𝛽)

ad
d/0 𝛽
a
𝛽
𝑝 𝑦 𝛽 = 𝑝(𝑥, 𝑧|𝛽)

統計モデルと推定

12

推定と幾何学的理解【おまけ】確率変数の一部は観測されないことが多い…．
確率変数は，𝒚 = (𝒙, 𝒛)のよう分割され，𝒙は観測できて𝒛は観測できない隠れ変数であるとする．（𝑺：確率分布空間）
構造推定

𝑺𝒚
D

M

理論モデルのパラメータを推定する

𝑺𝒚f𝑞b 𝑥 𝑟(𝑧|𝑥)

𝛽ad

𝛽ad/0

D

M

a∗

𝛽

一致する

𝑝 𝑦 𝛽 = 𝑝(𝑥, 𝑧|𝛽)
統計モデルと推定

𝛽a∗𝑞b 𝑥 f 𝑟(𝑧′|𝑥′)

ad
d/0 𝛽
a
𝛽
𝑝 𝑦′ 𝛽 = 𝑝(𝑥′, 𝑧′|𝛽)13

一般化推定方程式の背景
数学的には，一般化線形モデル(GLM)の拡張である擬似尤度ス
コア推定方程式の特殊な場合．
GLMの仮定仮定①：従属変数が指数分布族に従う．
仮定②：従属変数がそれぞれ独立である．尤度関数 L(𝜃|𝑦) = ∏65l0 𝑓(𝑦5 )e.g. 継時的，繰り返し測定，クラスタデータ
pdf

確率変数が指数分布に従う
AND
それぞれ独立統計モデルと推定

0 事象Aの確率変数

14

表記の整理・𝒚𝒊𝒕 ：従属変数𝑖 = 1, … , 𝐾：対象の数
𝑡 = 1, … , 𝑛5 ：観測回数
・𝒙𝒊𝒕 ：𝑝×1の共分散ベクトル・𝝁𝒊 = 𝑬(𝒚𝒊𝒕 |𝒙𝒊𝒕 )：期待値
・𝜷： 𝑝×1の回帰係数ベクトル
・𝒀𝒊 = 𝒚𝒊𝟏 , … , 𝒚𝒊𝒏𝒊
・𝑿𝒊 = 𝒙𝒊𝟏 , … , 𝒙𝒊𝒏𝒊𝑻𝑻：𝑛5 ×1の従属変数ベクトル：𝑛5 ×𝑝の𝑖(𝑖 = 1, … , 𝐾)番目の共分散ベクトル・𝜽𝒊𝒕 = 𝒉(𝜼𝒊𝒕 )：パラメータ・𝜼𝒊𝒕 = 𝒙𝒊𝒕 𝜷
統計モデルと推定

15

GLMと最尤推定
Step1. 確率変数𝒀𝒊 の確率分布を仮定．
指数分布族は一般に，

𝑓 𝑦5 = exp{{𝑦5 𝜃5 − 𝑎 𝜃5 + 𝑏(𝑦5 )}𝜙}

Step2. 確率変数𝒀𝒊 の期待値𝝁𝒊 の回帰モデルを設定．
𝜇5 = 𝐸 𝑦5d 𝑥5d = 𝑎f 𝜃5d = 𝑎f ℎ 𝜂5d= 𝑎′(ℎ(𝑥5d 𝛽))

Step3. 仮定した確率分布から尤度関数𝐋(𝜽|𝒚)を作り，これを} )を求める．
最大化する𝜷(𝜷𝜕→対数尤度関数を𝛽について偏微分して極値を求める． e.g. 𝜕𝜷 logL 𝜃 𝑦 = 0
𝑲

スコア関数

𝑼 𝑰 𝜷 = Œ 𝑿 𝒊 𝑻 ∆𝒊 𝑺 𝒊 = 𝟎
𝒊l𝟏

where

∆5 = 𝑑𝑖𝑎𝑔(𝑑𝜃5d ⁄𝑑𝜂5d )，𝑆5 = 𝑌5 − 𝑎′(𝜃5d )

一般化線形モデルと最尤推定

16} の性質最尤推定量𝜷性質①. 一致性
確率1で，𝛽a6 → 𝛽 ∗ (𝐾 → ∞)性質②. 漸近正規性
𝐾 𝛽a6 − 𝛽 ∗ → 𝑁 0, 𝐼(𝛽 ∗ )•0(𝐾 → ∞)

𝛽a6 ：推定値，𝛽∗ ：真値，𝐾：対象数
一般化線形モデルと最尤推定

17擬似尤度スコア推定方程式の着想(Wedderburn, 1974)

従属変数の真の分布が，仮定した分布でない場合にも妥当な統計的推論を行いたい．
e.g. ポアソン回帰モデルのover-dispersion問題

GLMの仮定の再確認
仮定①：従属変数が指数分布族に従う．←この仮定①を緩める
仮定②：従属変数がそれぞれ独立である．尤度関数 L(𝜃|𝑦) = ∏65l0 𝑓(𝑦5 )𝑌5 の従う確率分布が指数分布族であるとき平均と分散は，

𝑬 𝒚𝒊 = 𝝁𝒊 = 𝒂f 𝜽𝒊 ， 𝑽𝒂𝒓 𝒀𝒊 = 𝒂ff 𝜽𝒊 /∅で，任意の関数 𝑣(œ) によって関係は一意に決まる．

𝑽𝒂𝒓 𝒀𝒊 = 𝒗(𝝁𝒊 )/∅擬似尤度スコア方程式

18

擬似尤度
分散𝑽𝒂𝒓 𝒀𝒊 を扱うのはやめてある重み𝑽𝒊 として考えてみる．
𝑽𝒊 = 𝒗(𝝁𝒊 )/∅重み𝑽𝒊 を用いた擬似尤度スコア推定方程式を解いて得られる推} が，確率分布（分散）を仮定して計算した場合の最尤推
定量𝜷定量𝜷∗ と同じ性質（一致性・漸近正規性）を持てば妥当な統計的推論であるということ．
memo
𝑌5 の従う確率分布が指数分布族であるとき平均と分散は，

𝐸 𝑦5 = 𝜇5 = 𝑎f 𝜃5 ， 𝑉𝑎𝑟 𝑌5 = 𝑎ff 𝜃5 /∅で，任意の関数 𝑣(œ) によって関係は一意に決まる．

𝑉𝑎𝑟 𝑌5 = 𝑣(𝜇5 )/∅擬似尤度スコア方程式

19擬似尤度スコア推定方程式
𝑲

𝑼𝑸𝑳 𝜷 = Œ 𝑫𝒊 𝑻 𝑽𝒊 •𝟏 𝑺𝒊 = 𝟎
𝒊l𝟏

where

𝜕{𝑎f 5 (𝜃5 )}𝐷5 == 𝐴5 ∆5 𝑋5𝜕𝛽𝑉5 = 𝑣(𝜇5 )/∅𝑆5 = 𝑌5 − 𝑎′(𝜃5 )
𝐴5 = 𝑑𝑖𝑎𝑔{𝑣(𝜇5 )}memo
スコア推定方程式𝑲

𝑼𝑰 𝜷 = Œ 𝑿𝒊 𝑻 ∆𝒊 𝑺𝒊 = 𝟎

where∆5 = 𝑑𝑖𝑎𝑔(𝑑𝜃5d ⁄𝑑𝜂5d )，𝑆5 = 𝑌5 − 𝑎′(𝜃5d )𝒊l𝟏

擬似尤度スコア方程式20

漸近正規性の証明} を真値𝜷∗ のまわりでテーラー展開する．
𝑼𝑸𝑳 𝜷
𝑈¤¥ 𝛽a − 𝑈¤¥ 𝛽∗ =
1
𝐾

∴

𝑈¤¥ 𝛽

} − 𝜷∗𝑲 𝜷∗

𝜕𝑈¤¥ 𝛽 |¦l¦∗
𝜕𝛽

1 𝜕
= −
𝑈 𝛽 |¦l¦∗𝐾 𝜕𝛽 ¤¥

𝛽a − 𝛽∗ + 𝑂¨ (1)𝐾 𝛽a − 𝛽∗ + 𝑂¨ (1)

•𝟏
𝟏 𝝏
𝟏
∗
= −
𝑼𝑸𝑳 𝜷 |𝜷l𝜷
𝑼𝑸𝑳 𝜷∗ + 𝑶𝒑 (𝟏)
𝑲 𝝏𝜷
𝑲

=

𝑯•𝟏

where
𝐻= −

∗

𝐁

+ 𝑶𝒑 (𝟏)1 𝜕
1𝑈¤¥ 𝛽 |¦l¦∗ ，𝐵 =𝑈¤¥ 𝛽∗𝐾 𝜕𝛽
𝐾

擬似尤度スコア方程式21

漸近正規性の証明（続き）∴} − 𝜷∗
𝑲 𝜷

•𝟏
𝟏 𝝏
𝟏
∗
= −
𝑼𝑸𝑳 𝜷 |𝜷l𝜷
𝑼𝑸𝑳 𝜷∗ + 𝑶𝒑 (𝟏)
𝑲 𝝏𝜷
𝑲

𝑯•𝟏

=

∗

𝐁

+ 𝑶𝒑 (𝟏)

6

1 𝜕
𝜕 1
𝜕𝜇5 ² •0𝐻=−
𝑈¤¥ 𝛽 |¦l¦∗ = − limŒ
𝑉5 (𝑌5 − 𝜇5 ) = 𝑂¨ (1)6→± 𝜕𝛽 𝐾
𝐾 𝜕𝛽
𝜕𝛽
5l0

𝐵=

1
𝐾

𝑈¤¥ 𝛽 ∗ = lim

6→±（∵ 台数の法則）

6𝜕𝜇5 ² •0Œ
𝑉5 (𝑌5 − 𝜇5 ) ~𝑁𝑜𝑟𝑚𝑎𝑙 𝑑𝑖𝑠𝑡𝑟𝑖𝑏𝑢𝑡𝑖𝑜𝑛
𝜕𝛽
𝐾

1

5l0

（∵ 中心極限定理）∴} − 𝜷∗𝑲 𝜷~𝑵𝒐𝒓𝒎𝒂𝒍 𝒅𝒊𝒔𝒕𝒓𝒊𝒃𝒖𝒕𝒊𝒐𝒏

擬似尤度スコア方程式22一致性の証明：平均
平均
𝐸 𝐵 =𝐸𝐾 𝛽a − 𝛽 ∗

𝐸
1

6

Œ

𝑁 5l0

𝜕𝜇5
𝜕𝛽

∴

²

= 𝐸 HB + 𝑂¨ 1

𝑉5 •0 𝑌5 − 𝜇5𝑬=

} − 𝜷∗𝑲 𝜷= 𝑂¨ 1 ∗ 𝐸 𝐵 + 𝑂¨ (1)
1

6

Œ

𝐾 5l0

𝜕𝜇5 ² •0
𝑉5 𝐸 𝑌5 − 𝜇5
𝜕𝛽

=0

=𝟎

擬似尤度スコア方程式23一致性の証明：分散分散

𝑉𝑎𝑟

𝐾 𝛽a − 𝛽 ∗= 𝐸 { 𝐾 𝛽a − 𝛽 ∗ }{ 𝐾 𝛽a − 𝛽 ∗ }²
= 𝐸 𝑯•𝟏 𝑩𝑩𝑻 (𝑯•𝟏 )𝑻 +𝑂¨ (1)= 𝑯•𝟏 𝐸 𝑩𝑩𝑻 (𝑯•𝟏 )𝑻6

𝜕 1
𝜕𝜇5 ² •0𝑯=−
𝑉5 𝑌5 − 𝜇5
Œ
𝜕𝛽 𝐾
𝜕𝛽
5l0

𝜕
𝜕𝜇5 ² •0→−𝐸[
𝑉5 (𝑌5 − 𝜇5 )]𝜕𝛽𝜕𝛽𝜕 Â 𝜇5= −𝐸[
𝜕𝛽² 𝜕𝛽

²

²
𝜕𝜇
𝜕 •0
5
𝑉5 •0 𝑌5 − 𝜇5 −𝑉5 •0
𝑉5𝑉5 •0 𝑌5 − 𝜇5
𝜕𝛽𝜕𝛽𝜕𝜇5 ² •0 𝜕+
𝑉5 {𝑌 − 𝜇5 }]
𝜕𝛽
𝜕𝛽 5
𝜕𝜇5

𝜕𝜇5 ² •0
=
𝑉5
𝜕𝛽
𝜕𝛽

∴

6

1
𝜕𝜇5 ² •0 𝜕𝜇5
𝑯 = lim Œ
𝑉5
6→± 𝐾
𝜕𝛽
𝜕𝛽
5l0

擬似尤度スコア方程式24一致性の証明：分散（続き）
分散

𝐾 𝛽a − 𝛽 ∗𝑉𝑎𝑟= 𝐸 { 𝐾 𝛽a − 𝛽 ∗ }{ 𝐾 𝛽a − 𝛽 ∗ }²
= 𝐸 𝑯•𝟏 𝑩𝑩𝑻 (𝑯•𝟏 )𝑻 +𝑂¨ (1)= 𝑯•𝟏 𝐸 𝑩𝑩𝑻 (𝑯•𝟏 )𝑻= 𝑰𝟎 •𝟏 𝑰𝟏 𝑰𝟎 •𝟏} − 𝜷∗ ~ 𝑵(𝟎, 𝑰𝟎 •𝟏 𝑰𝟏 𝑰𝟎 •𝟏 )
𝑲 𝜷

∴

where

6

𝑌5 が指数分布族に従うときは

1
𝑰𝟎 = lim Œ6→± 𝐾

𝜕𝜇5 ² •0 𝜕𝜇5
𝑉5
𝜕𝛽
𝜕𝛽

𝑉5 = 𝑉𝑎𝑟(𝑌5 )なので𝐼Ç = 𝐼0
∴ 𝐾 𝛽Æ − 𝛽∗ ~ 𝑁(0, 𝐻•0 )1
𝑰𝟏 = lim Œ6→± 𝐾

𝜕𝜇5 ² •0𝜕𝜇5𝑉5 œ 𝑉𝑎𝑟(𝑌5 ) œ 𝑉5 •0
𝜕𝛽
𝜕𝛽

5l0
6

5l0

6

𝐸 𝐁𝐁

𝐓

1
= lim Œ6→± 𝐾

𝜕𝜇5 ² •0•0 𝜕𝜇5²
𝑉5 𝐸[ 𝑌5 − 𝜇5 𝑌5 − 𝜇5 ]𝑉5𝜕𝛽𝜕𝛽1= lim Œ6→± 𝐾

𝜕𝜇5 ² •0𝜕𝜇5𝑉5 œ 𝑉𝑎𝑟(𝑌5 ) œ 𝑉5 •0
𝜕𝛽
𝜕𝛽

5l0
6

5l0

擬似尤度スコア方程式

𝑞. 𝑒. 𝑑.

25

擬似尤度スコア推定方程式まとめ
・従属変数の真の分布が未知な場合にも，重み𝑽𝒊 と平均𝝁𝒊の関係（任意の関数 𝑣(œ)）さえ適当に与えてやれば，妥当な統計的推論を行える．
・特に，真の分布が指数分布であった場合(𝑉5 = 𝑉𝑎𝑟(𝑌5 ))の推定量は，一般化線形モデルの最尤推定量と一致する．
・推定方程式は，陽に解くことは難しいので，反復重み付
き二乗法など数値計算を行えば解が求まる．
e.g. Newton-Raphson method

擬似尤度スコア方程式

26

一般化推定方程式(GEE)の着想(Zeger, S. L., Liang, K. Y. , 1986)

従属変数が独立でない場合にも妥当な統計的推論を行いたい．
e.g. 土地利用など繰り返し観測データ

GLMの仮定の再確認
仮定①：従属変数が指数分布族に従う．
仮定②：従属変数がそれぞれ独立である．尤度関数 L(𝜃|𝑦) = ∏65l0 𝑓(𝑦5 )仮定②も緩める：𝑪𝒐𝒗(𝒀𝒊 )を考える．
従属変数𝑌5 が独立でない（相関構造を持つ）確率分布を見つけるのは難しいので，仮定①の制約条件を緩めた擬似尤度スコア
推定方程式の考え方を応用して，仮定②の制約条件を取り除くことを考える．
一般化推定方程式27

“working” correlation matrix
𝑽𝒊 = 𝑨𝟏⁄𝟐 𝑹𝒊 (𝜶)𝑨𝟏⁄𝟐 /∅“working” correlation matrix

where

𝑨𝒊 = 𝑑𝑖𝑎𝑔{𝑣(𝝁𝒊 )}𝑹𝒊 ：相関行列であるという要件を満たすn×n対称行列
𝛼： 𝑹を完全に特徴付ける𝑠×1ベクトル∅：スケールパラメータ6

擬似尤度スコア推定方程式

𝑈¤¥ 𝛽 = Œ 𝐷5 ² 𝑉5 •0 𝑆5 = 0
5l0

where

𝜕{𝑎f 5 (𝜃5 )}𝐷5 == 𝐴5 ∆5 𝑋5𝜕𝛽
𝑆5 = 𝑌5 − 𝑎′(𝜃5 )

一般化推定方程式

相関構造を表現できるようにする[パラメーターαを導入]↓
𝑉5 = 𝑣(𝜇5 )/∅𝐴5 = 𝑑𝑖𝑎𝑔{𝑣(𝜇5 )}28

一般化推定方程式：GEE
𝑲

𝑼𝑮𝑬𝑬 𝜷, 𝜶 = Œ 𝑫𝒊 𝑻 𝑽𝒊 •𝟏 𝑺𝒊 = 𝟎
𝒊l𝟏

where

𝜕{𝑎f 5 (𝜃5 )} 𝜕𝜇5
𝑫5 === 𝐴5 ∆5 𝑋5𝜕𝜷𝜕𝜷
𝑽5 = 𝑨𝟏⁄𝟐 𝑹𝒊 (𝜶)𝑨𝟏⁄𝟐 /∅𝑺𝒊 = 𝒀𝒊 − 𝝁𝒊 𝜷 = 𝒀𝒊 − 𝒂′(𝜽𝒊 )𝑨5 = 𝑑𝑖𝑎𝑔{𝑣(𝝁5 )}擬似尤度スコア推定方程式と同様．パラメーターαに依存することのみ異なる．
一般化推定方程式29

GEEの仮定と証明
(i)
(ii)
(iii)

𝐾 𝛼b − 𝛼 = 𝑂¨ 1𝑔𝑖𝑣𝑒𝑛 𝛽 𝑎𝑛𝑑 𝜙𝐾 𝜙Æ − 𝜙 = 𝑂¨ 1𝑔𝑖𝑣𝑒𝑛 𝛽𝜕𝛼b(𝛽, 𝜙)⁄𝜕𝜙 ≤ 𝐻 𝑌, 𝛽 𝑤ℎ𝑖𝑐ℎ 𝑖𝑠 𝑂¨ (1)6𝑈ÒÓÓ 𝛽, 𝛼 = Œ 𝐷5 ² 𝑉5 •0 𝑆5 = 0 は次のように書き直せる．
5l0𝛂, 𝜷, 𝝓依存𝑲} 𝜷
Ð 𝜷, 𝝓Œ 𝑼𝒊 𝜷, 𝜶=𝟎𝒊l𝟏

また，便宜上Ð {𝜷, 𝜙Æ(𝜷)} を定義しておく．
𝜶∗ 𝜷 = 𝜶
一般化推定方程式30

𝜶の漸近正規性の証明
Ð {𝜷, 𝜙Æ(𝜷)}を𝜶のまわりでテーラー展開する．
𝜶∗ 𝜷 = 𝜶∗∑65l0 𝑈5 {𝛽, 𝛼 (𝛽)}

𝐾

=
=

∑6
5l0 𝑈5 (𝛽, 𝛼)𝐾∑65l0 𝑈5 (𝛽, 𝛼)𝐾+∑65l0 𝜕𝑈5 (𝛽, 𝛼)/𝜕𝛼𝐾

𝛼 ∗ − 𝛼 + 𝑂¨ (1)∑65l0 𝜕𝑈5 (𝛽, 𝛼)/𝜕𝛼+𝐾 𝛼 ∗ − 𝛼 + 𝑂¨ (1)𝐾= 𝑷∗ + 𝑸∗ 𝑹∗ + 𝑂¨ (1)
where
∗

𝑷 =

∑6
5l0 𝑈5 (𝛽, 𝛼)𝐾Ð𝜕𝜶𝜷, 𝜙Æ 𝜷𝜕𝜙𝜙Æ − 𝜙Ð 𝜷, 𝜙Æ 𝜷 − 𝜶Ð 𝜷, 𝜙
𝜶
= lim
} →𝝓
𝝓
𝜙Æ − 𝜙𝑸∗ =1
lim Œ 𝑫𝒊 𝑻 𝑽𝒊 •𝟏 𝑺𝒊 = 𝑂¨ (1)
𝐾 6→±
5l0

Ð 𝜷, 𝜙Æ 𝜷𝑹∗ = 𝐾 𝜶= 𝐾

6Ð 𝜷, 𝜙 + 𝜶Ð 𝜷, 𝜙 − 𝜶−𝜶Ð𝜕𝜶Ð 𝜷, 𝜙 − 𝜶 = 𝑂¨ (1){𝜷, 𝜙Æ(𝜷)} + 𝜶𝜕𝜙

一般化推定方程式31

𝜶の漸近正規性の証明（続き）
Ð {𝜷, 𝜙Æ(𝜷)}を𝜶のまわりでテーラー展開する．
𝜶∗ 𝜷 = 𝜶∗∑65l0 𝑈5 {𝛽, 𝛼 (𝛽)}

𝐾
↑

=
=

漸近的に𝑷∗ に一致∑65l0 𝑈5 (𝛽, 𝛼)𝐾∑65l0 𝑈5 (𝛽, 𝛼)𝐾+∑65l0 𝜕𝑈5 (𝛽, 𝛼)/𝜕𝛼𝐾

𝛼 ∗ − 𝛼 + 𝑂¨ (1)∑65l0 𝜕𝑈5 (𝛽, 𝛼)/𝜕𝛼+𝐾 𝛼 ∗ − 𝛼 + 𝑂¨ (1)𝐾= 𝑷∗ + 𝑸∗ 𝑹∗ + 𝑂¨ (1)
where
∗

𝑷 =

∑6
5l0 𝑈5 (𝛽, 𝛼)
𝐾

6

𝑸∗ =

1
lim Œ 𝑫𝒊 𝑻 𝑽𝒊 •𝟏 𝑺𝒊 = 𝑂¨ (1)
𝐾 6→±
5l0

↑

漸近的に正規分布
Ð 𝜷, 𝜙Æ 𝜷𝑹∗ = 𝐾 𝜶= 𝐾

Ð 𝜷, 𝜙 + 𝜶Ð 𝜷, 𝜙 − 𝜶−𝜶Ð𝜕𝜶Ð 𝜷, 𝜙 − 𝜶 = 𝑂¨ (1){𝜷, 𝜙Æ(𝜷)} + 𝜶𝜕𝜙
一般化推定方程式𝜙Æ − 𝜙32

𝜶の一致性の証明：平均と分散∗ (𝛽)}
∑6
𝑈
{𝛽,
𝛼
5
5l0

𝐾

の平均と分散は，擬似尤度スコア推定

方程式の場合に導出したのと同様に導出出来る．

𝐸

𝑉𝑎𝑟∗ (𝜷)}
∑6
𝑈
{𝜷,
𝛼
5
5l0

∗ (𝜷)}
∑6
𝑈
{𝜷,
𝛼
5l0 5

𝐾

𝐾

=𝟎

6

= lim {Œ 𝐷5 ² 𝑉5 •0 œ 𝐶𝑜𝑣(𝑌5 ) œ 𝑉5 •0 𝐷5 ² }
6→±

5l0

一般化推定方程式33

𝜷の漸近正規性の証明} を真値𝜷∗ のまわりでテーラー展開することで，擬似尤𝑼𝑮𝑬𝑬 𝜷度スコア推定方程式の場合と同様に議論できて，} − 𝜷∗
𝑲 𝜷

•𝟏
𝟏 𝝏
𝟏
∗
= −
𝑼𝑮𝑬𝑬 𝜷, 𝜶 (𝜷) |𝜷l𝜷∗
𝑼𝑮𝑬𝑬 𝜷, 𝜶∗ (𝜷)|𝜷l𝜷∗ + 𝑶𝒑 (𝟏)
𝑲 𝝏𝜷
𝑲

∴
where

} − 𝜷∗ ~ 𝑵(𝟎, 𝑱𝟎 •𝟏 𝑱𝟏 𝑱𝟎 •𝟏 )𝑲 𝜷
6

1𝑱𝟎 = lim Œ6→± 𝐾

𝜕𝜇5 ² •0 𝜕𝜇5
𝑉5
𝜕𝛽
𝜕𝛽

1
𝑱𝟏 = lim Œ6→± 𝐾

𝜕𝜇5 ² •0𝜕𝜇5𝑉5 œ 𝐶𝑜𝑣(𝑌5 ) œ 𝑉5 •0
𝜕𝛽
𝜕𝛽

5l0
6

5l0

一般化推定方程式

𝑞. 𝑒. 𝑑.
34

一般化推定方程式の数値計算} 𝜷
Ð 𝜷, 𝝓
一般化推定方程式 ∑𝑲
𝒊l𝟏 𝑼𝒊 𝜷, 𝜶

= 𝟎 は陽に解くことは} と𝜶
Ð の二段階推定𝜷

できないので数値計算で解く．
Modified Fisher scoring algorithm•𝟏𝑲} 𝒋/𝟏 = 𝜷
} 𝒋 − Œ 𝑫𝒊 𝑻 (𝜷} 𝒋) œ 𝑽} 𝒋 ) œ 𝑫𝒋 (𝜷} 𝒋)Û𝒊 •𝟏 (𝜷
𝜷
𝒊l𝟏

where

𝑲

} 𝒋) œ 𝑽} 𝒋 ) œ 𝑺𝒊 (𝜷} 𝒋)Û𝒊 •𝟏 (𝜷œ Œ 𝑫𝒊 𝑻 (𝜷𝒊l𝟏

} (𝜷)}]
Û𝒊 •𝟏 𝜷 = 𝑽𝒊 [𝜷, 𝜶Ð {𝜷, 𝝓𝑽

define 𝐷 = 𝐷0 ² , … , 𝐷6 ²²

, 𝑆 = 𝑆0 ² , … , 𝑆6 ²

Ü5 ：𝑛𝐾×𝑛𝐾 𝑏𝑙𝑜𝑐𝑘 𝑑𝑖𝑎𝑔𝑜𝑛𝑎𝑙 𝑚𝑎𝑡𝑟𝑖𝑥
𝑉𝑡ℎ𝑒 𝑚𝑜𝑑𝑖𝑓𝑖𝑒𝑑 𝑑𝑒𝑝𝑒𝑛𝑑𝑒𝑛𝑡 𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒

𝒁 = 𝑫𝜷 − 𝑺

一般化推定方程式

35

一般化推定方程式まとめ・(回帰係数でなく分散に関心がある場合には使えないが）従属変数が独立でない場合にも，分散の仮定なしで妥当な統計的推論を行えるように，“working” correlation
matrixを導入することで一般化線形モデルを拡張したものである．

擬似尤度スコア方程式

36

参考文献
・ Shannon, C. E., A mathematical theory of communication, Part I, Part II. Bell Syst.
Tech. J., Vol.27, pp.623-656, 1948.
・Wedderburn, R. W., Quasi-likelihood functions, generalized linear models,
and the Gauss—Newton method. Biometrika, Vol.61, No.3, pp.439-447, 1974.
・Zeger, S. L., Liang, K. Y., Longitudinal data analysis for discrete and continuous outcomes. Biometrics, pp.121-130, 1986.
・Liang, K. Y., Zeger, S. L., Longitudinal data analysis using generalized
linear models. Biometrika, Vol.73, No.1, pp.13-22, 1986.
・柳本武美, 推定方程式に基づく推定 最尤法とモーメント法から. 応用統計学,
Vol.24, No.1, pp.1-12, 1995.
・池田思朗, 情報源の構造推定: 確率モデルの構造探索アルゴリズム, Doctoral
Dissertation, 1996.
・甘利俊一(2014). 情報幾何学の新展開. サイエンス社.
・Sopasakis, A. et al., Information metrics for improved traffic model fidelity
through sensitivity analysis and data assimilation. Transportation Research Part B:
Methodological, Vol.86, pp.1-18, 2016.
擬似尤度スコア方程式

37