忍者ブログ
×

[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。


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とするとプールした誤差項を使わず、そのまんま水準別で検定したのと同じ結果になる
Anova(aov(result~b, data=dat, subset=(a=="a1")), type=3, error=(lmres))
car:::Anova.III.lm, car:::Anova.II.lmを見比べてみると、タイプ3のときのAnova.III.lm指定した誤差項の情報(error.df, error.SS) がlinear.hypothesisの中で使われていないようだ。
というか、linear.hypothesisでerror.SSやerror.dfの情報を使っていないような…
car:::linear.hypothesis.lm
car:::linear.hypothesis.default

…分散分析ばっかりふえるのう…
PR
Comment
Trackback
Trackback URL

Comment form
Title
Color & Icon Vodafone絵文字 i-mode絵文字 Ezweb絵文字  
Comment
Name
Mail
URL
Password
プロフィール
HN:
tao
HP:
性別:
非公開
職業:
会社員
趣味:
アウトドア、自転車、ジョギング、英語学習
自己紹介:
・千葉在住のサラリーマンです。データ分析っぽいことが仕事。
・今年英検1級取得。今はTOEIC高得点を目指して勉強中。
・興味のあることは野球、アウトドア、英語学習、統計、プログラミング、PC関係などなど。
ブログ内検索
freead
順位表
プロ野球データFreak
セリーグ順位表
パリーグ順位表