計量経済学講義 第 25 回 R による計量経済分析 Part-1 2018 年 1 月 5 日 ( 金 )1 限 担当教員 : 唐渡 広志 研究室 : 経済学研究棟 4 階 432 号室 email: kkarato@eco.u-toyama.ac.jp website: http://www3.u-toyama.ac.jp/kkarato/ 1
講義の目的 より高度な計量経済分析を行うために総合 的な統計分析ソフト R の基本的な使い方を, 学びます keywords: コンソール画面, エディタ画面, スクリプト, データ ベクトル, クロス集計表, 独立性の検定, 順位相関係 数, 平均の差の検定, マン ホイットニーの U 検定 2
R とは フリーの統計解析ソフトウェア ( プログラミング言語 ) 大学の総合情報基盤センターにインストール済み プラットホーム :Linux, OS X, Windows ホームページ The R Project for Statistical Computing ( https://www.rproject.org/ ) ダウンロードとインストール The Comprehensive R Archive Network ( https://cran.rproject.org/welcome.html ) 使い方 ( 参考になるサイト ) http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html http://www.okadajp.org/rwiki/ 3
コンソール画面 コンソール画面 ここにコマンド ( 命令 ) を打ち込む 4
演算子と関数 (1) 基本演算子 + 足し算 引き算 * 掛け算 / 割り算 ^ 累乗 例. 以下の計算を実行 2 5 95 12 5 1 3 4 3 1 0.5 5
演算子と関数 (2) 2 sqrt(2) 1 e exp(1) log1 log(1) log 3 5 log(5,3) 6
演算子と関数 (3) 関数と代入 練習問題 (1) > x<- 10 > x^2 [1] 100 [1] 1 2 3 2 [2] e 0.03 1 > b<- sqrt(2) > round(b,2) [3] ln105 ln100 [1] 1.41 [4] 2 0.5x ただし, x 4 <- は代入を表す 7
エディタ (1) R の ファイル メニューの 新しいスクリプト をクリック エディタが起動 エディタに計算内容やコマンドを記述しておくと, エディタからそのコードを実行することができる エディタ画面 8
エディタ (2) エディタに記述した内容は保存しておく 例. 以下のコードをエディタに記述して保存する ( 保存するのは計算結果ではなく計算内容 ) b<- sqrt(2) round(b,2) ( コンソール画面ではなく ) エディタ画面をアクティブにした状態で ファイル メニューの 別名で保存 を選択 ファイル名を 0105.r として保存する ( 拡張子は r ) 〇〇.r は R のソースコード ファイルである ただし, 保存する場所に注意 ローカルのドライブではなく, ネットワークドライブか自分の USB メモリに保存すること 9
エディタ (3) コマンドの実行 実行したい部分を選択状態にして, 2 行まとめて実行 右クリック カーソル行または選択中の R コードを実行 または, ctrl + R のショートカットを利用する コンソール画面に結果が表示される > b<- sqrt(2) > round(b,2) カーソル行だけ実行 [1] 1.41 10
データ ベクトル (1) エディタに以下を記述 x<- c(2,4,6,8) x y<- c(7,4,12,5) x+y # xとyの和 x*y # xとyの積 データ x の作成 データ データ y の作成 x i y i x i y i i x i i の表示 の表示 の表示 : xi : yi 2,4,6,8 7,4,12,5 # 記号以下にコメント文を書くことができる Rでは# の部分を無視する c(2,4,6,8) のようにいくつかの値をひとまとめにしたものをデータ ベクトルとよぶ また,<- は代入を表している したがって, x<- c(2,4,6,8) はデータ ベクトル c(2,4,6,8) を x という文字に代入している x+y や x*y のように 長さ の等しい複数のデータ ベクトルを用いてデータ間の演算ができる 11
データ ベクトル (2) 連番データ a<- 1:10 a [1] 1 2 3 4 5 6 7 8 9 10 seq(1,10,2) [1] 1 3 5 7 9 A:B A, Bが整数のときは AからBまでの連番データを作成 seq A,B,m A, Bが整数のときは AからBまでのm個おきのデータを作成 ベクトルの 長さ length(x) [1] 4 length(a) [1] 10 length x データベクトルの長さ, データ x i の数 12
データ ベクトル (3) ベクトルの要素 x[3] [1] 6 x[x<5] [1] 2 4 x[4]<- 10 x [1] 2 4 6 10 x [ k] : xのk 番目の要素 x [ x j] : x j を満たす x 4 番目の値を10に上書き の要素を表示 ベクトル演算 > sum(x) # 合計 [1] 22 > mean(x) # 平均 [1] 5.5 > var(x) # 分散 [1] 11.66667 13
R の統計関数 14
散布図と回帰直線 (1) plot(x,y) eq<- lm(y~x) abline(eq) eq Call: lm(formula = y ~ x) Coefficients: (Intercept) x plot ( x, y) : 横軸に x, 縦軸に y の散布図を作成 y ~ x は回帰モデル y x u を意味する y ~ x の記述を formula( 式 ) とよぶ lm y ~ x は回帰モデルを最小 2乗法で推定するコマンド lm : linear model 最小 2乗法の計算結果 ( 切片と傾き ) が eq に代入されている abline eq : 計算された切片と傾きに基づいて, 回帰直線を散布図に書き込む 7.31429-0.05714 summary(eq) summary (eq) : 最小 2 乗推定の詳細な結果を表示する 15
散布図と回帰直線 (2) summary(eq) summary (eq) : 最小 2 乗推定の詳細な結果を表示する Call: lm(formula = y ~ x) Residuals: 残差 1 2 3 4-0.200-3.086 5.029-1.743 Coefficients: 推定値標準誤差 t 値 p 値 Estimate Std. Error t value Pr(> t ) (Intercept) 7.31429 4.59432 1.592 0.252 x -0.05714 0.73568-0.078 0.945 Residual standard error: 4.352 on 2 degrees of freedom 回帰の標準誤差と残差 2 乗和の自由度 Multiple R-squared: 0.003008, Adjusted R-squared: -0.4955 F-statistic: 0.006033 on 1 and 2 DF, p-value: 0.9452 決定係数, 自由度調整済み決定係数 ゼロスロープ検定の F 値, 分子 分母の自由度,F 検定の p 値 16
外部データの読み込み (1) http://www3.u-toyama.ac.jp/kkarato/2017/econometrics/ から data.csv をダウンロード 〇〇〇.docx, 〇〇〇.xlsx, 〇〇〇.txt, 〇〇〇.png のようにファイル名 〇〇〇 の後ろについている docx, xlsx, txt, png は拡張子とよばれている 拡張子はファイルの種類を識別するための文字列と思えばよい たいていのウィンドウズ PC では拡張子によって開くアプリケーションを定義づけている (docx はワード,xlsx はエクセルなど ) パソコンによっては拡張子が表示されないように設定されている その場合, コントロール パネルの エクスプローラーのオプション の 表示 において 登録されている拡張子は表示しない のチェックをはずしておく csv は comma-separated values は, いくつかのフィールドを区切り文字であるカンマ, で区切ったテキストデータおよびテキストファイル エクセルで開こうと思えば開けるが, テキスト エディター ( 例えば秀丸 ) やメモ帳でも開くことができる csv や txt などのカンマ区切りテキストファイルを R で読み込むには, read.csv コマンドを使う data<- read.csv( data.csv ) 読み込んだ data のことをデータ フレーム (data.frame) とよぶ 17
外部データの読み込み (2) data.csv は今作成しているソースコード ファイル 0105.r と同じ場所に置く ローカルドライブの ドキュメント には保存しない ネットワークドライブまたは自分のUSBメモリの適当なフォルダに保存しておく データ読み込みのコマンド data<- read.csv( data.csv ) attach(data) n<- dim(data)[1] str(data) summary(data) read. csv : csv ファイルを読み込んで データフレームを追加するコマンド クォーテーション (shift dim はデータの次元 ( 観測値の数 2) で括る ファイル名 ( 拡張子付き ) をダブル 読み込んだ内容を data というオブジェクト に代入する ( オブジェクト名は何でもよい ) attach は変数へのアクセスを簡単にするコマンド dim[1] は観測値の数を示す観測値の数を nとする str はデータの構造を一覧表示するコマンド summary はデータの要約統計量を示すコマンド 変数の数 ) を調べるコマンド 18
回帰分析 例. 配偶者間の年齢 配偶者の年齢 = a + b 本人の年齢 + u age.sp age.sf sf: self 本人 sp: spouse 配偶者 plot(age.sf, age.sp) eq<- lm(age.sf~age.sp) summary(eq) 配偶者の年齢 = b 1 + b 2 本人の年齢 + b 3 18 歳未満の子供の数 + u age.sp age.sf c18 eq<- lm(age.sf~age.sp+c18) summary(eq) 19
データ集約 (1) データ閲覧と編集 fix(data) 閲覧が終わったら, データエディタ画面を閉じる ( をクリック ) データの表示 block 居住している 地域ブロック の番号が表示される データの分布 x<- table(block) x barplot(x) table : データの度数分布を表示するコマンド barplot : 棒グラフを描画するコマンド 20
データ集約 (2) クロス集計表と独立性の検定 検証する仮説 H 0 : 居住している地域と性別の分布は独立である H 1 : 居住している地域と性別の分布は独立でない x<- table(block,sex) x chisq.test(x) block 地域ブロック sex 性別男性女性 1 2 北海道 東北 1 146 123 関東 2 304 288 中部 3 248 236 近畿 4 176 180 中国 四国 5 98 103 九州 6 120 140 tableコマンドを 2 つの変数に使うとクロス集計表が作成できる 6 2のクロス集計表なので, 独立性の検定の自由度は df = (6-1) (2-1)=5 検定統計量 ( カイ2 乗値 ) は X-squared = 4.1808,p 値は 0.5237 なので有意水準 5% で帰無仮説 H 0 を棄却できない (p 値が0.05 未満のときH 0 を棄却 ) 居住している地域と性別の分布は独立でない とは言えない地域によって男女比に偏りがあるサンプルにはなっていない 自由度 5の カイ 2乗分布 0 5 10 15 20 カイ 2乗 p 値 0.5237 21
データ集約 (3) 使用データの尺度と関係性の分析 質的データ量的データ名義尺度順序尺度間隔尺度非尺度 名義尺度クロス集計クロス集計 ( クロス集計 ) ( クロス集計 ) 質的データ 順序尺度 クロス集計順位相関係数 ( クロス集計 ) 順位相関係数 ( クロス集計 ) 順位相関係数 量的データ 間隔尺度 非尺度 ( クロス集計 ) ( 順位相関係数 ) 相関係数 ( クロス集計 ) ( 順位相関係数 ) 相関係数 ( クロス集計 ) ( 順位相関係数 ) 相関係数 22
データ集約 (4) スピアマンの順位相関係数 cor(city,wllive) cor(city,wllive,method="spearman") cor.test(city,wllive,method="spearman") city wllive 1 2 3 4 1 274 191 41 4 2 290 187 50 7 3 521 259 78 15 4 139 78 27 1 cor (x,y) は変数 x, y の相関係数を計算するコマンド method= spearman を追加すると, 順位相関係数を計算する cor.test(x,y) は相関係数がゼロであるという帰無仮説を検定するコマンド Spearman's rank correlation rho data: city and wllive S = 1735600000, p-value = 0.1567 alternative hypothesis: true rho is not equal to 0 sample estimates: rho -0.03046662 23
データの制限 (1) income.sf[age.sf==40] income.sf[age.sf<40] income.sf[sex==1] income.sf[sex!=1] 年齢が40歳に等しい人の年収階級年齢が40歳未満の人の年収階級性別が男性である人の年収階級性別が男性でない人の年収階級 演算子を使った条件式の書き方変数 == a : a に等しい変数 >= a : a 以上変数 <= a : a 以下変数 > a : a より大きい変数 < a : a より小さい変数!= a : a に等しくない income.sf[age.sf>=30 & age.sf<40] income.sf[edu.sf<30 age.sf>=60] 論理演算子 a & b : a かつ b a b : a または b 30 40 50 60 70 年齢が30歳以上 かつ 40歳未満の人の年収階級年齢が30歳未満 または 60歳以上の人の年収階級 24
データの制限 (2) x<- income.sf[age.sf<40] y<- sex[age.sf<40] table(x,y) 年齢が 40 歳未満の人の 年収階級と性別のクロス集計 x<- income.sf[sex==2 & job.sf==2] y<- edu.sf[sex==2 & job.sf==2] table(x,y) 性別が女性 かつ 就労形態が 常時雇用の一般従業者の 年収階級と学歴のクロス集計 25
平均の差の検定 平均の差の検定 (Welch s t-test) x<- age.sf[ownhouse==1 & income.hh==2] y<- age.sf[ownhouse==1 & income.hh==4] t.test(x,y) 対象となる母集団 x 居住形態が持ち家かつ世帯年収が 350-550 万円未満の人の年齢 y 居住形態が持ち家かつ世帯年収が 850 万円以上の人の年齢検証する仮説 H 0 : x と y の年齢に差がない H 1 : x と y の年齢に差がある t.test : x と y の母集団分散が異なる ( 不等分散 ) という条件下での t 検定 p 値が 0.05 未満のとき H 0 を棄却 26
中央値の差の検定 マン ホイットニーの U 検定 ( ウィルコクスンの順位和検定 ) x<- job.cont[sex==1 & work.sf==1 & edu.sf==3] y<- job.cont[sex==1 & work.sf==1 & edu.sf==5] wilcox.test(x,y) 対象となる母集団 x 性別が男性かつ仕事をしており, 学歴が高校卒の人の就労継続意向 y 性別が男性かつ仕事をしており, 学歴が大学 大学院の人の就労継続意向検証する仮説 H 0 : x と y の就労継続意向に差がない H 1 : x と y の就労継続意向に差がある wilcox.test : マン ホイットニーの U 検定 ( ウィルコクスンの順位和検定 ) 標準正規分布での t 検定 p 値が 0.05 未満のとき H 0 を棄却 27