三要因の分散分析(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 <- 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         125.89166   4  31.472915                
A         318.37176   1 318.371763 3.15076 0.15055
AxS       404.18356   4 101.045890                
B         179.77618   2  89.888091 0.92205 0.43617
BxS       779.89719   8  97.487149                
C          18.80278   3   6.267594 0.05300 0.98313
CxS      1419.11659  12 118.259716                
AxB        15.29180   2   7.645901 0.07374 0.92954
AxBxS     829.53652   8 103.692065                
AxC       424.39788   3 141.465959 1.73303 0.21336
AxCxS     979.55093  12  81.629244                
BxC      1124.59561   6 187.432602 2.24433 0.07347
BxCxS    2004.33581  24  83.513992                
AxBxC    1271.00070   6 211.833450 1.71934 0.15964
AxBxCxS  2956.95344  24 123.206393                
T       12851.70244 119 107.997499                

# print.SABC 関数がないときには空白欄に NA が表示されることがある

> summary(ans) # 推定周辺平均
        Mean        SE      LCL      UCL
a=1 52.03100 1.2835538 48.46728 55.59472
a=2 48.77333 0.7490903 46.69353 50.85314
b=1 49.67875 1.5267199 45.43990 53.91760
b=2 52.12575 1.6794315 47.46290 56.78860
b=3 49.40200 0.7140144 47.41958 51.38442
c=1 50.97000 1.7611640 46.08022 55.85978
c=2 49.93933 2.2108857 43.80093 56.07774
c=3 50.14400 1.3111927 46.50355 53.78445
c=4 50.55533 1.7793617 45.61503 55.49563


・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI

Made with Macintosh