このページについて

概要:このページでは、Rを使って計量ミクロする方法を解説します。計量ミクロ(ミクロ計量)については、北村先生が書かれた「ミクロ計量経済学入門」「ミクロ計量経済学とは何か」を参照してください(そのうち加筆します)。

親ページ:このページの親ページは、Rを使って計量経済分析です。

目次

Rawデータの整形

Excel 2007/Excel 2008

ミクロデータを利用しようとする場合、中にはすでに加工された状態で、かつSTATAやSPSS等ですぐに利用できるようバイナリファイルによって提供されているものもありますが(例えば、LSMS、DHS、PSID)、固定長のレコードデータを自分で加工しなければならない場合も少なくありません。

ほんの数年前まで、レコードが6万を超える大きなデータになるとExcelでは扱うことができなかったため、Perl、Ruby、Pythonといったスクリプト言語でテキストデータを切り分ける、もしくはレコード長のデータの加工に対応しているSASを利用するしかありませんでした。

しかし、Excel 2007(Macは2008)ではシートのサイズが最大で 104万8576行 × 1万6384列 と大きく改善されたこともあり、簡単に固定長のレコードに保存されたミクロデータを扱うことができるようになりました。

追記: 5万行×1000列あまりのデータを読み込んだところ、3万行までしか読み込めないというバグがあることが分かりました。特にエラーを表示するわけでもなく、本当にたちが悪いです。

Rでミクロ計量分析を行う際のデータ処理

以下では、Rでミクロ計量分析を行う際に必要不可欠なデータ処理について説明します。

Mergeについて

sexとincomeという二つのデータフレームがあるとし、この二つのデータを紐付けして一つのデータフレームを作成したいとします。ここでIDは一意であるとします。

IDsex
1Male
2Female
3Male


IDincome
1520
2430
4730

共通するものだけを残す

> merge(sex,income,by="ID")
IDsexincome
1Male520
2Female430

全てを残す

引数に all=TRUE を指定

> merge(sex,income,by="ID",all=TURE)
IDsexincome
1Male520
2Female430
3MaleNA
4NA730

左結合

引数に all.x=TRUE を指定

> merge(sex, income,by="ID", all.x=TURE)
IDsexincome
1Male520
2Female430
3MaleNA

右結合

引数に all.y=TRUE を指定

> merge(sex, income,by="ID", all.y=TURE)
IDsexincome
1Male520
2Female430
4NA730

ダミー変数の作成

Rでは説明変数に因子が含まれている場合、因子を自動でダミー変数に変換してくれるのでパッケージを利用するだけであればそのままでもかまいません。しかし、計量をSTATAやEviewsといった裏で何を計算しているのかわからない怪しいソフトウェアでなく、わざわざRでやろうとするみなさんの中には行列演算を行う場合も少なくないと思います。そんな場合は自分でダミー変数を作成する必要があります。ここでsexを、MaleFemaleという因子をもつベクトルとします。

> 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について

欠損値 (NA)について

R+Stataでミクロ計量分析

分析をはじめる前に

先にも述べましたが、PSIDやLSMSなどのデータは".csv"などの形式だけではなく、STATAのバイナリファイル(以下、".dta")でも提供されています。そこで、ここでは".dta"をRで読んで分析するための方法を紹介します。

".dta"をRに読み込ませよう

パッケージの読み込み

install.packages("foreign") # foreignをインストール
library(foreign) # Rで.dtaを読むためのパッケージを指定

データの読み込み

data<-read.dta("file name.dta") # file nameを変更してください
fix(data) # Stataの"browse"コマンドに相当。ただし、マイクロデータの読み込みには相応の時間がかかるので注意。

ヘックマンのρの導出プログラム Written by Akihiko NODA

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と大きく異なる点)

ミクロ計量のモデル

ロジットモデル(Logit Model),プロビットモデル(Probit Model)

パッケージの読み込み

 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)

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)

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)

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")

一般化順序ロジットモデル(Generalized Ordered Logit Model)

トービットモデル(Tobit Model)

打ち切りデータ(Truncated Data)

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")
添付ファイルの画像

検閲データ(Censored Data)

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)

ヘックマンの2段階推定

カウントデータ

サンプル選択問題

パネル

(;゚д゚) R言語でパネルデータ・モデルを推計せよ!

References

書籍

リンク


添付ファイル: fileologit.jpg 3676件 [詳細] filetobit.jpg 4327件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2013-10-27 (日) 23:35:36