読者です 読者をやめる 読者になる 読者になる

疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

R package'nnet'と'MASS'を用いて多項ロジスティック回帰分析を行う

今回はRのパッケージ'nnet'を用いて、多項ロジスティック回帰分析を実施する。

はじめに

よく医学の論文で目にするのが「ロジスティック回帰分析」である。このように何も定義していない場合は、結果変数が二値(Yes or No)である二項ロジスティック回帰分析であることが多い。非常に頻繁に活用し、オッズ比やROC曲線を引く場合にも必要な分析である。
今回はその拡張である多項ロジスティック回帰分析(結果変数が3つ以上の場合)を実施する。

カテゴリーに順番のない場合(名義変数)からスタート!

考え方としては3水準ある場合には、ひとつのレベルを基準にして、その他2つのレベルとの確率の比を求める。
Rではおなじみのirisデータを用いる。

nnetをインストールし、読み込む

install.packages("nnet")
library(nnet)

早速、分析する。

変数Speciesにはsentosa, versicolor, virginicaの3種類(順番に意味はない)が含まれている。ここではSpal.lengthとspeciesの関連を調査する。
この場合familyは必然的に多項分布を指定するため入力の必要はない。(通常のロジスティック回帰分析はbinomialでリンク関数はlogit)

multinom(Species ~ Sepal.Length, data =iris)

出力がこちら

# weights:  9 (4 variable)
# initial  value 164.791843 
# iter  10 value 91.337114
# iter  20 value 91.035008
# final  value 91.033971 
# converged
# Call:
# multinom(formula = Species ~ Sepal.Length, data = iris)
#
# Coefficients:
#                     (Intercept) Sepal.Length
# versicolor   -26.08339     4.816072
# virginica    -38.76786     6.847957
# 
# 
# Residual Deviance: 182.0679 
# AIC: 190.0679 

この結果でcoefficientsの部分に注目すると、どちらもSepal.Lengthの部分が正の値になっていることからSepal.Lengthが長くなれば、sentosaよりもその他2種になる可能性が高いことを示唆している。

しかし、この結果ではversicolorとvirginicaの比較は難しいが、上記のモデルより

{ \displaystyle
log\frac{P(virginica)}{P(versicolor)} = (a_2-a_1) + (b_2-b_1) Sepal.Length
}

と記述できるので、この場合{\beta_(virginica)}=6.847957が{\beta_(versicolor)}=4.816072よりも大きいことからSepal.Lengthが長くなるにつれて、virginicaが多い確率の方が高いことがわかった。

次にカテゴリー間に順序のある場合(順序変数)

今回はvcdというパッケージ内のhousingというデータセットを用いる。
(20160218現在、このデータセットはなかったのでhousingを模したデータを作成し、用いた)

summary(housing)
#     sat       cont         Freq       
# High  :24   High:36   Min.   : 6.127  
# Low   :24   Low :36   1st Qu.:26.622  
# Medium:24             Median :51.068  
#                                Mean   :48.389  
#                               3rd Qu.:70.541  
#                                Max.   :85.991  

家主の満足度を示す変数Satと住民間の強さを示す変数Contを用いる。
満足度を示すSatはlow, medium, highの3つの順番のあるグループに分けられる。
ここでは満足度を二つのグループに分けることを考える。
つまり、
{ \displaystyle
log\frac{P(low)}{P(high) + P(medium)} = \alpha_L + \beta Cont high
}

{ \displaystyle
log\frac{P(low)+ P(medium)}{P(high)} = \alpha_M + \beta Cont high
}

このようにどちらのモデルにもCont highが含まれているものを比例オッズモデルと呼ぶ。
この場合にはMASSパッケージに含まれるpolr関数を用いる。

polr(sat ~ cont, data = housing, weight = Freq)
# Call:
# polr(formula = sat ~ cont, data = housing, weights = Freq)
# 
# Coefficients:
# contLow 
# -0.70035 
# 
# Intercepts:
#   High|Low Low|Medium 
# -0.9746715  0.4876438 
# 
# Residual Deviance: 7530.293 
# AIC: 7536.293 

結果の解釈としては、 {\alpha_L}=-0.975、 {\alpha_M}=0.488、 {\beta}=-0.700となっており、この場合Cont highが負になっていることからContがLowよりもHighの時に満足度は低下していることになる。

参考文献

藤井良宜著、金明晢編『Rで学ぶデータサイエンス1 カテゴリカルデータ解析]』共立出版株式会社(2010)


20160218
RF