基礎統計学（明治大学グローバルビジネス研究科） 福地純一郎のページ

基礎統計学 講義ノート

2012

年

5

月

29

日バージョン

福地純一郎

[email protected]

第

1

章

統計学とは

自然現象、社会現象は多数の要因が複雑にからみあい、さまざまな出来事が生じています。

私たちは、これらの現象を客観的に理解するために、いろいろなデータを収集し分析していま

す。この授業では、統計学の基礎を理解し、Rを用いて実際にデータを分析することによって、

統計学の基礎を身につけます。各自の研究でデータ解析ができるようになることを目標とする。

現代の統計学は確率論という数学が基礎にあり、統計学の理論を深く理解するためには、一

定程度の数学の知識が必要で、多くの時間勉強する必要があります。しかしながら、この授業

では統計学の理論の説明は最小限にとどめ、実際に統計学を使ってみることを目的とします。

さまざまなデータを分析し、統計学が広い範囲で役にたつことを実感するでしょう。

1.1

データの種類

測定の尺度からデータの種類を分類してみる。実際に多く使われるデータは、比例尺度、順

序尺度、名義尺度である。比例尺度とは、長さ、重さ、広さ、金額などのような量を測る尺度

であり、比例尺度で表される変数は連続的な値をとるから連続変数である。順序尺度は成績の

「不可、可、良、優」のように順序に意味のあるような尺度である。名義尺度は、性別「男、女」

や、産業の分類「第一次、第二次、第三次産業」などのように複数カテゴリーに分割されるが、

本来数値的意味はない尺度である。順序尺度や名義尺度で表される変数をカテゴリカル変数と

いい、そのデータをカテゴリカルデータという。経済、経営では、連続データもカテゴリカル

データのどちらも重要である。この講義ノートでは、主に第3章から第6章まで連続データの

分析方法を学び、第7章でカテゴリカルデータの分析方法を学ぶ。

本講義ノート作成のために参考にした文献

「The R Tips」、舟尾暢男著、2009年、オーム社

「Rによるやさしい統計学」、山田剛史・杉澤武俊・村井潤一郎著、2008年、オーム社

第

2

章

R

の基礎

Rはデータ解析および統計分析用のフリーのソフトウエアです。世界中の人々がRの発展に

貢献しています。

2.1

R

の起動、四則演算

Rのアイコン（青いマーク）をダブルクリックし、Rを起動します。Rの操作は、プロンプ

ト（>）の右にRのコマンドを入力することによって行います。

足し算、引き算、掛け算、割り算はそれぞれ

+, -, *, /

を用います。たとえば2 + 4を計算するときは、

> 2+4

と入力し、[Enter]キーを押します。すると以下のように結果が得られます。

> 2+4

[1] 6

以下では、8−3, 3×4, 12÷3 を求めています。

> 8-3

[1] 5

> 3*4 [1] 12

> 12/3

[1] 4

平方根を求めるときは

sqrt()

を使います。

> sqrt(2)

[1] 1.414214

また、累乗を計算するときは ˆ を使います。たとえば、42を求めるときは

> 4ˆ2

[1] 16

とします。

2.2

R

の終了

Rを終了するときは、

> q()

とします。すると、作業スペースを保存しますか？と聞いてくるので、「いいえ」をクリックし

てください。

2.3

以前に計算した式を呼び出す

Rでは、一度計算した式を呼び出すことができます。キーボードの上矢印[↑]キーを押すと、

起動してから現在まで計算したコマンドを呼び出すことができます。

2.4

代入

前節では、Rによる四則演算の方法を学びました。実際に分析を行うときには、計算結果を

どこかに保存しておくと便利です。Rでは、計算結果を変数に代入することができます。たと

えば、sqrt(3) の結果を x という変数に代入するためには以下のコマンドを入力します。

> x<-sqrt(3)

x に入っている値を確認するためには、単に変数の名前を書けばよい。

> x

[1] 1.732051

課題 2.2

√

5の値を変数xに代入せよ。次にxの値を確認せよ。

オブジェクトの削除

Rで作業をしていると、すでに保存したオブジェクトが不要になることがある。以下のコマ

ンドで、作業スペースのオブジェクトを削除する。

> rm(list=ls())

2.5

スクリプトファイルによる実行

準備

USBメモリにrdataという名前のフォルダを作成してください。この授業で使うデータファ

イルはこのフォルダに保存します。https://sites.google.com/site/junfukuchistat/ から二つの

スクリプトファイルによるコマンドの実行

今までは、コマンドの実行を、直接プロンプトの右側に打ち込んで行ってきた。この方法だ

と、毎回コマンドを打ち込む必要があるし、入力したコマンドに誤りがあったときには修正が

しにくい。

実際のデータ分析をするときには、スクリプトファイルを使うと大変便利である。

スクリプトファイルによるコマンドの実行の手順

１．Rのメニューの「ファイル」の中の「新しいスクリプト」をクリック

２．「Rエディタ」というウインドウにRのコマンドを入力する。

３．入力したRのコマンドを反転させる（左クリックしてドラッグする）。

４．右クリックして、「カーソル行あるいは選択中のRコードを実行」をクリックするとRの

コマンドが実行される。

５．一連のRのコマンドを入力してから、最後に名前を付けてフォルダrdata に保存する。

この方法を使えば、利用したRのコードを保存できるし、一連のRのコマンドを一気に実行

することもできる。

課題 2.3 新しいスクリプトを開き以下のRコードを入力しなさい。

x<-sqrt(5) 5の平方根を求める

x

次に、この2行のコマンドを実行しなさい。

最後に、このスクリプトファイルを chap2.r という名前で保存せよ。

2.6節以降の2章のRのコードはすべて、chap2.rに入力し、実行してください。また、3

章以降も、Rのコードはすべてスクリプトに入力し、chap3.r, chap4.rなどのように章ご

とに名前を付けて、フォルダrdataに保存してください。

2.6

R

によるデータの扱い

実際のデータを扱う方法を説明する。最も単純な方法は、直接データを入力する方法である。

関数 cは、複数の数値を合わせてベクトルを作成する。たとえば、高校生３人のお小遣いの

データ6000,8000,12000をベクトル kozukaiとして定義するには以下のようなコマンドを入

力すればよい。kozukaiを定義した後、kozukaiと入力し、kozukaiの中身を表示させる。

> kozukai<-c(6000,8000,12000)

> kozukai

[1] 6000 8000 12000

Rでは、定義したベクトルに対して関数を適用し、統計分析をすることができる。お小遣いデー

タの平均を求めてみよう。データの平均を求めるRの関数はmeanである。

> mean(kozukai)

3人のお小遣いの平均が8666.7円であることがわかった。このように関数を多用することがR

の1つの特徴である。

ベクトルについての計算も簡単に行えます。ベクトルに対して、+, -, *, /を用いると要素ご

とにその計算を行います。

> kozukai

[1] 6000 8000 12000

> kozukai+1000

[1] 7000 9000 13000

> kozukai/2

[1] 3000 4000 6000

> kozukai

[1] 6000 8000 12000

> kozukai*2

[1] 12000 16000 24000

> kozukai/3

[1] 2000.000 2666.667 4000.000

> kozukai/2

[1] 3000 4000 6000

2.7

csv

ファイルからのデータの読み込み

csv

ファイルとは

直接データを入力する方法は単純であるが、データが大量にあるときには、不便である。通

常は、データをcsv ファイルという形式で保存しておき、そのファイルからRにデータを読み

込むこんだ後に分析を行う。

CSV形式とは、Comma Separated Value format の頭文字をとったもので、データをカンマ

「,」や改行で区切って並べたテキスト形式のことをいう。CSV形式は、多くのソフトウエアか

ら読み込むことができて便利です。

例として car.csv をエクセルで開くと以下のような形式になっている。

id age income maker

1 29 500 トヨタ

2 36 600 日産

3 32 400 トヨタ

4 45 1000 ホンダ

5 44 1200 ジャガー

6 47 800 ホンダ

7 28 400 なし

8 37 300 なし

9 37 600 ダイハツ

10 45 900 日産

RでCSVファイルを読む込むときは、car.csvのように横の並び1つ1つが個体、縦のなら

R

で

CSV

ファイルを読み込む

RでCSVファイルを読み込むのは関数 read.csv を使って簡単に行える。

> car<-read.csv("D:/rdata/car.csv")

上のコマンドで、car.csvのデータがオブジェクトcarに保存された。オブジェクト名を入

力すると、中身が表示される。

> car

id age income maker

1 1 29 500 トヨタ

2 2 36 600 日産

3 3 32 400 トヨタ

4 4 45 1000 ホンダ

5 5 44 1200 ジャガー

6 6 47 800 ホンダ

7 7 28 400 なし

8 8 37 300 なし

9 9 37 600 ダイハツ

10 10 45 900 日産

car には age という量的変数とmakerというカテゴリカル変数の両方が含まれている。Rで

は、上のcarのようなデータセットをデータフレームと呼んでいる。

データフレーム内の変数を指定するときは、car$age, car$maker のように、

データフレーム名$変数名

と書く。データフレーム名$変数名を入力すれば、その変数の中身が表示される。

> car$maker

[1] トヨタ 日産 トヨタ ホンダ ジャガー ホンダ なし なし

[9] ダイハツ 日産

Levels: ジャガー ダイハツ トヨタ なし ホンダ 日産

> car$age

[1] 29 30 21 45 56 51 28 22 25 45

car$ageの平均を求めるには、以下のように入力する。

> mean(car$age)

[1] 38

課題 2.4 Rを用いてcar$incomeの平均を求めなさい。

課題 2.5 （復習）以下の計算を、Rのスクリプトchap2.rに書き、Rのコードを実行し

て求めなさい。結果が正しいかどうか確認すること。

3_×5 34

56_÷7

√

2.8

関数の作成

これまで、見てきたようにRでは、関数の利用が重要である。実は、新しい関数を自分で作

成することができる。

関数定義の基本形式は以下である。

関数名 <- function(引数){

処理の一行目

処理の二行目

・・・・・

return(計算結果を返す)

}

この章で用いた関数

関数名 関数の意味

q Rを終了する

sqrt 平方根を求める

c 複数の数値からベクトルを作成する

mean データの平均を求める

第

3

章

1

変数のデータ分析

3.1

ヒストグラム

1変数データの分析で最も基本的なのが，データの分布の状況を見ることである．ヒストグ

ラムとは，データがある範囲を複数の区間（階級という）に分割し，各階級に属する観測値の

個数（度数という）を柱の面積で表したものである．図 3.1は，2010年7月−2011年2月の日

経平均株価の収益率1のヒストグラムである．ヒストグラムとは，各区間に含まれるデータの個

数（度数という）を柱の面積に比例するように描かれた図のことである．日経平均収益率の分

布がおおよそ左右対称であることがわかる．

Histogram of Nikkei$return

Nikkei$return

Frequency

−0.04 −0.03 −0.02 −0.01 0.00 0.01 0.02 0.03

0

5

10

15

20

25

30

図 3.1: 日経平均の収益率のヒストグラム

経済・経営データでは，分布が左右対称でないものも多い．図 3.1は，平成22年の23区住

居地域の公示価格データ（単位は円/m2）のヒストグラムである．20−80万円の範囲に多く

のデータがあるが，100万円を超えるデータもそれなりにあり，右側に伸びている．このよう

な分布は右に歪んだ分布という．経済データでは分布が右に歪んでいることがよくある（所得，

貯蓄など）．逆に左側にすそが伸びている場合，左に歪んだ分布という．

3.2

データの代表値−平均、メディアン−

ヒストグラムよりもさらにデータを要約して，データの特徴を表すいくつかの数値を求める

ことも多い．まず，1つの数値でデータを代表させようとする場合，これをデータの代表値と

Histogram of kouji_23ku$kakaku

kouji_23ku$kakaku

Frequency

0 500000 1000000 2000000 3000000

0

50

100

150

200

250

図 3.2: 23区住居地域の公示価格データのヒストグラム

いう．まず，データを

x1, x2. . . . , xn

で表そう．

3.2.1

平均

データx₁, x₂. . . . , x_nの平均は，

¯

x= 1

n

n X

i=1

xi

で定義される．平均は平均値とも呼ばれる．

平均の性質1

n X

i=1

(xi−x¯) = 0

平均の性質2

平均は外れ値に引っ張られやすい．

ただし，外れ値とは他のデータの値からかけ離れて大きい値，あるいは小さい値のことをいう．

例 3.1 以下の数値はある会社の従業員の年間所得（万円）が以下であるとする．

330,280,230,240,390,290,340,1580

このときx¯= 460（万円）と計算できる．この会社の従業員8人のうち7人は所得が平均460万

3.2.2

メディアン

メディアンは，データを小さいものから大きいものへ順番に並べた時に中央に位置する値で

ある．大きさの順に並べたデータを

x(1) < x(2) <· · ·< x(n)

で表す．メディアンは

M d=

(

x₍n+1

2 ) n

が奇数の場合,

x(n/2)+x(n/2+1) nが偶数の場合.

で定義される．メディアンは外れ値に影響を受けにくい特徴がある．上の例で、年間所得のメ

ディアンは？

モード，メディアン，平均の関係

右にゆがんだ分布ではモード<メディアン<平均，になる．左にゆがんだ分布では，平均<

メディアン<モード，の順番になる．

3.3

データの散らばりの尺度−分散，標準偏差−

平均はデータの一つの代表値であった。代表値とともに、データがどれだけ散らばっている

かを測る尺度も有用である。散らばりの尺度として用いられるのが分散(variance)で、

S2 ₌ 1

n

n X

i=1

(xi−x¯)2 (3.1)

で定義される。分散の正の平方根

S =√S2

は標準偏差(standard deviation)という．分散は

s2 = 1

n₋1

n X

i=1

(xi−x¯)2 (3.2)

で定義されることもある．n−1 で割る理由は後ほど説明される．s2を不偏分散ともいう．標

準偏差は

s=√s2

3.4

R

で

1

変数データ分析

まず、apartment.csvをRに読み込み、apartmentという名前のデータフレームを作り

ます。

> apartment<-read.csv("D:/rdata/apartment.csv")

上のコードでは，Dドライブにデータがあることを想定している．

> head(apartment) %関数 head はデータフレームの最初の7行を表示する関数。

area time years yachin shurui1 shurui2

1 7.4 14 39 1.7 1K

2 8.3 9 29 3.8 1K

3 10.0 13 31 3.8 1K

4 10.1 18 15 4.7 1K

5 11.8 15 19 4.0 1K

6 12.5 15 11 5.7 1K

このデータは203件分の賃貸アパートマンションのデータであり、areaは物件の面積(m2)、time

は駅までの時間距離(分)、years は建物の築年数(年)、yachin は家賃(万円)を意味する．

ここでは，area(面積)の分析を行う．まず，データがどのようにばらついているかを視覚的

に見るため、関数 hist()でデータのヒストグラムを作成する．

> hist(apartment$area)

Histogram of apartment$area

apartment$area

Frequency

0 20 40 60 80 100 120

0

10

20

30

40

50

60

図 3.3: 専有面積のヒストグラム

関数hist()で,hist(変数名, breaks=N)とすると階級の数がNに近くなるように階級

幅を選んでヒストグラムが描かれる．

データの平均，メディアン，分散，標準偏差はRの関数mean()，median()，var()，sd()

> mean(apartment$area)

[1] 29.43793

> median(apartment$area)

[1] 24

> var(apartment$area)

[1] 229.4022

> sd(apartment$area)

[1] 15.14603

ボックスプロットという図が用いられることもある．

> boxplot(apartment$area)

20

40

60

80

100

120

図 3.4: 専有面積のボックスプロット

関数名 関数の意味

hist データのヒストグラムを描く

mean データの平均を求める

median データのメディアンを求める

var データの分散を求める

sd データの標準偏差を求める

第

4

章

2

変数のデータ分析

4.1

散布図

> apartment<-read.csv("Z:/rdata/apartment.csv")

> plot(yachin˜area,data=apartment)

˜ を使うことに注意。

20 40 60 80 100 120

5

10

15

20

25

area

yachin

図 4.1: 家賃と専有面積の散布図

図4.1は家賃と専有面積の散布図である。家賃と専有面積はおおよそ右上がりの直線に近い

関係がある。このような関係を正の相関関係という。

家賃と築年数の間にはわずかであるが、右下がりの直線に近い関係がある。このような関係

を負の相関関係という。

課題

yachin と area の散布図を描くRのコードをスクリプトに書き、実行しなさい。

yachin と years の散布図を描くRのコードをスクリプトに書き、実行しなさい。

スクリプトファイルにchap4.r という名前をつけてフォルダrdata に保存しなさい。

0 10 20 30 40

5

10

15

20

25

years

yachin

図 4.2: 家賃と築年数の散布図

4.2

相関係数

相関関係を数値で測定する尺度として相関係数がある。データが(x₁, y₁),(x₂, y₂), . . . ,(x_n, y_n)

で与えられている時、x と y の相関係数は

r=

Pn

i=1(xi−x¯)(yi−y¯) pPn

i=1(xi−x¯)2 Pn

i=1(yi−y¯)2

で定義される。どのようなデータについても−1≤r ≤1である。Rでは関数corで相関係数

を求められる。

> cor(apartment$yachin,apartment$area)

[1] 0.8592066

>

第

5

章

回帰分析

5.1

単回帰分析

相関係数を用いて2変数のデータを分析するときは、二つの変数を対称的に扱った。ある変

数の動きを他の変数の一次式で説明する分析手法を回帰分析という。前章の図4.1を見ると、賃

貸アパート・マンションの専有面積と家賃の間には強い直線的な関係があることがわかる。こ

のような場合、家賃と専有面積の関係を表す近似的な一次式

yachin =a+b_×area

を求めると便利である。ここでは、area が yachin を説明する式を考えている。値を説明され

る方の変数を被説明変数、説明する方の変数を説明変数という。通常は、最小2乗法と呼ばれ

る方法によって、係数a と bを求める。a, bを求めることが回帰分析の第一歩である。Rでは

関数lmによって回帰分析を行う。関数lmは

lm(被説明変数名˜説明変数名,data=データフレーム名)

という形をとる。

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

> result1<-lm(yachin˜area,data=apartment)

> result1

Call:

lm(formula = yachin ˜ area, data = apartment)

Coefficients:

(Intercept) area

3.1553 0.2118

関数lmの結果をresult1に保存した。result1と入力して求めた係数がわかった。結果は

\

yachin = 3.1553 + 0.2118_×area (5.1)

と書いて表す。式(5.1)を回帰式という。次に、求めた回帰式を散布図上に描いてみよう。

> plot(yachin˜area,data=apartment)

> abline(result1)

20 40 60 80 100 120

5

10

15

20

25

area

yachin

図 5.1: 回帰直線

回帰式のあてはめ

専有面積と家賃の間には原因‐ 結果関係があり、さらに図??から、専有面積と家賃の間には

近似的に直線的な関係があることがわかる。このようなデータに対しては回帰分析が有効であ

る。家賃をyで、専有面積をx で表すとする。このときyとxの近似的な関係を表す一次関数

y=a+bx (5.2)

を求めることができれば便利である。このような一次関数を求めるために，実用上もっともよ

く使われるのが最小二乗法である。この方法によればa,b は次の公式で与えられる。

b=

n X

i=1

(xi −x¯)(yi−y¯) n

X

i=1

(xi−x¯)2

a= ¯y₋bx¯

結果は以下のようになる。

¯

y= 8, x¯= 20

b= 146

460 = 0.3174, a= 8−0.3174×20 = 1.652

得られた a, b を使った式

y=a+bx

付録：最小

2

乗法の考え方

( x

i

, y

i

)

a + b x

x

i

y

i

x y

a+bx

i

図 5.2: 最小二乗法の考え方

x₋y 平面上の直線は y =a+bx の形に書ける。ただし a, b は定数である。観測点 (x_i, y_i)

から垂直に下ろした線が直線y=a+bx とぶつかる点は(x_i, a+bx_i) である。したがってその

2点の差はy_i−a−bx_i となる。そこで、直線y=a+bx がデータ(x₁, y₁), . . . ,(x_n, y_n)からど

の程度「ずれ」ているかを測る尺度として、「ずれ」の2乗の和

S(a, b) =

n X

i=1

(yi−a−bxi)2

を考える。S(a, b) を最小にするa, b の値を求める方法を最小二乗法という。

決定係数

最小2乗法を使って推定された回帰モデルを

y=a+bx (5.3)

と書こう。この回帰モデルを用いて、 x_i の値から計算される y の値をyˆ_i と表す。すなわち

ˆ

yi =a+bxi, i= 1, . . . , n.

である。yˆ_i を y_i に対する回帰値という。当然、実際の観測値y_i と回帰値yˆ_i は一致するとは限

らない。y_i と yˆ_i の差

ei =yi−yˆi, i= 1, . . . , n (5.4)

を残差という。(5.4)式を変形すると

y i y

x y = α + β x

x i y

i

図 5.3: 回帰直線

と書ける。この式から、yˆ_i は y_i の値のうち回帰モデルで説明できる部分、e_i は回帰モデルで

説明できない部分である。データ全体で見たy_i の変動のうち yˆ_i の変動で説明できる部分が大

きければよい。まず、(5.5)式の両辺から y¯を引く。

yi−y¯= (ˆyi−y¯) +ei

次に、この式の両辺を2乗すると

(yi−y¯)2 = (ˆyi−y¯)2+ 2 (ˆyi−y¯)ei+e2i

を得る。これを、すべての y_i について合計すると以下のようになる。

n X

i=1

(yi−y¯)2

| {z }

yiの変動 =

n X

i=1

(ˆyi−y¯)2

| {z }

ˆ yiの変動

+

n X

i=1

e2_i

| {z }

残差の変動

(5.6)

(5.6)式は、y₁, y₂, . . . , y_n 変動が、回帰モデルで説明できる部分と説明できない部分とに分解で

きることを意味している。(5.6)式を、簡単に次のように書きなおす。

ST |{z}

全変動

= SR |{z}

回帰平方和

+ Se |{z}

残差平方和

(5.7)

ST に対して S_R の占める割合が大きければ大きいほど、推定された回帰モデルのデータに対す

る当てはまりがいいことを意味する。この割合をR2 で表す。

R2 = SR

ST

= ST −Se

ST

= 1₋ Se

ST

(5.8)

このR2 を推定された回帰モデル(5.3)の決定係数という。S_T, S_R,S_e はすべて正だから、

0≦R2 ≦1

5.2

回帰係数に関する推測

以上の節では、回帰直線を求めて、その当てはまりの良さ（決定係数）について説明した。

得られたデータは、いろいろな誤差や不確実性を含むので、得られた回帰式の係数は絶対的な

ものではない。より進んだ分析手法では、データの不確実性を考慮した上で、得られた回帰式

から何が結論できるかを調べる。

yi, xi それぞれは被説明変数、説明変数のi番目の観測値とする。y_i は x_i と誤差項 u_i とで

以下のように決定されるものと仮定する。

yi =α+βxi+ui, i= 1,2, . . . , n (5.9)

ただしu_i, i = 1,2, . . . , nは独立で、それぞれ正規分布 N(0, σ2) に従う。またα, βは未知の

定数である。このモデルを単回帰モデルという。β の最小2乗推定量をβˆで表す。上の仮定の

下で、

ˆ

β₋β σ_βˆ ∼

N(0, 1) (5.10)

であることを証明することができる。ただしσ_ˆ

β はβˆの標準偏差で，

σ_βˆ =

σ

s _n X

i=1

(xi−x¯)2

(5.11)

で与えられる．(5.10)式左辺のσ_ˆ

βの値は未知なので，その推定量で置き換える必要がある．誤

差項の分散σ2の推定量は残差e₁, . . . , e_n を用いて

ˆ

σ2 = 1

n₋2

n X

i=1

e2_i (5.12)

で与えられる．ˆσ2を用いれば，βˆの標準偏差の推定量は以下のようになる．

se( ˆβ) = s _n σˆ X

i=1

(xi−x¯)2

(5.13)

この推定量はβˆの標準誤差(standard error)とよばる．このとき

ˆ

β₋β

se( ˆβ) ∼t(n−2) (5.14)

が成り立つ．ただし，t(k) は自由度kのt分布を表す．この事実は，βに関する仮説検定に利

5.2.1

_β

についての仮説検定

回帰分析では、以下のような帰無仮説の仮説検定が重要である。

H0 :β = 0 (5.15)

対立仮説としては両側対立仮説 H₁ : β 6= 0 を設定する場合は両側検定と言われ、片側対立仮

説H₁ : β > 0 あるいは H₁ : β < 0を設定する場合は片側検定という。この仮説検定が重要な

のは、帰無仮説(5.15)が棄却されれば、β 6= 0 つまり、得られた回帰式は意味のあるものと結

論できるからである。この仮説検定の検定統計量は

t= βˆ

se( ˆβ)

である．両側検定の場合，この値の絶対値が十分大きければ、帰無仮説(5.15)が間違っている

証拠と考えられる．自由度n−2 の t分布の上側 100α%点をt_α(n−2)で表すものとする．t∗

をt検定統計量の実現値としよう．t∗はt値とよばれる．もし

|t∗_{|> t}

α/2(n−2) (5.16)

であれば，有意水準αでH₀を棄却する．(5.16)式が成り立たなければ，H₀は有意水準αでは

棄却されない．

H0を棄却できるかどうかは，P 値を用いて判断することもできる．P 値は両側検定の場合，

t検定統計量の絶対値が帰無仮説のもとで |t∗|以上になる確率のことである．P 値がα以下で

回帰分析の詳細な結果はsumary関数で見ることができる。

> summary(result1)

Call:

lm(formula = yachin ˜ area, data = apartment)

Residuals:

Min 1Q Median 3Q Max

-8.3304 -0.7731 -0.1072 0.5670 13.9522

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 3.155272 0.294394 10.72 _{<2e-16 ***}

area 0.211831 0.008897 23.81 _{<2e-16 ***}

---Signif. codes: 0‘ ***’ 0.001‘ **’ 0.01‘ *’ 0.05‘ .’ 0.1‘ ’ 1

Residual standard error: 1.915 on 201 degrees of freedom

Multiple R-squared: 0.7382, Adjusted R-squared: 0.7369

F-statistic: 566.9 on 1 and 201 DF, p-value: < 2.2e-16

上の関数lmの結果の読み方を説明する．Coefficients:の下の表記の説明

表 5.1:

Estimate 回帰パラメータの推定値

Std. Error 標準誤差

t value t値

Pr(>|t|) P値

Multiple R-squared 決定係数

Adjusted R-squared 自由度調整済み決定係数

上の結果をまとめると、以下のようになる。

回帰式

\

yachin=3.16 + 0.21area

(10.72) (23.81)

括弧内はt 値。決定係数R2 = 0.738専有面積は有意である。

5.3

役に立つテクニック

cv2<-read.csv("cv2.csv")

散布図の点にラベルを付ける。

plot(conv˜pop,data=cv2)

plot(conv˜pop,data=cv2)

text(cv2$pop,cv2$conv,labels=cv2$todofuken,pos=3)

2000 4000 6000 8000 10000 12000

0 1000 2000 3000 4000 5000 6000 pop conv ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

図 5.4: 都道府県別コンビニ数と人口

5.4

重回帰分析

前節までは、説明変数が一つの場合の回帰分析について述べてきた。説明変数が二つ以上あ

る場合の分析を重回帰分析という。

yi =α+β1x1i+β2x2i +βkxki+ui (5.17)

apartmentのデータで、家賃を被説明変数、面積、築年数、時間距離の三つを説明変数と

して重回帰分析を行ってみよう。

> apartment<-read.csv("c:/rdata/apartment.csv")

> regresult<-lm(yachin˜area+years+time,data=apartment)

> summary(regresult)

Call:

lm(formula = yachin ˜ area + years + time, data = apartment)

Residuals:

-5.75155 -0.78941 -0.05127 0.64746 13.05786

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 6.208969 0.440257 14.103 _{< 2e-16 ***}

area 0.208335 0.007386 28.208 _{< 2e-16 ***}

years -0.114126 0.013563 _{-8.414 7.64e-15 ***}

time -0.163914 0.029711 _{-5.517 1.06e-07 ***}

---Signif. codes: 0‘ ***’ 0.001‘ **’ 0.01‘ *’ 0.05‘ .’ 0.1‘ ’ 1

Residual standard error: 1.578 on 199 degrees of freedom

Multiple R-squared: 0.824, Adjusted R-squared: 0.8213

F-statistic: 310.6 on 3 and 199 DF, p-value: < 2.2e-16

>

課題(提出)

被説明変数をyachin, 説明変数をarea, yearsとし回帰分析を行いなさい。回帰式、t値、決定

係数R2を書き、解釈を書きなさい。

被説明変数をyachin/area,説明変数をarea, yearsとし回帰分析を行いなさい。回帰式、t値、

5.5

ダミー変数

回帰分析で名義変数を説明変数として使うこともできる。賃金構造基本調査から得られる年

齢別、学歴別の平均賃金データを見てみよう。

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

> plot(wage˜age,data=wagestat)

> plot(wage˜age,pch=as.integer(education),data=wagestat)

30 40 50 60

200

250

300

350

400

450

500

age

wage

図 5.5: 賃金カーブ

散布図から、賃金は年齢の関数であり、大学・大学院卒の賃金カーブが高卒の賃金カーブの

上に来ている。

> wagereg<-lm(wage˜age+ageˆ2+education,data=wagestat)

> summary(wagereg)

Call:

lm(formula = wage ˜ age + ageˆ2 + education, data = wagestat)

Residuals:

Min 1Q Median 3Q Max

-140.28 -39.31 15.45 48.92 90.57

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 117.546 53.854 _{2.183 0.043377 *}

education大学・大学院卒 128.070 31.623 4.050 0.000832 ***

---Signif. codes: 0‘ ***’ 0.001‘ **’ 0.01‘ *’ 0.05‘ .’ 0.1‘ ’ 1

Residual standard error: 70.71 on 17 degrees of freedom

Multiple R-squared: 0.6119, Adjusted R-squared: 0.5663

F-statistic: 13.4 on 2 and 17 DF, p-value: 0.0003204

レポート

提出期限： 2012年7月9日 (月)22時までにメールで提出する．7月10日授業時に紙に印刷し

たものを提出して下さい．

提出内容： 各自がデータを用意し，以下のように重回帰分析を行ってください

１．自分が興味を持つデータを用意する．ただし以下の条件を満たしていることが重要です．

（１）被説明変数、説明変数と考えられる変数があること．説明変数は２個以上．説明変数と

してダミー変数が含まれていてもいい．

（２）標本の大きさ（観測値の個数）はおおよそ２０以上であることが望ましい．標本の大き

さが小さいと推定の精度が低くなるため．

レポートの内容（量はA4で5枚以上）

1．分析の目的（分析後に考えても良い）

2．データの出所．変数の定義．

3．各変数を分析する．（ヒストグラム、代表値、散布図、相関係数など）

4．重回帰分析

4−1．二つ以上の回帰式を推定する．有意ではない説明変数が含まれていても構わない．

（異なる説明変数の組合せの回帰式．あるいは説明変数を変換してみる．）回帰分析の結果を書

く(Rの出力を貼り付けるのではなく、標準的な方法で結果を書く)．

4−2．結果の解釈．結論．

第

6

章

カテゴリカルデータの分析

いままでの章では、長さ、重さ、金額などのように数値で測ることができる変数とデータを

扱ってきた。アンケート調査などでは、性別、職業など複数のカテゴリーに分類される変数も

重要である。このような変数をカテゴリカル変数という。特に、名義尺度で測定された変数を

名義変数、順序尺度で測定された変数を順序変数という。

6.1

データの読み込み、度数分布表

keitai.csvは200人の社会人（男女は同数）について、契約している携帯電話の会社を調べた

結果（架空データ）である。まず、このデータをkeitai という名前のデータフレームとして読

み込む。

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

> head(keitai)

id sex carrier

1 1 女 docomo

2 2 女 au

3 3 女 docomo

4 4 男 au

5 5 男 docomo

6 6 男 docomo

カテゴリカルデータの度数分布表は関数 xtabsで書くことができる。ここでは、変数carrier

を集計し、オブジェクトx に保存する。

> x<-xtabs(˜carrier,data=keitai)

> x

carrier

au docomo softbank

61 103 36

次に集計結果を度数の大きさの順に並べ替えてみる。

> x<-x[order(x,decreasing=TRUE)]

> x

carrier

docomo au softbank

docomo

au softbank

図 6.1: 円グラフ

円グラフの描き方

> pie(x,radius=1,clockwise=TRUE)

6.2

2

元分割表

二つ以上のカテゴリカル変数の関係を分析するときには分割表（クロス表ともいう）を用い

る。関数 xtabsによって分割表を作成できる。

xtabs(˜変数1+変数2, data=データフレーム名)

というコマンドで変数1と変数2の分割表が作成される。データフレームkeitai を使って性別

と携帯会社の分割表を作成してみよう。

> y<-xtabs(˜sex+carrier,data=keitai)

> y

carrier

sex au docomo softbank

女 22 58 20

男 39 45 16

このような二つの変数の分割表を、2元分割表、あるいは簡単に2元表という。上のコマンド

では、2元表に集計した結果をオブジェクトyに保存した。関数prop.tableを用いれば、割

合で表示した分割表が作成できる。

> prop.table(y,margin=1)

carrier

sex au docomo softbank

女 0.22 0.58 0.20

男 0.39 0.45 0.16

6.3

独立性の検定

この節では、2元表に対して最もよく行われる仮説検定である独立性の検定を説明する。仮

説検定とは、データを発生した確率分布についての仮説が正しいのかどうかをデータを用いて

検証する方法であり、統計学のあらゆる応用で利用されている。まず、基準となる仮説を設定

する。この仮説を帰無仮説という。次に、帰無仮説が成り立たないときに成り立つ仮説を設定

する。この仮説を対立仮説という。

もう一度、性別と携帯会社の分割表を見てみよう。

au docomo softbank 合計

女 22 58 20 100

男 39 45 16 100

合計 61 103 36 200

「性別と携帯会社選択は独立である」という仮説を帰無仮説とする。このような仮説は独立

性の仮説という。独立性の仮説の下で、それぞれのセルの期待度数を求めると、以下のように

なる。

au docomo softbank 合計

女 30.5 51.5 18 100

男 30.5 51.5 18 100

合計 61 103 36 200

実際の度数と独立性の仮定の下での期待度数との乖離は

χ2 = (22−30.5)

2

30.5 +

(58₋51.5)2

51.5 +· · ·+

(16₋18)2

18 = 6.823

で測れる。この値が、帰無仮説の下では起こり得ないような大きい値であれば、帰無仮説が成り

立っていない証拠と考えられる。帰無仮説の下で起こる確率が非常に小さいような値が得られ

た時には帰無仮説を退ける。これが仮説検定の考え方である。小さい確率として通常用いられ

ているのが1% と 5% で、この確率を有意確率という。独立性の検定の結果にはP値(p-value)

と呼ばれる値が出力される。P値が0.05よりも小さければ、帰無仮説（独立性）は棄却される。

標本の大きさが大きいときには、帰無仮説の下でのχ2検定統計量の確率分布は、自由度2の

χ2

分布で近似される1。

Rでは、関数summaryを使えば、二つのカテゴリカル変数の独立性の検定が行うことがで

きる。

> summary(y)

Call: xtabs(formula = ˜sex + carrier, data = keitai)

Number of cases in table: 200

Number of factors: 2

Test for independence of all factors:

Chisq = 6.823, df = 2, p-value = 0.03299

>

summary(y)の出力の下の2行は、独立性のχ2検定の結果である。p-value は日本語でP値

と呼ばれる。P値とは、帰無仮説が正しいときに、得られた結果よりも極端な値が起こる確率

のことである。P値<0.05であれば有意水準0.05で帰無仮説は棄却される。携帯電話会社デー

タでは、性別と携帯会社選択が独立であるという帰無仮説は有意水準5%で棄却される。つま

り、性別によって携帯会社選択の確率が異なる、と結論することができる。

上のデータ例では、性別と携帯会社選択が独立であるという帰無仮説は有意水準5% で棄却

された。では、P値が0.05よりも大きいときはどうだろうか。P値が0.05 よりも大きいとき

は、帰無仮説は有意水準5%では棄却されない。したがって、二つの変数が独立ではない、と結

論できない。しかし、独立である、とも結論できない。

6.4

3

元分割表

ファイルalcohol.csvは、大学生100人について、性別、サークル活動、飲酒習慣につい

て調査したデータ（架空のデータ）である。

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

> head(alcohol)

id sex circle alcohol

1 1 男 あり あり

2 2 女 なし あり

3 3 男 あり あり

4 4 男 あり あり

5 5 男 あり あり

6 6 女 なし なし

このデータは三つの変数からなり、sexは性別を、circleは大学におけるサークル活動の有無を、

alcoholは飲酒習慣の有無を意味している。まず、circleと alcohol の2元表を作成してみよう。

> y<-xtabs(˜circle+alcohol,data=alcohol)

> y

alcohol

circle あり なし

あり 70 50

なし 30 50

> prop.table(y,margin=1)

alcohol

circle あり なし

あり 0.5833333 0.4166667

なし 0.3750000 0.6250000

circleと alcoholの2元表から、サークル活動をしている学生としていない学生では、飲酒習慣

の割合に差があることがわかる。

> summary(y)

Call: xtabs(formula = ˜circle + alcohol, data = alcohol)

Number of cases in table: 200

Number of factors: 2

Test for independence of all factors:

summary(y)の結果を見ると、χ2検定のP値は0.004 で0.01よりも小さく、サークル活動の 有無と飲酒習慣の有無との間には関連があると結論できる。三つのカテゴリカル変数のすべて

を使って分割表を作成するときは、以下のようにxtabsを用いる。

xtabs(˜変数1+変数2+変数3, data=データフレーム名)

このコマンドで、変数3で層別された2元表が作成される。このような分割表を3元表という。

以下で、性別で層別した2元表（3元表）を作成してみる。

> xtabs(˜circle+alcohol+sex,data=alcohol)

, , sex = 女

alcohol

circle あり なし

あり 10 30

なし 15 45

, , sex = 男

alcohol

circle あり なし

あり 60 20

なし 15 5

性別を考慮した3元表で見てみると、実はサークル活動の有無によっては飲酒習慣の割合は変

わらないことがわかる。この例のように、層別した分割表における変数の連関と層別しない分

割表での連関が異なってしまう現象をシンプソンのパラドックスという。

このデータでは、性別がサークル活動と飲酒習慣という両方の変数に影響を及ぼしているこ

とがシンプソンのパラドックスが生じている原因なのである。二つの変数の間に連関が見られ

る時も、二つの変数に影響を与える変数がないかどうかよく考える必要がある。

課題 6.1 性別と飲酒習慣の分割表を作り、ワードのファイルの書きなさい。独立性の検定の結

果も書く。

役に立つ方法

層別して分割表を求め、独立性の検定を行う。

サークル活動と飲酒習慣の分割表を男のデータだけで作り、独立性の検定を行うのは、以下

のようにすればよい。

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

> table<-xtabs(˜circle+alcohol+sex,data=alcohol)

>

> table[,,"男"]

alcohol

あり 60 20

なし 15 5

> chisq.test( table[,,"男"])

Pearson’s Chi-squared test with Yates’ continuity correction

data: table[, , "男"]

X-squared = 0.0833, df = 1, p-value = 0.7728

chisq.test(分割表)で、独立性の検定を行う。上の場合、帰無仮説は棄却されない。

課題 6.2 alcoholデータで、サークル活動と飲酒習慣の分割表を男のデータだけで独立性の

検定を行いなさい。

モザイクプロット

分割表を視覚的に表すグラフとしてモザイクプロットがある。

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

> table1<-xtabs(˜circle+alcohol,data=alcohol)

> mosaicplot(table1,col=TRUE)

table1

circle

alcohol

... ...

...

...

図 6.2: サークル活動と飲酒頻度のモザイクプロット

6.5

順序カテゴリカルデータの分析

順序尺度で測定された変数を順序変数という。たとえば、顧客満足度調査で得られる回答は

順序変数である。

２つの順序変数の分析も分割表で行う。下の表は、ある会社の従業員の所得の水準と仕事に

不満 どちら 満足

でもない

低位 8 3 1

中位 5 9 4

高位 2 5 10

この分割表のデータでは、所得水準が低位で仕事に不満な人、所得が高水準で仕事に満足な人

の割合が多いので、所得水準と仕事の満足度の間には、正の連関があるであろうことがわかる。

6.5.1

順序変数間の連関の尺度

順序変数のとる値には、大小関係しかないので、相関係数を計算することができない。二つ

の順序変数の連関の程度の尺度としてGoodmanとKruskalによるγ(ガンマ) の考え方を説明

する。

従業員満足度調査の分割表を見てみよう。このデータに含まれる従業員のペアを考え、一方

の従業員の所得の水準が他方よりも高く、かつ満足度も他方よりも高い水準であるとき、この

ペアをconcordantペアという。また、一方の従業員の所得の水準が他方よりも高いが、満足度

は他方よりも低い水準であるとき、このペアをdiscordantペアという。concordant pairs の個

数Cは

C = 8(9 + 4 + 5 + 10) + 3(4 + 10) + 5(5 + 10) + 9(10) = 224 + 42 + 75 + 90 = 431

と求められる。discordant pairsの個数Dは

D= 3(5 + 2) + 1(5 + 9 + 2 + 5) + 9(2) + 4(2 + 5) = 21 + 21 + 18 + 28 = 88

と求められる。C/(C+D) はconcordant pairs の割合であり、D/(C+D)はdiscordant pairs

の割り合いをを意味する。γ は

γ = C−D

C+D

で定義される。γ はconcordant pairsの割合からdiscordant pairsの割り合いを引いたものだか

ら、−1≤γ ≤1であり、この値が1に近ければ正の連関があり、−1に近ければ負の連関があ

ると言える。γ は

γ = C−D

C+D =

431₋88 431 + 88 =

343

以下の課題で、Rを使って分割表からγを求める。（関数定義については、第2.8節を参照

せよ。）

課題 6.3 （ 重 要 ）以 下 のRコ ー ド を ス ク リ プ ト に コ ピ ー し な さ い ．以 下 で 定 義 さ れ た 関 数

gamma2は、与えられた2元分割表からγの値を求める関数である。（青木繁伸教授のHPから 転載した）

gamma2 <- function(f)

{

R <- nrow(f)+2

C <- ncol(f)+2

g <- matrix(0, nr=R, nc=C)

g[2:(R-1), 2:(C-1)] <- f

cc <- dd <- 0

for (i in 2:(R-1)) {

for (j in 2:(C-1)) {

cc <- cc+g[i,j]*sum(g[1:(i-1), 1:(j-1)], g[(i+1):R, (j+1):C]) dd <- dd+g[i,j]*sum(g[1:(i-1), (j+1):C], g[(i+1):R, 1:(j-1)]) }

}

return((cc-dd)/(cc+dd))

}

CS.csvの中の2つの変数を選び、分割表を作成しなさい．さらにgamma を求めなさい．結