myscratchpad

実験デザイン用のパッケージ agricolae

# ラテン方格デザインでのランダマイズ
library(agricolae)
varieties<-c("perricholi","yungay","maria bonita","tomasa")
lsd <-design.lsd(varieties,number=1001,seed=23)
lsd # print field book. plotsは生物統計学の区画。心理学では参加者にあたる（と思う）
plots <-as.numeric(lsd[,1])
trt <-as.character(lsd[,4])
dim(plots)<-c(4,4)
dim(trt) <-c(4,4)
print(t(plots))
print(t(trt))

# 関数一覧
 http://rss.acs.unt.edu/Rdoc/library/agricolae/html/00Index.html
# 開発元
 http://tarwi.lamolina.edu.pe/~fmendiburu/
# CRAN Task View: Design of Experiments (DoE) & Analysis of Experimental Data
 # これは日本語化されると便利かもしれない
 http://cran.at.r-project.org/web/views/ExperimentalDesign.html

st <- "a,b,c,d,e"
st2 <- "f,g,h"
st3 <- "i,j,k,l"
lst <- list(st, st2, st3)
# これを3行6列の表にしたい。空白はNAでうめる
(sts <- lapply(lst, function(x) unlist(strsplit(x, ","))))
lens <- lapply(sts, length)

reslist <- list()
for (i in 1:length(sts)) {
lens <- length(sts[[i]])
ko <- 6
sta <- c(sts[[i]], rep(NA, 6-lens))
reslist[[i]] <- sta
}

res2 <- t(data.frame(reslist))
dimnames(res2) <- NULL
data.frame(res2)

# 行ごとにx内の要素に一致するものの数を数える
apply(res2, 1, function(a) sum(complete.cases(match(a,x))))
# 列ごとにx内の要素に一致するものの数を数える
apply(res2, 2, function(a) sum(complete.cases(match(a,x))))

x <- factor(c(c(10, 20, 30), "AA"))
x

x[x=="AA"] <- NA
mean(x) # 因子なのでエラー
x2 <- as.numeric(as.character(x))
x2
mean(x2, na.rm=T)

前のをちょっと変えた。
正答・誤答のベクトルと比較して、正答の反応時間だけ使うようにした。

nn <- 50
rts<- round(c((rnorm(nn, mean=300, sd=20)+(300*rexp(nn)))))
rsp <- sample(c("correct", "incorrect"), 50, replace=T, prob=c(0.80, 0.20))

elrt <- function(rts, rsp="allcorrect", cres="correct", cop=2.5, ifm=FALSE) {
rts0 <- rts
ifelse(rsp=="allcorrect", rsp<-rep(cres, length(rts)), rsp<- rsp)
cres <- cres
cop <- cop
## 正答以外の反応の添字
irind <- which(rsp!=cres)
rts1 <- rts0
rts1[irind] <- NA
## 正答以外の反応の反応時間を除く
rts2 <- na.omit(rts1)
rtl <- log10(rts2) # 常用対数に変換
meanv <- mean(rtl)
sdv <- sd(rtl)
cv.l <- meanv-(cop*sdv)
cv.u <- meanv+(cop*sdv)
rtl2 <- log10(rts1) # 常用対数に変換
lind <- which(rtl2 < cv.l)
uind <- which(rtl2 > cv.u)
    nart <- rts1
    nart[lind] <- 10^(cv.l) # 対数値を元に戻す
    nart[uind] <- 10^(cv.u)
    rsp2 <- rsp
    rsp2[lind] <- "lower"
    rsp2[uind] <- "upper"
    reslist <- list()
    reslist$nart <- nart
    reslist$logrt <- rtl2
    reslist$responses <- rsp2
    reslist$outliers.l <- rts[lind]
    reslist$outliers.u <- rts[uind]
    reslist$cop <- paste(cop, "SD", sep="")
    reslist$cv <- c(cv.l, length(rts[lind]), cv.u, length(rts[uind]))
    names(reslist$cv) <-  c("lower", "nofl", "upper",  "nofu")
ifelse(ifm==FALSE, return(rtl2), return(reslist))
}
elrt(rts) # 反応ベクトルを指定しないときは全て正答とする
elrt(rts, rsp=rsp, ifm=T)
 # 反応ベクトルがあるときはrspを指定する。指定がないときは全て正反応として扱われる。
 # nartは不正解をNA、外れ値をカットオフ値で置換した生RTを返す。長さは元のrtsベクトルと同じ
 # logrtは上の対数。ベクトルの長さは同じ。メインの分析対象はこれ
## 除去・置換後の反応をとりだす
x <- elrt(rts, rsp, ifm=T)
x[[1]] # rt
x[[2]] # log
x[[3]] # 反応


# 以下のデータフレームのように、各参加者IDと反応時間のベクトルがある場合

# 参加者個々にelrt関数を適用してみる
 # データ生成
id <- gl(3,10)
rts<- round(c((rnorm(30, mean=300, sd=20)+(300*rexp(30)))))
rts[30] <- 5000 # 意図的に外れ値を入れておく
rsp <- sample(c("correct", "incorrect"), 30, replace=T, prob=c(0.80, 0.20))
dat <- data.frame(id, rts, rsp)
dat # 仮データ

res <- list()
ids <- dat[,1]
lvs <- levels(dat[,1]) # id
for (i in 1:nlevels(dat[,1])) {
  rtx <- dat[,2][ids==lvs[i]]
  rsx <- as.character(dat[,3][ids==lvs[i]])
  res1 <- elrt(rtx, rsp=rsx, ifm=T)
  res2 <- res1$logrt # 対数値のみ返す
  res[[i]] <- res2
}
 # tapplyとかmapplyで上手くできるやり方あるのかな…

unlist(res) # ベクトルで返す
 # 参加者個々じゃないと
 elrt(rts, rsp=rsp)

# 全部をつなげてみよう
data.frame(dat, unlist(res), elrt(rts, rsp=rsp))

total sum of squares の求め方
dat <- data.frame(
s=factor(paste("s", 1:40, sep="")),
a=factor(c(rep("a1", 20), rep("a2", 20))),
b=factor(rep(c(rep("b1",5), rep("b2",5), rep("b3",5), rep("b4",5)),2)),
result = c(3,3,1,3,5, 4,3,4,5,7, 6,6,6,4,8, 5,7,8,7,9, 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,5, 2,3,3,4,6))
lmsst <- lm(result~1, dat)
sum(resid(lmsst)^2)
sum((dat$result-mean(dat$result))^2)