#行の中で、「#」以降は、コメントアウトで処理しません。

#====================================================

#ブートストラップを行うRのパッケージは、boot, simpleboot, bootstrapなどがあるが

#本章の分析はブートストラップ標本として平均波形を作り、そこからP300を抽出するという作業が

#あるため単純にこれらのパッケージを利用できない。したがって、ブートストラップの過程を書き下して

#いる。仕組みを学習するとよいでしょう。

#Cドライブにv09dataというフォルダを作り、v09c4ERP.csvを置いているという設定）

#データの読み込み

data=read.table("c:/v09data/v09c4ERP.csv",header=T,sep=",",row.names=1)

#標本サイズの指定。Rは関連データ、Uは非関連データを示す。

nR=20

nU=80

#オリジナルデータから関連データと非関連データを抽出

dataR=head(data,nR)

dataU=tail(data,nU)

#関連と非関連のオリジナルデータから平均波形を求める

avewaveOR=colMeans(dataR)

avewaveOU=colMeans(dataU)

#関連と非関連のオリジナルデータから2つのP300とP300振幅差を求める

P300OR=max(avewaveOR)

P300OU=max(avewaveOU)

P300OD=P300OR-P300OU

#乱数の種を指定、ブートストラップ標本サイズを1000、ブートストラップ標本を格納するベクトルを作成

set.seed(314)

B=1000

P300Dboot=numeric(B)

#ブートストラップ

for (b in 1:B){ #以下の操作を1000回(B回)くりかえす

iR=sample(1:nR,replace=TRUE) #1から20(nR)まで重複を許して乱数を生成

avewaveR=colMeans(dataR[iR,]) #抽出された関連データの平均波形を求める

P300R=max(avewaveR) #作成された平均波形からP300を抽出（注）

iU=sample(1:nU,replace=TRUE) #1から80(nU)まで重複を許して乱数を生成

avewaveU=colMeans(dataU[iU,]) #抽出された非関連データの平均波形を求める

P300U=max(avewaveU) #作成された平均波形からP300を抽出

P300D=P300R-P300U #P300振幅差を計算。Dはdifference

P300Dboot[b]=P300D #ブートストラップ統計量をベクトルに格納

}

#（注）CSV用のデータはあらかじめ300ms〜800msまでのデータを選んでいるので、得られた平均波形から

#最大値を選べば、その値がP300になる。

#ブートストラップ分布（図4-12に相当）の作成

hist(P300Dboot,freq=T,main="P300振幅差",breaks=40)

#ブートストラップ分布の偏り（ブートストラップ統計量の平均値）

mean(P300Dboot)

#パーセンタイル信頼区間(95%)

quantile(P300Dboot,probs=c(0.025,0.975))

#バイアス修正パーセンタイル信頼区間(BCa)法 (汪・桜井,2011,P126参照)

#偏り修正定数の計算

w=qnorm(sum(P300Dboot<=P300OD)/B)

#加速定数acの計算

thetaR=numeric(nR)

thetaU=numeric(nU)

for (i in 1:nR){thetaR[i] = max(colMeans(dataR[-i,]))-max(colMeans(dataU))}

for (i in 1:nU){thetaU[i] = max(colMeans(dataR))-max(colMeans(dataU[-i,]))}

VR=(nR-1)*(mean(thetaR)-thetaR)

VU=(nU-1)*(mean(thetaU)-thetaU)

nu=sum(VR^3)/nR^3+sum(VU^3)/nU^3

de=(sum(VR^2)/nR^2+sum(VU^2)/nU^2)^(3/2)

ac=nu/(6*de)

#BCa信頼区間(95%)の出力

d=w+qnorm(c(0.025,0.975))

alphaBCa=pnorm(w+d/(1-ac*d))

sort(P300Dboot)[round(B*alphaBCa)]

#参考文献 汪金芳・桜井裕仁(2011)ブートストラップ入門 共立出版

#====================================================

#2015.04.09 松田いづみ・荘島宏二郎