myscratchpad

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

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))

##　心理統計学の基礎p226 より。感謝

x1 <- c(12,12,7,17,14,9,10,13,15,12,12,15,11,14,17,17,16,15,15,10,12,9,12,12,19,11,14,15,15,15,16,15,12,10,11,12,15,13,15,12,12,12,13,17,13,11,14,16,12,12)
x2 <- c(2,2,2,3,2,2,3,3,3,1,3,3,2,2,4,2,4,3,4,2,2,1,2,2,4,2,3,2,3,3,2,3,2,2,3,1,2,3,2,2,2,3,3,3,2,3,2,4,2,2)
y <- c(6,11,11,13,13,10,10,15,11,11,16,14,10,13,12,15,16,14,14,8,13,12,12,11,16,9,12,13,13,14,12,15,8,12,11,6,12,15,9,13,9,11,14,12,13,9,11,14,16,8)


dat <- data.frame(x1, x2, y)
library(psych)
describe(dat)[2:4]
cor(dat)

r12 <- cor(x1,x2)
ry1 <- cor(x1,y)
ry2 <- cor(x2,y)

# 部分相関 (semipartial correlation, part correlation)
## x2からx1の影響を除き、そのx2とyとの相関を求める
ry21 <- (ry2 -ry1*r12)/sqrt(1-r12^2)
ry21
## x1からx2の影響を除き、そのx1とyとの相関を求める
ry12 <- (ry1 -ry2*r12)/sqrt(1-r12^2)
ry12
## 関数にしておこう
spcor <- function(x,y,el) {
  rxe <- cor(x,el)
  ryx <- cor(x,y)
  rye <- cor(el,y)
  cat("xからelの影響を除き、そのxとyの相関\n")
  return((ryx -rye*rxe)/sqrt(1-rxe^2))
}
spcor(x2, y, x1)
ry21

# 偏相関 (partial correlation)
## x2からx1の影響を除き、yからもx1の影響を除き、そのうえでx2とyとの相関をとる
pry21 <- (ry2-ry1*r12)/((sqrt(1-ry1^2))*(sqrt(1-r12^2)))
pry21
## 部分相関から偏相関を求める
ry21
pry21*(sqrt(1-ry1^2))
## 偏相関から部分相関を求める
pry21
ry21/sqrt(1-ry1^2)
## 関数にしておこう
pcor <- function(x,y,el) {
  rxe <- cor(x,el)
  ryx <- cor(x,y)
  rye <- cor(el,y)
  cat("xからelの影響を除き、yからelの影響を除き、そのxとyの相関\n")
  return((ryx-rye*rxe)/((sqrt(1-rye^2))*(sqrt(1-rxe^2))))
}
pcor(x2, y, x1)
pry21

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)

multilevelパッケージのsobel関数
http://davidakenny.net/dtt/mediate.htmのデータ。感謝

library(foreign)
dat <- data.frame(read.spss("http://davidakenny.net/dtt/morse_et_al.sav"))
names(dat) <- abbreviate(names(dat))
dmtr <- as.integer(dat$trtm=="Experimentals")

Treatment <- dmtr
DayHoused <- dat$hsng
HousingContacts <- dat$sr_h

library(psych)
print(corr.test(data.frame(Treatment, DayHoused, HousingContacts)), digits=3)

dat2 <- data.frame(Treatment, DayHoused, HousingContacts)
sapply(dat2, mean)
sd(dat2)

library(multilevel)
(sbl <- sobel(Treatment, HousingContacts, DayHoused))
(res <- rbind(c=sbl$Model.1[2,], a=sbl$Model.3[2,], b=sbl$Model.2[3,], "c'"=sbl$Model.2[2,]))

# 覚書
# 1. 予測変数が目的変数を予測する(6.558)
# 2. 予測変数が媒介変数を予測する(0.289)
# 3. 予測変数と媒介変数を投入して目的変数を予測する。このとき、媒介変数が目的変数を予測し(10.710) 、 予測変数、媒介変数を投入したとき、1. と比べて予測変数の係数が低い(3.468)


# sobel 関数の中身は、重回帰分析を3回やってるだけ
  summary(lm(DayHoused~Treatment)) # 1．　予測変数で目的変数を予測
  summary(lm(HousingContacts~Treatment)) # 2．　予測変数で媒介変数を予測
  summary(lm(DayHoused~Treatment+HousingContacts)) # 3．　予測変数と媒介変数を投入して目的変数を予測。予測変数 (Treatment) が有意じゃなくなって、媒介変数 (HousingContacts) が有意になればOK


#      Housing Contacts
#             /\
#            /  \
#           /    \
#          /      \
# 0.289*  /        \  10.710*
#        /          \
#       /            \
#      /              \
#     /                \
#Treatment-----------> Days Housed
#        3.468 (6.558*)


# Aroian, Goodmanのz値を出してみる。
# 参考: http://people.ku.edu/~preacher/sobel/sobel.htm

a <- sbl$Model.3[2,1] # HousingContacts~Treatment　モデルの、HousingContacts (媒介変数) に対するTreatment (予測変数) の回帰係数
b <- sbl$Model.2[3,1] #  DayHoused～Treatment+HousingContacts モデルの、HousingContacts (媒介変数) の回帰係数
a*b # いわゆる間接効果。これを以下の標準誤差で標準化し、z得点化して検定するのがsobel test
sga <- sbl$Model.3[2,2]
sgb <- sbl$Model.2[3,2]

sobeldn <- sqrt(a^2*sgb^2+b^2*sga^2)
aroiandn <- sqrt(a^2*sgb^2+b^2*sga^2+sga^2*sgb^2)
Goodmandn <- sqrt(a^2*sgb^2+b^2*sga^2-sga^2*sgb^2)

# Sobel
a*b/sobeldn
# Aroian
a*b/aroiandn
# Goodman
a*b/Goodmandn

PSYCH-R Home Page
Rと心理学のメーリングリスト

Edinburgh Psychology R-users
エジンバラ大のRwiki (?)
sem用のOpenMxというのがあるらしい


Learning R for Researchers in Psychology
上記リンクのあった元記事

最近リンク貼るばっかだな…