#======================================================
#	7 ツリーモデルとrandomForest
#======================================================
#======================================================
#	7.1 ツリーモデル
#======================================================
#パッケージのロード
library(rpart)         #ツリーモデルのため
library(rpart.plot)    #プロットのための追加機能
library(MASS)          #stepAICのため

#---------------<回帰ツリー>---------------
#---------------<適用例：自動車>---------------
summary(car.test.frame)

#燃費を重量で推定
car.rpt <- rpart(Mileage ~ Weight, data=car.test.frame)
car.rpt #当てはめ結果の出力
#描画
#図7.1
prp(car.rpt, type=2, extra=1, digits=0)

#残差平方和（RSS）
(mu <- mean(car.test.frame$Mileage)) #平均
(RSS.1 <- sum((car.test.frame$Mileage - mu)^2)) #残差平方和

#分岐の例
mu.left <- mean(car.test.frame$Mileage
		[car.test.frame$Weight>=2800])        #左の集団の平均
(RSS.left <- sum((car.test.frame$Mileage
	[car.test.frame$Weight>=2800] - mu.left)^2)) #左の集団の残差平方和
mu.right <- mean(car.test.frame$Mileage
		[car.test.frame$Weight<2800])         #右の集団の平均
(RSS.right <- sum((car.test.frame$Mileage
	[car.test.frame$Weight<2800] - mu.right)^2)) #右の集団の残差平方和
(RSS.left + RSS.right) #合計
 
#最適な分岐
#中間点
mid <- (sort(car.test.frame$Weight)[1:59] + 
	sort(car.test.frame$Weight)[2:60])/2
#初期化
RSS.all <- numeric(59)

#各中間点ごとに、残差平方和の合計を計算
for(i in 1:59) {
    mu.left <- mean(car.test.frame$Mileage
		[car.test.frame$Weight>=mid[i]])
    RSS.left <- sum((car.test.frame$Mileage
	[car.test.frame$Weight>=mid[i]] - mu.left)^2)
    mu.right <- mean(car.test.frame$Mileage
		[car.test.frame$Weight<mid[i]])
    RSS.right <- sum((car.test.frame$Mileage
	[car.test.frame$Weight<mid[i]] - mu.right)^2)
    RSS.all[i] <- RSS.left + RSS.right		
}
mid[which.min(RSS.all)] #最適な分岐点
#描画
#図7.2
plot(mid, RSS.all, type="b", xlab="分岐点", ylab="残差平方和の合計")
abline(v = 2567.5, lty=2, col="red")

#複雑度パラメータ
car.rpt$frame

#検算
(1354.583-(361.2 + 186.93333))/1354.583

#予測（ツリーモデル）
pred <- predict(car.rpt)
ord <- order(car.test.frame$Weight, pred)
plot(Mileage ~ Weight, data=car.test.frame, pch=16)
lines(car.test.frame$Weight[ord], pred[ord])

#---------------<複数の連続な予測変数>---------------
#燃費を重量と排気量で推定
car.rpt2 <- rpart(Mileage ~ Weight + Disp., data=car.test.frame)
car.rpt2 #当てはめ結果の出力
#描画
#図7.3
prp(car.rpt2, type=2, extra=1, digits=0)

#---------------<カテゴリカルな予測変数を含む場合>---------------
#燃費を排気量と国別で推定
car.rpt3 <- rpart(Mileage ~ Disp. + Country, data=car.test.frame)
car.rpt3 #当てはめ結果の出力
#描画
#図7.4
prp(car.rpt3, type=2, extra=1, digits=0)

#---------------<適用例：Boston住宅価格>---------------
#データの確認
library(mlbench)
data(BostonHousing)

#サンプリング
set.seed(456)
n <- nrow(BostonHousing)  #観測値の数
n.train <- ceiling(n/2)   #観測値の数の約50％
train <- sample(1:n, n.train)

#当てはめ（ツリーモデル）
boston.rpt0 <- rpart(medv ~ ., data=BostonHousing[train, ],
        control=rpart.control(cp=0.0001))
boston.rpt0

#cpの検討
#一覧表
printcp(boston.rpt0)
(23299.870 - (8890.188 + 2880.716))/23299.870
(8890.188 + 2880.716)/23299.870
0.27293 + 0.055522
#描画
#図7.5
plotcp(boston.rpt0)

#当てはめ（ツリーモデル２）
boston.rpt1 <- rpart(medv ~ ., data=BostonHousing[train, ],
        control=rpart.control(cp=0.0148954))

#描画
#図7.6
prp(boston.rpt1, type=2, extra=1, digits=0)

#推定結果
boston.rpt1
printcp(boston.rpt1)

#予測と平均平方誤差
pred.rpt0 <- predict(boston.rpt0, newdata=BostonHousing[-train, ])
mse.rpt0 <- sum((BostonHousing$medv[-train] - pred.rpt0)^2)/
        (n - n.train)
pred.rpt1 <- predict(boston.rpt1, newdata=BostonHousing[-train, ])
mse.rpt1 <- sum((BostonHousing$medv[-train] - pred.rpt1)^2)/
        (n - n.train)

#当てはめ（線形モデル）
boston.lm0 <- lm(medv ~ . ^2, data=BostonHousing[train, ])

#変数選択
boston.lm1 <- stepAIC(boston.lm0)

#要約
summary(boston.lm1)

#予測と平均平方誤差
pred.lm <- predict(boston.lm1, newdata=BostonHousing[-train, ])
mse.lm <- sum((BostonHousing$medv[-train] - pred.lm)^2)/
        (n - n.train)

#比較
#ツリーモデル
#剪定前
mse.rpt0
#剪定後
mse.rpt1
#線形モデル
mse.lm

#---------------<分類ツリー>---------------

#複雑度又は不純度
#図7.7
p <- seq(0, 1, length=200)
gini <- 1 - (p^2 + (1-p)^2)　#Gini係数
info <- -p*log(p)-(1-p)*log(1-p) #情報エントロピー
plot(p, gini, ylab="Impurity", type="l", ylim=c(0, 0.7))
lines(p, info, lty=2)
legend("topleft", c("Gini", "Information"), lty=1:2)

#----------<適用例：乳がん>---------

#ツリーモデル
#データ
data(BreastCancer)　#データのロード
summary(BreastCancer)
BC <- na.omit(BreastCancer)　#NAを除去
BC <- BC[ ,-1]　#Idを除去

#サンプリング
set.seed(456)
n <- nrow(BC)
train <- sample(1:n, ceiling(n/2))

#モデルの当てはめ
BC.rpt <- rpart(Class ~ ., data=BC[train, ])
BC.rpt
printcp(BC.rpt)
#図7.8
prp(BC.rpt, type=2, extra=1, digits=0)

#予測
BC.pred <- predict(BC.rpt, newdata=BC[-train, ], type="class")
(ans.rpt <- table(BC$Class[-train], BC.pred))
error.rpt <- (sum(ans.rpt) - sum(diag(ans.rpt)))/sum(ans.rpt)

#決して実行しないでください。
#BC.glm0 <- glm(Class ~ (Cl.thickness + Cell.size + Cell.shape + Marg.adhesion
#	+Epith.c.size + Bare.nuclei + Bl.cromatin + Normal.nucleoli
#	+ Mitoses)^2, data=BC[train, ], family=binomial)

#こちらを実行する。
BC.glm0 <- glm(Class ~ ., data=BC[train, ], family=binomial, 
    maxit=50)
BC.glm1 <- stepAIC(BC.glm0)
summary(BC.glm1)
anova(BC.glm1, test="Chi")
BC.glm.pred <- predict(BC.glm1, newdata=BC[-train, ], type="response")
(ans.glm <- table(BC$Class[-train], BC.glm.pred>0.5))
error.glm <- (sum(ans.glm) - sum(diag(ans.glm)))/sum(ans.glm)

#比較
#ツリーモデル
error.rpt
#線形モデル
error.glm

#======================================================
#	7.2 randomForest
#======================================================
#パッケージのロード
library(randomForest)

#復元抽出
set.seed(456)
sample(1:10, 5, replace=T)

#---------------<回帰問題のためのrandomForest>---------------
#---------------<適用例：自動車>---------------

#サンプリングの準備
car.all <- car.test.frame[ ,-3] #NAがある変数を除く
n <- dim(car.all)[1]　          #サンプル数全体
n.train <- n - ceiling(n/2)　   #訓練データ

#ツリーモデル
mse.car.rpt <- numeric(100) #数値ベクトルの初期化
set.seed(456)
for (i in 1:100) { #100回ループ
       train <- sample(1:n, n.train)　#サンプリング
       car.rpt <- rpart(Mileage ~ ., data=car.all[train, ], 
	        control = rpart.control(minsplit=5, cp=0.0001))
       car.rpt.pred <- predict(car.rpt, newdata=car.all[-train, ])
       mse.car.rpt[i] <- sum((car.all$Mileage[-train] - car.rpt.pred)^2)/
               (n - n.train)
}
 
#線形モデル
mse.car.lm1 <- numeric(100) #数値ベクトルの初期化
set.seed(456)
for (i in 1:100) { #100回ループ
       train <- sample(1:n, n.train)　#サンプリング
       car.lm0 <- lm(Mileage ~ (.)^2, data=car.all[train,  -c(2, 4)])
       car.lm1 <- stepAIC(car.lm0, trace=F)
       car.lm1.pred <- predict(car.lm1, newdata=car.all[-train, ])
       mse.car.lm1[i] <- sum((car.all$Mileage[-train] - car.lm1.pred)^2)/
              (n - n.train)
}

#randomForest
mse.car.rf <- numeric(100) #数値ベクトルの初期化
set.seed(456)
for (i in 1:100) { #100回ループ
       train <- sample(1:n, n.train)　#サンプリング
       car.rf <- randomForest(Mileage ~ ., data=car.all[train, ], mtry=6)
       car.rf.pred <- predict(car.rf, newdata=car.all[-train, ])
       mse.car.rf[i] <- sum((car.all$Mileage[-train] - car.rf.pred)^2)/
              (n - n.train)
}

#MSE比較
#ツリーモデル
mean(mse.car.rpt)

#線形モデル
mean(mse.car.lm1)

#randomForest
mean(mse.car.rf)

#箱ひげ図
#図7.9
mse <- c(mse.car.rpt, mse.car.lm1, mse.car.rf)
model <- factor(rep(c("rpart", "lm", "randomForest"), each=100), 
	levels=c("rpart", "lm", "randomForest"))
boxplot(log(mse) ~ model, ylab="log(MSE)", col="springgreen")

#randomForestの結果表示
importance(car.rf)

#図7.10
plot(car.rf)

#---------------<分類問題のためのrandomForest>---------------
#---------------<適用例：乳がん>---------------

#サンプリング
set.seed(456)
n <- nrow(BC)
train <- sample(1:n, ceiling(n/2))

#モデルの当てはめ
x <- BC[ ,1:9]
y <- BC[ ,10]
BC.rf <- randomForest(x, y, data=BC[train, ], mtry=9)

#予測
new <- cbind(x, y)
BC.pred.rf <- predict(BC.rf, newdata=new[-train, ])
(ans.rf <- table(BC$Class[-train], BC.pred.rf))
#誤判別率
error.rf <- (sum(ans.rf) - sum(diag(ans.rf)))/sum(ans.rf)
error.rf

importance(BC.rf)

#---------------<適用例：DNA>---------------
#---------------<tree>---------------

#サンプリング
set.seed(456)
n <- nrow(dna)
train <- sample(1:n, ceiling(n/2))

#モデルの当てはめ
dna.rpt <- rpart(Class ~ ., data = dna[train, ])
printcp(dna.rpt)

#予測
dna.pred <- predict(dna.rpt, newdata=dna[-train, ], type="class")
(ans.rpt <- table(dna$Class[-train], dna.pred))
error.rpt <- (sum(ans.rpt) - sum(diag(ans.rpt)))/sum(ans.rpt)
error.rpt      #誤判別率

#---------------<randomForest>---------------

#推定用データの構成
x <- dna[ ,-41]
y <- dna$Class

#モデルの当てはめ
dna.rf <- randomForest(x, y, data=dna[train, ])

#予測
new <- cbind(x, y)
dna.pred.rf <- predict(dna.rf, newdata=new[-train, ])
(ans.rf <- table(dna$Class[-train], dna.pred.rf))

#誤判別率
error.rf <- (sum(ans.rf) - sum(diag(ans.rf)))/sum(ans.rf)
error.rf
