生物統計学

11月11日の授業について

今日の授業用のファイルと補足のページです。

マウスのデータ(エクセル形式)

マウスのデータ(CSV形式)

マウスのデータで一元配置の分散分析を行う手順

作業ディレクトリの設定後、

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%で有意に差があると言える。