本講義では、Rを用いて GSのための予測モデル構築を行う方法について説明する。
必要なパッケージ
require(here)
require(gaston)
require(glmnet)
require(BGLR)
require(ranger)
require(elmNNRcpp)例として利用するデータ
ここでは、Rを用いてGSモデル構築を行う方法を解説する。 データ例として、Zhaoら（2011; Nature Communications 2:467）がイネ遺伝資源のGWASに用いた表現型データと、 McCouchら（2016; Nature Communication 7:10532）が用いたHigh Density Rice Array (HDRA)を用いてジェノタイピングされたSNP遺伝子型データ（いずれも、Rice Diversityからダウンロードできる）を用いる。 後者のデータは、70万SNPsのデータのうち、他のマーカーとの冗長性が無く、遺伝子型データの欠測率が低く（< 1%）、かつ、マイナーアリル頻度（minor allele frequency: MAF）が高い（> 5%）、38,769 SNPsのデータを用いる。なお、欠測値については、Beagle 5.1（Blowning et al. 2018）を用いて補完をしてある。
データの読み込み
ゲノムワイドな一塩基多型（SNP）データの読み込み
まず、vcf 形式のマーカー遺伝子型データ(McCouch ら 2016) を読み込む。vcf 形式のデータの読み込みには、gastonパッケージに含まれる関数read.vcfを用いる。
bm <- read.vcf("data3/HDRA-G6-4-RDP1-RDP2-NIAS.AGCT_MAF005MIS001NR_imputed.vcf.gz") ## ped stats and snps stats have been set.
## 'p' has been set.
## 'mu' and 'sigma' have been set.bm # show the summary of the file## A bed.matrix with 1568 individuals and 38769 markers.
## snps stats are set
## ped stats are setイネ遺伝資源アクセッション1,568系統についての、38,769 SNPsの遺伝子型データが含まれている。
系統/品種の属性情報と表現型データの読み込み
次に、おなじ遺伝資源について得られた表現型データを読み込む。
# read phenotypic data
pheno <- read.csv("data3/RiceDiversityPheno4GWASGS.csv", row.names = 1)
head(pheno) # show the data summary## Flowering.time.at.Arkansas Flowering.time.at.Faridpur
## NSFTV1@86f75d2b.0 75.08333 64
## NSFTV3@5ef1be74.0 89.50000 66
## NSFTV4@81d03b86.0 94.50000 67
## NSFTV5@5533f406.0 87.50000 70
## NSFTV6@0d125c0e.0 89.08333 73
## NSFTV7@e37be9e5.0 105.00000 NA
## Flowering.time.at.Aberdeen FT.ratio.of.Arkansas.Aberdeen
## NSFTV1@86f75d2b.0 81 0.9269547
## NSFTV3@5ef1be74.0 83 1.0783133
## NSFTV4@81d03b86.0 93 1.0161290
## NSFTV5@5533f406.0 108 0.8101852
## NSFTV6@0d125c0e.0 101 0.8820132
## NSFTV7@e37be9e5.0 158 0.6645570
## FT.ratio.of.Faridpur.Aberdeen Culm.habit Leaf.pubescence
## NSFTV1@86f75d2b.0 0.7901235 4.0 1
## NSFTV3@5ef1be74.0 0.7951807 7.5 0
## NSFTV4@81d03b86.0 0.7204301 6.0 1
## NSFTV5@5533f406.0 0.6481481 3.5 1
## NSFTV6@0d125c0e.0 0.7227723 6.0 1
## NSFTV7@e37be9e5.0 NA 3.0 1
## Flag.leaf.length Flag.leaf.width Awn.presence
## NSFTV1@86f75d2b.0 28.37500 1.2833333 0
## NSFTV3@5ef1be74.0 39.00833 1.0000000 0
## NSFTV4@81d03b86.0 27.68333 1.5166667 0
## NSFTV5@5533f406.0 30.41667 0.8916667 0
## NSFTV6@0d125c0e.0 36.90833 1.7500000 0
## NSFTV7@e37be9e5.0 36.99000 1.5500000 1
## Panicle.number.per.plant Plant.height Panicle.length
## NSFTV1@86f75d2b.0 3.068053 110.9167 20.48182
## NSFTV3@5ef1be74.0 4.051785 143.5000 26.83333
## NSFTV4@81d03b86.0 3.124565 128.0833 23.53333
## NSFTV5@5533f406.0 3.697178 153.7500 28.96364
## NSFTV6@0d125c0e.0 2.857428 148.3333 30.91667
## NSFTV7@e37be9e5.0 2.859021 119.6000 24.62500
## Primary.panicle.branch.number Seed.number.per.panicle
## NSFTV1@86f75d2b.0 9.272727 4.785975
## NSFTV3@5ef1be74.0 9.166667 4.439706
## NSFTV4@81d03b86.0 8.666667 5.079709
## NSFTV5@5533f406.0 6.363636 4.523960
## NSFTV6@0d125c0e.0 11.166667 5.538646
## NSFTV7@e37be9e5.0 13.333333 5.338019
## Florets.per.panicle Panicle.fertility Seed.length Seed.width
## NSFTV1@86f75d2b.0 4.914658 0.879 8.064117 3.685183
## NSFTV3@5ef1be74.0 4.727388 0.750 7.705383 2.951458
## NSFTV4@81d03b86.0 5.146437 0.935 8.237575 2.926367
## NSFTV5@5533f406.0 4.730921 0.813 9.709375 2.382100
## NSFTV6@0d125c0e.0 5.676468 0.871 7.118183 3.281892
## NSFTV7@e37be9e5.0 5.435903 0.907 7.347255 2.528891
## Seed.volume Seed.surface.area Brown.rice.seed.length
## NSFTV1@86f75d2b.0 2.587448 3.914120 5.794542
## NSFTV3@5ef1be74.0 2.111145 3.658709 5.603658
## NSFTV4@81d03b86.0 2.154759 3.713856 6.120342
## NSFTV5@5533f406.0 1.940221 3.700054 7.089883
## NSFTV6@0d125c0e.0 2.202360 3.657370 5.136483
## NSFTV7@e37be9e5.0 1.774465 3.474791 5.485345
## Brown.rice.seed.width Brown.rice.surface.area
## NSFTV1@86f75d2b.0 3.113958 3.511152
## NSFTV3@5ef1be74.0 2.566450 3.287373
## NSFTV4@81d03b86.0 2.466983 3.334126
## NSFTV5@5533f406.0 2.044608 3.297792
## NSFTV6@0d125c0e.0 2.803858 3.284328
## NSFTV7@e37be9e5.0 2.188882 3.101569
## Brown.rice.volume Seed.length.width.ratio
## NSFTV1@86f75d2b.0 7.737358 2.188
## NSFTV3@5ef1be74.0 5.105350 2.611
## NSFTV4@81d03b86.0 5.126958 2.815
## NSFTV5@5533f406.0 4.085133 4.076
## NSFTV6@0d125c0e.0 5.555183 2.169
## NSFTV7@e37be9e5.0 3.591045 2.905
## Brown.rice.length.width.ratio Seed.color Pericarp.color
## NSFTV1@86f75d2b.0 1.861 0 0
## NSFTV3@5ef1be74.0 2.183 0 0
## NSFTV4@81d03b86.0 2.481 0 0
## NSFTV5@5533f406.0 3.468 0 0
## NSFTV6@0d125c0e.0 1.832 0 0
## NSFTV7@e37be9e5.0 2.506 0 0
## Straighthead.suseptability Blast.resistance Amylose.content
## NSFTV1@86f75d2b.0 4.833333 8 15.61333
## NSFTV3@5ef1be74.0 7.831667 4 23.26000
## NSFTV4@81d03b86.0 NA 3 23.12000
## NSFTV5@5533f406.0 8.333333 5 19.32333
## NSFTV6@0d125c0e.0 8.166667 3 23.24000
## NSFTV7@e37be9e5.0 8.585000 2 20.18667
## Alkali.spreading.value Protein.content
## NSFTV1@86f75d2b.0 6.083333 8.45
## NSFTV3@5ef1be74.0 5.638889 9.20
## NSFTV4@81d03b86.0 5.527778 8.00
## NSFTV5@5533f406.0 6.027778 9.60
## NSFTV6@0d125c0e.0 5.444444 8.50
## NSFTV7@e37be9e5.0 6.000000 11.75
## Year07Flowering.time.at.Arkansas
## NSFTV1@86f75d2b.0 73.0
## NSFTV3@5ef1be74.0 88.0
## NSFTV4@81d03b86.0 105.5
## NSFTV5@5533f406.0 86.5
## NSFTV6@0d125c0e.0 85.5
## NSFTV7@e37be9e5.0 105.0
## Year06Flowering.time.at.Arkansas
## NSFTV1@86f75d2b.0 77.16667
## NSFTV3@5ef1be74.0 91.00000
## NSFTV4@81d03b86.0 83.50000
## NSFTV5@5533f406.0 88.50000
## NSFTV6@0d125c0e.0 92.66667
## NSFTV7@e37be9e5.0 NAphenoには、36形質（一部は、形質と試験地あるいは年次との組合せ）のデータが含まれている。
ここでは、形質例として種子の長さと幅の比（Seed.length.width.ratio）を用いる。種子長のデータを抜き出した後、各要素（各系統に対応）に名前を付ける。その後、na.omit関数を用いて欠測値を除く。
y <- pheno$Seed.length.width.ratio # select seed length width ratio as a target trait
names(y) <- rownames(pheno) # naming the elements of y
y <- na.omit(y) # remove missing elements例として解析する表現型
GSは、GWASでは有意なマーカーが検出されないような形質で特に有効である。ここでは、対象形質を、「個体あたりの稔性のある頴花数」として解析を行う。これを被説明変数\(y\)とする。
y <- pheno$Panicle.number.per.plant * pheno$Florets.per.panicle * pheno$Panicle.fertility
names(y) <- rownames(pheno) # naming the elements of y with the rownames
y <- na.omit(y) # remove NA entriesGWAS
準備した\(y\)について、GWASをおこなってみる。
bm.wp <- select.inds(bm, id %in% names(y))
bm.wp@ped$pheno <- y[bm.wp@ped$id]
gwa <- association.test(bm.wp, method = "lmm", test = "wald", eigenK = eigen(GRM(bm.wp)))
fdr <- p.adjust(gwa$p, method = "BH")
col <- rep("black", nrow(gwa))
col[gwa$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa, pch = 20, col = col)
rm(gwa, bm.wp) # save memory space for RStudio Cloud environmentFDR 5%で有意なSNPは無く、マンハッタンプロットでも顕著なピークは見られない。
ゲノミック予測モデルの構築
ここから、ゲノミックセレクション(genomic selection: GS)で用いられるゲノミック予測(genomic prediction: GP)のためのモデル構築を行ってみる。ここでは、以下の4つの方法を用いて予測モデルを構築して、その精度を比較してみる。
- 正則化線形回帰（リッジ回帰、lasso）(RR, LA)
- ゲノミックBLUP（best linear unbiased predictor: 最良線形不偏予測量）（GB）
- random forest (RF)
- extreme learning machine (EL)
ゲノムデータの取り出し（行列への変換）
MAFが10%よりも大きなSNPにしぼり、その後、ゲノムデータを説明変数の行列\(X\)として取り出す。なお、RStudio Cloudを用いた解析の場合、限られた計算能力のために、全SNP（30,219）を入れた解析を行うと時間を要するだけでなく、時折異常終了することがあった。そこで、マーカー数を1/10に減らすこととする。減らし方は、30個に1個のマーカーを規則的にとることとする。
X.all <- as.matrix(select.snps(bm, maf > 0.1))
dim(X.all)## [1] 1568 30219(sel <- c(T, rep(F, 29)))## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSEX <- X.all[, sel]
dim(X)## [1] 1568 1008rm(bm, X.all) # save memory space説明変数\(X\)について、表現型データ\(y\)がある系統と、無い系統に分ける。最終的には、前者のデータから後者を予測する。後述するように、こうした予測を行うことによって、表現型データの無い遺伝資源についても、有望な系統のスクリーニングを行うことが可能となる。
X.wp <- X[names(y), ] # X with phenotypes (wp)
X.wop <- X[!(rownames(X) %in% names(y)), ] # X without phenotypes (wop)
c(nrow(X.wp), nrow(X.wop))## [1] 356 1212rm(X) # save memory space表現型データ\(y\)を持つ系統が356系統、持たない系統が1212系統ある。
次に、表現型データのある356系統のSNP遺伝子型データX.wpから、多型の無いマーカーを除いておく。多型の有無のチェックには、分散の計算を用いている。分散が0でなければ多型があると判断できる。
poly <- apply(X.wp, 2, var) > 0
X.wp <- X.wp[, poly]
dim(X.wp)## [1] 356 1008次に、表現型データの無い1212系統のSNP遺伝子型データ（X.wop）についても、同じSNPのセットに揃えておく。
X.wop <- X.wop[, poly]正則化線形回帰（RR, LA）
まずは、glmnetパッケージを用いて正則化線形回帰を行う。glmnetでは、alphaというパラメータがあり、これを0に設定するとRR、1に設定するとLASSOを用いてモデル構築が行われる。まずは、alphaを0にして、リッジ回帰を行う。
alpha <- 0 # the value 0 corresponds to RRglmnetパッケージのcv.glmnet関数を用いてリッジ回帰のモデルを構築する。cv.glmnet関数は、回帰係数の大きさに対して課すペナルティの大きさを調整するパラメータを、交差検証（cross-validation）によって最適化してくれる。関数名の最初のcvは交差検証を意味している。説明変数\(X\)の標準化は、ここでは行わない（0, 1, 2と基準化されているので）。
model.ridge <- cv.glmnet(X.wp, y, alpha = alpha, standardize = F)次に、得られたモデルをもとに、マーカー遺伝子型データから表現型データを計算する。ただし、ここでは、モデル構築に用いた系統のSNP遺伝子型データを入力にするため、モデルに基づく「予測」（prediction）ではなく、モデルの「あてはめ」（fitting）である。関数名predictに惑わされないよう注意する。
y.fit <- predict(model.ridge, newx = X.wp, s = "lambda.min") あてはめられた\(y\)と観察値の\(y\)について比較する。両者の相関係数を計算する。
plot(y, y.fit)
abline(0, 1)
cor(y, y.fit)## 1
## [1,] 0.7877959なお、入力される説明変数が多い場合には、あてはまりの良いモデルを構築するのは容易である。したがって、あてはまりの良さは、モデルの予測能力を反映しない。モデルの予測能力を評価するためには、モデル構築に用いなかったデータに対する「予測」を行い、その精度を評価する必要がある。そこで、本当はモデル構築に用いることができるデータの一部を、モデル構築から除いておき、除いておいたデータを予測して、その予測値を観察された値と比較する。除いておいたデータは、「本当はモデル構築に用いることができるデータ」であり、表現型データ（\(y\)）があるデータである。したがって、予測値を計算して、「本当は観察されていた\(y\)の値」との比較が可能である。
\(n\)-fold交差検証とよばれる方法では、無作為にデータを\(n\)グループ分割し、そのうち1グループをモデル構築から除いておき、残りの\(n-1\)クラスで予測モデルを構築する。そして、除いておいた1グループについて、予測を行い、実測値との比較を行う。では、無作為に1〜\(n\)のIDを系統にふってみよう。ここでは、\(n=3\)とする。無作為にIDを並べ替えるにはsample関数を用いる。
n.fold <- 3
head(id <- 1:length(y) %% n.fold + 1, 30)## [1] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1head(cv.id <- sample(id), 30) # randomly reorder## [1] 2 3 2 2 2 1 1 1 1 3 3 3 3 1 1 2 2 3 3 3 3 2 2 3 2 1 2 2 1 3次に、無作為に割り振られたIDが1番の系統を除いて、モデル構築用データを作成してみる。
y.train <- y[cv.id != 1]
X.wp.train <- X.wp[cv.id != 1, ]データを分割したことにより、多型がなくなるマーカーが出現するかもしれない。そこでこれらを改めて取り除いておく。
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]訓練用データでモデルを構築し、1分割目の予測を行なう。
model <- cv.glmnet(X.wp.train, y.train, alpha = alpha, standardize = F)
y.pred <- predict(model, newx = X.wp[cv.id == 1, poly], s = "lambda.min")1分割目について、観察値と予測値を比較する。観察値と予測値の相関係数を精度の指標として計算する。
plot(y[cv.id == 1], y.pred)
abline(0, 1)
cor(y[cv.id == 1], y.pred)## 1
## [1,] 0.65274実際には、全分割について、同様の計算を繰り返し、それらをまとめて精度評価する。そこで、以上の操作を、for文を用いて1〜3の分割について繰り返す。
# repeat this analysis for #1 - #5 folds
y.pred.rr <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y[cv.id != i]
X.wp.train <- X.wp[cv.id != i, ]
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]
model <- cv.glmnet(X.wp.train, y.train, alpha = alpha, standardize = F)
y.pred.rr[cv.id == i] <- predict(model, newx = X.wp[cv.id == i, poly], s = "lambda.min")
}## [1] 1
## [1] 2
## [1] 3最後に全分割の予測値と観察値を合わせて比較する。相関係数も精度の指標として計算する。
plot(y, y.pred.rr)
abline(0, 1)
cor(y, y.pred.rr)## [1] 0.6923296リッジ回帰で推定された各SNPマーカーの効果（回帰係数）を図示してみる。全データを用いてモデルを構築し、それをもとにSNPマーカーの効果を図示する。効果（回帰係数）は、coef関数で取り出せる。
beta.ridge <- coef(model.ridge, s = "lambda.min")[-1] # model.ridge has been obtained above
plot(abs(beta.ridge), main = "Ridge", ylab = "Coefficients", xlab = "Position (bp)")
推定効果をプロットしてみると、特に目立って効果の大きなSNPsがあるわけではなく、多くのSNPがある程度の遺伝効果を持っていることが分かる。
最後に、表現型データがある356系統から得られた予測モデルを、表現型データの無い1212系統に適用して、後者の表現型値を予測する。
# prediction of phenotypes of genotypes without phenotypic records
y.pred.wop <- predict(model.ridge, newx = X.wop, s = "lambda.min")
conf <- hist(c(y.pred.rr, y.pred.wop), col = "#ff0000ff", border = "#ffffff")
hist(y.pred.rr, col = "#0000ffff", border = "#ffffff", breaks = conf$breaks, add = T)
上図は、モデル作成に用いられた表現型データをもつ356系統の交差検証による予測値（青）と、表現型データの無い1212系統の予測値（赤）のヒストグラムである。ここで、個体あたりの稔性頴花数が多いほうが良いとすると、このような予測から、そうした系統が1212系統にどのくらい含まれているかが分かる。
例えば、この値が16を超える系統が1212系統中60以上あることが分かる。
sum(y.pred.wop > 16)## [1] 30最近、多数の遺伝資源についてゲノムデータが解析・公開されるようになってきている。例えば、イネでは遺伝資源3,000系統について全ゲノム配列のリシークエンスが行われ、そのデータが既に公開されている（The 3,000 Rice Genome Project 2014, Gigasicence 3:7; Sun et al. 2017, Nucl. Aci. Res. 45:597-605）。また、今回例として用いているデータも1,500系統以上について70万SNPsの遺伝子型データが解析され、公開されている。こうしたデータを活用することにより、表現型データのスクリーニングが難しいために、死蔵されることが多かった遺伝資源を簡易にスクリーニングできるようになる。したがって、GSを用いることで、有用な遺伝資源を迅速に発見し、育種に効率的に活用する道が開ける（Yu et al. 2016, Nature Plants 2: 16150; Tanaka and Iwata 2017, Theor. Appl. Genet. 131:93-105）。
次に、LASSOを用いた正則化線形回帰を行う。LASSOを用いる場合には、alphaを1にする。なお、alphaを1にする以外は、上で行った計算と全く同じである。
alpha = 1 # do the same thing with alpha = 1
y.pred.la <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y[cv.id != i]
X.wp.train <- X.wp[cv.id != i, ]
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]
model <- cv.glmnet(X.wp.train, y.train, alpha = alpha, standardize = F)
y.pred.la[cv.id == i] <- predict(model, newx = X.wp[cv.id == i, poly], s = "lambda.min")
}## [1] 1
## [1] 2
## [1] 3plot(y, y.pred.la)
abline(0,1)
cor(y, y.pred.la) ## [1] 0.6501848全データを用いて推定されたSNP効果（回帰係数）の絶対値を図示する。
model.lasso <- cv.glmnet(X.wp, y, alpha = alpha, standardize = F)
beta.lasso <- coef(model.lasso, s = "lambda.min")[-1]
plot(abs(beta.lasso), main = "LASSO", ylab = "Coefficients", xlab = "Position (bp)")
リッジ回帰に比べ、効果をもつSNPsの数が少ないことが分かる。また、効果をもつ場合は、リッジ回帰に比べ、効果の大きさが大きいことも分かる。しかし、LASSOを用いた場合でも、複数のSNPsに重み付けされており、この形質に関わる遺伝子の数が少なくないことが示唆される。
最後に、先と同じように、表現型データが無い1212系統の予測値を計算する。
y.pred.wop <- predict(model.lasso, newx = X.wop, s = "lambda.min")
conf <- hist(c(y.pred.la, y.pred.wop), col = "#ff0000ff", border = "#ffffff")
hist(y.pred.la, col = "#0000ffff", border = "#ffffff", breaks = conf$breaks, add = T)
sum(y.pred.wop > 16)## [1] 51用いたモデル化手法とモデルが異なるため、リッジ回帰の場合と異なる結果が得られる。実際に適用する場合には、複数の手法の中から、それぞれの形質に合った手法を経験的に選んで用いるのが良い。どの形質にも万能な手法は無い（データ科学の世界では、これを「No free lunch定理」とよぶ）。
ゲノミックBLUP（GB）
ゲノミックBLUPとよばれるモデルを用いて予測モデルを構築する。そのためには、ゲノム関係行列とよばれるゲノムワイドマーカーの多型に基づいて計算される系統間の遺伝的関係を予め計算しておく必要がある。
grm <- tcrossprod(scale(X.wp)) / ncol(X.wp)計算しておいたゲノム関係行列を用いて、予測モデルを構築する。まずは、全データに対して当てはめる。
eta <- list(list(K = grm, model = "RKHS"))
model.gb <- BGLR(y, ETA = eta, verbose = F)
y.fit <- model.gb$yHat
plot(y, y.fit)
abline(0, 1)
cor(y, y.fit)## [1] 0.8889414次に、交差検証を行う。
y.pred.gb <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y
y.train[cv.id == i] <- NA
model <- BGLR(y.train, ETA = eta, verbose = F)
y.pred.gb[cv.id == i] <- model$yHat[cv.id == i]
}## [1] 1
## [1] 2
## [1] 3plot(y, y.pred.gb)
abline(0,1)
cor(y, y.pred.gb) ## [1] 0.6985989機械学習 Random Forest（RF）
ここでは、機械学習法のひとつ、ランダムフォレスト（RF）を用いて、予測モデルを作成する。関数rangerを用いる。
まずは、全データに対して当てはめる。
df <- data.frame(y, X.wp)
model.rf <- ranger(y ~ ., data = df)
y.fit <- predict(model.rf, data = df)$predictions
plot(y, y.fit)
abline(0, 1)
cor(y, y.fit)## [1] 0.9588809少し、当てはめ過ぎ（over-fitting）のようにも見える。
交差検証を行う。
y.pred.rf <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y[cv.id != i]
X.wp.train <- X.wp[cv.id != i, ]
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]
df <- data.frame(y = y.train, X.wp.train)
model <- ranger(y ~ ., data = df)
y.pred.rf[cv.id == i] <- predict(model, data = data.frame(X.wp[cv.id == i, poly]))$predictions
}## [1] 1
## [1] 2
## [1] 3plot(y, y.pred.rf)
abline(0,1)
cor(y, y.pred.rf) ## [1] 0.7394781予測精度は良い。
機械学習 Extreme Learning Machine (EL)
最後に、階層的ニューラルネットワークの一つであるExtreme Learning Machine（EL）を用いて予測モデルを作成する。
まずは、全てのデータに当てはめる。
model.elm <- elm_train(X.wp, as.matrix(y), nhid = 5000, actfun = 'sig')
y.fit <- elm_predict(model.elm, newdata = X.wp)[,1]
plot(y, y.fit)
abline(0,1)
cor(y, y.fit)## [1] 0.9995071モデルは、完全に観察データに当てはまってしまっている。
次に、交差検証を行い、未知のデータにおける予測精度を評価する。
y.pred.el <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y[cv.id != i]
X.wp.train <- X.wp[cv.id != i, ]
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]
model <- elm_train(X.wp.train, as.matrix(y.train), nhid = 5000, actfun = 'sig')
y.pred.el[cv.id == i] <- elm_predict(model, newdata = X.wp[cv.id == i, poly])[,1]
}## [1] 1
## [1] 2
## [1] 3plot(y, y.pred.el)
abline(0,1)
cor(y, y.pred.el) ## [1] 0.7098838予測精度の比較
最後に、5つの手法の予測値と実測値の比較を行う。5つの手法間の関係も分かるように対散布図を作成する。
df <- data.frame(Obs = y, RR = y.pred.rr, LA = y.pred.la, GB = y.pred.gb, RF = y.pred.rf, EL = y.pred.el)
pairs(df)
予測精度をまとめておく。
(r <- cor(df)[1,-1])## RR LA GB RF EL
## 0.6923296 0.6501848 0.6985989 0.7394781 0.7098838barplot(r)
今回のデータではRFやGBが比較的良い予測精度を示している。
なお、今回は３分割交差検証を一度しか行っていない。確率的変動により、実行するごとに精度は変化する。実際には、複数回の分割検証を行って精度比較を行う必要がある。