myscratchpad

pairwise.t.test関数はt値が出力されないのでt値と自由度を出力するようにした。
信頼区間を出してもいいが、p値を調整した信頼区間に何か意味があるのかよくわからないので勘弁してやろう

source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275849015")

## サンプル
dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9, 10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))

ptest(dat$result, dat$a)
## 参考に、pairwise.t.testで
pairwise.t.test(dat$result, dat$a)
## paired=Tにしたとき、あるいはプールされた標準偏差を用いない場合は単にt検定を繰り返しているだけなのでこっちでやる。

## 参考に、青木先生の関数。感謝
source("http://aoki2.si.gunma-u.ac.jp/R/src/Bonferroni.R", encoding="euc-jp")
Bonferroni(dat$result, dat$a,  "Holm")
round(p.adjust(Bonferroni(dat$result, dat$a,  "Holm")$result2[,2], "holm"), 5)

source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275792703")
tobj # 中身の表示
## 要素名はabbreviate関数で短くしてある。
## row names were found from a short variable and have been discarded という警告が出るが気にしない。…いや、できれば消したいがやり方がわからない
### 関係ないが、この日本語訳は"列名は短い変数から見つけられ、捨て去られました"となる。たぶん"列名"ではなく"行名"の誤訳だろう。

## サンプル
x1 <- rnorm(10)
x2 <- rnorm(10)
res <- t.test(x1, x2)
res2 <- t.test(x1, x2, paired=T)
tobj(res)
tobj(res2)

## t検定繰り返しサンプル
dat <- data.frame(a1=c(9,7,8,8,12,11,8,13), a2=c(6,5,6,3,6,7,10,9), a3=c(10,13,8,13,12,14,14,16), a4=c(9,11,13,14,16,12,15,14))
dat

cbn <- combn(names(dat), 2) # 組合せ
tresdat <- tobj(t.test(rnorm(10), rnorm(10), paired=T)); tresdat[1,] <- NA # 結果格納用データフレーム

for (i in 1:ncol(cbn)) {
tres <- t.test(dat[,cbn[1,i]], dat[,cbn[2,i]], paired=T)
tob <- tobj(tres)
tresdat[i,] <- tob
}
rownames(tresdat) <- apply(cbn, 2, function(x) paste(x, collapse=","))
tresdat
## Holm法でp値の調整をする
pholms <- p.adjust(tresdat[,"p.vl"], "holm")
tresdat <- data.frame(tresdat, pholms)
colnames(tresdat) <- c(colnames(tresdat[-ncol(tresdat)]), "Pr(>|t|).Holm")
print.anova(tresdat) ## 有意の星をつける。文字型の変数は01にされる


source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275792703")
## 縦型データで繰り返し、対応なしにしてみる
dat <- data.frame(a1=c(9,7,8,8,12,11,8,13), a2=c(6,5,6,3,6,7,10,9), a3=c(10,13,8,13,12,14,14,16), a4=c(9,11,13,14,16,12,15,14))
dat
datl <- stack(dat)
val <- datl$value
idv <- datl$ind
cbn <- combn(levels(idv), 2)
tresdat <- tobj(t.test(rnorm(10), rnorm(10), paired=F)); tresdat[1,] <- NA # 結果格納用データフレーム
for (i in 1:ncol(cbn)) {
tres <- t.test(val[which(idv==cbn[1,i])], val[which(idv==cbn[2,i])], paired=F)
tob <- tobj(tres)
tresdat[i,] <- tob
}
rownames(tresdat) <- apply(cbn, 2, function(x) paste(x, collapse=","))
tresdat
## Holm法でp値の調整をする
pholms <- p.adjust(tresdat[,"p.vl"], "holm")
(tresdat <- data.frame(tresdat, pholms))
colnames(tresdat) <- c(colnames(tresdat[-ncol(tresdat)]), "Pr(>|t|).Holm")
print.anova(tresdat) ## 有意の星をつける。文字型の変数は01にされる

grp <- factor(rep(c(51,55,57,59,61), each=3))
val <- c(12.2,18.8,18.2,22.2,20.5,14.6, 20.8,19.5,26.3, 26.4,32.6,31.3, 24.5,21.2,22.4)

## Leveneの等分散性の検定。こっちよりもBartlett検定のほうがいいらしい。
## car, lawstat, s20xに同じ関数がある
library(car)
leveneTest(val, grp)
leveneTest(lm(val~grp)) # lm関数のオブジェクトも受け取れる。便利
levene.test(val, grp) # もうすぐなくなるらしい

library(lawstat)
levene.test(val, grp) # デフォルト
levene.test(val, grp, "median") # デフォルトと同じ
levene.test(val, grp, "mean") # spssと同じ出力

# s20xパッケージにも同じ関数がある

## Bartlettの等分散性の検定。Rの組み込み関数
bartlett.test(val, grp)

## 2要因参加者間も同じ。以下のデータを参考。感謝
http://academic.udayton.edu/gregelvers/psy216/SPSS/2wayanovabs.htm

cls <- factor(rep(c("dst", "lct"), each=10))
gpa <- rep(factor(rep(c("high", "low"), each=5)), 2)
pnt <- c(332,380,371,366,354,259.5,302.5,296,349,309,354.67,353.5,304,365,339,306,339,353,351,333)

library(psych)
describe.by(pnt, list(cls, gpa))

fct <- factor(paste(cls, gpa, sep=""))
library(lawstat)
levene.test(pnt, fct)
levene.test(pnt, fct, "mean")
bartlett.test(pnt, fct)

この記事の続き
以下参考文献。
# 帰無仮説有意性検定 (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("http://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

dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9, 10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))

lmres <- lm(result~a, dat)
smry <- summary(lmres)
smry$sigma^2 # 残差の平均平方和
smry$sigma^2 * lmres$df.residual
anova(lmres)

# 被験者間しか使えない