post Image
『Rによるデータサイエンス』をcaretで再現する

Rによるデータサイエンス-データ解析の基礎から最新手法まで

にある各実行をRのcaretパッケージで試します。

対象外

  • caretは教師あり学習の判別と予測がメインなので、教師なし学習や生存分析、時系列モデルなどは対象外。
  • 対象外:第1-6章、第11-12章、第17章

全体のメモ

  • lmなど通常の関数で作成したモデルに対応するcaret使用のモデルは $finalModel となる。
  • ただ一部項目は異なるようです。
    • 例 lm
      • cars.lm<-lm(dist ~ speed, data = cars)
      • cars.caret.lm<-train(dist ~ speed, data = cars, method = ‘lm’)
      • cars.lm に対応するのは cars.caret.lm$finalModel
        • 予想結果等は一致する
        • モデルの統計値等の項目が異なる場合があるので、図示や統計値抽出の場合は要注意

第7章 線形回帰分析

7.2 線形単回帰分析 lm

library(caret)

# モデル作成
cars.lm<-train(dist ~ speed, data = cars, method = 'lm')
summary(cars.lm)

# 残差
residuals(cars.lm)
cars.lm$finalModel$residuals

# 回帰係数
round(coefficients(cars.lm$finalModel), 2)
round(cars.lm$finalModel$coefficients, 2)

# 予測
予測値<-predict(cars.lm)
残差<-residuals(cars.lm)
data.frame(cars, 予測値, 残差)

# グラフ
plot(cars)
abline(cars.lm$finalModel, lwd = 2)

par(oma = c(0, 0, 4, 0))
par(mfrow = c(2, 2))
plot(cars.lm$finalModel)
par(mfrow = c(1, 1))
par(oma = c(0, 0, 0, 0))

7.3 線形重回帰分析 lm

# モデル作成
air.lm<-train(Ozone ~ Solar.R + Wind + Temp, data = airquality, method = 'lm')
summary(air.lm)

# 回帰係数
round(coefficients(air.lm$finalModel), 2)

# 相互作用モデル
air.lm2<-train(Ozone ~ (Solar.R + Wind + Temp)^2, data = airquality, method = 'lm')
summary(air.lm2)

# 変数・モデルの選択
#air.lm3<-step(air.lm2$finalModel) はエラーとなる

第8章 非線形回帰分析

8.2 ロジスティック回帰、8.3 多項式回帰

  • 非線形回帰nlsは対応なし。

8.4 一般化線形モデル

air.glm1<-train(Ozone ~ Solar.R + Wind + Temp, data = airquality, method = 'glm', family = 'gaussian')
air.glm2<-train(Ozone ~ Solar.R + Wind + Temp, data = airquality, method = 'glm', family = 'Gamma')
AIC(air.glm1$finalModel)
AIC(air.glm2$finalModel)
  • 一般化線形モデルによるロジスティック回帰
年度<-c(1966:1984)
普及率<-c(0.003, 0.016, 0.054, 0.139, 0.263, 0.423, 0.611, 0.758, 0.859, 0.903, 0.937, 0.954, 0.978, 0.978, 0.982, 0.985, 0.989, 0.988, 0.992)
TV<-data.frame(年度, 普及率)

# モデル作成
fm5<-train(普及率 ~ 年度, data = TV, method = 'glm', family = 'binomial')
plot(年度, 普及率, type = 'l')
lines(年度, fitted(fm5), lty = 2, col = 'red', lwd = 2)
legend(locator(1), c('実測値', '予測値'), col = 1:2, lty = 1:2)

# 予測
predict(fm5$finalModel, type = 'response')

names(fm5); names(fm5$finalModel)  # どちらもglmで作成したモデルとは少し異なる

予測結果に直接影響ないが、caretを用いているため、fm5もfm5$finalModelもglmから直接求めたモデルとは結果の項目が異なっている。

8.5 平滑化回帰と加法モデル

  • 平滑化回帰と加法モデル
    比較のための多項式モデルは下記で実行可能だが、平滑化回帰は対応なし。
df<-data.frame(x1, y1)
lmp2<-fitted(train(y1 ~ poly(x1, 2), data = df, method = 'lm'))
lines(x1, lmp2, lty = 4, col = 4, lwd = 2)
lmp5<-fitted(train(y1 ~ poly(x1, 5), data = df, method = 'lm'))
lines(x1, lmp5, lty = 2, col = 2, lwd = 2)
lmp8<-fitted(train(y1 ~ poly(x1, 8), data = df, method = 'lm'))
lines(x1, lmp8)
legend(locator(1), c('poly(x1, 2)', 'poly(x1, 5)', 'poly(x1, 8)'), col = c(4, 2, 1), lty = c(4, 2, 1), lwd = 2)
  • 加法モデル
    • 平滑化スプラインの対応がないので、gamの値がテキストと異なる。
airq2<-na.exclude(airquality[, 1:4])
airq.gam<-train(Ozone ~ Solar.R + Wind + Temp, data = airq2, method = 'gam')
AIC(airq.gam$finalModel)
airq.glm<-train(Ozone ~ Solar.R + Wind + Temp, data = airq2, method = 'glm')
AIC(airq.glm$finalModel)

第9章 線形判別分析

9.3 ケーススタディ(lda)

iris.lab<-c(rep('S', 50), rep('C', 50), rep('V', 50))
iris1<-data.frame(iris[, 1:4], Species = iris.lab)
even.n<-2 * (1:75) - 1
iris.train<-iris1[even.n, ]
iris.test<-iris1[-even.n, ]

Z.lda<-train(Species ~., data = iris.train, method = 'lda')
Z.lda$finalModel

apply(Z.lda$finalModel$means%*%Z.lda$finalModel$scaling, 2, mean)

table(iris.train[, 5], predict(Z.lda))

#plot(Z.lda$finalModel$fitted, dimen = 1) プロットできない
#plot(Z.lda$finalModel$fitted, dimen = 2) プロットできない
Z<-predict(Z.lda$finalModel, iris.train[, -5])
plot(Z$x, type = 'n', ylim = c(-8, 8))
text(Z$x, labels = iris.train$Species)

Y<-predict(Z.lda$finalModel, iris.test[, -5])
table(iris.test[, 5], Y$class)

plot(Y$x, type = 'n')
text(Y$x, labels = iris.test$Species)

図9.2に対応するものはプロット出来ない。ldaとcaret/ldaで出力形式が違うためで、探せば見つかると思うが分からず
図9.3も同様に直接はプロットできなかった。図9.3(a)はpredictを用いてほぼ同等の図がプロット出来た。

9.4 交差確認

iris.CV<-train(Species ~., data = iris, method = 'lda', trControl = trainControl(method = 'LOOCV'))
(lda.tab<-table(iris[, 5], predict(iris.CV)))

sum(lda.tab[row(lda.tab) == col(lda.tab)]) / sum(lda.tab)
sum(lda.tab[row(lda.tab) != col(lda.tab)]) / sum(lda.tab)

第10章 非線形判別分析

10.2 判別関数による判別分析 (qda)

前章で作成したiris.train, iris.testデータを用いる

Z.qda<-train(Species ~., data = iris.train, method = 'qda')
table(iris.train[, 5], predict(Z.qda))

Y.qda<-predict(Z.qda, iris.test[, -5])
table(iris.test[, 5], Y.qda)

10.3 距離による判別分析

  • 対応なし

10.4 多数決による判別分析 (k最近傍法 k-Nearest Neighbor)

iris.knn<-train(Species ~., data = iris.train, method = 'knn', tuneGrid = expand.grid(k = 5))
iris.test.knn<-predict(iris.knn, iris.test[, -5])
table(iris.test.knn, iris.test[, 5])

# 図示
pairs(iris[, 1:4], pch = as.character(iris.test[ ,5]), col = c(3, 2)[(iris.test$Species != iris.test.knn) + 1])

Rplot.png

# 交差確認
iris.cv<-train(Species ~., data = iris, method = 'knn', trControl = trainControl(method = 'LOOCV'), tuneGrid = expand.grid(k = 5))
table(iris[, 5], predict(iris.cv))
  • 今回はk=5と指定されたが、kの最適値を求めるには下記のような感じ(1<=k<=10の範囲で探している)
knnGrid <- expand.grid(k = c(1:10))
iris.knn<-train(Species ~., data = iris.train, method = 'knn', tuneGrid = knnGrid)
iris.test.knn<-predict(iris.knn, iris.test[, -5])
table(iris.test.knn, iris.test[, 5])

10.5 ベイズ判定法

library(mlbench); data(Glass); G<-Glass[, c(1:5, 10)]

# NaiveBayesはデフォルトではカーネル法を用いないがcaretのnbはカーネル法を用いるのでオプションでFALSEにする
m1<-train(Type ~., data = G, method = 'nb', tuneGrid = data.frame(fL = 0, usekernel = F))
m1.p<-predict(m1)
tem1<-table(G$Type, m1.p)
1 - sum(diag(tem1)) / sum(tem1)

# カーネル法を用いる 
m1<-train(Type ~., data = G, method = 'nb')
m1.p<-predict(m1)
tem1<-table(G$Type, m1.p)
1 - sum(diag(tem1)) / sum(tem1)

第13章 樹木モデル

分類木

set.seed(20)
iris.rp<-train(Species ~., data = iris, method = 'rpart')
print(iris.rp$finalModel, digit = 1)

テキストの木と少し異なる。ただ、rpartで実行してもこの結果と同じ。テキストではmvpartを用いているので異なるか?ただ、mvpartはもうRでサポートされていない。

plot(iris.rp$finalModel, uniform = T, branch = 0.6, margin = 0.05)
text(iris.rp$finalModel, use.n = T, all = T)

木のプロットもmvpartを用いているテキストの図13.4とは異なるが、rpartのみの図と一致しているので、問題なし。

  • 木の剪定
printcp(iris.rp$finalModel)

iris.rp1<-prune(iris.rp$finalModel, cp = 0.02)
plot(iris.rp1, uniform = T, branch = 0.6, margin = 0.05)
text(iris.rp1, use.n = T)

mvpart使用より最初から枝が少ないので、変化なし。

# plotcp(iris.rp$finalModel) は動かず、、、

caretならtuneGridで行くのだろう。例えばこんなか

(iris.rp.test<-train(Species ~., data = iris, method = 'rpart', tuneGrid = expand.grid(cp = seq(0.01, 0.50, 0.01))))
  • 判別
even.n<-2 * (1:75) - 1
iris.train<-iris[even.n, ]
iris.test<-iris[-even.n, ]
set.seed(20)
iris.rp2<-train(Species ~., data = iris.train, method = 'rpart')
#plotcpは動かないので略
iris.rp3<-predict(iris.rp2, iris.test[, -5])
table(iris.test[, 5], iris.rp3)

回帰木

cars.rp<-train(dist ~ speed, data = cars, method = 'rpart')
cars.rp$finalModel

#printcpは略

cars.rp1<-prune(cars.rp$finalModel, cp = 0.044)
plot(cars.rp1, uniform = T, margin = 0.05)
text(cars.rp1, use.n = T)

多変量回帰木

  • mvpartのデータセットspiderが必要であることと、目的変数が複数である回帰はあまり頻度が高くないと思われ、略

13.4 捕遣と注釈

  • treeはなし
  • J48
iris.j48<-train(Species ~., data = iris, method = 'J48')
iris.j48$finalModel

library(party)
plot(iris.j48$finalModel)
plot(iris.j48$finalModel, type = c('extended'))

テキストでは図13.14の(a), (b)が出力されるとなっているが、extendedの指定に関わらず同じ図で、テキストとも表示がやや異なる。RWeka/J48でもcaretのj48でも同じ図。
Rplot01.png

  • ctree
cars.ctr<-train(dist ~., data = cars, method = 'ctree')
plot(cars.ctr$finalModel)
t.style<-node_hist(cars.ctr$finalModel, ymax = 0.06, xscale = c(0, 150), col = 'red', fill = hsv(0.6, 0.5, 1))
plot(cars.ctr$finalModel, terminal_panel = t.style)

第14章 ニューラルネットワーク

even.n<-2 * (1:75)
iris.train<-iris[even.n, ]
iris.test<-iris[-even.n, ]

iris.nnet<-train(Species ~., data = iris.train, method = 'nnet', tuneGrid = data.frame(size = 3, decay = 0.1))
iris.nnetp<-predict(iris.nnet, iris.test[, -5])
table(iris.test[, 5], iris.nnetp)

lvq

以下の通りに一応lvqは実施できるがlvqinitに対応するものや、lvq1, lvq2, lvq3, olvq1の別に対応するものはないようだ。

cd1<-train(Species ~., data = iris.train, method = 'lvq')
test.re<-predict(cd1, iris.test[, -5])
table(iris.test[, 5], test.re)

第15章 カーネル法とサポートベクターマシン

  • kpcaはなし。

判別

  • ksvm
library(kernlab); data(spam)
set.seed(50)
tr.num<-sample(4601, 2500)
spam.train<-spam[tr.num, ]
spam.test<-spam[-tr.num, ]

#set.seed(50) #テキストと同様にここで初期化をすると値がずれる。SVMの本質ではないのでどうでもよいがここでは初期化をなくした。

(spam.svm<-train(type ~., data = spam.train, method = 'svmRadial', trControl = 
trainControl(method = 'CV', number = 3)))
Support Vector Machines with Radial Basis Function Kernel 

2500 samples
  57 predictor
   2 classes: 'nonspam', 'spam' 

No pre-processing
Resampling: Cross-Validated (3 fold) 
Summary of sample sizes: 1667, 1666, 1667 
Resampling results across tuning parameters:

  C     Accuracy   Kappa    
  0.25  0.9044001  0.7960315
  0.50  0.9152011  0.8200693
  1.00  0.9215993  0.8345351

Tuning parameter 'sigma' was held constant at a value
 of 0.02769061
Accuracy was used to select the optimal model using 
 the largest value.
The final values used for the model were sigma =
 0.02769061 and C = 1. 
spam.pre<-predict(spam.svm, spam.test[, -58])
(spam.tab<-table(spam.test[, 58], spam.pre))
         spam.pre
          nonspam spam
  nonspam    1224   58
  spam         95  724
1 - sum(diag(spam.tab)) / sum(spam.tab)
0.07282247

テキストと値が異なるがksvmでも同じ値となる。

y<-as.matrix(iris[51:150, 5])
iris1<-data.frame(iris[51:150, 3:4], y)
set.seed(0)
ir.ksvm<-train(y ~., data = iris1, method = 'svmRadial')
#plot(ir.ksvm, data = iris1[, 1:2]) #動かず、 ir.ksvm$finalModelだと作図されるが違う図
table(iris1$y, predict(ir.ksvm, iris1[, 1:2]))
             versicolor virginica
  versicolor         47         3
  virginica           2        48

テキストと値が異なるがksvmでも同じ値となる。

回帰

x1<-seq(-10, 10, 0.1); set.seed(10)
y1<-50 * sin(x1) + x1^2 + 10 * rnorm(length(x1), 0, 1)
xy.svm<-train(y1 ~ x1, data = data.frame(x1, y1), method = 'svmRadial')
sy.pre<-predict(xy.svm, x1)
plot(x1, y1, type = 'l')
lines(x1, sy.pre, col = 'red', lty = 2)
legend(locator(1), c('実測値', '予測値'), lty = c(1, 2), col = c(1, 2))

epsilonを指定するメソッドがcaretのsvmにはないので略。

第16章 集団学習

テキストでは adabag パッケージを使用している。caretでは adabag パッケージを使用するものとして、AdaBag, AdaBoost.M1 があり、前者はバギングとブースティングの両方として紹介され、後者はブースティングとして紹介されている。ここでは前者をバギング、後者をブースティングとして示します。

バギング

library(mlbench)
data(BreastCancer)
even.n<-2 * (1:349)
BC.train<-BreastCancer[even.n, -1]
BC.test<-BreastCancer[-even.n, -1]
for(i in 1:9) {
  BC.train[, i]<-as.integer(BC.train[, i])
  BC.test[, i]<-as.integer(BC.test[, i])
}

# BC.train には欠損値があり、bagging ではそのまま扱えるようだが、下記ではそのまま扱えないので欠損値を埋める処理をする。
set.seed(20)
tmp.pre<-preProcess(BC.train, method = 'bagImpute')
BC.train<-predict(tmp.pre, BC.train)

set.seed(20)
BC.ba<-train(Class ~., data = BC.train, method = 'AdaBag')

# 同様にBC.test  にも欠損値があり埋める処理をする。
set.seed(20)
tmp.pre<-preProcess(BC.test, method = 'bagImpute')
BC.test<-predict(tmp.pre, BC.test)

BC.bap<-predict(BC.ba, BC.test)
BC.bap<-data.frame(BC.bap); names(BC.bap)<-'class'
(ta.ba<-table(BC.test[, 10], BC.bap$class))
            benign malignant
  benign       216        11
  malignant      3       120

テキストと値が異なってしまっている。欠損値の補完の問題だろうか、、、

1 - sum(diag(tb.ba)) / sum(tb.ba)
[1] 0.04

ブースティング

set.seed(20)
BC.ad<-train(Class ~., data = BC.train, method = 'AdaBoost.M1')

BC.adp<-predict(BC.ad, newdata = BC.test[, -10])
BC.adp<-data.frame(BC.adp); names(BC.adp)<-'class'
(tb.ad<-table(BC.adp$class, BC.test[, 10]))
            benign malignant
  benign       219         3
  malignant      8       120

テキストと値が異なってしまっている。欠損値の補完の問題かもしれない。

1 - sum(diag(tb.ad)) / sum(tb.ad)
[1] 0.03142857

ランダムフォレスト

data(BreastCancer)
even.n<-2 * (1:349)
BC.train<-BreastCancer[even.n, -1]
BC.test<-BreastCancer[-even.n, -1]
for(i in 1:9) {
  BC.train[, i]<-as.integer(BC.train[, i])
  BC.test[, i]<-as.integer(BC.test[, i])
}

BC.train<-na.omit(BC.train)

set.seed(20)
BC.rf<-train(Class ~., data = BC.train, method = 'rf')
print(BC.rf$finalModel)
Call:
 randomForest(x = x, y = y, mtry = param$mtry) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 2

        OOB estimate of  error rate: 3.85%
Confusion matrix:
          benign malignant class.error
benign       214         7  0.03167421
malignant      6       111  0.05128205

テキストとは値が異なってしまっている。テキストとは異なり、 mtry = param$mtry が指定されてしまっているためか。パラメータの指定でどうにかなるかも知れないが、ここではスルーします。

summary(BC.rf$finalModel)
                Length Class      Mode     
call               4   -none-     call     
type               1   -none-     character
predicted        338   factor     numeric  
err.rate        1500   -none-     numeric  
confusion          6   -none-     numeric  
votes            676   matrix     numeric  
oob.times        338   -none-     numeric  
classes            2   -none-     character
importance         9   -none-     numeric  
importanceSD       0   -none-     NULL     
localImportance    0   -none-     NULL     
proximity          0   -none-     NULL     
ntree              1   -none-     numeric  
mtry               1   -none-     numeric  
forest            14   -none-     list     
y                338   factor     numeric  
test               0   -none-     NULL     
inbag              0   -none-     NULL     
xNames             9   -none-     character
problemType        1   -none-     character
tuneValue          1   data.frame list     
obsLevels          2   -none-     character

caret使用時は、BC.rf$finalModel が randomForest のモデルに該当する以下も同様。

BC.rf$finalModel$type
[1] "classification"
plot(BC.rf$finalModel)
BC.rf$finalModel$importance
                MeanDecreaseGini
Cl.thickness          10.3031669
Cell.size             31.0735166
Cell.shape            25.1238751
Marg.adhesion         10.2548693
Epith.c.size          13.5882463
Bare.nuclei           20.6794913
Bl.cromatin           20.7562703
Normal.nucleoli       19.0439173
Mitoses                0.9920097
varImpPlot(BC.rf$finalModel)
BC.rfp<-predict(BC.rf, BC.test[, -10])
(BC.rft<-table(BC.test[, 10], BC.rfp))
           BC.rfp
            benign malignant
  benign       219         8
  malignant      3       120
(1 - sum(diag(BC.rft)) / sum(BC.rft))
[1] 0.03142857

『 機械学習 』Article List
Category List

Eye Catch Image
Read More

Androidに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

AWSに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Bitcoinに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

CentOSに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

dockerに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

GitHubに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Goに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Javaに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

JavaScriptに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Laravelに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Pythonに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Rubyに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Scalaに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Swiftに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Unityに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Vue.jsに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

Wordpressに関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。

Eye Catch Image
Read More

機械学習に関する現役のエンジニアのノウハウ・トレンドのトピックなど技術的な情報を提供しています。コード・プログラムの丁寧な解説をはじめ、初心者にもわかりやすいように写真や動画を多く使用しています。