Knowledge As Practice

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

論文「顧客収益性の統計的分析」を読む

最近、統計分析の学習をしていなかったので、勉強になる論文を読んでみました。マルチレベル分析の経営学での応用事例を検索していて見つけたものです。偶然というか、第1著者は私の専門職学位論文の分析パートをチェックしてくださった先生、第3著者は専門職大学院のゼミの先生です*1

 
論文タイトルは「顧客収益性のと統計的分析 ─管理会計研究のマルチレベル分析の適用可能性─」*2です。内容は、階層線形モデル(HLM、マルチレベル分析)を管理会計研究・実務へ活かす可能性を探ったものです。戦略論、管理会計研究の既存研究を概観し、シミュレーションを行っています。

 
特に勉強になったのは、

  • ある階層のデータが少なくとも5を下回ると結果の誤差は無視できなくなるほど大きくなる。
  • ある階層のデータが少なくとも10あれば、HLM は頑健な結果を得られる。
  • 分散の極端な偏りに対しても、頑健な結果が得られる。

というものです。

 
HLM を使う利点として「多階層におよぶ影響がひとつの利益指標に集約されている時にHLMを用いると、それらがどのような階層からの影響力が大きいのかを説明できる」と述べています。例として、先行研究レビューと顧客別利益を支店・担当者・顧客・時間という4つの階層でシミュレーションを行っています。

 
シミュレーションで使用された R のパッケージはnlmeです。glmmMLlmerTest しか使ったことがないので、このパッケージを知ったことも勉強になりました。マルチレベル分析については、心理学や生態学の例を書籍やネットで見ることが多かったので、経営学のフィールドにいる自分としては、こういう論文はとても身近に感じます。

 
なお、マルチレベル分析の組織マネジメントへの応用としては、次の本があります。マーケティング・サイエンスの世界では、高度な分析がたくさん行われているようですが*3、その他の経営学の分野ではまだまだ少ないように思います。

*1:私の専門フィールドはマーケティング論(の特にサービス研究)ですが、ちょっと大人の事情で管理会計の先生のゼミに所属していたのです。

*2:新井康平, 大浦啓輔, 加登豊. (2014). 顧客収益性の統計的分析: 管理会計研究へのマルチレベル分析の適用可能性. 原価計算研究, 38(2), 78-88.

*3:学会に行っても、ほとんど理解できずに帰ってきてしまいます…。

サービス・マーケティング3つのパラダイム

イタリアのマーケティング学会誌(たぶん)の Editorial*1 が勉強になったので、自分用メモ。著者は Evert Gummesson。日本では、ガメソンとかグメソンとかグマソンとか呼ばれています。たぶん、いちばん近い発音はグマソンかと。ストックホルムビジネススクールの先生です。

 
PDFファイルはこちら↓
http://www.francoangeli.it/Area_RivistePDF/getArticolo.ashx?idArticolo=45509

 
自分の研究に関係がありそうなものをネット上で物色してたどり着きました。イタリア語はよくわからないのですが、どうもサービス・マーケティングの特集号のようです。この Editorial はサービス・マーケティング研究の歴史を3つにわけて簡単に紹介しています*2

 
その3つは次のとおりです。

  • パラダイム1 1970年代まで サービスが認知される、すべてはモノ・製造業の時代。
  • パラダイム2 1970年代〜2000年まで モノとサービスの違いを追求した時代
  • パラダイム3 2000年以降 共有・相互依存およびシステムアプローチの時代

 
現代のサービス研究に重要なのは3つ目のパラダイムです。このパラダイムには、さらに3つの柱があると著者は言っています。

  1. S-Dロジック
  2. サービス・サイエンス
  3. 複雑性の研究(Complexity Studies)

 
S-Dロジックはすでに当ブログで詳しく説明しています。サービス・サイエンスはIBMのお商売です*3。私には3つ目の「複雑性の研究」が勉強になりました。マーケティングの複雑性は、ケーススタディ・ネットワーク理論・システム理論で研究されている、と。システム理論による「VSAアプローチ」という方法があることも初めて知りました。いまはそれぞれ何のことだかよくわかりませんが、参考文献がちゃんと紹介されているので、後で調べてみることにしますå。

 
自分の研究フィールドは零細業ですから、すぐに役立てられるかわからないですが、ちょうど研究の幅を広げたいなと思っていたところに読んだ論考なので、勉強になりました。

*1:Gummesson, E. (2012). Editoriale. The three service marketing paradigms: which one are you guided by?. Mercati e Competitività.

*2:サービス・マーケの歴史はFisk et al.の有名な論文がありますが、それより簡潔で読みやすいです。

*3:こんなこと言うと怒られるかも…。

今年のデータ分析の学習はじめはドキュメント作成からスタート

今年は Reproducible research を意識したいです(下のスライドが参考になります)。昨年はデータと分析結果がとっ散らかって、自分が分析した手順・内容を忘れしまい、何度も同じことを繰り返すことがありましたので…。

www.slideshare.net

 
そこで、いつもは図書館で借りて読んでいる*1『ドキュメント・プレゼンテーション生成』を買ってきて、この週末に読みました。

ドキュメント・プレゼンテーション生成 (シリーズ Useful R 9)

ドキュメント・プレゼンテーション生成 (シリーズ Useful R 9)

 
恒例というか私の環境が悪いせいで、本どおりにやってもうまくいかず、PDF 出力で苦労しました。↓がうんともすんとも動かない。

knit2pdf("minimal.Rnw", compiler = "lualatex", encoding="UTF8")

 
しかたがないので、Rstudio で PDF 出力をしていけばいいかとあきらめました。Rstudio だとちゃんと PDF 出力できました。他にもbrowseURL("sample.html")とすると、勝手に Rstudio が立ち上がるので困ったり、ITエンジニア・プログラマ属性がある方には簡単に解決できそうなこともヒーヒー言いながら、読了。

 
なお、フルパスでブラウザのありかを指定すればOKです(以下は私の場合)。

browseURL("sample.html", browser="C:\\Program Files (x86)\\Mozilla Firefox\\firefox.exe")

 
その後「R Markdownで楽々レポートづくり」を読んで復習しました。ということで、今年の内輪的な勉強会は R Markdown でぜんぶやろうと思います。

gihyo.jp

*1:もし、大阪市中央図書館に行って、お目当てのR関連の本が見つからなければ、だいたい私が借りています。

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値が出るので、だいたい有意かどうかはわかります。

Rによるカテゴリカルデータ分析事例(2) ~(1) の補足:等分散性を確認していたか?~

前回のエントリーを書いたあと、現在は手に入りにくい本『Rによる統計解析の基礎』(私は買いました)の最新版『Rによる保健医療データ解析演習』を確認し、さらに「マイナーだけど最強の統計的検定 Brunner-Munzel 検定」を読むと、ウィルコクソンの順位和検定(マン・ホイットニーのU検定)*1クラスカル・ウォリス検定には等分散性の確認が必要のようです。

Rによる統計解析の基礎 (Computer in Education and Research)

Rによる統計解析の基礎 (Computer in Education and Research)

 
『Rによる保健医療データ解析演習』は↓の「ダウンロード」セクションでゲットできます。とても勉強になります。
http://minato.sip21c.org/msb/index.html

 
d.hatena.ne.jp

念のため、分散の同質性を確認します。この確認にはフリグナー・キリーン検定を行います。デフォルトで使えるfligner.testか {coin}パッケージの fligner_test を使います。せっかくなので、前回と同様、{coin} パッケージを使ってみます。データは前回のエントリーと同じものです。帰無仮説は「各群(性別または年齢層)の分散は同じ」。有意水準を5%とします。

 
次の分析結果から、帰無仮説を棄却することはできませんでした。よって、前回のエントリーの分析を行ってOKです。

> fligner_test(cs ~ gender, data =d)

    Asymptotic Two-Sample Fligner-Killeen Test

data:  cs by gender (0, 1)
Z = -1.0815, p-value = 0.2795
alternative hypothesis: true ratio of scales is not equal to 1

> fligner_test(cs ~ age_seg, data =d)

    Asymptotic K-Sample Fligner-Killeen Test

data:  cs by age_seg (1, 2, 3, 4, 5)
chi-squared = 5.6545, df = 4, p-value = 0.2265

> fligner_test(x ~ gender, data =d)

    Asymptotic Two-Sample Fligner-Killeen Test

data:  x by gender (0, 1)
Z = 0.00056251, p-value = 0.9996
alternative hypothesis: true ratio of scales is not equal to 1

> fligner_test(x ~ age_seg, data =d)

    Asymptotic K-Sample Fligner-Killeen Test

data:  x by age_seg (1, 2, 3, 4, 5)
chi-squared = 6.177, df = 4, p-value = 0.1863

*1:等分散性が確認できれなければ、 Brunner-Munzel 検定を使います。

Rによるカテゴリカルデータ分析事例(1)  {coin}パッケージによるウィルコクソンの順位和検定とクラスカル・ウォリス検定

先日(2015年12月26日土曜日)に開催された Kobe.R 第23回で発表を行いました。当日はちょっと風邪気味だったため、あまり準備ができず、テキトーな内容になってしまいました…。そのため、ここで補足をしたいと思います。

kobexr.doorkeeper.jp

 
発表タイトルは「Rによるカテゴリカルデータ分析事例」。現在、私が分析を進めているデータを使った内容です。扱ったデータの構成は

  • 性別(0が男性、1が女性の名義尺度)
  • 年齢層(1が20代、2が30代、3が40代、4が50代、5が60代の名義尺度)
  • あるサービスを受けた満足度(1~7の順序尺度)
  • 変数x(1~5の順序尺度)

というもの。実際のデータセットはもっと多くて50個ほど変数があります。サンプルサイズは500です。今回使用するデータは↓においてあるので、興味がある方はやってみてください。 https://app.box.com/s/0jikm4cfeu4iopdhthfkyqlrw4z7ibed

 
分析によって調べたいことは次の3つです。

  1. 満足度および変数xは性別によって異なるか。
  2. 満足度および変数xは年齢層によって異なるか。
  3. 変数xは満足度に正の影響を与えているか(変数xが高いほど、満足度も高いか)。

 
今回は1と2を扱います。3は次回のエントリーに回します。では、やってみます。有意水準は5%とします。

 

0.使用するデータの準備と確認

データの構成は次のようになっています。

> d <- blogyou.data.20151229
> str(d)
'data.frame':   500 obs. of  4 variables:
 $ gender : int  0 0 0 0 0 1 1 0 0 1 ...
 $ age_seg: int  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 ...

 
性別と年齢層を因子型に変えて、もう1度構成を確認します。

d$gender <- as.factor(d$gender)
d$age_seg <- as.factor(d$age_seg)
> str(d)
'data.frame':   500 obs. of  4 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 ...

 

1.満足度および変数xは性別によって異なるか?

性別は2値の名義尺度、満足度と変数xは順序尺度です。まずはクロス表を作ってデータを眺めてみます。

 
●満足度と性別の関係

> xtabs(~gender + cs, data=d)
1 2 3 4 5 6 7
男性(0) 1 3 6 38 66 59 35
女性(1) 4 1 13 41 90 80 63

 
ん~、どうでしょう? 女性のほうが高い満足度にある人が多いような…。図示してみます。

> library(lattice)
> histogram(~factor(cs)|gender, data = d)
## cs は離散値のため、factor() で因子型に直しています。以下、x も同様です。

f:id:hikaru1122:20151229181954p:plain

では、実際に性別によって満足度の差があるか検定します。上図より正規分布に従っているようには見えないので、ノンパラメトリック検定のウィルコクソンの順位和検定(マン・ホイットニーのU検定)を行います。帰無仮説は「性別によって差がない(差が0)」、対立仮説が「性別によって差がある(差が0でない)」です。

 
さて、ここで{coin}パッケージを使います。ウィルコクソンの順位和検定を行うには、wilcox.test や {exactRankTests} の wilcox.exact を使うのですが、{exactRankTests}をインストールして使おうとしたら「{coin}を使え!」とコンソール画面で怒られたので、{coin} を使うことにします。

## {coin}パッケージがない場合は、インストールしてください。
library(coin)
> kekka <- wilcox_test(cs ~ gender, distribution = "exact", data = d)
## distribution = "exact"とすると、パソコンによっては結果が出るまで時間がかかります。
> kekka
    Exact Wilcoxon-Mann-Whitney Test
data:  cs by gender (0, 1)
Z = -1.0953, p-value = 0.2737
alternative hypothesis: true mu is not equal to 0

 
帰無仮説は棄却できませんでした。したがって、性別によって満足度に違いがあるとは言えませんでした。次に、変数xと性別の関係と調べます。

 
●変数xと性別の関係
先ほどと同じ検定を行いますが、その前にクロス表とヒストグラムで確認します。

> xtabs(~gender + x, data =d)
1 2 3 4 5
男性(0) 52 30 51 53 22
女性(1) 47 56 63 74 52
histogram(~factor(x)|gender, data =d)

f:id:hikaru1122:20151229184129p:plain

女性は特に変数xが5の人が多いように見えます。検定します。

> kekka2 <- wilcox_test(x ~ gender, distribution = "exact", data = d)
> kekka2

    Exact Wilcoxon-Mann-Whitney Test

data:  x by gender (0, 1)
Z = -2.2005, p-value = 0.02769
alternative hypothesis: true mu is not equal to 0

変数x は性別によって違いがあることがわかりました。

 

2.満足度および変数xは年齢層によって異なるか。

年齢層は5つの水準がある名義変数です。年齢層と満足度および変数xの関係をクロス表とヒストグラムで確認します。

> xtabs(~age_seg + cs, data =d)
1 2 3 4 5 6 7
20代(1) 0 1 4 18 32 25 20
30代(2) 4 2 3 18 34 22 17
40代(3) 0 1 6 19 30 25 19
50代(4) 0 0 0 14 31 37 18
60代(5) 1 0 6 10 29 30 24
> xtabs(~age_seg + x, data =d) 
1 2 3 4 5
20代(1) 24 23 17 19 17
30代(2) 20 22 24 19 15
40代(3) 25 10 24 29 12
50代(4) 18 16 21 29 16
60代(5) 12 15 28 31 14
histogram(~factor(cs)|age_seg, data =d)
histogram(~factor(x)|age_seg, data =d)

f:id:hikaru1122:20151229184314p:plain f:id:hikaru1122:20151229184322p:plain

どうも年齢層による違いは満足度、変数xともに見られそうにありません。検定します。水準が3つ以上なので、先ほどのウィルコクソンの順位和検定(マン・ホイットニーのU検定)ではなく、クラスカル・ウォリス検定を行います。

## kruskal_test() は {coin} パッケージにあります。
kekka3 <- kruskal_test(cs ~ age_seg, data = d)
kekka4 <- kruskal_test(x ~ age_seg, data = d)

 
満足度および変数xの年齢層による違いは次の結果から差があるとはいえませんでした。

> kekka3

    Approximative Kruskal-Wallis Test

data:  cs by age_seg (1, 2, 3, 4, 5)
chi-squared = 8.0373, p-value = 0.0895

> kekka4

    Asymptotic Kruskal-Wallis Test

data:  x by age_seg (1, 2, 3, 4, 5)
chi-squared = 5.4405, df = 4, p-value = 0.245

次回に続きます。今回、参考にしたものはこちらの2冊です。

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

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

人文・社会科学のためのカテゴリカル・データ解析入門

人文・社会科学のためのカテゴリカル・データ解析入門

Rでダミー変数をつくる。makedummies 関数で。

●2015年12月28日 追記

makedummies 関数の作者の方から、コメントで情報をいただきました。 以下のように、 install_github で 簡単に makedummies が使えるようになっています。

## devtools がなければインストール→ install.packages("devtools")
library(devtools)
install_github("toshi-ara/makedummies")

●2016年1月3日 追記

makedummies がさらに便利になりました。詳細は↓です。
http://hikaru1122.hatenadiary.jp/entry/2015/12/30/213032

以上、追記おわり。

 

今回のモチベーション

最近、ずっとカテゴリカルデータをいじっていて、ダミー変数を作る機会が出てきました。Rでダミー変数を作るのには、パッケージがあったはずだと思い、調べていたら、こちらの2つがみつかりました。

 
dummyVars についてはこちら。 estrellita.hatenablog.com

 
dummy.data.frame についてはこちら。 qiita.com

 

makedummies 関数

ただ、これらはダミー変数を作るときに基準となる変数もそのまま作ります。なので、基準となる変数のダミーはいりません。勝手に基準となるダミーを作らない関数がないかと探ししたところ、ありました。

 
それがこちらの makedummies 関数です。 github.com

 
作ってくださった方に感謝です。データフレームしか受け付けないので、その点だけ注意すればダミー変数が簡単につくれます。github のソースをコピーして、適当なファイルに保存して、適当なフォルダに入れておきます。今回は、makedummies.R というファイルにして、getwd() で表示されるフォルダに入れています。

source("makedummies.R")
> d <- read.delim("clipboard")
> str(d)
'data.frame':   18 obs. of  2 variables:
 $ y: int  5 7 4 6 6 6 6 6 7 5 ...
 $ x: int  1 2 3 1 2 3 1 2 3 1 ...
> head(d)
  y x
1 5 1
2 7 2
3 4 3
4 6 1
5 6 2
6 6 3

 
x のダミー変数を作りたいと思います。x が2のときに「1」となるダミー、x が3のときに「1」となるダミーを作ります。1は基準なのでダミーはいりません。ダミー変数としたいものは、因子型である必要があるので、因子型に直します。

> d$x <- factor(d$x)
> str(d)
'data.frame':   18 obs. of  2 variables:
 $ y: int  5 7 4 6 6 6 6 6 7 5 ...
 $ x: Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...

 
順序付きのカテゴリカル変数でもOKです。あとは makedummies 関数に入れるだけ。自動的に変数名がつきます。無事、基準となる変数はダミーが作られていません。ちなみに、複数の因子型変数がある場合、全部ダミー変数を作ります。

> d2 <- makedummies(d)
> head(d2)
  y x_2 x_3
1 5   0   0
2 7   1   0
3 4   0   1
4 6   0   0
5 6   1   0
6 6   0   1

私には重宝する関数です。改めて感謝。

 
ちなみに、カテゴリカルデータの分析はこちらの2冊が勉強なりました。

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

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

人文・社会科学のためのカテゴリカル・データ解析入門

人文・社会科学のためのカテゴリカル・データ解析入門

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