RStudio を利用して R による機械学習 (分類) を行う

このファイルは R Markdown で作成されてます
今日行った R でのデータ解析について、講義の終わりに R Markdown でレポートを作成してみましょう

探索的データ解析　　

データの視覚化や相関性の解析など、予測モデルの構築に必要な探索的データ解析を行います
Rstudio を開いてデータを読み込んでみましょう

1 データの読み込み

setwd("/Users/ohmoriohmori/Desktop/アグリバイオ/講義関連/2019年度/ionome/")
Con <- read.csv(("C1C2.csv"), row.names=1) # row.names=1 は 1列目を行の名前にするという意味
Dro <- read.csv(("D1D2.csv"), row.names=1)
Y <- rbind(Con, Dro) # Con と Dro の行を束ねる
dim(Y) # Yの行列の大きさを確認

## [1] 384  24

head(Y) # Yの先頭部分 6 行を表示

tail(Y) # Yの末尾部分 6 行を表示

Y[1:5,1:5] # Yの[1から5行, 1から5列] を表示

2 データを視る

2-1 オブジェクトの内容を情報付きで表示する

str(Y) # 読み込まれたデータ型の確認

## 'data.frame':    384 obs. of  24 variables:
##  $ BlockID: Factor w/ 2 levels "C","D": 1 1 1 1 1 1 1 1 1 1 ...
##  $ PlotID : int  301 301 301 301 302 302 302 302 303 303 ...
##  $ IndID  : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ As     : num  224 282 173 264 372 ...
##  $ Br     : num  0.683 0.293 0.451 0.594 0.115 ...
##  $ Ca     : num  1001242 773615 833767 1172171 437903 ...
##  $ Cd     : num  0.232 0 0.279 0.293 0 ...
##  $ Cl     : num  3.02 2.41 3.78 1.94 1.74 ...
##  $ Co     : num  0.236 0.301 0.446 0.207 0.456 ...
##  $ Cs     : num  0.00614 0.00646 0.01581 0.01306 0.02404 ...
##  $ Cu     : num  0.185 0.162 0.289 0.21 0.338 ...
##  $ Fe     : num  1.79 2.4 2.46 1.61 2.11 ...
##  $ K      : num  15.5 13.2 21.1 13.7 41.4 ...
##  $ Mn     : num  2.19 2.37 1.79 2.56 0.78 ...
##  $ Mo     : num  0.0408 0.0311 0.0206 0.0393 0.026 ...
##  $ Ni     : num  0.372 0.341 0.563 0.362 0.662 ...
##  $ P      : num  8156 5242 2269 2691 7270 ...
##  $ Rb     : num  0.0633 0 0.0147 0 0.0467 ...
##  $ S      : num  25185 11879 8546 10291 6761 ...
##  $ Si     : num  2179 7687 6359 11248 867 ...
##  $ Sr     : num  1.07 0.676 1.139 0.839 0.422 ...
##  $ Zn     : num  0.309 0.162 0.204 0.175 0.251 ...
##  $ name   : Factor w/ 4 levels "AKASAYA","Houjaku Kuwazu",..: 4 4 4 4 1 1 1 1 2 2 ...
##  $ AFW    : num  47 54.7 36.2 43.8 68.8 ...

2-2 基本統計量を知る

summary(Y$Ca) # summary 関数で Y の Ca 列のデータを要約

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    8459  322861  439444  490948  662835 1198897      13

summary(Y[Y$BlockID=="C", "Ca"]) # Y の[ Y の BlockID 列が C の行, Ca の列] を指定して要約

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    8459  461296  650713  612376  767536 1172171       3

summary(Y[Y$BlockID=="D", "Ca"]) # 同様に、BlockID 列が D の行の Ca 列を要約

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  121923  272418  347136  364850  428560 1198897      10

2-3 boxplot でデータを視る

Control のみの箱ヒゲ図

par(family = "HiraKakuProN-W3") # 日本語表示のためにフォントを指定
boxplot(Y[Y$BlockID=="C", "Ca"], # Y の[ Y の BlockID 列が C の行, Ca の列]を対象に指定
        main = "Control の葉中 Ca", # 図のタイトルを指定
        col=2, # 図に使う色を指定、番号は 1 から順に「黒，赤，緑，青，水色，紫，黄，灰」となっている
        varwidth=TRUE, # 箱の幅を各項目のサンプルサイズの平方根に比例させて表示
        notch=TRUE, # ノッチ(切れ目)で各群の中央値の95%信頼区間を表示
        outline=TRUE) # 外れ値を表示

Dorought のみの箱ヒゲ図

par(family = "HiraKakuProN-W3")
boxplot(Y[Y$BlockID=="D", "Ca"], 
        main = "Dorought の葉中 Ca", 
        col=3, 
        varwidth=TRUE, 
        notch=TRUE, 
        outline=TRUE)

Control と Dorought の箱ひげ図による比較

library(beeswarm)
dim(Y)

## [1] 384  24

Y$BlockID # Y の BlockID 列を表示して、情報を確認

##   [1] C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
##  [36] C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
##  [71] C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
## [106] C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
## [141] C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
## [176] C C C C C C C C C C C C C C C C C D D D D D D D D D D D D D D D D D D
## [211] D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D
## [246] D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D
## [281] D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D
## [316] D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D
## [351] D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D
## Levels: C D

#
par(family = "HiraKakuProN-W3")
boxplot(Ca ~ BlockID, data = Y, # Y のデータで、Ca を BlockID をグループにして箱ひげ図を描く
        col = 2:3, # C と D にそれぞれ色を指定 col = c(2,3) と書いてもいい
        main = "Control と Dorought の葉中 Ca の比較",
        varwidth=TRUE, 
        notch=TRUE,
        outline=TRUE
        )
# Box plot に beaswarm を重ね書きして、データの分布を確認
beeswarm(Ca ~ BlockID, data = Y,
         col = "orange", 
         method = "square", # 点の配置方法を指定、他には swarm, center, hex, など 
         pch = 20, # プロット点の形を指定
         cex = 0.5, # 点の大きさを指定、標準は 1
         add = TRUE) # TRUE で図を重ね描きする

for文で作業を繰り返し、全元素についての箱ひげ図を作る

head(Y)

# 箱ひげ図を描きたい列の名前を、traitlist に代入、c() はベクトルを作成する基本的な関数
traitlist <-  c("As", "Br" , "Ca", "Cd", "Cl",  "Co" , "Cs" , "Cu", "Fe", "K", "Mn" ,"Mo", "Ni" ,"P", "Rb", "S", "Si" ,"Sr", "Zn", "AFW")
traitlist # 間違いなく代入されたか確認

##  [1] "As"  "Br"  "Ca"  "Cd"  "Cl"  "Co"  "Cs"  "Cu"  "Fe"  "K"   "Mn" 
## [12] "Mo"  "Ni"  "P"   "Rb"  "S"   "Si"  "Sr"  "Zn"  "AFW"

# 作画の準備
par(family = "HiraKakuProN-W3")
par(mfrow=c(2,5)) # 画面を 2 行 5 列に分割する、1 画面に戻す場合は mfrow=c(1,1) と指定する
# for 文で traitlist にある列の名前について箱ひげ図の作画を繰り返す
  for(i in traitlist){
    boxplot(Y[,i] ~ Y$BlockID,
            col = 2:3, 
            main = paste0("葉中", i, "の比較"), # paste0()で文字の結合
            varwidth=TRUE, 
            notch=TRUE,
            outline=TRUE)
    beeswarm(Y[,i] ~ Y$BlockID,
             col = "orange",
             method = "square", 
             pch = 20, 
             cex=0.4, 
             add = TRUE)
  }

2-4 ヒストグラムでデータを視る

Control と Dorought を分けずにヒストグラム作成

hist(na.omit(Y$P), # na.omit() で NA を含む行を削除、Y の P 列を対象に指定
     freq = FALSE, # FALSE の場合確率密度、TRUE の場合頻度を表示
     breaks = "freedman-diaconis", # 階級数と階級幅を指定、"freedman-diaconis"は適切な階級幅を計算するアルゴリズム
     xlab = "P", # X軸ラベルを指定
     ylab = "Density", # Y軸ラベルを指定
     main = "All Histgram",
     col = "#ff00ff40") # HTMLカラーコードで色を指定 (https://www.colordic.org) ff00ff は magenta、後ろの 40 は透過度
#
lines(density(na.omit(Y$P)), # 密度曲線を追加 
      col = "orange", 
      lwd = 2) # 線の太さを指定
#
rug(na.omit(Y$P), # ラグプロットを追加
    col = "#ff00ff40") 
#
abline(v = median(na.omit(Y$P)), # 中央値を表示
       lty = 1, # 線の種類を指定
       col = 2, # 中央値を赤色で表示
       lwd = 1) 
#
abline(v = mean(na.omit(Y$P)), # 平均値を表示
       lty = 1,
       col = 4, # 平均値を青色で表示
       lwd = 1)

Control のみのヒストグラム

hist(na.omit(Y[Y$BlockID=="C", "P"]), 
     freq = FALSE,
     breaks = "freedman-diaconis",
     xlab = "P", 
     ylab = "Density",
     main = "Control Histgram",
     col = "#FF000040")
#
lines(density(na.omit(Y[Y$BlockID=="C", "P"])), 
      col = "orange", 
      lwd = 2)
#
rug(na.omit(Y[Y$BlockID=="C", "P"]),
    col = "#FF000040")
#
abline(v = median(na.omit(Y[Y$BlockID=="C", "P"])), #中央値を表示
       lty = 1, 
       col = 2, 
       lwd = 1)

Dorought のみのヒストグラム

hist(na.omit(Y[Y$BlockID=="D", "P"]), 
     freq = FALSE,
     breaks = "freedman-diaconis",
     xlab = "P", 
     ylab = "Density", 
     main = "Dorought Histgram",
     col = "#0000ff40")
#
lines(density(na.omit(Y[Y$BlockID=="D", "P"])), 
      col = "orange", 
      lwd = 2)
#
rug(na.omit(Y[Y$BlockID=="D", "P"]), 
    col = "#0000ff40")
#
abline(v = median(na.omit(Y[Y$BlockID=="D", "P"])), #中央値を表示
       lty = 1, 
       col = 2, 
       lwd = 1)

Control と Dorought のヒストグラムを重ね書き

hist(na.omit(Y[Y$BlockID=="C", "P"]), 
     xlim = c(0, 23000), # データの要約を見て、X軸の範囲を指定
     freq = FALSE, 
     breaks = seq(0,23000,800), # 0 ~ 23000 の範囲を 800 ずつ分ける
     col = "#FF000040", 
     border = "#ff00ff", #棒の枠線の色を指定
     xlab = "P", 
     ylab = "Density",
     main = "Histgram")
hist(na.omit(Y[Y$BlockID=="D", "P"]), 
     xlim = c(0, 23000), # 図を重ねるので上にあわせる
     freq = FALSE, 
     breaks = seq(0,23000,800), # 図を重ねるので上にあわせる
     col = "#0000ff40", 
     border = "#0000ff", 
     add = TRUE)
lines(density(na.omit(Y[Y$BlockID=="C", "P"])), 
      col = 2, 
      lwd = 1)
lines(density(na.omit(Y[Y$BlockID=="D", "P"])), 
      col = 4, 
      lwd = 1)
abline(v = median(na.omit(Y[Y$BlockID=="C", "P"])), 
       lty = 1,
       col = 2,
       lwd = 1) 
abline(v = median(na.omit(Y[Y$BlockID=="D", "P"])), 
       lty = 1,
       col = 4,
       lwd = 1) 
rug(na.omit(Y[Y$BlockID=="C", "P"]), 
    col = "#FF000040", 
    lwd = 2)
rug(na.omit(Y[Y$BlockID=="D", "P"]), 
    col = "#0000ff40", 
    lwd = 2)
#図に凡例を加える
legend("topright", #凡例の位置を指定
       pch = c(16,16),
       col = c("#FF000080", "#0000ff80"), 
       legend = c("Control","Drought"), #凡例の文字を指定
       text.col = c("#FF000080", "#0000ff80"), #凡例の文字の色を指定
       cex = 2, #凡例文字の大きさ
       pt.cex = 1, #マークの大きさ
       bty = 'n', #凡例の枠のスタイル、"o" だと四方を枠線で囲う
       inset = c(0.15, -0.01), # インセットの位置を微調整 c(横位置, 縦位置)
       x.intersp=0.3, # インセットの間隔を微調整
       y.intersp=0.7)

2-5 データの全体を視る

tabplot で変数をソートしてグラフ化し、全体を眺める

library(tabplot)

## Loading required package: bit

## Attaching package bit

## package:bit (c) 2008-2012 Jens Oehlschlaegel (GPL-2)

## creators: bit bitwhich

## coercion: as.logical as.integer as.bit as.bitwhich which

## operator: ! & | xor != ==

## querying: print length any all min max range sum summary

## bit access: length<- [ [<- [[ [[<-

## for more help type ?bit

## 
## Attaching package: 'bit'

## The following object is masked from 'package:base':
## 
##     xor

## Loading required package: ff

## Attaching package ff

## - getOption("fftempdir")=="/var/folders/7z/0_f2n_ts72g995z53qk2l3hm0000gn/T//RtmpbA46h0"

## - getOption("ffextension")=="ff"

## - getOption("ffdrop")==TRUE

## - getOption("fffinonexit")==TRUE

## - getOption("ffpagesize")==65536

## - getOption("ffcaching")=="mmnoflush"  -- consider "ffeachflush" if your system stalls on large writes

## - getOption("ffbatchbytes")==16777216 -- consider a different value for tuning your system

## - getOption("ffmaxbytes")==536870912 -- consider a different value for tuning your system

## 
## Attaching package: 'ff'

## The following objects are masked from 'package:bit':
## 
##     clone, clone.default, clone.list

## The following objects are masked from 'package:utils':
## 
##     write.csv, write.csv2

## The following objects are masked from 'package:base':
## 
##     is.factor, is.ordered

## Loading required package: ffbase

## 
## Attaching package: 'ffbase'

## The following objects are masked from 'package:ff':
## 
##     [.ff, [.ffdf, [<-.ff, [<-.ffdf

## The following objects are masked from 'package:base':
## 
##     %in%, table

#
head(Y) # データの確認

# tabplot に使いたい列の名前を、trait に代入、c() はベクトルを作成する基本的な関数
trait <-  c("As", "Br" , "Ca", "Cd", "Cl",  "Co" , "Cs" , "Cu", "Fe", "K", "Mn" ,"Mo", "Ni" ,"P", "Rb", "S", "Si" ,"Sr", "Zn", "AFW","BlockID")
trait # 間違いなく代入されたか確認

##  [1] "As"      "Br"      "Ca"      "Cd"      "Cl"      "Co"      "Cs"     
##  [8] "Cu"      "Fe"      "K"       "Mn"      "Mo"      "Ni"      "P"      
## [15] "Rb"      "S"       "Si"      "Sr"      "Zn"      "AFW"     "BlockID"

#
tableplot (dat = Y[,c(trait)], 
           from = 0, # データの範囲 (下限)、0 ~ 100 の範囲で指定
           to = 100, # データの範囲 (上限)、0 ~ 100 の範囲で指定
           #select = c(Ca, Sr, K, Rb, BlockID), # 視覚化する変数を指定できる、何も指定しなければの全変数が対象
           sortCol = "Ca", # Ca 列でソートする
           nBins = 50, # 分割数を指定、 100 だと 100 分割 (1%区切)
           colorNA_num = "red", # NA データの色を指定
           scales = "lin" #スケールの指定、lin, log もしくは auto を選ぶ
           #, subset_string = "BlockID == 'C'" #フィルターを指定することもできる
           )

2-6 データの相関を視る

元素間の (総当たりの) 相関係数 (r) を計算する

# 相関係数の計算に使いたい列の名前を、trait に代入
trait <-  c("As", "Br" , "Ca", "Cd", "Cl",  "Co" , "Cs" , "Cu", "Fe", "K", "Mn" ,"Mo", "Ni" ,"P", "Rb", "S", "Si" ,"Sr", "Zn", "AFW")
trait # 間違いなく代入されたか確認

##  [1] "As"  "Br"  "Ca"  "Cd"  "Cl"  "Co"  "Cs"  "Cu"  "Fe"  "K"   "Mn" 
## [12] "Mo"  "Ni"  "P"   "Rb"  "S"   "Si"  "Sr"  "Zn"  "AFW"

#
cor(Y[,c(trait)], # cor()関数で、相関係数の計算 
    use="pairwise.complete.obs", # 二変数の組合せごとに欠損値を含むサンプルを除く方法で計算
    method="pearson") # pearsonの積率相関係数を計算

##               As           Br          Ca           Cd           Cl
## As   1.000000000  0.079055846  0.31028397 -0.085291424 -0.004461164
## Br   0.079055846  1.000000000  0.28563117  0.182568163  0.608964088
## Ca   0.310283974  0.285631167  1.00000000 -0.091959969 -0.079284891
## Cd  -0.085291424  0.182568163 -0.09195997  1.000000000  0.208301016
## Cl  -0.004461164  0.608964088 -0.07928489  0.208301016  1.000000000
## Co  -0.140466236  0.077739926 -0.36014110  0.284941233  0.223468911
## Cs  -0.070842857 -0.006479914 -0.13887026  0.056394069  0.061161902
## Cu  -0.258349530  0.151716467 -0.37937879  0.320919156  0.206771165
## Fe  -0.045472279  0.130252625 -0.28559458  0.252928616  0.310674034
## K   -0.253804764 -0.031555635 -0.60613800  0.236611908  0.359272079
## Mn   0.033266591  0.353104077  0.08407369  0.126698245  0.311433649
## Mo   0.047202485 -0.180226339  0.15791573 -0.151197935 -0.321776332
## Ni  -0.280951144  0.018625018 -0.56475596  0.305277180  0.215409001
## P   -0.057246295 -0.225625538 -0.13962123 -0.012615317 -0.376270468
## Rb  -0.221602275 -0.429805266 -0.42965734 -0.001295872 -0.323585323
## S    0.092916858  0.064490856  0.26234864  0.023603569 -0.145339030
## Si   0.165108845  0.181073404  0.43096392 -0.015870763 -0.025643421
## Sr   0.220462538 -0.068831805  0.74035922 -0.150042982 -0.255373050
## Zn  -0.257828344  0.022977478 -0.35930212  0.109664910 -0.033310137
## AFW  0.193867395 -0.360205380  0.39136474 -0.259844797 -0.443291866
##              Co           Cs         Cu          Fe           K
## As  -0.14046624 -0.070842857 -0.2583495 -0.04547228 -0.25380476
## Br   0.07773993 -0.006479914  0.1517165  0.13025262 -0.03155563
## Ca  -0.36014110 -0.138870256 -0.3793788 -0.28559458 -0.60613800
## Cd   0.28494123  0.056394069  0.3209192  0.25292862  0.23661191
## Cl   0.22346891  0.061161902  0.2067712  0.31067403  0.35927208
## Co   1.00000000  0.334750706  0.5379913  0.59242129  0.67939617
## Cs   0.33475071  1.000000000  0.1788428  0.13527170  0.27644699
## Cu   0.53799132  0.178842764  1.0000000  0.38185434  0.44714867
## Fe   0.59242129  0.135271704  0.3818543  1.00000000  0.43116818
## K    0.67939617  0.276446995  0.4471487  0.43116818  1.00000000
## Mn   0.25983741  0.076306737  0.1220785  0.23380626  0.16352330
## Mo  -0.37012858 -0.138666068 -0.2564617 -0.29495047 -0.32794912
## Ni   0.72171882  0.175457987  0.6886155  0.47012005  0.71270756
## P    0.22580609  0.085526587  0.1856352 -0.04833730  0.19520343
## Rb   0.14145519 -0.002361797  0.1324541 -0.02825120  0.34492098
## S    0.09765117 -0.020442175  0.1339574  0.02925047 -0.01612103
## Si  -0.13558564 -0.106440986 -0.1567966 -0.05328718 -0.26761261
## Sr  -0.39870718 -0.075382067 -0.4759750 -0.36439970 -0.56804205
## Zn   0.26334753  0.079583374  0.2544933  0.06970183  0.26540395
## AFW -0.42187749 -0.030901697 -0.4206691 -0.45584923 -0.47989737
##              Mn           Mo          Ni           P           Rb
## As   0.03326659  0.047202485 -0.28095114 -0.05724630 -0.221602275
## Br   0.35310408 -0.180226339  0.01862502 -0.22562554 -0.429805266
## Ca   0.08407369  0.157915734 -0.56475596 -0.13962123 -0.429657336
## Cd   0.12669824 -0.151197935  0.30527718 -0.01261532 -0.001295872
## Cl   0.31143365 -0.321776332  0.21540900 -0.37627047 -0.323585323
## Co   0.25983741 -0.370128579  0.72171882  0.22580609  0.141455186
## Cs   0.07630674 -0.138666068  0.17545799  0.08552659 -0.002361797
## Cu   0.12207852 -0.256461695  0.68861553  0.18563520  0.132454123
## Fe   0.23380626 -0.294950467  0.47012005 -0.04833730 -0.028251205
## K    0.16352330 -0.327949117  0.71270756  0.19520343  0.344920979
## Mn   1.00000000 -0.213793917  0.31441734 -0.11110945 -0.215615947
## Mo  -0.21379392  1.000000000 -0.32876296  0.07671938  0.129341560
## Ni   0.31441734 -0.328762958  1.00000000  0.22335853  0.336191285
## P   -0.11110945  0.076719379  0.22335853  1.00000000  0.487504051
## Rb  -0.21561595  0.129341560  0.33619128  0.48750405  1.000000000
## S    0.02922530 -0.035973096 -0.01985881  0.45766989  0.015099686
## Si   0.20640562  0.001952701 -0.25538983 -0.09792299 -0.225426232
## Sr  -0.15166473  0.205009322 -0.66263109 -0.13804555 -0.281508663
## Zn   0.15131201 -0.084045159  0.36164971  0.28229732  0.298881592
## AFW -0.28422885  0.314748873 -0.65500728  0.13319383 -0.062386169
##               S           Si          Sr          Zn         AFW
## As   0.09291686  0.165108845  0.22046254 -0.25782834  0.19386740
## Br   0.06449086  0.181073404 -0.06883180  0.02297748 -0.36020538
## Ca   0.26234864  0.430963921  0.74035922 -0.35930212  0.39136474
## Cd   0.02360357 -0.015870763 -0.15004298  0.10966491 -0.25984480
## Cl  -0.14533903 -0.025643421 -0.25537305 -0.03331014 -0.44329187
## Co   0.09765117 -0.135585641 -0.39870718  0.26334753 -0.42187749
## Cs  -0.02044218 -0.106440986 -0.07538207  0.07958337 -0.03090170
## Cu   0.13395738 -0.156796638 -0.47597504  0.25449326 -0.42066909
## Fe   0.02925047 -0.053287176 -0.36439970  0.06970183 -0.45584923
## K   -0.01612103 -0.267612611 -0.56804205  0.26540395 -0.47989737
## Mn   0.02922530  0.206405621 -0.15166473  0.15131201 -0.28422885
## Mo  -0.03597310  0.001952701  0.20500932 -0.08404516  0.31474887
## Ni  -0.01985881 -0.255389830 -0.66263109  0.36164971 -0.65500728
## P    0.45766989 -0.097922991 -0.13804555  0.28229732  0.13319383
## Rb   0.01509969 -0.225426232 -0.28150866  0.29888159 -0.06238617
## S    1.00000000  0.081203017  0.10364855  0.03217804  0.19829644
## Si   0.08120302  1.000000000  0.36301471 -0.14116449  0.19310561
## Sr   0.10364855  0.363014710  1.00000000 -0.33980912  0.60798874
## Zn   0.03217804 -0.141164495 -0.33980912  1.00000000 -0.29779054
## AFW  0.19829644  0.193105610  0.60798874 -0.29779054  1.00000000

#
cor_res <- cor(Y[,c(trait)], # 結果を cor_res に代入
               use="pairwise.complete.obs", 
               method="pearson") 
#
write.csv(cor_res,"corrilation.csv") #cor_res を"corrilation.csv"としてワーキングディレクトリに保存

視覚的に相関係数を示す

library(corrplot)

## corrplot 0.84 loaded

corrplot(cor(Y[,c(trait)], 
            use="pairwise.complete.obs", 
            method="pearson"))

C、D の条件ごとに相関係数 (r) を計算する

control での元素間の相関係数

cor(Y[Y$BlockID=="C",c(trait)], 
    use="pairwise.complete.obs", 
    method="pearson")

##              As          Br          Ca          Cd          Cl
## As   1.00000000  0.15712576  0.17508633 -0.09317700  0.14292614
## Br   0.15712576  1.00000000  0.60478305  0.12654990  0.69723462
## Ca   0.17508633  0.60478305  1.00000000  0.04312002  0.28187640
## Cd  -0.09317700  0.12654990  0.04312002  1.00000000  0.17552335
## Cl   0.14292614  0.69723462  0.28187640  0.17552335  1.00000000
## Co  -0.02906399  0.01099533 -0.07824669  0.21964733  0.21343595
## Cs  -0.04686634 -0.05482708 -0.17046613  0.17106087  0.15620169
## Cu  -0.18831684  0.02995971 -0.14741992  0.27202555  0.14832132
## Fe   0.02797170  0.03573507 -0.14193811  0.19000250  0.29847144
## K   -0.16344175 -0.25979413 -0.45608155  0.12551178 -0.01855352
## Mn   0.18530015  0.60470657  0.36418839  0.08786325  0.48824725
## Mo  -0.03833350 -0.21658091 -0.09214804 -0.17809119 -0.41685069
## Ni  -0.09897902 -0.04736216 -0.17251391  0.28547868  0.06745389
## P    0.04464147 -0.04736569 -0.08430054  0.03769769 -0.07333995
## Rb  -0.19366076 -0.37882998 -0.40378128  0.01345047 -0.39604440
## S    0.09055270  0.11461631  0.26936129  0.06082767 -0.03609165
## Si   0.09407283  0.33861448  0.33042895  0.05761119  0.19406888
## Sr  -0.01296876  0.21257648  0.57014143  0.06544653  0.07533913
## Zn  -0.18053749 -0.07040990 -0.20245548  0.14445964 -0.09838000
## AFW -0.03778921 -0.41545378 -0.31257789 -0.19065198 -0.22100808
##               Co          Cs           Cu           Fe           K
## As  -0.029063994 -0.04686634 -0.188316838  0.027971703 -0.16344175
## Br   0.010995334 -0.05482708  0.029959710  0.035735073 -0.25979413
## Ca  -0.078246689 -0.17046613 -0.147419919 -0.141938106 -0.45608155
## Cd   0.219647329  0.17106087  0.272025551  0.190002504  0.12551178
## Cl   0.213435952  0.15620169  0.148321324  0.298471445 -0.01855352
## Co   1.000000000  0.36556623  0.450920527  0.713635042  0.60424601
## Cs   0.365566228  1.00000000  0.244222213  0.322152860  0.30870665
## Cu   0.450920527  0.24422221  1.000000000  0.342245143  0.31326712
## Fe   0.713635042  0.32215286  0.342245143  1.000000000  0.45869582
## K    0.604246009  0.30870665  0.313267122  0.458695822  1.00000000
## Mn   0.185434607  0.02807283 -0.003910801  0.335894180 -0.04888500
## Mo  -0.345918125 -0.10746925 -0.207239006 -0.373309835 -0.19273395
## Ni   0.800616131  0.29928899  0.648759282  0.563219628  0.67576844
## P    0.393565742  0.06791941  0.253115902  0.252147843  0.51900753
## Rb   0.114434760 -0.02516006  0.122510217  0.002124609  0.50798303
## S    0.319723842 -0.03197932  0.211083275  0.283916248  0.24980876
## Si   0.007712699 -0.12221613 -0.000987608  0.057289991 -0.15818743
## Sr  -0.065747829 -0.07955780 -0.168336221 -0.183280263 -0.38105729
## Zn   0.111191727  0.02560148  0.025585739  0.035345128  0.34891722
## AFW -0.088565957  0.01357191  0.050499193 -0.109881632  0.01028965
##               Mn          Mo          Ni             P           Rb
## As   0.185300152 -0.03833350 -0.09897902  0.0446414720 -0.193660760
## Br   0.604706573 -0.21658091 -0.04736216 -0.0473656857 -0.378829982
## Ca   0.364188386 -0.09214804 -0.17251391 -0.0843005421 -0.403781284
## Cd   0.087863251 -0.17809119  0.28547868  0.0376976909  0.013450474
## Cl   0.488247251 -0.41685069  0.06745389 -0.0733399463 -0.396044400
## Co   0.185434607 -0.34591812  0.80061613  0.3935657424  0.114434760
## Cs   0.028072833 -0.10746925  0.29928899  0.0679194072 -0.025160058
## Cu  -0.003910801 -0.20723901  0.64875928  0.2531159022  0.122510217
## Fe   0.335894180 -0.37330984  0.56321963  0.2521478434  0.002124609
## K   -0.048884996 -0.19273395  0.67576844  0.5190075254  0.507983026
## Mn   1.000000000 -0.27012893  0.08472920  0.0095712213 -0.363583892
## Mo  -0.270128926  1.00000000 -0.25435374 -0.0322532594  0.209937541
## Ni   0.084729203 -0.25435374  1.00000000  0.5585685444  0.286068627
## P    0.009571221 -0.03225326  0.55856854  1.0000000000  0.444347063
## Rb  -0.363583892  0.20993754  0.28606863  0.4443470633  1.000000000
## S    0.182841486 -0.12868508  0.37982053  0.5560949923  0.117628331
## Si   0.400862131 -0.11994754 -0.03771163 -0.0987686432 -0.188951209
## Sr  -0.034029299 -0.05282475 -0.21155788 -0.2806124736 -0.246019428
## Zn   0.021730187  0.02987598  0.13284933  0.1772627344  0.350167675
## AFW -0.255522361  0.16286266 -0.02664745 -0.0002624772  0.100516487
##               S           Si          Sr          Zn           AFW
## As   0.09055270  0.094072831 -0.01296876 -0.18053749 -0.0377892065
## Br   0.11461631  0.338614478  0.21257648 -0.07040990 -0.4154537833
## Ca   0.26936129  0.330428954  0.57014143 -0.20245548 -0.3125778930
## Cd   0.06082767  0.057611191  0.06544653  0.14445964 -0.1906519800
## Cl  -0.03609165  0.194068882  0.07533913 -0.09838000 -0.2210080758
## Co   0.31972384  0.007712699 -0.06574783  0.11119173 -0.0885659575
## Cs  -0.03197932 -0.122216129 -0.07955780  0.02560148  0.0135719063
## Cu   0.21108328 -0.000987608 -0.16833622  0.02558574  0.0504991930
## Fe   0.28391625  0.057289991 -0.18328026  0.03534513 -0.1098816323
## K    0.24980876 -0.158187425 -0.38105729  0.34891722  0.0102896521
## Mn   0.18284149  0.400862131 -0.03402930  0.02173019 -0.2555223605
## Mo  -0.12868508 -0.119947544 -0.05282475  0.02987598  0.1628626613
## Ni   0.37982053 -0.037711625 -0.21155788  0.13284933 -0.0266474505
## P    0.55609499 -0.098768643 -0.28061247  0.17726273 -0.0002624772
## Rb   0.11762833 -0.188951209 -0.24601943  0.35016767  0.1005164871
## S    1.00000000  0.027898071 -0.08126895  0.04678525 -0.0116477156
## Si   0.02789807  1.000000000  0.21078896 -0.01147076 -0.1475816355
## Sr  -0.08126895  0.210788957  1.00000000 -0.06862429 -0.2128820049
## Zn   0.04678525 -0.011470760 -0.06862429  1.00000000 -0.1085389946
## AFW -0.01164772 -0.147581636 -0.21288200 -0.10853899  1.0000000000

Drought での元素間の相関係数

cor(Y[Y$BlockID=="D",c(trait)], 
    use="pairwise.complete.obs", 
    method="pearson")

##               As          Br          Ca           Cd          Cl
## As   1.000000000  0.15466928  0.29716566  0.047541013  0.15915035
## Br   0.154669280  1.00000000  0.47939313  0.130304138  0.62944458
## Ca   0.297165656  0.47939313  1.00000000  0.062395284  0.30958848
## Cd   0.047541013  0.13030414  0.06239528  1.000000000  0.11871006
## Cl   0.159150350  0.62944458  0.30958848  0.118710058  1.00000000
## Co  -0.038895902 -0.11999487 -0.25418612  0.193937869 -0.03594459
## Cs  -0.084911181  0.01901011 -0.08520072 -0.073666160  0.02030309
## Cu  -0.124931916  0.01402888 -0.11982053  0.213052884 -0.08350576
## Fe   0.143917349 -0.01192630  0.08482170  0.150962025  0.10064763
## K   -0.134824375 -0.17765526 -0.43218629  0.138485561  0.22178449
## Mn   0.009875022  0.07722788  0.19386327  0.063084676  0.19722366
## Mo  -0.005130894  0.01135691  0.15173384  0.005802612 -0.22070574
## Ni  -0.220851290 -0.44350469 -0.40740622  0.144353421 -0.27087585
## P   -0.226420991 -0.32648452 -0.43570324 -0.005651794 -0.46762918
## Rb  -0.211731020 -0.59084931 -0.51118542 -0.073169462 -0.52626107
## S   -0.053191246  0.15400266 -0.02400238  0.105884012 -0.08631745
## Si   0.187095641  0.27664979  0.52690086  0.102280473  0.28432269
## Sr   0.278725845  0.27425915  0.69357309  0.027309000  0.35674213
## Zn  -0.220385391 -0.07483540 -0.26176649 -0.043462302 -0.26696066
## AFW -0.120742032 -0.01885354 -0.25410780  0.006254089 -0.29339384
##              Co          Cs            Cu           Fe           K
## As  -0.03889590 -0.08491118 -1.249319e-01  0.143917349 -0.13482438
## Br  -0.11999487  0.01901011  1.402888e-02 -0.011926303 -0.17765526
## Ca  -0.25418612 -0.08520072 -1.198205e-01  0.084821699 -0.43218629
## Cd   0.19393787 -0.07366616  2.130529e-01  0.150962025  0.13848556
## Cl  -0.03594459  0.02030309 -8.350576e-02  0.100647627  0.22178449
## Co   1.00000000  0.33744274  3.401234e-01  0.302111246  0.55755916
## Cs   0.33744274  1.00000000  1.091197e-01 -0.025318985  0.29423882
## Cu   0.34012339  0.10911972  1.000000e+00  0.078317010  0.16654039
## Fe   0.30211125 -0.02531898  7.831701e-02  1.000000000  0.08963233
## K    0.55755916  0.29423882  1.665404e-01  0.089632329  1.00000000
## Mn   0.16502617  0.10401461  7.699431e-05  0.021312199  0.10236204
## Mo  -0.18117848 -0.16408060 -3.940536e-02 -0.029327077 -0.23800645
## Ni   0.53345841  0.14636383  4.705195e-01  0.008287647  0.46699815
## P    0.24646948  0.11270239  2.988758e-01 -0.141558740  0.15252000
## Rb   0.07643369  0.01239774  4.004776e-02 -0.178473448  0.16778484
## S    0.07007703  0.02243498  4.118957e-01  0.026562767 -0.03465589
## Si  -0.06019284 -0.05722048 -6.788422e-02  0.315596447 -0.18283385
## Sr  -0.16412368 -0.03149428 -1.851668e-01  0.289691385 -0.15584776
## Zn   0.15371955  0.10690535  1.666320e-01 -0.187726152 -0.07012555
## AFW  0.16406038  0.06097175  1.931415e-01 -0.215222064 -0.04219111
##                Mn           Mo           Ni            P          Rb
## As   9.875022e-03 -0.005130894 -0.220851290 -0.226420991 -0.21173102
## Br   7.722788e-02  0.011356912 -0.443504686 -0.326484525 -0.59084931
## Ca   1.938633e-01  0.151733845 -0.407406220 -0.435703242 -0.51118542
## Cd   6.308468e-02  0.005802612  0.144353421 -0.005651794 -0.07316946
## Cl   1.972237e-01 -0.220705735 -0.270875847 -0.467629179 -0.52626107
## Co   1.650262e-01 -0.181178480  0.533458415  0.246469477  0.07643369
## Cs   1.040146e-01 -0.164080602  0.146363825  0.112702389  0.01239774
## Cu   7.699431e-05 -0.039405360  0.470519513  0.298875751  0.04004776
## Fe   2.131220e-02 -0.029327077  0.008287647 -0.141558740 -0.17847345
## K    1.023620e-01 -0.238006449  0.466998146  0.152520004  0.16778484
## Mn   1.000000e+00 -0.044263624  0.278621262 -0.149889133 -0.16452745
## Mo  -4.426362e-02  1.000000000 -0.114943801  0.128679653  0.12563303
## Ni   2.786213e-01 -0.114943801  1.000000000  0.385878724  0.44887194
## P   -1.498891e-01  0.128679653  0.385878724  1.000000000  0.56982013
## Rb  -1.645275e-01  0.125633025  0.448871943  0.569820131  1.00000000
## S   -1.817025e-02 -0.065220362  0.061737869  0.384419907 -0.07135124
## Si   2.739788e-01  0.004811533 -0.205955847 -0.376288454 -0.38749596
## Sr   2.520957e-01  0.045951217 -0.340159023 -0.584364831 -0.49314414
## Zn   1.097997e-01  0.002271636  0.201936646  0.430377112  0.22092269
## AFW  9.117986e-03  0.053987603  0.130675450  0.618049580  0.20504274
##               S           Si          Sr           Zn          AFW
## As  -0.05319125  0.187095641  0.27872584 -0.220385391 -0.120742032
## Br   0.15400266  0.276649793  0.27425915 -0.074835402 -0.018853540
## Ca  -0.02400238  0.526900860  0.69357309 -0.261766494 -0.254107800
## Cd   0.10588401  0.102280473  0.02730900 -0.043462302  0.006254089
## Cl  -0.08631745  0.284322693  0.35674213 -0.266960661 -0.293393839
## Co   0.07007703 -0.060192841 -0.16412368  0.153719551  0.164060378
## Cs   0.02243498 -0.057220478 -0.03149428  0.106905349  0.060971750
## Cu   0.41189566 -0.067884220 -0.18516683  0.166632038  0.193141456
## Fe   0.02656277  0.315596447  0.28969139 -0.187726152 -0.215222064
## K   -0.03465589 -0.182833849 -0.15584776 -0.070125547 -0.042191111
## Mn  -0.01817025  0.273978838  0.25209573  0.109799716  0.009117986
## Mo  -0.06522036  0.004811533  0.04595122  0.002271636  0.053987603
## Ni   0.06173787 -0.205955847 -0.34015902  0.201936646  0.130675450
## P    0.38441991 -0.376288454 -0.58436483  0.430377112  0.618049580
## Rb  -0.07135124 -0.387495963 -0.49314414  0.220922693  0.205042740
## S    1.00000000 -0.012977733 -0.15900218  0.200008674  0.228285254
## Si  -0.01297773  1.000000000  0.43395076 -0.209464779 -0.153321638
## Sr  -0.15900218  0.433950765  1.00000000 -0.321219456 -0.466225959
## Zn   0.20000867 -0.209464779 -0.32121946  1.000000000  0.402560515
## AFW  0.22828525 -0.153321638 -0.46622596  0.402560515  1.000000000

Control と Dorought それぞれの相関係数を図示する

corrplot(cor(Y[Y$BlockID=="C",c(trait)], 
             use="pairwise.complete.obs", 
             method="pearson") )

corrplot(cor(Y[Y$BlockID=="D",c(trait)], 
             use="pairwise.complete.obs", 
             method="pearson") )

相関係数を図中に表示する

library(corrplot)
corrplot.mixed(cor(Y[Y$BlockID=="C",c(trait)], 
                   use="pairwise.complete.obs", 
                   method="pearson"), 
               　　number.cex=0.5,
              　　 tl.cex=0.5)

元素間の相関を、1 対 1 の散布図を作成して視る

plot(Y[Y$BlockID=="C",c("Ca")], Y[Y$BlockID=="C",c("Sr")]) # 散布図の作成
#
(cor(Y[Y$BlockID=="C",c("Ca")] , 
     Y[Y$BlockID=="C",c("Sr")], 
     use="pairwise.complete.obs",  
     method="pearson")) #相関係数の計算と表示

## [1] 0.5701414

# sprintf()関数で、計算した相関係数の小数点以下２桁までを cor に代入
cor <- sprintf("%.2f",cor(Y[Y$BlockID=="C",c("Ca")] , 
                          Y[Y$BlockID=="C",c("Sr")], 
                          use="pairwise.complete.obs",  
                          method="pearson")) 
#図に凡例を加える
legend("bottomright", #凡例の位置を指定
       legend=paste0("r = ",cor),
       bty="n")

cor.text() 関数で相関があるかどうかを検定 (相関検定)

cor.test(Y[Y$BlockID=="C",c("Ca")], Y[Y$BlockID=="C",c("Sr")]) #無相関かどうかの検定

## 
##  Pearson's product-moment correlation
## 
## data:  Y[Y$BlockID == "C", c("Ca")] and Y[Y$BlockID == "C", c("Sr")]
## t = 9.4901, df = 187, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4652735 0.6592262
## sample estimates:
##       cor 
## 0.5701414

Control と Dorought での１対１の相関を同時に視る

plot(Y$Ca, Y$Sr, 
     col=c("red","green")[Y$BlockID])
summary(lm(Y[Y$BlockID=="C",c("Sr")]~ Y[Y$BlockID=="C",c("Ca")])) # lm () 関数による単回帰分析, BlockID=="C"

## 
## Call:
## lm(formula = Y[Y$BlockID == "C", c("Sr")] ~ Y[Y$BlockID == "C", 
##     c("Ca")])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4624 -0.1555 -0.0205  0.1420  0.5414 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  2.394e-01  4.996e-02   4.792 3.36e-06 ***
## Y[Y$BlockID == "C", c("Ca")] 7.314e-07  7.707e-08   9.490  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2252 on 187 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.3251, Adjusted R-squared:  0.3215 
## F-statistic: 90.06 on 1 and 187 DF,  p-value: < 2.2e-16

summary(lm(Y[Y$BlockID=="D",c("Sr")]~ Y[Y$BlockID=="D",c("Ca")])) # lm () 関数による単回帰分析, BlockID=="D"

## 
## Call:
## lm(formula = Y[Y$BlockID == "D", c("Sr")] ~ Y[Y$BlockID == "D", 
##     c("Ca")])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.23387 -0.05836 -0.01277  0.04173  0.27986 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  -8.442e-03  1.655e-02   -0.51    0.611    
## Y[Y$BlockID == "D", c("Ca")]  5.436e-07  4.208e-08   12.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08321 on 180 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.481,  Adjusted R-squared:  0.4782 
## F-statistic: 166.8 on 1 and 180 DF,  p-value: < 2.2e-16

abline(lm(Y[Y$BlockID=="C",c("Sr")]~ Y[Y$BlockID=="C",c("Ca")]), # 推定回帰直線を追加, BlockID=="C"
       col="red")
abline(lm(Y[Y$BlockID=="D",c("Sr")]~ Y[Y$BlockID=="D",c("Ca")]), # 推定回帰直線を追加, BlockID=="D"
       col="green")
cor1 <- sprintf("%.2f",cor(Y[Y$BlockID=="C",c("Ca")] ,
                           Y[Y$BlockID=="C",c("Sr")], 
                           use="pairwise.complete.obs", 
                           method="pearson")) 
cor2 <- sprintf("%.2f",cor(Y[Y$BlockID=="D",c("Ca")] , 
                           Y[Y$BlockID=="D",c("Sr")], 
                           use="pairwise.complete.obs", 
                           method="pearson")) 
legend("bottomright",
       legend = c(paste0("Control: r = ",cor1), 
                  paste0("Drought: r = ",cor2)), 
                  bty="n", 
                  text.col = c("red", "green"))

R で機械学習　

3 機械学習 (randomForest)

3-1 randomForest で分類

library(randomForest)

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

head(Y) # データの確認

#
Dependent  <- "BlockID" # C か D かを分類したいので Dependent に BlockID を代入 (目的変数、従属変数と呼ばれる)
Independent  <- c("As", "Br", "Ca", "Cd", "Cl", "Co", "Cs", "Cu", "Fe", "K", "Mn", "Mo", "Ni", "P", "Rb", "S", "Si", "Sr", "Zn") # 説明変数 (独立変数) にしたいデータ名を取り出す、今回は 19 元素全てを使って C か D かの分類モデルを作る
X2 <- na.omit(Y[,c(Dependent,Independent)]) # Yから必要な行・列だけのデータを選んで X2 に代入
head(X2) # 間違っていないか確認

#元素のデータで C か D かを予測する randomforest のモデルを作る
set.seed(1000)
ionome.rf.class <- randomForest(BlockID ~
                          As+Br+Ca+Cd+Cl+Co+Cs+Cu+Fe+K+Mn+Mo+Ni+P+Rb+S+Si+Sr+Zn,
                          data = X2,
                          proximity = TRUE, # 個々のデータの近接性を計算、MDS plot での結果の可視化に使用
                          importance = TRUE # 特徴量加工による重要度も計算・出力
                          ) 
#
ionome.rf.class # モデルの評価結果を表示

## 
## Call:
##  randomForest(formula = BlockID ~ As + Br + Ca + Cd + Cl + Co +      Cs + Cu + Fe + K + Mn + Mo + Ni + P + Rb + S + Si + Sr +      Zn, data = X2, proximity = TRUE, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 2.43%
## Confusion matrix:
##     C   D class.error
## C 184   5  0.02645503
## D   4 178  0.02197802

3-2 分類の結果を多次元尺度法 (MDS plot) で視る

proximityに基づいて、個体間の類似度を多次元尺度法 (MDS) で視覚化

MDS.result <- MDSplot(ionome.rf.class, #ランダムフォレストのモデル
                      X2$BlockID # ランダムフォレストで分類した因子
                      )

## Warning in RColorBrewer::brewer.pal(nlevs, "Set1"): minimal value for n is 3, returning requested palette with 3 different levels

誤分類されているものだけ色を変えて視る
正解データと判別結果を照合

TF <- ionome.rf.class$y == ionome.rf.class$predicted # 同じものが TRUE で、異なるものが FALSE で表示される
TF # 結果の確認

##   [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [12]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
##  [23]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [34]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [45]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [56]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [67]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [78]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [89]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [100]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [111]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [122]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [133]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [144]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [155]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [166]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [177]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
## [188]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [199]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [210]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [221]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [232]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [243]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [254]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [265]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [276]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [287]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [298]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE FALSE
## [309]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [320]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [331]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [342]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [353]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [364]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

正しく分類されているか否かの論理値（TRUE/FALSE）を数値（1/0）に変換

TF.n <- as.numeric(TF) # TF を数字型に変更して TF.n に代入
TF.n # 結果の確認

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1
##  [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [176] 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [211] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [246] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [281] 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 0 1 1 1 1 1 1 1
## [316] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [351] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

MDS plot に反映させる

plot(MDS.result$points, #MDS.result 内の point データを指定
     pch = as.numeric(X2$BlockID),  # マーカー指定の番号として、数値型にした X2 の BlockID 列 (1か2になっている) を指定
     col = TF.n + 2) # 色指定の番号として TF.n (0か1) に +2 した色を指定、+1でも+3でもいい (TFで色分けされることになる)
#誤分類されている番号を取得
n <- grep(FALSE, TF) # TFの中で FALSE のものの番号 (=元のサンプルデータの行数) を取得して、n に代入
#誤分類されているものの番号を表示
text(MDS.result$points[n, ] + 0.03, # [n, ]で誤分類されている番号の行のみ対象に指定、+ 0.03 でマーカーとずらしてテキストを表示
     as.character(n), # n に代入した番号を文字として利用
     cex = 1)

3-3 説明変数の重要度を視る

分類での説明変数の重要度（影響度）を表示

importance(ionome.rf.class)

##             C         D MeanDecreaseAccuracy MeanDecreaseGini
## As  0.3599485  2.966985            2.5556239         1.057988
## Br 12.7736285 14.778444           17.0057861         8.484295
## Ca 10.7062046 15.478205           17.5428014        18.054834
## Cd  6.4932883  1.088168            6.1424112         2.046616
## Cl 18.4822033 16.379894           22.3149927        11.294549
## Co  6.5735882  4.709777            7.9662462         3.585676
## Cs  2.7300523 -1.762537            0.4001371         1.659457
## Cu  6.1353206  7.798747            9.6228660         7.745024
## Fe 10.1906990  9.598195           13.4773391        12.472941
## K   8.6991423  6.102422           10.4476086        10.905235
## Mn  4.3711452  5.365930            6.6339661         2.270898
## Mo  2.1331008  2.119479            2.8565734         1.669175
## Ni 22.1595128 26.755922           33.5019392        39.733353
## P  12.5678506 11.952058           15.6330340         5.154011
## Rb  3.2110027  3.436116            4.7941015         1.636970
## S   7.9169847  7.874162           10.6948104         4.089192
## Si  4.1541032  7.178963            8.3175226         2.997612
## Sr 26.1361241 28.211966           35.2205644        45.594647
## Zn  2.0567738  6.500600            6.6807379         4.426212

importance <- as.data.frame(ionome.rf.class$importance)

重要度を図示

#特徴量加工による重要度(MeanDecreaseAccuracy)
#ジニ係数による重要度(MeanDecreaseGini)
varImpPlot(ionome.rf.class)

3-4 分類の結果を partialPlot で視る

partialPlot では、縦軸の値が高いほど、そのサンプルが第1群（Control）である可能性が高く
逆に縦軸の値が低いほど、第2群（Dorought）である可能性が高いということを表しています。
横軸は、各説明変数 (元素) の値です。

n <- 10 # 上位 10 の説明変数を選んで視る
Independent #説明変数を確認

##  [1] "As" "Br" "Ca" "Cd" "Cl" "Co" "Cs" "Cu" "Fe" "K"  "Mn" "Mo" "Ni" "P" 
## [15] "Rb" "S"  "Si" "Sr" "Zn"

x <- as.matrix(X2[rownames(X2),Independent])
head(x) # 間違いないか確認

##       As        Br        Ca        Cd       Cl        Co          Cs
## 1 224.49 0.6832959 1001242.0 0.2321977 3.020088 0.2357462 0.006141176
## 2 282.00 0.2932211  773615.0 0.0000000 2.408051 0.3007813 0.006456006
## 3 172.58 0.4511103  833767.4 0.2792589 3.780524 0.4457043 0.015805737
## 4 264.24 0.5942183 1172170.8 0.2934162 1.937571 0.2068109 0.013059519
## 5 371.57 0.1154127  437902.6 0.0000000 1.744277 0.4555889 0.024038441
## 6 214.99 0.1077018  359460.6 0.1700242 2.096367 0.3158041 0.018823619
##          Cu       Fe        K        Mn         Mo        Ni       P
## 1 0.1854331 1.793914 15.51089 2.1882190 0.04080477 0.3719638 8156.19
## 2 0.1623305 2.395069 13.17087 2.3727699 0.03105810 0.3407753 5241.50
## 3 0.2892326 2.462523 21.12123 1.7866403 0.02061806 0.5631262 2269.07
## 4 0.2097509 1.611651 13.71322 2.5631044 0.03931035 0.3620368 2690.66
## 5 0.3380146 2.114425 41.40882 0.7797614 0.02599244 0.6624682 7269.51
## 6 0.2121268 2.372219 39.30810 1.1544799 0.04730723 0.3925130 2622.33
##           Rb        S       Si        Sr        Zn
## 1 0.06334190 25185.47  2179.44 1.0696411 0.3093640
## 2 0.00000000 11879.18  7687.41 0.6761877 0.1620859
## 3 0.01471262  8545.92  6359.48 1.1392773 0.2041940
## 4 0.00000000 10290.78 11247.67 0.8385686 0.1746725
## 5 0.04671583  6760.89   867.46 0.4219804 0.2506585
## 6 0.06387170 12238.02  2768.05 0.3629122 0.4155615

# 寄与度の高い順に説明変数の番号 (1からnまで) を取得
rk <- order(importance$MeanDecreaseGini, decreasing = TRUE)[1 : n] 
rk

##  [1] 18 13  3  9  5 10  2  8 14 19

#
par(family = "HiraKakuProN-W3")
par(mfrow = c(2, 5)) 
for(i in rk){
  partialPlot(ionome.rf.class, 
              x, 
              Independent[i],
              main = paste0("葉中の",Independent[i]), 
              xlab = Independent[i], 
              ylab = "Partial Dependency", 
              col = "red")
}

箱ひげ図で傾向を確認

rk

##  [1] 18 13  3  9  5 10  2  8 14 19

colnames(x)

##  [1] "As" "Br" "Ca" "Cd" "Cl" "Co" "Cs" "Cu" "Fe" "K"  "Mn" "Mo" "Ni" "P" 
## [15] "Rb" "S"  "Si" "Sr" "Zn"

par(mfrow = c(2, 5)) 
par(family = "HiraKakuProN-W3")
for(i in rk){
  boxplot(x[, i] ~ X2$BlockID, 
          main = paste0("葉中の",colnames(x)[i]), 
          col = 2:3)
}

3-5 学習の収束状況を視る

plot(ionome.rf.class)

3-6 分類の交差検定 (クロスバリデーション)

下準備１(データの再確認をしよう)

library(randomForest)
setwd("/Users/ohmoriohmori/Desktop/アグリバイオ/講義関連/2019年度/ionome/")
Con <- read.csv(("C1C2.csv"), row.names=1)
Dro <- read.csv(("D1D2.csv"), row.names=1)
Y <- rbind(Con, Dro)
head(Y)

dim(Y)

## [1] 384  24

#
Dependent <-"BlockID"
Independent <- c("As", "Br" , "Ca", "Cd", "Cl",  "Co" , "Cs" , "Cu", "Fe", "K", "Mn" ,"Mo", "Ni" ,"P", "Rb", "S", "Si" ,"Sr", "Zn")
X2 <- na.omit(Y[,c(Dependent,Independent)])
head(X2) # BlockID とionome (元素) データだけになっているか確認

dim(X2) # NA削除後のデータ数を確認

## [1] 371  20

下準備２モデルに使用する x (Independent、説明変数、独立変数) と y (Dependent、目的変数、従属変数) を作る

x <- as.matrix(X2[rownames(X2),Independent])
head(x)# ionome データだけになっているか確認

##       As        Br        Ca        Cd       Cl        Co          Cs
## 1 224.49 0.6832959 1001242.0 0.2321977 3.020088 0.2357462 0.006141176
## 2 282.00 0.2932211  773615.0 0.0000000 2.408051 0.3007813 0.006456006
## 3 172.58 0.4511103  833767.4 0.2792589 3.780524 0.4457043 0.015805737
## 4 264.24 0.5942183 1172170.8 0.2934162 1.937571 0.2068109 0.013059519
## 5 371.57 0.1154127  437902.6 0.0000000 1.744277 0.4555889 0.024038441
## 6 214.99 0.1077018  359460.6 0.1700242 2.096367 0.3158041 0.018823619
##          Cu       Fe        K        Mn         Mo        Ni       P
## 1 0.1854331 1.793914 15.51089 2.1882190 0.04080477 0.3719638 8156.19
## 2 0.1623305 2.395069 13.17087 2.3727699 0.03105810 0.3407753 5241.50
## 3 0.2892326 2.462523 21.12123 1.7866403 0.02061806 0.5631262 2269.07
## 4 0.2097509 1.611651 13.71322 2.5631044 0.03931035 0.3620368 2690.66
## 5 0.3380146 2.114425 41.40882 0.7797614 0.02599244 0.6624682 7269.51
## 6 0.2121268 2.372219 39.30810 1.1544799 0.04730723 0.3925130 2622.33
##           Rb        S       Si        Sr        Zn
## 1 0.06334190 25185.47  2179.44 1.0696411 0.3093640
## 2 0.00000000 11879.18  7687.41 0.6761877 0.1620859
## 3 0.01471262  8545.92  6359.48 1.1392773 0.2041940
## 4 0.00000000 10290.78 11247.67 0.8385686 0.1746725
## 5 0.04671583  6760.89   867.46 0.4219804 0.2506585
## 6 0.06387170 12238.02  2768.05 0.3629122 0.4155615

dim(x)

## [1] 371  19

y <- X2[,Dependent]
length(y) # x の行数とあっているか確認

## [1] 371

sample関数を用いて無作為にデータ全体をトレーニングセットとテストデータに分割する

sample関数の利用例
a <- 1:10
sample(a)  # sample関数で無作為に並べ替えをする
sample(a)  # sample関数は実行する毎に異なる並びになる

n.fold <- 10   # ここを10にすると10分割交差検証になる
id <- 1:length(y) %% n.fold + 1 # %%は剰余の計算(割った余り)、+1 して 0 表示をなくしている
id

##   [1]  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4
##  [24]  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7
##  [47]  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10
##  [70]  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3
##  [93]  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6
## [116]  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9
## [139] 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2
## [162]  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5
## [185]  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8
## [208]  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1
## [231]  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4
## [254]  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7
## [277]  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10
## [300]  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3
## [323]  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6
## [346]  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9
## [369] 10  1  2

cv.id <- sample(id) # 振られた番号を並び替える
cv.id

##   [1]  3  1  6  6  9  3  8  6  4  6  8  9  2 10  3  2  2  2  1  8  4  5  9
##  [24]  1  7  7 10  3  3  3  8 10  2 10  4  4  7  7  6  5  6  4 10  2  1  3
##  [47]  6  1  6  7  6  3  1 10  6  8  8  2  9  8  6  2  7  3  5  8 10  6  4
##  [70]  4 10  3  4  1  7  4  7  4  9  3  4  3  5  9  7  3  4  2  8  1  3  1
##  [93]  5  9 10  4  9  9  5 10  7  1  2  2 10  6  5  2  1  5  6  7  9  6  4
## [116]  3  5 10  8  8  9  8  3  9  5  4  2  1  8  3  3  8  8  6  9  6 10  9
## [139]  9  6  7  7  7  2  6  2  4  4  4  7  7  9  5  3 10  1  2 10  1  2  4
## [162]  7  7  6  5  8  9  3  3  5 10  1  8  2  5 10  6 10  4  8  8  9  7  1
## [185]  6  2 10  8 10  7  3  9  5  7  5  9  8  1  4  6  6  3  9  4  6  8  7
## [208]  8  3  9  7  4  4  1  2  9  2  2  1  7  3  3  1  4  8  5  1  9  6  5
## [231] 10  6  8  2  3  1  6  2  8  5  3  5  3  9  1  6  1  4  1  9  8  7  6
## [254]  4 10  5  2  9  4  9  2  8  9  1  1  8 10  9  7  8  2  7  7 10  4  5
## [277]  6  2 10  6  4  2  1 10  1  1 10  3  5  7  6  3 10  3  3  5  9 10  2
## [300]  2  5  3  2 10  1 10  6  7  9  9  7  9  4  5  9  2  5  7  7  3  5  4
## [323] 10  5  8  4  1  7  4  5  7  2  5  5  5 10  5  7  2  4  7 10  4  2  1
## [346]  8  1  1  3  5  1  6  2  8  4 10  1 10  3  9  8 10  6  6  8  8  5  9
## [369]  2  5  8

y.train <- y[cv.id != 1] #1分割目を除いたものをトレーニングセット、除かれた1分割をテストデータとする
X.train <- x[cv.id != 1, ] 
length(y.train)

## [1] 334

dim(X.train)

## [1] 334  19

分類のモデルを作る

set.seed(1000)
Result <- randomForest(y = y.train, x = X.train, importance=TRUE)
y.pred  <- predict(Result, x[cv.id == 1, ]) #テストデータの予測
(result <- table(y.pred, X2[cv.id == 1, ]$BlockID)) #()で括って結果のテーブルを表示

##       
## y.pred  C  D
##      C 16  1
##      D  0 20

(accuracy_prediction = sum(diag(result)) / sum(result)) #精度の計算（行列の対角合計 / 行列の合計）

## [1] 0.972973

varImpPlot(Result) #重要度グラフ表示

plot(Result) #学習の収束状況をプロットする

以上を、1分割から10分割まで繰り返し行う

y.pred <- rep(NA, length(y))  # y (予測データ) の入れ物を予め準備
all_imp <- NULL  # 10分割クロスバリデーション中 (10回分) の importance を保存するための入れ物も用意
#
for(i in 1:n.fold) {
  print(i)
  y.train <- y[cv.id != i]
  X.train <- x[cv.id != i, ]
　set.seed(1000) 
  Result <- randomForest(y = y.train, x = X.train, importance = TRUE, proximity = TRUE)
  y.pred[cv.id == i] <- predict(Result, x[cv.id==i, ])
  #
  imp <- importance(Result) # 各モデルの importance をimp に代入
  all_imp <- rbind(all_imp, imp) # 10回分を列で結合
}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10

(result <- table(y.pred, X2$BlockID))

##       
## y.pred   C   D
##      1 184   4
##      2   5 178

(accuracy_prediction <- sum(diag(result)) / sum(result))

## [1] 0.9757412

write.csv(all_imp, "class_vari_importance.csv") # 10 回分の importance を csvファイルに出力

イオノーム情報を用いた機械学習 (分類)

Yoshihro Ohmori

2019 年 5 月 21 日