データ穴リストのブログ

文系へっぽこリサーチャーから穴リストへジョブチェンジ中。渋谷の片隅でニヤニヤしております。・・・動いたら、ね。

Rで線形判別、の前に線形回帰、で、つまずく。

さて、本題に入ります。

 

が、何を書こうか決めてなかった。

とりあえず最近*1仕事でやってることにするか。

 

最近、「製品の上市前にある程度ウケるかどうかチェックしようよ」みたいな流れになってまして、

ようやく予測的なお話ができる・・・!

 

 

・・・んですが、元のデータはユーザビリティテストの値。

そして学習できるデータのある過去製品は…たぶん10製品くらい。

すくネェェエェェェエエエ

 

………と、世の中で活躍されているデータサイエンティストさんたち(呼称に思うところのある方々もいらっしゃると思いますが)は驚かれるかもしれませんが、

…世の中こんなもんですよね?

前時代的と言われても泣かない。*2

 

とりあえず線形判別分析

まぁデータの話はさておき、

正直初めてなんですよね〜、判別分析って。

 

本で読んでても大体

「うんうん、回帰モデルと一緒でしょ?

 まぁ機械学習的な発想も(しらんけど)この辺から来てんだよね。」

くらいの理解でおりました

(そして大概読み飛ばしてました)。

だって今まで予測関係全く必要なかったんだもの。

 思えば座学を受けてもその辺は結構スルーだった気がする。

 

 

というわけで、「いまさら聞けない判別分析」、やってみます。

 

で、まず判別分析って判別式は書くけどそのモデルの評価ってどこでやってるの?

ってところでまずつまづきました。はやい!

 

だって、誤判別率やらクロスバリデーションやらって話は、モデルができてからの話でしょ?

そのモデルの統計的有意性については一言も言ってないよね?(え、違うのかな?)

 

の前に、線形回帰

じゃぁってんで、とりあえず回帰から行ってみよう、と。

ここの変数が有意かどうか、モデルが有意かどうか確認してから判別してみましょうかと。

 

しかし判別分析のフローって、こういうので適切なんだろうか???

そもそもそういう、モデルの探索方法だとか、探索まではさておき、分析手法の組み合わせ方とかってほんと説明少ないですよね〜。

 

調査あるあるな「因子分析⇒クラスタ化」も、分析手法の解説本にはそれぞれ書いてあるけど、フローとしては書かれていなかった(まぁそのフローも批判されてるけど)。

 

そんなわけで、とりあえず線形回帰をやってみる。

 コードをハイライトする、みたいなやり方がよく分からないので、とりあえず分けとく。

>|r|
#まずは線形回帰モデルにしてみる。変数は色々投入#####
ach.model1 <- lm(income ~ obj1 * obj2, data = train.data)
ach.model2 <- lm(income ~ obj1 + obj2, data = train.data)
ach.model3 <- lm(income ~ obj1 + obj2 + obj4, data = train.data)
ach.model4 <- lm(income ~ obj1 * obj2 + obj4, data = train.data)
ach.model5 <- lm(income ~ obj1 + obj2 + obj4 + obj3 , data = train.data)
ach.model6 <- lm(income ~ obj1 * obj2 + obj4 + obj3 , data = train.data) #なんかこれだめ
ach.model7 <- lm(income ~ obj1 + obj3, data = train.data)
ach.model8 <- lm(income ~ obj1 * obj3, data = train.data)
ach.model9 <- lm(income ~ obj1 * obj3 + obj2, data = train.data)
ach.model10 <- lm(income ~ obj1 * obj3 * obj2, data = train.data)
ach.model11 <- lm(income ~ obj1 * obj3 + obj2 -1 , data = train.data)

とりあえず上のはあまり参考にはなりません。

色々投入してみたというだけ。

 

#モデルを1個ずつ確認。個々に見るのはこれで良いが、多いといちいちめんどくさい。
summary(ach.model1)
summary(ach.model2)
summary(ach.model3)
summary(ach.model4)
summary(ach.model5)
summary(ach.model6)
summary(ach.model7)
summary(ach.model8)
summary(ach.model9)
summary(ach.model10)
summary(ach.model11)
summary(ach.model12)
summary(ach.model13)

 

で、とりあえず各独立変数、あるいは予測因子、

それぞれを色んなパターンで投入して、各係数のp値を確認。

 

こんな結果が出た。

  • F検定が有意でないモデル
  • 定数項(α)が有意でないモデル
  • モデルの一部が部分的に有意でないモデル(にほんごがおかしい)
  • 諸々有意だけど決定係数が異常に高いモデル
  • 定数項(α)を除くと有意になるモデル

とりあえず、まず言いたい、

モデルいっぱいあって見づらい

 

色んなパターン試したけど、何回もやってると「あれ、コレ前に試した時どうだったんだっけ?」とか、思い出せなくなる。。。

世のアナリストさん達は少なくとも10分前に走らせたサマリーくらい覚えておられるのだろうか。

残念ながら僕の脳内メモリはそんなに優秀ではない。

 

こまけぇことはいいんだよ!

とりあえず、F検定が有意でないやつ、

決定係数が非常に小さいやつ、

どうだったか概観できるようにしときたい。

 

ということで、こんな風にしてみた。

 

# 上の処理を一括で行なってデータフレームにしておく#####
rm(model.fit)
model.fit <- as.data.frame(NULL)
rr <- as.vector(NULL)
adjrr <- as.vector(NULL)
p.val <- as.vector(NULL)
for(i in 1:13){ #最大値はモデルの数で。
samari <- summary(get(paste("ach.model",i,sep=""))) #モデルを一つずつ投入

rr <- c(rr, samari$r.squared) #決定係数
adjrr <- c(adjrr, samari$adj.r.squared) #自由度調整済み決定係数
p.val <- c(p.val, round((1-pf(samari$fstatistic["value"],samari$fstatistic["numdf"],samari$fstatistic["dendf"])), digits = 3))

model.fit <- cbind(rr,adjrr,p.val)
}
#できた! けど、なんか変なValueってのがいてデータフレームになってない。
class(model.fit)

model.fit <- as.data.frame(model.fit) #データフレームにそのままする。そのままにしておくと怒られる。
class(model.fit)

row.names(model.fit) <- NULL #行番号「Value」なるものを消す。

model.fit #できた。 

 

僕がようやくRをやる気になれた本、

「みんなのR」には、著者が作って自画自賛しておるcoefplotパッケージで各回帰係数を複数モデルで比較するためのmultiplot(obj1,obj2,...objn)でモデルの比較するよー

と書いてあるが、

#回帰係数coefficientを確認
require(coefplot)
coefplot(ach.model1)
coefplot(ach.model2)
coefplot(ach.model3)

multiplot(ach.model1,ach.model2,ach.model3,ach.model4,ach.model5,ach.model6,ach.model7,ach.model8,ach.model9,ach.model10,ach.model11,ach.model12,ach.model13) 

 

これ自体はすごく便利だと思う。

これで分析はかなり捗るだろうことは想像に難くない*3

でもこの前段階、ゴミみたいなモデルをゴミだと覚えておきたいの。

いや、覚えたくないから覚えさせておきたいの。

 

この辺り、もう少し読み進めれば適切な対応が出来るのだろうか。。。

 

さて、それはさておき、候補になりそうなモデルが出来た。

 

次は正規性の確認とかをプロットして、、、

というのは書くの面倒なので、省略します。

その辺はどっかのテキストで見てね!

 #4つのプロット。意味はよく分からんがこの辺はテキストに書いてあった。

par(mfrow=c(2,2))
plot(ach.model8)
par(mfrow=c(1,1)) #戻しとく。これを書いてくれないのはいけない。 

 

いよいよ判別

 

判別にはMASSパッケージを使います。

これはもう特に特徴ありません。丸コピです。

 

 require(MASS)

#判別式を作成。事前確率、グループ平均、係数がそれぞれ出る。
(ach.lda1 <- lda(check ~ lgi.cho_goo * lgi.repfeel.over3, data = train.data))
(ach.lda3 <- lda(check ~ lgi.cho_goo + lgi.repfeel.over3 + lgi.photogenic, data = train.data))
(ach.lda8 <- lda(check ~ lgi.cho_goo * lgi.no1, data = train.data))
(ach.lda9 <- lda(check ~ lgi.cho_goo * lgi.no1 + lgi.repfeel.over3, data = train.data))
(ach.lda11 <- lda(check ~ lgi.cho_goo * lgi.no1 + lgi.repfeel.over3 -1 , data = train.data))
(ach.lda12 <- lda(check ~ lgi.cho_goo * lgi.no1-1, data = train.data)) 

 

モデルが出来たら、予測。これも特に考えることはない。

まぁ、線形判別でなければこの辺でどの分布に当て込むか、methodの指定かなにかをするんだろうけど。

今はとりあえずのカタチを作ることが先決。

話はそれからだ。

 

#判別モデルをデータに当てはめ。予測結果を計算
(ach.Result1 <- predict(ach.lda1,train.data))
(ach.Result3 <- predict(ach.lda3,train.data)) #これだけ文句を言われた(1~4までのxの値がNA?)
(ach.Result8 <- predict(ach.lda8,train.data))
(ach.Result9 <- predict(ach.lda9,train.data))
(ach.Result11 <- predict(ach.lda11,train.data))
(ach.Result12 <- predict(ach.lda12,train.data)) 

それで、確認。

#予測結果と元データの結果を比較
table(train.data$check,ach.Result1$class) #これだけ誤判別が発生。他は完全に判別?
table(train.data$check,ach.Result3$class)
table(train.data$check,ach.Result8$class)
table(train.data$check,ach.Result9$class)
table(train.data$check,ach.Result11$class)
table(train.data$check,ach.Result12$class) 

さてさて、何か色々出てきますが、

説明読んでも忘れるから書いとこう。

 

さぁ予測が出来た!

バンザイ!

おしまい!

 

とはいかない!

 

我々の目的は何だったか?

そう、「なんかそれっぽい絵を作ってそれらしいことを言ってニヤける(ついでにそれなりに褒めてもらう)こと」

でしたね?

 

というわけで、次はビジュアライズです。

 

が、長くなったので後編につづく。 

 

 

*1:とかいいつつ、とりあえず記事書いてから外に出せるコードに直してって暇を見つけながらーと思ったら2か月くらい放置してしまったorz

*2:開発自体にかなりの時間を要する製品で、かつ製品間で統一したデータを取り出したのがそこそこ最近ってことで。。。無理やり合わせることもまぁできなくもないけど。。。

*3:ただね、、、デスクトップだとちゃんと動くんだけど、ノートだとエラーが出てきちんと動いてくれないんだよね。。作業環境は合わせているはずなんだけど。。勉強会したメンバーでも分かれてたから、何か原因があるんだと思うけど。デバッグしてみないといけないのかな。。