Rによるカテゴリカルデータ分析事例(3) ~{rms}パッケージによる順序ロジスティック回帰~
●2016年1月3日 追記
本エントリー中で使っている makedummies() がさらに便利に使えるようになっています。
アップデート版の使い方はコメント欄を参照してください。
以上、追記終わり。
----
前回までに、あるサービスを受けた500人の顧客の満足度および変数xは、性別や年齢層というデモグラフィック的な要素によって違いを受けるとは言えないことがわかりました。唯一、変数xは性別によって違いがあることがわかりました。
前回のエントリーはこちら↓
hikaru1122.hatenadiary.jp
変数xというとわかりにくいので、顧客のサービス参加度(以下、参加度)としておきます。本当はもっと複雑で、Kobe.R では分析そっちのけでこの議論が少しあったのですが、ここでは簡潔に「参加度」としておきます。要は、あるサービスに関わろうとする度合いが強い顧客ほど満足度が高いのではないか、ということを調べます。
知りたいこと
顧客のサービス参加度は満足度(cs)に正の影響を与えているかどうか。
→参加度が高いほど、満足度も高くなるのでは?
分析の手順
1.まず相関をチェック。
2.順序ロジスティック回帰分析で、参加度の満足度への影響を調べる。
使用するデータ
使うデータは前回のエントリーと同じです。ダウンロードは↓から。
https://app.box.com/s/0jikm4cfeu4iopdhthfkyqlrw4z7ibed
実証分析
0.データの確認 まずクロス表で参加度と満足度の状況を確認します。たしかに、参加度が大きいほど、満足度も高い人が多そうです。
> xtabs(~x + cs, data =d)
満足度1 | 2 | 3 | 4 | 5 | 6 | 7 | |
---|---|---|---|---|---|---|---|
参加度1 | 2 | 2 | 4 | 51 | 26 | 12 | 2 |
2 | 1 | 1 | 7 | 14 | 38 | 15 | 10 |
3 | 0 | 0 | 4 | 6 | 61 | 34 | 9 |
4 | 0 | 1 | 2 | 6 | 28 | 63 | 27 |
5 | 2 | 0 | 2 | 2 | 3 | 15 | 50 |
1.相関係数をチェック
回帰分析の前に、相関係数を調べておきたいと思います。特に相関がなければ、回帰分析をする必要がないと考えるからです。参加度は1~5の順序尺度、満足度は1~7の順序尺度です。連続値ではないので、厳密には普通の相関係数(ピアソンの積率相関係数)で調べることはできません。
5つや7つの順序尺度なら、連続値とみなして別にいいんじゃないかという意見もあります*1。私もそうしたいです、楽だし。でも、今回は学びのために、少し厳密になってみます。
3つ以上カテゴリを持つ順序カテゴリカル変数には、ポリコリック相関係数を調べます。cs とxを順序カテゴリカル変数に変換しておき、データの構成を確認します。
> d$cs2 <- ordered(d$cs) > d$x2 <- ordered(d$x) > str(d) 'data.frame': 500 obs. of 6 variables: $ gender : Factor w/ 2 levels "0","1": 1 1 1 1 1 2 2 1 1 2 ... $ age_seg: Factor w/ 5 levels "1","2","3","4",..: 3 4 4 4 3 1 2 4 3 3 ... $ cs : int 5 7 4 6 6 6 5 6 7 5 ... $ x : int 3 2 1 1 3 4 2 5 2 3 ... $ cs2 : Ord.factor w/ 7 levels "1"<"2"<"3"<"4"<..: 5 7 4 6 6 6 5 6 7 5 ... $ x2 : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 3 2 1 1 3 4 2 5 2 3 ...
ポリコリック相関係数を出すには、{polycor}パッケージのpolychor()
を使います。
> library(polycor) > polychor(d$cs2, d$x2) [1] 0.5913401
まずまずの値でしょうか。hetcor()
だとp値を知ることができます。出てくる数値は同じく、0.59です。今回は有意(p値 1.665e-12)のようです。
> hetcor(d$cs2, d$x2) Two-Step Estimates Correlations/Type of Correlation: d$cs2 d.x2 d$cs2 1 Polychoric d.x2 0.5913 1 Standard Errors: d$cs2 d.x2 0.02967 Levels: 0.02967 n = 500 P-values for Tests of Bivariate Normality: d$cs2 d.x2 1.665e-12
2.順序ロジスティック回帰
正の相関があることが確認できました。顧客の参加度(1~5)がサービスの満足度(1~7)に影響を与えているかどうかを調べます。どちらも間隔尺度なので、本来は普通の回帰分析を使うことはできません*2。応答変数が順序カテゴリカル変数の場合、順序ロジスティック回帰が適切です。説明変数もカテゴリカルデータです。ダミー変数を作って、分析をしてみることにします。
まずデータの準備です。少し回りくどいですが、初めにcs2 ではなく、cs を抜き出したのは、makedummies()
が順序カテゴリカル変数をすべてダミー変数化してしまうためです*3。顧客の参加度のみダミー変数を作りたいため、順序尺度ではない cs を抜き出しました。
> library(dplyr) ## データdから cs と x2 を抜き出します。 > d2 <- select(d, cs, x2) ## 顧客の参加度x2 のダミー変数を作ります。 > library(makedummies) > d2 <- makedummies(d2) ## 最後に 満足度 cs を順序カテゴリカルにします。 > d2$cs <- ordered(d2$cs) > str(d2) 'data.frame': 500 obs. of 5 variables: $ cs : Ord.factor w/ 7 levels "1"<"2"<"3"<"4"<..: 5 7 4 6 6 6 5 6 7 5 ... $ x2_2: num 0 1 0 0 0 0 1 0 1 0 ... $ x2_3: num 1 0 0 0 1 0 0 0 0 1 ... $ x2_4: num 0 0 0 0 0 1 0 0 0 0 ... $ x2_5: num 0 0 0 0 0 0 0 1 0 0 ...
ダミー変数のベースとなるのは「顧客の参加度 1」です。では、順序ロジスティック回帰分析を行います。使用するパッケージは{rms}で、関数はorm()
です*4。
> orm(cs ~ ., data = d2) ## ↑の ~. は cs 以外の全部の変数を使う指示の省略形です。 Logistic (Proportional Odds) Ordinal Regression Model orm(formula = cs ~ ., data = d2) Frequencies of Responses 1 2 3 4 5 6 7 5 4 19 79 156 139 98 Model Likelihood Discrimination Rank Discrim. Ratio Test Indexes Indexes Obs 500 LR chi2 213.58 R2 0.364 rho 0.569 Unique Y 7 d.f. 4 g 1.555 Median Y 5 Pr(> chi2) <0.0001 gr 4.736 max |deriv| 1e-05 Score chi2 215.96 |Pr(Y>=median)-0.5| 0.302 Pr(> chi2) <0.0001 Coef S.E. Wald Z Pr(>|Z|) y>=2 3.4341 0.4644 7.39 <0.0001 y>=3 2.8328 0.3564 7.95 <0.0001 y>=4 1.6431 0.2283 7.20 <0.0001 y>=5 -0.1027 0.1860 -0.55 0.5809 y>=6 -1.9960 0.2191 -9.11 <0.0001 y>=7 -3.8202 0.2575 -14.83 <0.0001 x2_2 1.1229 0.2786 4.03 <0.0001 x2_3 1.6875 0.2598 6.50 <0.0001 x2_4 2.7115 0.2729 9.94 <0.0001 x2_5 4.4556 0.3537 12.60 <0.0001
係数 x2_2 ~ x2_5 を確認すると、顧客の参加度1をベースにして、2~5と参加度が上がるごとに、係数が増えています。よって、顧客の参加度は満足度に正の影響を与えていると言えるので、顧客の参加を促す施策を行ったほうがよいことがわかります。
今回、参考にしたもの
- 作者: 藤井良宜,金明哲
- 出版社/メーカー: 共立出版
- 発売日: 2010/04/22
- メディア: 単行本
- クリック: 13回
- この商品を含むブログ (3件) を見る
- 作者: 太郎丸博
- 出版社/メーカー: ナカニシヤ出版
- 発売日: 2005/07
- メディア: 単行本
- 購入: 1人 クリック: 26回
- この商品を含むブログ (12件) を見る
- 作者: 山本勲
- 出版社/メーカー: 中央経済社
- 発売日: 2015/10/21
- メディア: 単行本
- この商品を含むブログ (1件) を見る
相関係数の補足
なお、ノンパラメトリックなスピアマンの順位相関係数を知りたい場合は、次のようにします。ポリコリック相関係数のほうが高い値が出やすいように思います。今回も少し高かったです(今回のポリコリック相関係数は0.59)。
> cor(d$cs, d$x, method = "spearman") [1] 0.5690423
順序ロジスティック回帰の補足
順序ロジスティック回帰は{MASS}パッケージのpolr()
がよく例として出されます。以下、polr()
でやってみます。なお、polr()
の結果はp値が出力されません*5。そのため、{texreg}パッケージのscreenreg()
でp値を確認します。
> library(MASS) > kekka <- polr(cs ~ ., data = d2) ## ↑の ~. は cs 以外の全部の変数を使う指示の省略形です。 > library(texreg) > screenreg(kekka) ## カッコ内は標準誤差です。 Re-fitting to get Hessian =========================== Model 1 --------------------------- x2_2 1.12 *** (0.28) x2_3 1.69 *** (0.26) x2_4 2.71 *** (0.27) x2_5 4.46 *** (0.35) --------------------------- AIC 1345.59 BIC 1387.73 Log Likelihood -662.79 Deviance 1325.59 Num. obs. 500 =========================== *** p < 0.001, ** p < 0.01, * p < 0.05
結果は同じです。こっちのほうがシンプルな結果が出るので、やはりpolr()
が便利なのでしょうか。