共分散分析 Last modified: Jun 29, 2004
目的
共分散分析を行う
使用法
covar.test(dat, cp1, cp2, cp3)
引数
dat データ行列（行がケース，列が変数）
cp1 独立変数の列番号
cp2 従属変数の列番号
cp3 群変数の列番号
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/covar_test.R", encoding="euc-jp")
# 共分散分析
covar.test <- function( dat, # データ行列
cp1, # 独立変数の列番号
cp2, # 従属変数の列番号
cp3) # 群変数の列番号
{
dat <- subset(dat, complete.cases(dat[,c(cp1, cp2, cp3)])) # 欠損値を持つケースを除く
x <- dat[,cp1] # 独立変数
y <- dat[,cp2] # 従属変数
g <- dat[,cp3] # 群変数
nj <- table(g) # 各群の例数
n <- sum(nj) # 全例数
k <- length(nj) # 群の個数
# 独立変数について
mx <- mean(x) # 全体の平均値
mxj <- tapply(x, g, mean) # 各群の平均値
sstx <- (n-1)*var(x) # 全変動
ssbx <- sum(nj*(mxj-mx)^2) # 群間変動
sswx <- sstx-ssbx # 群内変動
# 従属変数について
my <- mean(y) # 全体の平均値
myj <- tapply(y, g, mean) # 各群の平均値
ssty <- (n-1)*var(y) # 全変動
ssby <- sum(nj*(myj-my)^2) # 群間変動
sswy <- ssty-ssby # 群内変動
spt <- (n-1)*cov(x, y) # 共変動
spb <- sum(nj*(mxj-mx)*(myj-my)) # 群間共変動
spw <- spt-spb # 群内共変動
ss.wy <- sswy-spw^2/sswx # 全群に共通な傾きと各群ごとの切片を持つ回帰直線からの変動
ss.ty <- ssty-spt^2/sstx # 全データに基づく回帰直線からの変動
ss.by <- ss.ty-ss.wy # 各群の回帰の差に起因する変動
hensa.x <- x-mxj[g] # 各群の平均値からの偏差
hensa.y <- y-myj[g] # 各群の平均値からの偏差
xy <- hensa.x*hensa.y
xx <- hensa.x^2
numerator <- tapply(xy, g, sum)
denominator <- tapply(xx, g, sum)
a <- numerator/denominator # 各群のデータに基づく回帰直線の傾き
b <- myj-a*mxj # 各群のデータに基づく回帰直線の切片
predict.y <- a[g]*x+b[g] # 各群ごとのデータに基づく回帰直線による予測値
ss.wyj <- sum((y-predict.y)^2) # 各群ごとのデータに基づく回帰直線からの変動
df.r <- k-1
df.e <- n-2*k
df.t <- n-k-1
diff.reg <- ss.wy-ss.wyj # 各群の回帰の差による推定誤差平方和
ms.r <- diff.reg/df.r # 各群の回帰の差による推定誤差分散
ms.e <- ss.wyj/df.e # 各群の推定誤差の和による推定誤差分散
ms.w <- ss.wy/df.t # 平均回帰に基づく推定誤差による推定誤差分散
f <- ms.r/ms.e # 検定統計量
p <- pf(f, df.r, df.e, lower.tail=FALSE) # P 値
anova.table <- data.frame( # 分散分析表としてまとめる
"SS"=c(diff.reg, ss.wyj, ss.wy),
"d.f."=c(df.r, df.e, df.t),
"MS"=c(ms.r, ms.e, ms.w)
)
rownames(anova.table) <- c("group x slope", "residual", "total")# 行の名前
test.result <- c("F value"=f, "d.f.1"=df.r, "d.f.2"=df.e, "P value"=p)
if (p <= 0.05) { # 各群における回帰直線が等しくないとき，
part2 <- anova.table2 <- test.result2 <- NULL # それ以上の検定は行わない
}
else { # 各群における回帰直線の傾きが等しくないとはいえないとき
part2 <- "H0: 共変量で調整した平均値は同じである" # 帰無仮説
ms.by <- ss.by/df.r
ms.wy <- ss.wy/df.t
anova.table2 <- data.frame( # 分散分析表としてまとめる
"SS"=c(ss.by, ss.wy, ss.ty),
"d.f."=c(df.r, df.t, n-2),
"MS"=c(ms.by, ms.wy, ss.ty/(n-2))
)
rownames(anova.table2) <- c("effect & group", "residual", "total")# 行の名前
f2 <- ms.by/ms.wy # 検定統計量
p2 <- pf(f2, df.r, df.t, lower.tail=FALSE) # P 値
test.result2 <- c("F value"=f2, "d.f.1"=df.r, "d.f.2"=df.t, "P value"=p2)
}
return(list( part1="H0: 各群の回帰直線の傾きは同じである",
result1.1=anova.table,
result1.2=test.result,
part2=part2,
result2.1=anova.table2,
result2.2=test.result2))
}
使用例
dat <- matrix(c( # データ行列例（ファイルから読んでも良い）
5.9, 1.87, 1, # 1列目が独立変数，2列目が従属変数，3列目が群を表す変数
5.7, 1.19, 1,
5.7, 1.22, 1,
6.7, 1.46, 1,
6.5, 1.35, 1,
6.2, 1.16, 1,
6.3, 1.62, 1,
6.7, 2.00, 1,
5.6, 1.28, 2,
5.3, 1.06, 2,
5.5, 1.17, 2,
5.6, 1.74, 2,
5.5, 1.13, 2,
5.3, 1.18, 2,
5.2, 1.03, 2,
5.6, 1.23, 2,
5.5, 0.90, 2,
5.5, 1.24, 2
), ncol=3, byrow=TRUE)
covar.test(dat, 1, 2, 3)
出力結果例
$part1
[1] "H0: All slopes are equal."
$result1.1
SS d.f. MS
group x slope 0.03576607 1 0.03576607
residual 0.91878786 14 0.06562770
total 0.95455392 15 0.06363693
$result1.2
F value d.f.1 d.f.2 P value
0.5449843 1.0000000 14.0000000 0.4725688
$part2
[1] "H0: Adjusted means are equal"
$result2.1
SS d.f. MS
effect & group 7.909420e-06 1 7.909420e-06
residual 9.545539e-01 15 6.363693e-02
total 9.545618e-01 16 5.966011e-02
$result2.2
F value d.f.1 d.f.2 P value
1.242898e-04 1.000000e+00 1.500000e+01 9.912519e-01
解説ページ
直前のページへ戻る E-mail to Shigenobu AOKI
