草薙の研究ログ

英語の先生をやってます。

Rで反応時間データの基礎的処理まとめ

これで外国語教育でやるような処理は大体できると思う。ご自由にどうぞ。

library(retimes)	#ない場合はインストールすること
library(MASS)
library(ks)

x<-rexgauss(1000,300,200,500)	#数値例の生成

#基礎
x<-x[!is.na(x)]	#欠損の除外(ここにはない)
summary(x)	#5数要約+平均
skew(x)	#歪度
quantile(x,seq(0,1,.1))	#10%点ずつ
quantile(x,seq(0,1,.05))	#5%点ずつ
data.frame("Counts"=hist(x)$counts,"Breaks"=hist(x)$breaks[-1])	#度数分布表

#外れ値
y1<-x[x<mean(x)+2*sd(x)]	#2SD以上を除去
length(y1)/length(x)	#除去率
y2<-x[x<mean(x)+3*sd(x)]	#3SD以上を除去
y3<-ifelse(x>mean(x)+2*sd(x),mean(x),x)	#2SD以上を平均値に置き換え
y4<-ifelse(x>mean(x)+3*sd(x),mean(x),x)	#3SD以上を平均値に置き換え
y5<-ifelse(x>mean(x)+2*sd(x),mean(x)*2+sd(x),x)	#2SD以上を2SDに置き換え
y6<-ifelse(x>mean(x)+3*sd(x),mean(x)*3+sd(x),x)	#3SD以上を3SDに置き換え

#変換
y6<-log(x,10)	#対数変換(底を10)
y7<-log(x,2)	#対数変換(底を2)
y8<-1/x	#逆数変換

#最尤法による分布へのフィット
m1<-timefit(x)	#Ex-Gaussian分布へのフィット
m1@par	#パラミター
m1@logLik	#対数尤度
AIC(m1)	#AIC
BIC(m1)	#BIC
m1<-timefit(x,iter=100)	#100個のブートストラップ標本
m1@bootPar	#ブートストラップ標本ごとの推定値
apply(m1@bootPar,2,quantile,c(.025,.975))	#パーセンタイル法による信頼区間
apply(m1@bootPar,2,sd)	#ブートストラップ標準誤差

m2<-fitdistr(x,densfun="normal")	#正規分布へのフィット
m3<-fitdistr(x,densfun="lognormal")	#対数正規分布へのフィット
m4<-fitdistr(x,densfun="gamma")	#ガンマ分布へのフィット
m5<-fitdistr(x,densfun="weibull")	#ワイブル分布へのフィット

m2$estimate	#推定値
m2$sd	#標準誤差
m2$loglik	#対数尤度
AIC(m2)	#AIC
BIC(m2)	#BIC


#モーメント法によるEx-Gaussian分布へのフィット
mexgauss(x)

#可視化
boxplot(x,range=F,horizontal=F)	#箱ひげ図
k1<-kde(x)	#カーネル
plot(k1)	#図示
hist(x)	#ヒストグラム


#推定した母数をもつEx-Gaussian分布の確率密度曲線
q<-seq(0,max(x)*1.2,1)
plot(q,dexgauss(q,m1@par[[1]],m1@par[[2]],m1@par[[3]]),type="l",ylab="p",xlab="RT")


#2017/4/6追記

ある読者の方から,外れ値の除去を行列でやるにはどうしたらよいか?という質問をいただきました。ありがとうございます。たくさんやり方はありますが,以下のようにすれば最も単純ではないでしょうか。

#このような実験では最も一般的なデータ形式,
#つまり,行に項目,列に被験者の情報が入っている行列
#またはデータフレームdatがあるとして,

hazure<-function(x){
x[x<mean(x)+2*sd(x)]
}

apply(dat,2,hazure)

#これでOK
#これで除去したもののみを返しているが,
#これではデータの整理がつかないので,
#この外れ値除去後に外れ値の平均を取る。
#つまり,外れ値を除去した被験者平均。

hazure.mean<-function(x){
mean(x[x<mean(x)+2*sd(x)])
}

apply(dat,2,hazure.mean)

これでどうでしょうか?