2008年6月アーカイブ

母集団の平均が等しいと仮定できる場合

bigclass.csvを例にとり、男女間で身長の平均に差があるかどうかを検定しましょう。まず、データを読み込みます。

> class <- read.csv("bigclass.csv")

男子と女子のデータをclassから抽出し、それぞれclass_male, class_femaleに代入します。

> class_male <- subset(class,性別=="M")

> class_female <- subset(class,性別=="F")

class_maleの内容を表示してみます。

> class_male

名前 年齢 性別 身長 体重

6 TIM 12 M 60 84

7 JAMES 12 M 61 128

8 ROBERT 12 M 51 79

12 JOHN 13 M 65 98

13 JOE 13 M 63 105

14 MICHAEL 13 M 58 95

準備が出来たので、男女間の身長の差の検定を行いましょう。平均の差の検定にはt.testという関数を利用します。母集団の分散が等しいと仮定しているので、var.equal=Tをオプションとして追加します。また、ここでは両側検定を行いましょう。両側検定を行うにはalternative="two.sided"というオプションを追加します。

> t.test(x=class_male$身長,y=class_female$身長,var.equal=T, alternative="two.sided")

Two Sample t-test

data: class_male$身長 and class_female$身長

t = 2.3687, df = 38, p-value = 0.02304

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

0.4390238 5.6013802

sample estimates:

mean of x mean of y

63.90909 60.88889

p-valueに注目しましょう。0.02304なので、有意水準を5%とすると0.023<0.05より、帰無仮説を棄却し、対立仮説が正しいと判断します。つまり、「男女の母集団の平均身長には差がある」と判断されます。

次に、対立仮説を「男子の平均身長が女子よりも高い」とする片側検定を行いましょう。

今回はalternative="greater"というオプションを追加します。

> t.test(x=class_male$身長,y=class_female$身長,var.equal=T, alternative="greater")

Two Sample t-test

data: class_male$身長 and class_female$身長

t = 2.3687, df = 38, p-value = 0.01152

alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:

0.8705471 Inf

sample estimates:

mean of x mean of y

63.90909 60.88889

p-valueが0.0152なので、5％の有意水準で帰無仮説が棄却され、対立仮説が正しいと判断されます。つまり、「男子の平均身長が女子よりも高い」ということになります。

母集団の分散が等しいと仮定できない場合

は生の落花生と、炒った落花生をそれぞれ10匹のネズミに与えたときのタンパク質の摂取量です。生の落花生といった落花生でタンパク質の平均摂取量が異なるかどうか検定します。

> peanuts <- read.csv("peanuts.csv")

> peanuts

生 炒り

1 61 55

2 60 54

3 56 47

4 63 59

5 56 51

6 63 61

7 59 57

8 56 54

9 44 62

10 61 58

仮説検定にはt.testを利用しますが、今回はvar.equal=Fというオプションを追加します。

> t.test(x=peanuts$生,y=peanuts$炒り,var.equal=F, alternative="two.sided")

Welch Two Sample t-test

data: peanuts$生 and peanuts$炒り

t = 0.9185, df = 17.347, p-value = 0.371

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-2.716617 6.916617

sample estimates:

mean of x mean of y

57.9 55.8

p-valueが0.371なので、有意水準を5%とすると、0.371>0.05より、帰無仮説を棄却できません。つまり、生の落花生と、炒った落花生でも、ネズミが摂取するタンパク質の平均に差はないと判断されます。帰無仮説を「生の落花生の方がタンパク質摂取量の平均値が大きい」とする片側検定は次のようになります。

> t.test(x=peanuts$生,y=peanuts$炒り,var.equal=F, alternative="greater")

Welch Two Sample t-test

data: peanuts$生 and peanuts$炒り

t = 0.9185, df = 17.347, p-value = 0.1855

alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:

-1.872925 Inf

sample estimates:

mean of x mean of y

57.9 55.8

この結果から、有意水準を5％とするならば、片側検定と同じ結果になることが分かります。

2つの標本に対応がある場合

10人の被験者の体温を電子式体温計と水銀式体温計で測った結果が

に保存されています。

> temp <- read.csv("temperature.csv")

> temp

電子式 水銀式

1 37.21 37.02

2 37.12 36.88

3 37.04 37.19

4 37.45 37.27

5 36.93 36.83

6 36.96 36.96

7 37.23 36.92

8 37.15 37.23

9 37.22 37.02

10 37.05 36.78

これは対応のあるデータなので、t.testにpaired=Tというオプションを付加します。

> t.test(x=temp$電子式,y=temp$水銀式, paired=T,alternative="two.sided")

Paired t-test

data: temp$電子式 and temp$水銀式

t = 2.5765, df = 9, p-value = 0.02987

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

0.01537252 0.23662748

sample estimates:

mean of the differences

0.126

p-valueから、5％の有意水準で帰無仮説「2つの種類の体温計で測った体温の平均値に差がない」が棄却され、「平均値に差がある」と判断されます。

片側検定も同様に行えます。

> t.test(x=temp$電子式,y=temp$水銀式, paired=T,alternative="greater")

Paired t-test

data: temp$電子式 and temp$水銀式

t = 2.5765, df = 9, p-value = 0.01493

alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:

0.0363543 Inf

sample estimates:

mean of the differences

0.126

練習問題

1. ある地域には2つのレストランがある。ある1年間の営業日の中から無作為に12日を選び、両レストランの売り上げを記録したデータが

である。このデータから両レストランの平均売上に差があるといえるか検定しなさい。

2. 東京電力と東芝の1959年1月から2001年12月までの月間株価収益率が、

に保存されている。このデータをつかって、2つの銘柄の収益率の平均に差があるかどうかを有意水準5％で検定しなさい。

クロス集計表の作成

まず、データをRに読み込みます。

> car <- read.csv("carsurvey.csv")

データの中身を確認してみましょう。

> car

性別 既未婚 年齢 生産国 サイズ タイプ

1 男性 既婚 34 米国 大型 ファミリー

2 男性 未婚 36 日本 小型 スポーツ

3 男性 既婚 23 日本 小型 ファミリー

4 男性 未婚 29 米国 大型 ファミリー

このデータから、既婚者・未婚者で車の好みが異なるかどうかを検証してみましょう。つまり、「既未婚」の列と、「タイプ」の列を使ってクロス集計表とモザイク図を作成します。

クロス集計表を作成するにはtableという関数を利用します。

> marital_cartype <- table(car$既未婚,car$タイプ)

代入した変数の中身を確認すると、クロス集計表が表示されます。

> marital_cartype

スポーツ ファミリー ワーク

既婚 45 119 32

未婚 55 36 16

モザイク図の作成

次にこのクロス集計表からモザイク図を作成します。Mosaicplotという関数を利用します。

> mosaicplot(marital_cartype,main="既未婚と車の種類")

独立性の仮説検定

結婚の有無と車の選好が独立（無関係）かどうかを仮説検定しましょう。帰無仮説と対立仮説はそれぞれ、

帰無仮説：結婚の有無と車の選好が独立。

帰無仮説：結婚の有無と車の選好が独立ではない。

です。

> chisq.test(marital_cartype)

Pearson's Chi-squared test

data: marital_cartype

X-squared = 26.9629, df = 2, p-value = 1.397e-06

ここで注意しなければならないところは「p-value」です。この値が有意水準よりも小さければ帰無仮説が棄却されます。もし有意水準を一般的な5%にすると1.397e-06<0.05より、帰無仮説が棄却されます(1.397e-06はを意味します)。つまり結婚の有無と車の選好は独立ではないと考えられます。

「X-squared」と「df」はそれぞれ検定統計値と自由度を表します。これらも仮説検定では重要な情報ですが、今回はこれには触れません。

練習問題

「生産国」と「サイズ」、「既未婚」と「サイズ」についてもそれぞれ独立性の検定を行って下さい。

1. 前期（1990年～96年まで）
   

1. 入力Y範囲：B1からB8
2. 入力X範囲：C1からC8
3. 一覧の出力先：A18
4. 「ラベル」にチェック

1. 入力Y範囲：B9からB16
2. 入力X範囲：C9からC16
3. 一覧の出力先：A37
4. 「ラベル」にチェックしない

1. 入力Y範囲：B1からB16
2. 入力X範囲：C1からC16
3. 一覧の出力先：A56
4. 「ラベル」にチェック

1. 入力Y範囲：B1からB11
2. 入力X範囲：C1からD11
3. 「ラベル」にチェック

1. 入力Y範囲：C1からC16
2. 入力X範囲：B1からB16
3. 「ラベル」にチェック

• 残差の標準誤差は、「回帰統計」中の「標準誤差」に表示されています。
• 回帰係数の標準誤差は「切片」および「GDP]と、「標準誤差」が交差するセルに表示されています。
• ｔ値は「切片」および「GDP]と、「ｔ」が交差するセルに表示されています。

1. 分析ツールの「回帰分析」のメニューで以下の設定をした回帰分析を行って下さい。
1. 入力Y範囲：B1からB18
2. 入力X範囲：C1からD18
3. 「ラベル」にチェック
2. さらに、以下の設定をした回帰分析を行って下さい。
1. 入力Y範囲：B1からB18
2. 入力X範囲：C1からE18
3. 「ラベル」にチェック

1. 年齢を横軸に、給与額を縦軸にとる散布図を描いてください。
2. 給与額を被説明変数、給与額とその2乗値を説明変数にする回帰分析を行います。
1. セルC2に「=B2^2」と入力し、C13まで下方向にコピーします。
2. 分析ツールの「回帰分析」を利用し、
1. 入力Y範囲：A1からA13
2. 入力X範囲：B1からC13
3. 「ラベル」にチェック

1. 入力Y範囲：B1からB11
2. 入力X範囲：C1からD11
3. 「ラベル」にチェック

各都道府県のデータを使って、各変数間の散布図の描画と相関係数の計算を行いましょう。まず、データを読み込みます。

> pref <- read.csv("prefectures.csv")

prefと入力すると、prefに代入されたデータが表示されます。

県コ..ド 県名 面積 人口 持家比率 地方交付税 地方税収 商店数 電力消費量 自動車保有

1 1 北海道 83408 5666 54 694 543 63 8623 2065

2 2 青森 9606 1507 72 234 113 22 2128 459

3 3 岩手 15274 1429 73 239 118 21 2088 464

4 4 宮城 7284 2287 61 159 234 30 3468 812

5 5 秋田 11613 1227 80 216 101 20 1733 404

以下では各県の「面積」～「自動車保有」の列だけを利用して分析を進めていくので、必要な列だけを取り出します。

> pref2 <- pref[3:10]

上の命令は変数prefの3列目から10列目をとりだし、変数pref2に代入せよという意味です。

> pref2

面積 人口 持家比率 地方交付税 地方税収 商店数 電力消費量 自動車保有

1 83408 5666 54 694 543 63 8623 2065

2 9606 1507 72 234 113 22 2128 459

3 15274 1429 73 239 118 21 2088 464

4 7284 2287 61 159 234 30 3468 812

5 11613 1227 80 216 101 20 1733 404

各変数間の散布図を作成します。pairsという関数を利用します。

> pairs(pref[3:10],gap=0) #gap=0 は各散布図間の空白を0にする

各変数間の相関係数を計算します。相関係数を計算する関数はcorです。

> cor(pref2)

面積 人口 持家比率 地方交付税 地方税収 商店数 電力消費量 自動車保有

面積 1.00000000 0.1121399 -0.10981581 0.86357637 0.01554832 0.09766679 0.04847636 0.1749191

人口 0.11213987 1.0000000 -0.80391731 -0.15413683 0.72394462 0.97563958 0.99193483 0.9689834

持家比率 -0.10981581 -0.8039173 1.00000000 0.05355635 -0.50281434 -0.81125287 -0.80544938 -0.7503827

地方交付税 0.86357637 -0.1541368 0.05355635 1.00000000 -0.18031202 -0.15466839 -0.22705848 -0.1073727

地方税収 0.01554832 0.7239446 -0.50281434 -0.18031202 1.00000000 0.76653324 0.75725724 0.6397272

商店数 0.09766679 0.9756396 -0.81125287 -0.15466839 0.76653324 1.00000000 0.98409466 0.9264742

電力消費量 0.04847636 0.9919348 -0.80544938 -0.22705848 0.75725724 0.98409466 1.00000000 0.9407229

自動車保有 0.17491910 0.9689834 -0.75038274 -0.10737274 0.63972724 0.92647417 0.94072291 1.0000000

まず、bigclass2.csvをRで読み込みます。

class <- read.csv("bigclass2.csv")

男子と女子のデータをclassから抽出し、それぞれclass_male, class_femaleに代入します。

class_male <- subset(class,性別=="M")

class_female <- subset(class,性別=="F")

class_maleの内容を表示してみます。

class_male

名前 年齢 性別 身長 体重

6 TIM 12 M 60 84

7 JAMES 12 M 61 128

8 ROBERT 12 M 51 79

12 JOHN 13 M 65 98

13 JOE 13 M 63 105

14 MICHAEL 13 M 58 95

男子の身長のヒストグラムと箱ひげ図を表示します。

par(mfrow=c(2,1))

boxplot(class_male$身長,horizontal=T,main="男子身長の箱ひげ図",xlab="身長")

hist(class_male$身長,main="男子身長のヒストグラム",xlab="身長")

同様に女子についても表示してみます。

par(mfrow=c(2,1))

boxplot(class_female$身長,horizontal=T,main="女子身長の箱ひげ図",xlab="身長")

hist(class_female$身長,main="女子身長のヒストグラム",xlab="身長")

基本統計量も見てみましょう。

mean(class_male$身長)

[1] 63.90909

mean(class_female$身長)

[1] 60.88889

sd(class_male$身長)

[1] 4.308453

sd(class_female$身長)

[1] 3.61189

これより男性の身長の方が女性の身長に比べ平均身長が高く、ばらつきも大きいことが分かります。

pref <- read.csv("prefectures.csv")

次に以下の3行を入力します。

par(mfrow=c(2,1)) #グラフ画面を2分割（2行1列）

boxplot(pref$商店数,horizontal=T,main="商店数の箱ひげ図")

hist(pref$商店数,xlab="商店数",main="商店数のヒストグラム")

以下のグラフが出力されます。

ヒストグラムより、分布が右にゆがんでいることがわかります。また、箱ひげ図から外れ値と考えられるものが５つあります。どの件が外れ値になっているのかを、データを商店数でソートすることで調べてみましょう。

pref[order(pref$商店数,decreasing=T),]

県コ..ド 県名 面積 人口 持家比率 地方交付税 地方税収 商店数 電力消費量

13 13 東京 2183 11573 40 0 3860 144 23746

27 27 大阪 1884 8543 48 39 1098 112 16986

23 23 愛知 5147 6715 58 16 921 80 11961

14 14 神奈川 2412 8104 52 33 887 72 13974

28 28 兵庫 8382 5458 60 262 544 70 9623

40 40 福岡 4966 4849 53 250 434 65 8090

この結果から、５つの外れ値は東京、大阪、愛知、神奈川、兵庫だということが分かります。つまり、人口が多い都市にはたくさんの商店があるという、ごく当たり前の結果が現れています。

次に、人口1000人あたりの商店数を計算することで人口規模の要因を除去してヒストグラムを作成しましょう。次の式を入力してください。

pref$千人当たり商店数 <- pref$商店数 / pref$人口 * 1000

これは、読み込んだデータprefの「商店数」の列を「人口」の列で除したうえで1000を乗じ、その結果を「千人あたり商店数」として付け加えることを意味します。

次に「千人あたり商店数」の箱ひげ図とヒストグラムを作成しましょう。

par(mfrow=c(2,1))

boxplot(pref$千人当たり商店数,horizontal=T,main="千人当たり商店数の箱ひげ図")

hist(pref$千人当たり商店数,xlab="千人当たり商店数",main="千人当たり商店数のヒストグラム")

千人あたりの商店数でデータを並び替えてみましょう。

pref[order(pref$千人当たり商店数,decreasing=T),]

県コ..ド 県名 面積 人口 持家比率 地方交付税 地方税収 商店数 電力消費量 自動車保有 千人当たり商店数

39 39 高知 7104 827 67 183 60 15 1447 264 18.137848

16 16 富山 4246 1124 80 135 124 20 1914 456 17.793594

32 32 島根 6626 775 76 181 65 13 1265 257 16.774194

36 36 徳島 4143 837 70 156 73 14 1514 294 16.726404

47 47 沖縄 2264 1267 56 189 72 21 2169 454 16.574586

30 30 和歌山 4722 1095 73 160 97 18 2056 358 16.438356

5 5 秋田 11613 1227 80 216 101 20 1733 404 16.299919

この結果から何が読み取れるでしょうか？考えてみて下さい。

score <- read.csv("classscore.csv")

score

名前 第1回得点 第2回得点

1 KATIE 67 80

2 LOUISE 48 82

3 JANE 53 78

4 JACLYN 60 76

5 LILLIE 29 83

この結果からKATIEの成績は第1回目に比べ第2回目の方が良かったとは言えません。なぜなら、クラス全体の平均得点や点数のばらつきを考慮しなければならないからです。

第1回目と第2回目の得点のヒストグラムを表示してみましょう。

par(mfcol=c(2,1))

hist(score$第1回得点,main="第1回得点",xlab="得点")

hist(score$第2回得点,main="第2回得点",xlab="得点")

第1回目の方がばらつきが大きく、平均点が低いことが分かります。

次に、第1回目の得点と第2回目の得点をそれぞれ基準化します。Rにはデータを基準化するための関数scaleが用意されているので、これを利用します。

score$z1 <- scale(score$第1回得点)

score$z2 <- scale(score$第2回得点)

scoreを表すると基準化されたデータがz1、z2として表示されます。

score

名前 第1回得点 第2回得点 z1 z2

1 KATIE 67 80 1.853667896 0.104549

2 LOUISE 48 82 -0.097561468 0.522745

3 JANE 53 78 0.415919943 -0.313647

4 JACLYN 60 76 1.134793920 -0.731843

5 LILLIE 29 83 -2.048790833 0.731843

この結果から、KATIEは第1回より第2回の成績が悪いことが分かります。

ちなみに基準化得点を10倍し、50を加えたものがいわゆる「偏差値」と呼ばれるものです。偏差値の計算をしてみましょう。

score$偏差値1 <- score$z1 * 10 + 50

score$偏差値2 <- score$z2 * 10 + 50

score

名前 第1回得点 第2回得点 z1 z2 偏差値1 偏差値2

1 KATIE 67 80 1.853667896 0.104549 68.53668 51.04549

2 LOUISE 48 82 -0.097561468 0.522745 49.02439 55.22745

3 JANE 53 78 0.415919943 -0.313647 54.15920 46.86353