このページについて |
概要:このページでは、Rを使って計量ミクロする方法を解説します。計量ミクロ(ミクロ計量)については、北村先生が書かれた「ミクロ計量経済学入門」や「ミクロ計量経済学とは何か」を参照してください(そのうち加筆します)。
親ページ:このページの親ページは、Rを使って計量経済分析です。
目次 |
ミクロデータを利用しようとする場合、中にはすでに加工された状態で、かつSTATAやSPSS等ですぐに利用できるようバイナリファイルによって提供されているものもありますが(例えば、LSMS、DHS、PSID)、固定長のレコードデータを自分で加工しなければならない場合も少なくありません。
ほんの数年前まで、レコードが6万を超える大きなデータになるとExcelでは扱うことができなかったため、Perl、Ruby、Pythonといったスクリプト言語でテキストデータを切り分ける、もしくはレコード長のデータの加工に対応しているSASを利用するしかありませんでした。
しかし、Excel 2007(Macは2008)ではシートのサイズが最大で 104万8576行 × 1万6384列 と大きく改善されたこともあり、簡単に固定長のレコードに保存されたミクロデータを扱うことができるようになりました。
追記: 5万行×1000列あまりのデータを読み込んだところ、3万行までしか読み込めないというバグがあることが分かりました。特にエラーを表示するわけでもなく、本当にたちが悪いです。
以下では、Rでミクロ計量分析を行う際に必要不可欠なデータ処理について説明します。
sexとincomeという二つのデータフレームがあるとし、この二つのデータを紐付けして一つのデータフレームを作成したいとします。ここでIDは一意であるとします。
ID | sex |
1 | Male |
2 | Female |
3 | Male |
ID | income |
1 | 520 |
2 | 430 |
4 | 730 |
> merge(sex,income,by="ID")
ID | sex | income |
1 | Male | 520 |
2 | Female | 430 |
引数に all=TRUE を指定
> merge(sex,income,by="ID",all=TURE)
ID | sex | income |
1 | Male | 520 |
2 | Female | 430 |
3 | Male | NA |
4 | NA | 730 |
引数に all.x=TRUE を指定
> merge(sex, income,by="ID", all.x=TURE)
ID | sex | income |
1 | Male | 520 |
2 | Female | 430 |
3 | Male | NA |
引数に all.y=TRUE を指定
> merge(sex, income,by="ID", all.y=TURE)
ID | sex | income |
1 | Male | 520 |
2 | Female | 430 |
4 | NA | 730 |
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)
先にも述べましたが、PSIDやLSMSなどのデータは".csv"などの形式だけではなく、STATAのバイナリファイル(以下、".dta")でも提供されています。そこで、ここでは".dta"をRで読んで分析するための方法を紹介します。
パッケージの読み込み
install.packages("foreign") # foreignをインストール library(foreign) # Rで.dtaを読むためのパッケージを指定
データの読み込み
data<-read.dta("file name.dta") # file nameを変更してください fix(data) # Stataの"browse"コマンドに相当。ただし、マイクロデータの読み込みには相応の時間がかかるので注意。
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
パッケージの読み込み
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")
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)
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)
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)
> newdata1<-data.frame(pared=c(0,1), public=mean(mydata$public), gpa=mean(mydata$gpa)) > newdata1$predicted<-predict(ologit,newdata=newdata1,type="response")
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")
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)