myscratchpad

by関数とAnovaのerror引数を使ってプールした誤差項による検定を行う
こちらのページを参考にした。感謝。卒論でやったという。最近の卒論はレベルが高いのう
データはテクニカルブックp95
2要因とも対応がない場合

# データ入力
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))

# 合計を出す
aggregate(dat[4], list(dat[,3], dat[,2]), sum)

# lmとcarのAnovaで分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(result~a*b, data=dat)
library(car)
Anovares <- Anova(lmres, type=3)

# a要因についてプールした誤差項で検定
by(dat, dat$a, function(x) Anova(lm(result~b, data=x), type=2, error=lmres))
# b要因についてプールした誤差項で検定
by(dat, dat$b, function(x) Anova(lm(result~a, data=x), type=2, error=lmres))
## セルサイズが等しいため、テクニカルブックp101どおりの結果になる
## by関数がこんなに便利とは、Anovaで誤差項が指定できるとは知らなかった

## 問題
上記の分析で、type=3とするとプールした誤差項を使わず、そのまんま水準別で検定したのと同じ結果になる というか、linear.hypothesisでerror.SSやerror.dfの情報を使っていないような…
car:::linear.hypothesis.lm
car:::linear.hypothesis.default

…分散分析ばっかりふえるのう…

#　対応のない1要因
# データ生成。テクニカルブックp. 87
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))

# 平均と合計
aggregate(dat[2], list(dat[,1]), mean)
aggregate(dat[2], list(dat[,1]), sum)

# lmで分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(result~a, dat)
# carのAnovaで分析
library(car)
Anovares <- Anova(lmres, type=3)
Anovares

# aovと比較
summary(aov(result~a, dat))


#　対応のある1要因
# データ生成。テクニカルブックp. 92
dat <- data.frame(s=factor(1:8), 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))

# 平均と合計
sapply(dat[2:5], mean)
sapply(dat[2:5], sum)

# 分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(a1, a2, a3, a4)~1, dat)
# carのAnovaで分析
library(car)
afact <- factor(c("a1", "a2", "a3", "a4"))
idat <- data.frame(afact)
Anovares <- Anova(lmres, idata=idat, idesign=~afact, type=3)
summary(Anovares)
(AnovaTable <- Anovastat(Anovares))

# aovと比較
library(reshape)
dat2 <- melt(dat, id="s")
names(dat2) <- c("s", "a", "result")
summary(aov(result~a+s+Error(a:s), dat2))

carパッケージ、Anova関数の分散分析表とかはprintで結果を出力しているだけ。
データフレームとか行列になってると、後からF値や自由度を取り出して効果量なんかを出しやすいので関数 "Anovastat" を作ってみた。
少し改良して"Astat"にした。
といっても、car:::summary.Anova.mlmやcar:::print.linear.hypothesis.mlmを適当に書き換えただけ。Rjpwikiの記事に感謝。

## 関数の読み込み
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1264349365")
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1273264153")
# 我ながらきたないコードだなぁ…

## 使い方

# 分析サンプル
library(car)
hsb2 <- read.table('http://www.ats.ucla.edu/stat/R/faq/hsb2.csv', header=T, sep=",")
result1 <- Anova(lm(cbind(write,read)~female+math+science+socst, hsb2))


## Anovastatで結果の取り出し
Tableres <- Anovastat(result1)
Tableres <- Astat(result1)
Tableres
# manovaはWilksをデフォルトで出力する。単変量Anovaには偏イータ二乗、f二乗も出力する。
Tableres$Anova # 単変量分散分析の結果。混合要因のときのみ出力される (上の分析例では出力されない) 。通常はUAnova (単変量分散分析表) 、Mauchly, GG(Greenhouse-Geisser) 、HF (Huynh-Feldt) がmatrixで出力される。
# Tableres$Anova$UAnovaとかで取り出す
Tableres$Manova # 多変量分散分析の結果
Tableres$Manova[[2]] # とかやれば個々の分散分析表がデータフレームで取り出せる

## 被験者間要因のみのときはそもそもデータフレームなので上の関数は必要ない
result2 <- Anova(lm(write~female*math, hsb2))
class(result2)

各水準間の差得点を条件間で比較するだけ
入戸野先生のホームページの情報、pdfが大変分かりやすかった。感謝

# データの準備。データ解析テクニカルブックのp. 116より
dat <- data.frame(
s=factor(c("s1", "s2", "s3", "s4", "s5")),
a1b1=c(3, 3, 1, 3, 5), a1b2=c(4, 3, 4, 5, 7), a1b3=c(6, 6, 6, 4, 8), a1b4=c(5, 7, 8, 7, 9),
a2b1=c(3, 5, 2, 4, 6), a2b2=c(2, 6, 3, 6, 4), a2b3=c(3, 2, 3, 6, 5), a2b4=c(2, 3, 3, 4, 6)
)
dat
sapply(dat[2:9], sum)
sapply(dat[2:9], mean)

options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(a1b1,a1b2,a1b3,a1b4,a2b1,a2b2,a2b3,a2b4)~1, dat)
library(car)
afact <- factor(c(rep("a1", 4), rep("a2", 4)))
bfact <- factor(c(rep(c("b1","b2","b3","b4"),2)))
idat <- data.frame(afact, bfact)
Anovares <- Anova(lmres, idata=idat, idesign=~afact*bfact, type=3)
summary(Anovares)

# 交互作用対比の準備
data1 <- dat[2:5]
data2 <- dat[6:9]
commat <- combn(x=1:4, 2) # combn関数で組み合わせをつくる
combnN <- ncol(commat)
c1 <- commat[1,]
c2 <- commat[2,]

# t検定結果格納用のオブジェクトを用意する
tres <- data.frame(matrix(1, 6, 9))
colnames(tres) <- c("d1mean", "d2mean", "d1sd", "d2sd", "t.value", "df", "p.value", "CIlower", "CIupper")
tresam <- data.frame(matrix(1,6,2)) # 文字列用
colnames(tresam) <- c("alternative", "method")
rwnames <- letters[1:6] # 行名用
# t検定繰り返し
for (i in 1:combnN) {
d11 <- data1[c1[i]]
d12 <- data1[c2[i]]
d21 <- data2[c1[i]]
d22 <- data2[c2[i]]
d1 <- d11[,1]-d12[,1]
d2 <- d21[,1]-d22[,1]
tobject <- t.test(d1, d2, alternative="two.sided", paired=T)
tres.vct <- c(mean(d1), mean(d2), sd(d1), sd(d2), tobject[[1]], tobject[[2]], tobject[[3]], tobject[[4]][1], tobject[[4]][2])
tres[i,] <- tres.vct
tresam[i,] <- c(tobject[[7]], tobject[[8]])
cntname <- paste(names(d11), "-", names(d12), " vs. ", names(d21), "-", names(d22), sep="")
rwnames[i] <- cntname
}

rownames(tres) <- rwnames
tres <- data.frame(tres, tresam)
# t検定結果
tres
# holm法でp値の調整
(p.holm <- p.adjust(tres["p.value"][,1], "holm"))
(tresp <- data.frame(tres, p.holm=p.holm))

## 入戸野先生の解説は有意水準を=.05として、これをp値の順位で割って有意水準自体を調整している。
hprank <- c(3,5,6,1,4,2) # 差の小ささ。最も差が小さいのは3番目ということ
hp <- rep(0.05, 6)
hp/hprank # これとtres["p.value"][,1]を比較する
(tresp2 <- data.frame(tres, p.holm=p.holm, p.ntn=hp/hprank))
## なので、有意水準を.05とし、tres["p.value"][,1]を順位倍して.05を超えるかどうか見ても同じ。p.adjustでやってるのはこれ。ただし、第n+1ステップが第nステップの有意水準より小さくならないように調整する。小さくなった場合はn+1の水準を使う


## 別のやり方
dat2 <- dat[c(2,5,6,9)]
options()$contrasts

lmres2 <- lm(cbind(a1b1, a1b4, a2b1, a2b4)~1, dat2)
afact2 <- factor(c("a1", "a1", "a2", "a2"))
bfact2 <- factor(c("b1", "b4", "b1", "b4"))
idat2 <- data.frame(afact2, bfact2)
idat2
library(car)
Anovares2 <- Anova(lmres2, idata=idat2, idesign=~afact2*bfact2, type=3)
Ares2 <- Astat(Anovares2)
Ares2$unv
tres[3,]

## グラフにしてみる
pmat <- matrix(sapply(dat2, mean), 2)
dimnames(pmat) <- list(c("b1", "b4"), c("a1", "a2"))

palette(grey.colors(2))
barplot(pmat, beside=T)
legend("top", rownames(pmat), fill = 1:2)
win.graph()
barplot(t(pmat), beside=T)
legend("top", colnames(pmat), fill = 1:2)

#aovの結果
Error: a:s
          Df  Sum Sq Mean Sq F value    Pr(>F)   
a          3 217.500  72.500 21.4437 1.347e-06 ***
s          7  77.000  11.000  3.2535   0.01682 * 
Residuals 21  71.000   3.381 

# Anovaの結果
                SS num Df Error SS den Df         F       Pr(>F)
(Intercept) 3444.5      1       77      7 313.13636 4.537639e-07
afact        217.5      3       71     21  21.44366 1.347309e-06

Anova関数の被験者要因は
平方和: (Intercept) のError SS = 77
平均平方和: (Intercept) のError SS/ (Intercept) のden DF 77/7 = 11
誤差平方和: 要因のError SS/要因のden Df = 71/21 = 3.380952
F値: 平均平方和/誤差平方和 = 11/3.380952 = 3.253521
要因自由度: Interceptのden Df = 7
誤差自由度: 要因のden Df = 21


## aov結果のとりだし
aovrestb <- summary(aov(result ~ a + s + Error(a:s), RB))
## aovrestb[[1]][[1]]["Sum Sq"]