量的アプローチによる価値共創の研究発表＠第6回知識共創フォーラム in 金沢
先週（2016年3月12〜13日），金沢で開催された第6回知識共創フォーラムで発表をしてきました。相馬・清水（2015）*1に大きく影響を受けた分析の結果発表でした。もともと潜在ランク理論はどこかで知っていて，自分の研究に使えそうだなぁと思っていたところに相馬・清水（2015）が出たので，その時はとても驚きました。
発表タイトルは「価値共創の促進に意味はあるか？」です。情報学系な人には参考にならないと思いますが，経営学系で量的な研究をされる方（本当に少ないと思うけれど…）の少しでも参考になれば，と思います。
発表に使ったスライドは下です。
www.slideshare.net
分析の方法・ツールなど
分析の方法・ツールなどは次のとおりです。過去，本ブログや Kobe.R で発表してきたまとめ的な分析を行いました。
| 分析の方法 | 使用したツール |
|---|---|
| 潜在ランク理論 | exametrica |
| ウィルコクソンの順位和検定 | R, wlicox.test() |
| クラスカル・ウォリス検定 | R, kruskal.test() |
| 線形連関の検定 | R, coin::lbl_test() |
| ベイズ推定による多項ロジスティック回帰 | R, brms::brm() |
| ロジスティック回帰（普通の） | R, glm() |
今回の発表はケーススタディが多い（ような気がする）共創研究に量的アプローチを行った例です。分析に使ったアンケートは1〜7のリッカート尺度なので，独立性の検定にはノンパラメトリック検定を使っています。
ベイズ推定による多項ロジスティック回帰
brms パッケージを使用しました。以下のコードで一発です。下でも書きましたが，brm()のオプションはまだ理解していません。
library(brms) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) fit1 <- brm(qol ~ kyousou, data = d, family = categorical)
使用したデータはこちらに置いておきました。気になる方はお試しください。
分析でわからないところ
誰にも聞けなかった（というか，周りにわかる人がいない）ため，自信がない・わからないところは次のとおりです。
最後の2つの分析（多項ロジスティックと普通のロジスティック回帰）では，アンケート結果を間隔尺度として扱いました。尺度の扱いに一貫性がないのはどうかなぁ，と思っていますが判断がつきません…。
さらに，多項ロジスティック回帰はアンケート結果を名義尺度にしてしまいました。これも適切かどうか…。1〜7のうちどの選択肢が選ばれるかという観点から多項ロジスティック回帰でやってみようとした次第です。
多項ロジスティック回帰はベイズ推定を使ってやってみましたが，オプションは初期値ままです。収束はしたので大丈夫かとは思いましたが，適切な初期設定を知りたいところです。
発表を終えて
実は，今回が初めての口頭発表でした。聞いてくれた方のうち，おもしろいと思ってくれた方はいたようでホッとしました。どのへんがおもしろかったか，ちゃんと確認しておけばよかった。また今回の発表では因果関係を示したわけではないので，因果を調べられるようになりたいものです。
*1:相馬敏彦・清水裕士(2015)「ワンランク上のブランド・コミットメントはどう形成されるのか？ ─顧客の潜在ランクへの分類と拡張版投資モデルのブランドの適用─」『マーケティングジャーナル』第35巻第3号，75-94頁．
その有意差に意味あるの？ ～ノンパラメトリック検定と効果量の出し方～
Kobe.R 第24回で発表をしてきました。今回は来週日曜日に金沢で発表する内容のデータ分析部分（特にRの操作）をメインに切り出したものを発表しました。
↓使用したスライドはこちらにあります↓
RPubs - Rによるノンパラメトリック検定と効果量の出し方
slideshare はこちらです。Rstuido だと Rpubs にも PDF にも簡単出力できるので便利ですね。
トピックとしては
- wilcox.test（マイ・ホイットニーのU検定）は formula を使ったほうがわかりやすいこと。
- 効果量 cliff の d は orddom パッケージの delta で確認できること。
- R の発表は Rmd が便利だよ。
の３つでした。
使用したデータはこちらからダウンロードできるようにしました。どうぞお試しください。
Dropbox - data2.txt
@kazutan先生にご参加いただき、Rmd を使って発表できたいい経験ができました。というか、急遽、Rmarkdown について発表をお願いいたしました。急なお願いにも対応してくださり、感謝です。
【追記 2016年3月6日】
今回のデータのように順位にタイがある場合のウィルコクソンの順位和検定は coin パッケージの wilcox_test() がよりよいようです。でも、結果はほとんど変わりません。
【普通のウィルコクソンの順位和検定】
wilcox.test(rank ~ gender, data = d, correct = F)
Wilcoxon rank sum test
data: rank by factor(gender)
W = 27532, p-value = 0.07193
alternative hypothesis: true location shift is not equal to 0
【coin パッケージのウィルコクソンの順位和検定*1】
wilcox_test(rank ~ factor(gender), data = d)
Asymptotic Wilcoxon-Mann-Whitney Test
data: rank by factor(gender) (0, 1)
Z = -1.7995, p-value = 0.07193
alternative hypothesis: true mu is not equal to 0
【より正確にするならこう】
wilcox_test(rank ~ factor(gender), data = d, distribution = "exact")
Exact Wilcoxon-Mann-Whitney Test
data: rank by factor(gender) (0, 1)
Z = -1.7995, p-value = 0.07194
alternative hypothesis: true mu is not equal to 0
*1:gender を因子型にしないと動かない。
論文「顧客収益性の統計的分析」を読む
最近、統計分析の学習をしていなかったので、勉強になる論文を読んでみました。マルチレベル分析の経営学での応用事例を検索していて見つけたものです。偶然というか、第1著者は私の専門職学位論文の分析パートをチェックしてくださった先生、第3著者は専門職大学院のゼミの先生です*1。
論文タイトルは「顧客収益性のと統計的分析 ─管理会計研究のマルチレベル分析の適用可能性─」*2です。内容は、階層線形モデル（HLM、マルチレベル分析）を管理会計研究・実務へ活かす可能性を探ったものです。戦略論、管理会計研究の既存研究を概観し、シミュレーションを行っています。
特に勉強になったのは、
- ある階層のデータが少なくとも5を下回ると結果の誤差は無視できなくなるほど大きくなる。
- ある階層のデータが少なくとも10あれば、HLM は頑健な結果を得られる。
- 分散の極端な偏りに対しても、頑健な結果が得られる。
というものです。
HLM を使う利点として「多階層におよぶ影響がひとつの利益指標に集約されている時にHLMを用いると、それらがどのような階層からの影響力が大きいのかを説明できる」と述べています。例として、先行研究レビューと顧客別利益を支店・担当者・顧客・時間という4つの階層でシミュレーションを行っています。
シミュレーションで使用された R のパッケージはnlmeです。glmmML と lmerTest しか使ったことがないので、このパッケージを知ったことも勉強になりました。マルチレベル分析については、心理学や生態学の例を書籍やネットで見ることが多かったので、経営学のフィールドにいる自分としては、こういう論文はとても身近に感じます。
なお、マルチレベル分析の組織マネジメントへの応用としては、次の本があります。マーケティング・サイエンスの世界では、高度な分析がたくさん行われているようですが*3、その他の経営学の分野ではまだまだ少ないように思います。
学習を促す組織文化 -- マルチレベル・アプローチによる実証分析
- 作者: 北居明
- 出版社/メーカー: 有斐閣
- 発売日: 2014/03/05
- メディア: 単行本
- この商品を含むブログ (1件) を見る
サービス・マーケティング3つのパラダイム
イタリアのマーケティング学会誌（たぶん）の Editorial*1 が勉強になったので、自分用メモ。著者は Evert Gummesson。日本では、ガメソンとかグメソンとかグマソンとか呼ばれています。たぶん、いちばん近い発音はグマソンかと。ストックホルム・ビジネススクールの先生です。
PDFファイルはこちら↓
http://www.francoangeli.it/Area_RivistePDF/getArticolo.ashx?idArticolo=45509
自分の研究に関係がありそうなものをネット上で物色してたどり着きました。イタリア語はよくわからないのですが、どうもサービス・マーケティングの特集号のようです。この Editorial はサービス・マーケティング研究の歴史を3つにわけて簡単に紹介しています*2。
その3つは次のとおりです。
- パラダイム１ 1970年代まで サービスが認知される、すべてはモノ・製造業の時代。
- パラダイム２ 1970年代〜2000年まで モノとサービスの違いを追求した時代
- パラダイム３ 2000年以降 共有・相互依存およびシステムアプローチの時代
現代のサービス研究に重要なのは3つ目のパラダイムです。このパラダイムには、さらに3つの柱があると著者は言っています。
- S-Dロジック
- サービス・サイエンス
- 複雑性の研究(Complexity Studies)
S-Dロジックはすでに当ブログで詳しく説明しています。サービス・サイエンスはIBMのお商売です*3。私には3つ目の「複雑性の研究」が勉強になりました。マーケティングの複雑性は、ケーススタディ・ネットワーク理論・システム理論で研究されている、と。システム理論による「VSAアプローチ」という方法があることも初めて知りました。いまはそれぞれ何のことだかよくわかりませんが、参考文献がちゃんと紹介されているので、後で調べてみることにしますå。
自分の研究フィールドは零細業ですから、すぐに役立てられるかわからないですが、ちょうど研究の幅を広げたいなと思っていたところに読んだ論考なので、勉強になりました。
今年のデータ分析の学習はじめはドキュメント作成からスタート
今年は Reproducible research を意識したいです（下のスライドが参考になります）。昨年はデータと分析結果がとっ散らかって、自分が分析した手順・内容を忘れしまい、何度も同じことを繰り返すことがありましたので…。
www.slideshare.net
そこで、いつもは図書館で借りて読んでいる*1『ドキュメント・プレゼンテーション生成』を買ってきて、この週末に読みました。
ドキュメント・プレゼンテーション生成 (シリーズ Useful R 9)
- 作者: 高橋康介,金明哲
- 出版社/メーカー: 共立出版
- 発売日: 2014/06/24
- メディア: 単行本
- この商品を含むブログ (2件) を見る
恒例というか私の環境が悪いせいで、本どおりにやってもうまくいかず、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 でぜんぶやろうと思います。
Rによるカテゴリカルデータ分析事例(3) ～{rms}パッケージによる順序ロジスティック回帰～
●2016年1月3日 追記
本エントリー中で使っている makedummies() がさらに便利に使えるようになっています。
アップデート版の使い方はコメント欄を参照してください。
以上、追記終わり。
----
前回までに、あるサービスを受けた500人の顧客の満足度および変数ｘは、性別や年齢層というデモグラフィック的な要素によって違いを受けるとは言えないことがわかりました。唯一、変数ｘは性別によって違いがあることがわかりました。
前回のエントリーはこちら↓
hikaru1122.hatenadiary.jp
変数ｘというとわかりにくいので、顧客のサービス参加度（以下、参加度）としておきます。本当はもっと複雑で、Kobe.R では分析そっちのけでこの議論が少しあったのですが、ここでは簡潔に「参加度」としておきます。要は、あるサービスに関わろうとする度合いが強い顧客ほど満足度が高いのではないか、ということを調べます。
知りたいこと
顧客のサービス参加度は満足度（cs）に正の影響を与えているかどうか。
→参加度が高いほど、満足度も高くなるのでは？
分析の手順
１．まず相関をチェック。
２．順序ロジスティック回帰分析で、参加度の満足度への影響を調べる。
使用するデータ
使うデータは前回のエントリーと同じです。ダウンロードは↓から。
https://app.box.com/s/0jikm4cfeu4iopdhthfkyqlrw4z7ibed
実証分析
０．データの確認 まずクロス表で参加度と満足度の状況を確認します。たしかに、参加度が大きいほど、満足度も高い人が多そうです。
> 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～5の順序尺度、満足度は1～7の順序尺度です。連続値ではないので、厳密には普通の相関係数（ピアソンの積率相関係数）で調べることはできません。
5つや7つの順序尺度なら、連続値とみなして別にいいんじゃないかという意見もあります*1。私もそうしたいです、楽だし。でも、今回は学びのために、少し厳密になってみます。
3つ以上カテゴリを持つ順序カテゴリカル変数には、ポリコリック相関係数を調べます。cs とｘを順序カテゴリカル変数に変換しておき、データの構成を確認します。
> 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
２．順序ロジスティック回帰
正の相関があることが確認できました。顧客の参加度（1～5）がサービスの満足度（1～7）に影響を与えているかどうかを調べます。どちらも間隔尺度なので、本来は普通の回帰分析を使うことはできません*2。応答変数が順序カテゴリカル変数の場合、順序ロジスティック回帰が適切です。説明変数もカテゴリカルデータです。ダミー変数を作って、分析をしてみることにします。
まずデータの準備です。少し回りくどいですが、初めにcs2 ではなく、cs を抜き出したのは、makedummies() が順序カテゴリカル変数をすべてダミー変数化してしまうためです*3。顧客の参加度のみダミー変数を作りたいため、順序尺度ではない cs を抜き出しました。
> library(dplyr) ## データｄから 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()が便利なのでしょうか。
Rによるカテゴリカルデータ分析事例(2) ～(1) の補足：等分散性を確認していたか？～
前回のエントリーを書いたあと、現在は手に入りにくい本『Rによる統計解析の基礎』（私は買いました）の最新版『Rによる保健医療データ解析演習』を確認し、さらに「マイナーだけど最強の統計的検定 Brunner-Munzel 検定」を読むと、ウィルコクソンの順位和検定（マン・ホイットニーのU検定）*1とクラスカル・ウォリス検定には等分散性の確認が必要のようです。
Rによる統計解析の基礎 (Computer in Education and Research)
- 作者: 中澤港
- 出版社/メーカー: ピアソンエデュケーション
- 発売日: 2003/10
- メディア: 単行本
- 購入: 2人 クリック: 28回
- この商品を含むブログ (7件) を見る
『Rによる保健医療データ解析演習』は↓の「ダウンロード」セクションでゲットできます。とても勉強になります。
http://minato.sip21c.org/msb/index.html
念のため、分散の同質性を確認します。この確認にはフリグナー・キリーン検定を行います。デフォルトで使える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 検定を使います。
