## -*- mode: R -*-
## 2019年10月11日(金)
##前回との変更点は、tree→rpartにしたこと。

## ##参照ファイルは、
## forPrint5rpart.pdf
## command5rpart.txt
## ・source
## ・read.table
## splice.All <- read.table("splice.dat",header=TRUE)

####################################
######  Treeパッケージの準備  ###### 
####################################
install.packages("rpart", repos="http://cran.ism.ac.jp/")
library(rpart)


############################
######  決定木の説明  ######
############################
setwd("C:\\Users\\student\\Desktop")
##splicingのデータ。
splice.All <- read.table("splice.dat",header=TRUE)
dim(splice.All)##3181個のデータ、塩基の数は、19。先頭はラベル({EI,IE,N})
splice.All[20,]
## 1番目の要素がlabel。全部で19塩基。
##   label s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19
##20    EI  T  C  T  C  C  A  A  G  G   T   G   A   G   A   T   C   A   C   C

summary(splice.All)
##  label     s1      s2      s3      s4      s5      ... s18     s19    
##  EI: 764   A:660   A:636   A:651   A:666   A:766   ... A:698   A:719  
##  IE: 766   C:964   C:951   C:944   C:889   C:882   ... C:816   C:789  
##  N :1651   G:701   G:704   G:703   G:657   G:854   ... G:924   G:913  
##            T:856   T:890   T:883   T:969   T:679   ... T:743   T:760  
       
##手始めに、ラベルが{EI}もしくは{N}を選んでやってみよう。|はorの事。
splice <- splice.All[splice.All$label=="EI"|splice.All$label=="N",]

dim(splice)##2415個のデータ、20塩基の数は、19。先頭はラベル({EI,N})
summary(splice)
##大体の様子が判る。
##  label     s1      s2      s3      s4      ...  s18     s19    
##  EI: 764   A:619   A:571   A:614   A:611   ...  A:509   A:525  
##  IE:   0   C:641   C:606   C:588   C:582   ...  C:618   C:607  
##  N :1651   G:614   G:649   G:659   G:603   ...  G:707   G:686  
##            T:541   T:589   T:554   T:619   ...  T:581   T:597  

##spliceのデータは、label,s1,s2,....,s18,s19 (1~19番目の塩基)の情報。
splice.tree <- rpart(label ~ ., data=splice,
                     parms = list(split = "information"))
##「spliceデータで、labelを予測してね」。「.」は、 label以外の属性値を
##用いて予測する「当て物関数」での、お約束。
##つまり、labelを、(1~19番目の塩基)の情報に基づいて予測してね。
splice.tree
## n= 2415 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
##  1) root 2415 764 N (0.316356108 0.000000000 0.683643892)  
##    2) s10=T 1177 420 EI (0.643160578 0.000000000 0.356839422)  
##      4) s9=G 855  98 EI (0.885380117 0.000000000 0.114619883)  
##        8) s13=G 665  15 EI (0.977443609 0.000000000 0.022556391) *
##        9) s13=A,C,T 190  83 EI (0.563157895 0.000000000 0.436842105)  
##         18) s8=G 127  23 EI (0.818897638 0.000000000 0.181102362) *
##         19) s8=A,C,T 63   3 N (0.047619048 0.000000000 0.952380952) *
##      5) s9=A,C,T 322   0 N (0.000000000 0.000000000 1.000000000) *
##    3) s10=A,C,G 1238   7 N (0.005654281 0.000000000 0.994345719) *
plot(splice.tree);text(splice.tree, use.n = TRUE, all = TRUE)

##EIと分類するのは、
##  8) s13=G 665  15 EI (0.977443609 0.000000000 0.022556391) *
## 18) s8=G 127  23 EI (0.818897638 0.000000000 0.181102362) *
path.rpart(splice.tree, node = c(8,18))
##  node number: 8 
##    root
##    s10=T
##    s9=G
##    s13=G
## 
##  node number: 18 
##    root
##    s10=T
##    s9=G
##    s13=A,C,T←Gじゃなくても
##    s8=G であれば、EIだよ。

##予測と実際の違い。newdataと言いつつ、学習データをつかっている。
splice.prd <- predict(splice.tree, newdata=splice, type="class")
##予測のラベルと本当のラベルの
table(splice.prd,splice$label)
## splice.prd   EI   IE    N
##         EI  754    0   38
##         IE    0    0    0
##         N    10    0 1613

diag(table(splice.prd,splice$label))##対角成分を取り出した。

##正解の数
sum(diag(table(splice.prd,splice$label)))
##全数
sum(table(splice.prd,splice$label))
##正解率は、98%
sum(diag(table(splice.prd,splice$label)))/sum(table(splice.prd,splice$label))

#############################
##全部のラベルでやってみよう。
#############################
splice.All <- read.table("splice.dat",header=TRUE)
splice <- splice.All##上書きして、全ラベル({EI,IE,N})使う。
splice.tree <- rpart(label ~ ., data=splice,
                     parms = list(split = "information"))

plot(splice.tree, branch=0.9, margin=0.04)##ちょっとだけ、枝の角度を調整する。
text(splice.tree, use.n=T)##Nodeに含まれる、{EI,IE,N}の数。
splice.tree
## n= 3181 前のよりもちょっと複雑。
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
##  1) root 3181 1530 N (0.240176045 0.240804778 0.519019176)  
##    2) s8=G 1818 1054 IE (0.342684268 0.420242024 0.237073707)  
##      4) s10=T 952  335 EI (0.648109244 0.261554622 0.090336134)  
##        8) s9=G 750  133 EI (0.822666667 0.144000000 0.033333333)  
##         16) s13=G 537   24 EI (0.955307263 0.040968343 0.003724395) *
##         17) s13=A,C,T 213  109 EI (0.488262911 0.403755869 0.107981221)  
##           34) s11=A 90   13 EI (0.855555556 0.088888889 0.055555556) *
##           35) s11=C,G,T 123   45 IE (0.219512195 0.634146341 0.146341463) *
##        9) s9=A,C,T 202   61 IE (0.000000000 0.698019802 0.301980198)  
##         18) s7=A 154   13 IE (0.000000000 0.915584416 0.084415584) *
##         19) s7=C,G,T 48    0 N (0.000000000 0.000000000 1.000000000) *
##      5) s10=A,C,G 866  351 IE (0.006928406 0.594688222 0.398383372)  
##       10) s7=A 622  109 IE (0.004823151 0.824758842 0.170418006)  
##         20) s6=C,T 543   45 IE (0.003683241 0.917127072 0.079189687) *
##         21) s6=A,G 79   16 N (0.012658228 0.189873418 0.797468354) *
##       11) s7=C,G,T 244    5 N (0.012295082 0.008196721 0.979508197) *
##    3) s8=A,C,T 1363  143 N (0.103448276 0.001467351 0.895084373)  
##      6) s13=G 421  138 N (0.327790974 0.000000000 0.672209026)  
##       12) s10=T 213   76 EI (0.643192488 0.000000000 0.356807512)  
##         24) s9=G 150   13 EI (0.913333333 0.000000000 0.086666667) *
##         25) s9=A,C,T 63    0 N (0.000000000 0.000000000 1.000000000) *
##       13) s10=A,C,G 208    1 N (0.004807692 0.000000000 0.995192308) *
##      7) s13=A,C,T 942    5 N (0.003184713 0.002123142 0.994692144) *

##予測
splice.prd <- predict(splice.tree, newdata=splice, type="class")

##予測のラベルと本当のラベルの
table(splice.prd,splice$label)
## splice.prd   EI   IE    N
##         EI  727   30   20
##         IE   29  717   74
##         N     8   19 1557

diag(table(splice.prd,splice$label))##対角成分を取り出した。
##正解率は、94%。(前のは2種類を当てるので簡単だったので98%)
sum(diag(table(splice.prd,splice$label)))/sum(table(splice.prd,splice$label))

##非常に細かいルールを作る。
splice.tree <- rpart(label ~ ., data=splice,
                     parms = list(split = "information"),
                     control = rpart.control(
                       minsplit=0,    # 分岐の最少データ数(小さいと細かく分岐)
                       cp=0,          # モデルの複雑さ(小さいと細かく分岐)
                       maxdepth = 30)) # 木の最大深さ(大きいと細かく分岐)

splice.prd <- predict(splice.tree, newdata=splice, type="class")
table(splice.prd,splice$label)
##正解率は、99%。お任せでは、94%
sum(diag(table(splice.prd,splice$label)))/sum(table(splice.prd,splice$label))

plot(splice.tree, branch=0.9, margin=0.04);text(splice.tree)



###些細なルールを作ってしまうと正しく判定できないことに陥る。
##7割のデータを選んで、学習。3割のデータで、テスト。
length(splice.All[,1])##全データの数

n <- length(splice.All[,1])
lrn <- sort(sample(n, round(n*0.7)))##7割の学習データを選ぶ。

splice <- splice.All[lrn,]
splice.tree <- rpart(label ~ ., data=splice,
                     parms = list(split = "information"),
                     control = rpart.control(minsplit=0, cp=0,
                       maxdepth = 30))
splice.prd <- predict(splice.tree, newdata=splice, type="class")

plot(splice.tree); text(splice.tree)

##3割のデータでテスト
splice.prd <- predict(splice.tree, newdata=splice.All[-lrn,], type="class")
##予測と本当のラベル
(tbl <- table(splice.prd,(splice.All[-lrn,])$label))

##正解率は、93.7%。
##「さっきの正解率(学習データが最大どこまで予測できるか)」は、99%。
sum(diag(tbl))/sum(tbl)
##つまり、未学習のデータを持ってきたときに予測率が悪くなる。

##ちょっとだけ、[-lrn]の使い方。
seq(10,100,by=10)
##5番目、3番目、8番目を選ぶ。
(seq(10,100,by=10))[c(5,3,8)]
##5番目、3番目、8番目を含まない物全部。
(seq(10,100,by=10))[-c(5,3,8)]
##つまり、splice.All[lrn,]番目を学習に用いて、splice.All[-lrn,]を検査に使う

##植木屋さんにお任せ。Trivialな枝から、切っていって、
##未学習データに対するエラーがどの様に推移するか。
splice.tree <- rpart(label ~ ., data=splice.All,
                     parms = list(split = "information"),
                     control = rpart.control(minsplit=0, cp=0,
                       maxdepth = 30))

plotcp(splice.tree)
printcp(splice.tree)

##12の辺りのCP(複雑度)は、0.00294118なので、、、
##            CP nsplit  rel error  xerror      xstd
## 1  0.22908497      0 1.00000000 1.00000 0.0184181
## 2  0.15490196      2 0.54183007 0.54183 0.0161816
## .....
## 8  0.00294118     12 0.10980392 0.11111 0.0082910

splice.tree <- rpart(label ~ ., data=splice.All,
                     parms = list(split = "information"),
                     control = rpart.control(minsplit=0, cp=0.00294118,
                       maxdepth = 30))

plot(splice.tree, branch=0.9, margin=0.04);text(splice.tree)
##ってのが、よさそう。

splice.tree
## 16) s13=G 537   24 EI (0.955307263 0.040968343 0.003724395) *
## 24) s9=G 150   13 EI (0.913333333 0.000000000 0.086666667) *
## 34) s11=A 90   13 EI (0.855555556 0.088888889 0.055555556) *
##   
## 18) s7=A 154   13 IE (0.000000000 0.915584416 0.084415584) *
## 20) s6=C,T 543   45 IE (0.003683241 0.917127072 0.079189687) *
## 70) s6=C,T 98   22 IE (0.183673469 0.775510204 0.040816327) *
##   
##  7) s13=A,C,T 942    5 N (0.003184713 0.002123142 0.994692144) *
## 11) s7=C,G,T 244    5 N (0.012295082 0.008196721 0.979508197) *
## 13) s10=A,C,G 208    1 N (0.004807692 0.000000000 0.995192308) *
## 19) s7=C,G,T 48    0 N (0.000000000 0.000000000 1.000000000) *
## 21) s6=A,G 79   16 N (0.012658228 0.189873418 0.797468354) *
## 25) s9=A,C,T 63    0 N (0.000000000 0.000000000 1.000000000) *
## 71) s6=A,G 25   11 N (0.360000000 0.080000000 0.560000000) *

path.rpart(splice.tree, node = c(16,24,34))
##EIと分類されるルールは、、、
##  node number: 16   node number: 24   node number: 34  
##    s8=G              s8=A,C,T          s8=G           
##    s9=G              s9=G              s9=G           
##    s10=T             s10=T             s10=T          
##                                        s11=A          
##    s13=G             s13=G             s13=A,C,T 

path.rpart(splice.tree, node = c(18,20,70))
##IEと分類されるルールは、、、
##  node number: 18  node number: 20  node number: 70 
##                       s6=C,T           s6=C,T   
##    s7=A               s7=A                      
##    s8=G               s8=G          　 s8=G     
##    s9=A,C,T                            s9=G     
##    s10=T              s10=A,C,G        s10=T    
##                                        s11=C,G,T
##                                        s13=A,C,T