Rで計量ミクロ
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
|
ログイン
]
開始行:
|このページについて|
''概要'':このページでは、Rを使って計量ミクロする方法を解...
''親ページ'':このページの親ページは、[[Rを使って計量経済...
|目次|
#contents
----
*Rawデータの整形 [#zd164cc4]
**Excel 2007/Excel 2008 [#w81225f3]
ミクロデータを利用しようとする場合、中にはすでに加工され...
ほんの数年前まで、レコードが6万を超える大きなデータになる...
しかし、Excel 2007（Macは2008）ではシートのサイズが最大で...
COLOR(RED){追記: 5万行×1000列あまりのデータを読み込んだと...
*Rでミクロ計量分析を行う際のデータ処理 [#b05e228f]
以下では、Rでミクロ計量分析を行う際に必要不可欠なデータ処...
**Mergeについて [#l5b6b07d]
sexとincomeという二つのデータフレームがあるとし、この二つ...
|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では説明変数に因子が含まれている場合、因子を自動でダミー...
> levels(sex)
[1] "Female" "Male"
そしてこのsexという因子ベクトルからfemaleというダミー変数...
ダミー変数を作成する際は、''ifelse''という関数を使います...
> female <- ifelse(sex=="Female",1,0)
**Appendについて [#o3b04512]
**欠損値 (NA)について [#h3a53568]
*R＋Stataでミクロ計量分析 [#m0eee53d]
**分析をはじめる前に [#g9ce678e]
先にも述べましたが、PSIDやLSMSなどのデータは".csv"などの...
**".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...
Stata10ではヘックマンのρを直接推定していません（Stata10 R...
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)*...
rho
rhosd
**欠損値について（Stataと大きく異なる点） [#oeadafa8]
*ミクロ計量のモデル [#u79bd7f2]
**ロジットモデル（Logit Model），プロビットモデル（Probit...
パッケージの読み込み
install.packages("multcomp")
library("multcomp")
データの読み込み
mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/...
attach(mydata)
glmの引数としてfamily=binomial(link="logit")を指定。プロ...
result <- glm(admit~gre+gpa+topnotch, family=binomial(l...
summary(result)
Call:
glm(formula = admit ~ gre + gpa + topnotch, family = bin...
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...
(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/...
attach(mydata)
library(VGAM)
resl <- vglm(brand~female+age, family=multinomial())
summary(resl)
**多項プロビットモデル（Multinomial Probit Model） [#l242...
mnp{MNP}
mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/...
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/...
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 ...
400 3e-11 24.18 3 0 ...
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$...
> newdata1$predicted<-predict(ologit,newdata=newdata1,ty...
**一般化順序ロジットモデル（Generalized Ordered Logit Mod...
**トービットモデル（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", "R...
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 v...
## 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 = y...
fit = vglm(y ~ x, cnormal1(zero=2), trace=TRUE, extra=ex...
**ヘックマンの２段階推定 [#a9aaf328]
**カウントデータ [#x70cfbf1]
**サンプル選択問題 [#kc48fd67]
**パネル [#z128a1af]
[[(；ﾟдﾟ) R言語でパネルデータ・モデルを推計せよ！:http://...
*References [#d6520dd2]
**書籍 [#ae2b6356]
-[[浅野・中村(2000):http://www.amazon.co.jp/gp/product/46...
-[[Greene(2003):http://www.amazon.co.jp/gp/product/013110...
-[[Hamilton(1994):http://www.amazon.co.jp/gp/product/0691...
-[[北村(2005):http://www.amazon.co.jp/gp/product/40000991...
-[[Wooldridge(2002):http://www.amazon.co.jp/Econometric-A...
-[[講座 ミクロ統計分析 2 - ミクロ統計の集計解析と技法:htt...
**リンク [#vb5d7ae1]
-[[Econometrics in R:http://cran.r-project.org/doc/contri...
-[[CRAN Task View: Computational Econometrics:http://cran...
-[[UCLAのサイト:http://www.ats.ucla.edu/stat/r/]]
-[[パッケージ一覧:http://cran.md.tsukuba.ac.jp/src/contri...
終了行:
|このページについて|
''概要'':このページでは、Rを使って計量ミクロする方法を解...
''親ページ'':このページの親ページは、[[Rを使って計量経済...
|目次|
#contents
----
*Rawデータの整形 [#zd164cc4]
**Excel 2007/Excel 2008 [#w81225f3]
ミクロデータを利用しようとする場合、中にはすでに加工され...
ほんの数年前まで、レコードが6万を超える大きなデータになる...
しかし、Excel 2007（Macは2008）ではシートのサイズが最大で...
COLOR(RED){追記: 5万行×1000列あまりのデータを読み込んだと...
*Rでミクロ計量分析を行う際のデータ処理 [#b05e228f]
以下では、Rでミクロ計量分析を行う際に必要不可欠なデータ処...
**Mergeについて [#l5b6b07d]
sexとincomeという二つのデータフレームがあるとし、この二つ...
|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では説明変数に因子が含まれている場合、因子を自動でダミー...
> levels(sex)
[1] "Female" "Male"
そしてこのsexという因子ベクトルからfemaleというダミー変数...
ダミー変数を作成する際は、''ifelse''という関数を使います...
> female <- ifelse(sex=="Female",1,0)
**Appendについて [#o3b04512]
**欠損値 (NA)について [#h3a53568]
*R＋Stataでミクロ計量分析 [#m0eee53d]
**分析をはじめる前に [#g9ce678e]
先にも述べましたが、PSIDやLSMSなどのデータは".csv"などの...
**".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...
Stata10ではヘックマンのρを直接推定していません（Stata10 R...
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)*...
rho
rhosd
**欠損値について（Stataと大きく異なる点） [#oeadafa8]
*ミクロ計量のモデル [#u79bd7f2]
**ロジットモデル（Logit Model），プロビットモデル（Probit...
パッケージの読み込み
install.packages("multcomp")
library("multcomp")
データの読み込み
mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/...
attach(mydata)
glmの引数としてfamily=binomial(link="logit")を指定。プロ...
result <- glm(admit~gre+gpa+topnotch, family=binomial(l...
summary(result)
Call:
glm(formula = admit ~ gre + gpa + topnotch, family = bin...
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...
(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/...
attach(mydata)
library(VGAM)
resl <- vglm(brand~female+age, family=multinomial())
summary(resl)
**多項プロビットモデル（Multinomial Probit Model） [#l242...
mnp{MNP}
mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/...
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/...
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 ...
400 3e-11 24.18 3 0 ...
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$...
> newdata1$predicted<-predict(ologit,newdata=newdata1,ty...
**一般化順序ロジットモデル（Generalized Ordered Logit Mod...
**トービットモデル（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", "R...
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 v...
## 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 = y...
fit = vglm(y ~ x, cnormal1(zero=2), trace=TRUE, extra=ex...
**ヘックマンの２段階推定 [#a9aaf328]
**カウントデータ [#x70cfbf1]
**サンプル選択問題 [#kc48fd67]
**パネル [#z128a1af]
[[(；ﾟдﾟ) R言語でパネルデータ・モデルを推計せよ！:http://...
*References [#d6520dd2]
**書籍 [#ae2b6356]
-[[浅野・中村(2000):http://www.amazon.co.jp/gp/product/46...
-[[Greene(2003):http://www.amazon.co.jp/gp/product/013110...
-[[Hamilton(1994):http://www.amazon.co.jp/gp/product/0691...
-[[北村(2005):http://www.amazon.co.jp/gp/product/40000991...
-[[Wooldridge(2002):http://www.amazon.co.jp/Econometric-A...
-[[講座 ミクロ統計分析 2 - ミクロ統計の集計解析と技法:htt...
**リンク [#vb5d7ae1]
-[[Econometrics in R:http://cran.r-project.org/doc/contri...
-[[CRAN Task View: Computational Econometrics:http://cran...
-[[UCLAのサイト:http://www.ats.ucla.edu/stat/r/]]
-[[パッケージ一覧:http://cran.md.tsukuba.ac.jp/src/contri...
ページ名:
![[PukiWiki]](image/econwiki.png "[PukiWiki]")
