![[PukiWiki]](image/econwiki.png "[PukiWiki]")
|このページについて| ''概要'':このページでは、Rを使って計量ミクロする方法を解説します。計量ミクロ（ミクロ計量）については、北村先生が書かれた[[「ミクロ計量経済学入門」:http://www.ier.hit-u.ac.jp/~kitamura/PDF/A217.pdf]]や[[「ミクロ計量経済学とは何か」:http://www.ier.hit-u.ac.jp/~kitamura/PDF/A214.pdf]]を参照してください（そのうち加筆します）。 ''親ページ'':このページの親ページは、[[Rを使って計量経済分析]]です。 |目次| #contents ---- *Rawデータの整形 [#zd164cc4] **Excel 2007/Excel 2008 [#w81225f3] ミクロデータを利用しようとする場合、中にはすでに加工された状態で、かつSTATAやSPSS等ですぐに利用できるようバイナリファイルによって提供されているものもありますが（例えば、LSMS、DHS、PSID）、固定長のレコードデータを自分で加工しなければならない場合も少なくありません。 ほんの数年前まで、レコードが6万を超える大きなデータになるとExcelでは扱うことができなかったため、Perl、Ruby、Pythonといったスクリプト言語でテキストデータを切り分ける、もしくはレコード長のデータの加工に対応しているSASを利用するしかありませんでした。 しかし、Excel 2007（Macは2008）ではシートのサイズが最大で 104万8576行 × 1万6384列 と大きく改善されたこともあり、簡単に固定長のレコードに保存されたミクロデータを扱うことができるようになりました。 COLOR(RED){追記: 5万行×1000列あまりのデータを読み込んだところ、3万行までしか読み込めないというバグがあることが分かりました。特にエラーを表示するわけでもなく、本当にたちが悪いです。} *Rでミクロ計量分析を行う際のデータ処理 [#b05e228f] 以下では、Rでミクロ計量分析を行う際に必要不可欠なデータ処理について説明します。 **Mergeについて [#l5b6b07d] sexとincomeという二つのデータフレームがあるとし、この二つのデータを紐付けして一つのデータフレームを作成したいとします。ここでIDは一意であるとします。 |ID|sex|h |1|Male| |2|Female| |3|Male| ~ |ID|income|h |1|520| |2|430| |4|730| ***共通するものだけを残す [#b8ad211d] > merge(sex,income,by="ID") |ID|sex|income|h |1|Male|520| |2|Female|430| ***全てを残す [#uc846344] 引数に ''all=TRUE'' を指定 > merge(sex,income,by="ID",all=TURE) |ID|sex|income|h |1|Male|520| |2|Female|430| |3|Male|COLOR(red):NA| |4|COLOR(red):NA|730| ***左結合 [#m009b9f1] 引数に ''all.x=TRUE'' を指定 > merge(sex, income,by="ID", all.x=TURE) |ID|sex|income|h |1|Male|520| |2|Female|430| |3|Male|COLOR(red):NA| ***右結合 [#he747f9e] 引数に ''all.y=TRUE'' を指定 > merge(sex, income,by="ID", all.y=TURE) |ID|sex|income|h |1|Male|520| |2|Female|430| |4|COLOR(red):NA|730| **ダミー変数の作成 [#ya134c3c] Rでは説明変数に因子が含まれている場合、因子を自動でダミー変数に変換してくれるのでパッケージを利用するだけであればそのままでもかまいません。しかし、計量をSTATAやEviewsといった裏で何を計算しているのかわからない怪しいソフトウェアでなく、わざわざRでやろうとするみなさんの中には行列演算を行う場合も少なくないと思います。そんな場合は自分でダミー変数を作成する必要があります。ここでsexを、''Male''と''Female''という因子をもつベクトルとします。 > levels(sex) [1] "Female" "Male" そしてこのsexという因子ベクトルからfemaleというダミー変数を作りたいとします。 ダミー変数を作成する際は、''ifelse''という関数を使います。R使いなら、間違ってもforで回してifで条件分岐をしようなんて考えてはいけません。ifelse関数は、ifelse(test, yes, no)となっており、testがTUREならyesを、FALSEならnoとなるようなベクトルなり行列を返してくれる関数です。この場合、testにはFemaleならTRUEを、それ以外ならFALSEとなる条件式を入れ、yesに1を、noに0を入れればいいことになります。 > female <- ifelse(sex=="Female",1,0) **Appendについて [#o3b04512] **欠損値 (NA)について [#h3a53568] *R＋Stataでミクロ計量分析 [#m0eee53d] **分析をはじめる前に [#g9ce678e] 先にも述べましたが、PSIDやLSMSなどのデータは".csv"などの形式だけではなく、STATAのバイナリファイル（以下、".dta"）でも提供されています。そこで、ここでは".dta"をRで読んで分析するための方法を紹介します。 **".dta"をRに読み込ませよう [#rde9dba3] パッケージの読み込み install.packages("foreign") # foreignをインストール library(foreign) # Rで.dtaを読むためのパッケージを指定 データの読み込み data<-read.dta("file name.dta") # file nameを変更してください fix(data) # Stataの"browse"コマンドに相当。ただし、マイクロデータの読み込みには相応の時間がかかるので注意。 **ヘックマンのρの導出プログラム [[Written by Akihiko NODA:http://at-noda.com]] [#m1b1baf0] Stata10ではヘックマンのρを直接推定していません（Stata10 Reference A-H pp.557を参照）。よって、論文に掲載する際にはちょっとした数値計算が必要となります。以下では、デルタ法（J.Wooldridge(2002) pp.44を参照）を用いてヘックマンのρを計算します。 x<- # insert the "coef" of "/athrho" y<- # insert the "sd" of "/athrho" rho<-(exp(x*2)-1)/(1+exp(x*2)) rhosd<-((exp(2*x)*2/(1+exp(2*x))-(exp(2*x)-1)*(exp(2*x)*2)/(1+exp(2*x))^2)*y^2*(exp(2*x)*2/(1+exp(2*x))-(exp(2*x)-1)*(exp(2*x)*2)/(1+exp(2*x))^2))^(1/2) rho rhosd **欠損値について（Stataと大きく異なる点） [#oeadafa8] *ミクロ計量のモデル [#u79bd7f2] **ロジットモデル（Logit Model），プロビットモデル（Probit Model） [#be87ed2d] パッケージの読み込み install.packages("multcomp") library("multcomp") データの読み込み mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/dae/logit.csv")) attach(mydata) glmの引数としてfamily=binomial(link="logit")を指定。プロビットの場合は引数としてfamily=binomial(link="probit")を指定 result <- glm(admit~gre+gpa+topnotch, family=binomial(link="logit")) summary(result) Call: glm(formula = admit ~ gre + gpa + topnotch, family = binomial(link = "logit")) Deviance Residuals: Min 1Q Median 3Q Max -1.3905 -0.8836 -0.7137 1.2745 1.9572 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.600814 1.096379 -4.196 2.71e-05 *** gre 0.002477 0.001070 2.314 0.0207 * gpa 0.667556 0.325259 2.052 0.0401 * topnotch 0.437224 0.291853 1.498 0.1341 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 499.98 on 399 degrees of freedom Residual deviance: 478.13 on 396 degrees of freedom AIC: 486.13 Number of Fisher Scoring iterations: 4 尤度比検定 anova(result, test="Chisq") **多項ロジットモデル（Multinomial Logit Model） [#eb3830c7] mlogit{mlogit} multinomial{VGAM} mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/dae/mlogit.csv")) attach(mydata) library(VGAM) resl <- vglm(brand~female+age, family=multinomial()) summary(resl) **多項プロビットモデル（Multinomial Probit Model） [#l2429e72] mnp{MNP} mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/dae/mlogit.csv")) attach(mydata) library(MNP) resl <- mnp(brand~female+age) summary(resl) **順序ロジットモデル（Ordered Logit Model） [#q18f041b] lrm{Design} library(Design) mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/dae/ologit.csv")) ddist <- datadist(mydata[2:ncol(mydata)]) options(datadist="ddist") fit <- lrm(apply~.,data=mydata) fit Logistic Regression Model lrm(formula = apply ~ ., data = mydata) Frequencies of Responses 0 1 2 220 140 40 Obs Max Deriv Model L.R. d.f. P C Dxy Gamma Tau-a R2 Brier 400 3e-11 24.18 3 0 0.605 0.21 0.212 0.119 0.07 0.235 Coef S.E. Wald Z P y>=1 -2.20332 0.7795 -2.83 0.0047 y>=2 -4.29877 0.8043 -5.34 0.0000 pared 1.04766 0.2658 3.94 0.0001 public -0.05868 0.2979 -0.20 0.8438 gpa 0.61575 0.2606 2.36 0.0182 # オッズ比のグラフ表示 > sum.fit <- summary(fit) > plot(sum.fit) #ref(ologit.jpg,left,nowrap,50%,添付ファイルの画像) > newdata1<-data.frame(pared=c(0,1), public=mean(mydata$public), gpa=mean(mydata$gpa)) > newdata1$predicted<-predict(ologit,newdata=newdata1,type="response") **一般化順序ロジットモデル（Generalized Ordered Logit Model） [#q18f041b] **トービットモデル（Tobit Model）[#f0fa0d66] *** 打ち切りデータ（Truncated Data） [#w89bd93e] tobit{VGAM} n <- 1000 x <- seq(-1, 1, len=n) f <- function(x) 1 + 4*x ystar <- f(x) + rnorm(n) y = pmax(ystar, 0) fit <- vglm(y ~ x, tobit(Lower=0), trace=TRUE) summary(fit) plot(x, y, main="Tobit model", las=1) legend(-0.9, 3, c("Truth", "Estimate"), col=c("Blue", "Red"), lwd=2) lines(x, f(x), col="blue", lwd=2) lines(x, fitted(fit), col="red", lwd=2, lty="dashed") #ref(tobit.jpg,left,nowrap,50%,添付ファイルの画像) *** 検閲データ（Censored Data） [#dccb382f] cnormal1{VGAM} n = 1000 x = runif(n) ystar = rnorm(n, mean=100 + 15 * x, sd=exp(3)) # True values ## Not run: hist(ystar) L = runif(n, 80, 90) # Lower censoring points U = runif(n, 130, 140) # Upper censoring points y = pmax(L, ystar) # Left censored y = pmin(U, y) # Right censored extra = list(leftcensored = ystar < L, rightcensored = ystar > U) fit = vglm(y ~ x, cnormal1(zero=2), trace=TRUE, extra=extra) **ヘックマンの２段階推定 [#a9aaf328] **カウントデータ [#x70cfbf1] **サンプル選択問題 [#kc48fd67] **パネル [#z128a1af] [[(；ﾟдﾟ) R言語でパネルデータ・モデルを推計せよ！:http://uncorrelated.no-ip.com/cgi-bin/view.cgi/20090414/T]] *References [#d6520dd2] **書籍 [#ae2b6356] -[[浅野・中村(2000):http://www.amazon.co.jp/gp/product/4641160805/sr=8-1/qid=1155206048/ref=sr_1_1/250-7742157-1058652?ie=UTF8&s=gateway]] -[[Greene(2003):http://www.amazon.co.jp/gp/product/0131108492/sr=8-1/qid=1155281794/ref=pd_bowtega_1/250-7742157-1058652?ie=UTF8&s=gateway]] -[[Hamilton(1994):http://www.amazon.co.jp/gp/product/0691042896/ref=pd_sim_gw_4/250-7742157-1058652?ie=UTF8]] -[[北村(2005):http://www.amazon.co.jp/gp/product/4000099116/250-7513282-3114617?v=glance&n=465392]] -[[Wooldridge(2002):http://www.amazon.co.jp/Econometric-Analysis-Cross-Section-Panel-Data/dp/0262232197/ref=pd_bxgy_b_img_b/503-7013025-2018353?ie=UTF8]] -[[講座 ミクロ統計分析 2 - ミクロ統計の集計解析と技法:http://www.nippyo.co.jp/book/1490.html]] **リンク [#vb5d7ae1] -[[Econometrics in R:http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf]] -[[CRAN Task View: Computational Econometrics:http://cran.md.tsukuba.ac.jp/src/contrib/Views/Econometrics.html]] -[[UCLAのサイト:http://www.ats.ucla.edu/stat/r/]] -[[パッケージ一覧:http://cran.md.tsukuba.ac.jp/src/contrib/PACKAGES.html]]
