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)
これでどうでしょうか?