目次

統計工学実験 第3週

連絡

最終課題

準備

準備。 

Sys.setenv("http_proxy"="http://130.153.8.66:8080/")
install.packages(c("mvtnorm", "mvpart"), dependencies=TRUE)
library(mvtnorm)
library(mvpart)
library(MASS)

データがない人は、次のコードも実行する必要がある。 

tic.learn <- read.table("http://kdd.ics.uci.edu/databases/tic/ticdata2000.txt")
tic.eval <- read.table("http://kdd.ics.uci.edu/databases/tic/ticeval2000.txt")
tic.test <- read.table("http://kdd.ics.uci.edu/databases/tic/tictgts2000.txt")
tic.eval <- cbind(tic.eval, tic.test)
colnames(tic.eval)[86] <- "V86"
rm(tic.test)

さらに先週のレポート課題と同様に、V1とV5のグループ化も済ませるとよいかもしれない。 あるいは今週の課題で改めて、グループ化を見直してもよい。

クラスタリング

データの層別を行う方法の総称。生成した層のことを、クラスタ(群)という。大きく階層クラスタリングと非階層クラスタリングに分かれる。

階層クラスタリングは、レコード間の距離を全ての組み合わせについて算出して作成した行列を距離行列から出発する。 まずは最も近いレコード同士をグループ(群)にまとめる。その結果、2レコードは1つのクラスタに属するが、残りのレコードは点のままである。 そのクラスタと残りのレコードの間の距離行列を算出し、最も近いレコード同士、最も近いレコードとグループ、もしくは最も近いグループ同士を再びグループに
まとめる。 これを繰り返すことで、クラスタ(群)分けを得るのが、階層クラスタリングである。 階層クラスタリングには、レコード同士の距離だけでなく、レコードとグループの距離、グループ同士の距離、を定める必要がある。 またクラスタを生成する過程をグラフに表したものを、デンドログラムと呼ぶ。

それに対して、クラスタ(群)ごとの群内のばらつきが最小になるように、目的関数を定め、離散最適化法を用いてクラスタを生成するのが、非階層クラスタリング
である。 代表的な非階層クラスタリングのひとつである$k$-means法では、各群のユークリッド距離の平均が最小になるように、クラスタ平均を付置する。

k-means法によるクラスタリング

まず楕円状に分布するデータ(n=1000)を擬似乱数を用いて生成する。 

test.data <- rmvnorm(1000,mean=c(1,1),sigma=matrix(c(1,0.5,0.5,1),ncol=2))
plot(test.data)

これにk-means法をかけてみる。k-means法では、クラスタ数kは分析者が指定するものであり、ここでは2としてみる。さらに最小値を求めるために、ランダムに1000の初期値から出発(nstart=1000)して、最も適したクラスタ分けを採用するように指定してある。 

test.kmeans <- kmeans(test.data, centers=2, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster)

最初のグラフに色分けがついたはずで、これがクラスタである。 クラスタはデータのみから生成する層であり、これによる層別に技術的な背景はない。 クラスタ数を3, 4, 5と変えてみると、2の場合とは異なった分割が得られることが見て取れるはずである。

4枚のグラフをまとめて描くようにしてみた。 

par(mfrow=c(2,2))
test.kmeans <- kmeans(test.data, centers=2, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=2")
test.kmeans <- kmeans(test.data, centers=3, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=3")
test.kmeans <- kmeans(test.data, centers=4, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=4")
test.kmeans <- kmeans(test.data, centers=5, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=5")
par(mfrow=c(1,1))

同じことを、目視でもほぼ2群に分けられそうなデータ(n=2000)でも行ってみる。 

test.data <- rbind(rmvnorm(1000,mean=c(2.5,2.5),sigma=matrix(c(1,0.5,0.5,1),ncol=2)),
                   rmvnorm(1000,mean=c(-2.5,-2.5),sigma=matrix(c(1,0.5,0.5,1),ncol=2)))
plot(test.data)

そもそも別々に生成した2群のデータを結合しているので、k=2が一番自然なクラスタ分けを得る。 グラフを描かせるとやはりk=2の場合が最もらしいことが確認できる。 それでも3以上でもk-means法がクラスタを生成する。 

par(mfrow=c(2,2))
test.kmeans <- kmeans(test.data, centers=2, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=2")
test.kmeans <- kmeans(test.data, centers=3, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=3")
test.kmeans <- kmeans(test.data, centers=4, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=4")
test.kmeans <- kmeans(test.data, centers=5, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=5")
par(mfrow=c(1,1))

戯れにk=20としてみると、2つの群それぞれを、細かく群分けする様子が見て取れる。 

test.kmeans <- kmeans(test.data, centers=20, nstart=1000)
print(test.kmeans)
plot(test.data, col=test.kmeans$cluster, sub="k=20")

これはベクトル量子化と呼ばれる、多次元信号(画像など)の量子化法(デジタル信号への変換)と数理的には同等である。 また個々のクラスタごとに予測モデルを構築することは、学習理論における分割統治の考え方に合致している。

分割統治

この考え方をより進めると、下記のような手法に辿り着く。

層別

k-means法による層別

k-means法は、数値変数でしか用いられないため、データは「tic.learn[, c(1:85)]」と数値変数のフィールドを指定してkmeans()関数を実行する。 

tic.learn.kmeans <- kmeans(tic.learn[, c(1:85)], centers=2, nstart=1000)

k-means法の結果を適用するには、レコードごとに「centers」の座標とのユークリッド距離を算出し、一番小さくなる層にその変数を割り当てる。 

date()
tic.learn$cluster <- rep(0, dim(tic.learn)[1] )
for( i in c(1:dim(tic.learn)[1]) ) {
  tic.kmeans.dist <- rep(0, max(tic.learn.kmeans$cluster) )
  for( j in c(1:max(tic.learn.kmeans$cluster)) ) {
    tic.kmeans.dist[j] <- sum( (tic.learn[i,c(1:85)]-tic.learn.kmeans$center[j,])^2 )
  }
  tic.learn$cluster[i] <- sort.list(tic.kmeans.dist)[1]
}
date()

念のため、kmeansの結果と比較して、計算に誤りがないことを確認するためにクロス集計を行う。この結果が対角であれば、計算はあっている。 

table(tic.learn$cluster,tic.learn.kmeans$cluster)

centers=に指定する値は、いろいろ変えてみるといい。

変数の値を直接指定した層別(V1の例)

例えばV1で層別する場合、まずその変数を集計して、何種類の値があるかを確認する。 

table(tic.learn$V1)

V1は1から41の値を取ることが確認できたら、次のように、clusterという変数をtic.learn内に作り、まず0で埋める。 そして、V1の値ごとに層の値を指定するコードを書き、Rに実行させる。下記は、先週に行った決定木分析を踏まえた層別の例である。 

tic.learn$cluster <- rep(0, dim(tic.learn)[1])
tic.learn$cluster[tic.learn$V1==1] <- 1
tic.learn$cluster[tic.learn$V1==2] <- 2
tic.learn$cluster[tic.learn$V1==3] <- 1
tic.learn$cluster[tic.learn$V1==4] <- 2
tic.learn$cluster[tic.learn$V1==5] <- 2
tic.learn$cluster[tic.learn$V1==6] <- 1
tic.learn$cluster[tic.learn$V1==7] <- 2
tic.learn$cluster[tic.learn$V1==8] <- 2
tic.learn$cluster[tic.learn$V1==9] <- 2
tic.learn$cluster[tic.learn$V1==10] <- 2
tic.learn$cluster[tic.learn$V1==11] <- 2
tic.learn$cluster[tic.learn$V1==12] <- 1
tic.learn$cluster[tic.learn$V1==13] <- 2
tic.learn$cluster[tic.learn$V1==14] <- 2
tic.learn$cluster[tic.learn$V1==15] <- 2
tic.learn$cluster[tic.learn$V1==16] <- 2
tic.learn$cluster[tic.learn$V1==17] <- 2
tic.learn$cluster[tic.learn$V1==18] <- 2
tic.learn$cluster[tic.learn$V1==19] <- 2
tic.learn$cluster[tic.learn$V1==20] <- 2
tic.learn$cluster[tic.learn$V1==21] <- 2
tic.learn$cluster[tic.learn$V1==22] <- 2
tic.learn$cluster[tic.learn$V1==23] <- 2
tic.learn$cluster[tic.learn$V1==24] <- 2
tic.learn$cluster[tic.learn$V1==25] <- 2
tic.learn$cluster[tic.learn$V1==26] <- 2
tic.learn$cluster[tic.learn$V1==27] <- 2
tic.learn$cluster[tic.learn$V1==28] <- 2
tic.learn$cluster[tic.learn$V1==29] <- 2
tic.learn$cluster[tic.learn$V1==30] <- 2
tic.learn$cluster[tic.learn$V1==31] <- 2
tic.learn$cluster[tic.learn$V1==32] <- 2
tic.learn$cluster[tic.learn$V1==33] <- 2
tic.learn$cluster[tic.learn$V1==34] <- 2
tic.learn$cluster[tic.learn$V1==35] <- 2
tic.learn$cluster[tic.learn$V1==36] <- 2
tic.learn$cluster[tic.learn$V1==37] <- 2
tic.learn$cluster[tic.learn$V1==38] <- 2
tic.learn$cluster[tic.learn$V1==39] <- 2
tic.learn$cluster[tic.learn$V1==40] <- 2
tic.learn$cluster[tic.learn$V1==41] <- 2

最後に、0の値が残っていないかを確認するために、集計する。 

table(tic.learn$cluster)

変数の値を直接指定した層別(V2とV4の組み合わせの例)

例えばV2とV4の組み合わせで、下記のように層別することを考える。

V2/V4123456
1111133
2222233
3333333
4333333
5333333
6333333
7333333
8333333
9333333
10333333

この場合も、まずその変数を集計して、何種類の値があるかを確認する。 

table(tic.learn$V2, tic.learn$V4)

すべての組み合わせを書くのは煩雑なので、最初にcluster変数は3で埋めてしまう。 

tic.learn$cluster <- rep(3, dim(tic.learn)[1])

clusterが1になる箇所と2になる箇所だけ、書く。 

tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==1] <- 1
tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==2] <- 1
tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==3] <- 1
tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==4] <- 1
tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==1] <- 2
tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==2] <- 2
tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==3] <- 2
tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==4] <- 2

最後に、確認のために、集計する。 

table(tic.learn$cluster)

検証用データの層別

検証用データに層別変数 cluster を加えるには、上のコードの「tic.learn」をすべて「tic.eval」に変えればよい。

層別後の解析

層別後の解析

層ごとに回帰分析を適用したり、決定木を適用したりする。 

for( i in sort(unique(tic.learn$cluster)) ) {
  tic.learn.cluster <- tic.learn[tic.learn$cluster == i,]
  tic.working.lm <- lm(V86~V2+V3+V4, data=tic.learn.cluster)
  print(summary(tic.working.lm))
}

全体に回帰分析を行うのに較べて、残差の標準偏差 (Residual standard error) が小さくなっている層と、大きくなっている層があることを確認できるだろうか。 

print(summary(lm(V86~V2+V3+V4, data=tic.learn)))

考えたルールに基づく対象限定 (層別なし)

各変数に閾値を設けてルールを生成したとする。 たとえば、「V47が5.5以上かつV44が1未満」または「V47が5.5以上かつV1が{1,3,6,8,12,20}のどれか」、というルールは 次のように記す。 

(tic.learn$V47>5.5 & tic.learn$V44<1) | (tic.learn$V47>5.5 & (tic.learn$V1==1 |tic.learn$V1==3 | tic.learn$V1==6 | tic.learn$V1==8 | tic.learn$V1==12 | tic.learn$V1==20) )

「&」が「かつ(AND)」、「|」が「または(OR)」である。

このルールを検証用データに適用するには、 

tic.learn$visit <- (tic.learn$V47>5.5 & tic.learn$V44<1) | (tic.learn$V47>5.5 & (tic.learn$V1==1 |tic.learn$V1==3 | tic.learn$V1==6 | tic.learn$V1==8 | tic.learn$V1==12 | tic.learn$V1==20) )

と、訪問するか否かを二値(TRUE, FALSE)で表すオブジェクトを生成する。 このモデルに予測に基づいた訪問の成果を検証するには、訪問対象のリストtic.learn$visitと正解V86のクロス集計を行えばよい。 <code> table(tic.learn$visit) table(tic.learn$visit, tic.learn$V86) </code>

参考までに、上の訪問ルールを評価用のデータtic.evalに適用した場合、次のようになる。 

table(tic.eval$visit)
FALSE  TRUE 
 3029   971 

table(tic.eval$visit, tic.eval$V86)
tic.eval.visit    0    1
         FALSE 2878  151
         TRUE   884   87

ここでは、訪問対象に884+87=971人を選定し、そのうちの87人が実際に契約してくれる人だったことになる。 契約率は87/971=8.96%。また誤判別率は 

(884+151)/4000

で25.9%となる。

回帰分析に基づく対象限定 (層別なし)

まず、すべての変数を用いて回帰分析を行う。 

tic.lm <- lm(V86~V1gr+V2+V3+V4+V5gr+V6+V7+V8+V9+V10
            +V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
            +V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
            +V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
            +V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
            +V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
            +V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
            +V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
            +V81+V82+V83+V84+V85, data=tic.learn)
print(tic.lm)
summary(tic.lm)
plot(tic.lm)

上で「object 'V1gr' not found」などとエラーが表示されたら、先週のデータの読み込みからやり直し。 この結果をtic.lmに収めたら、次はAICで変数選択をさせる。 

tic.lm.aic <- stepAIC(tic.lm)

推定したモデルに基づいて、個別のレコードごとの契約の予測を立てる。その予測値が0.05以上なら訪問してみることとした場合、次のようになる。 

tic.learn$visit <- predict(tic.lm.aic, newdata=tic.learn)>0.05

集計してみると。 

table(tic.learn$visit)
table(tic.learn$visit, tic.learn$V86)

回帰分析に基づく対象限定 (層別あり)

学習したモデルに基づいて、訪問対象を狭めるには、predict()という関数を用いて、訪問対象か否かというリストを作成する。 まずクラスタ1について、設定まで調整したモデルを、学習用データ(tic.learn)から得る。(rpart関数に指定してあるモデルは適当なので、各自の得たものを用いること) 

tic.lm.1 <- stepAIC(lm(V86~., data=tic.learn[tic.learn$cluster==1, ]))

次に、このモデル(ここではtic.rpart.1)を、そのまま学習用データ(tic.learn)に適用して、契約してくれるか否かの予測を行う。 この際、下に記した0.05という閾値は調整(増やしたり減らしたり)が必要かもしれない。 

tic.learn$visit <- rep(0, dim(tic.learn)[1])
tic.learn$visit[tic.learn$cluster==1] <- predict(tic.lm.1, newdata=tic.learn[tic.learn$cluster==1,])>0.05

同様のことを、クラスタごとに、すべてのクラスタについて行う。(ここではk=2) 

tic.lm.2 <- stepAIC(lm(V86~., data=tic.learn[tic.learn$cluster==2, ]))
tic.learn$visit[tic.learn$cluster==2] <- predict(tic.lm.2, newdata=tic.learn[tic.learn$cluster==2,])>0.05
tic.lm.3 <- stepAIC(lm(V86~., data=tic.learn[tic.learn$cluster==3, ]))
tic.learn$visit[tic.learn$cluster==3] <- predict(tic.lm.3, newdata=tic.learn[tic.learn$cluster==3,])>0.05
...

この層別解析の結果をまとめる。 

table(1*(predict(tic.lm.1)>0.05), tic.learn$V86[tic.learn$cluster==1])
table(1*(predict(tic.lm.2)>0.05), tic.learn$V86[tic.learn$cluster==2])
table(1*(predict(tic.lm.3)>0.05), tic.learn$V86[tic.learn$cluster==3])
...

全体の集計は、上のテーブルを足しても良いし、集計し直しても良い。 

table(tic.learn$visit, tic.learn$V86)

決定木に基づく対象限定 (層別なし)

学習したモデルに基づいて、訪問対象を狭めるには、predict()という関数を用いて、訪問対象か否かというリストを作成する。 まず、設定まで調整したモデルを、学習用データ(tic.learn)から得る。 

tic.rpart <- rpart(V86~., data=tic.learn, control=c(cp=0.005))

次に、このモデル(ここではtic.rpart)を、そのまま学習用データ(tic.learn)に適用して、契約してくれるか否かの予測を行う。 この際、下に記した0.05という閾値は、確率が0.05以上なら訪問する、という意味なので、調整(増やしたり減らしたり)が必要かもしれない。 

tic.learn$visit <- predict(tic.rpart, newdata=tic.learn)[,2]>0.05

このモデルに予測に基づいた訪問の成果を検証するには、訪問対象のリストtic.visitと検証用データの正解V86のクロス集計を行えばよい。 

table(tic.learn$visit)
table(tic.learn$visit, tic.learn$V86

参考までに、評価用データに適用した場合は、次のようになる。 

tic.eval$visit <- predict(tic.rpart, newdata=tic.eval)[,2]>0.05
table(tic.eval$visit)

tic.eval$visit 
FALSE  TRUE
 2389  1611

table(tic.eval$visit, tic.eval$V86)

tic.eval$visit    0    1
         FALSE 2310   79
         TRUE  1452  159

ここでは、訪問対象に1452+159=1611人を選定し、そのうちの159人が実際に契約してくれる人だったことになる。契約率は159/1452=11.0%。 また誤判別率は 

(79+1452)/4000

で38.275%となる。

決定木に基づく対象限定 (層別あり)

学習したモデルに基づいて、訪問対象を狭めるには、predict()という関数を用いて、訪問対象か否かというリストを作成する。 まずクラスタ1について、設定まで調整したモデルを、学習用データ(tic.learn)から得る。(rpart関数に指定してあるモデルは適当なので、各自の得たものを用いること) 

tic.rpart.1 <- rpart(V86~., data=tic.learn[tic.learn$cluster==1, ], control=c(cp=0.005))

あるいは 

tic.rpart.1 <- rpart(V86~., data=subset(tic.learn, (cluster==1) ), control=c(cp=0.005))

とsubset()関数を用いて、データの一部を用いるように指定する。

次に、このモデル(ここではtic.rpart.1)を、そのまま学習用データ(tic.learn)に適用して、契約してくれるか否かの予測を行う。 この際、下に記した0.05という閾値は調整(増やしたり減らしたり)が必要かもしれない。 

tic.learn$visit <- rep(0, dim(tic.learn)[1])
tic.learn$visit[tic.learn$cluster==1] <- predict(tic.rpart.1, newdata=tic.learn[tic.learn$cluster==1,])>0.05

同様のことを、クラスタごとに、すべてのクラスタについて行う。(ここではk=2) 

tic.rpart.2 <- rpart(V86~., data=tic.learn[tic.learn$cluster==2, ], control=c(cp=0.005))
tic.learn$visit[tic.learn$cluster==2] <- predict(tic.rpart.2, newdata=tic.learn[tic.learn$cluster==2,])>0.05
tic.rpart.3 <- rpart(V86~., data=tic.learn[tic.learn$cluster==3, ], control=c(cp=0.005))
tic.learn$visit[tic.learn$cluster==3] <- predict(tic.rpart.3, newdata=tic.learn[tic.learn$cluster==3,])>0.05
...

この層別解析の結果をまとめる。 

table(1*(predict(tic.rpart.1)>0.05), tic.learn$V86[tic.learn$cluster==1])
table(1*(predict(tic.rpart.2)>0.05), tic.learn$V86[tic.learn$cluster==2])
table(1*(predict(tic.rpart.3)>0.05), tic.learn$V86[tic.learn$cluster==3])
...

全体の集計は、上のテーブルを足しても良いし、集計し直しても良い。 

table(tic.learn$visit, tic.learn$V86)

ここも層別なしの解析と比較しておく。

検証

以上のように、層別も加えて分析して作り上げたモデルが、目標を達成しているかどうか、を検証用データ (tic.eval) で確認する。

レポート提出について

レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること

項目指定
提出期限実験実施の翌週の月曜日の午後7時0分まで (でも今回は12月25日から冬休みなので、12月25日の午後7時でいい)
提出方法電子メールに添付 (宛先は配付資料に記載)
ファイル形式Wordファイル (LaTeXで作成する場合は、dvipdfmxでPDFに変換すること)
メールの件名統計実験3レポート提出(XXXXXXX)
レポートファイルの名称統計実験3_XXXXXXX.doc あるいは 統計実験3_XXXXXXX.docx
提出部数レポートは各自1通ずつ。レポートの表紙に、共同実験者の学籍番号と氏名を記すこと。

実験当日の修正

自分のMacだと問題が発生せず、実験室のPCだとRが落ちる不具合があります。メモリが少ないためかもしれません。Rが勝手に終了する場合には、kmeansのコードは次のように nstart を削除し、iter.max を追加する修正をお願いします。

k-means法による層別 (iter.max追加、nstart削除)

k-means法は、数値変数でしか用いられないため、データは「tic.learn[, c(1:85)]」と数値変数のフィールドを指定してkmeans()関数を実行する。 

tic.learn.kmeans <- kmeans(tic.learn[, c(1:85)], centers=2, iter.max=100)

k-means法の結果を適用するには、レコードごとに「centers」の座標とのユークリッド距離を算出し、一番小さくなる層にその変数を割り当てる。 

date()
tic.learn$cluster <- rep(0, dim(tic.learn)[1] )
for( i in c(1:dim(tic.learn)[1]) ) {
  tic.kmeans.dist <- rep(0, max(tic.learn.kmeans$cluster) )
  for( j in c(1:max(tic.learn.kmeans$cluster)) ) {
    tic.kmeans.dist[j] <- sum( (tic.learn[i,c(1:85)]-tic.learn.kmeans$center[j,])^2 )
  }
  tic.learn$cluster[i] <- sort.list(tic.kmeans.dist)[1]
}
date()

念のため、kmeansの結果と比較して、計算に誤りがないことを確認するためにクロス集計を行う。この結果が対角であれば、計算はあっている。 

table(tic.learn$cluster,tic.learn.kmeans$cluster)

centers=に指定する値は、いろいろ変えてみるといい。