Proc GAMの紹介

東京大学大学院医学系研究科生物統計学
田中司朗*
松山裕

発表のアウトライン


イントロダクションProc GAMの特徴



散布図の平滑化

GAM （generalized additive model, 一般化加法モデル）


平滑化手法（LOESS, 平滑化スプラインなど）

推定アルゴリズム 自由度の設定
Proc GAMの文法・出力
 解析事例: Kyphosis（脊柱後弯症）データ
Proc GAMで何が出来るか?
散布図の平滑化
 柔軟な交絡調整 結果変数の予測 用量反応関係の把握 ...


Proc GAMの特徴


多機能な回帰分析を行うプロシジャ

結果変数: 連続データ, 2値データ, 計数データ...



平均構造のノンパラメトリックなモデル化

結果のグラフィカルな提示に有用
 比較してみると...
Proc LOESS Proc GENMOD
Proc GAMと他のプロシジャとの比較
結果変数の型
連続データ

パラメトリック平均構造のモデル

カテゴリカルデータ

LOESS
GENMODGAM
LOESSGENMOD
GAM

GAM

GAM

ノンパラメトリック歴史的な流れ

70年代: 一般化線型モデル


カテゴリカルデータのモデルなどを統一的に表現

80年代: ノンパラメトリック回帰の発展
探索的データ解析の一手法として
 散布図の平滑化: 連続量




GAM
Hastie T, Tibshirani R. Generalized Additive Models,
Chapman and Hall: London; 1990.
 古典的なノンパラメトリック回帰をカテゴリカルデータの場合に拡張Melanoma（黒色腫）データ（Houghton D, et al）
コネチカット州がん登録データベースによる記述疫学
1936年から1972年までの悪性黒色腫発生の推移
 結果変数: 10万人あたりの年齢調整発生数
 説明変数: 年代Melanomaデータの散布図
発生数/10万人

年代

回帰分析の出力発生数/10万人

年代

直線では捉えきれない曲線的な傾向
発生数/10万人

年代回帰分析: 線型性の仮定


Anscombe 1973の例
回帰分析の出力は同じ
 回帰係数 残差平方和...
 グラフを正しく要約しているのは(a)だけ
 (b): 非線型 (c), (d): 外れ値
ノンパラメトリック: 線型性の仮定を緩和


単回帰モデルyとxの関係を一次関数として特定
 推定する対象は回帰係数（β0 , β1）yi = β0 + β1xi + εi

ノンパラメトリック回帰のモデル
s(x)は特定しない（何らかの非線型の関数でよい）
 関数の形をデータから推定する


yi= s(xi) + εi

単純なノンパラメトリック回帰（移動平均）
発生数/10万人
k個
この点の予測値を推定したいとき次の点に移る
前後k個の点に枠を作る
次の点に移る
この範囲の平均を取る

年代

移動平均の出力発生数/10万人

なんとなく曲線的な傾向が年代移動平均

重み付き平均の形で表される

w y∑sˆ( x ) =ni =1jijik⎧ 1 （xiがx jのk近傍に含まれるとき）
wij = ⎨
⎩0 （含まれないとき）Proc GAMでは散布図の平滑化が可能
 より性能のよい平滑化手法を採用
LOESS

平滑化スプライン

GAM（一般化加法モデル）
古典的な散布図平滑化の手法をカテゴリカルデータへ拡張 まずは一般化線型モデルについて...
一般化線型モデル


線型モデルやカテゴリカルデータ解析のモデルを統一的に表現

ロジスティック回帰, Poisson回帰...

2つの構成要素から成る
サンプリングモデル
 平均構造のモデル（GAMではこの部分を拡張）サンプリングモデル


結果変数の従う分布は指数型分布族


正規分布, 二項分布, Poisson分布などを含む

平均構造のモデル


結果変数Yiの平均と説明変数Xiとの関係


ただし、E(yi)=µi, g(µi)はリンク関数

g(µi)=Xiβ
線型回帰: g(µi)=µi
 ロジスティック回帰: g(µi)= log[µi/(1-µi)]GAM: 平均構造に加法モデルを仮定


結果変数Yiの平均と説明変数Xiとの関係


ただし、E(yi)=µi, g(µi)はリンク関数

g ( µi ) = s0 + ∑ j =1 s j ( xij )p

sj(x)は特定しない（何らかの非線型の関数でよい）関数の形をデータから推定する

発表のアウトライン


イントロダクションProc GAMの特徴



散布図の平滑化

GAM （generalized additive model, 一般化加法モデル）


いくつかの平滑化手法（LOESS, 平滑化スプラインなど）

推定アルゴリズム 自由度の設定
Proc GAMの文法・出力
 解析事例: Kyphosis（脊柱後弯症）データ
LOESS（Locally Weighted Scatterplot Smoother）


移動平均では

推定する点の近くには1, 遠くには0の重みを与えて平均を取ったらうまくいった

LOESSの2つの工夫
最小二乗法による推定
 よりよい重み 局所重み付き（Locally Weighted）最小二乗推定LOESSの最小二乗推定

以下の重み付き残差平方和を最小にする推定量
wijは重み（後述）

RSS = ∑i =1 wij [ yi − [ β 0 j + β1 j ( xi − x j )]]n2

単回帰のモデルを当てはめる

LOESSの重み

数式で表すと複雑だが...
3 3⎧ ⎡⎤⎛⎞|−|
xx32ij⎪⎪ ⎢1 − ⎜⎟ ⎥ （xiがx jのk近傍に含まれるとき）
wij = ⎨ 5 ⎢ ⎜ max | x − x | ⎟ ⎥ij ⎠⎦⎪ ⎣ ⎝⎪⎩ 0（含まれないとき）

LOESSの重み

重みをグラフにしてみると


近くには大きい重み、遠くには小さい重みk個の点が含まれる
バンド幅

重みのイメージ（移動平均）
発生数/10万人

この範囲の平均を取る

年代

重みのイメージ（LOESS）
発生数/10万人

重み付き最小二乗推定

拡大

年代

局所での推定のイメージ（移動平均）

平均をとる
s(xi)を推定

局所での推定のイメージ（LOESS: 単回帰）

単回帰を当てはめる

s(xi)を推定

LOESSの出力発生数/10万人

年代移動平均はk個の点の平均

k=1のときs(xi)の推定量は、その点自身（自由度はn）



バイアスはないが、推定精度が悪い

k=nのときs(xi)の推定量は、全体平均（自由度は1）



推定精度はよいが、バイアスがある

精度とバイアスのトレードオフはバンド幅で決まる

k=1のとき（バイアスはないが精度が悪い）
発生数/10万人

年代

k=nのとき（精度はよいがバイアスがある）
発生数/10万人

年代

バンド幅の意味

バンド幅を決めるkは平滑化パラメータ
平滑化の程度を規定する
 推定精度とバイアスのトレードオフ
 解析者が指定するチューニングパラメータ発表のアウトライン



イントロダクションProc GAMの特徴



散布図の平滑化

GAM （generalized additive model, 一般化加法モデル）


いくつかの平滑化手法（LOESS, 平滑化スプラインなど）

推定アルゴリズム 自由度の設定
Proc GAMの文法・出力
 解析事例: Kyphosis（脊柱後弯症）データ
三次のスプライン関数


区分多項式の一種
範囲を適当に区分する 区間ごとに三次関数を当てはめる
 ただし, 区間の繋ぎ目（節点）では連続で滑らか


三次関数
節点

平滑化スプライン以下のペナルティ付き残差平方和を最小とする推定量

RSS = ∑i =1[ yi − s ( xi )] +n2λ22[s''(x)]dx∫

xの実現値を節点とした三次の自然スプラインが解平滑化スプライン以下のペナルティ付き残差平方和を最小とする推定量

RSS = ∑i =1[ yi − s ( xi )] +n2

λ22[s''(x)]dx∫

xの実現値を節点とした三次の自然スプラインが解λは平滑化パラメータ平滑化スプラインの出力（λ≒∞ ）発生数/10万人

年代

平滑化スプラインの出力（λ≒0）
発生数/10万人

年代

平滑化スプライン以下のペナルティ付き残差平方和を最小とする推定量

RSS = ∑i =1[ yi − s ( xi )] +n2

λ22[s''(x)]dx∫

第一項はデータへの当てはまりをよくする基準
 第二項は曲線の変動を小さくする基準λは2つのバランスを規定する

自由度


平滑化パラメータは手法により異なるLOESSではk 平滑化スプラインではλ

自由度により統一的に表現

実質的なパラメータの数に変換自由度

回帰分析では, 最小二乗推定量は行列の形で表せる


Aはハット行列yˆ = X ( X t X ) −1 X t y= Ay

自由度（パラメータの数）は行列Aのランク

単回帰における帰無仮説「β1=0」のt検定


モデルの自由度は2（切片と傾き）
残差の自由度はn-2

自由度LOESS, 平滑化スプライン推定量もこの形に表せる


Aは平滑化パラメータと一対一に対応

sˆ = Ay

自由度
LOESSではAtAのトレース（対角要素の和） 平滑化スプラインではAのトレース

こう決めておけば平滑化パラメータを統一的に扱える自由度の設定

解析者が指定

扱っているデータ解析の背景からデータから推定GCV（一般化クロスバリデーション）法以下の基準を最小にする自由度を選ぶn∑i =1[ yi − sˆ( xi )]2
nGCV =[n − trace( A)]

2

自由度5
発生数/10万人

年代自由度7
発生数/10万人

年代

自由度9
発生数/10万人

年代

自由度11
発生数/10万人

年代

自由度13
発生数/10万人

年代

自由度15
発生数/10万人

年代

自由度17
発生数/10万人

年代

自由度19
発生数/10万人

年代

自由度21
発生数/10万人

年代

自由度23
発生数/10万人

年代

自由度25
発生数/10万人

年代

重ねてみると発生数/10万人

年代

説明変数が複数ある場合は?
1. 多次元の関数を考える


実質的なパラメータ数が加速的に増加

yi = s( xi1, xi 2 ,..., xip ) + ε i

2. 加法モデルの制約を入れる
1つ1つの説明変数の効果は足し算で効く 交互作用はない

yi = s0 + ∑ j =1 s j ( xij ) + ε ip

薄板スプライン
二変量を扱える平滑化手法
 以下のペナルティ付き残差平方和を最小とする推定量
222222⎡2⎛ ∂ s ⎞ ⎛∂ s⎞ ⎤λ ⎛⎜ ∂ s ⎞⎟n⎟⎟ + ⎜ 2 ⎟ ⎥ dx1dx2
RSS = ∑i =1[ yi − s ( x1i , x2i )] + ∫ ⎢⎜ 2 ⎟ + 2⎜⎜
2 ⎢⎝ ∂x1 ⎠∂x1∂x2 ⎠ ⎜⎝ ∂x2 ⎟⎠ ⎥⎝⎦⎣

Proc GAMで指定可能推定アルゴリズム


二つとも推定は最小二乗法
LOESSだと局所重み付き...
 平滑化スプラインだとペナルティ付き...


最小二乗解を求めるには?
一度に全ての説明変数について推定するのは大変 ここで加法モデルが効いてくるバックフィッティングアルゴリズム


特徴
説明変数ごとに解を更新する
 順番に残差に対して当てはめを行う

例えば, k番目の関数skの更新ではsk以外に関する残差を求める
 Rikを結果変数, xikのみを説明変数として解を求める
Rik = yi − sˆ0 − ∑ j =1 sˆ j ( xij ) − ∑ j = k +1 sˆ j ( xij )k −1p

バックフィッティングアルゴリズム

ステップ1（初期値を与える）
ˆs0 (1) = Y , sˆ1(1) = sˆ2 (1) = ... = sˆ p (1) = 0



ステップ2（反復計算）(1)( 2)
R1にx1を当てはめて, sˆ1 をsˆ1 に更新
(1)

( 2)

R2にx2を当てはめて, sˆ2 をsˆ2 に更新
M
全ての説明変数が終わ ったら次のループへ

ステップ3（収束の判定）
残差平方和がこれ以上減少しなくなったらストップ

Proc GAMの文法
PROC GAM < option > ;CLASS variables ;MODEL response = covariates < /options > ;SCORE data=dataset out=dataset ;OUTPUT <out=dataset> <keyword> ;BY variables ;ID variables ;FREQ variable ;Proc GAMの文法
PROC GAM < option > ;CLASS variables ;MODEL response = covariates < /options > ;SCORE （予測用のデータセットを指定）
OUTPUT <OUT=dataset> <keyword> ;BY variables ;ID （ID用の変数を指定）
FREQ variable ;MODELステートメント
MODEL response = covariates < /options > ;当てはめるモデルを指定する
 平滑化手法
 自由度
 結果変数の分布

MODELステートメント
MODEL response = SPLINE (covariate1, DF=3)
PARAM (covariate2, covariate3)/ DIST=BINOMIAL;
平滑化手法
 自由度
 結果変数の分布
MODELステートメント
平滑化手法


自由度




PARAM, LOESS, SPLINE, SPLINE2（薄板スプライン）
直接指定, オプション（METHOD=GCV）

結果変数の分布GAUSSIAN, BINOMIAL, BINARY, GAMMA,IGAUSSIAN, POISSON
 リンク関数はカノニカルリンク
OUTPUTステートメント
OUTPUT <OUT=dataset> <keyword> ; 予測値・信頼区間・回帰診断統計量などを出力
 キーワード: 変数名
PRED : 予測値
 UCLM : 上側信頼区間, LCLMだと下側 ADIAG : テコ比統計量
 RESID : 残差統計量
 STD: 標準偏差 ALL
: 上記の変数全てを出力
OUTPUTステートメント
OUTPUT OUT=dataset ALL;
 予測値・信頼区間・回帰診断統計量などを出力
 キーワード: 変数名
PRED : P̲response, P̲covariate1
 UCLM : UCLM̲covariate1
 ADIAG : ADIAG̲covariate1
 RESID : R̲covariate1: STD̲covariate1 （スプライン, LOESSではSTDP̲） STD
ALL: 上記の変数全てを出力

平滑化する説明変数の扱いについて一次の項と二次以上に分けて出力 例えば, 年齢の効果は
回帰係数×Age + P̲Age 一次の効果: 回帰係数×Age 二次以上の効果: P̲Age

上側信頼区間回帰係数の上側信頼区間×Age + UCLM̲Age

GCV基準により選択された自由度
1を引いて出力Kyphosis（脊柱後弯症）データ（Bell D, et al）
目的: Kyphosisの予後因子の探索 対象：椎弓切除術を受けた子供83人
 エンドポイント: Kyphosisの有無


予後因子


Age: 年齢

その他の変数 StartVert: （手術を脊椎のどの部位から行ったか） NumVert: （手術を行った部位の数）

Kyphosisデータの入力
data Kyphosis;
input Age StartVert NumVert Kyphosis @@;datalines;
71 5 3 0
158 14 3 0 128 5 4 1
2 1 5 0
1 15 4 0 1 16 2 0
61 17 2 0
37 16 3 0 113 16 2 0
59 12 6 1
82 14 5 1 148 16 3 0
18 2 5 0
1 12 4 0 243 8 8 0・・・;

2値データだと…
Kyphosisの有無

年齢（月）SASのプログラム


結果変数




Kyphosis

説明変数
平滑化（自由度3）: Age

線型:StartVert, NumVert

入出力データセット
解析用データセット: Kyphosis
 出力データセット: Kyphosis̲out2
 予測用入力データセット: Kyphosis̲in 予測用出力データセット: Kyphosis̲out1
SASのプログラム
PROC GAM DATA=Kyphosis;
MODEL Kyphosis = PARAM(StartVert NumVert)
SPLINE(Age, DF=3)/ DIST=BINOMIAL;SCORE DATA=Kyphosis̲in out=Kyphosis̲out1;
OUTPUT OUT=Kyphosis̲out2 ALL;
RUN;

SASの出力
大切なのはこの3つ

Parameter Estimates
パラメトリックで指定した変数の回帰係数の出力 平滑化した変数（Linear(Age)）は一次の部分だけFit Summary of Smoothing Components
平滑化した変数の自由度やGCV基準の当てはまり
 自由度は2


一次の項は除いて考えるAnalysis of Deviance（デビアンスを用いた分散分析）
平滑化した変数の検定
 Ageは有意OUTPUTステートメントの出力

SCOREステートメントの出力

Proc GENMODで解析した結果
StartVert, NumVertの結果は同様の結果
 Ageは有意ではない


GENMODの結果をグラフにするとKyphosisの有無

年齢（月）GAMの結果をグラフにするとKyphosisの有無*破線は95%信頼区間

年齢（月）

自由度2
Kyphosisの有無

年齢（月）

自由度5
Kyphosisの有無

年齢（月）GCV基準を用いてみるとKyphosisの有無*GCV基準による自由度は2.27

年齢（月）まとめ: 発表内容


ノンパラメトリック回帰

散布図の平滑化

GAMの紹介

平滑化手法（LOESS, 平滑化スプラインなど）

推定アルゴリズム 自由度の設定
Proc GAMの文法・出力
 Kyphosisデータ解析の紹介
まとめ: 注意点
信頼区間・自由度などをグラフに示す
 自由度の設定


平滑化の程度によって印象がかなり異なる
 いくつかの値を試してみる


平滑化した変数の出力


一次の項と二次以上に分けて出力

参考文献




マニュアル・教科書
 SAS Institute. SAS/STAT® User’s Guide Version9.1, SAS Institute, Cary,NC, U.S.A.; 2004.
 Hastie T, Tibshirani R. Generalized Additive Models, Chapman and Hall:
London; 1990.
 Ruppert D, Wand MP, Carroll RJ. Semiparametric regression, Cambridge
University Press: Cambridge; 2003.
原著
 Houghton AN, Flannery J, Viola MV. Malignant Melanoma in Connecticut
and Denmark. International Journal of Cancer 1980; 25: 95-104.
 Bell D, et al. Spinal Deformation Following Multi-Level Thoracic andLumbar Laminectomy in Children. 1989; Submitted for publication.
 Anscombe. Graphs in Statistical Analysis. American Statistician 1973; 27:
17-21