第
1
章
統計学とは
自然現象、社会現象は多数の要因が複雑にからみあい、さまざまな出来事が生じています。
私たちは、これらの現象を客観的に理解するために、いろいろなデータを収集し分析していま
す。この授業では、統計学の基礎を理解しRを用いて実際にデータを分析することによって、
統計学の基礎を身につけます。各自の研究でデータ解析ができるようになることを目標とする。
現代の統計学は確率論という数学が基礎にあり、統計学の理論を深く理解するためには、一
定程度の数学の知識が必要で、多くの時間勉強する必要があります。しかしながら、この授業
では統計学の理論の説明は最小限にとどめ、実際に統計学を使ってみることを目的とします。
本講義についての情報は,私の作成するHP https://sites.google.com/site/junfukuchistat/
に書きます。ここから, 授業(2011年度) → 統計学 と進むとこの授業に関する情報が見られ
ます。講義ノートの最新版は、授業期間中このページからダウンロードできます。
課題が講義ノートの中にいくつかありますが、これらの課題は必ずやってみてください。
課題(提出)は所定の期日まで提出してください。毎回、提出用の課題が出題されます。課
題提出先は
(福地のメールアドレス)です。
授業方法、成績の評価の方法:教科書と講義ノートを用いて進める。課題、小テスト、レポー
ト。毎回課題があるので、出席が重要です。課題はRを使用するもの。小テストは2回。
1.1
データの種類
測定の尺度からデータの種類を分類してみる。実際に多く使われるデータは、比例尺度、順
序尺度、名義尺度である。比例尺度とは、長さ、重さ、広さ、金額などのような量を測る尺度で
あり、比例尺度で表される変数は連続的な値をとるから連続変数とよばれる。順序尺度は成績
の「不可、可、良、優」のように順序に意味のあるような尺度である。名義尺度は、性別「男、
女」や、産業の分類「第一次、第二次、第三次産業」などのように複数カテゴリーに分割され
るが、本来数値的意味はない尺度である。順序尺度や名義尺度で表される変数をカテゴリカル
変数といい、そのデータをカテゴリカルデータという。経済、経営では、連続データもカテゴ
リカルデータのどちらも重要である。
本講義ノート作成のために参考にした文献
「The R Tips」、舟尾暢男著、2009年、オーム社
「Rによるやさしい統計学」、山田剛史・杉澤武俊・村井潤一郎著、2008年、オーム社
10月13日(木)準備
1.ブラウザーを起動し、福地のHPを開いてください。福地のHPは
https://sites.google.com/site/junfukuchistat/
です。このページを大学からも家からも見られるようにしてください。
2.デスクトップにrdataという名前のフォルダを作る。
2.福地のHPの「2011年授業」→「統計学」のページに行き、統計学の講義ノートをダウン
ロードし、フォルダrdataに保存します。
3.次に、二つのcsvファイル(apartment.csv, car.csv)をダウンロードし、フォルダ rdataに
保存します。
4.Rを起動します(スタートボタン→すべてののプログラム→R)。
5.Rのメニューから,ファイル→フォルダの変更を選び,desktop のフォルダrdataを指定
します。
第
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
スクリプトファイルによる実行
準備
スクリプトファイルによるコマンドの実行
今までは、コマンドの実行を、直接プロンプトの右側に打ち込んで行ってきた。この方法だ
と、毎回コマンドを打ち込む必要があるし、入力したコマンドに誤りがあったときには修正が
しにくい。
スクリプトファイルによるコマンドの実行の手順
1.Rのメニューの「ファイル」の中の「新しいスクリプト」をクリック
2.「Rエディタ」というウインドウにRのコマンドを入力する。
3.入力したRのコマンドを反転させる(左クリックしてドラッグする)。
4.右クリックして、「カーソル行あるいは選択中のRコードを実行」をクリックするとRの
コマンドが実行される。
5.一連の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などのように章ごとに名
前を付けて、デスクトップに保存してください。
2.6
R
によるデータの扱い
実際のデータを扱う方法を説明する。最も単純な方法は、直接データを入力する方法である。
関数cは、複数の数値を合わせてベクトルを作成する。たとえば、高校生3人のお小遣いのデー
タ 6000,8000,12000をベクトル kozukai として定義するには以下のようなコマンドを入力す
ればよい。kozukaiを定義した後、kozukaiと入力し、kozukaiの中身を表示させる。
> kozukai<-c(6000,8000,12000) > kozukai
[1] 6000 8000 12000
Rでは、定義したベクトルに対して関数を適用し、統計分析をすることができる。お小遣いデー
タの平均を求めてみよう。データの平均を求めるRの関数はmeanである。
> mean(kozukai) [1] 8666.667
3人のお小遣いの平均が8666.7円であることがわかった。このように関数を多用することがR
の1つの特徴である。
ベクトルについての計算も簡単に行えます。ベクトルに対して、+, -, *, /を用いると要素ご
とにその計算を行います。
> kozukai
> 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つが個体、縦のなら
びが変数になるように入力しておかなければならない。
課題 2.4 ファイルcar.csv を開き、中身を見なさい。変数はいくつありますか?観測値 はいくつありますか?
R
で
CSV
ファイルを読み込む
> car<-read.csv("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.5 Rを用いてcar$incomeの平均を求めなさい。
課題 2.6 (復習)以下の計算を、Rのスクリプトchap2.rに書き、Rのコードを実行し
て求めなさい。結果が正しいかどうか確認すること。
3×5 34
56÷7
√
2.8
関数の作成
これまで、見てきたようにRでは、関数の利用が重要である。実は、新しい関数を自分で作
成することができる。
関数定義の基本形式は以下である。
関数名 <- function(引数){
処理の一行目
処理の二行目
・・・・・
return(計算結果を返す)
}
この章で用いた関数
関数名 関数の意味
q Rを終了する
sqrt 平方根を求める
c 複数の数値からベクトルを作成する
mean データの平均を求める
第
3
章
1
変数のデータ分析
まず、apartment.csvをRに読み込み、apartmentという名前のデータフレームを作ります。
> apartment<-read.csv("apartment.csv")
> 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 は家賃(万円)を意味します。
3.1
ヒストグラム
ここでは、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
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.2: 日経平均収益率のヒストグラム
経済・経営データでは,分布が左右対称でないものも多い.図3.1 は,平成22年の23区住
居地域の公示価格データ(単位は円/m2)のヒストグラムである.20−80万円の範囲に多く
のデータがあるが,100万円を超えるデータもそれなりにあり,右側に伸びている.このよう
な分布は右に歪んだ分布という.経済データでは分布が右に歪んでいることがよくある(所得,
貯蓄など).逆に左側にすそが伸びている場合,左に歪んだ分布という.
Histogram of kouji_23ku$kakaku
kouji_23ku$kakaku
Frequency
0 500000 1000000 2000000 3000000
0
50
100
150
200
250
図 3.3: 23区住居地域の公示価格データのヒストグラム
3.2
データの代表値−平均、メディアン−
ヒストグラムよりもさらにデータを要約して,データの特徴を表すいくつかの数値を求める
ことも多い.まず,1つの数値でデータを代表させようとする場合,これをデータの代表値と
いう.まず,データを
x1, x2. . . . , xn
で表そう.
3.2.1
平均
データx1, x2. . . . , xnの平均は,
¯
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万
円よりも小さい.平均が1つの大きな値1580万円に引っ張られてしまっていることがわかる.
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)
=1
n
(x1−x¯)2+ (x2 −x¯)2+· · ·+ (xn−x¯)2 (3.2)
で定義される。分散の正の平方根
S =√S2
を標準偏差(standard deviation)という.
(xi−x¯)は偏差と呼ばれ,
Pn
i=1(xi−x¯)2は偏差2乗和と呼ばれる.
分散は
s2 = 1
n−1
n
X
i=1
(xi−x¯)2 (3.3)
で定義されることもある.n−1 で割る理由は後ほど説明される.s2を不偏分散ともいう.標
準偏差は
s=√s2
で定義される.Rの関数var(),sd()で得られるのはs2,s である.
関数meanはデータの平均を求める。
> mean(apartment$area) [1] 29.43793
> var(apartment$area) [1] 229.4022
分散の正の平方根S =
√
S2は標準偏差と呼ばれ,Rでは関数sdで求められる。
> sd(apartment$area) [1] 15.14603
関数名 関数の意味
hist データのヒストグラムを描く
mean データの平均を求める
var データの分散を求める
第
4
章
2
変数のデータ分析
Rのコードはchap4.rに入力して、保存してください。
4.1
散布図
> 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
相関係数
相関関係を数値で測定する尺度として相関係数がある。データが(x1, y1),(x2, y2), . . . ,(xn, yn)
で与えられている時、x と y の相関係数は
r=
n
P
i=1
(xi−x¯)(yi−y¯)
r n P
i=1
(xi−x¯)2 n
P
i=1
(yi−y¯)2
(4.1)
で定義される。どのようなデータについても−1 ≤ r ≤ 1である。Rでは関数corで相関係数
を求められる。
> cor(apartment$yachin,apartment$area) [1] 0.8592066
>
4.3
散布図にラベルを貼る
cv2<-read.csv("cv2.csv")head(cv2)
xと yはベクトルであるとする.zは文字列の配列とする.text(x,y,labels=z,pos=3)は,
x座標,y座標がx, yで指定される点にzのラベルを描くことを意味する.散布図の点にラベ
ルを描く.
plot(conv~pop,data=cv2)
text(cv2$pop,cv2$conv,labels=cv2$todofuken,pos=3)
図 4.3: 都道府県別コンビニ数と人口
課題(提出)
cv2のデータを用いて一人当たりgppとdivorceの散布図を描きなさい.県名のラベルを貼るこ
第
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
(5.3)
a= ¯y−bx¯ (5.4)
結果は以下のようになる。
¯
y= 8, x¯= 20
b= 146
460 = 0.3174, a= 8−0.3174×20 = 1.652
得られた a, b を使った式
y=a+bx (5.5)
付録:最小
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 は定数である。観測点 (xi, yi)
から垂直に下ろした線が直線y=a+bx とぶつかる点は(xi, a+bxi) である。したがってその
2点の差はyi−a−bxi となる。そこで、直線y =a+bx がデータ(x1, y1),(x2, y2), . . . ,(xn, yn)
からどの程度「ずれ」ているかを測る尺度として、「ずれ」の2乗の和
S(a, b) =
n
X
i=1
(yi−a−bxi)2
を考える。S(a, b) を最小にするa, b の値を求める方法を最小二乗法という。
決定係数
最小2乗法を使って推定された回帰モデルを
y=a+bx (5.6)
と書こう。この回帰モデルを用いて、 xi の値から計算される y の値をyˆi と表す。すなわち
ˆ
yi =a+bxi, i= 1,2, . . . , n. (5.7)
である。yˆi を yi に対する回帰値という。当然、実際の観測値yi と回帰値yˆi は一致するとは限
らない。yi と yˆi の差
ei =yi−yˆi, i= 1,2, . . . , n (5.8)
を残差という。(5.8)式を変形すると
y
i
y
x y = α + β x
x
i
y
i
図 5.3: 回帰直線
と書ける。この式から、yˆi は yi の値のうち回帰モデルで説明できる部分、ei は回帰モデルで
説明できない部分である。データ全体で見たyi の変動のうち yˆi の変動で説明できる部分が大
きければよい。まず、(5.9)式の両辺から y¯を引く。
yi−y¯= (ˆyi−y¯) +ei (5.10)
次に、この式の両辺を2乗すると
(yi−y¯)2 = (ˆyi−y¯)2+ 2 (ˆyi−y¯)ei+e2i (5.11)
を得る。これを、すべての yi について合計すると以下のようになる。
n
X
i=1
(yi−y¯)2
| {z }
yiの変動 =
n
X
i=1
(ˆyi−y¯)2
| {z }
ˆ yiの変動
+
n
X
i=1
e2i
| {z }
残差の変動
(5.12)
(5.12)式は、y1, y2, . . . , yn 変動が、回帰モデルで説明できる部分と説明できない部分とに分解
できることを意味している。(5.12)式を、簡単に次のように書きなおす。
ST
|{z}
全変動
= SR
|{z}
回帰平方和
+ Se
|{z}
残差平方和
(5.13)
ST に対して SR の占める割合が大きければ大きいほど、推定された回帰モデルのデータに対す
る当てはまりがいいことを意味する。この割合をR2 で表す。
R2 = SR
ST
= ST −Se
ST
= 1− Se
ST
(5.14)
このR2 を推定された回帰モデル(5.6)の決定係数という。ST, SR,Se はすべて正だから、
0≦R2 ≦1
5.2
回帰係数に関する推測
以上の節では、回帰直線を求めて、その当てはまりの良さ(決定係数)について説明した。
得られたデータは、いろいろな誤差や不確実性を含むので、得られた回帰式の係数は絶対的な
ものではない。より進んだ分析手法では、データの不確実性を考慮した上で、得られた回帰式
から何が結論できるかを調べる。
yi, xi それぞれは被説明変数、説明変数のi番目の観測値とする。yi は xi と誤差項 ui とで
以下のように決定されるものと仮定する。
yi =α+βxi+ui, i= 1,2, . . . , n (5.15)
ただしui, i = 1,2, . . . , nは独立で、それぞれ正規分布 N(0, σ2) に従う。またα, βは未知の
定数である。このモデルを単回帰モデルという。β の最小2乗推定量をβˆで表す。上の仮定の
下で、
ˆ
β−β σβˆ ∼
N(0, 1) (5.16)
であることを証明することができる。ただしσˆ
β はβˆの標準偏差で,
σβˆ =
σ
s n X
i=1
(xi−x¯)2
(5.17)
で与えられる.(5.16)式左辺のσˆ
βの値は未知なので,その推定量で置き換える必要がある.誤
差項の分散σ2の推定量は残差e1, . . . , en を用いて
ˆ
σ2 = 1
n−2
n
X
i=1
e2i (5.18)
で与えられる.ˆσ2を用いれば,βˆの標準偏差の推定量は以下のようになる.
se( ˆβ) = s n σˆ X
i=1
(xi−x¯)2
(5.19)
この推定量はβˆの標準誤差(standard error)とよばる.このとき
ˆ
β−β
se( ˆβ) ∼t(n−2) (5.20)
が成り立つ.ただし,t(k) は自由度kのt分布を表す.この事実は,βに関する仮説検定に利
用される.
5.2.1
β
についての仮説検定
回帰分析では、以下のような帰無仮説の仮説検定が重要である。
対立仮説としては両側対立仮説 H1 : β 6= 0 を設定する場合は両側検定と言われ、片側対立仮 説H1 : β > 0 あるいは H1 : β < 0を設定する場合は片側検定という。この仮説検定が重要な
のは、帰無仮説(5.21)が棄却されれば、β 6= 0 つまり、得られた回帰式は意味のあるものと結
論できるからである。この仮説検定の検定統計量は
t= βˆ
se( ˆβ) (5.22)
である.両側検定の場合,この値の絶対値が十分大きければ、帰無仮説(5.21)が間違っている
証拠と考えられる.自由度n−2 の t分布の上側 100α%点をtα(n−2)で表すものとする.t∗
をt検定統計量の実現値としよう.t∗はt値とよばれる.もし
|t∗|> t
α/2(n−2) (5.23)
であれば,有意水準αでH0を棄却する.(5.23)式が成り立たなければ,H0は有意水準αでは
棄却されない.
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専有面積は有意である。
課題(メールで提出,件名は東大1027氏名とする.10月27日(木)に出すこと.)
教科書3.7節の操作を36ページablineまで行いなさい.ただし、被説明変数をgsmとし,説
明変数をpopとすること.
1.回帰分析の結果を書く.その際,教科書35ページのように被説明変数にはハットを付け
る.回帰式、t値、 決定係数R2を書き、解釈を書きなさい。(ワードで書く)
5.3
見せかけの相関
5.4
重回帰分析
前節までは、説明変数が一つの場合の回帰分析について述べてきた。説明変数が二つ以上あ
る場合の分析を重回帰分析という。被説明変数がy,説明変数がk個ありx1, x2, . . . , xkで表す
とする.これらの変数のi番目の個体についての観測値を(yi, x1i, x2i, . . . , xki)で表すとする.以
下のモデルを考える.
yi =α+β1x1i+β2x2i+· · ·+βkxki+ui, (i= 1,2, . . . , n) (5.24)
このモデルを重回帰モデルという.
apartmentのデータで、家賃を被説明変数、面積、築年数、時間距離の三つを説明変数とし
て重回帰分析を行ってみよう。
> apartment<-read.csv("apartment.csv")
> regresult<-lm(yachin~area+years+time,data=apartment) > summary(regresult)
Call:
lm(formula = yachin ~ area + years + time, data = apartment)
Residuals:
Min 1Q Median 3Q Max -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=6.21 + 0.21area −0.114years −0.164time
(14.1) (28.2) (−8.41) (−5.51)
課題(レポート準備)(メール本文に書く.件名は東大(レポート準備)氏名とする.締切り11
月24日(木)18時.)
レポート(提出締切りは1月上旬)では、各自が用意したデータを統計的手法で分析する。基本
は重回帰分析とする(研究に関連して分析したいデータがある人は他の分析方法でよい)。(1)
レポートで分析を行う計画のデータ(分析対象とする変数名)を具体的に述べなさい.データ
のソース(省庁名、統計の名前、サイトの名前.実験やアンケート調査の場合は、その概要な
ど)を書くこと。三つ以上の変数(性別や実験条件などカテゴリカルな変数も含めてよい)を
扱うこと.回帰分析の場合は,どの変数を被説明変数とするかも書くこと.
(2)(1)であげたデータをどのような目的で分析するか説明しなさい。
注意:この課題は、レポートの分析の準備を行うことが目的です。実際のレポート作成時に使
用データが変更になっても構いません.どのよう分析手法を用いたらよいのか不明であるとき
には、口頭またはメールで質問して下さい。
課題(ワードに書く.メール件名は,東大1117氏名とする.締切り11月17日(木)18時.)
1. todofuken_new.csvを授業のサイトからダウンロードし、Zドライブのフォルダrdataに保
存する.
2. Rでこのデータを読み込み,データフレームにtodofuken_newと名前を付ける.
3. crimerate2006を被説明変数、unempとpopdensityを説明変数として回帰分析を行う.ただ
し,crimerate2006は2006年の犯罪率,unempは完全失業率,popdensityは人口密度である。
4.教科書47ページ(c)と同様の形式で、回帰式、t値、自由度調整済み決定係数を書く.分析
結果を「∼は水準∼で有意である.∼は有意ではない」、のように説明する.
5.4.1
変数を変換したものを説明変数とする場合
plot(crime2006/pop~I(pop/area),data=todofuken)
crimereg<-lm(crime2006/pop~unemp+I(pop/area),data=todofuken) summary(crimereg)
5.5
ダミー変数
教科書43ページ、4.6節を読んでください.
性別、雇用されているかどうか、既婚かどうかなどの質的情報を説明変数として用いる時は、
ダミー変数を使う.カテゴリーがk個からなる情報のときにはk−1個のダミー変数を定義す
ればよい.たとえば、学歴として中学校卒、高校卒、大学卒、大学院卒の四つのカテゴリーの
情報を使いたいときには、これらのカテゴリーのうち三つのカテゴリーについて(たとえば高
校卒、大学卒、大学院卒)ダミー変数を作り,これらを説明変数に用いる.
課題(ワードに書く.メール件名は,東大1201氏名とする.締切り12月1日(木)18時.)
1. 教科書のサイトhttps://sites.google.com/site/econometricsr/
のコード、データダウンロードのページを見る.
2. denryoku.csv と6章_回帰分析(時系列データ).rをダウンロードし、Zドライブのフォ
ルダrdataに保存する.
3. スクリプト6章_回帰分析(時系列データ).rのコードを一行ずつ実行する.ただし,
# 残差を保存する
# 残差をプロットする
の部分はとばす.時系列プロットをワードに貼り付ける.
曜日祝日効果の回帰分析の所まで実行したら終了.教科書86ページの結果を同じことを確認す
る.回帰分析の結果(回帰式,t値,決定係数,各説明変数が有意かどうか.)を書き、メール
で提出.
4. (csvファイル作成の練習)各自がレポートで使用するデータをネット等から、コピーし、
csv ファイルに保存する.変数は2つ入れる.以下のことに注意.
(1)データと変数名のみをcsvファイルに貼りつける.
(2)変数名は、英小文字で.短めでわかりやすい名前が望ましい.
(3)数値にはカンマ区切り入れない
できたファイルも添付してメール.
12月8日課題以下の作業を自分のPCで行う。
1.以下のサイトに書いてある手順で、USBメモリにRをインストールする.
http://humansci.let.hokudai.ac.jp/m/terao/stat/r/usb_install/usb_install.html
(このサイトの手順に逐次したがってインストールする.)
2. USBメモリのRgui.exe をダブルクリックし、Rが起動することを確かめる.
3.パッケージlmtest, FinTS, varsをインストールする.
手順
パッケージ→パッケージのインストールをクリック
サイトを選択(Japanのサイトを選ぶ。)
lmtestを選択
インストール終了まで待つ
他の二つも同様に行う.