Knowledge As Practice

JAIST(東京)で Transformative Service Research に取り組んでる社会人大学院生の研究・勉強メモ

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と参加度が上がるごとに、係数が増えています。よって、顧客の参加度は満足度に正の影響を与えていると言えるので、顧客の参加を促す施策を行ったほうがよいことがわかります。

 

今回、参考にしたもの

カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)

カテゴリカルデータ解析 (Rで学ぶデータサイエンス 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()が便利なのでしょうか。

*1:書籍によっては「連続値をみなして分析OK」と書いてあるものがあります。

*2:そんなに結果は変わらないし、私はいいかなと思うのですが…。

*3:makedummies についてはこちらを参照してください。

hikaru1122.hatenadiary.jp

*4:lrm関数を用いても、同じことができます。

*5:とはいえ、summaryでt値が出るので、だいたい有意かどうかはわかります。

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 - 継承 4.0 国際 ライセンスの下に提供されています。