| ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データであ... |
|
| 投稿者：tosh 2018/05/09(Wed) 11:39 No. 22537 |
ロジスティック回帰分析を行う際、説明変数として用いたい変数が、３種類以上の値を含む名義尺度データであった場合に、それを複数のダミー変数に置き換える方法が、いろいろあることについて悩んでいます。
インターネットで調べられた限りですが、おもに次の(1)〜(3)の方法があるようでした。
これらはどのように使い分ければよいでしょうか。
-----------------
説明変数Ａに、"x","y","z" の３種類の名義が含まれている場合、
(1)ダミー変数 { d1, d2, d3 } を作り、変数Ａが、
"x" なら { 1, 0, 0 }
"y" なら { 0, 1, 0 }
"z" なら { 0, 0, 1 }
に置き換える方法。
(2)ダミー変数 { d1, d2 } を作り、変数Ａが、
"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { 0, 0 }
に置き換える方法。
(3)ダミー変数 { d1, d2 } を作り、変数Ａが、
"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { -1, -1 }
に置き換える方法。
-----------------
上に挙げたうち、どれか１つの方法を紹介しているサイトはいくつか見つかるのですが、それぞれの方法の比較や、どう使い分ければよいか、という点についての情報がなかなか見つからず、質問させていただきました。
また、この問題が体系的に解説されている書籍などで、もしおすすめのものがありましたら、ご紹介いただければ幸いです（それぞれの方法の結果の解釈、オッズ比への逆変換の方法、ステップワイズ法で変数を増減する場合にダミー変数はどう扱うのか、などについてもよく知りたいと思っています）。 |
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：後医は名医 2018/05/09(Wed) 12:52 No. 22538 |
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：tosh 2018/05/10(Thu) 12:31 No. 22541 |
ご返信ありがとうございます。
多重共線性は、互いに相関が強い複数の変数が説明変数の中に含まれていると、起こる問題ですよね。
（１）の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…？
（追記）
すみません、次のように修正しました
修正後「 d1 が 1 なら d2 + d3 は必ず0に決定されてしまうので、」
修正前「 d1 が 1 なら d2 とd3 は必ず0に決定されてしまうので、」 |
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：青木繁伸 2018/05/10(Thu) 22:50 No. 22542 |
> １）の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…？
そういうこと
R では，変数Ａが名義尺度であるということを
FA = factor(A)
としておけば，
glm(result ~ FA)
でちゃんと処理してくれます。
変数 A は 2 個のダミー変数で表現できます。
3 個以上の変数を使っても，R では 3 個目からは分析に使われません。他の統計ソフトではエラーメッセージを吐いて死んでしまうこともあるでしょう。
独立な2変数変数であれば，どんな値であっても(0 と 1 でなくても -6.23123 と 578.21234 でもなんでも），その変数に対する係数は見かけは違いますが，その係数を使って計算される予測値は全く同じになります（そうでないと困ります。よい結果が出るような変数値のセットがある！なんてことになると，客観性がなくなります）。
ではなぜ，（3）のような変数を作らない！！かというと，出てくる値が3種類になるからです。0と1だと計算も簡単だしわかりやすい。
それにしても 0,1,-1 を与えるというのは，どんなメリットがあるのかなあ？ |
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：tosh 2018/05/11(Fri) 16:38 No. 22543 |
ご返信、ありがとうございます。
(3)の意義については、この方法が載っていたpdf資料の内容を後でアップいたします。
（私のパソコンのDLフォルダの整理が悪く、見つかり次第）
通常は (2) の置き換えを行い、(1)はダメ、ということですね。
さっそく、R で、名義尺度の変数を factor型 にし、glm メソッド にかけてみました。
（スクリプトは、別に投稿します）
このスクリプトは次のようなロジスティック回帰分析を行うつもりで書きました。
・データセット: birthwt 低体重出生とそのリスク因子のデータ
・目的変数: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
・説明変数: 母親の年齢、母親の体重、母親の人種、...etc
母親の人種 race {1:白人 2:黒人 3:その他 } が名義尺度なので、スクリプトの中であらかじめ次のように変換しています。
# race を factor型 に変換する
DF1$race <- factor(DF1$race)
glmメソッド の summary をみると、以下に引用したように race2, race3 という２つの変数が生成されているのを確認しましたが、これらは race { 1:白人 2:黒人 3:その他 } をどのように変換した変数なのかを確認する方法がわかりません。
---------------------------------------------------------------
# 結果表示
summary(out1.glm)
Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF1)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8946 -0.8212 -0.5316 0.9818 2.2125
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.48062 1.19689 0.402 0.68801
age -0.02955 0.03703 -0.798 0.42489
lwt -0.03397 0.01524 -2.229 0.02580 *
race2 1.27226 0.52736 2.413 0.01584 *
race3 0.88050 0.44078 1.998 0.04576 *
smoke 0.93885 0.40215 2.335 0.01957 *
ptl 0.54334 0.34540 1.573 0.11571
ht 1.86330 0.69753 2.671 0.00756 **
ui 0.76765 0.45932 1.671 0.09467 .
ftv 0.06530 0.17239 0.379 0.70484
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.28
Number of Fisher Scoring iterations: 4
---------------------------------------------------------------
ここでの {race2, race3} は、たとえばですが、内部的に
race=1 なら { 1, 0 }
race=2 なら { 0, 1 }
race=3 なら { 0, 0 }
のように変換されているのだと思いますが、それをどのように確認できるのでしょうか？
str(out1.glm)
としてみたりもしましたが… |
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：tosh 2018/05/11(Fri) 16:40 No. 22544 |
以下、スクリプトです。
----------------------
###################################################
#
# ロジスティック回帰分析例
#
# https://www.gixo.jp/blog/3200/
# より一部改変
#
library(MASS)
str(birthwt)
# *** 目的変数 ***
# low:新生児の体重が2.5kg未満か否か(2.5kg未満が1)
#
# *** 説明変数 ***
# age:母親の年齢
# lwt:母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke:母親の喫煙有無(喫煙有りが1)
# ptl:母親の早産経験の有無(経験有りが1)
# ht:母親の高血圧の有無(有りが1)
# ui:母親の子宮神経過敏の有(有りが1)
# ftv:妊娠第１期における医師の訪問回数
#
# *** 使わない ***
# bwt:新生児の体重(グラム単位)
# ■ 前処理
# bwtは使わないので除く
DF1 <- birthwt[,1:9]
# lwt をポンド単位からグラム単位に変換する
DF1$lwt <- DF1$lwt * 0.454
# race を factor型 に変換する
DF1$race <- factor(DF1$race)
str(DF1)
# ■ ロジスティック回帰
out1.glm <- glm(
low ~ . ,
family = binomial(link = "logit"),
data = DF1)
# 注: モデル式 low ~ . は、目的変数に low、説明変数に low 以外の全ての変数を使うという意味
# 結果表示
summary(out1.glm)
# 各説明変数の係数を指数変換し、オッズ比を表示
exp(out1.glm$coefficients) |
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：青木繁伸 2018/05/11(Fri) 21:30 No. 22545 |
例としてあげられたデータ birthwt$race は整数型で 1, 2, 3 をとる変数でしたが，factor(birthwt$race) で factor型の変数にすると，1, 2, 3 というレベルを持つ factor 変数になります。> DF1$race
[1] 2 3 1 1 1 3 1 3 1 1 3 3 3 3 1 1 2 1 3 1 3 1 1 3 3 1 1 1 2 2 2 1
途中省略
[161] 1 3 3 3 2 1 3 3 1 1 2 2 2 3 3 1 1 1 1 2 3 3 1 3 1 3 3 2 1
Levels: 1 2 3 ← ここ
このような factor 変数を多変量解析などに使用すると，一番最初の Level が基準値 0 になります。つまり，必要なダミー変数において，取る値の全てが0
一般的には，変数を factor 変数に変換すると，"辞書順の最初の値が基準値 0 になります。
詳しくは factor 関数のオンラインヘルプを参照してもらいたいですが，ちょっと見以下のような結果になりますが，訳ワカメに陥らないようにしっかり理解してください。
> sex = c("male", "female", "U.K.", "male", "male", "female")
> factor(sex, labels=c(1,2,0))
[1] 2 1 0 2 2 1
Levels: 1 2 0
> factor(sex, levels=c("male", "female", "U.K."))
[1] male female U.K. male male female
Levels: male female U.K.
> factor(sex, levels=c("male", "female", "U.K."), labels=c(1,2,0))
[1] 1 2 0 1 1 2
Levels: 1 2 0
さて，このようなことから，R が作っているダミー変数を自分で作るとすると以下のようになります> DF2 = DF1
> DF2$race = NULL # factor である race を削除
> Race = birthwt$race # 元のデータは整数型 1, 2, 3
> DF2$Race2 = (Race == 2)+0 # Race が 2 のときのみ整数値 1，それ以外なら 0
> DF2$Race3 = (Race == 3)+0 # Race が 3 のときのみ整数値 1，それ以外なら 0
> head(DF2) # Race = 1 のときは DF2$Race2 = DF2$Race3 = 0
low age lwt smoke ptl ht ui ftv Race2 Race3
85 0 19 82.628 0 0 0 1 0 1 0
86 0 33 70.370 0 0 0 0 3 0 1
87 0 20 47.670 1 0 0 0 1 0 0
88 0 21 49.032 1 0 0 1 2 0 0
89 0 18 48.578 1 0 0 1 0 0 0
91 0 21 56.296 0 0 0 0 0 0 1
> out2.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF2)
> summary(out2.glm)
Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF2)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8946 -0.8212 -0.5316 0.9818 2.2125
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.48062 1.19689 0.402 0.68801
age -0.02955 0.03703 -0.798 0.42489
lwt -0.03397 0.01524 -2.229 0.02580
smoke 0.93885 0.40215 2.335 0.01957
ptl 0.54334 0.34540 1.573 0.11571
ht 1.86330 0.69753 2.671 0.00756
ui 0.76765 0.45932 1.671 0.09467
ftv 0.06530 0.17239 0.379 0.70484
Race2 1.27226 0.52736 2.413 0.01584
Race3 0.88050 0.44078 1.998 0.04576
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.28
Number of Fisher Scoring iterations: 4 ということで，DF1 を使ったときと同じ結果になりました。
なお，reorder という関数がありますが，この存在を知ると，きっと悪夢を見ます（^_^;）
複雑怪奇な変換規則を理解するよりは，自分で「適切な」ダミー変数を作成する方が幸せかも知れません。
|
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：tosh 2018/05/12(Sat) 04:06 No. 22547 |
factor関数をテストしてみると、自動的に定まる Levels の順序が独特ですね…。
> factor(c(100,1,3,1,-5))
[1] 100 1 3 1 -5
Levels: -5 1 3 100 → 数値は昇順なのはまあ分かりますが、
> factor(c("b","a","A","c","A","b","C"))
[1] b a A c A b C
Levels: a A b c C → 文字列は辞書順だが小文字が先…？
最初のLevelがどの値になるか、いまいち確信が持てないので、ご説明にあるように、名義尺度をfactor型に変換する際は、 Levels パラメータ と Labels パラメータをちゃんと与えて、最初のLevelになる値を明示すればいい、と考え、次のように進めました。
> DF3 = DF1
> DF3$race = NULL
> DF3$race = factor(birthwt$race,levels=c(1,2,3),labels=c("白人","黒人","その他"))
> out3.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF3)
> summary(out3.glm)
Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF3)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4806 1.1969 0.40 0.6880
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race黒人 1.2723 0.5274 2.41 0.0158 *
raceその他 0.8805 0.4408 2.00 0.0458 *
この結果が、お示しいただいた out.glm2 の結果と一致するのを見て、やっと気付いたのですが、out1.glm に含まれるダミー変数 race2 , race3 の "2"や"3"は、2:黒人、3:その他 のラベルなんですね。
（No.22543を投稿した時点では、通し番号のようなものかな、と誤解していました）
ダミー変数 { race黒人, raceその他 } は、おそらく
race=1 について { 0, 0 }
race=2 について { 1, 0 }
race=3 について { 0, 1 }
という対応になっていて、
> exp(out3.glm$coefficients)
(Intercept) age lwt smoke ptl
1.62 0.97 0.97 2.56 1.72
ht ui ftv race黒人 raceその他
6.44 2.15 1.07 3.57 2.41
から、
・人種が「黒人」であることは、白人である場合と較べ、低体重児の出生リスクが 3.57倍になる
・人種が「その他」 〃 2.41倍になる
という解釈に行き着くのか、と。
|
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：tosh 2018/05/14(Mon) 19:48 No. 22550 |
投稿No. 22537に挙げた（３）の方法
（ダミー変数を -1, 0, 1 などとおく方法）
についてですが、掲載されているのは、ある大学病院で開催された「ＥＢＭ・統計セミナー」という研修会のPDF資料でした（現在はDLできなくなっていました）。
この資料には、次の２つの変換方法が紹介されており、
ア）基準とするカテゴリ（reference group）を仮定する場合
イ）カテゴリ全体の平均をreferenceとする場合
-------------------
ア）基準とするカテゴリ（reference group）を仮定する場合
例 名義尺度データが｛ "薬剤なし", "薬剤Ｐ", "薬剤Ｑ" } のような場合
ダミー変数 D1, D2 を
"薬剤なし" = {0,0}
"薬剤Ｐ" = {1,0}
"薬剤Ｑ" = {0,1}
とおく。ロジスティック回帰の結果、D1, D2 について、
b1 : D1 に対応する回帰係数
b2 : D2 に対応する回帰係数
が推定され、
exp(b1) は"薬剤なし"に対する"薬剤Ｐ"のオッズ比、
exp(b2) は"薬剤なし"に対する"薬剤Ｑ"のオッズ比、
という解釈となる。
--------------------
イ）カテゴリ全体の平均をreferenceとする場合
例 名義尺度データが｛ "地域Ａ", "地域Ｂ", "地域Ｃ" ｝のような場合
ダミー変数 D1, D2 を
"地域Ａ" = {1,0}
"地域Ｂ" = {0,1}
"地域Ｃ" = {-1,-1}
とおく。
b1: "地域Ａ"の回帰係数
b2: "地域Ｂ"の回帰係数
b3: "地域Ｃ"の回帰係数
とし、モデル上で、
b1 + b2 + b3 = 0
という制約を置いたうえで推定する。
計算上は、b1, b2 しか推定されないが、b3 = - b1 - b2 として計算可能である。
exp(b1)をとればオッズ比が得られる。
--------------------
引用ここまでです。
最後の一文にある、「exp(b1)をとれば」は、「exp(-b1-b2)をとれば」ではないのかな…？ という疑問がありますが、原文ママです。
なお、この資料の末尾の参考文献には次が挙げられていました。
青木繁伸（2009）「Ｒによる統計解析」
大橋靖雄（2011）「医師のための臨床統計学基礎編」
丹後俊郎（1996）「路地すてぃく回帰分析SASを利用した統計解析の実際」
豊田秀樹（1992）「原因をさぐる統計学」
服部環・海保博之（1996）「Ｑ＆Ａ心理データ解析」
森實敏夫（2004）「入門 医療統計学」 |
| Re: ロジスティック回帰に用いる説明変数が、３種類（以上）の値を含む名義尺度データ... |
| 投稿者：tosh 2018/05/14(Mon) 20:12 No. 22551 |
No.22545でお示しいただいた、自分でダミー変数を作る方法を真似して、（３）の方法でダミー変数を作成する処理を追加してみました。
###################################################
#
# ロジスティック回帰分析例 2
# https://www.gixo.jp/blog/3200/
# より一部改変
library(MASS)
str(birthwt)
#
# low: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
# age: 母親の年齢
# lwt: 母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke: 母親の喫煙有無(喫煙有りが1)
# ptl: 母親の早産経験の有無(経験有りが1)
# ht: 母親の高血圧の有無(有りが1)
# ui: 母親の子宮神経過敏の有(有りが1)
# ftv: 妊娠第１期における医師の訪問回数
# bwt: 新生児の体重(グラム単位)
options(digits = 4)
### DF1 ダミー変数を、Rで自動生成する
# bwtは使わないので除く
DF1 <- birthwt[,1:9]
# lwtをポンド単位からg単位に変換する
DF1$lwt <- DF1$lwt * 0.454
# raceをfactor型に変換する
DF1$race <- factor(DF1$race,levels = c(1,2,3),labels = c("白人","黒人","その他"))
str(DF1)
# ロジスティック回帰
out1.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF1)
summary(out1.glm)
exp(out1.glm$coefficients)
#### DF4 ダミー変数を、白人{1,0} 黒人{0,1} その他{-1,-1} とする方法で自作する
DF4 <- DF1
# ダミー変数に変換する
race = birthwt$race
DF4$race白人 = ifelse(race==1, 1, ifelse(race==2, 0, -1))
DF4$race黒人 = ifelse(race==1, 0, ifelse(race==2, 1, -1))
DF4$race = NULL
str(DF4)
# ロジスティック回帰
out4.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF4)
summary(out4.glm)
exp(out4.glm$coefficients)
# "raceその他" の偏回帰係数、オッズ比
b3 = - coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"]
names(b3) = c("raceその他_coef")
b3_exp = exp(- coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"])
names(b3_exp) = c("raceその他_exp(coef)")
-----------------
スクリプトここまでです。
以下、結果を貼ります。
-----------------
> summary(out4.glm)
Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF4)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1982 1.1509 1.04 0.2978
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race白人 -0.7176 0.2699 -2.66 0.0079 **
race黒人 0.5547 0.3232 1.72 0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.3
Number of Fisher Scoring iterations: 4
> exp(out4.glm$coefficients)
(Intercept) age lwt smoke ptl ht ui ftv race白人
3.3142 0.9709 0.9666 2.5570 1.7217 6.4450 2.1547 1.0675 0.4879
race黒人
1.7414
> b3
raceその他_coef
0.1629
> b3_exp
raceその他_exp(coef)
1.177
------------
結果は以上です。
資料の説明からすると、
「低体重児の出生確率は、全人種の平均出生確率を基準にすると、
・白人の場合、0.4879倍
・黒人の場合、1.7414倍
・その他の人種の場合、1.177倍」
というような解釈になるのかな、と思いましたが…。
これ以上の詳細は、資料を作った方に直接質問するのが筋かと思いますが、具体的にやってみたことの報告のつもりで投稿しました。 |
|
