今日の授業用のファイルと補足のページです。
作業ディレクトリの設定後、
mouse <- read.table("mouse.csv", sep=",", header=T)
m <- mean(mouse$Weight)
ma <- by(mouse$Weight, mouse$Food, mean)
sst <- sum((mouse$Weight-m)^2) # 総変動（確認のため、実際の分散分析には不要）
ssa <- sum((ma[mouse$Food]-m)^2) # 級（水準）間変動
sse <- sum((mouse$Weight-ma[mouse$Food])^2) # 級（水準）内変動（誤差変動）
pha <- 3-1 # 級（水準）間変動の自由度、水準数-1
pht <- 9-1 # 総変動の自由度、全データ数-1
phe <- pht-pha # 級（水準）内変動の自由度
msa <- ssa/pha # 級間の平均平方
mse <- sse/phe # 級内の平均平方
f0 <- msa/mse # F値
qf(0.01, pha, phe, lower.tail=F) # 棄却域、F値がこれより大きければ帰無仮説は棄却される。
Rでは以前に話したようにベクトルというデータ型を多用します。
petunia <- c(10.3, 11.5, 12.2, 11.3, 10.8, 10.9, 10.5, 10.9, 10.4)
としたとき、petuniaは長さが9の数値型のベクトルとなります。Rではベクトルの演算が簡単にでき、
petunia - 3 とすると、それぞれの値から3を引いたベクトルが返ります。
m <- mean(petunia) としたとき
(petunia - m)^2 で各データの偏差平方がベクトルとして得られ、sum((petunia - m)^2) とすることで偏差平方の和、つまり全変動が得られます。
ベクトルは文字列でも可能で、
species <- c("axillaris", "axillaris", "axillaris", "parodii", "parodii", "parodii", "subandina", "subandina", "subandina") のように作れます。factor関数を使うと文字型ベクトルを因子として扱うことができます。
sp <- factor(species)
> summary(species)
Length Class Mode
9 character character
> summary(sp)
axillaris parodii subandina
3 3 3
のように、違いがわかります。
ma <- by(petunia, sp, mean) としたとき、
maはby型のデータであり、spの水準の情報を内部に持っているため、ma[1]あるいはma["axillaris"]のように値を取り出すことができます。さらに、ma[sp] とするとspの各水準に対する水準の平均がベクトルとして得られます。
> ma[sp] sp axillaris axillaris axillaris parodii parodii parodii subandina subandina subandina 11.33333 11.33333 11.33333 11.00000 11.00000 11.00000 10.60000 10.60000 10.60000
これからmを引いて２乗したものを合計してやれば級間（水準間）変動になります。
ssa <- sum((ma[sp]-m)^2)
また、各データからこれを引いて２乗したものを合計してやれば級内（水準内）変動になります。
sse <- sum((petunia-ma[sp])^2)
これらをそれぞれの自由度で割った値が級間平均平方MSAと級内平均平方MSEであり、MSA ÷ MSEが求める検定統計量F0で、この値が棄却域に入れば帰無仮説が棄却できます。自由度については全変動に対する自由度φTが全データ数 - 1、水準間変動に対する自由度φAが水準の数 - 1、とわかりやすいので、φEは
φT = φA + φE の関係を使って計算するのがいいでしょう。