本章で使用する東大社研・若年パネル調査（JLPS-Y）wave1のオープンデータ（u001）は以下のURLで公開されていますので、取得してください。
https://csrda.iss.u-tokyo.ac.jp/infrastructure/urd/
分析結果の解釈については、書籍をご参照ください。以下はWindows10、R version 4.2.1で実行しています。
データを読み込むにあたって、作業場所の指定とそこにデータが置いてあることが前提になります。 R本体を操作している場合は、「ファイル」→「ディレクトリの変更」でデータの置いてある場所を指定するのが簡単です。 RStudioを操作している場合は、例えばデスクトップの「R」という名前のフォルダにデータがあるとすると
setwd("C:/Users/ユーザー名/Desktop/R")を最初に実行するのがよいでしょう。ユーザー名は各自異なるので注意です。
使用するデータを読み込み、表示します。なお、R version 4.1以前で読み込む際には
data1 <- read.csv("ファイル名.csv", fileEncoding = "UTF-8-BOM")のように、エンコードのオプションを指定する必要があります。version 4.2以降では必要ありませんので、以下ではオプションを指定していません。
data1 <- read.csv("u001.csv")本章で使用するパッケージを読み込みますが、パッケージはインストールしておく必要があります。 以下のようにコマンドでインストールするか、RやRStudioからクリックでインストールすることもできます。
install.packages(“stargazer”, dependencies = TRUE)パッケージの機能が使えるように読み込みます。
library(dplyr)##
## 次のパッケージを付け加えます: 'dplyr'## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter, lag## 以下のオブジェクトは 'package:base' からマスクされています:
##
## intersect, setdiff, setequal, unionlibrary(car)## 要求されたパッケージ carData をロード中です##
## 次のパッケージを付け加えます: 'car'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## recodelibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(mfx)## 要求されたパッケージ sandwich をロード中です## 要求されたパッケージ lmtest をロード中です## 要求されたパッケージ zoo をロード中です##
## 次のパッケージを付け加えます: 'zoo'## 以下のオブジェクトは 'package:base' からマスクされています:
##
## as.Date, as.Date.numeric## 要求されたパッケージ MASS をロード中です##
## 次のパッケージを付け加えます: 'MASS'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## select## 要求されたパッケージ betareg をロード中ですlibrary(MASS)
library(erer)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo使用する変数の準備をします。
data1 <- data1 %>% mutate(
fem=recode(sex,'1=0;2=1'),
marr=recode(ZQ50,'1=0;2:4=1'),
mar=recode(ZQ50,'2=1;1=0;3:4=0'),
age=ifelse(mbirth<=4,2007-ybirth,2007-ybirth-1),
edu3=recode(ZQ23A,'1:2=1;3:4=2;5:6=3;7:9=NA'),
ZQ23Ar=recode(ZQ23A,'1:2=1;3:4=2;5:6=3;7:9=NA'),
edu=as.factor(ifelse(ZQ24<=2, ZQ23Ar, NA)),
gend=as.factor(recode(ZQ39A,'1:2=1;3=2;4:5=3;6:9=NA')),
wtyp=as.factor(recode(JC_1,'1:2=1;3:8=2;10=3;11:99=NA')),
lsat=as.factor(recode(ZQ30D,'1=5;2=4;3=3;4=2;5=1;9=NA')),
lsatq=recode(ZQ30D,'1=5;2=4;3=3;4=2;5=1;9=NA'))リコードができているか確認します。
性別
table(data1$sex,data1$fem,useNA = "ifany")##
## 0 1
## 1 499 0
## 2 0 501既婚
table(data1$ZQ50,data1$marr,useNA = "ifany")##
## 0 1
## 1 637 0
## 2 0 345
## 4 0 18有配偶
table(data1$ZQ50,data1$mar,useNA = "ifany")##
## 0 1
## 1 637 0
## 2 0 345
## 4 18 0学歴
table(data1$ZQ24,data1$edu,useNA = "ifany")##
## 1 2 3 <NA>
## 1 220 278 271 3
## 2 26 15 22 0
## 3 0 0 0 144
## 9 0 0 0 21table(data1$ZQ23A,data1$edu,useNA = "ifany")##
## 1 2 3 <NA>
## 1 12 0 0 2
## 2 234 0 0 7
## 3 0 171 0 20
## 4 0 122 0 9
## 5 0 0 264 104
## 6 0 0 29 15
## 7 0 0 0 2
## 9 0 0 0 9生活満足度
table(data1$ZQ30D,data1$lsat,useNA = "ifany")##
## 1 2 3 4 5 <NA>
## 1 0 0 0 0 163 0
## 2 0 0 0 431 0 0
## 3 0 0 257 0 0 0
## 4 0 103 0 0 0 0
## 5 37 0 0 0 0 0
## 9 0 0 0 0 0 9性別役割意識
table(data1$ZQ39A,data1$gend,useNA = "ifany")##
## 1 2 3 <NA>
## 1 47 0 0 0
## 2 211 0 0 0
## 3 0 260 0 0
## 4 0 0 182 0
## 5 0 0 283 0
## 6 0 0 0 10
## 9 0 0 0 7現在の仕事
table(data1$JC_1,data1$wtyp,useNA = "ifany")##
## 1 2 3 <NA>
## 1 7 0 0 0
## 2 489 0 0 0
## 3 0 150 0 0
## 4 0 31 0 0
## 5 0 5 0 0
## 6 0 18 0 0
## 7 0 24 0 0
## 8 0 6 0 0
## 10 0 0 126 0
## 11 0 0 0 48
## 12 0 0 0 92
## 99 0 0 0 4ロジットモデルを推定します。
logit <- glm(marr~age+fem+edu+wtyp+gend, family=binomial(link=logit), data=data1)
summary(logit)##
## Call:
## glm(formula = marr ~ age + fem + edu + wtyp + gend, family = binomial(link = logit),
## data = data1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2188 -0.8401 -0.4370 0.9628 2.2546
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.20751 0.73982 -11.094 < 2e-16 ***
## age 0.27882 0.02406 11.588 < 2e-16 ***
## fem 0.41360 0.18383 2.250 0.02446 *
## edu2 -0.42951 0.20962 -2.049 0.04046 *
## edu3 -0.54938 0.20914 -2.627 0.00862 **
## wtyp2 -0.16154 0.19900 -0.812 0.41694
## wtyp3 1.29774 0.25974 4.996 5.84e-07 ***
## gend2 -0.19767 0.22517 -0.878 0.38000
## gend3 -0.23543 0.20269 -1.162 0.24541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1120.21 on 817 degrees of freedom
## Residual deviance: 892.23 on 809 degrees of freedom
## ( 182 個の観測値が欠損のため削除されました )
## AIC: 910.23
##
## Number of Fisher Scoring iterations: 4マクファデンの疑似決定係数を計算します。
1-(logit$deviance/logit$null.deviance)## [1] 0.2035181プロビットモデルを推定します。
probit <- glm(marr~age+fem+edu+wtyp+gend, family=binomial(link=probit), data=data1)
summary(probit)##
## Call:
## glm(formula = marr ~ age + fem + edu + wtyp + gend, family = binomial(link = probit),
## data = data1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2628 -0.8482 -0.4097 0.9670 2.3221
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.00032 0.42100 -11.877 < 2e-16 ***
## age 0.16950 0.01368 12.394 < 2e-16 ***
## fem 0.24778 0.10925 2.268 0.02332 *
## edu2 -0.24717 0.12467 -1.983 0.04740 *
## edu3 -0.32356 0.12478 -2.593 0.00951 **
## wtyp2 -0.09080 0.11826 -0.768 0.44260
## wtyp3 0.78092 0.15180 5.145 2.68e-07 ***
## gend2 -0.12018 0.13402 -0.897 0.36989
## gend3 -0.14600 0.12065 -1.210 0.22624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1120.21 on 817 degrees of freedom
## Residual deviance: 889.94 on 809 degrees of freedom
## ( 182 個の観測値が欠損のため削除されました )
## AIC: 907.94
##
## Number of Fisher Scoring iterations: 5マクファデンの疑似決定係数を計算します。
1-(probit$deviance/probit$null.deviance)## [1] 0.2055625線形確率モデルを推定します。
lpm <- lm(marr~age+fem+edu+wtyp+gend, data=data1)
summary(lpm)##
## Call:
## lm(formula = marr ~ age + fem + edu + wtyp + gend, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.96830 -0.33399 -0.05329 0.38846 0.97117
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.091854 0.118617 -9.205 < 2e-16 ***
## age 0.054351 0.003872 14.037 < 2e-16 ***
## fem 0.075950 0.033837 2.245 0.02507 *
## edu2 -0.085842 0.038726 -2.217 0.02692 *
## edu3 -0.116948 0.039096 -2.991 0.00286 **
## wtyp2 -0.028089 0.036592 -0.768 0.44292
## wtyp3 0.248055 0.046287 5.359 1.09e-07 ***
## gend2 -0.035842 0.041590 -0.862 0.38906
## gend3 -0.043555 0.037519 -1.161 0.24603
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4318 on 809 degrees of freedom
## ( 182 個の観測値が欠損のため削除されました )
## Multiple R-squared: 0.2497, Adjusted R-squared: 0.2423
## F-statistic: 33.65 on 8 and 809 DF, p-value: < 2.2e-16表１２－１ 結婚経験の有無に関する推定結果
stargazer(logit, probit, lpm, type="text",star.cutoffs=c(0.1,0.05,0.01),keep.stat = c("n","ll","rsq"))##
## ============================================
## Dependent variable:
## -----------------------------
## marr
## logistic probit OLS
## (1) (2) (3)
## --------------------------------------------
## age 0.279*** 0.170*** 0.054***
## (0.024) (0.014) (0.004)
##
## fem 0.414** 0.248** 0.076**
## (0.184) (0.109) (0.034)
##
## edu2 -0.430** -0.247** -0.086**
## (0.210) (0.125) (0.039)
##
## edu3 -0.549*** -0.324*** -0.117***
## (0.209) (0.125) (0.039)
##
## wtyp2 -0.162 -0.091 -0.028
## (0.199) (0.118) (0.037)
##
## wtyp3 1.298*** 0.781*** 0.248***
## (0.260) (0.152) (0.046)
##
## gend2 -0.198 -0.120 -0.036
## (0.225) (0.134) (0.042)
##
## gend3 -0.235 -0.146 -0.044
## (0.203) (0.121) (0.038)
##
## Constant -8.208*** -5.000*** -1.092***
## (0.740) (0.421) (0.119)
##
## --------------------------------------------
## Observations 818 818 818
## R2 0.250
## Log Likelihood -446.115 -444.970
## ============================================
## Note: *p<0.1; **p<0.05; ***p<0.01表１２－２
平均的な限界効果
logitmfx(marr~age+fem+edu+wtyp+gend, data=data1, atmean=FALSE) ## Call:
## logitmfx(formula = marr ~ age + fem + edu + wtyp + gend, data = data1,
## atmean = FALSE)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## age 0.0511405 0.0060835 8.4064 < 2.2e-16 ***
## fem 0.0759803 0.0335095 2.2674 0.023364 *
## edu2 -0.0778624 0.0371848 -2.0939 0.036266 *
## edu3 -0.1009618 0.0379023 -2.6637 0.007728 **
## wtyp2 -0.0296003 0.0363835 -0.8136 0.415896
## wtyp3 0.2465307 0.0471572 5.2279 1.715e-07 ***
## gend2 -0.0360236 0.0406958 -0.8852 0.376054
## gend3 -0.0432318 0.0371564 -1.1635 0.244622
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "fem" "edu2" "edu3" "wtyp2" "wtyp3" "gend2" "gend3"平均での限界効果
logitmfx(marr~age+fem+edu+wtyp+gend, data=data1, atmean=TRUE) ## Call:
## logitmfx(formula = marr ~ age + fem + edu + wtyp + gend, data = data1,
## atmean = TRUE)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## age 0.0676135 0.0057313 11.7971 < 2.2e-16 ***
## fem 0.0998524 0.0440362 2.2675 0.023359 *
## edu2 -0.1026049 0.0491060 -2.0895 0.036666 *
## edu3 -0.1304636 0.0483281 -2.6995 0.006944 **
## wtyp2 -0.0388863 0.0475030 -0.8186 0.413011
## wtyp3 0.3127506 0.0570833 5.4788 4.281e-08 ***
## gend2 -0.0475023 0.0535470 -0.8871 0.375017
## gend3 -0.0569446 0.0488466 -1.1658 0.243702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "fem" "edu2" "edu3" "wtyp2" "wtyp3" "gend2" "gend3"オッズ比
logitor(marr~age+fem+edu+wtyp+gend, data=data1) ## Call:
## logitor(formula = marr ~ age + fem + edu + wtyp + gend, data = data1)
##
## Odds Ratio:
## OddsRatio Std. Err. z P>|z|
## age 1.321563 0.031798 11.5878 < 2.2e-16 ***
## fem 1.512250 0.277997 2.2499 0.024455 *
## edu2 0.650828 0.136426 -2.0490 0.040462 *
## edu3 0.577310 0.120737 -2.6269 0.008617 **
## wtyp2 0.850834 0.169317 -0.8117 0.416938
## wtyp3 3.661019 0.950903 4.9964 5.842e-07 ***
## gend2 0.820639 0.184780 -0.8779 0.380002
## gend3 0.790228 0.160169 -1.1616 0.245413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1順序のある３値以上の質的従属変数の回帰
線形モデルを推定します。
lm <- lm(lsatq~age+fem+mar+edu+wtyp, data=data1)
summary(lm)##
## Call:
## lm(formula = lsatq ~ age + fem + mar + edu + wtyp, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0806 -0.5705 0.1548 0.6149 1.9618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.872674 0.270837 14.299 < 2e-16 ***
## age -0.023841 0.009398 -2.537 0.011375 *
## fem 0.221707 0.074024 2.995 0.002826 **
## mar 0.535325 0.077356 6.920 9.1e-12 ***
## edu2 0.261437 0.084933 3.078 0.002152 **
## edu3 0.246689 0.084763 2.910 0.003708 **
## wtyp2 -0.291406 0.080600 -3.615 0.000318 ***
## wtyp3 -0.283835 0.103557 -2.741 0.006262 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9573 on 819 degrees of freedom
## ( 173 個の観測値が欠損のため削除されました )
## Multiple R-squared: 0.09269, Adjusted R-squared: 0.08493
## F-statistic: 11.95 on 7 and 819 DF, p-value: 1.442e-14順序ロジットモデルを推定します。
ologit <- polr(lsat~age+fem+mar+edu+wtyp, data=data1, method = "logistic")
summary(ologit)##
## Re-fitting to get Hessian## Call:
## polr(formula = lsat ~ age + fem + mar + edu + wtyp, data = data1,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## age -0.04276 0.01835 -2.330
## fem 0.43123 0.14372 3.000
## mar 1.03753 0.15533 6.680
## edu2 0.50194 0.16581 3.027
## edu3 0.47593 0.16467 2.890
## wtyp2 -0.53462 0.15517 -3.445
## wtyp3 -0.56579 0.20313 -2.785
##
## Intercepts:
## Value Std. Error t value
## 1|2 -3.8654 0.5579 -6.9287
## 2|3 -2.3692 0.5339 -4.4374
## 3|4 -0.8926 0.5279 -1.6907
## 4|5 1.3020 0.5308 2.4529
##
## Residual Deviance: 2177.539
## AIC: 2199.539
## ( 173 個の観測値が欠損のため削除されました )data1$used <- TRUE;
data1$used[na.action(ologit)] <- FALSE;
table(data1$used)##
## FALSE TRUE
## 173 827ologitn <- polr(lsat~1, data=data1, subset=used==1, method = "logistic")
summary(ologitn)##
## Re-fitting to get Hessian## Call:
## polr(formula = lsat ~ 1, data = data1, subset = used == 1, method = "logistic")
##
## No coefficients
##
## Intercepts:
## Value Std. Error t value
## 1|2 -3.2456 0.1831 -17.7289
## 2|3 -1.7932 0.0994 -18.0357
## 3|4 -0.4095 0.0710 -5.7668
## 4|5 1.6521 0.0947 17.4540
##
## Residual Deviance: 2256.12
## AIC: 2264.12マクファデンの疑似決定係数を計算します。
1-(logLik(ologit)/logLik(ologitn))## 'log Lik.' 0.0348303 (df=11)表１２－３ 生活満足度の推定結果
stargazer(ologit, lm, type="text",star.cutoffs=c(0.1,0.05,0.01),keep.stat = c("n","rsq"))##
## =========================================
## Dependent variable:
## ----------------------------
## lsat lsatq
## ordered OLS
## logistic
## (1) (2)
## -----------------------------------------
## age -0.043** -0.024**
## (0.018) (0.009)
##
## fem 0.431*** 0.222***
## (0.144) (0.074)
##
## mar 1.038*** 0.535***
## (0.155) (0.077)
##
## edu2 0.502*** 0.261***
## (0.166) (0.085)
##
## edu3 0.476*** 0.247***
## (0.165) (0.085)
##
## wtyp2 -0.535*** -0.291***
## (0.155) (0.081)
##
## wtyp3 -0.566*** -0.284***
## (0.203) (0.104)
##
## Constant 3.873***
## (0.271)
##
## -----------------------------------------
## Observations 827 827
## R2 0.093
## =========================================
## Note: *p<0.1; **p<0.05; ***p<0.01ologitやstargazerでは閾値の有意性が表示できないので、念のため以下で確認します。
ctable <- coef(summary(ologit))##
## Re-fitting to get Hessianp <- pnorm(abs(ctable[, "t value"]), lower.tail=FALSE)*2
(ctable <- cbind(ctable, "p value"=p))## Value Std. Error t value p value
## age -0.0427615 0.0183502 -2.330302 1.979020e-02
## fem 0.4312341 0.1437242 3.000428 2.696009e-03
## mar 1.0375291 0.1553288 6.679568 2.396470e-11
## edu2 0.5019418 0.1658078 3.027251 2.467887e-03
## edu3 0.4759306 0.1646740 2.890138 3.850731e-03
## wtyp2 -0.5346155 0.1551677 -3.445405 5.702034e-04
## wtyp3 -0.5657916 0.2031318 -2.785342 5.347129e-03
## 1|2 -3.8654154 0.5578869 -6.928672 4.248080e-12
## 2|3 -2.3691896 0.5339114 -4.437421 9.104315e-06
## 3|4 -0.8926044 0.5279479 -1.690706 9.089303e-02
## 4|5 1.3019746 0.5307890 2.452904 1.417080e-02表１２－４ 生活満足度に対する各属性のオッズ比
exp(coef(ologit))## age fem mar edu2 edu3 wtyp2 wtyp3
## 0.9581399 1.5391558 2.8222349 1.6519259 1.6095113 0.5858945 0.5679104表１２－５ 生活満足度に対する平均での限界効果
me_ol <- ocME(w = ologit, digits=4)##
## Re-fitting to get Hessianme_ol## effect.1 effect.2 effect.3 effect.4 effect.5
## age 0.0013 0.0035 0.0054 -0.0049 -0.0053
## fem -0.0136 -0.0353 -0.0539 0.0490 0.0537
## mar -0.0307 -0.0799 -0.1274 0.0987 0.1394
## edu2 -0.0147 -0.0389 -0.0637 0.0511 0.0662
## edu3 -0.0140 -0.0370 -0.0604 0.0488 0.0626
## wtyp2 0.0189 0.0472 0.0638 -0.0685 -0.0613
## wtyp3 0.0214 0.0522 0.0650 -0.0772 -0.0615