今日の授業用のファイルと補足のページです。
作業ディレクトリの設定後、
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を引いて2乗したものを合計してやれば級間(水準間)変動になります。
ssa <- sum((ma[sp]-m)^2)
また、各データからこれを引いて2乗したものを合計してやれば級内(水準内)変動になります。
sse <- sum((petunia-ma[sp])^2)
これらをそれぞれの自由度で割った値が級間平均平方MSAと級内平均平方MSEであり、MSA ÷ MSEが求める検定統計量F0で、この値が棄却域に入れば帰無仮説が棄却できます。自由度については全変動に対する自由度φTが全データ数 - 1、水準間変動に対する自由度φAが水準の数 - 1、とわかりやすいので、φEは
φT = φA + φE の関係を使って計算するのがいいでしょう。
トウモロコシの草丈 (cm) に関して二つの処理区より下記のデータを得た。A区とB区の草丈に違いがあると言えるか。Rを用いて5%水準で実際にt検定を行いレポートとして提出して下さい。
A区 308, 319, 313, 304, 296, 302, 307
B区 236, 260, 245, 275, 251, 245
Rについては、t.test関数を使わずにRプログラムの授業用解説ページのt検定の項目を参考にして偏差平方和とサンプル数からt0の値を計算し、qt関数を用いて5%水準の棄却限界値を求め、帰無仮説が棄却できるかどうかを判断して下さい。
レポートの形式はワードまたは適当なテキストエディタを使ってRの出力をコピーペーストし、そこに必要な文章を書き込んで下さい。
提出期限: 2016年11月22日(火)
提出先: Moodleまたは電子メールでhkokubun@faculty.chiba-u.jpまで
サラブレッド種とアラブ種の体高の平均値に差があるか、t検定によって検定する。
帰無仮説 H0: サラブレッドとアラブの体高に差はない
対立仮説 H1: サラブレッドとアラブの体高に差はある
それぞれのデータを変数に付値する。
> sara<-c(160,168,158,165,161) > arab<-c(153,155,149,152)
t0を計算する。
> m1<-mean(sara) > m2<-mean(arab) > ss1<-sum((sara-m1)^2) > ss2<-sum((arab-m2)^2) ・・・
自由度 (n1-1)+(n2-1) = 7の有意水準5%の両側検定における棄却限界値を求める。
> qt(0.025,7,lower.tail=F) [1] 2.364624
t0の値が棄却域に入るので帰無仮説H0は棄却できる。従って、サラブレッド種とアラブ種の体高の平均は有意水準5%で有意に差があると言える。