Rでのコントラストについてのメモをいくつかまとめる(未整理)
y = u + b1 x1 + b2 x2 + b3 x3 + b4 x4
というモデルが想定される場合,五つのパラメータを想定する必要はない。 treatment mean は b1 を 0 と見なし,x1 の平均を u とし,他の b は,u との差と考える。
これに対して Helmert Contrasts は,
全体平均を最初のパラメータとし 最初と二つ目の平均の平均と,最初の平均との差を第二のパラメータ 最初と二つ目の平均の平均と,最初と二つ目さらに三つ目の平均の平均との差を第三のパラメータ 最初と二つ目さらに三つ目の平均の平均と,全体平均の差を第四のパラメータ
Brian S. Everitt, Torsten Hothorn: A Handbook of Statistical Analyses using R
# p.75 下
# treatment contrasts を理解する
comp <- read.table( "http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/competition.txt", header = TRUE) attach(comp) comp.aov <- aov(biomass ~ clipping) summary(comp.aov ) summary.lm(comp.aov)
# 上の出力の内容を確認
(means <- tapply(biomass, clipping, mean)) means - means[1] # control 群の平均を引く.これが「効果」
# 標準誤差 Std.Err の意味
# Intercept の誤差は,control 群の平均の誤差
# プールされた平均誤差分散を control の繰り返し数で割って平方根を取る
sqrt(4961/6)
# 2行目以降の標準誤差は平均の「差」の誤差
# 水準が独立であれば,平均の差の分散は,二つの水準の分散の和
sqrt(2 * 4916/6) detach(comp)
# Fox p.143 より
Baumann <- read.table("http://150.59.18.68/Baumann.txt") attach(Baumann) baumann.aov1 <- aov(post.test.3 ~ group) summary.lm(baumann.aov1) # 通常のコントラストによる解析 # これは Basel と他の水準の差 contrasts(group) <- matrix(c(1,-0.5,-0.5, 0,-1,1), 3,2) contrasts(group) # 直交するコントラストを作成 baumann.aov2 <- aov(post.test.3 ~ group) summary.lm(baumann.aov2)
# 標準的方法と最新の方法の間には差があるが
# 最新の二つの方法の間には差は無い