東京大学 | 大学院医学系研究科 | 国際保健学専攻 | 人類生態学教室TOP | 2nd
社会統計学第10回
「回帰分析，重回帰分析，ロジスティック回帰分析」（2001年11月29日）
トップ | 更新情報 | 研究と教育 | 業績 | 計算機 | 写真 | 枕草子 | 著者 | 目安箱 | 書評 | 社会統計学目次
最終更新：
2001年12月7日 金曜日 13時57分
全部を一括して読む | 前回へ
講義概要
- 回帰と重回帰：
- ▼1984年から1993年までのプロ野球の１試合平均入場者数（単位：千人）の推移は下表のようになっている。
| 年次 | 84 | 85 | 86 | 87 | 88 | 89 | 90 | 91 | 92 | 93 |
| セリーグ | 28 | 29 | 29 | 31 | 31 | 31 | 31 | 32 | 35 | 34 |
| パリーグ | 13 | 12 | 16 | 18 | 21 | 23 | 22 | 24 | 24 | 24 |
- ▼これを見ると，途中から伸び悩んでいるが，ある程度直線的に増加しているように見える。このように，ほぼ比例関係にある量的なデータ群をうまく代表する直線を求めるには，「各データの点から直線までのずれの大きさの合計」を最小にすればよい。この場合，「年次」が予め決まっている値であるのに対して，入場者数は測定値であり，誤差を含む可能性があるので，「データ点から直線までのずれ」を評価するには，データ点と直線の最短距離よりも，データ点から直線に垂直に下ろした線分の長さの二乗和を使う方がよい。この，ずれを最小にする直線を，「年次を独立変数，入場者数を従属変数とする回帰直線」と呼ぶ（計算方法は，数学的には検量線の求め方（pdf形式）を参照されたい）。
- Rでは，nenji<-c(84:93); central<-c(28,29,29,31,31,31,31,32,35,34)として，plot(central‾nenji); z<-lm(central‾nenji); abline(z)とすれば上の図が描かれる。なお，plot(nenji,central)とplot(central‾nenji)は同義である。
- ▼一般に，２組あるいはそれ以上の対で観測されるデータが与えられていて，そのうちの１つの変数（従属変数，または目的変数という）の変動の様子が残りの変数（独立変数，または説明変数という）の変動の仕方と，ある一定の関係を保ちながら変化していることがある。この場合，従属変数の変動を捉えるのに独立変数の変動に関する情報が利用できる。この「一定の関係」を求めることを回帰分析という。関係が直線のとき線型回帰分析といい，lm()が線型回帰分析を行う関数である。
- ▼独立変数が複数あって，それらによって同時に従属変数の変動を説明しようとするとき，とくに重回帰分析というが，同じ関数で実行できる。例えば，独立変数が年齢(変数名ageとする)と体重(変数名wtとする)，従属変数が収縮期血圧(変数名SBPとする)であるなら，その重回帰分析をするには，lm(SBP‾age+wt)とすればよい。
- ロジスティック回帰
- ▼前回問題とした，３つ以上のカテゴリカル変数の関係を同時に調べるのにも，この重回帰の考え方が使える。しかし，カテゴリカル変数の場合，値の範囲が有限であるため，単純な直線関係にはならない。そこで，一般化線型モデルを用いる。関数名はglm()である。独立変数に量的な変数が含まれ，従属変数がカテゴリカル変数の場合，カテゴリカル変数があるカテゴリになる確率pを，log(p/(1-p))という形で無限小から無限大までの範囲をとるように変換して（ロジット変換と呼ぶ），一般化線型モデルを適用するとできる分析がロジスティック回帰分析である。
- 例：
- ▼前回から計算練習に使っているuser.dat（エコポイント・チェックを使って作った仮のデータ）でQ4とQ10の関係には年齢と性別の影響が両方あったので，年齢を量的変数として，NQ4がTRUE（いつもまたは大体買い物袋を持参している）の人はFALSEの人に比べて，どれくらいNQ10がTRUE（いつもまたは大体マイカーを避けて公共交通を利用している）の割合が高いかを考える。
- ▼NQ4がTRUEである割合をpとすると，log(p/(1-p))=β0+β1AGE+β2SEX+β3NQ10というモデルが考えられる。TRUEを1，FALSEを0とすると，他の変数を調整した上でのNQ10の効果は，
- NQ10がTRUEの場合のNQ4がTRUEである割合をp1とすれば，
- log(p1/(1-p1))=β0+β1AGE+β2SEX+β3
- NQ10がFALSEの場合のNQ4がTRUEである割合をp0とすれば，
- log(p0/(1-p0))=β0+β1AGE+β2SEX
- となり，辺々引けば
- log(p1/(1-p1))-log(p0/(1-p0))=β3
- 結局この左辺は，log[(p1/(1-p1))/(p0/(1-p0))]となる。これは対数オッズ比である。つまり，AGEとSEXの影響を調整した上でNQ10がTRUEである人はFALSEの人に比べて，exp(β3)倍，NQ4がTRUEになりやすいと言える。普通の言葉で書き直せば，「性別と年齢の影響を調整した上で，いつもまたは大体マイカーを避けて公共交通を利用している人は，そうでない人に比べて，いつもまたは大体買い物袋を持参する傾向がexp(β3)倍大きい」ということである。なお，β0〜β3は，一般化線型モデルの係数として推定される。
- Rでのプログラム：
- result<-glm(NQ4‾NQ10+AGE+SEX, family=binomial(logit)); summary(result)
- 出力：
- Call:
- glm(formula = NQ4 ‾ NQ10 + AGE + SEX, family = binomial(logit))
- Deviance Residuals:
| Min | 1Q | Median | 3Q | Max |
| -0.9472 | -0.7143 | -0.5803 | -0.4478 | 1.9311 |
- Coefficients:
| Estimate | Std. Error | z value | Pr(>|z|) |
| (Intercept) | -2.2496 | 0.6354 | -3.540 | 0.000399 *** |
| NQ10 | 0.2969 | 0.5291 | 0.561 | 0.574639 |
| AGE | 0.2767 | 0.3169 | 0.873 | 0.382463 |
| SEX | 0.7372 | 0.5109 | 1.443 | 0.149070 |
- ---
- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
- (Dispersion parameter for binomial family taken to be 1)
- Null deviance: 100.08 on 99 degrees of freedom
- Residual deviance: 96.64 on 96 degrees of freedom
- AIC: 104.64
- Number of Fisher Scoring iterations: 3
- ▼とりあえずは，Coefficients:のEstimateにβ0〜β3が得られ，Pr(>|z|)にそれぞれが0である確率が得られている，と判断すればよい。
- ▼モデル自体の説明力に関しては，AIC（赤池の情報量基準）を見る。一般に，未知パラメータ数kの適当なモデルMを仮定したときの最大対数尤度をMLL(M)と書くと，
- AIC(M)=-2MLL(M)+2k
- と定義される。AICは，モデルの当てはまりの悪さを示す値であるから，複数のモデルのAICを比べた場合は，AICが小さい方を採択するわけである。
フォロー
- ロジスティック回帰で出てきたlogは，高校数学で出てきたのと同じ？
- ▼その通り，対数のことです。統計学でlogが出てきて底が明示されていない場合は，自然対数と考えていいです。
- 一般化線型モデルとは？
- ▼線型モデルを一般化したものです。線型モデルは，普通の量的変数の間の直線回帰を指します。これを一般化すれば，t検定，分散分析，共分散分析，回帰分析，判別分析，正準相関分析などをすべて共通の数学モデルで扱うことができ，それを一般化線型モデルといいます。英語ではGeneral Linear Modelというので，SASでもPROC GLMですが，Rでの関数名もglmです。
- ▼基本的には，Y=β0+βX+εという形で表されます。Yが従属変数群，Xが独立変数群，β0が切片群，βが係数群です。
- Rの日本語版のようなものはないのですか？
- ▼Rには日本語版はありません。ただ，統計用語は英語の方が曖昧さがありませんし，Rのマニュアルはかなり日本語訳されていますから，現状でもそれほど問題はないと思います。
- あるデータをもっているときどんなふうに分析したらいいということはどうやって決める？
- ▼一般的な統計的分析の手順としては，第１回で示したことの繰り返しになりますが，もう少し詳しく書くと，
- （研究全体の）目的を明確にして調査を行う
- データ化（エディティング，コーディング，データ入力）
- データの図示（幹葉表示など）（第３回を参照）
- 代表値（平均値，中央値，分散，四分位偏差など）の計算（必要なら区間推定）（第４回，第５回を参照）
- 作業仮説を立てる
- 複数のグループが予め決めてあるデザインの場合は，グループ間で各変数の値に違いがあるかどうかを検定する。注目する変数がカテゴリカルな変数なら，クロス集計表を作ってフィッシャーの正確な確率を計算するとか，複数カテゴリ間での比率の差の検定を行うし（第６回を参照），量的な変数なら，グループが２つだけの場合はt検定や順位和検定や符号付順位検定などを行うし（第11回を参照），グループが複数のときは一元配置分散分析と多重比較を行う（第12回を参照）。
- 複数の変数間の関連があるかどうかを検定する。
- 従属関係がない場合は，カテゴリカル変数ならクロス集計表を作って独立性のカイ二乗検定をしたり，フィッシャーの正確な確率を求めたり（第7回を参照），一致度の指標や関連性の指標を求めたり（第8回を参照），それらを攪乱要因によって層別したりするし（第9回を参照），量的な変数ならピアソンの積率相関係数を求めたり，スピアマンやケンドールの順位相関係数を求めたりする（第11回を参照）。
- 従属関係がある場合は，変数がカテゴリカル変数であるか量的な変数であるかによって，使える分析方法は異なる。簡単にまとめると，以下のようになる（なお，ライブラリの関数を使うには，library(ライブラリ名)とする）。
| 従属変数 | 独立変数 | 分析法（Rでの関数名） | 参照回など |
|---|
| カテゴリ変数１つ | カテゴリ変数または量的変数，１つまたは複数 | ロジスティック回帰分析(glm)，判別分析(MASSライブラリのldaまたはmdaライブラリのmda) | 第10回（ロジスティック回帰のみ） |
| 量的変数１つ | カテゴリ変数１つ | 一元配置分散分析(oneway.testまたはkruskal.test)と多重比較(pairwise.t.test) | 第12回 |
| 量的変数１つ | カテゴリ変数複数 | 多元配置分散分析(anova) | 第13回 |
| 量的変数１つ | カテゴリ変数１つ（クラス変数）と量的変数１つまたは複数 | 共分散分析(lmまたはaovで実行できる) | 第12回 |
| 量的変数１つ | 量的変数複数 | 重回帰分析(lm) | 第10回 |
| 量的変数１つ | カテゴリ変数複数 | ダミー変数を使った重回帰分析(lm) | 第13回 |
| 量的変数複数 | カテゴリ変数複数 | 多変量分散分析(manova) | 触れない |
| 量的変数複数 | 量的変数複数 | 正準相関分析(mvaライブラリのcancor) | 触れない |
- 可能なデータについては，次のような，より進んだ解析を実行する。
- 複数の変数間の関係を，少数の直交する軸に集約して見たい場合，因子分析(mvaライブラリのfactanal)や主成分分析(mvaライブラリのprcomp)を実行する。
- データ間の類縁関係を調べたいとき，クラスタ分析(mvaライブラリを使い，データ間の距離行列を計算するにはdistを用い，クラスタの連結をするにはhclustやkmeans，デンドログラムを書かせたいときはplot.dendrogramなど)を実行する。
- 打ち切りのある間隔データを分析したいときは生存時間解析(survivalライブラリを使い，Survで生存時間を示すオブジェクトを生成し，カプランマイヤ法ならsurvsumやsurvfit，加速モデルならsurvreg，コックス回帰ならcoxph，条件付きロジスティック回帰ならclogitなど)を実行する。
- 時系列で起こる事象の変化について分析したい場合は，ある時点での状態が前の時点での状態から独立でないので，時系列解析と呼ばれる特別な手法を用いる。RではtsライブラリでARIMAモデル(arima0)などがサポートされている。
- 複数の変数と潜在変数を含めたネットワークを分析したい場合は，パス解析や共分散構造解析などの手法がある。
- という感じになります。上の表に出てきた，量的な変数に関連する部分を，この講義の残りの回で説明する予定です。
全部を一括して読む | 次回へ