三要因の分散分析（SABC タイプ；RBFpqr デザイン；被検者内計画） Last modified: Jun 06, 2014
目的
三要因の分散分析（SABC タイプ；RBFpqr デザイン；被検者内計画）を行う。
使用法
SABC(data)
引数
data 4次元配列のデータ（使用例を参照のこと）
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/SABC.R", encoding="euc-jp")
SABC <- function(tbl) {
d <- dim(tbl)
n <- d[4]
r <- d[1]
q <- d[2]
p <- d[3]
X <- sum(tbl)^2/n/p/q/r
A <- sum(apply(tbl, 3, sum)^2)/n/q/r
B <- sum(apply(tbl, 2, sum)^2)/n/p/r
C <- sum(apply(tbl, 1, sum)^2)/n/p/q
S <- sum(apply(tbl, 4, sum)^2)/p/q/r
AB <- sum(apply(tbl, c(3, 2), sum)^2)/n/r
AC <- sum(apply(tbl, c(3, 1), sum)^2)/n/q
BC <- sum(apply(tbl, c(2, 1), sum)^2)/n/p
AS <- sum(apply(tbl, c(4, 3), sum)^2)/q/r
BS <- sum(apply(tbl, c(4, 2), sum)^2)/p/r
CS <- sum(apply(tbl, c(4, 1), sum)^2)/p/q
ABS <- sum(apply(tbl, c(4, 2, 3), sum)^2)/r
ACS <- sum(apply(tbl, c(4, 1, 3), sum)^2)/q
BCS <- sum(apply(tbl, c(4, 2, 1), sum)^2)/p
ABC <- sum(apply(tbl, c(3, 1, 2), sum)^2)/n
ABCS <- sum(tbl^2)
SS.A <- A - X
SS.B <- B - X
SS.C <- C - X
SS.AB <- AB - A - B + X
SS.AC <- AC - A - C + X
SS.BC <- BC - B - C + X
SS.ABC <- ABC - AB - AC - BC + A + B + C - X
SS.T <- ABCS - X
SS.S <- S - X
SS.AS <- AS - A - S + X
SS.BS <- BS - B - S + X
SS.CS <- CS - C - S + X
SS.ABS <- ABS - AB - AS - BS + A + B + S - X
SS.ACS <- ACS - AC - AS - CS + A + C + S - X
SS.BCS <- BCS - BC - BS - CS + B + C + S - X
SS.ABCS <- ABCS - ABC - ABS - ACS - BCS +
AB + AC + BC + AS + BS + CS - A - B - C - S + X
SS <- c(SS.S, SS.A, SS.AS, SS.B, SS.BS, SS.C, SS.CS, SS.AB,
SS.ABS, SS.AC, SS.ACS, SS.BC, SS.BCS, SS.ABC, SS.ABCS, SS.T)
n1 <- n - 1
p1 <- p - 1
q1 <- q - 1
r1 <- r - 1
df <- c(n1, p1, p1 * n1, q1, q1 * n1, r1, r1 * n1, p1 * q1,
p1 * q1 * n1, p1 * r1, p1 * r1 * n1, q1 * r1, q1 * r1 * n1,
p1 * q1 * r1, p1 * q1 * r1 * n1, n * p * q * r - 1)
MS <- SS/df
suf <- 1:7 * 2
p.value <- F.value <- rep(NA, 16)
F.value[suf] <- MS[suf]/MS[suf + 1]
p.value[suf] <- pf(F.value[suf], df[suf], df[suf + 1], lower.tail = FALSE)
ANOVA.table <- data.frame(SS, df, MS, F.value, p.value)
rownames(ANOVA.table) <- c("S", "A", "AxS", "B", "BxS", "C", "CxS",
"AxB", "AxBxS", "AxC", "AxCxS", "BxC", "BxCxS", "AxBxC", "AxBxCxS", "T")
result <- list(ANOVA.table=ANOVA.table, tbl=tbl)
class(result) <- "SABC"
return(result)
}
# ANOVA 表の print メソッド（SABC 関数が返すオブジェクト）
print.SABC <- function(obj) {
x <- obj$ANOVA.table
printf <- function(x, fmt) if (is.na(x)) "" else sprintf(fmt, x)
x[,4] <- sapply(x[,4], printf, "%.5f")
x[,5] <- sapply(x[,5], printf, "%.5f")
print.data.frame(x)
}
# 推定周辺平均を出力する summary メソッド（SABC 関数が返すオブジェクト）
summary.SABC <- function(obj) {
tbl <- obj$tbl
d <- dim(tbl)
n <- d[4]
r <- d[1]
q <- d[2]
p <- d[3]
means.a <- apply(tbl, c(4, 3), sum)/q/r
means.b <- apply(tbl, c(4, 2), sum)/p/r
means.c <- apply(tbl, c(4, 1), sum)/p/q
mean <- c(colMeans(means.a), colMeans(means.b), colMeans(means.c))
SE <- c(apply(means.a, 2, sd),
apply(means.b, 2, sd),
apply(means.c, 2, sd)) / sqrt(n)
t95 <- qt(0.975, n - 1)
EMM <- data.frame(Mean = mean, SE = SE,
LCL = mean - t95 * SE,
UCL = mean + t95 * SE)
rownames(EMM) <- c(paste("a", 1:p, sep="="),
paste("b", 1:q, sep="="),
paste("c", 1:r, sep="="))
print(EMM)
}
使用例
SABCタイプ（被検者内計画）のデータが以下のようにまとめられているとする。
要因 A は 2 水準，要因 B は 3 水準，，要因 C は 4 水準をもち，各被験者は要因 A，要因 B，要因 C のすべての水準の組み合わせ（例では 2 × 3 × 4 = 24 通り）についてデータを採取される。
---------------------------------------------------------------------------------------------------------------
A1 A2
------------------------------------------------- -------------------------------------------------
B1 B2 B3 B1 B2 B3
--------------- --------------- --------------- --------------- --------------- ---------------
C1 C2 C3 C4 C1 C2 C3 C4 C1 C2 C3 C4 C1 C2 C3 C4 C1 C2 C3 C4 C1 C2 C3 C4
--------------- --------------- --------------- --------------- --------------- ---------------
被検者1 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24
被検者2 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34 x35 x36 x37 x38 x39 x40 x41 x42 x43 x44 x45 x46 x47 x48
被検者3 x49 x50 x51 x52 x53 x54 x55 x56 x57 x58 x59 x60 x61 x62 x63 x64 x65 x66 x67 x68 x69 x70 x71 x72
被検者4 x37 x38 x39 x40 x41 x42 x43 x44 x45 x46 x47 x48 x37 x38 x39 x40 x41 x42 x43 x44 x45 x46 x47 x48
被検者5 x49 x50 x51 x52 x53 x54 x55 x56 x57 x58 x59 x60 x61 x62 x63 x64 x65 x66 x67 x68 x69 x70 x71 x72
被検者6 x73 x74 x75 x76 x77 x78 x79 x80 x81 x82 x83 x84 x85 x86 x87 x88 x89 x90 x91 x92 x93 x94 x95 x96
被験者7 x97 x98 …
:
---------------------------------------------------------------------------------------------------------------
以上のようなデータを，R の配列として以下のように付値する。
配列名を data とすると
data <- array(c(x1, x2, x3, ..., x96, ...), dim=c(4, 3, 2, 7))
dim=c(4, 3, 2, 7) の3つの数値は順に，要因 C の水準数，要因 B の水準数，要因 A の水準数，被験者数。
出力結果例
# 森敏昭，吉田寿夫「心理学のためのデータ解析テクニカルブック」北大路書房，P. 152 の例
> dat <- c(
+ 2,5,9, 6,3,6, 1,2,1, 5,5,5,
+ 6,7,10, 6,6,7, 2,1,5, 3,6,5,
+ 5,9,13, 8,5,5, 5,4,3, 4,6,9,
+ 7,9,14, 10,8,6, 2,5,5, 6,7,7)
> tbl <- array(dat, dim = c(3, 2, 2, 4))
> ans <- SABC(tbl)
> ans # 分散分析表
SS df MS F.value p.value
S 60.333333 3 20.1111111
A 96.333333 1 96.3333333 51.00000 0.00565
AxS 5.666667 3 1.8888889
B 3.000000 1 3.0000000 1.42105 0.31893
BxS 6.333333 3 2.1111111
C 33.500000 2 16.7500000 60.30000 0.00011
CxS 1.666667 6 0.2777778
AxB 56.333333 1 56.3333333 101.40000 0.00209
AxBxS 1.666667 3 0.5555556
AxC 6.166667 2 3.0833333 2.92105 0.13007
AxCxS 6.333333 6 1.0555556
BxC 24.500000 2 12.2500000 5.37805 0.04591
BxCxS 13.666667 6 2.2777778
AxBxC 41.166667 2 20.5833333 6.07377 0.03614
AxBxCxS 20.333333 6 3.3888889
T 377.000000 47 8.0212766
# print.SABC 関数がないときには空白欄に NA が表示されることがある
> summary(ans) # 推定周辺平均
Mean SE LCL UCL
a=1 7.166667 0.7905694 4.650722 9.682611
a=2 4.333333 0.5400617 2.614616 6.052051
b=1 5.500000 0.8193267 2.892537 8.107463
b=2 6.000000 0.5046084 4.394111 7.605889
c=1 4.875000 0.6166104 2.912671 6.837329
c=2 5.500000 0.7430231 3.135369 7.864631
c=3 6.875000 0.5994789 4.967190 8.782810
直前のページへ戻る E-mail to Shigenobu AOKI
