myscratchpad

この記事の続き
以下参考文献。
# 帰無仮説有意性検定 (Null Hypothesis Significance Testing, NHST) なんかやめてグラフで信頼区間を示せばいいじゃない、という内容。この信頼区間の計算にプールされた誤差項などを使う。
# 対比についてはそのうち勉強しよう

Masson, M. E. J., & Loftus, G. R. (2003). Using confidence intervals for graphically based data interpretation. Canadian Journal of Experimental Psychology, 57, 203-220.
http://web.uvic.ca/psyc/masson/ML.pdf
Loftus, G. R., & Masson, M. E. J. (1994). Using confidence intervals in within-subject designs. Psychonomic Bulletin & Review, 1, 476-490.
http://web.uvic.ca/psyc/masson/LM.pdf
## 別に関係ないのだが、この二つはほとんど内容が同じだ。こういうのは二重投稿にならないのだろうか

## Masson & Loftus (2003) の例をやってみる

library(car)
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1273264153")

dat <- data.frame(
incong=c(784,853,622,954,634,751,918,894),
cong=c(632,702,598,873,600,729,877,801),
neutral=c(651,689,606,855,595,740,893,822))


# 参加者間1要因3水準, Figure 1
sdat <- stack(dat)
names(sdat) <- c("rts", "cnd")
lmres.b <- lm(rts~cnd, sdat)
Anovares.b <- Anova(lmres.b)
Ares.b <- Astat(Anovares.b)
Ares.b$unv
## 信頼区間の算出
(msew <- Ares.b$unv[2,1]/Ares.b$unv[2,2])
(nj <- 8) # 各群の人数。違う場合はそれぞれ変える
(tcr <- qt(0.975, 21)) # mseの自由度をt値の自由度にする
(ci <- sqrt(msew/nj)*tcr) # 信頼区間。p205
# 作図
mv <- tapply(sdat[,1], sdat[,2], mean)
mv <- mv[c(2,1,3)]; names(mv) <- c("Incong.", "Cong.", "Neutral")
(ul <- mv+ci) # 上側
(ll <- mv-ci) # 下側
bp <- barplot(mv, ylim=c(0,1100), xlab="Group", ylab="Response Latency")
arrows(x0=bp, y0=ul, x1=bp, y1=ll, angle=90, code=3)


# 参加者内1要因3水準 Figure2
tfact <- factor(c("incong", "cong", "neutral"))
idat <- data.frame(tfact)
lmres.w <- lm(cbind(incong, cong, neutral)~1, dat)
Anovares.w <- Anova(lmres.w, idata=idat, idesign=~tfact)
(Ares.w <- Astat(Anovares.w))
## 信頼区間の算出
(msesc <- Ares.w$unv[2,3]/Ares.w$unv[2,4])
(nj <- 8) # 各群の人数。違う場合はそれぞれ変える
(tcr <- qt(0.975, 14)) # mseの自由度をt値の自由度にする
(ci <- sqrt(msesc/nj)*tcr) # 信頼区間。p207。mseを参加者数で割る。mse自体が平方和を自由度で割ったものだが、さらに参加者数で割るのは、サンプリング分布の区間推定ということ (?)
mean(dat) + ci
mean(dat) - ci

# 分散分析の仮定
## 要約
#参加者間デザインで分散の等質性が崩れていた場合、方程式 (1) に基づき普通の信頼区間を出す
#参加者内デザインの場合、球面性仮定を見て、崩れていた場合はプールされた誤差項を使った信頼区間は適切ではない (分散分析自体も適切ではない) 。実際的にはイプシロンを見る。イプシロンが0.75未満はF値の修正を行う。信頼区間についてはプールされた誤差項ではなく調べたい対をもってきてその誤差項で個々に分析し、信頼区間を出す。対なので、1条件につき2つの信頼区間が算出されることもある (Figure3) 。

## 参加者内1要因3水準データで
# 球面性とイプシロンを確認する p207-
Ares.w
# incongruentとcongruent
tfactic <- factor(c("incong", "cong"))
idatic <- data.frame(tfactic)
lmres.ic <- lm(cbind(incong, cong)~1, dat)
Anovares.ic <- Anova(lmres.ic, idata=idatic, idesign=~tfactic)
(Ares.ic <- Astat(Anovares.ic))
## 信頼区間の算出
(msescic <- Ares.ic$unv[2,3]/Ares.ic$unv[2,4])
(nj <- 8) # 各群の人数。違う場合はそれぞれ変える
(tcr <- qt(0.975, 7)) # mseの自由度をt値の自由度にする
sqrt(msescic/nj)*tcr # 信頼区間。p208

### 多要因デザイン p208- ##
## 多要因デザインの各条件の信頼区間は基本的に交互作用項のmseを用いる。交互作用について正確に図示するには差得点の信頼区間を対比を用いて別途算出する必要がある…がそんなグラフは一般的ではないので今回は勘弁してやろう
## 参加者間多要因デザイン (各要因が2水準のとき) p208-
# 要約 群内平均平方和を使って方程式2とほぼ同じようにやる。分散が2：1みたいに等質じゃないときは個々に信頼区間を出す。Table 4が分散分析の結果
# 仮想データ (ローデータはない) Table3. p209. 完全参加者間2x2デザイン。従属変数は正答率。各セル12、48人の参加者。
(table3 <- data.frame(factorA=c("a1", "a1", "a2", "a2"), factorB=c("b1", "b2", "b1", "b2"), means=c(0.51, 0.50, 0.58, 0.61), sds=c(0.08, 0.09, 0.11, 0.08), n=rep(12, 4)))
## 信頼区間の算出 # 分散分析表から
(msew <- 0.009)
n <- 12
(tcr <- qt(0.975, 44)) # t分布の値がちょっと違うのはなんでかなー
sqrt(msew/n)*tcr
# 対比を使おう

## 参加者内多要因デザイン p 211-
## 要約　全ての誤差項を足し、全ての自由度を足した値で割って分子となる平均平方和を出す。
# 仮想データ (ローデータはない) Table4. p211. 完全参加者内2x2デザイン。従属変数は反応時間。12人の参加者。
# Table4のデータを参加者
(table4 <- data.frame(factorA=c("a1", "a1", "a2", "a2"), factorB=c("b1", "b2", "b1", "b2"), means=c(588, 504,525,478), sds=c(69,78,90,67), nob=rep(12, 4)))
## 信頼区間の算出 (psd文末誤植修正後)
msesab <- (4665+3672+7045)/(11+11+11)
n <- 12 # 観測値の数。参加者12人、完全被験者内なので12個
tcr <- qt(0.975, 33) # msesabの分母にした自由度を全て足す。交互作用項の自由度ではない
sqrt(msesab/n)*tcr

# 混合要因デザイン　p214
table6 <- data.frame(
df = c(1,38,1,1,38,79),
ss = c(23018,318868,43665,3605,9974,399130),
ms = c(23018,8391,43665,3605,262,NA),
f = c(2.74, NA, 166.35, 13.73, NA,NA),
row.names=c("rp", "within", "prime", "rpxp", "sxrpxp", "total")
)
table6
## 参加者間要因の信頼区間の算出は方程式2に基づく
(msew <- 8391)
(nj <- 40) # 原文では(2)20 なんで2をかける。参加者内水準が2水準なので1人の参加者につき観測値が2つあるから
(tcr <- qt(0.975, 38))
sqrt(msew/nj)*tcr
## 参加者内要因の信頼区間
(msesc <- 262)
(n <- 20)
(tcr <- qt(0.975, 38))
sqrt(msesc/n)*tcr

# 3水準以上の要因がある場合　p216
## 要約: 交互作用が複雑になるが、全ての対比の信頼区間を報告する必要はない。理論的に重要なところだけでOK。被験者内要因がある場合、差得点を算出できるので上記混合デザインのように差得点の信頼区間を出せばよい。
## 信頼区間の算出。Figure8の左に使う
msesc <- (0.615+1.084+0.675)/(23+46+46)
msesc
n <- 24
n
tcr <- qt(0.975, (23+46+46))
sqrt(msesc/n)*tcr

# 仮想データを生成して同じ事をやってみる。乱数なので結果は違う
# Table7 仮想データ生成
dat <- data.frame(
i.sem = rnorm(24, mean=0.66, sd=0.11),
i.nsm = rnorm(24, mean=0.51, sd=0.16),
i.new = rnorm(24, mean=0.30, sd=0.11),
e.sem = rnorm(24, mean=0.32, sd=0.16),
e.nsm = rnorm(24, mean=0.47, sd=0.20),
e.new = rnorm(24, mean=0.28, sd=0.14)
)

library(car)
iefact <- factor(c("inc", "inc", "inc", "exc", "exc", "exc"))
smfact <- factor(c("sem", "nsm", "new", "sem", "nsm", "new"))
idat <- data.frame(iefact, smfact)
lmres <- lm(cbind(i.sem, i.nsm, i.new, e.sem, e.nsm, e.new)~1, dat)
Anovares <- Anova(lmres, idata=idat, idesign=~iefact*smfact)
Astat(Anovares)
tb <- Astat(Anovares)$unv
tb

## 信頼区間の算出
msesc <- sum(tb[2:4, 3])/sum(tb[2:4, 4])
msesc
n <- 24
n
tcr <- qt(0.975, sum(tb[2:4, 4]))
sqrt(msesc/n)*tcr